Rapid evolution fuels transcriptional plasticity to ocean acidification

Abstract Ocean acidification (OA) is postulated to affect the physiology, behavior, and life‐history of marine species, but potential for acclimation or adaptation to elevated pCO2 in wild populations remains largely untested. We measured brain transcriptomes of six coral reef fish species at a natural volcanic CO2 seep and an adjacent control reef in Papua New Guinea. We show that elevated pCO2 induced common molecular responses related to circadian rhythm and immune system but different magnitudes of molecular response across the six species. Notably, elevated transcriptional plasticity was associated with core circadian genes affecting the regulation of intracellular pH and neural activity in Acanthochromis polyacanthus. Gene expression patterns were reversible in this species as evidenced upon reduction of CO2 following a natural storm‐event. Compared with other species, Ac. polyacanthus has a more rapid evolutionary rate and more positively selected genes in key functions under the influence of elevated CO2, thus fueling increased transcriptional plasticity. Our study reveals the basis to variable gene expression changes across species, with some species possessing evolved molecular toolkits to cope with future OA.


| INTRODUC TI ON
Global ocean surface pH is projected to decline at a rate of approximately 0.02 pH units per decade with the ongoing uptake of anthropogenic atmospheric CO 2 by the oceans , a process termed ocean acidification (OA). There is concern as to how marine life will respond and whether adaptation to this rapid acidification is possible. Higher partial pressures of carbon dioxide (pCO 2 ) in seawater cause blood arterial pCO 2 to increase in fish and other marine animals (Brauner et al., 2019). Fishes can prevent acidosis and maintain their acid-base balance at higher pCO 2 levels by compensatory ion exchange (Brauner et al., 2019), but this can alter energy budgets (Lefevre, 2019) and impact other biochemical processes (Brauner et al., 2019). Therefore, although some fish appear to be relatively unaffected by projected future CO 2 levels (Clark et al., 2020;Munday et al., 2009bMunday et al., , 2011Munday et al., , 2019a, a decade of laboratory experiments indicate that predicted OA conditions could affect some marine fishes' physiological performance, growth, survival (Hannan et al., 2020(Hannan et al., , 2021Munday et al., 2019a;Nagelkerken et al., 2021), and behaviors (e.g., olfactory, auditory, vision, behavioral lateralization, learning, activity, anxiety, boldness, foraging behavior, and homing ability) (Goldenberg et al., 2018;Jiahuan et al., 2018;Munday et al., 2019a;Nagelkerken & Munday, 2016;Paula et al., 2019aPaula et al., , 2019bPorteus et al., 2018).
The reported behavioral impairments in elevated pCO 2 conditions have been associated with altered function of GABA A (γ-aminobutyric acid type A) receptors-the major inhibitory neurotransmitter in the vertebrate brain essential for normal brain function, neuronal activity, information processing, and plasticity through maintaining the polarity and amplitude of Cl − fluxes (Bhat et al., 2010;Heubl et al., 2017;Heuer et al., 2016Heuer et al., , 2019Lai et al., 2015;Nilsson et al., 2012;Schunter et al., 2019). This indicates that these impairments have a neural origin (Heuer et al., 2019;Paula et al., 2019a). Furthermore, other physiological responses to elevated pCO 2 , such as ion exchange and acid-base regulation, are also controlled by the brain (Wang et al., 2020). Consequently, changes in gene expression in the brain caused by elevated pCO 2 can be used to understand the mechanisms underpinning neural effects of elevated pCO 2 in marine organisms and the likely impacts on their behavior, physiology, and adaptive capacity.
In fishes tested in laboratory experiments, gene expression patterns indicate that neural responses to elevated pCO 2 vary considerably among species (Hamilton et al., 2017;Lai et al., 2016Lai et al., , 2017. Nonetheless, some recurring expression changes linked to ion transporters, including the GABA receptors (Schunter et al., 2018;Williams et al., 2019), have been identified from the transcriptional architecture of Acanthochromis polyacanthus brains (Schunter et al., 2016(Schunter et al., , 2018 and Oncorhynchus kisutch olfactory tissues (Williams et al., 2019). Notably, Ac. polyacanthus exhibited a shift in expression of core circadian genes (CCGs) after elevated pCO 2 exposure (Schunter et al., 2016). While these controlled experiments provide valuable insights into the effects of elevated pCO 2 on brain function in fish, they lack the interface with other natural factors such as predators, competitors, water currents, and food supply that are present in the wild (Langdon & Atkinson, 2005;Riebesell & Gattuso, 2015). Therefore, the patterns of gene expression observed in laboratory experiments under elevated pCO 2 may not fully capture the response of fish to future ocean acidification (OA) conditions in nature.
Volcanic CO 2 seeps, such as the Upa-Upasina Reef in Milne Bay Province (Normanby Island, Papua New Guinea), can be used as natural laboratories where CO 2 rises from the substratum and acidifies the surrounding seawater to levels similar to, or sometimes beyond, the projections for OA by the end of this century (Fabricius et al., 2011;Hall-Spencer et al., 2008;Munday et al., 2014). These CO 2 seeps provide a unique opportunity to investigate the longerterm effects of OA conditions on fishes in their natural habitat and how this varies among species. While previous studies found that elevated pCO 2 at CO 2 seeps either altered behavior (e.g., reduced predator-avoidance) and physiology in fishes (Munday et al., 2014;Nagelkerken et al., , 2017, or did not affect predator recognition (Cattano et al., 2017), the direct effect on brain function in wild fishes has not been tested. Here, we sequenced the transcriptomes of the brains of adult fish of five damselfish species (Ac. polyacanthus, Amblyglyphidodon curacao, Dascyllus aruanus, Pomacentrus adelus, and Pomacentrus moluccensis) and one cardinalfish species (Ostorhinchus compressus) from a reef within the Upa-Upasina CO 2 seep in Papua New Guinea (pH 7.77, pCO 2 846 µatm (Fabricius et al., 2011)) and an adjacent reef (500 m distance) with ambient pCO 2 (pH 8.01, pCO 2 443 µatm (Fabricius et al., 2011); Figure 1), including an additional collection of Ac. polyacanthus individuals from the CO 2 seep during a storm that flushed ambient pCO 2 water over the CO 2 seeps site reducing the seep pCO 2 . Based on 14,634 orthologous genes across the six fish species, we investigated the effects of natural in situ exposure to elevated pCO 2 by evaluating genetic and brain gene expression variation to detect similarities and disparities in the molecular and functional responses to OA. We further evaluated the evolutionary drivers of adaptive potential in these coral reef fishes to future OA.

| Molecular response across species from the CO 2 seep
The six coral reef fish species exhibited varying degrees of transcriptional changes in response to elevated pCO 2 (Figure 2a). In fish collected at the CO 2 seep, compared with the nearby control reef, Ac. polyacanthus experienced transcriptional changes that were an order of magnitude higher than all other species (1679 differentially expressed genes; DEGs, Table S1), followed by D. aruanus (169 DEGs, Table S1). In contrast, Am. curacao only displayed one DEG (C-X-C chemokine receptor type 4). Furthermore, no significant genetic sequence variations were detected between individuals of each species collected at the CO 2 seep and the adjacent control reef ( Figure  S1a,b). The largest standing genetic variation within species was also observed in Ac. polyacanthus and D. aruanus (i.e., 918 and 923 singlenucleotide polymorphisms [SNPs], respectively; Table S2). The large variations in how the brain of these six coral reef fishes responded to elevated pCO 2 resulted in no genes commonly differentially expressed across all species and a maximum of four DEGs shared among three species. Sixty-five genes exhibited common differential expression between the CO 2 seep and control reef for two or more species (Table S3, Figure S2), which were involved in important biological processes, such as the response to stimulus (e.g., HSPB1, HSP47, SUN1), immune response (e.g., IRF1, IRF4, NFIL3), and circadian rhythm (CR) (e.g., PER1, CLOCK, BMAL1) (Figure 2b), indicating the fundamental responses to elevated pCO 2 across species in a natural setting.
CRs are near-24-h oscillations found in nearly all aspects of physiological processes in the vertebrate brain and body (Dunlap & Loros, 2016;Logan & McClung, 2019), of which 3-10% of transcripts display rhythmical expression (Staels, 2006). These molecular rhythms are generated by a transcriptional-translational feedback loop in CCGs (Logan & McClung, 2019). Environmental perturbations can reset CR in vertebrates, which is characterized by changes in expression of CCGs (Logan & McClung, 2019). Elevated pCO 2 exposure has been shown to induce expression changes in the CR pathway in Ac. polyacanthus (Schunter et al., 2016) as well as an anemonefish, Amphiprion percula, in laboratory settings (Schunter et al., 2021). Three species in this study-D. aruanus, P. moluccensis, and Ac. polyacanthus-displayed significant expression differences in CCGs between the CO 2 seep and control reef ( Figure 2c, Table S4).
However, these CCGs were modulated differently in Ac. polyacanthus, with more downstream changes compared with P. moluccensis and D. aruanus. The difference is supported by gene interaction networks, which quantify the interaction between the CCGs and other DEGs in the GO term "CR" (GO:0007623) and showed no interaction for P. moluccensis and D. aruanus ( Figure S3). Of 12 DEGs (Table S4) that directly interacted with the CCGs in Ac. polyacanthus, some genes regulate the CCGs activity via ubiquitination (UBE3A), phosphorylation (GSK3B), alternative splicing (ANM5), mediating nicotinamide adenine dinucleotide (NAD + ) biosynthesis (NAMPT), coactivating cAMP response element binding protein (CRTC1), or binding to either the promoter (KLF10) or untranslated region (RBM4). The transcriptional changes of CCGs in Ac. polyacanthus, therefore, generate a more significant impact on downstream genes. F I G U R E 1 Schematic illustrating samples collected for six species of coral reef fishes and the molecular and bioinformatic analyses pipeline. Wild adult fishes were collected from both the CO 2 seep (average pCO 2 846 µatm (Fabricius et al., 2011)) and the adjacent control reef (separated by approximately 500 m, average pCO 2 443 µatm (Fabricius et al., 2011)), and brains were dissected. RNA extraction and sequencing protocols were performed, and expressed orthologous genes among all species were identified to examine differential gene expression and potential genetic variances between individuals from the CO 2 seep and control reefs and investigate disparity of protein sequences between species In the main network (1438 genes, 9621 interactions) of all DEGs of Ac. polyacanthus, 632 and 1307 DEGs are connected to the CCGs by two or three nodes, respectively. This emphasizes the importance of CR in the transcriptional adjustments that occur in the brain in response to elevated pCO 2 and downstream genes are either directly or indirectly changed following the expression shift in CCGs.
An additional differential response between fish from the CO 2 seep and control reef that was common across the five species and more significant for three species (i.e., Ac. polyacanthus, D. F I G U R E 2 Molecular responses of six coral reef fish species from the natural CO 2 seep and nearby control reef. (a) number of differentially expressed genes (DEGs) between the CO 2 seep and control reef for each species. The height of each bar represents the number of DEGs in each species. Each color associated to a fish species indicates the specific DEGs for the species, and the black part shows the overlapping DEGs between this and at least one another species. (b) Functional enrichment of DEGs for each species, the size of the circle or triangle indicates the ratio of DEGs per function to total DEGs per species, the triangle indicates DEGs that were overrepresented in the respective function, red to blue coloration ranged from highly significantly overrepresented to not significant. Three more specific functions (T-helper cell differentiation, T cell activation, B cell activation) within immune response were displayed separately to show the extraordinarily high proportion of DEGs related to immune response in O. compressus. (c) Simplified schematic of the molecular transcriptional-translational feedback loops of core circadian DEGs. A red background indicates upregulation and grey indicates downregulation at the CO 2 seep. Red lines indicate positive regulation, and gray lines indicate negative regulation of the downstream genes. (d) The regulation of immune response DEGs in three species under elevated pCO 2 . Ac. polyacanthus is in red, D. aruanus is in black, and O. compressus is in green, with each dot representing one immune response related gene and its upregulation (Log2Foldchange >0) or downregulation (Log2Foldchange <0) at the CO 2 seep with some key genes highlighted by name  were previously observed across various oxygen uptake rate metrics (Hannan et al., 2021). Nonetheless, previous laboratory experiments where fishes were exposed to elevated pCO 2 also indicated a change in the expression of genes involved in the immune response (Machado et al., 2020;de Souza et al., 2014). Hence, immune regulation appears to be an important function involved in the response to elevated pCO 2 across numerous species of reef fishes.

| Distinct regulation of intracellular pH and ion transport in Acanthochromis polyacanthus
A key challenge for fishes when coping with elevated pCO 2 is maintaining blood and tissue pH homeostasis to avoid acidosis. Ac.
polyacanthus displayed 25 DEGs (Table S6) Table S6). These core pH i regulation genes are involved in the extrusion of intracellular H + or uptake of extracellular HCO 3 − during acid-base regulation in many fishes (e.g., rainbow trout, dogfish, zebrafish (Brauner et al., 2019)). In addition, CAH1 and CAHZ are precursors to, or directly related to carbonic anhydrase (CA), the enzyme that catalyzes the rapid conversion of CO 2 to H + and HCO 3 − and the reverse. CAH1 was upregulated in both Ac. polyacanthus and D. aruanus from the CO 2 seep. This pattern has been found in other fishes, albeit upon exposure to extremely high CO 2 levels (Georgalis et al., 2006;Perry et al., 2010). However, the CA-dependent mechanism may be in place during a generalized acidosis-one that would ensue with exposure to elevated environmental CO 2 -to maintain or even enhance O 2 delivery to tissues (Rummer & Brauner, 2011;Rummer et al., 2013).
Furthermore, cardiac β1-adrenergic receptors (ADRB1), which are part of the stress response and could help relieve impairments from potential hypercapnic tachycardia (Miller et al., 2014), were significantly expressed in Ac. polyacanthus from the CO 2 seep. Likewise, hemoglobin (HBA, HBA1, HBB), which is the primary protein responsible for O 2 transport (Rummer & Brauner, 2015;Rummer et al., 2013), was also significantly expressed at elevated levels in Ac. polyacanthus from the CO 2 seep. Because a generalized acidosis reduces the affinity and carrying capacity of hemoglobin for O 2 , increased expression of various hemoglobin genes as well as those related with CA, as mentioned above, could counter the effect of a generalized acidosis such that O 2 uptake and delivery are safeguarded (Harter & Brauner, 2017;Munday et al., 2019b;Rummer & Brauner, 2011).
What is interesting is that these responses related to acid-base balance and O 2 delivery were observed exclusively in Ac. polyacanthus. As such, Ac. polyacanthus is sensitive to elevated pCO 2, with regulations to pH i that may provide this species advantages to cope with elevated pCO 2 . In fact, while other pathways found here have been previously reported in Ac. polyacanthus in controlled laboratory studies (Schunter et al., 2016(Schunter et al., , 2018, the adjustments related to intracellular pH have not previously been shown, which reveals a key mechanism when fish are exposed to elevated pCO 2 in the wild. Neural signal transductions are mediated through ion channels and pumps to passively or actively push ions in and out of cells (Gadsby, 2009). Ac. polyacanthus sampled from the CO 2 seep exhibited significant expression changes in a larger number of ion transporters compared with individuals from the control reef ( Figure 3b; Figure S4, Table S7), a pattern not found for other species (Figure 3b; Figure S4). Interestingly, Ac. polyacanthus from the CO 2 seep displayed overall repressed expression for Ca 2+ /K + transporters (Table   S7). Of these transporters, voltage-dependent calcium channels and KCC2 are involved in the GABAergic pathway. This pattern, along with the attenuated expression of adenylyl cyclase (ADCY2), one GABA A receptor (GBRB3), and two GABA B receptors (GABR1, GABR2), implicates that Ac. polyacanthus displayed a repression of the GABAergic pathway when exposed to elevated pCO 2 (Table S8, Figure S5). Besides, elevated pCO 2 induced Ac. polyacanthus with a decreased expression in genes related to neurotransmitter signals, such as five evolutionarily conserved protein families (RIMS, UN13b, RIMB2, LIPA2, RB6I2) (Sudhof, 2012) and two glutamate receptors (NMDE4, NMDE3A) (Brassai et al., 2015). Hence, Ac. polyacanthus, when faced with long-term elevated pCO 2 in the wild, use a large repertoire of neuronal regulators that potentially affect neural signal transduction in synapses.
The ion and neuronal regulatory activity in Ac. polyacanthus might display rhythmic expression patterns driven by the CCGs.

| Plasticity in transcriptional changes of Acanthochromis polyacanthus at the CO 2 seep
Flexibility in gene expression adjustments for Ac. polyacanthus was observed in the response to short-term exposure to ambient pCO 2 levels at the seep in our study. While our samples were collected under stable environmental conditions (i.e., flat sea, little to no wind), a storm event lasting around 24 h provided the opportunistic collection of Ac. polyacanthus from the CO 2 seep (hereafter, "storm" individuals). The wind during storm events, which by itself can decrease surface ocean pCO 2 by 150-200 µatm (Massaro et al., 2012), can flush the oceanic ambient pCO 2 waters into the CO 2 seep, leading to a dilution of the high pCO 2 waters. Hence, this storm event allowed us to also investigate the effects of a short-term reduction in pCO 2 on fish that have been exposed to chronic high pCO 2 in the wild. Storm individuals showed a similar expression pattern to control reef individuals (97 DEGs; Table S1), but larger differences compared with CO 2 seep individuals (3212 DEGs, Table S1; Figure   S6). A common 1399 genes (Table S3) were differentially expressed between control reef and storm individuals compared with CO 2 seep individuals ( Figure S6). Moreover, if some genes of individuals from the CO 2 seep exhibited upregulation, then the storm individuals had these genes downregulated to a lower expression in comparison with control and vice versa (Figure 4a), suggesting an "overcompensated" pattern.80.3% of these common DEGs, such as the CCGs ( Figure S7), revealed this "over-compensated" pattern, which was also indicated from the increased number of DEGs enriched in a wide range of functions during the storm ( Figure S8), such as the CR, pH i regulation, ion transport, immune responses, and GABAergic pathways (Tables S4-S8; Figure 4b). The GABAergic pathway showed decreased expression in individuals from the CO 2 seep compared with individuals from the control reef, but elevated expression during the storm ( Figure S5). Hence, with high pCO 2 water from the seep becoming mixed and replaced by ambient pCO 2 oceanic water during the storm for 24 h, molecular adjustments in Ac. polyacanthus exposed to chronically elevated pCO 2 were counteracted or even over-compensated after acute pCO 2 reduction. This species has been found in previous studies to be potentially adapted to fluctuating CO 2 , as they display better swimming performance and a greater aerobic scope than individuals from stable elevated pCO 2 conditions (Hannan et al., 2020), together reflecting that Ac. polyacanthus can flexibly cope with variations pCO 2 . In addition, Ac. polyacanthus has been shown to exhibit larger molecular responses on short-term exposure to elevated pCO 2 (4 days) compared with chronic exposure (Schunter et al., 2018;Tsang et al., 2020). As such, the present study F I G U R E 3 Regulation of pH i , Ca 2+ , and K + transport in response to elevated pCO 2 in four coral reef fish species. Genes with significant upregulation at the CO 2 seep have a red background and downregulation is indicated in blue. (a) Core acid-extruder and uptake genes involved in acid-base homeostasis with differential expression for the different species. (b) Overall downregulation of calcium and potassium transporters across the cell membrane in Ac. polyacanthus Na + also reveals that acute changes in CO 2 levels can cause larger transcriptional changes in the brain compared with chronic exposure to elevated pCO 2 in the wild.

| Disparity between species: Evolutionary rate and positive selection
Changes in gene expression are expected to drive evolutionary differences across species (King & Wilson, 1975;Wray, 2007), yielding divergent acclimation or adaptation potential to elevated pCO 2 . By applying Expression Variance and Evolution (EVE) models (Rohlfs & Nielsen, 2015) to the phylogeny of the six species ( Figure S9), 183 genes were detected to have significantly diverged expression patterns among all species (Table S9). Of these, ATP7A and APBA2, two genes that are both involved in neural signal transduction, showed notable divergences between Ac. polyacanthus and the other species ( Figure S10). Evolutionary differences are also partly driven by protein sequence divergences across species (King & Wilson, 1975;Wray, 2007). Hence, based on the phylogenetic tree ( Figure S9), we used nonsynonymous and synonymous ratios (d N /d S ) as indicators of protein sequence divergences to estimate gene evolutionary rate and positive selection for each species. Ac. polyacanthus exhibited a more rapid evolutionary rate than other species (Figure 5a, Table S10). Moreover, Ac. polyacanthus has a more rapid evolutionary rate also in genes related to the CR (GO:0007623), regulation of pH i (GO:0051453), and ion transmembrane transport (GO:0034220; Figure 5b, Table   F I G U R E 4 Over-compensated expression patterns to acute environmental changes during a storm in Acanthochromis polyacanthus individuals at the CO 2 seep. (a) Two trends of expression changes are detected in 1399 common differentially expressed genes (DEGs) between storm versus CO 2 seep individuals, and control versus CO 2 seep individuals. For illustration purposes, a subset (right: 52 genes; left: 51genes) of these common DEGs with normalized reads counts <= 200 were selected. Gene expression increased (left) or decreased (right) from control reef to CO 2 seep but over-compensated from CO 2 seep to storm by a downregulation (left) or upregulation (right). (b) The number of DEGs increased from CO 2 seep to storm in core circadian rhythm genes, ion transporters, and those related to pH i regulation, immune response and the GABAergic pathway

(b)
Total Core CR genes Immune response pHi regulation Ion transporters GABAergic pathway S11), which also exhibit large expression differences between the CO 2 seep and control site. This species, therefore, displays a genetic adaptive system with evolved expression patterns in key functions to cope with elevated pCO 2 .
Positive selection is an important source of evolutionary innovation and a major force behind the divergence of species (Kosiol et al., 2008) and can, therefore, reveal adaptive potential. By the comparisons between protein sequences of all six species in this study, Ac. polyacanthus displayed more positively selected genes (PSGs) (1251 PSGs, Table S12) (Figure 5d, Figure S12). The most notable gene in this network was cyclic AMP-responsive element-binding protein 1 (CREB1) with four sequential nonsynonymous mutations at the binding domain ( Figure S13). CREB1 is a well-known transcription factor that drives the transcription and translation of clock genes (Colwell, 2011;O'Neill et al., 2008), which interact with DEGs including CCGs (BMAL1, CRY1, PER1, NFIL3), an ion transporter gene (KCMA1), and GABAergic pathway genes (ADCY2, ADCY9; Figure S12). Due to this, the circadian output of Ac. polyacanthus exposed to pCO 2 is pro-  In fact, Ac. polyacanthus is known to have strong population genetic structure and even possible divergence within the species (Miller-Sims et al., 2008;Planeset al., 2001) possibly owing to the direct development of juveniles and lack of a larval dispersal phase.
While we can only hypothesize, it is possible that this life-history trait leaves Ac. polyacanthus with enhanced potential to adapt to local conditions more than species with a dispersive larval phase and as a species therefore experiences elevated evolutionary rates and positive selection, which in turn provides this species with more transcriptional plasticity to respond to environmental changes. In addition to low evolutionary rate, the other fish species with higher dispersal may cope with environmental changes more frequently than Ac. polyacanthus, therefore these species are not as sensitive to environmental changes as Ac. polyacanthus and exhibited lower transcriptional plasticity in response to elevated pCO 2 .
Environmental conditions, including pCO 2 levels, can change rapidly in nature, which is not often recreated under laboratory conditions. Here, we show that a tropical fish species, Ac. polyacanthus, displayed a means to react to elevated CO 2, as well as variability in environmental CO 2 at a natural CO 2 seep, through high transcrip-

| Sampling and sample processing
Adult coral reef fishes were collected at a shallow coral reef The median pH and pCO 2 of seawater at seeps were 7.77 units and 843 µatm, comparing with a higher pH (8.01 units) and lower pCO 2 (443 µatm) at control sites. There are no significant difference of the temperature and salinity between the CO 2 seeps and control sites (Fabricius et al., 2011). This study was conducted in accordance with James Cook University animal ethics permit A2534. Six selected fish species (130 individuals) were sampled across multiple days, with each species being collected from the CO 2 seep and control reefs at the same time of day (Table S14). The six species were selected due to their abundance at the control reef and CO 2 seep and include five damselfishes, Ac. polyacanthus (control: 11, CO 2 seep: 7; storm: 9), Am. curacao (control: 9, CO 2 seep: 10), D. aruanus (control: 10, CO 2 seep: 10), Pomacentrus adelus (control: 11, CO 2 seep: 11),

| Sequence processing and transcriptome assembly
To only work with high-quality sequences we assessed the raw read quality with FastQC v0.11.8 by removing adapter sequences and low-quality sections and reads with  (Table S15). To address potential contaminates from the sequence reads we removed potential archaea, bacteria, fungi, and viral reads with Kraken v2.0.7-beta (Wood & Salzberg, 2014) (--confidence 0.7), which filtered out around 5 million reads (  et al., 2016) were applied to assess their completeness and quality scores, which indicated that they include almost single copy orthologues in the Actinopterygii database (Table S16), and good contig ratio in orthologous genes were much higher (82.9%-93.8%) than de novo assemblies and ORFs (Table S17).

| Orthologous gene detection and annotation
To be able to compare the molecular response across the six species, we detected orthologous genes based on the ORF protein sequences using the default parameters in OrthoFinder v2.3.3 (Emms & Kelly, 2015). 116,503 initial orthologous groups among the six species were then filtered to contain at least one transcript for each species (remaining 17,564 orthogroups). We further filtered out groups with the smallest E-value was kept. Thirdly, For BLAST hits per species that were assigned to the same protein, we selected the one with highest number of mapped reads for each species. Finally, we removed the orthogroups with blast hits to the same protein in all six species if their average E-value of the six hits was larger than 1E-30.
When more than one protein was detected in the same orthogroup, the one with largest total bit score was retained. As a result, 14,635 orthologous groups (Table S18) that include only one representative transcript per species were chosen as the final set of orthologous genes for the further analysis. The final set of orthogroups were then annotated using the UniProt Blast results that represent each orthologous group and GO terms and InterPro pathways were added with OmicsBox (Conesa et al., 2005).

| Gene expression analyses
To identify gene expression differences between samples collected at CO 2 seep and control reefs, the high-quality reads were mapped against the reference transcript for each species with Bowtie 2 v2.3.5.1 (Langmead et al., 2009)  compressus with name including "B cell," "T cell," "T-helper" (Table   S6) as immune-related functions, all DEGs of each species underlying these functions were extracted as the DEGs related to immune functions. And the genes with name including "calcium", "sodium", "potassium", and "chloride" were regarded as the preliminary ion transport-related genes.
To understand gene expression shifts in the six species through a phylogenetic perspective, the analysis of EVE model (Rohlfs & Nielsen, 2015), which models gene expression as a quantitative trait across a phylogeny to estimate expression variance within versus among species, was performed based on a maximum likelihood (ML) tree of the six species. To prepare the ML tree for EVE model, the species-specific protein sequence in each orthologous group were aligned, respectively, using Clustal Omega v1.2.4 (Sievers et al., 2011) (-t Protein), and the gaps within alignments were removed by Gblocks v0.91b (Talavera & Castresana, 2007) (-t=p). The orthologous group would be excluded when any sequences was shorter than 50 codons, which left 14,456 of the total orthologous groups. To construct the ML tree RAxML-NG v0.9.0git (Kozlov et al., 2019) was run based on concatenate sequences (7,096,195 base pairs, 83.61% invariant sites) of the orthologous groups. The best ML tree was selected based on 50 random and 50 parsimony-based starting tree topologies and the LG +G4m model (LG substitution matrix with discrete GAMMA N categories, ML estimate of alpha) (Bernal et al., 2020 as the lineage-specific divergence genes (higher differences between lineages) and plastic genes (higher differences in expression within their own lineage than between lineages).

| SNP calling, population structure, and outlier loci detection
To identify whether the potential genetic variances exist between samples from CO 2 seep and control reefs, SNP calling was performed among all six species. Samples were individually mapped using Bowtie2 against the references (as described above). The reads per sample were sorted, assigned to new read-group, marked duplicate, and reordered by PICARD in the GENOME ANALYSIS TOOLKIT (GATK v3.8-1-0) (DePristo et al., 2011), which was used for SNP calling and SNP filtering ( Figure S1c). After splitting the reads that contain Ns in their cigar string, the SNPs within samples were detected using HaplotypeCaller in GATK. At first, the genotypes were filtered based on the Phred-scaled quality score using a threshold of 30. All resulting SNPs were merged together by "-T CombineVariants -genotypeMergeOptions UNIQUIFY." Next, potential SNPs were removed if there were more than three SNPs in 35 base pair window size. We also filtered out potential SNPs if any of the following con-  (Berdan et al., 2015) were retained. To compare the sequence variance between individuals from CO 2 seep and control site, VCFtools (v0.1.13) was used to keep biallelic SNPs with "PASS" flag and remove SNPs with low F ST (0 or "non"), low MAF (< 0.1, to exclude recent mutations that are uninformative loci) (Roesti et al., 2012) and high heterozygosity (large than 0.8, to remove potential paralogs) (Berdan et al., 2015).
Once the high-quality SNPs were obtained between samples from CO 2 seep and control site, ADMIXTURE v1.3.0 (Alexander & Lange, 2011) was performed to estimate individual ancestry from the SNPs data sets by simulated K (the priori number of groups) from 1 to 10.
Putative outlier loci were detected using a Bayesian approach implemented in BayeScan 2.0 (Foll & Gaggiotti, 2008) under two different assumptions that the samples from CO 2 seep and control site were one population or two different populations.

| Estimation of the evolutionary rate and positively selected genes
The ML tree was also applied to estimate the evolutionary rate of our six species. To estimate lineage-specific evolutionary rates and positively selected genes (PSGs), the nucleotide sequence of each orthologous gene per species were aligned to the representative protein by GeneWise (Birney & Durbin, 2000) (Du et al., 2020). The MCMC process was run for 5,000,000 steps and sampled every 5000 steps. MCMCtree suggested that the divergence time between Ac. polyacanthus and Am. curacao in our study is approximately 24.4 Mya ( Figure   S16), which is similar to previous estimates (McCord et al., 2021).
Secondly, to construct the phylogenetic tree of all individuals in our study, we assembled the sequences of 50 randomly selected genes from 2686 single-copy orthologues referring the methods in Kuang et al. (2018). We kept the longest ORFs estimated by TransDecoder v3.0.1 (Haas & Papanicolaou, 2018) for each assembled gene, and the concatenated ORFs sequences (45,343 bp) of the final 31 single-copy orthologs (Table S19) that were captured by all 130 individuals were aligned by MUSCLE v.3.8.31 (Edgar, 2004)

| Gene interaction network analysis
To further understand the interaction, and relative importance of DEGs and PSGs across six coral fish species, the gene interaction network was estimated using Cytoscape v3.7.0 (Shannon et al., 2003) to draw edges between genes with functional interactions reported for humans in the STRING database (Szklarczyk et al., 2016). Gene interaction network is selected if its nodes are more than three.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.