Brain transcriptome of gobies inhabiting natural CO2 seeps reveal acclimation strategies to long‐term acidification

Abstract Ocean acidification (OA) is known to affect the physiology, survival, behaviour and fitness of various fish species with repercussions at the population, community and ecosystem levels. Some fish species, however, seem to acclimate rapidly to OA conditions and even thrive in acidified environments. The molecular mechanisms that enable species to successfully inhabit high CO2 environments have not been fully elucidated especially in wild fish populations. Here, we used the natural CO2 seep in Vulcano Island, Italy to study the effects of elevated CO2 exposure on the brain transcriptome of the anemone goby, a species with high population density in the CO2 seep and investigate their potential for acclimation. Compared to fish from environments with ambient CO2, gobies living in the CO2 seep showed differences in the expression of transcripts involved in ion transport and pH homeostasis, cellular stress, immune response, circadian rhythm and metabolism. We also found evidence of potential adaptive mechanisms to restore the functioning of GABAergic pathways, whose activity can be affected by exposure to elevated CO2 levels. Our findings indicate that gobies living in the CO2 seep may be capable of mitigating CO2‐induced oxidative stress and maintaining physiological pH while meeting the consequent increased energetic costs. The conspicuous difference in the expression of core circadian rhythm transcripts could provide an adaptive advantage by increasing the flexibility of physiological processes in elevated CO2 conditions thereby facilitating acclimation. Our results show potential molecular processes of acclimation to elevated CO2 in gobies enabling them to thrive in the acidified waters of Vulcano Island.

