Upwelling-level acidification and pH/pCO2 variability moderate effects of ocean acidification on brain gene expression in the temperate surfperch, Embiotoca jacksoni.

Acidification-induced changes in neurological function have been documented in several tropical marine fishes. Here, we investigate whether similar patterns of neurological impacts are observed in a temperate Pacific fish that naturally experiences regular and often large shifts in environmental pH/pCO2 . In two laboratory experiments, we tested the effect of acidification, as well as pH/pCO2 variability, on gene expression in the brain tissue of a common temperate kelp forest/estuarine fish, Embiotoca jacksoni. Experiment 1 employed static pH treatments (target pH = 7.85/7.30), while Experiment 2 incorporated two variable treatments that oscillated around corresponding static treatments with the same mean (target pH = 7.85/7.70) in an eight-day cycle (amplitude ± 0.15). We found that patterns of global gene expression differed across pH level treatments. Additionally, we identified differential expression of specific genes and enrichment of specific gene sets (GSEA) in comparisons of static pH treatments and in comparisons of static and variable pH treatments of the same mean pH. Importantly, we found that pH/pCO2 variability decreased the number of differentially expressed genes detected between high and low pH treatments, and that interindividual variability in gene expression was greater in variable treatments than static treatments. These results provide important confirmation of neurological impacts of acidification in a temperate fish species and, critically, that natural environmental variability may mediate the impacts of ocean acidification.

INTRODUCTION

Ocean acidification (OA; here defined as both increased ocean pCO2 and decreased pH) has been identified as a major threat to marine species (Kroeker et al., 2013). Several studies have documented changes in neurological functioning, including altered cognition, sensory function and behaviour, in marine fish (e.g., Domenici et al., 2012; Hamilton et al., 2013; Hamilton et al., 2017; Munday et al., 2010; Pistevos et al., 2015), raising concerns about neurological impacts leading to changes in the strength of species interactions (e.g., predation). In contrast, more recent work has questioned the generality and replicability of such impacts across studies and species (Clark et al., 2020a). Additionally, much of the evidence of neurological impacts comes from studies of a few tropical reef species under static pH/pCO2 regimes (Nagelkerken & Munday, 2016). Physiological evidence indicates these neurological impacts may be the result of a hypercapnia‐driven reversal of electrochemical gradients in GABAergic neurons. This has been hypothesized to result from internal acid–base balance processes that lead to an accumulation of intracellular [HCO3 −] and/or a decrease in extracellular [Cl−] (Heuer & Grosell, 2014; Nilsson et al., 2012). This shift in ion concentrations is thought to cause neuron depolarization upon GABAA receptor activation rather than the hyperpolarization expected under nonacidified conditions, reversing the functional nature of these neurons from inhibitory to excitatory and presumably causing the observed shifts in cognition and behaviour (Heuer & Grosell, 2014; Nilsson et al., 2012; Schunter et al., 2019). Given this body of evidence, altered neurological function may be a major pathway through which changing seawater carbonate chemistry will impact fitness in marine fish. Continued work elucidating the molecular mechanisms underlying these changes is therefore critical for moving the field forward.

Many coastal ecosystems experience significant environmental variability over a range of temporal scales, including fluctuations in seawater pH/pCO2 (Chan et al., 2017; Hofmann et al., 2011; Kang et al., 2022; Kroeker et al., 2020). In upwelling regions, where deeper, more acidic water is brought to the ocean surface, pH can vary by half a unit over a period of weeks (Hirsh et al., 2020; Hofmann et al., 2011). In seagrass beds, pH can vary by a whole unit over a period of hours to weeks due to fluctuations in photosynthesis and respiration and tidal movement (Duarte et al., 2013; Hofmann et al., 2011). These fluctuations often reach or exceed predictions for the mean future ocean pH under OA (Gruber et al., 2012; Hauri et al., 2013; Takeshita et al., 2015). In previous studies, exposure to low pH/high pCO2 seawater has affected indicators of fish neurological function anywhere from 2–12 days after exposure has ceased (Hamilton et al., 2013; Munday et al., 2010), but it is unclear what duration of exposure elicits these effects. How temporal environmental variability moderates fish responses to low pH/high pCO2 remains critically understudied, leaving many unanswered questions about how OA may realistically affect populations in nature (but see Jarrold et al., 2017; Jarrold & Munday, 2019; Schunter et al., 2021). Investigations into the effects of pH/pCO2 variability are important for accurate prediction of the severity of impacts acidification will have on natural populations and ecosystems (Kroeker et al., 2020). For example, if variability dampens the effects documented in studies using only static pH/pCO2 treatments (as seen in Jarrold et al., 2017 and Jarrold & Munday, 2019, where diel pCO2 fluctuations ameliorated impairments in behaviour and growth seen under static decreases in pCO2), we may be overestimating the effects of OA, and overlooking an important role that variability may play as a provider of temporal refuge. Conversely, if variability exacerbates negative impacts of acidification, acting as an additional stressor, we may be underestimating the potential impact of acidification on natural populations.

We expect physiological responses to OA to be reflected in the gene expression of the affected organism (Griffiths et al., 2019; Hamilton et al., 2017). In particular, we expect changes in brain gene expression to be associated with shifts in neurological and cognitive function (Schunter et al., 2016). Given the proposed mechanism of OA‐induced cognitive impairment described above, we expect experimental acidification to impact expression in genes related to the maintenance of homeostasis and neuronal signalling, such as ion transporter and signal receptor genes, and those involved in the GABAergic signalling pathway. Changes in expression in these gene categories have been noted in spiny damselfish (Schunter et al., 2018) and three‐spined stickleback (Lai et al., 2016), but this has not yet been investigated in a temperate reef species with an evolutionary history of exposure to fluctuating pCO2.

Here, we present two experimental studies of the effects of acidification on brain gene expression in a common temperate reef fish, the black surfperch (Embiotoca jacksoni). Surfperches make up a large proportion of fish biomass on California rocky reefs (Laur & Ebeling, 1983) and support an immensely popular recreational fishery. E. jacksoni is found in both upwelling reef systems and estuarine seagrass ecosystems, and therefore has an evolutionary history of exposure to variable pH conditions that are often more extreme than those experienced by tropical reef fish (Hofmann et al., 2011). A few studies have investigated the effects of acidification on temperate reef fishes (e.g., Cline et al., 2020; Hamilton et al., 2017; Kwan et al., 2017), but surfperches are unique because they exhibit viviparity and no pelagic larval phase, with young born as developed juveniles. Additionally, E. jacksoni has limited adult dispersal (Bernardi, 2000, 2005; Hixon, 1981). These two life‐history traits increase the likelihood of adaptation to local environmental conditions in E. jacksoni, which may lead to divergent effects of acidification in this species compared to other temperate fish (e.g., greater physiological adaptation to OA in populations that have historically experienced local acidification).

Experiment 1 was designed to determine the presence and extent of any impacts of acidification on E. jacksoni brain gene expression and used a static and more extreme acidified treatment (pH ~7.30). In Experiment 2, we used a less extreme static treatment and incorporated two variable treatments with different mean pH levels to mimic upwelling‐scale pH variability. This experiment was designed to test the potential role of temporal pH/pCO2 variability in mediating any neurological effects of acidification. Together, the results of these experiments provide a more comprehensive understanding of the impacts of acidification on marine organisms, particularly in dynamic, temperate ecosystems.

MATERIALS AND METHODS

Collections and acclimation

We collected young‐of‐the‐year E. jacksoni from Elkhorn Slough (Monterey County, CA) using a beach seine. Collected fish were placed in coolers and driven back to UCSC‐CSC, where they were kept in outdoor flow‐through containers until the start of each experiment. For exact dates of collections, acclimation periods, and experimental manipulations see Table S1.

