Development and validation of molecular biomarkers for the green-lipped mussel (Perna canaliculus)

ABSTRACT

Globally, there is a move towards using local, native species for ecotoxicological risk assessments. Anthropogenic stressors from urban, agricultural, and industrial activities can impact the health of receiving ecosystems. Biomarkers can provide valuable insights as early warning signals of the potential environmental impacts of stressors. The aim of this study was to develop biomarkers in the green-lipped mussel (Perna canaliculus), a potential bioindicator of environmental health for coastal marine ecosystems in New Zealand. Reverse transcription quantitative polymerase chain reaction (RT-qPCR) assays targeting the expression of genes involved in oxidative stress, xenobiotic transfer, membrane transportation, cellular and DNA response/repair, and endocrine disruption were developed and validated for P. canaliculus. We found significant modulation of genes associated with oxidative stress, xenobiotic transfer, membrane transport, endocrine disruption, and genotoxicity in P. canaliculus following 48-hour exposures to copper and benzo[α]pyrene. These results demonstrate the potential of P. canaliculus as a bioindicator species for environmental risk assessment. The gene expression assays showed potential as early indicators of exposure to the chemicals tested but require additional validation to assess their ability to predict effects at higher levels of biological organisation.

KEYWORDS:

• Ecotoxicology
• copper
• benzo[α]pyrene
• RT-qPCR
• gene expression
• New Zealand

Introduction

Integration of gene expression analysis has become a valuable approach to complement chemical and ecological methods to assess the biological effects of pollution. Its inclusion into environmental risk assessment frameworks can achieve more effective biomonitoring to manage and protect ecosystems (Snape et al. 2004; Rudén et al. 2017). Reverse transcription quantitative polymerase chain reaction (RT-qPCR), the most commonly used technique to measure gene expression, is a fast and low-cost method recognised for its efficacy, sensitivity, and reproducibility (Bustin et al. 2009). Gene expression data are used in ecotoxicology to better understand mechanisms of toxicity, identify possible biological responses to contaminant stress, and predict potential adverse effects of toxicants linked to higher levels of biological organisation (Wood et al. 2013). RT-qPCR has been predominantly used in laboratory settings on model organisms. However, there has been a shift to integrate non-standard test species into regulatory frameworks due to the challenges of characterising the hazards of contaminants of emerging concern (Rudén et al. 2017). Despite a rapidly growing availability of nucleotide and protein sequences, there is still limited data for ‘non-model’ organisms (Martyniuk and Simmons 2016).

Bivalve mollusks, particularly mussels, are globally used as sentinel organisms to monitor the impacts of contaminants in coastal ecosystems (Goldberg 1975; Beyer et al. 2017; Strehse and Maser 2020). Mussels fulfil almost all the requirements as a valuable sentinel species, as: (i) their global distribution allows comparison among a wide range of different ecosystems, (ii) they tolerate environmental fluctuations (such as varying salinity, temperature and oxygen levels), (iii) their sessile nature allows for reliable assessments of local conditions, and (iv) as filter feeders, they are indiscriminately exposed to contaminants that can bioaccumulate in their tissues (Bayne et al. 1976; Strehse and Maser 2020). Their use in ecotoxicology is well established and they have been incorporated into monitoring programs such as the U.S. Mussel Watch program (Goldberg 1975), Biological Effects of Environmental Pollution in Marine Coastal Ecosystem (EU BEEP project) (Lehtonen et al. 2006), and the European BioMar initiative (Narbonne et al. 1999). Mytilus spp. (blue mussels) are commonly used for biomonitoring and ecotoxicological hazard assessments worldwide using a wide array of biomarkers to determine sublethal effects of contaminants (e.g. ROS production, biochemical markers, molecular biomarkers) (Monserrat et al. 2007; Beyer et al. 2017; Martínez-Gómez et al. 2017; Strehse and Maser 2020) whereas other mussel species with less prevalent distributions have not been fully explored as potential bioindicators of environmental health.

The green-lipped mussel (Perna canaliculus, Gmelin 1791) is a New Zealand endemic species of cultural importance to Māori as a taonga, or treasured species. It is commercially important as the country's most significant aquaculture species (Webb et al. 2020). Historically, New Zealand has experienced relatively low levels of contamination compared to the rest of the world, but it is a developed country and the expanding population combined with intensification of urban, agricultural, and industrial activities are increasing contaminant loadings into the receiving environment. Current risk assessment by direct toxicity assessment (DTA) have started incorporating tests using local, native species. For instance, a recent study has compared the sensitivity of early life stages of Mytilus galloprovincialis and P. canaliculus and found that both species responded comparably to triclosan exposure (Rolton et al. 2022). M galloprovincialis coexists in intertidal and shallow subtidal zones across New Zealand but P. canaliculus has a broader ecological niche making it a more relevant receptor organism. It can thrive in soft sediment, which may lead to P. canaliculus being exposed to sediment-bound contaminants that M. galloprovincialis may not interact with (Morton and Miller 1968; Gardner 2000). While P. canaliculus has been previously used in studies assessing contaminant uptake (Chandurvelan et al. 2012; Webb et al. 2020), its use in ecotoxicology is yet to be fully realised. Investigating the application of P. canaliculus as a bioindicator species would provide New Zealand with a novel endemic bivalve species for application in environmental risk assessment processes (Breitholtz et al. 2006).