Most knowledge about these impacts of ocean acidification (OA) on teleosts, however, comes from short-term laboratory-based studies and mostly only accounts for short-term effects of OA. Accurate assessment of the future biological impacts of OA on species and populations requires consideration of their potential for acclimation and adaptive evolution (Sunday et al., 2014). The effects of elevated CO 2 are also known to be variable both within and across species.
Inter-individual variability in CO 2 tolerance has been reported both in laboratory studies and in the wild (Cattano et al., 2018;Lehmann et al., 2022;Monroe et al., 2021;Munday, McCormick, et al., 2013;Paula et al., 2019;Welch & Munday, 2017). The underlying reason for the observed variable effect is not fully known but could be the result of specific genetic and epigenetic signatures or environmental factors or a combination of those. Environmental factors such as nutrient gradients, species composition, habitat modifications and predation risk influence the effects of elevated CO 2 on various aspects of fish physiology and behaviour (Ferrari et al., 2017;Mirasole et al., 2019;Nagelkerken et al., 2015;Zunino et al., 2019). Studies on wild fish populations that naturally inhabit high CO 2 environments provide a unique opportunity for integrating the direct and indirect effects of elevated CO 2 exposure on fish and investigate the molecular mechanisms enabling acclimation.
Volcanic CO 2 seeps in the ocean, where the emission of CO 2rich fluids locally reduce pH, have been used as natural analogues to study the effects of long-term exposure of fish to elevated CO 2 in their natural environment. Shifts in habitat, trophic diversity and fish assemblages at CO 2 seeps have been reported which could obscure some of the negative effects of OA on fish (Fabricius et al., 2011;Hall-Spencer et al., 2008;Mirasole et al., 2019;Nagelkerken et al., 2015;Vizzini et al., 2017). For example, a decrease in predator abundance at the CO 2 seep could offset the increased mortality of some fish species as a consequence of sensory and behavioural impairment resulting from exposure to elevated CO 2 levels .
Habitat shifts are another factor that can work antagonistically with CO 2 to influence species behaviour. Gobies from the CO 2 seeps in Vulcano Island, Italy living in more open habitats showed reduced predator escape behaviour, however, the same was not observed in gobies living in more sheltered sites indicating the influence of habitat on species response to elevated CO 2 exposure (Nagelkerken et al., 2015). Similarly, habitat shifts resulting in net resource enrichment at CO 2 seeps seems to benefit certain fish species by providing sufficient food to meet the increased energy demands of living in a high CO 2 environment (Connell et al., 2013;Mirasole et al., 2019;Nagelkerken et al., 2015). It is therefore important to consider the indirect effects of OA while making inferences about the adaptive potential of species to future OA conditions and volcanic CO 2 seeps provide the perfect opportunity to do so.
Prolonged exposure to elevated CO 2 at the seeps may enable fish to acclimate through physiological and behavioural compensation resulting in long-term survival and subsequent adaptation.
Indeed studies using transcriptomics to determine gene expression patterns in the brain and gonads of fish from natural CO 2 seeps have revealed plasticity and adaptive selection to elevated CO 2 (Kang et al., 2022;Petit-Marty et al., 2021). These findings suggest that species that are continuously exposed to elevated CO 2 levels can develop mechanisms to mitigate the negative effects of CO 2 exposure and have the potential to adapt to acidified environments. However, given the relatively small spatial range of CO 2 seeps compared with the pelagic larval dispersal distance of most marine species, the population in the CO 2 seep will mostly come from a heterogeneous larval pool. The recently settled larvae in the seep must therefore rapidly respond to the acidified environment via phenotypic plasticity, resulting in physiological adjustments to ensure survival. These plastic responses established during development could be crucial to allow organisms to acclimate and eventually lead to adaption to environments with elevated CO 2 levels (Munday, Warner, et al., 2013;Sunday et al., 2014). As the CO 2 levels continue to rise developmental plasticity could be a crucial process facilitating the establishment of adaptive processes across species enabling them to survive in the increasingly acidifying oceans of the future.
While several controlled laboratory experiments (Esbaugh, 2017;Heuer & Grosell, 2014;Ishimatsu et al., 2005Ishimatsu et al., , 2008Kelly & Hofmann, 2013;Lehmann et al., 2022;Strader et al., 2020;Tresguerres & Hamilton, 2017) have investigated the effects of elevated CO 2 on fish physiology, behaviour and molecular processes, we are only starting to learn about the effects of CO 2 on brain function in wild fish species. Here, we sequenced and analysed the brain transcriptome of the anemone goby, Gobius incognitus to identify gene expression patterns underlying acclimation to elevated CO 2 in the wild population. This benthic fish species has a very limited home range making it an ideal model to study the effects of elevated CO 2 exposure in situ. Previous studies have reported potentially negative physiological effects such as changes in otolith morphology and skeletal composition in Gobius from the natural CO 2 seeps at Vulcano Island (Mirasole et al., 2017(Mirasole et al., , 2021, however, they still maintain higher population densities at the seeps (Nagelkerken et al., 2015). This could mean that gobies have developed alternative strategies/mechanisms to offset the negative effects of elevated CO 2 exposure. In fact, swimming activity and predator perception of Gobius incognitus living in natural CO 2 seeps in Vulcano Island were not affected by elevated CO 2 indicating higher tolerance of this species to changing ocean conditions (Spatafora et al., 2022). By examining brain gene expression patterns of the anemone goby from the natural CO 2 seeps in Vulcano Island, Italy, we aim to identify molecular signatures of acclimation that enable this species to mitigate the negative effects of elevated CO 2 and thrive in the acidified waters of Vulcano Island. pH 8.19 ± 0.03) located around 600 m from the low pH site (Figure 1).
The control site was similar to the low pH site in terms of orientation (southeast), depth and habitat (sandy with beds of Cymodocea nodosa seagrass). A total of eight individuals were sampled at each site at the same time of day (between 11:00 am and 12:00 pm). The total length of all fish was measured in the field, immediately after which they were euthanized on ice and the brain tissues were dissected, preserved in RNAlater, and stored at −80°C until further processing.
The age of all 16 individuals was estimated using the Von Bertalanffy growth function (VBGF; Von Bertalanffy, 1938) which is based on a bioenergetic expression of fish growth. Size-at-age data can be used to fit the parameters of a growth function, such as the VBGF. In our case, we extrapolated the age of our analysed Gobius specimens using the growth data from the work of Sasal et al. (1996). On average all the sampled fish were 4 cm long and 2 years old (Table S1).

| RNA extraction and sequencing
Frozen whole brains were homogenized in RLT-lysis buffer (Qiagen) for 30 s at 30 Hz/s using Qiagen TissueLyzer II with 4-5 single-use silicon beads per sample and total RNA was extracted using the Qiagen RNeasy Micro Kit. On-column DNase digestion was performed using RNase-Free DNase1 (Qiagen) following the manufacturer's instructions, to remove any DNA contamination. RNA concentration and quality were determined using Qubit and Agilent 2100 Bioanalyzer, respectively (mean RIN = 7.79), following which RNA-Seq libraries were built using the KAPA-stranded mRNA-Seq kit with 100 ng total RNA as input from each of the 16 individuals.
The samples were subsequently sequenced on an Illumina NovaSeq 6000 to get 151 bp, paired-end reads at the Centre for PanorOmics Sciences, The University of Hong Kong. An average of 33 ± 2.8 million read pairs were obtained from the 16 RNA libraries sequenced (Table S2).

| Sequence processing and primary de novo transcriptome assembly
The quality of all the sequencing reads was determined using FastQC v0.11.5 (Andrews, 2010). To ensure the high quality of all reads, we trimmed adapters and sequences with an average Phredscaled quality score below 30 using fastp v0.20.0 (Chen et al., 2018) with a 4-base sliding window approach. Sequences less than 40 bp in length after trimming were removed. The trimmed sequences were further filtered to remove potential contaminants, ambiguous bases and rRNA sequences. Potential contaminant sequences were removed using kraken v2.0.8-beta (Wood & Salzberg, 2014 ambiguous base calls (Ns), the longest sub-sequence of at least 40 bp without any ambiguous bases was extracted and retained if it was at least half the original sequence length and had an average Phred-scaled quality score above 30 using the 'filter_illumina' script from DRAP (Cabau et al., 2017). Reads originating from ribosomal RNA (rRNA) were identified and removed by mapping the sequences to the SSUParc (SILVA_138_SSUParc_tax_silva.full_metadata. gz) and LSUParc data sets (SILVA_132_LSUParc.full_metadata.gz) from SILVA (Quast et al., 2013) using Bowtie2 v2.4.1 (Langmead & Salzberg, 2012) in the very sensitive, local mode. An average of 26 ± 2.5 million adapter-free, filtered, quality-trimmed and decontaminated paired reads were finally retained after the stringent filtering process from all 16 samples (Table S2). These high-quality reads were subsequently merged and assembled de novo using Trinity v2.8.4  within DRAP v1.92 (parameters: -c '-x'; --no-trim; -m bwa; --no-rate) using the pooled-assembly strategy resulting in a primary transcriptome assembly consisting of 633,517 contigs.

| De novo transcriptome assembly optimization
Candidate coding regions (open reading frames, ORFs) within the transcript contigs assembled by DRAP were identified using TransDecoder v5.5.0  and only contigs with ORFs of at least 100 amino acids long were retained to prevent retention of false-positive ORFs as done in previous studies on CO 2 seeps (Kang et al., 2022;Petit-Marty et al., 2021). These contigs were then blasted against the UniProtKB/Swiss-Prot Release 2021_01 database (Apweiler, 2004) using the BLASTX algorithm (Blast+ v2.6.0; Altschul et al., 1990) with an e-value cut-off 1e-05 and retaining only one blast hit per contig (--max_target_seqs = 1). For each assembled contig only the best ORF was retained prioritized first based on blast homology and then based on ORF length (--retain_ blastp_hits; --single_best_only) after which 113,895 contigs were retained. The quality of the selected contigs was further assessed using the runAssesment module within DRAP v1.92 and 68 putative chimeric contigs were removed. The high-quality RNA-seq reads were then re-mapped to the filtered transcriptome assembly using Bowtie2 v2.4.1 in the sensitive mode and disallowed any discordant and unpaired alignments but without any upper limit on the number of alignments. To identify and cluster individually assembled contigs belonging to the same gene, the alignment results were used by Corset v1.09 (Davidson & Oshlack, 2014) to hierarchically cluster contigs based on the proportion of multi-mapped reads and the expression patterns. For each of the Corset-generated clusters, the longest transcript sequence was extracted using the fetchClus-terSeqs.py script from Corset-tools (https://github.com/Adamt arant o/Corset-tools) and those clusters with no aligned reads was removed from the de novo transcriptome assembly resulting in 44,614 high-quality contigs retained.
To further validate the transcriptome assembly, blast search (TBLASTX, e-value = 1e-05, max_target_seqs = 1) was performed against the complete cDNA sequence of Neogobius melanostomus (round goby) obtained from Ensembl (release-104; Howe et al., 2021), the phylogenetically closest species with a sequenced genome to Gobius incognitus (Thacker & Roje, 2011). To further reduce transcript redundancy and identify potential alternative transcripts that failed to cluster in the previous steps, the blast results were filtered to retain only the longest transcript for each blast hit to get the final de novo transcriptome assembly. Transcripts with no blast hit to Neogobius melanostomus were also retained in the final assembly as these could potentially be transcripts that are specific to our target species, Gobius incognitus. The final transcriptome assembly consisted of 43,349 non-redundant transcripts. The basic statistics of the transcriptome assembly determined using rnaQUAST v2.1.0 (Bushmanova et al., 2016), 'stats.sh' script from BBMap v38.87 (Bushnell, 2014) and 'contig_ExN50_statistic.pl' from Trinity accessory scripts (https://github.com/trini tyrna seq/trini tyrna seq/ tree/maste r/util/misc), are shown in Table S3 and Figure S1.

| Transcriptome annotation and evaluation of assembly completeness
The de novo assembled transcriptome was compared against the UniProtKB/Swiss-Prot Release 2021_01 database using BLASTX (evalue = 1e-05 and HSP length cut-off = 33) resulting in 43. 6% (24,454) of the assembled transcripts with blast hits to known proteins. Gene Ontology (GO) mapping and annotation were then carried out within OmicsBox v1.4.11 (Götz et al., 2008). InterProScan was used to identify the protein domains within the transcript sequences, and the InterProScan GOs were merged with the GO annotations resulting in a total of 24,875 (42.6%) annotated transcripts. For contigs that did not blast to the Swiss-Prot database, the Ensembl gene ID obtained by blasting against Neogobius melanostomus was retained.
Combining the annotation results from OmicsBox and the additional 9362 transcripts with blast hits to the cDNA sequence of Neogobius melanostomus, 78.9% (34,237) of the transcripts in the final de novo assembly has a functional annotation (Tables S3 and S4).
To incorporate the transcript expression levels in evaluating the contiguity of the de novo assembly, the ExN50 (expression-informed N50) metric was determined. ExN50 (expression-informed N50) represents the N50 value for the most highly expressed transcripts representing x% of the total normalized expression data and is more appropriate for transcriptome assemblies than the conventional N50 statistic (Galachyants et al., 2019). A peak of the ExN50 value at higher × percentiles indicates the effectiveness of the de novo assembly pipeline to assemble longer transcripts (Sahraeian et al., 2017). The expression levels of the transcripts were quantified using RSEM v1.3.1 (Li & Dewey, 2014) and the ExN50 value was calculated across different expression percentiles using the TMM normalized expression values. The Gobius incognitus reference transcriptome assembled here has a sufficiently high coverage of the longer transcripts indicated by a peak of ExN50 value at E86 (Table S5 and Figure S2). The level of completeness of the transcriptome assembly was determined using BUSCO v4.1.3 (Manni et al., 2021) using the Actinopterygii_odb10 and Eukaryota_odb10 sets of Benchmarking Universal Single-Copy Orthologs (BUSCO) which revealed high representation of both the Actinopterygian core genes and the Eukaryota core genes. Specifically, 76.4% and 91.4% of Actinopterygii core genes and Eukaryota core genes, respectively, were recovered completely. An additional 7.7% of Actinopterygii core genes and 7.1% of Eukaryota core genes were recovered partially ( Figure S3) indicating sufficiently good assembly and annotation completeness, similar to other de novo assemblies (Pomianowski et al., 2021;Wong et al., 2019).

| Differential expression analysis
The high-quality, filtered RNA-seq reads were re-mapped to the final transcriptome assembly using Bowtie2 v2.4.1 in the sensitive mode, disallowing any discordant and unpaired alignments and using additional parameters (--dpad 0 --gbar 99,999,999 --mp 1,1 --np 1 --score-min L,0,-0.1) to ensure the resulting alignment files are compatible with RSEM, which was used to quantify the transcript expression levels. On average 72% of the RNA-Seq reads mapped back to the assembled transcriptome (Table S6) indicating the sufficiently high quality of the de novo transcriptome. Principal component analysis (PCA) was performed to determine the overall transcript expression patterns and identify any outlier samples in R v4.1.1 (R Core Team, 2021) following which differential expression analysis was done using DESeq2 v1.32.0 (Love et al., 2014). The transcripts were considered to be significantly differentially expressed between the control and CO 2 seep sites if the FDR corrected p-value (p adj ) was less than 0.05. Additionally, to remove further false positives, we set other filters such as absolute log2 fold change (abs log2FC) greater than 0.3 and baseMean greater than 10 as done in previous studies on brain fish transcriptomics (Schunter et al., 2016(Schunter et al., , 2018 (Foll & Gaggiotti, 2008) with the assumption that the individuals from the control and CO 2 seep sites belong to two different populations.  Table S7). The majority of these (68%) were expressed at higher levels in fish from the CO 2 seep site compared to the control site ( Figure S4).

| Functional analysis of DE transcripts
Functional enrichment analysis identified 434 GO terms to be enriched among the DE transcripts (Table S8). The transcripts associated with each enriched GO term was further categorized into six broader functional groups (acid-base and ion homeostasis, neural activity and related effects, circadian rhythm, metabolism, cellular stress response and immune response) based on their functional description from the NCBI gene database (https://www.ncbi.nlm.nih. gov/gene/) and UniProt knowledgebase (UniProtKB; https://www. unipr ot.org/; Figure 3 and Table S9).
The five most significant DE transcripts (with the lowest p adj values) were all associated with the circadian rhythm. These include two core circadian genes Period circadian protein homolog 3 (PER3) and Orphan nuclear receptor RVR (NR1D2), two nuclear hormone receptor genes involved in transcriptional regulation of core circadian genes Nuclear receptor ROR-alpha A (RORA) and Nuclear receptor ROR-beta (RORB) and a transcriptional regulator known to regulate circadian rhythm genes Nuclear factor interleukin-3-regulated protein (NFIL3; Figure S4). In fact, all core clock transcripts were DE between the control and CO 2 seep site; specifically, there was an overall downregulation of repressors of the circadian rhythm and upregulation of the activators (Figure 4). Furthermore, NOCT a circadian deadenylase that functions as a post-transcriptional modulator in circadian control of metabolic genes was upregulated in individuals from the CO 2 seep (Table S9).
Various transcripts involved in cellular stress response was also DE in our study. In fact, this group had the highest number of DE transcripts majority of which were upregulated in individuals from the CO 2 seep (Figure 4). This included transcripts known to play a role in the production of reactive oxygen species (ROS), transcripts  (Table S9).
Ion channels and transporters were another major functional group that was differentially regulated between fish from F I G U R E 2 Overall transcript expression pattern across the 16 samples from the natural CO 2 seep and control site at Vulcano Island. Principal component analysis (PCA) was performed using the regularized log-transformed (rlog) counts of all expressed transcripts. The circles represent the individuals from the control site, and the triangles represent the individuals from the low pH site. The blue and grey region is the 95% confidence interval of the samples from the CO 2 seep and control sites, respectively. F I G U R E 3 Functional groups associated with the differentially expressed transcripts. Each circle represents one transcript, and the colour represents the FDR corrected p-value (p adj ) from the differential expression analysis. the CO 2 seep and control sites. Notably, transcripts encoding the Na + /K + transporting ATPase (ATP1A1) and carbonic anhydrase 10 (CA10), which play critical roles in ion regulation and maintaining pH homeostasis, were upregulated in fish from the CO 2 seep site ( Figure 3). Additionally, various transcripts involved in actin cytoskeleton organisation were also DE (Table S9 and Figure 4). This

F I G U R E 4
Relative expression levels of transcripts with enriched functions of fish from control and CO 2 seep sites. Each row represents a transcript, and the row z-score is calculated from the log2 transformed transcripts per million (TPM) values for each transcript. The transcripts are clustered based on expression levels. could indirectly affect membrane transport as the cytoskeleton closely interacts with the plasma membrane. Changes in ionic composition across neuronal membranes can affect the activity of neurotransmitters such as GABA and glutamate and we found transcripts involved in degradation of GABA (SLC6A13 and GABAT) and a transcript encoding the metabotropic glutamate receptor 4 (GRM4) to be differentially expressed in our study (Table S9 and Figure 4).
Additionally, transcripts encoding chloride transporters and claudins, whose expression is affected by neural activity, were also DE (Table S9 and Figure 4).
Several transcripts involved in metabolism were also found to be differentially regulated between fish from control and CO 2 seep sites. Transcripts involved in key metabolic processes that result in increased ATP production such as glycolysis, gluconeogenesis, TCA cycle, mitochondrial beta-oxidation, mitochondrial respiratory chain, and lipid and amino acid metabolism were expressed at higher levels in gobies from the CO 2 seep site. Interestingly, transcripts involved in fatty acid synthesis, an ATP-consuming process were downregulated in individuals from the CO 2 seep (Table S9 and Figure 4).

| Outlier loci detection
We evaluated the genetic background of the fish living in the CO 2 seep and control site. A total of 8629 high-quality SNPs were obtained across all 16 individuals from the CO 2 seep and control sites.
There was no genetic structuring found among the analyzed samples and no outlier was detected ( Figure S5). Therefore, there was no evidence of genetic divergence between the individuals from the control and CO 2 seep site. Gobies have a pelagic larval phase and hence the individuals from the CO 2 seep and control sites come from a heterogeneous pool of larvae and there is no evidence of postsettlement processes between the sites.

| DISCUSS ION
Natural CO 2 seeps provide an excellent opportunity to study the effects of long-term exposure of fish to elevated CO 2 in their natural environment and determine their adaptive potential while also accounting for various indirect effects of ocean acidification (OA) such as habitat shifts and changes in species assemblages. Here, we investigated changes in the brain transcriptome of the anemone goby from a natural CO 2 seep in Vulcano Island to identify molecular processes enabling acclimation to elevated CO 2 in the wild. Previous studies at Vulcano Island have found limited behavioural effects of CO 2 on this fish species and increased population size at the seep suggesting that this species has higher tolerance to elevated CO 2 conditions possibly through differential regulation of gene expression enabling physiological and behavioural adjustments (Nagelkerken et al., 2015;Spatafora et al., 2022).
We found large differences in overall brain transcript expression levels between individuals from CO 2 seep and ambient control sites.
The expression patterns were more variable among individuals from the CO 2 seep compared to those from the control site which is not due to the age or length of the fish. Additionally, the SNP analysis did not reveal any post-settlement selection processes which could be driving the increased variability of transcript expression levels in fish from the CO 2 seep compared to those from the control site. Several previous studies have reported variability in species response to elevated CO 2 (Cattano et al., 2018;Monroe et al., 2021;Munday, McCormick, et al., 2013;Paula et al., 2019;Welch & Munday, 2017), which could be due to variation in acclimatization strategies or epigenetic differences. The variation within this population could also result from inter-individual differences in capacity for plasticity when faced with environmental change (Forsman, 2015). Given that developmental plasticity is potentially the primary mechanism driving local adaptation within the CO 2 seep, it is highly likely that the observed variation in transcript expression is resulting from postsettlement plastic responses. Further research is, however, needed to fully understand the underlying causes of individual variation in CO 2 tolerance especially since intraspecific variability can be heritable (Welch & Munday, 2017) and hence is a key factor in influencing adaptive capacity of populations to future global change (Lehmann et al., 2022;Munday, McCormick, et al., 2013).
The majority of differentially expressed (DE) transcripts exhibited elevated expression in fish from the CO 2 seep similar to what has been observed in gonads of the common triplefin from the CO 2 seep in New Zealand, another species known to have increased population size in acidified environments (Petit-Marty et al., 2021).
The gobies as well as the common triplefin seem to be able to elicit acclimation responses to elevated CO 2 by differentially regulating transcripts associated with crucial functions which potentially enables them to maintain cellular homeostasis and thrive in acidified waters at CO 2 seeps showing a more general adaptive capacity of species that thrive in these environments. Specifically, for the gobies, 2.3% of the brain transcriptome involved in key functions such as acid-base and ion homeostasis, neurological function, circadian rhythm, metabolism, cellular stress response and immune response were differentially expressed (DE) between fish from the CO 2 seep and control sites. Cellular Stress Response was the most gene-rich response in gobies at the CO 2 seep. We observed DE of various transcripts associated with redox-sensitive pathways such as NF-kB and MAPK signalling as well as transcripts known to induce an increase in the intracellular concentration of reactive oxygen species (ROS). Therefore, gobies inhabiting the CO 2 seep could experience changes in cellular redox balance resulting in oxidative stress (Poli et al., 2012;Zhang et al., 2016). They, however, also seem to be capable of counteracting stress-induced damage by upregulating transcripts involved in cellular detoxification and redox regulation. Oxidative stress in fish can trigger cell damage and apoptosis (Pizzino et al., 2017;Strader et al., 2020) and also cause conformational changes and oxidative damage to all macromolecules including DNA, proteins and lipids thereby affecting the transmission of genetic information and resulting in intracellular accumulation of unfolded and misfolded proteins (Juan et al., 2021) which is cytotoxic (Hetz, 2012). Differential regulation of transcripts involved in DNA damage response and DNA repair, unfolded protein response (UPR), chaperones such as the heat shock protein (HSP) family observed here suggest that gobies are capable of minimizing oxidative damage on macromolecules caused by ROS. High constitutive expression of several universally conserved genes in the cellular stress response machinery in fish from the CO 2 seep indicates an adaptive response to increase cellular resistance to oxidative stress-induced cell damage.
Transcripts involved in innate and adaptive immune response and the complement system were found to be upregulated in gobies from the CO 2 seep site. Elevated levels of cellular stress can result in increased infections from opportunistic pathogens and hence activation of immune responses could be a protective mechanism to prevent potential infections in the face of increased physiological stress in environments with elevated CO 2 levels (de Souza et al., 2016).
Immune response genes are commonly differentially regulated in fish exposed to elevated CO 2 (de Souza et al., 2014(de Souza et al., , 2016Machado et al., 2020) revealing it as an important function that is regulated under elevated CO 2 conditions. Additionally, two hypoxia-inducible transcripts, HIF1A and EGLN2, were upregulated in individuals from the CO 2 seep. While the primary function of these transcripts is to regulate cellular response to hypoxic conditions, their expression could also be enhanced in response to cellular stress. Increased mRNA expression of HIF1A in response to hypercapnia has been previously reported under stable oxygen levels (Machado et al., 2020).
These results indicate that gobies living in the naturally acidified waters of Vulcano Island are able to elicit physiological compensatory mechanisms to counter the cellular stress induced by prolonged exposure to elevated CO 2 levels.
Living in CO 2 seeps also induced expression differences in various transcripts associated with ion transport and conserving pH balance. In particular, carbonic anhydrase (CA), an enzyme which plays a key role in CO 2 excretion, acid-base balance and ion regulation was upregulated in gobies from the CO 2 seep, which could facilitate increased CO 2 excretion thereby preventing acidosis and also enhance transport of O 2 to the brain (Rummer & Brauner, 2011).
Additionally, transcripts encoding potassium channel transporters (KCNAB3, KCNK10 and KCNJ10), and an Na + /K + -transporting ATPase (ATP1A1) were differentially expressed (DE) in the gobies. Acid-base balance and ion homeostasis in fish have closely linked processes as the transport of H + and HCO 3 − ions are coupled to K + , Na + and Cl − transport (Gilmour & Perry, 2009). Regulation of ion transport is therefore key to ensure proper cellular functioning, maintain pH homeostasis and prevent respiratory acidosis which is a common effect of elevated CO 2 exposure in fish (Heuer & Grosell, 2014;Ishimatsu et al., 2005). DE of transcripts encoding potassium channel proteins in our study could potentially be involved in enabling gobies living in the CO 2 seep to sense ambient CO 2 levels and in turn trigger downstream homeostatic responses to restore internal acid-base status as potassium channels are involved in CO 2 chemoreception and sensing acid-base disturbance (Qin et al., 2010). Additionally, the increased expression of ATP1A1 in fish from the CO 2 seep may play an important role in pH homeostasis as the sodium and chloride gradient created by this enzyme affects the activity of other key secondary transporters involved in pH regulation (Casey et al., 2009;Singh Jaggi et al., 2015;Tresguerres et al., 2020). Transcripts involved in actin cytoskeleton organization were also DE which combined with DE of Ras GTPases, which are involved in remodelling the cytoskeleton, could imply changes in membrane-cytoskeleton interactions in individuals from the CO 2 seep, in turn affecting membrane transport (de Curtis & Meldolesi, 2012). Therefore, gobies inhabiting the naturally acidified waters of Vulcano Island can efficiently regulate various ion transport mechanisms to buffer any intracellular pH changes caused by elevated environmental CO 2 and thereby maintain cellular homeostasis.
We found subtle evidence of regulation of GABAergic signalling in gobies from the CO 2 seep which could be a consequence of changes in ion channel activity as neural signalling pathways are reliant on ion transport and stable ion gradients. Two transcripts, SLC6A13 and GABAT regulating the turnover rate of extracellular γaminobutyric acid (GABA), which is a major neurotransmitter in the brain, were upregulated in gobies from the seep. Given that GABA is not prone to enzymatic degradation, increased expression of these transcripts could be an adaptive response to regulate GABAergic signalling which is known to be altered in high CO 2 environments as a result of pH buffering processes (Heuer et al., 2016;Hübner & Holthoff, 2013;Kang et al., 2022;Nilsson et al., 2012;Schunter et al., 2019). Furthermore, upregulation of the glutamate metabotropic receptor 4 (GRM4) in fish from the CO 2 seep could be a compensatory response to reduce presynaptic GABA release (Bobeck et al., 2014;Flor et al., 1995). Additionally, transcripts involved in the transport of chloride ions, whose concentration gradient influences GABA activity, were DE, which could potentially regulate GABAergic neurotransmission (Jentsch et al., 2002). Particularly, the downregulation of SLC12A1 in the gobies could be an adaptive response to restore any changes in GABA function (Grosell, 2011;Hübner & Holthoff, 2013) similar to what has been observed with the potassium chloride co-transporter 2 (KCC2) gene in the damselfish Acanthochromis polyacanthus exposed to elevated CO 2 (Schunter et al., 2018). Therefore, gobies inhabiting the naturally acidified waters of Vulcano Island possibly use alternative strategies to ensure the proper functioning of neural signalling pathways reliant on the maintenance of stable ionic gradients which are altered by pH buffering processes.
The circadian rhythm is a key pathway regulating nearly all physiological processes and all core transcripts in this pathway showed differential expression. CO 2 -induced changes in the circadian rhythm have been repeatedly observed across various species and are possibly linked to changes in GABA receptor activity (Lee et al., 2021). The circadian pathway enables organisms to anticipate and respond to environmental changes which is evolutionarily advantageous (Ayyar & Sukumaran, 2021;Dmitriev & Mangel, 2000;Prokkola & Nikinmaa, 2018). Changes in the expression of circadian genes in response to elevated CO 2 , have been suggested to result in an amplitude or phase shift of the circadian clock providing an adaptive advantage by increasing flexibility of physiological processes (Kang et al., 2022;Lee et al., 2021;Schunter et al., 2016Schunter et al., , 2021. We found evidence of potential circadian regulation of metabolism in gobies in response to elevated environmental CO 2 . In addition, NOCT, an immediate early gene with deadenylase activity was upregulated in gobies from CO 2 seep. NOCT functions as a mediator of post-transcriptional circadian regulation of metabolic processes and can also respond to changes in the external environment (Stubblefield et al., 2012). Increased expression of NOCT in individuals from the CO 2 seep combined with DE of all core circadian genes suggests rhythmic regulation of energy metabolism, potentially enabling gobies to meet the increased energetic demands of living in an environment with higher than ambient levels of CO 2 . Given the ubiquity of the circadian clock, there is a need for further research to understand the biological implications of the effect of elevated CO 2 on the circadian rhythm and its impacts on neuronal signal transduction and other downstream physiological processes.
In addition to the above-mentioned rhythmic regulation of metabolism, we found an overall increase in the expression of transcripts involved in ATP production in gobies from the CO 2 seep.
Interestingly, transcripts involved in fatty acid biosynthesis, an anabolic process, were downregulated. The constant need to maintain acid-base balance and homeostasis while minimizing oxidative stress in a high CO 2 environment can be energetically expensive (Cattano et al., 2018;Kelly & Hofmann, 2013;Strader et al., 2020).
Decreased expression of transcripts involved in fatty acid synthesis could be an alternative strategy of gobies living in the seep to facilitate allocation of acetyl-CoA, a key metabolite that intersects with glycolysis, TCA cycle, fatty acid beta-oxidation and fatty acid synthesis (Currais et al., 2019), for use in energy-producing metabolic pathways. This in addition to the upregulation of SLC16A1 and SLC25A22, which are also involved in the transport of key metabolic intermediates, suggests that gobies selectively regulate metabolic pathways to ensure increase energy production. Additionally, all transcripts encoding hexose/glucose transporters, with the exception of SLC45A1 were also upregulated in gobies from the CO 2 seep which could facilitate increased availability of glucose in the brain. This is particularly important as glucose is the primary substrate for energy production in the brain. Despite being a glucose transporter, SLC45A1 was selectively downregulated probably because it is a sugar-H + symporter (Shimokawa et al., 2002) and hence facilitates H + import into neuronal cells thereby increasing intracellular pH which could be detrimental. Therefore, gobies have potentially acclimated to a high CO 2 environment by selective regulation of genes involved in metabolism to increased energy production. Changes in the energy budgets of fish living in acidified waters could result in physiological trade-offs (Kelly & Hofmann, 2013). However, the increased metabolic capacity of gobies observed here combined with CO 2 -driven habitat shifts resulting in the enrichment of key trophic resources at the CO 2 seep (Nagelkerken et al., 2015) potentially enables them to meet increased energy requirements. This aids in maintaining pH homeostasis and mitigating any potential cellular damage due to elevated CO 2 levels without a significant trade-off in the energy available for other biological processes. In fact, the CO 2 vents have a higher abundance of gobies compared to nearby control sites (Mirasole et al., 2019;Nagelkerken et al., 2015;Spatafora et al., 2022). At other CO 2 seeps at White Island, New Zealand, a similar habitat shift with net resource enrichment can be found and the common triplefin F. lapillum also exhibit a higher abundance with a similar increase in metabolic capacity (Nagelkerken et al., 2015;Petit-Marty et al., 2021). Therefore, gobies similar to the common triplefin, might benefit from the habitat shifts at CO 2 seeps and show an adaptive response to prolonged CO 2 exposure with an overall increase in metabolism.
Our findings show molecular signatures of acclimation in the anemone goby to the acidified waters of Vulcano Island. The gobies effectively regulate disturbances in acid-base balance, maintain cellular homeostasis and mitigate OA-induced oxidative stress. Gobies also showed potential adaptive mechanisms to maintain the proper functioning of neurological signalling pathways, which could be disrupted as a consequence of altered ionic composition due to pH buffering processes. Although these processes can be energetically expensive potentially resulting in trade-offs, the overall increase in metabolism seen in gobies from the CO 2 seep suggests that they meet the increased energy demands associated with living in an environment with higher than ambient CO 2 concentrations. This could be further facilitated by net trophic resource enrichment at the CO 2 seep. The increased population density of gobies at the CO 2 seep in Vulcano Island reported in previous studies combined with evidence of acclimation to elevated CO 2 environments (potentially mediated by developmental plasticity) found in our study indicates that gobies can effectively mitigate the behavioural effects of elevated CO 2 exposure. Given that a wide range of species show developmental plasticity (Munday, McCormick, et al., 2013;Munday, Warner, et al., 2013) this could be a vital process in decreasing species vulnerability and result in the establishment of sufficient standing phenotypic variation within populations which could become targets for selection to act upon in the future as OA continues to occur at a whole ocean scale.

ACK N O WLE D G E M ENTS
We thank Dr. Juan Diego Gaitan-Espitia and Dr. Kang Jingliang for providing guidance and feedback on the de novo transcriptome assembly.

FU N D I N G I N FO R M ATI O N
This project was funded by start-up funds from The University of Hong Kong to CS.

CO N FLI C T O F I NTE R E S T S TATE M E NT
All authors declare that they have no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
RNA-seq raw sequences and the de novo assembled transcriptome assemblies have been deposited in NCBI under BioProject PRJNA869880.