Experimental design

We conducted two separate experiments with similar methods in November 2015 and September 2017 at University of California, Santa Cruz's Coastal Science Campus (UCSC‐CSC). Both experiments treated E. jacksoni juveniles in outdoor flow‐through seawater systems. In Experiment 1 (2015), we set target treatments at pH 7.85, representing a common current upwelling condition along the coast of Central California, and pH 7.30, representing a current extreme estuarine event or future extreme upwelling event (Chan et al., 2017; Hofmann et al., 2011; Lowe et al., 2019; Takeshita et al., 2015). Both treatments in this experiment held the target pH constant (static) over the course of the experiment. Five randomly assigned juvenile E. jacksoni were distributed across two replicate tanks at pH 7.30 and four were distributed across two replicate tanks at pH 7.85 (opaque 200 L plastic drums). Seawater pH treatments were replicated only at the level of holding tanks. Replicates were then brought down to experimental pH levels over 7 days. Tissue sampling was conducted after 23 days of treatment (Table S1).

In Experiment 2 (2017), we incorporated upwelling‐scale pH variability into two of the treatments, and target pH levels were set at more conservative levels. Two static pH treatments were set at target pH levels of 7.85 and 7.70, approximating conservative present and future reef conditions during the upwelling season (Chan et al., 2017; Takeshita et al., 2015). For each static treatment, there was a corresponding variable treatment that oscillated around the same mean pH as the static treatment with an amplitude of ±0.15 pH and a period of 8 days (Figure 1), approximating a typical upwelling pattern (Hofmann et al., 2011). An additional treatment, hereafter referred to as “ambient”, had a static target pH of 8.00. However, because our pH control system was not capable of increasing pH above that of the incoming seawater, periodic natural decreases in the pH of the input seawater below 8.00 caused this treatment to exhibit an intermediate level of variability between that of the static and variable treatments (Figure S1). We used 10 header buckets (two per treatment) to create the five different pH treatments, with two replicate tanks per header. We randomly assigned six juvenile E. jacksoni to replicate tanks (translucent 61 L plastic containers; 6 individuals × 4 replicate tanks × 5 pH treatments) and allowed them to acclimate at ambient pH for 2–3 days. We then allowed the pH of each treatment to slowly approach its starting pH (target pH 8.00, 7.85, or 7.70) over a period of 2 days. After an additional 4–5 days, the variable treatments began their programmed oscillations (Figure S1). Due to logistical restrictions, the treatments were separated into two groups (pH 7.85 and ambient treatments, pH 7.70 treatments) that were staggered in their timing by 1 day (Figure S1). Using a custom‐built LabView program, set points for the variable treatments were changed throughout the experiment at intervals of 0.003125 pH/h to create 8‐day cycles. During this experiment, fish were removed from their treatment tanks on two occasions to conduct behavioural assays, after which they were returned to their treatment tanks. Because we were met with logistical challenges that precluded the proper execution of these trials, these data were not analysed. Eight days (1 full cycle of the variable treatments) were allowed to elapse between the last trial and tissue sampling, which was conducted after 22 days (Figure S1; Table S1).

FIGURE 1

Experiment design and data analysis pipeline for Experiments 1 and 2.

pH control system and sampling of seawater chemistry

Seawater pH was manipulated using a custom‐built feedback control system. Two large sumps received a continuous flow of ambient seawater. One of these sumps (“low pH”) was continuously bubbled with CO2 gas, while the other (“ambient”) was left untreated. Lines from both sumps fed seawater into header buckets at varying rates to create pH treatments. The pH of each bucket was continuously measured by Honeywell Durafet II sensors connected to Honeywell Universal Dual Analysers (UDAs; see Kapsenberg et al., 2017). Seawater pH in each header was controlled through a feedback system, where a solenoid valve determined the flow of low pH seawater to the header to either increase or decrease pH. The mixed treatment water from each header then flowed out into two replicate holding tanks. We oxygenated and mixed seawater in each header using air pumps/stones and/or water pumps (Experiment 2 only).

Prior to beginning each experiment, we calibrated the Durafet sensors from the header buckets using equimolar Tris buffer (DelValls & Dickson, 1998) obtained from the Dickson Laboratory (Scripps Institution of Oceanography). In Experiment 1, discrete water samples were taken from the replicate tanks at seven time points and used for characterization of carbonate chemistry via spectrophotometric pH analysis and open cell total alkalinity titration (Dickson et al., 2007). In Experiment 2, samples were taken from headers at five time points and used for post hoc calibration of Durafet pH measurements. Using a handheld sensor (YSI), we also measured temperature, pH, dissolved oxygen, and salinity in each replicate tank daily and, in Experiment 2, in each header as well to allow for calibration of YSI pH measurements to calibrated Durafet measurements. See Tables 1 and 2 for measured and calculated seawater parameters.

TABLE 1

Carbonate chemistry and environmental parameters for treatment containers in Experiment 1

Treatment	pHT (spec)	pCO2 (μatm)	Ω	TA (μmol/kg)	Temp (°C)	Salinity (ppt)	pHT (YSI)
Target pH 7.85	7.88 ± 0.02	599 ± 36	1.50 ± 0.08	2193 ± 59	13.0 ± 0.5	33.8 ± 0.1	7.89 ± 0.04
Target pH 7.30	7.35 ± 0.06	2204 ± 333	0.48 ± 0.07	2212 ± 20	12.2 ± 0.8	33.8 ± 0.1	7.36 ± 0.14

Note: Aragonite saturation state (Ω) and pCO2 were calculated with the R package seacarb (Gattuso et al., 2021) using the spectrophotometric pH and total alkalinity (TA) values from discrete bottle samples, and salinity and temperature values from YSI readings. All values are means ± SD. Mean pHT (spec) and TA were calculated from bottle samples taken at seven time points across the experiment. Mean pHT (YSI) was calculated from daily readings that were calibrated using the discrete bottle samples.

TABLE 2

Carbonate chemistry and environmental parameters for the headers of each treatment in Experiment 2

Treatment	pHT (Durafet) (hourly)	pCO2 (μatm)	Ω	TA (μmol/kg) (bottle samples)	Temp (°C) (Durafet)	Salinity (ppt) (YSI)	pHT (YSI) (daily)
Ambient	8.00 ± 0.04	452 ± 55	2.36 ± 0.23	2266 ± 3	17.5 ± 0.06	34.3 ± 0.1	8.01 ± 0.08
Target pH 7.85 ‐ Static	7.90 ± 0.01	586 ± 14	1.94 ± 0.07	2268 ± 3	17.6 ± 0.06	34.3 ± 0.1	7.90 ± 0.05
Target pH 7.85 ‐ Variable	7.89 ± 0.08	614 ± 136	1.93 ± 0.34	2268 ± 4	17.5 ± 0.06	34.3 ± 0.1	7.88 ± 0.11
Target pH 7.70 ‐ Static	7.76 ± 0.04	848 ± 64	1.46 ± 0.17	2268 ± 5	17.6 ± 0.06	34.3 ± 0.1	7.75 ± 0.08
Target pH 7.70 ‐ Variable	7.76 ± 0.09	870 ± 199	1.47 ± 0.30	2267 ± 4	17.6 ± 0.06	34.3 ± 0.1	7.74 ± 0.10