The aim of this study was to develop molecular tools in P. canaliculus for assessment of contaminant impacts in marine ecosystems. Development and validation of biomarkers associated with generic ecotoxicological responses is a key first step for assessing whether P. canaliculus can be used to characterise a broader range of chemical contaminants likely present in marine ecosystems. We designed assays that targeted a suite of genes encoding for oxidative stress responses, xenobiotic transfer, membrane transport, cellular and DNA response/repair, and endocrine disruption mechanisms in P. canaliculus, using RT-qPCR. The genes selected have previously been measured in other mussels from the Mytilidae family (e.g. Chatel et al. 2018). After the initial primer development, P. canaliculus was exposed to common environmental contaminants, benzo[α]pyrene (B[α]P) and copper (Cu²⁺), as they are well established in the literature, persistent and widespread urban contaminants, and previous studies demonstrated they induce genotoxicity, oxidative stress, enzymatic activity, and endocrine disruption in other mussel species (e.g. Banni et al. 2017; Vernon and Jha 2019). Gills and digestive gland were selected to assess the validity of the developed primers. Gills were chosen because they are the first line of exposure and defense to contaminants in the dissolved phase, serving as the interface between the mussel and the aquatic environment. The digestive gland was selected because filtered and ingested particulate matters accumulate here, making it a dominant organ for metabolism and biotransformation activities (Chipman et al. 1991). The development and validation of biomarkers associated with responses to B[α]P and Cu²⁺ are critical first steps for assessing whether P. canaliculus can be used to characterise a broader range of chemical contaminants present in coastal marine zones.

Materials and methods

Preparation of test solutions