Note: Aragonite saturation state (Ω) and pCO2 were calculated with the R package seacarb (Gattuso et al., 2021) using the Durafet pH and temperature values, average TA values from discrete bottle samples, and salinity values from YSI readings. All values are means ± SD. Mean pHT (Durafet) values were calculated using hourly averaged pH readings (from headers) that were calibrated using discrete (bottle) water samples and include only the time period of the first two pH cycles (1/2 September–17/18 September). Mean pHT (YSI) values were calculated using daily readings (from replicate containers) that were calibrated using bottle‐calibrated Durafet values (taken simultaneously) from the 1/2 September–17/18 September date range and include YSI readings from the entire length of the experiment (1/2 September–23/24 September).

Experimental considerations

A heat wave struck Santa Cruz during Experiment 2. This added stress may have contributed to the juvenile mortality observed during this experiment (48 out of 120 fish died, unrelated to treatment; ANOVA, F = 2.088, p = .133). Additionally, some of the Durafet sensors (4 of 10) used to control the pH in the headers experienced heavy fouling by microalgae toward the end of Experiment 2. This probably led to artificially high pH measurements for about 2 h around midday due to photosynthesis, and thus a corresponding over‐correction by the pH control system. Because of this issue, the pH of certain tanks was probably lower during midday than their respective set points. To better understand the scale of this over‐correction, we conducted a test 11 days after tissue sampling, in which all Durafets were placed in the same header with no active pH control (Figure S2). Though the effect of fouling on recorded pH appeared to strengthen over the 11 days since tissue sampling, this post‐hoc test revealed variability in the impact across headers, and a relatively even distribution of fouling across treatments (Figure S2). The greatest spike during this test occurred in one of the two ambient headers with a magnitude of ~0.5 pH units, but this treatment was not included in most of our analyses, and thus we believe it does not affect our conclusions. Examination of experiment Durafet readings from the ambient header (which, due to its high set point and limited pH control, displayed the true extent of the pH spikes) revealed that significant spikes (deviation of ~0.05 pH units or greater) in this most affected treatment only began occurring approximately 3 days before tissue dissection. Because our test indicated that the other headers were affected to a much smaller degree (Figure S2), the other treatments probably did not experience midday spikes of >0.05 pH units for any significant duration prior to the end of the experiment. To prevent the inclusion of spurious pH data points in the characterization of the experimental treatments, we used the continuous Durafet pH and temperature data for the dates 1–17 September, after which we used pH and temperature data from daily YSI readings taken from each header. This shift in sampling frequency probably explains much of the apparent increased variability of the pH treatments after 17 September (Figure S1), as the YSI data represents only a daily snapshot of the pH of each header. The pH sensor within the YSI is also functionally different from those within Durafets (Martz et al., 2010). Finally, outside of the daily spikes, the Durafet pH data collected 17–24 September showed no obvious departure from the precision seen earlier in the experiment. We therefore contend that, apart from the midday overcorrections experienced at the end of the experiment, the true precision and variability of the pH treatments was unchanged after 17 September, and any apparent changes reflect only a change in the pH‐sensing instrument used.

Fish care and handling

This experiment was run under the approval of UCSC IACUC project proposals BERNG1312 and KROEK1503_A2. We performed system checks at least daily and fed fish frozen shrimp every day (Experiment 1) or a mix of frozen brine shrimp, Spirulina brine shrimp, and mysis shrimp every other day (Experiment 2). Tanks were cleaned and excess food removed approximately seven times per week (Experiment 1) and at least once per week (Experiment 2). To minimize stress, a shelter was placed in each replicate. To reduce heat and sun exposure, shade cloth was kept over the top of the replicate tanks whenever water monitoring, cleaning, or feeding was not occurring.

Tissue sampling

At the end of each experiment, we dissected tissue from four individuals from each treatment. Individuals were dissected one at a time, with all dissections taking less than 10 min from the time of fish removal from its tank. Brain and lateral muscle tissue were dissected and sequenced in Experiment 1 for use in the transcriptome assembly, but only brain tissue was sequenced in Experiment 2. In Experiment 2, the whole brain was dissected at the approximate time when the variable pH treatments were crossing (in the ascending direction) their target mean pH levels (Figure S1). Only brain gene expression analysis will be further discussed here. Tissue was stored in screw‐cap tubes and flash frozen in liquid nitrogen. We stored all tissue samples at −80°C until RNA extraction.

RNA extraction and library preparation

Dissected whole brains were arbitrarily subsampled and homogenized using a Qiagen TissueLyser. A discussion of the potential effects of subsampling are included below in the “Interindividual variability in gene expression” section of the Discussion. We extracted RNA using the Qiagen RNeasy Mini extraction kit. RNA quality and quantity were assessed using a NanoDrop spectrophotometer and Qubit fluorometer. RNA was stored in DEPC‐treated water at −80°C. cDNA libraries were prepared from 1 μg of total RNA using the New England Biolabs NEBNext Ultra II RNA Library Prep kit. Prepared libraries were sequenced on an Illumina HiSeq 4000 (150 bp SE) at the QB3 Vincent J. Coates Genomics Sequencing Laboratory at the University of California, Berkeley.

Read processing and transcriptome assembly

We removed adapters and trimmed/removed low quality reads using the Trimmomatic software (v0.36; Bolger et al., 2014; parameters = LEADING:2 TRAILING:2 SLIDINGWINDOW:4:2 MINLEN:25) and quality checked the trimmed sequences using FastQC (version 0.11.7; Andrews, 2010). We used the trimmed reads from all sequenced samples from both experiments to assemble a brain/muscle tissue combined transcriptome for E. jacksoni using the genome‐guided TopHat/Cuffmerge/Cufflinks pipeline (default parameters; TopHat version 2.1.1, Cufflinks version 2.2.1; Trapnell et al., 2012). This pipeline creates separate assemblies for each sample, which are then merged. A draft, scaffold‐level E. jacksoni genome assembly was used as the reference (see Supporting Information Materials). We annotated the transcriptome assembly by running a blastx query (e‐value cutoff = 1e‐3; NCBI, Altschul et al., 1990) against the SwissProt database (uniprot_sprot.dat.gz downloaded 25 April 2020; The UniProt Consortium, 2021).

Multivariate analysis of gene expression

In Experiment 1, we sequenced transcripts from both muscle and brain tissue, but only brain gene expression will be discussed here. To characterize global gene expression of individuals, trimmed reads were aligned and quantified into gene‐level expression data using bowtie (version 1.2.3; Langmead et al., 2009) and RSEM (version 1.3.3; Li & Dewey, 2011) within the Trinity software package (version 2.9.1; Haas et al., 2013). Raw read counts were then filtered to remove genes with low expression using the default parameters of the filterByExpr function (min.count = 10, min.total.count = 15, large.n = 10, min.prop = 0.7) in the R package, edgeR (version 3.34.0; R Core Team, 2021; Robinson et al., 2010). We normalized the read counts using the TPM method, as implemented by the calcNormFactors function in the edgeR package, then log2‐transformed the data using the cpm function (prior.count = 2). The transformed data were dimensionally reduced through multidimensional scaling (metric MDS in Experiment 1, nMDS in Experiment 2) using Manhattan distances, as implemented through the wcmdscale and metaMDS functions in the vegan package for R (version 2.5.7; Oksanen et al., 2020).

To test whether global gene expression profiles differed among treatments, we ran a permutational multivariate analysis of variance (PERMANOVA; Anderson, 2001, 2017) on the transformed expression data using the adonis function of the vegan package (method = “manhattan”, perm = 1,000,000). In Experiment 1, the sole model factor was pH level (7.85, 7.30). In Experiment 2, the model was run with two factors: pH level (7.85, 7.70) and pH variability (static, variable). For Experiment 2, pairwise comparisons between treatments were conducted using the pairwise.adonis function of the pairwiseAdonis package (sim.method = “manhattan”, perm = 1,000,000) (Martinez Arbizu, 2020).

Differential gene expression analysis

Using the gene‐level counts matrix created by RSEM, we identified differentially expressed genes (DEGs) between all pairwise treatment comparisons using the edgeR package, as implemented through Trinity. To buffer against false positives and noise due to the experimental conditions described above, we used a conservative FDR cutoff value of 0.001 (‐P parameter) and a fold‐change cutoff of 1.5 (‐C parameter) to create the final list of DEGs for each treatment comparison. We then repeated MDS procedures and PERMANOVA tests as described above, using only these DEGs.

Functional enrichment analysis

To identify gene sets (groups of functionally related genes) that were significantly enriched in a given treatment comparison (e.g., pH 7.85 vs. pH 7.30) we used the threshold‐free analytical method, gene set enrichment analysis (GSEA; Subramanian & Tamayo et al., 2005) as implemented through the FGSEA package for R (Korotkevich et al., 2021). Given a ranked list of genes derived from differential expression analysis, this method yields a list of gene sets from user‐supplied gene set databases ‐ in this case gene ontology (GO; The Gene Ontology Consortium, 2020), KEGG (Kanehisa & Goto, 2000), and the MSigDB Hallmark collection (Liberzon et al., 2015) ‐ that are enriched among up‐ and downregulated genes. Enrichment analysis was completed for Experiments 1 and 2 separately, and the resulting enriched gene sets from analogous treatment comparisons were then contrasted across experiments to identify commonly enriched gene sets.

RESULTS

Seawater pH manipulation

Mean seawater pH levels were maintained near their target set points in each experiment. Seawater parameters and information on how they were calculated are given in Tables 1 and 2. Note that for consistency, we continue to use target pH levels to refer to each treatment.

Sequencing and transcriptome assembly

The TopHat/Cufflinks pipeline yielded a transcriptome made up of 71,933 assembled transcripts grouped into 39,258 putative genes. Aligning the trimmed reads back to the assembled transcriptome resulted in an 82.98% alignment rate. A blastx search of the transcriptome against the SwissProt database revealed that 8836 transcripts represented nearly full‐length transcripts (>80% alignment coverage), and that 16,574 proteins were represented in the transcriptome at some level of alignment coverage. We used BUSCO (version 4.0.6; Simão et al., 2015) to quantify the completeness of our transcriptome and found that of the 3640 BUSCO orthologs in the Actinopterygii data set, 76.9% were found complete (41.2% single‐copy, 35.7% duplicated), 6.5% were found fragmented, and 16.6% were missing. For more statistics on the assembly, see Table S2.

Total sequenced reads per sample are listed in Tables S3 and S4. Of the 39,258 putative genes in the assembled transcriptome, 22,961 (Experiment 1) and 33,597 (Experiment 2) remained in each dataset after filtering for genes with low expression.

Multivariate analyses of global gene expression

Single‐factor PERMANOVA analysis identified a strong and significant effect of pH level on global gene expression in both Experiment 1 (single‐factor; r 2 = 0.811, F = 25.766, p = .029) and Experiment 2 (two‐factor, ambient excluded; r 2 = 0.159, F = 2.53, p = .021), with pH explaining 81 and 16% of the observed variation, respectively (Tables S5 and S6). In Experiment 2, we did not detect an effect of pH variability on global gene expression (r 2 = 0.037, F = 0.59, p = .890) or an interaction of pH level and variability (r 2 = 0.052, F = 0.84, p = .524). Pairwise comparisons of all treatments (including ambient) revealed two comparisons were nearly significantly different (Table S7): ambient versus 7.70 static (r 2 = 0.230, F = 1.79, p = .086) and 7.85 static versus 7.70 static (r 2 = 0.292, F = 2.47, p = .057). We visualized the differences in global gene expression patterns between treatments using MDS (Figures S3 and S4).

We found 10,656 DEGs between the treatments in Experiment 1 (Figure 2). In Experiment 2, we found a total of 200 DEGs across all treatment comparisons (Table 3; Figure S5). The 7.85 static versus 7.70 static comparison produced the majority of DEGs (159) in this experiment. The 7.85 static versus 7.70 variable and 7.85 variable versus 7.70 static comparisons produced 11 DEGs each and six genes each were differentially expressed in the static versus variable comparisons of both the 7.85 and 7.70 pH levels (one gene was consistently differentially expressed across the two comparisons).

FIGURE 2

Heatmap of gene expression profiles for each individual in Experiment 1. Each column represents an individual fish, and each row represents a differentially expressed gene. Yellow colours represent upregulation in a given treatment and purple colours represent downregulation. Brighter hues represent larger differences in relative gene expression across the treatments.

TABLE 3

Number of DEGs detected across all treatment comparisons in Experiment 2

	Ambient	pH 7.70 static	pH 7.70 variable	pH 7.85 static
Ambient
pH 7.70 static	5
pH 7.70 variable	8	6
pH 7.85 static	3	159	11
pH 7.85 variable	6	11	9	6

Since we did not expect transcriptome‐wide shifts in gene expression across pH variability treatments in Experiment 2 (Figure S4), and were instead interested in how the expression of acidification response genes was affected by the introduction of environmental variability, we repeated our multivariate analyses of gene expression for only the DEG subset of Experiment 2. For consistency, analogous analyses were also performed for the Experiment 1 DEG subset (Table S8; Figure S6).

DEG expression differed among pH levels in Experiment 2 (r 2 = 0.388, F = 10.77, p = .004), with pH explaining 39% of the observed variation (Figure 3). We did not detect an effect of pH variability (r 2 = 0.041, F = 1.15 p = .291). The interaction of pH level and variability, however, was marginally significant (r 2 = 0.139, F = 3.85, p = .052; Table S9). Pairwise comparisons of all treatments revealed a significant difference between the pH 7.85 and 7.70 static treatments (p = .029) and nearly‐significant differences between the ambient and 7.85 static treatments (p = .057) and between the ambient and 7.70 static treatments (p = .086; Table S10). The comparison of 7.85 static versus 7.85 variable was also nearly significant (p = .057), but the comparison of 7.70 static versus 7.70 variable was less so (p = .171).

FIGURE 3

nMDS plot of DEG expression in Experiment 2. Points represent single individuals. Ellipses are 95% confidence ellipses.

Analysis of within‐treatment variances

To test the effect of pH variability on within‐treatment variability in gene expression, we calculated the variance of normalized gene expression for each DEG within each treatment in Experiment 2 (excluding ambient) and averaged the variance across all DEGs (Figure 4). For each pH level (7.85, 7.70), we then calculated F‐ratios by dividing the mean variance of the variable treatment by that of the static treatment. We log‐transformed these variances and compared the distributions to a t‐distribution centred at 0 using one‐tailed t‐tests (Table S11). We found that at both pH levels, the average variance in DEG expression was greater in the variable pH treatment than in the static pH treatment (pH 7.85: p = .0001; pH 7.70: p = .03487).

FIGURE 4

Box plot of within‐treatment variances in Experiment 2 (DEGs only, outliers removed for clarity). Diamonds mark the mean for each treatment. Notches represent a roughly 95% confidence interval around the median. Removed points lie outside of 1.5 times the IQR of each hinge.