A stock solution of soluble copper (Copper (II) sulfate pentahydrate, CAS # 7758-99-8) was prepared in ultrapure water at 1 g L⁻¹. A stock solution of B[α]P (CAS # 50-32-8) (1 g L⁻¹) was prepared in dimethyl sulfoxide (DMSO) (CAS # 67-68-5) as previous research has indicated DMSO induced limited effects on molecular endpoints in mussels and is compliant with OECD test guidelines (Chatel et al. 2012; Guo et al. 2020).

Mussel collection and laboratory maintenance

Early adult P. canaliculus (average shell length and weight ± standard deviation: 5.13 cm ± 0.46, 10.3 g ± 2.16) were collected from aquaculture farms in Pelorus Sound, New Zealand. Mussels were transported to a laboratory flow through system with flow rate greater than 1 L/h/organism (mASTM 2014) and acclimated for two weeks in natural seawater (filtered through a 0.22 µm filter and sterilised with UV-light), maintained at 33 ± 1 practical salinity units (PSU ± standard deviation), 17 °C, pH 8, dissolved oxygen above 90%, and with a 12:12 h photoperiod. Mussels were fed 1 L (6 × 10⁴ cells mL⁻¹) of a microalgae blend of 50:50 Chaetoceros muelleri and Tisochrysis lutea three times a week. Two days before the start of exposure, mussels were not fed to allow for depuration and reduce algal biomass in the digestive gland.

Experimental design

For Cu²⁺ exposure, 48 mussels (3 individuals × 4 contaminant concentrations × 4 replicates per condition) were selected and pooled in 3 L jars (1 mussel/L). They were exposed to four test conditions: (i) seawater control, (ii) 25 μg L⁻¹ Cu²⁺, (iii) 50 μg L⁻¹ Cu²⁺, and (iv) 100 μg L⁻¹ Cu²⁺. For a separate B[α]P experiment with an identical set-up, mussels were exposed to (i) seawater control, (ii) 15 μg L⁻¹ B[α]P, (iii) 30 μg L⁻¹ B[α]P, (iv) 60 μg L⁻¹ B[α]P, and (v) a DMSO solvent control (final concentration 0.01% in all treatments except for seawater controls). These test concentrations of Cu²⁺ and B[α]P were selected as previous research demonstrated responses over these concentration ranges in other mussel species (Chatel et al. 2012, 2018; Banni et al. 2017). Jars were aerated throughout the experiment using glass Pasteur pipettes attached to an air pump, filled with sterilised seawater, and maintained under lab conditions consistent with the main culture. All exposures lasted 48 h, during which time the mussels were not fed; a complete water change and chemical renewal for all test conditions were carried out after 24 h. We selected 48 h as our sampling time point after a preliminary time duration experiment. Briefly, control and exposed mussels were sampled at 24, 48, and 72 h and gene expression measured (Supp Figures 1 and 2). Immediately following 48 h of exposure, mussels were removed from the treatment solutions and the digestive glands and gills were dissected. Tissue samples were weighed and snap-frozen in liquid nitrogen. Dissected tissue was homogenised in liquid nitrogen using a Mixer Mill MM 400 (Retsch) and stored in a −80 °C freezer before RNA extraction.

RNA extraction and cDNA synthesis

Frozen tissue powder was removed from −80 °C storage and 1 mL of TRI Reagent® (Sigma-Aldrich, Auckland, NZ) added immediately to arrest RNA degradation. Total RNA was then extracted following the TRI Reagent® manufacturer’s protocol, with the addition of a second ethanol wash step to ensure unwanted organic compounds (i.e. trizol) had been removed from the RNA sample. Residual genomic DNA was removed by treating the total RNA with Turbo DNAse, following the manufacturer's protocol (Invitrogen, Auckland, NZ). RNA concentration and purity were measured using a NanoPhotometer® N60/N50 (Implen GmbH, Munich, Germany). First-strand cDNA synthesis was conducted using 5 μg of total RNA extract with random hexamers, according to the SuperScript™ IV First-Strand Synthesis System (Invitrogen, Auckland, NZ) protocol. Samples were stored at −20 °C until RT-qPCR analyses.

RT-qPCR primer design

To generate DNA sequences for the gene regions of interest, degenerate primers were first designed by aligning previously published sequences from the subfamily Mytilinae obtained from GenBank® (www.ncbi.nlm.nih.gov/genbank/) in Geneious Prime version 2020.2. Several degenerate primers were developed for each target gene, and their efficacy was tested in an Eppendorf Mastercycler Nexus. The primer combinations were tested on an amplification gradient from 50–60°C, and the conditions were as follows an initial denaturation at 95 °C for 5 min, followed by 39 cycles of denaturation at 94 °C for 30 s, annealing on a gradient of 50–60°C for 30 s, and extension at 72 °C. The extension time is dependent on length of amplicon (calculated based on 60 s per kilobase). Successful PCR products were purified and sent for Sanger sequencing in both directions by an external contractor (Genetic Analysis Services, University of Otago, Dunedin, NZ). Forward and reverse sequences were aligned using Geneious Prime and conflicts (i.e. ambiguous nucleotide base calls) were resolved by manual inspection. After removing low-quality or ambiguous DNA sequence reads, the remaining sequences were compared to existing sequences in GenBank using the BLAST online software (http://blast.ncbi.nlm.nih.gov/Blast.cgi) to ensure the correct target gene had been amplified. Additional information on degenerate primer design can be found in the supplemental material.

After confirming the DNA sequences of gene regions of interest following end-point PCR, RT-qPCR primers were designed using both Geneious Prime and the PrimerQuest® Tool (Integrated DNA Technologies, Inc. 2020) with the following criteria: primer size between 18 and 24 base pairs, GC content between 40% and 60%, amplicon size from 90 to 130 base pairs, and primer annealing temperatures in the 55–62 °C range. Following the Minimum Information for Publication of Quantitative Real-Time PCR Experiments (MIQE) guidelines, RT-qPCR assays were validated by assessing optimal annealing temperatures of the primers, optimal primer concentration, and reaction efficiency using decreasing ten-fold dilutions of cDNA (from 25 ng L⁻¹ to 0.0025 ng L⁻¹) (Bustin et al. 2009; Taylor et al. 2010; Bustin and Huggett 2017) on a Rotor-Gene Q machine (QIAGEN GmbH, Hilden, Germany) using a 72-sample rotor. Following the amplification program, melt curve analyses were performed from 55 to 95 °C to determine the specificity of the PCR amplicon.

RT-qPCR and reference gene selection

To investigate the efficacy of the developed primers (Table 1) to amplify the target biomarker, cDNA from both the Cu²⁺ and B[α]P exposure experiments was assessed by RT-qPCR. Reactions were performed in 20 μL volumes containing 10 μL PowerUp™ SYBR™ Green Master Mix, 7 μL DNase-free water, 1 μL cDNA diluted fourfold, and 1 μL of each primer diluted in the final reaction concentration to 0.5 μM except for actin (0.5 μM forward primer, 0.25 μM reverse primer) and Na/K ATPase specific primers (0.25 μM forward primer, 0.5 μM reverse primer), which performed better at these dilutions. Thermal cycling conditions consisted of an initial denaturation at 95 °C for 2 min, followed by 40 cycles of denaturation at 95 °C for 30 s, annealing for 30 s (see Table 1 for primer-specific temperatures), and extension at 72 °C for 15 s. Each gene was analyzed in triplicate and included no template control reactions. Data were reported as normalised relative expression (fold change) to control samples. This was done by recording the Ct values and calculating the ΔΔCt value as 2^–ΔΔCt (ΔΔCt = ΔCt (treated sample) – ΔCt (seawater control sample) and ΔCt = Ct (gene of interest) – Ct (housekeeping gene)). 18S ribosomal RNA (rRNA), 28S rRNA, actin, and elongation factor genes were selected as potential reference genes. To assess the stability of the candidate reference genes across the test, two different algorithms were used, normFinder (Andersen et al. 2004) and BestKeeper (Pfaffl et al. 2004).

Table 1. Primers used in the RT-qPCR assays with their respective amplicon size, annealing temperatures, and amplification efficiency.

Gene symbol	Gene name	Gene Function	Accession number	Primer sequences (5’-3’)	Amplicon size (bp)	Annealing temp. (°C)	Amplification efficiency
18S	18 S ribosomal RNA	Reference gene	MK419107	CGA ACG AGA CTC TAG CCT ACT A (F) GAA CAT CTA AGG GCA TCA CAG A (R)	129	60	0.98
28S	28 S ribosomal RNA	Reference gene	MK419107	CGA ATG AAG CGG GAA GAG CC (F) CGA TCA GGA GGA CTT GGG CA (R)	115	55	1.04
ACT	Actin	Reference gene	OK318446	CAG CCA TGT ACG TCG CCA TC (F) TTG GTA CGG TGT GGG TGA CA (R)	105	62	0.91
EF1	Elongation factor 1	Reference gene	OK318444	GGA TGG GCT ATT GAG AGG AAA G (F) GGT TTG TCT GTT GGT CTT TGT G (R)	98	62	0.97
MRP	Multidrug resistance-associated protein	Membrane transport	OK318439	GTG GCG GAC AGA AAC AGA GAG (F) CCA ACG TGG GAA TCT ACA GCG (R)	103	55	1.08
P-gp	P-glycoprotein	Membrane transport	OK318437	TGC CCA CAG ACT ATC GAC CAT C (F) AGT TCA TCG TGT GTA CCT CTC TCG (R)	90	55	1.08
Na/K ATPase	Sodium/potassium transporting ATPase	Membrane transport	OK318435	ACT GGT CAG CGG GCA AGA AG (F) TCC ACT TGA CGA CAA CTT CAG GG (R)	101	55	0.95
MT10	Metallothionein 10	Metal binding	OK318440	GGA GAT ATG CCT ACA CCT TGC (F) CAC AAT CGC TGC CAC ATT TAC (R)	118	55	1.01
HSP70	Heat shock protein 70	Generic stress	OK318441	CAG CCT GGT GTG TTG ATC CA (F) TGC TGG TGG GAT TCC TGT CA (R)	99	55	1.05
ER1	Estrogen receptor 1	Endocrine disruption	OK318443	ATG ATC CGT CAG AAC CAT ATC C (F) GTG AAC CTG GAA GTC CCA TTA (R)	117	62	0.98
CAT	Catalase	Oxidative stress	OK318445	CTT TGC CGA AGT TGA ACA GAT AG (F) TGG GTA TCA GAG TAG GAG AAG AG (R)	111	58	1.03
SOD	Superoxide dismutase	Oxidative stress	OK318434	GTG ACG GTG CTG TTA CAG GAA (F) GTT CAC CAG GAG CTA AGC CAG (R)	99	62	1.01
GST	Glutathione S-transferase	Xenobiotic transfer	OK318442	CTG GAC AGT TTC CCT AAC CTA AA (F) ATT CCC ATT GAT TGG CAT GTT C (R)	111	62	0.99
p53	p53 tumour suppressor-like protein	Cellular response/repair	OK318438	GGG TGT GTC ATC CGA GCC AT (F) CTC CTT GGA CGT GGC ATG GT (R)	93	60	1.02
Rad51	DNA repair protein Rad51	Cellular response/repair	OK318436	GGG TGG TGG AGA AGG AAA AGC (F) ACA CAT CAC TGC CAG ACA AAC C (R)	107	62	1.01

Confirmation of chemical exposure concentrations

To verify Cu²⁺ concentrations, a 100 mL aliquot was taken from each exposure solution, placed in sterile bottles provided by a commercial testing facility, Hill Laboratories, and stored at 4 °C until being sent for analysis. Upon receipt at the testing laboratory, water extracts were filtered (0.45 μm membrane filter) to isolate dissolved phase copper and the pH was adjusted to <2 with nitric acid for preservation. Water samples were digested in acid and residues of copper were measured by Inductively Coupled Plasma-Mass Spectrometry using a Perkin Elmer NexION 300 ICP-MS or 2000 ICP-MS using a recognised standard method.

To verify B[α]P concentrations, 5 mL aliquots were sampled from the exposure treatment and transferred to 7 mL amber glass EPA vials. Dichloromethane (DCM; 2 mL) was added to sterilise and stabilise residues of B[α]P in the aqueous solutions. Vials were firmly shaken to cause a primary extraction of B[α]P. The vials of partially extracted solutions were stored at 4 °C in the dark until complete extraction. Before total extraction, the solution vials were brought to room temperature and spiked with a surrogate standard of perylene-d₁₂ at a concentration of 10 or 50 μg L⁻¹, respectively, for lower and higher concentration treatments. The solutions were extracted in four separate batches and a quality assurance solvent blank was included in each batch. The vials of solution were extracted with the original aliquot of DCM by vortex mixing for 10 s, followed by standing to allow the separation of the aqueous and organic phases. The lower layer of DCM was carefully removed using a glass Pasteur pipette and transferred into a tapered glass test tube. A fresh 2 mL aliquot of DCM was added to the original sample solution, extracted, and transferred to the tapered test tube. This extraction procedure was repeated for a final third time.

The DCM sample extracts were each dried through a glass column containing 3 g of anhydrous sodium sulphate and collected in a glass round bottom flask. The tapered test tubes were washed with three individual 1 mL aliquots of DCM that were similarly passed through the drying column. Isooctane (0.5 mL) was added to the round bottom flasks and the DCM removed by rotary evaporation (Bucchi, water bath 40°C, condenser −11°C, vacuum 180 mb). An internal standard solution of benzo[α]pyrene-d₁₂ was added to the round bottom flasks. The contents were transferred to GC vials, which were capped and stored under refrigeration before analysis by gas chromatography mass-spectrometry.

Residues of B[α]P were analyzed using an Agilent 6890 gas chromatograph and 5975 MSD operating in electron impact mode and 70 eV, with source and quadrupole temperatures of 300 and 200°C, respectively, and a transfer line temperature of 300°C. Aliquots (1 μL) of calibration standards and sample extracts were injected into a split/less injector operating in pulsed splitless mode (40 psi for 1 min followed by purge flow of 50 mL min⁻¹) and a temperature of 300 ^oC. Sample components were separated on an Agilent J&W DB-5MS ultra inert column (30 m x 0.25 mm x 0.25 μm) using helium as carrier gas at a constant flow rate of 1.5 mL min⁻¹. Data acquisition was in synchronous SIM/Scan mode between 50 and 450 amu with a solvent delay of 18 mins. Mass specific ions were acquired for the internal and surrogate standards and benzo[a]pyrene. Data quantitation was completed using Agilent Chemstation data analysis software and internal standard quantitation. Calibration curves were prepared using eight calibration standards spanning concentrations from 2.5 to 500 μg L⁻¹.

Statistical analyses

Statistical analyses were performed using RStudio Team (2020), incorporating the dunn.test package. The distribution of datasets was checked for normality using a Shapiro–Wilk test. Following a significant Kruskal–Wallis test, the Dunn test was applied to compare data between treatments. Statistical differences were accepted when p ≤ 0.05. Additionally, partial least squares–discriminant analyses (PLS-DA) were completed using mixOmics, ade4 and RVAideMemoire packages to understand the expression of which genes discriminate treatment exposures (Jaumot et al. 2015). Genes were assigned Variable Importance on Projection (VIP) scores that rank their significance in the PLS model.

Results

Gene expression analysis

NormFinder (Supp Figure 3) and BestKeeper identified EF1 as the most stable reference gene while ACT was the second most stable. Gene expression results were then calculated relative to control samples and EF1 as the reference gene. Results for gene expression for both the Cu²⁺ and B[α]P conditions are summarised in Table 2. Specific figures for gene expression in the gills and digestive gland for both exposures can be seen in the supplemental material (Supp Figures 4–7).

Table 2. Gene expression results relative to control samples and reference gene EF1 expressed as P-values calculated using Kruskal-Wallis and Dunn tests for copper and benzo[α]pyrene. *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001.

Tissue	Copper (μg L^−1)
	Digestive gland				Gills
Gene	Solvent Control	25 μg L^−1	50 μg L^−1	100 μg L^−1	Solvent control	25 μg L^−1	50 μg L^−1	100 μg L^−1
CAT	–	0.28	0.30	0.25	–	0.47	0.03* (downregulated)	0.02* (downregulated)
SOD	–	0.13	0.09	0.08	–	0.07	0.15	0.006** (downregulated)
GST	–	0.09	0.33	0.19	–	0.44	0.21	0.08
MRP	–	0.30	0.44	0.005** (downregulated)	–	0.15	0.01** (downregulated)	0.005** (downregulated)
P-gp	–	0.47	0.47	0.02* (downregulated)	–	0.50	0.08	0.003** (downregulated)
Na/K ATPase	–	0.30	0.04* (upregulated)	0.01** (upregulated)	–	0.50	0.04* (downregulated)	0.007** (downregulated)
p53	–	0.30	0.07	0.001*** (downregulated)	–	0.23	0.01** (downregulated)	0.0005*** (downregulated)
Rad51	–	0.19	0.47	0.04* (downregulated)	–	0.13	0.06	0.009** (downregulated)
MT10	–	0.15	0.01** (upregulated)	0.01** (upregulated)	–	0.12	0.02* (upregulated)	0.08
HSP70	–	0.41	0.01** (upregulated)	0.09	–	0.19	0.19	0.19
ER1	–	0.21	0.05* (downregulated)	0.13	–	0.30	0.25	0.38
	Benzo[α]pyrene (μg L^−1)
		15 μg L^−1	30 μg L^−1	60 μg L^−1		15 μg L^−1	30 μg L^−1	60 μg L^−1
CAT	0.48	0.34	0.25	0.12	0.38	0.14	0.48	0.13
SOD	0.16	0.02* (upregulated)	a0.07	0.31	0.36	0.24	0.20	0.14
GST	0.19	0.27	0.29	0.10	0.20	0.50	0.20	0.50
MRP	0.19	0.01** (upregulated)	0.02* (upregulated)	0.05* (upregulated)	0.12	0.24	0.15	0.38
P-gp	0.22	0.04* (upregulated)	0.0008*** (upregulated)	0.0006*** (upregulated)	0.12	0.06	0.09	0.09
Na/K ATPase	0.41	0.06	0.15	0.17	0.45	0.32	0.45	0.36
p53	0.22	0.40	0.36	0.18	0.08	0.24	0.07	0.07
Rad51	0.18	0.02* (upregulated)	0.29	0.48	0.40	0.39	0.13	0.39
MT10	0.19	0.04* (downregulated)	0.004** (downregulated)	0.07	0.03* (upregulated)	0.27	0.45	0.37
HSP70	0.30	0.25	0.42	0.28	0.25	0.22	0.27	0.17
ER1	0.20	0.18	0.004** (upregulated)	0.04* (upregulated)	0.18	0.42	0.08	0.39

Copper exposure

Target genes in the gills and digestive gland showed no significant modulation (i.e. altered expression; Table 2) at 25 µg L⁻¹ Cu²⁺. However, both oxidative stress genes were significantly downregulated in gill tissue exposed to higher concentrations; CAT was downregulated at 50 and 100 µg L⁻¹ Cu²⁺ and SOD at 100 µg L⁻¹ Cu²⁺. GST mRNA levels were not significantly modulated and only a trend of downregulation was observed. The expression of genes encoding for three membrane transporters (MRP, pgp, and Na/K ATPase) was significantly downregulated at 50 and 100 μg L⁻¹ Cu²⁺ except for pgp, which was only significantly downregulated at 100 µg L⁻¹. There was also significant downregulated in the expression of cell cycle control genes; p53 at 50 and 100 μg L⁻¹ Cu²⁺ and Rad51 at 100 μg L⁻¹ Cu²⁺. MT10 was significantly upregulated at 50 μg L⁻¹ Cu²⁺. Expression of HSP70 and ER1 genes in the gills did not vary significantly across the different Cu²⁺ concentrations.

No significant modulation of expression was observed within the digestive gland for the oxidative stress genes CAT and SOD. As with the gills, MRP and P-gp were significantly downregulated at 100 μg L⁻¹ Cu²⁺, whereas Na/K ATPase was significantly upregulated at 50 and 100 μg L⁻¹ Cu²⁺. p53 and Rad51 showed significant downregulated at 100 μg L⁻¹ Cu²⁺. MT10 was significantly upregulated at both 50 and 100 μg L⁻¹ Cu²⁺. Significant activation of HSP70 and significant inhibition of ER1 were observed for exposure to 50 μg L⁻¹ Cu²⁺.

Benzo[a]pyrene exposure

mRNA levels in gills were not significantly changed except for one solvent control. In the digestive gland, significant changes in mRNA levels were observed in various genes from different modes of action. SOD and Rad51 showed significant upregulated at 15 µg L⁻¹ compared to the control. Significant activation was observed for MRP and P-gp at concentrations of 15, 30, and 60 μg L⁻¹. ER1 was also significantly upregulated at 30 and 60 μg L⁻¹ and MT10 exhibited significant downregulated at 15 and 30 µg L⁻¹.

Partial least squares-discriminant analysis (PLS-DA)

PLS-DA analysis was selected over PCA analysis due to it being a supervised method that reduces dimensionality and shows discrimination between variables with full awareness of treatments. PLS-DA can also be used for classification and selection of features. Additionally, PLS-DAs are a useful tool to aid in biomarkers and to make predictions using gene expression data (e.g. Nguyen and Rocke 2002; Christin et al. 2013).

PLS-DA models were run separately for expression data collected from the gills and digestive gland, repeated for Cu²⁺ and B[α]P exposures. The PLS-DA results were plotted on a factorial plane using correlation ellipses, representing the 95% confidence intervals of each treatment, and the genes used to generate the model were plotted in a correlation circle (Figures 1 and 2). VIP scores were calculated for all conditions (Table 3). Orientation of the vectors in Figure 2 describes the relationship between genes and how each gene describes the axes in Figure 1. The length of the vectors indicates the importance of the gene in describing the axes; the longer the arrow is from the centre the more important it is in describing the axis it is directed toward.

Figure 1. Partial Least Squares-Discriminant Analysis (PLS-DA) analysis for copper exposure at 25, 50, and 100 μg L⁻¹ in the A, digestive gland and B, the gills and for benzo[α]pyrene exposure at 15, 30, and 60 μg L⁻¹ in the C, digestive gland and the D, gills plotted on a factorial plane. Ellipses represent the 95% confidence intervals. Genes with a Variable Importance on Projection (VIP) score below 0.8 were removed from the model. Vectors showing the importance of each gene in determining the patterns shown are provided in Figure 2.

Figure 2. Partial Least Squares-Discriminant Analysis (PLS-DA) correlation circles for copper exposure in the A, digestive gland and the B, gills, and for benzo[α]pyrene exposure in the C, digestive gland and the D, gills plotted on a factorial plane. Orientation of the vectors describes the relationship between genes and how each gene describes the axes.

Table 3. Variable Importance on Projection (VIP) scores for all genes for each contaminant in gills and digestive gland, ordered from highest to lowest score.

Copper digestive gland		Copper gills		Benzo[α]pyrene digestive gland		Benzo[α]pyrene gills
Gene	VIP score	Gene	VIP score	Gene	VIP score	Gene	VIP scores
p53	1.37	SOD	1.50	P-gp	1.69	p53	1.41
MRP	1.29	p53	1.25	ER1	1.56	CAT	1.36
MT10	1.23	CAT	1.17	MRP	1.02	MRP	1.26
HSP70	1.14	MRP	1.11	MT10	0.98	SOD	1.08
Na/K ATPase	1.08	P-gp	1.08	Rad51	0.95	Rad51	0.99
Rad51	1.02	Na/K ATPase	1.04	CAT	0.92	P-gp	0.88
GST	0.92	Rad51	0.96	SOD	0.90	GST	0.87
P-gp	0.87	MT10	0.93	Na/K ATPase	0.69	HSP70	0.84
ER1	0.76	GST	0.62	HSP70	0.57	Na/K ATPase	0.83
CAT	0.45	ER1	0.30	p53	0.49	ER1	0.57
SOD	0.23	HSP70	0.25	GST	0.32	MT10	0.42

Higher values indicate the gene is more important in group discrimination. VIP values below 0.8 were removed from the PLS-DA analysis.

After copper exposure, gills and digestive glands analyses showed that almost all genes analyzed were important in describing variation between test conditions. p53, Rad51, MT10, MRP, P-gp, Na/K ATPase appear important in discriminating treatments in both gills and digestive gland (Table 3). SOD was the most important gene in the gills responding to Cu²⁺ exposure but the least important in the digestive gland. HSP70 and GST in the digestive gland were responsive to Cu²⁺ but not in the gills; conversely, CAT expression in the gills was responsive to Cu²⁺ exposure, but not in the digestive gland. However, ER1 was not crucial in discriminating Cu²⁺ concentrations in both gill and digestive gland tissues.

After B[α]P exposure, the expression of most genes in the gills discriminated between exposure concentrations, but in the digestive gland, fewer genes played an important role in discriminating treatments. The expressions of MRP, P-gp, CAT, SOD, and Rad51 were important in both the gills and digestive gland for discriminating between B[α]P concentrations; p53 was the most important gene in the gills for discriminating treatments but not important for the digestive gland whereas ER1 was the second most important in the digestive gland for discriminating treatments but not important in the gills.

Contaminant exposure concentrations

Verified concentrations of Cu²⁺ were averaged for replicate samples and are as follows: 22.7, 44.3, and 93 μg L⁻¹ Cu²⁺ (for nominal Cu²⁺ concentrations of 25, 50 and 100 μg L⁻¹, respectively). Verified concentrations of B[α]P were averaged between the replicates and are as follows: 13.85, 26.2, and 55.1 μg L⁻¹ B[α]P (for nominal B[α]P concentrations of 15, 30 and 60 μg L⁻¹, respectively). The mean difference between nominal and measured concentrations of B[α]P was less than 10% which is well within the degree of total experimental variability to be expected from an experiment of this type. Consequently, the nominal concentration of B[α]P is reported hereon in this paper.

The method detection limit for B[α]P, determined using a minimal signal-to-noise ratio of 3:1, was 0.01 μg L⁻¹. No detectable residues of B[α]P were measured in any of the quality assurance solvent blanks. Low concentrations of B[α]P were detected in two of the three replicate control and solvent control treatments on day zero at concentrations of 0.034 mg L⁻¹ (n = 2) and 0.044 mg L⁻¹ (n = 2), respectively. No detectable residues of B[α]P were measured in any of the triplicate control and solvent control treatments on day one.

The mean and 95% confidence intervals obtained for the recovery of the surrogate standard perylene-d₁₂ at the low (10 μg L⁻¹) and high (50 μg L⁻¹) spike rates were 99.2 ± 1.9% (n = 36) and 101.7 ± 1.4% (n = 12), respectively, and 99.8 ± 1.5% (n = 48) overall. The consistently high rate of recovery of the surrogate standard perylene-d₁₂ combined with the narrow 95% confidence intervals demonstrated the robustness of the sample extraction and analysis method.

Discussion

Contamination is present in New Zealand’s urban coastal ecosystems at concentrations similar to many other industrialised nations, but there is limited information on their effects on exposed biota (Stewart et al. 2014; Boehler et al. 2017). While biomonitoring has typically focused on apical endpoints (observable outcomes in a whole organism such as alterations to life-history traits and changes to population structure), it is advantageous to establish toxicogenomic approaches that can potentially detect effects at an early stage and facilitate the prediction of adverse outcomes at higher levels of biological organisation. This can allow for the development of tools for future risk assessments and guidelines for contaminants of concern. Establishing molecular tools requires: (1) the identification of candidate species that can be considered sentinel species of environmental health, and (2) molecular endpoints to investigate mechanistic pathways of toxicity. This study aimed to develop and validate P. canaliculus as a bioindicator species for ecotoxicological studies by developing and validating molecular biomarkers, using various gene targets historically used in ecotoxicology, targeting two tissue types, gills and digestive glands. The digestive gland and gills were selected in this study due to their essential roles in both detoxification of environmental contaminants and as a defensive barrier filtering suspended matter and pollutants. Within this study the reference chemicals were not measured in mussel tissues as the exposures were only for 48 h which we previously confirmed is adequate time to induce a biomarker response in mussel (Baettig et al. 2023). The development of biomarkers for use in field studies is more effective if dose–response relationships are adequately developed first under laboratory conditions, however future work should assess bioaccumulation potential of contaminants of emerging concern (Adams et al. 2001).

Copper exposure

Cu²⁺ is well established as a reference toxicant. Despite being well studied, data surrounding potential mechanisms of toxicity are not completely characterised in mussel species (Nguyen et al. 2018). The response of the oxidative stress markers, SOD and CAT, in P. canaliculus after Cu²⁺ exposure demonstrated similar results to the freshwater clam (Corbicula fluminea) where significant inhibition was observed in the gills following exposure to 10 and 50 μg L⁻¹ Cu²⁺ while a concentration of 50 μg L⁻¹ produced significant activation in the digestive gland (Bigot et al. 2011). HSP70 expression in the digestive gland, a biomarker of generic stress and thought to precede Cu²⁺ induced oxidative damage, was significantly upregulated in the present study, which is consistent with previous studies on M. galloprovincialis, M. edulis, and M. coruscus (Chatel et al. 2018; Xu et al. 2018; Franzellitti et al. 2020). MT10, a metallothionein isoform, is a metal-binding protein well documented in the literature to sequester metals and play a role in cellular antioxidative defense (Monserrat et al. 2007; Vergani et al. 2007). MT10 was significantly upregulated in this study in the digestive gland at 50 and 100 μg L⁻¹ Cu²⁺. Overall, these results suggest P. canaliculus response to Cu²⁺ exposure measured in classical biomarkers is comparable to previously evaluated mussel species.

Membrane transporters, particularly ATP-binding cassette (ABC) transporters, play an important role in defending aquatic organisms as a first line of cellular defense (Smital et al. 2004; Luckenbach and Epel 2008). In our study, the expression of P-gp and MRP were significantly downregulated in both the gills and digestive gland. These ABC transporters have been shown to respond to Cu²⁺ exposure in bivalves species, Crassostrea angulate, M. edulis, and M. galloprovincialis (Shi et al. 2015; Chatel et al. 2018; Franzellitti et al. 2020), suggesting that P. canaliculus may display equivalent responses to other mollusk species. Na/K ATPase, a gene crucial for translocating ions across membranes and maintaining ion homeostasis within the cell, was significantly induced in the digestive gland but was significantly downregulated in the gills following exposure to Cu²⁺. Cu²⁺ is well documented to negatively impact Na⁺ homeostasis and inhibit Na/K ATPase activity suggesting activation in the digestive gland may be a compensatory mechanism to counteract the disruption of ion balance by Cu²⁺ before the cell becomes overwhelmed (Monserrat et al. 2007). Freshwater mussels (Lampsilis siliquoidea) chronically exposed to Cu²⁺ (28 d) and blue mussels (M. edulis) exposed for 24 h showed decreased Na/K ATPase activity and Na/K ATPase expression in the gills (Jorge et al. 2013; Chatel et al. 2018), respectively, which is consistent with the findings of our study.

We observed inhibitions of Rad51 and p53 in both the gills and digestive gland. Tumour suppressor protein p53 plays an important role in cell-cycle arrest, apoptosis, the repair of genotoxic damage, and cell survival and regulation of oxidative stress (Monteiro et al. 2009). Rad51 protein facilitates homologous recombination repair of DNA double-strand breaks, protects genome integrity, and helps cells avoid apoptosis (Monteiro et al. 2009). One potential explanation for these results is that DNA damage has occurred that could not be repaired, leading to a cell death. This would need to be verified with additional sublethal endpoints to demonstrate DNA damage, but the observation in the present study is consistent with previous research (Monteiro et al. 2009; Chatel et al. 2018; Vernon and Jha 2019).

Gills are the first tissue exposed to environmental contaminants, making them susceptible to modulation of excretion pathways and DNA response and repair pathways (Smital et al. 2004; Luckenbach and Epel 2008). The significant inhibitions of MRP, P-gp, p53, and Rad51 observed, suggest an overwhelming of the gills. In the digestive gland, significant activation of the metal detoxifying gene MT10, membrane transporter Na/K ATPase, and the generic stress marker HSP70 was observed in response to exposure to Cu²⁺. In contrast, inhibition of DNA response and repair genes and ABC transporters was observed. Overall, our results suggest P. canaliculus response to Cu²⁺ exposure is comparable to those reported for other mussel species. Still, we cannot conclude categorically due to the limited number of studies that have measured gene expression in mussels exposed to Cu²⁺ (Chatel et al. 2018).

It can be challenging to establish meaningful conclusions on ecotoxicity using only pairwise analysis of gene expression therefore, integrative modelling to identify mechanistic trends can be useful to better understand how organisms respond to stressors (Jaumot et al. 2015). Multivariate analyses such as PLS-DA provide a holistic assessment of how the genes respond and can be used to discriminate between treatments as well as identify relationships between genes (Jaumot et al. 2015). After Cu²⁺ exposure, test conditions were separated in both tissue types, suggesting that response profiles differ. However, the overlap observed in the digestive gland of mussels exposed to the control and the lowest Cu²⁺ concentration suggested minimal differences between the two treatments. When looking at the relationship between genes, we see that DNA repair genes and ABC transporters all responded consistently, suggesting that exposure to Cu²⁺ could potentially affect membrane transporters and genotoxicity in P. canaliculus.

Benzo[α]pyrene exposure

Polycyclic aromatic hydrocarbons (PAHs) are one of the most widespread and persistent organic contaminants in the environment, and B[α]P is widely recognised for its genotoxicity and carcinogenicity (Banni et al. 2017). B[α]P produces reprotoxic and epigenetic effects at low concentrations, making it a contaminant of environmental concern (Guo et al. 2020). B[α]P has been demonstrated in marine bivalves to induce genotoxicity, oxidative stress, and endocrine disruption (Banni et al. 2017).

Phase I and II enzymes and oxidative stress markers play an important role in detoxifying B[α]P (Chatel et al. 2012). Previous research on the freshwater zebra mussel (Dreissena polymorpha) exposed to 10 μg L⁻¹ demonstrated oxidative stress and impact on membrane transporters within 12 h of exposure to B[α]P but no response was observed at 100 μg L⁻¹ (Chatel et al. 2012). This suggests that P. canaliculus has a similar response profile to B[α]P, further validating it as a potential sentinel species for New Zealand.

The most notable response observed following exposure to B[α]P was in the ABC transport genes, which were significantly upregulated at all concentrations in the digestive gland. Previous studies have shown that PAHs can be transported out of cells by multixenobiotic resistance-mediated processes and upregulation of P-gp in multiple taxa has been demonstrated when exposed to B[α]P (Chatel et al. 2012; Ferreira et al. 2014; Guo et al. 2020) at concentrations ranging from 2 to 280 μg L⁻¹.

Multivariate analysis of B[α]P exposure in the digestive gland of P. canaliculus demonstrated that DMSO solvent control overlapped with the control, demonstrating no impact of the solvent control on gene response. However, in the multivariate analysis of B[α]P exposure in the gills a clear separation between the control and the DMSO solvent control can be seen, indicating that DMSO may be the dominant descriptor in genotoxicity. Additionally, a significant effect was observed for gills exposed to the DMSO solvent control. This was unexpected and suggests the use of DMSO in the preparation of the B[α]P treatment doses can potentially mask the biological response to a contaminant. The final concentration of 0.01% (v/v) DMSO used in our experiments is recommended by the OECD (2019), but these test guidelines were not designed to assess responses at lower levels of biological organisation, such as gene expression. The interactions of solvents are often overlooked in ecotoxicity studies. Our results suggest that further investigations into their effects are required to clarify their contributions and determine the appropriate concentrations for gene expression studies to limit this potential confounding factor. However, the use of the biomarkers remains valid in a field situation where no solvents are present.

The results of our study suggest P. canaliculus response to B[α]P was limited compared to what has been reported in previous bivalve studies. This could potentially be explained by different exposure times or by the noise from the solvent, which should be further investigated. Regardless, only a few studies have investigated the toxicity of B[α]P through gene expression in mussels. While comparisons between individual genes may not have led to statistically significant results, the strength of molecular endpoints in ecotoxicity is to better understand mechanisms of toxicity of targeted contaminants following integrative, multimarker approaches, as well as to highlight how these effects can lead to adverse outcomes at higher levels of biological organisation.

Conclusion

After exposure to two reference contaminants, modulation of genes in P. canaliculus suggested similarities in response profiles to other bivalves, but differences were also evident. Exposure duration plays a role in biological response variations and without standardised guidelines it can be difficult selecting an optimal sampling time. Future studies should investigate chronological toxicity to establish standardised guidelines. The present study developed novel RT-qPCR primers, advancing knowledge on molecular tools for P. canaliculus and elucidates potential modes of action for legacy contaminants. Molecular tools can provide robust insights into mechanistic response to anthropogenic activity, but it is important acknowledge that some of the target biomarkers allow for a characterisation of generic stress whereas other primers may be more specific to the target contaminant. While molecular biomarkers can be useful tools, the present study highlights the need for integrative analysis to provide a holistic assessment of the organismal response to contamination. Linking the molecular biomarkers developed in the present study to additional endpoints such as protein-based assays, energetic reserves, and other endpoints at higher levels of biological organisation can improve understanding of the predictive capabilities of molecular endpoints and help establish adverse outcome pathways (Ankley et al. 2010). Incorporating multiple endpoints at different levels of biological organisation are necessary to fully elucidate biological response profiles and develop predictive ecotoxicological platforms. The present study aimed to develop molecular endpoints for P. canaliculus for use in ecotoxicity studies using a multimarker approach. However, their use is not limited to ecotoxicity studies as they can be applied to address other physiological investigations.

Acknowledgements

We would like to thank Sanford Ltd. for collecting and providing us with mussels from their aquaculture farms.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research was funded through the New Zealand Ministry of Business Innovation and Employment (MBIE) Endeavour grant: Managing the Risk of Emerging Contaminants (CAWX1708).