In Experiment 1, gene set enrichment analysis using FGSEA revealed 240 enriched gene sets among the upregulated genes and 343 enriched gene sets among the downregulated genes in the pH 7.30 treatment compared to the 7.85 treatment (FDR <0.05; Table S12). In Experiment 2, 61 gene sets were enriched among upregulated genes and 71 among downregulated genes in the 7.70 static treatment compared to the 7.85 static treatment (Table S13). At the 7.85 pH level, we found 44 enriched gene sets among upregulated genes and 202 among downregulated genes in the variable treatment compared to the static treatment (Table S14). At the 7.70 pH level, 115 and 22 gene sets were enriched among the upregulated and downregulated genes, respectively, in the variable treatment compared to the static treatment (Table S15). To aid interpretation, the enriched gene sets for the 7.85/7.30 comparison in Experiment 1 and the 7.85 static/7.70 static comparison in Experiment 2 were further collapsed into clusters of gene sets (using a gene set similarity coefficient) using the AutoAnnotate and clusterMaker2 applications for the Cytoscape software platform. These clusters were manually summarized based on their constituent gene sets (Tables 4 and 5).

TABLE 4

Summary of upregulated and downregulated gene set clusters in Experiment 1

Upregulated in pH 7.30 treatment		Downregulated in pH 7.30 treatment
Categorical cluster	Number of gene sets in each cluster	Categorical cluster	Number of gene sets in each cluster
Mitochondrion, aerobic respiration, mRNA export from nucleus	44	Transmembrane ion transport, regulation of synaptic signalling, ligand‐gated ion channel activity, behaviour, cognition and sensory perception	90
RNA metabolism, processing, splicing, modification, tRNA biosynthesis; ribosome biogenesis	41	Regulation of nervous system development and growth	60
Translation and protein localization	39	Synaptic vesicle membrane, regulation of clathrin‐dependent endocytosis	22
Muscle development	22	Axo‐dendritic transport	20
Organic acid catabolism	15	Synaptic membrane and synapse	19
Muscle contraction and adaptation, myogenesis	14	G protein‐coupled receptor signalling	15
Energy reserve and carbohydrate metabolic process	10	Exocytosis and secretion	14
Proteolysis, mRNA catabolism, negative regulation of cell cycle G2/M phase transition	10	Central nervous system development	12
Peroxisomal organization and transport, protein localization to organelle	8	Regulation of pH and iron ion transport	9
Innate immune response	6	Aminoglycan and glycoprotein metabolic process	8
Telomere maintenance via lengthening and organization	6	Calcium‐dependent phospholipid binding and cell–cell adhesion	8
RNA polymerase II	5	Dopamine secretion and transport	7
Protein modification by small protein conjugation or removal	3	Axon, distal axon and terminal bouton	6
Actin filament binding	2	Dendritic tree and neuron spine	6
Alpha actinin binding	2	GTPase activator activity	6
Cytoplasmic stress granule	2	Positive regulation of MAPK cascade	6
DNA polymerase activity	2	Receptor localization to synapse	6
Mitochondrial matrix and nucleoid	2	Regulation of vesicle fusion	6
Ribosome binding	2	Dendrite membrane	5
RNA helicase activity	2	Ephrin receptor signalling pathway	5
adipogenesis	1	Extrinsic component of cytoplasmic side of plasma membrane	5
ADP binding	1	Microtubule polymerization	5
Allograft rejection	1	Regulation of protein localization to membrane	5
Androgen response	1	Synaptic vesicle transport and localization	5
Cell substrate junction	1	Glycosphingolipid biosynthetic process	4
Cysteine‐type endopeptidase activity	1	Cortical Actin cytoskeleton	3
Fatty acid metabolism	1	Regulation of cell shape	3
Ficolin‐1‐rich granule lumen	1	Vascular transport	3
General transcription initiation factor binding	1	Intrinsic component of Golgi membrane	2
Interferon alpha response	1	Long term depression and vascular smooth muscle contraction	2
Lysine degradation	1	Negative regulation of secretion & transport	2
MYC targets version 1 (Hallmark)	1	Neuron apoptotic process	2
MYC targets version 2 (Hallmark)	1	Regulation of amyloid precursor protein catabolic process	2
Platelet morphogenesis	1	Regulation of neurotransmitter receptor activity	2
Positive regulation mitotic cell cycle	1	Regulation of small GTPase‐mediated signal transduction	2
Receptor signalling pathway via STAT	1	Response to catecholamine	2
rRNA binding	1	Synaptic vesicle recycling	2
Sarcolemma	1	Vesicle docking	2
Sarcoplasm	1	Amyotrophic lateral sclerosis	1
Starch & sucrose metabolism	1	Anchored component of membrane	1
Viral myocarditis	1	Cyclic nucleotide‐mediated signalling	1
		developmental maturation	1
		Endocytosis	1
		Gap junction	1
		Genes upregulated by KRAS activation	1
		Kinesin binding	1
		Long‐term potentiation	1
		Neuron migration	1
		Perinuclear region of cytoplasm	1
		Phosphoprotein binding	1
		Phosphoric diester hydrolase activity	1
		Protein serine threonine kinase inhibitor activity	1
		Regulation of neuron differentiation	1
		Renal system process	1
		Tau protein binding	1

Note: Enriched gene sets (GO, KEGG, hallmark) were clustered by similarity using the AutoAnnotate and clusterMaker2 applications for the Cytoscape software platform. Clusters were then manually examined and named. See Table S12 for the full list of enriched gene sets in this experiment.

TABLE 5

Summary of upregulated and downregulated gene set clusters in Experiment 2 (comparison of static treatments

Upregulated in pH 7.70 treatment		Downregulated in pH 7.70 treatment
Categorical cluster	Number of gene sets in each cluster	Categorical cluster	Number of gene sets in each cluster
RNA processing & splicing, histone methyltransferase complex	26	Immune response	34
Epigenetic regulation of gene expression and chromatin organization	8	Lymphocyte proliferation, differentiation and activation	26
DNA repair, recombination and replication	7	Secretory granule and myeloid leucocyte mediated immunity	13
mRNA export from nucleus	6	Endothelial cell migration and blood vessel morphogenesis	10
E‐box binding	4	JAK–STAT signalling pathway	9
Ubiquitin‐mediated proteolysis	4	Neuropeptide/G protein‐coupled receptor signalling pathway	7
Ubiquitin ligase complex	4	Cellular ion homeostasis	6
RNA phosphodiester bond hydrolysis	3	Positive regulation of MAPK cascade	6
Gene silencing	2	Cell–cell junction assembly	5
Nuclear speck	2	Developmental growth involved in morphogenesis	5
A band	1	Regulation of cytoskeleton and supramolecular fibre organization	5
Cell cortex region	1	Wound healing and regulation of body fluid levels	5
Inositol phosphate‐mediated signalling	1	Leucocyte migration and regulation of chemotaxis	4
Regulation of long‐term synaptic potentiation	1	External side of plasma membrane	3
Single‐stranded RNA binding	1	Leading edge membrane	3
Structural constituent of cytoskeleton	1	Plasma membrane signalling receptor complex	3
Transcription coregulator activity	1	Positive regulation of phagocytosis	3
		Protein complex involved in cell adhesion and integrin‐mediated signalling pathway	3
		Regulation of cytokine production	3
		Cilium movement and cell motility	2
		Collagen‐containing extracellular matrix	2
		Endocytic vesicle	2
		Positive regulation of cell‐substrate adhesion	2
		Receptor‐mediated endocytosis	2
		Regulation of peptidyl‐tyrosine phosphorylation	2
		Guanyl nucleotide binding	1
		Calcium ion binding	1
		Allograft rejection	1
		Membrane microdomain	1
		Complement system	1
		Positive regulation of cell population proliferation	1
		Smooth muscle contraction	1
		Superoxide metabolic process	1
		Response to organophosphorus	1
		Ras protein signal transduction	1
		Response to dopamine	1
		Inflammatory response	1
		Odontogenesis	1
		Coagulation	1
		Leucocyte transendothelial migration	1
		Pigment granule	1
		Cell adhesion molecule binding	1
		Ciliary plasm	1

Note: Enriched gene sets (GO, KEGG, hallmark) were clustered by similarity using the AutoAnnotate and clusterMaker2 applications for the Cytoscape software platform. Clusters were then manually examined and named. See Table S13 for the full list of enriched gene sets in this experiment.

To assess consistency of response to acidification across experiments, we determined the overlap in enriched gene sets between the pH 7.85/pH 7.30 comparison of Experiment 1 and the pH 7.85 static/pH 7.70 static comparison of Experiment 2 (Figure 5, Table 6). To assess consistency in expression response across the two static versus variable treatment comparisons of Experiment 2 (7.85 static/7.85 variable, 7.70 static/7.70 variable), we determined the overlap in enriched gene sets between these comparisons (Figure 6, Table S16).

FIGURE 5

Overlapping enriched gene sets across both experiments. “Up” and “down” refer to gene sets that were up‐ or downregulated in the lower pH treatment relative to the higher pH treatment in each experiment (i.e., pH 7.85 treatments are treated as baseline in both cases). Only static treatments are included for Experiment 2.

TABLE 6

Overlapping enriched gene sets between high pH versus low pH comparisons in Experiments 1 and 2 (upregulated vs. downregulated gene sets in both experiments)

Upregulated in both experiments	Downregulated in both experiments
GOBP ribonucleoprotein complex biogenesis	GOBP MAPK cascade
GOBP ncRNA processing	GOCC side of membrane
GOCC ribonucleoprotein complex	GOBP receptor mediated endocytosis
GOMF catalytic activity acting on RNA	GOCC cell surface
GOBP translational termination	GOBP positive regulation of protein kinase activity
GOBP RNA export from nucleus	GOBP positive regulation of MAPK cascade
GOBP RNA processing	GOCC cell leading edge
GOBP RNA phosphodiester bond hydrolysis	GOBP cell–cell junction organization
GOBP mRNA export from nucleus	GOBP cell–cell adhesion
GOBP nuclear export	GOBP endocytosis
GOBP mRNA metabolic process	GOBP exocytosis
GOCC U2 type spliceosomal complex	GOBP cell–cell junction assembly
GOBP RNA 3′‐end processing	GOBP cell growth
GOBP nucleic acid phosphodiester bond hydrolysis	GOBP taxis
KEGG spliceosome	GOCC secretory granule membrane
GOCC transferase complex	GOBP regulation of anatomical structure morphogenesis
GOBP protein modification by small protein conjugation	GOCC secretory vesicle
GOBP RNA localization	GOMF neuropeptide receptor activity
GOCC spliceosomal complex	GOBP cell junction assembly
GOBP protein modification by small protein conjugation or removal	KEGG cell adhesion molecules cams
GOBP RNA splicing	GOCC plasma membrane protein complex
GOBP mRNA processing	GOMF calcium ion binding
GOCC nuclear protein‐containing complex	GOCC cell projection membrane
GOCC intracellular protein‐containing complex	GOCC plasma membrane signalling receptor complex
	GOBP cell–cell adhesion via plasma membrane adhesion molecules
	GOBP developmental growth involved in morphogenesis
	GOBP developmental cell growth
	GOBP neuropeptide signalling pathway
	GOBP adenylate cyclase inhibiting G protein‐coupled receptor signalling pathway
	GOCC leading edge membrane
	GOCC vesicle membrane
	GOMF G protein‐coupled receptor activity
	GOCC receptor complex
	KEGG neuroactive ligand receptor interaction
	GOBP G protein‐coupled receptor signalling pathway
	GOMF molecular transducer activity

FIGURE 6

Overlapping enriched gene sets across static versus variable comparisons in Experiment 2. “Up” and “down” refer to gene sets that were up‐ or downregulated in the “variable” treatment relative to the “static” treatment for a given pH level (7.85 or 7.70; i.e., static treatments are treated as baseline).

DISCUSSION

Impacts of static acidification

Numerous studies have demonstrated impaired behaviour and sensory function in fish and other marine organisms when exposed to low pH/high pCO2 (Domenici et al., 2012; Hamilton et al., 2017, Hamilton et al., 2013; Munday et al., 2010; Pistevos et al., 2015, but see Clark et al., 2020a, 2020b). Though the mechanisms behind these changes are still poorly understood, significant effects of low pH/increased pCO2 on brain gene expression have been documented in a few marine fish species that demonstrate associated impairments in behaviour (Lai et al., 2016; Schunter et al., 2016, 2018, 2021). In this study we tested whether brain gene expression is similarly impacted in a temperate reef fish that experiences prolonged periods of natural acidification. Across both experiments presented here, we found that global gene expression was significantly affected by acidification (high vs. low pH). Comparing results across experiments, the number of detected DEGs increased with more extreme acidification, as did the number of enriched gene sets. A similar increase in DEGs with increased intensity of acidification has also been reported in the olfactory bulb of coho salmon (Williams et al., 2019). This marked increase in effect size indicates that further acidification past the already‐low pH of 7.70 can have a substantial additional impact on the physiology of marine fish. This pattern may have important implications for the management of marine ecosystems and the services they provide as our global society struggles to control CO2 emissions.

Although a greater number of gene sets were enriched in Experiment 1 than in the comparison of the static treatments of Experiment 2, similar enrichment themes emerged. In both experiments, static acidification led to the upregulation of gene sets related to turnover in the proteome and transcriptome that may reflect ongoing physiological adaptation to altered environmental conditions (Tables 4 and 5). Additionally, static acidification in both experiments led to the downregulation of gene sets related to the MAPK cascade, G protein‐coupled receptor signalling pathways, plasma membrane components, secretory vesicles and granules, neuroactive ligand‐receptor interaction, and calcium ion binding, indicating a general reduction in cell signalling, including neuroactive signalling, in response to high pCO2. In general, this is the opposite of the response seen in similar gene sets in spiny damselfish (Acanthochromis polyacanthus) (Schunter et al., 2018) and the olfactory bulb of coho salmon (Williams et al., 2019). Schunter et al. (2019) proposed that high pCO2‐induced changes in electrochemical gradients across GABAergic neuron membranes may initiate a “vicious cycle” of feedbacks and ultimately an increase in excitatory activity in the brain that may explain behavioural changes seen in other species. If this is indeed the case, the downregulation of gene sets related to neuroactive signalling seen here may represent a species‐specific adaptive response aimed at combating maladaptive runaway excitation in acidic waters (see discussion of GABAA receptor related genes below). Finally, downregulation of gene sets related to growth and morphogenesis, cell–cell adhesion, and the cytoskeleton indicate potential disruption of cell growth and development due to increased cellular stress. Similar themes of upregulated transcription and cellular stress response have also been documented in the muscle tissue of Pacific rockfish (Hamilton et al., 2017).

We also identified divergent sets of genes enriched between the moderate (Experiment 2, target pH 7.70) and extreme (Experiment 1, target pH 7.30) acidification treatments, indicative of a potential threshold effect as static pH decreases. In comparison to the static acidification in Experiment 2, static acidification in Experiment 1 resulted in the up‐ and downregulation of additional gene sets related to metabolic processes (Table 4). These changes may again indicate further shifts to the synthesis of stress response proteins, or to isoforms that are better suited to an altered cellular environment. Because, at least in humans, there can be interaction/crosstalk between cellular stress response pathways and the innate immune system signalling pathways (Muralidharan & Mandrekar, 2013), the upregulation of an additional six gene sets related to the innate immune response may further indicate increased cellular stress. The acidification in Experiment 1 also resulted in the downregulation of broad categories of gene sets related to basic neurological functions, behaviour, and cognition, which supports the hypothesis that acidification can lead to behavioural impairment in this species though, as mentioned above, the specific mechanisms through which OA induced alterations in neurobiology might impact fish behaviour are still not well understood (Tresguerres & Hamilton, 2017).

In Experiment 2, additional gene sets related to the regulation of gene expression (including epigenetic regulation) were upregulated, again indicating a systemic shift in gene expression and response to cellular stress. We also found a unique downregulation of a large number of gene sets related to immune response, which may in part reflect the external conditions of this experiment, particularly the unusually warm ambient temperatures. The combination of physical and chemical stressors may have led to the suppression of the immune system in fish in the acidified treatment. Suppression or dysregulation of immune function is a well‐established response to stress (Dhabhar, 2014), and heat stress‐induced immunosuppression, specifically, has been noted across various animal systems (Nardone et al., 2010).

The biological themes of the enriched gene sets in both experiments are consistent with enriched categories identified in previous studies in tropical reef fish (e.g., Schunter et al., 2016, 2018) and salmon (Williams et al., 2019). Interestingly, however, the pattern of enrichment in E. jacksoni under acidified conditions is generally opposite to the enrichment pattern found by Schunter et al. (2018) when comparing acute or developmentally (together: cis‐generationally) exposed spiny damselfish to control (untreated) individuals, but closely resembles the pattern of gene set enrichment that Schunter et al. found when comparing transgenerationally exposed A. polyacanthus to those that were developmentally exposed to acidified conditions (not the control treatment). The contrast of our results may reflect the transgenerational and evolutionary exposure history of E. jacksoni populations to naturally acidic environments. While A. polyacanthus on coral reefs may experience diurnal pCO2 fluctuations on the scale of ±50–150 μatm (Schunter et al., 2021 and references therein), E. jacksoni in upwelling regions are likely to regularly experience prolonged increases in pCO2 (days to weeks) from as low as ~300 to >1000 μatm (Chavez et al., 2018; Donham et al., 2022). Those living in estuaries can experience even greater shifts in carbonate chemistry over even shorter timescales (Duarte et al., 2013; Hofmann et al., 2011). It is therefore likely that the juvenile E. jacksoni in our experiments were transgenerationally exposed to acidified conditions in situ. Additionally, while both species lack a pelagic larval stage, A. polyacanthus is a substrate spawner and E. jacksoni is a live‐bearing species. This means that the E. jacksoni used for this experiment may also have developmentally experienced their mothers' natural environmental exposures prior to their birth. Because of this potential in‐situ transgenerational exposure, our experimental design may be more comparable to the transgenerational treatment used by Schunter et al. (2018).

As mentioned above, it has been suggested that the cause of previously documented acidification‐induced behavioural changes in fish may not only be due to a reversal of electrochemical gradients that flip the nature of GABAergic neurons from inhibitory to excitatory (Nilsson et al., 2012), but also due to a positive feedback cycle that may develop as a response to this increase in excitatory activity in the brain (Schunter et al., 2018, 2019). This proposed response consists of an increase in GABA release and in the abundance of GABAA receptors, which under nonacidified conditions would serve to reduce overactivity in the brain, but under acidified conditions probably act to exacerbate the overactivity. Some previous studies in fish have seen changes in expression consistent with this response, such as increased expression of GABAA receptor subunits and transporter genes (e.g., Lai et al., 2016; Schunter et al., 2016, 2018). However, the fish in Experiment 1 showed the opposite response in GABA‐related genes. In Experiment 1, GABAA receptor subunit isoforms α (1–6), β (1–3), γ (1–3), ρ (2), and π were all downregulated in the pH 7.30 treatment, along with many other GABA signalling genes, including glutamate decarboxylases gad1 and gad2 (Lai et al., 2016) and gabarapl2. Interestingly, a similar general downregulation of GABAergic signalling pathways was recently noted in A. polyacanthus at CO2 seeps, but not in other reef fish species (Kang et al., 2022). A study on Pacific coho salmon (Williams et al., 2019) also found no changes in GABAA receptor subunit expression in the olfactory bulb under increased pCO2, but did find an increase in the expression of a GABAB receptor subunit (gabbr2), which was instead downregulated in E. jacksoni in our Experiment 1. Williams et al. (2019) also found significant changes in the expression of other genes associated with GABA signalling, including downregulation of the slc6a13 gene involved in GABA uptake, which we also saw downregulated in E. jacksoni in Experiment 1 (see Table S17 for all differentially expressed GABA‐related genes). These divergent responses between E. jacksoni, Pacific salmon, and tropical fish species could represent species‐specific adaptation to differing environmental conditions. In the case of E. jacksoni, which frequently experiences periods of high pCO2, the downregulation of GABA‐related genes under high pCO2 may be an adaptation that prevents or interrupts the excitatory positive feedback cycle proposed by Schunter et al. (2019). Previous studies have also noted opposite responses in gene expression across species of the same taxa (Kang et al., 2022; Strader et al., 2020), and even across populations of the same species (Goncalves et al., 2016), but the extent of the role that transgenerational effects play in creating divergent responses is still unclear (but see Goncalves et al., 2016; Schunter et al., 2018). Importantly, however, our seemingly species‐specific results may indicate that E. jacksoni is preadapted to acidified conditions, whether through long‐term local adaptation or transgenerational plasticity. Because of its limited adult dispersal and lack of a pelagic larval phase, E. jacksoni may be more likely to be genetically adapted to local conditions than other species (Warner, 1997), and its live‐bearing reproduction may also facilitate adaptation through maternal effects. Kang et al. (2022) recently proposed a similar hypothesis to explain why A. polyacanthus (which also lacks a pelagic larval stage) differed from other co‐occurring damselfish species in its molecular response to elevated pCO2.

Interestingly, the response of GABA‐related genes to acidification varied between Experiments 1 and 2 (which used different levels of acidification). In response to the more moderate static acidification in Experiment 2, E. jacksoni showed an upregulation of two subunits of the GABAA receptor (gabra6 and gabrb3), which were instead downregulated in Experiment 1. Interestingly, the gabra6 subunit is also upregulated in spiny damselfish transgenerationally exposed to high pCO2 when compared to those that were only developmentally exposed, an effect opposite to that seen in the expression of other GABAA subunits in the same experiment (Schunter et al., 2018). No other GABA‐related genes were significantly affected by this treatment. In Experiment 1, the greater magnitude change in pH resulted in an opposite and much broader response of GABA‐related genes. These conflicting responses in the transcription of GABAA receptor subunits and other GABAA‐related genes indicate that in addition to varying across species, the response of the GABA signalling pathway to acidification/high pCO2 may also depend on the magnitude of the environmental change. Further study is needed to determine how the divergent transcriptomic response of E. jacksoni seen in our experiments translates to behaviour and overall fitness, how the magnitude of any emergent effects compare to those observed in other species, and the role transgenerational exposure plays in E. jacksoni response to acidification.

In both Experiments 1 and 2, an additional group of gene sets related to muscle tissue were identified as enriched in the acidified treatment. In Experiment 1, this included the upregulation of gene sets related to muscle development, contraction, and adaptation and muscle cell components, as well as the downregulation of the vascular smooth muscle contraction GO gene set. In Experiment 2, the A band GO gene set was upregulated, while the smooth muscle contraction gene set was downregulated. While smooth muscle is present in blood vessels in the brain, it is possible that the identification of some of these pathways (such as those related to striated muscle tissue) as enriched is in part due to the misannotation of genes in this nonmodel species to orthologous reference genes. Alternatively, because our transcriptome was assembled using both brain and muscle tissue, it is possible that some brain transcripts were misaligned to muscle‐exclusive reference transcripts during differential expression analysis.

Impacts of pH variability

Overall, we found that variability in pH moderated the differential gene expression seen under static acidification. pH variability decreased the number of DEGs detected by the edgeR analysis between the pH 7.85 and pH 7.70 treatments in Experiment 2 (from 159 genes when treatments were static to nine genes when both treatments were variable). This aligns with two previous studies that found that effects of pH on fish gene expression and behaviour were diminished by the incorporation of diel pH fluctuations (Jarrold et al., 2017; Schunter et al., 2021).

Functional enrichment analysis revealed many up‐ and downregulated gene sets between static and variable treatments at each pH level, though there were more enriched gene sets in the more moderate pH 7.85 comparison (329 gene sets) than the pH 7.70 comparison (139 gene sets). This difference may represent an acidification threshold nearer to the 7.85 treatment, where the majority of transcriptional adaptation to acidification is activated. Such a threshold effect in gene expression patterns has also been observed in the gill tissue of spider crabs exposed to two levels of acidification (Harms et al., 2014), as well as in the muscle tissue of blue rockfish (Sebastes mystinus; Hamilton et al., 2017), and thresholds in OA response have been noted across taxa (Bednaršek et al., 2021; Castillo et al., 2014; Wittmann & Pörtner, 2013). Interestingly, although 33 gene sets were commonly enriched across the pH 7.85 and 7.70 static‐variable comparisons, the majority of them (30) were enriched in opposite directions depending on the pH level, with variability at pH 7.85 eliciting a directional response mirroring that of static acidification, and variability at pH 7.70 eliciting the opposite response (Figure 6; Table S16). For example, at the 7.85 pH level, variability led to a downregulation of gene sets related to morphogenesis, development, cell differentiation, exocytosis, cell–cell adhesion, molecular transducer activity, and leucocyte mediated immunity, while variability at pH 7.70 led to upregulation in these gene sets compared to the static treatment. These contrasting responses indicate that pH variability can have opposing effects on brain physiology depending on the underlying mean pH level. This interactive effect of acidification and variability may again reflect a threshold in the neural response of fish to acidification. It may be that at more moderate pH levels, variability exacerbates the negative effects of acidification by temporarily dropping the pH further below the average, but under more extreme acidification, perhaps past a biological tipping point, any negative effects of further acidification introduced by temporary oscillations may be outweighed by the temporary relief provided by the upswing of the oscillations above the mean pH.

It is important to note that our interpretation of these results could be limited by the scope of our experimental design. In Experiment 2, we sampled tissue from individuals in each treatment when the variable treatments were increasing in pH and intersecting their corresponding static treatments. While this design keeps the pH at the time of sampling consistent between the static and variable treatments, it assesses expression at only a single time point, and therefore does not account for likely divergent expression patterns at different positions in the pH cycles of the variable treatments. Additional experiments are necessary to determine if and how gene expression differs in E. jacksoni depending on the trajectory and value of the pH at the time of sampling.

Interindividual variability in gene expression

A particularly striking finding from our experiments is the observation that gene expression variability across individuals was greater in the variable pH treatments of Experiment 2 than in the static treatments (Figure 4). This pattern indicates that the environmental variability introduced by the pH oscillations may be revealing significant “cryptic variation” (Rutherford, 2000, 2003; Rutherford & Lindquist, 1998) in the transcriptomic response of E. jacksoni to acidification. In the context of climate change, such phenotypic variation, if beneficial and heritable, could represent potential adaptive variation on which selection may act, allowing populations to adapt to ongoing changes in environmental conditions (Rutherford & Lindquist, 1998; Rutherford, 2000, 2003; Queitsch et al., 2002; reviewed in Ghalambor et al., 2007).

Patterns of expression across individuals within static treatments, and across functionally related genes within individuals, were notably consistent. This consistency provides evidence of a conserved stress response as described by Kültz (2005), and may again reflect a biochemical “switch” type response, activated at a certain environmental threshold. This idea is further supported in Experiment 2 by the similarity of expression profiles of some individuals in the pH 7.85 variable treatment to the expression profiles exhibited by those in the pH 7.85 static treatment, while others in the variable treatment exhibited expression profiles similar to those in the pH 7.70 static treatment (Figure S5; Figure 3).

Because we did not use the whole brain, and instead arbitrarily subsampled brain tissue from each individual, some of the interindividual variability in expression profiles may be the result of variability in the exact section(s) of the brain that was sampled for each individual. Conversely, it is possible that this sampling method could introduce treatment‐level bias in the brain region sampled that could lead to misleading signals of differential expression between treatments. However, expression profiles within static treatments were remarkably consistent across individuals, especially in Experiment 1 (Figure 2), indicating a low probability of sampling bias, and all individuals were subsampled in an arbitrary manner by a single researcher for each experiment. We therefore maintain that alternative sampling methods would have been unlikely to change the major patterns and conclusions presented here.

CONCLUSIONS

Overall, our results indicate that both acidification and pH/pCO2 variability can have significant impacts on the brain gene expression of a nearshore temperate fish species. Given recent debate regarding the generality of neurological impacts of OA on marine fish (Clark et al., 2020a, 2020b; Munday et al., 2020), our study provides evidence of neurological impacts, even in a species with a high likelihood for local adaptation to naturally low pH/pCO2. We found a significant effect of acidification on global gene expression in E. jacksoni brain tissue, and that the transcriptomic response was similar to a previous experiment that compared transgenerationally exposed tropical damselfish to individuals that were developmentally exposed (Schunter et al., 2018). These results suggest that the E. jacksoni in our experiments were exposed to ecologically relevant pH/pCO2 variability in situ, which may have influenced their response to acidification in the laboratory. Additionally, our results demonstrate that the incorporation of upwelling‐scale pH variability into acidification treatments has a substantial impact on the number of DE genes detected between moderate and low levels of acidification, indicating that temporal pH variability can moderate the impacts of acidification. Interestingly, we also found that the direction of the effect of variability on gene expression in certain genes depended on the degree of acidification. These opposing patterns of gene expression indicate that the impact of pH variability on fish brain physiology may be context‐dependent, perhaps serving as an additional stressor at more moderate levels of acidification, but as an ameliorating factor when the mean pH is more extreme. Finally, we observed significant variation in gene expression across individuals, and found that upwelling‐scale pH variability revealed additional cryptic phenotypic variation. This finding indicates that studies employing only static treatments may underestimate standing genetic variation in traits related to the response of fish to acidification. This cryptic variation may provide additional genetic variation on which selection may act and therefore increase the likelihood of successful adaptation of fish populations to acidification. In summary, our results emphasize the importance of considering environmental variability in global change experiments and demonstrate that a species with an evolutionary history of exposure to acidified and variable conditions exhibits a distinctive transcriptomic response in gene sets similar to those affected in species that have shown behavioural impairment.

AUTHOR CONTRIBUTIONS

JAT, KJK, GB, CAL, and YT conceptualized and designed the experiments. JAT and YT designed and tested the pH‐manipulation system. GB carried out and led the sequencing project for Experiment 1. JAT carried out Experiment 2 and led the sequencing project with laboratory assistance from CAL. GCL led the reference genome sequencing project. JAT conducted all data processing and analyses with input from KJK, GB, CAL, and YT. JAT drafted the paper. All authors reviewed and edited the drafts.

CONFLICT OF INTEREST

The authors declare no competing interests.

Appendix S1

Click here for additional data file.

Tables S12‐S15

Click here for additional data file.