The time course of molecular acclimation to seawater in a euryhaline fish

The Arabian pupfish, Aphanius dispar, is a euryhaline fish inhabiting both inland nearly-freshwater desert ponds and highly saline Red Sea coastal lagoons of the Arabian Peninsula. Desert ponds and coastal lagoons, located respectively upstream and at the mouths of dry riverbeds (“wadies”), have been found to potentially become connected during periods of intense rainfall, which could allow the fish to migrate between these different habitats. Flash floods would therefore flush Arabian pupfish out to sea, requiring a rapid acclimation to a greater than 40 ppt change in salinity. To investigate the molecular pathways of salinity acclimation during such events, a Red Sea coastal lagoon and a desert pond population were sampled, with the latter exposed to a rapid increase in water salinity. Changes in branchial gene expression were investigated via genome-wide transcriptome measurements over time from 6 h to 21 days. The two natural populations displayed basal differences in genes related to ion transport, osmoregulation and immune system functions. These mechanisms were also differentially regulated in seawater transferred fish, revealing their crucial role in long-term adaptation. Other processes were only transiently activated shortly after the salinity exposure, including cellular stress response mechanisms, such as molecular chaperone synthesis and apoptosis. Tissue remodelling processes were also identified as transient, but took place later in the timeline, suggesting their importance to long-term acclimation as they likely equip the fish with lasting adaptations to their new environment. The alterations in branchial functional pathways displayed by Arabian pupfish in response to salinity increases are diverse. These reveal a large toolkit of molecular processes important for adaptation to hyperosmolarity that allow for successful colonization to a wide variety of different habitats.

Salinity is one of the main abiotic factors shaping the distribution and habitat preference of fish and other aquatic taxa. Teleosts have evolved different strategies to maintain osmotic homeostasis depending on water ion concentrations. In seawater, fish drink copiously to combat water lost via osmosis, concurrently suffering passive intake of high quantities of salts. Specialized mitochondria-rich cells of the gill epithelium, called ionocytes, are then responsible for active excess ion secretion 1,2 . In freshwater, fish face ion loss and passive osmotic water intake. To compensate for the consequent dilution of their body fluids, they actively uptake ions from the surrounding medium through specialized ionocytes 1,2 . Given the profound differences in osmoregulatory mechanisms in fresh versus seawater, changes in salinity represent a significant challenge for fish, and the induced osmotic stress can interfere with physiological homeostasis and impairment of biological processes 3 . As a consequence, most teleost fishes are restricted to limited habitat salinities (stenohaline fish). Other species, termed euryhaline fish, have evolved osmoregulatory plasticity in the organs involved in maintenance of osmotic balance 4 . In particular plastic modifications to the gill epithelium, the main tissue responsible for osmoregulation in fish 2 , grant the ability to live in wider salinity ranges and exploit larger habitat diversity. www.nature.com/scientificreports/ The fish were transported to the King Abdullah University of Science and Technology (KAUST) Coastal and Marine Resources Core Lab and maintained in closed-system tanks. The conditions in the aquaria resembled those of the collection sites with 8.32 ± 0.03 pH, 10L:14D photoperiod, and salinity at 1.92 ppt ± 0.01 SE for the desert pond fish and 42.68 ppt ± 0.09 SE for the coastal lagoon fish (Suppl. Table 2). Salinity was achieved using seawater, diluted for the desert pond fish with dechlorinated tap water. In order to avoid potential confounding effects due to temperature differences, tank water was kept at 26 °C, an intermediate value between the two collection site ones (31.6 °C desert pond; 20.9 °C coastal lagoon). Fish were acclimated to such holding conditions for more than eight months before the salinity change experiment. All fish were kept at low densities, with maximum seven individuals per 80 L tank, in six replicate tanks per treatment, and fed once a day ad libitum with commercial pelleted feed. In October 2016, the osmotic challenge experiment started: fish from the desert pond population were directly transferred to seawater, in order to mimic the change in salinity they might experience during potential colonization of coastal lagoons of the Red Sea (Fig. 2). Temperature, photoperiod and pH were kept constant during the entire duration of the experiment in order to avoid any confounding effects. No mortality was recorded throughout the whole experiment. Fish were sampled pretransfer (0 h) from both populations (1.9 ppt and 43 ppt), and at five post-transfer time points (43 ppt; 6 h, 24 h, 72 h, 7 days, 21 days), always at the same time of the day to avoid potential confounding effects (except for the 6 h time point). At each sampling event, fish from different replicate tank were rapidly euthanized using an overdose of tricaine methanesulphonate  www.nature.com/scientificreports/ (MS-222, MP Biomedicals), sexed, and length was measured (Suppl.  58 . Moreover, to evaluate the completeness of the assembly, and to control for potential loss of core genes during the filtering process, BUSCO v3.02 59 was run after every filtering step, using the provided Actinopterygii set of 4584 Benchmarking Universal Single-Copy Orthologs. Annotation, differential gene expression and gene ontology analysis. Transcriptome annotation was performed first by BLAST searches of the TransDecoder predicted ORFs and the untranslated transcripts, using blastp and blastx algorithms respectively (evalue ≤ 1e −5 ), against the UniProtKb and the NCBI non-redundant (nr) databases. Where conflicts were found, the following order of priority was observed: UniProtKb/Swiss-Prot, NCBI nr, and UniProtKb/TrEMBL. BLAST results were then loaded in OmicsBox v1.4.11 60 , where Gene Ontology (GO), InterPro, and EggNOG 61 annotations were additionally performed using Blast2GO 62 .
To test for gene differential expression (DE) between the two source populations, as well as along the salinity acclimation time series, reads from each sample were first quantified using Salmon v1.1.0 63 in mapping-based mode against the de novo assembled transcriptome. Transcript abundance estimates were then summarized at gene level and imported in DESeq2 v1.26.0 64 using the package tximport v1.14.2 65 in R v3.6.1. To visually check for possible overall patterns, outliers and batch effects, raw counts were transformed using the vst function from DESeq2 package, and principal component analyses (PCAs) were run on the 500 most variable genes across samples using plotPCA function of the same package. Ellipses were drawn using ggplot2 stat_ellipse function in R. Based on PCA results and Cook's distances calculated by the DESeq function, three samples (DP_Rep1, 24h_Rep4, 21d_Rep5) were labelled as outliers and excluded from the DE analysis. Differentially expressed genes (DEGs) were identified from the remaining 30 samples (3-6 samples/time-point) running pairwise comparisons of control populations and post-transfer time points using the contrast function of DESeq2 with shrunken log2 Fold Change (log2FC) estimates by apeglm 66 . The threshold cut-offs chosen to identify significant DEGs were |log2FC| ≥ 1, False Discovery Rate (FDR) adjusted p-value 67 < 0.05 (Wald test), and a mean expression of > 10 reads (baseMean), as similarly done in other studies on osmotic responses in fish [68][69][70] .
To identify groups of genes showing comparable trends, such as rapid or longer-term responses to the salinity challenge, DEGs revealing similar expression patterns across time points were clustered. Additionally, Impul-seDE2 v1.10.0 71 , a Bioconductor R package specifically designed for time series data, was employed in case-only mode to discern steadily increasing or decreasing expression trajectories from transiently up-or down-regulated genes (FDR adjusted p-value 67 < 0.05). For all clustering purposes, Red Sea population samples were treated as if they belonged to an additional time-point, the last in the acclimation timeline (long-term acclimation beyond 3 weeks).
Functional enrichment analyses were performed in OmicsBox v1.4.11 60 , with the Fisher's Exact Test (FDR < 0.05) after removing duplicated annotations from the differentially expressed genes (DEGs) and identified gene cluster sets to find over-represented GO terms between controls and post-transfer time points.

Results
Transcriptome assembly and annotation. The first de novo transcriptome assembly for Aphanius dispar was created from a total of 1.38 billion raw reads with an average of 24.5 million reads per sample after trimming and error correction steps (Suppl. Table 3). These high-quality reads were used to de novo assemble the Arabian pupfish gill transcriptome, which resulted in 650,824 contigs with a 97.2% overall mapping rate, an N50 of 1275 bp, and 86.2% BUSCO completeness score using the Actinopterygii database (Suppl. Gene expression differences in natural populations. The two sampled populations, desert pond (DP) and the highly saline Red Sea (RS) coastal lagoon, exhibited a clear branchial gene expression separation even after months of acclimation to the same holding conditions, salinity excluded (time 0; Fig. 3). Such differences in gene expression were partly conserved throughout the salinity challenge (Suppl. Fig. 1). The pairwise comparison of fish in their original salinities resulted in 552 differentially expressed genes (DEGs; Suppl. Table 5), which is the largest expression difference of the experiment (Fig. 4). Several functions were differentially regulated in the two populations, in particular ion transport and immune system. Full names for gene symbols can be found in Suppl. Table 6. Ion transport related terms represented 70% of total enriched GO categories between the populations at time 0 (Suppl. Table 7). The RS individuals in particular upregulated specific seawater ion transporters, such as the Na + /K + /Cl − transporter (SLC12A2), the cystic fibrosis transmembrane conductance regulator (CFTR), and the transient receptor potential (TRP) cation channel subfamily V member 1 (TRPV1), together with genes . Several Na + /K + -transporting ATPase (NKA) subunits (ATP1A1, ATP1A3) were also differentially expressed between the two populations, and elevated expression in the nearly-freshwater DP fish was found for a different set of ion channels, including members of the transient receptor potential (TRP) cation channel superfamily (TRPM2, TRPM5, TRPM7, TRPV4), chloride channels like the chloride channel protein 2 (CLCN2), the inward rectifier potassium channel 2 (KCNJ2), and the intracellular channel inositol 1,4,5-trisphosphate receptor type 1 (ITPR1), together with its modulator, the calcium binding protein 1 (CABP1). Aquaporin 3 (AQP3) was also upregulated in DP samples. The two populations also showed differential regulation of the immune system. C-C motif chemokine ligands 3 and 25 (CCL3, CCL25) were upregulated in the RS population, while members of the C-X-C motif chemokine family (CXCL8, CXCL11.6, CXCL14) were upregulated in the DP fish, enriching the "chemokine activity" function. Moreover, C-C chemokine receptors (CCR4, CCR6) were upregulated in the gills of the DP population, while several genes coding for components of acquired immunity were upregulated in the seawater RS individuals, like immunoglobulins and major histocompatibility complex proteins (IGKC, IGHM, IGL1, IGKV4-1, MR1, H2-EA).
Time course of the molecular responses to the salinity challenge. The selected experimental timeline and the clustering analyses of differentially expressed genes along the acclimation window ( Fig. 4) allowed for the discovery of distinctively timed processes in the Arabian pupfish branchial response to the abrupt increase in salinity. While some DEGs were only transiently differentially regulated along the acclimation timeline (Fig. 5a), 231 and 269 genes steadily increased or decreased, respectively ( Fig. 5b; Suppl. Table 8), and in particular the downregulated genes exhibited enrichment in cell cycle related terms (Suppl. Table 9). Through the investigation of each time point, mechanisms typical of a short-term stress response were uncovered at 6 and 24 h post-transfer, while the 72 h and 7 day time points revealed cell cycle arrest and tissue remodelling events, and the last time point (21 days post-transfer) was characterized by longer-term acclimation processes, many resembling the Red Sea population transcriptional profile (Fig. 6). www.nature.com/scientificreports/ Just six hours after the start of the hyperosmotic challenge, 75 branchial genes showed a statistically significant change in expression compared to pre-transfer controls (Suppl. Table 10). This immediate reaction was based on a wide array of different functions (Suppl. Table 11). The acute osmotic stress caused the onset of cell signalling cascades, with the differential expression of prolactin receptor (PRLR), the stress responsive hydroxysteroid 11-beta dehydrogenase 2 (HSD11B2), which modulates intracellular cortisol levels, and serum/glucocorticoid regulated kinase 1 (SGK1), important in the cellular stress and DNA damage responses. Furthermore, the osmotic stress response was found to entail ion homeostasis pathways, through the regulation of ion transport and ion channel activity (CFTR, KCNJ2, TCAF, WNK2), and organic osmolyte synthesis and transport (GLUL, IMPA1, ISYNA1-B, SLC6A20, SLC5A7, ABCG20). In particular, ISYNA1-B was transiently upregulated between 6 and 24 h and, in the same interval, IMPA1 showed very high levels of expression, while being upregulated across the whole timeline (Suppl. Fig. 2a). Another important function found at 6 h post-transfer was tissue modification, by means of cellular proliferation and differentiation (GPM6B, NR4A3, NHSL1), and keratinization (CNFN, DSC1). Genes involved in the immune response were also differentially expressed at 6 h, and class I histocompatibility antigen, F10 alpha chain (HA1F) was upregulated post-transfer. Upregulation of genes implicated in lipid (ALOXE3, ALOX15B, CPT1A) and glucose (GAD1, SLC2A8) metabolism and transport was also found. Although a few DEGs were found to be related to circadian rhythm processes (NFIL3, NR1D1, NR1D2, PER2), this result is most likely due to the sampling time occurring at a different moment of the day, rather than an effect of the salinity change. For the same reason, the possibility that other DEGs in this time point might be responding to the difference in time of day rather than salinity cannot be excluded.  Table 12). "Gland development" was among the enriched GO terms (Suppl . Table 13), and included genes related to secretion, such as sodium and chloride channels (CFTR, CLCN2, SLC12A2), cell proliferation (E2F7, FGF7, FOXM1), hormonal response to stress actors (CRHR1, PRLR), cytoskeleton and extra-cellular matrix organization genes (DIAPH3, SOX9). "Lipid metabolic process" was also enriched at 24 h, with genes related to cell membrane glycosphingolipid and glycerophospholipid biosynthesis (B4GALNT1, PLAAT4) and genes involved in lipid transport, such as carnitine palmitoyltransferase 1A (CPT1A), known to be part of the mitochondrial fatty acid oxidation pathway, and already upregulated at 6 h. The increased energetic demand was additionally manifested by the upregulation of NADH:ubiquinone oxidoreductase complex assembly factor 4 (NDUFAF4), and solute carrier family 24 member 48 (SLC25A48), among others. Several ion transport actors and osmotic stress response regulators still played a role at 24 h as did various transcription factors, repressors, and regulators (Suppl. Table 11). Some genes involved in the DNA damage response pathway also still showed differential expression. Notably, several genes involved in keratinization (CNFN, K1C1, KRT13, S100A11, TGM5) were upregulated compared to time 0, while there was a downregulation of genes involved in the extracellular matrix degradation (ADAMTS5, COL23A1, EFEMP2, PLOD2, THSD4). Tissue remodelling processes happening in transferred fish were additionally reflected by the differential expression of genes involved in the regulation of apoptosis (Suppl.  www.nature.com/scientificreports/ C (CYC), solute carrier family 25 member 6 (SLC25A6) and voltage-dependent anion channel 2 (VDAC2), all thought to be involved in the mitochondrial apoptotic pathway. Clustering revealed that most of the above-mentioned genes exhibited a transiently up-or down-regulated pattern in the first 24 h, being differentially expressed compared to pre-transfer fish only in these early time points (Fig. 5a; Suppl. Fig. 2; Suppl. Table 14). The main differences between the 6 and 24 h post-transfer time points (109 DEGs; Suppl. Table 15), as portrayed by the many related enriched GO terms (Suppl.  Table 18), and the strongest signal of this contrast was the downregulation of more than 30 genes (Suppl. Table 11) involved in mitotic cell cycle and cell population proliferation (Suppl . Table 19), including genes fundamental for G1/S and G2/M phase transitions of the cell cycle. Concurrently, genes involved in DNA replication and repair showed downregulation compared to time 0 (Suppl. Table 11). Branchial tissue modifications were still happening at 72 h with DEGs involved in extracellular matrix degradation including: Ca 2+ -activated cysteine protease calpains (CAPN2, CAPN5, CAPN8), and matrix metallopeptidase 8 (MMP8), known to play a role in the breakdown of the extracellular matrix during tissue remodelling, as well as several glycoproteins with cell adhesion functions, such as CEA cell adhesion molecules (CEACAM5, CEACAM6). Genes belonging to the glycosaminoglycan-or mucopolysaccharides-biosynthesis pathway were upregulated at 72 h compared to pre-transfer conditions (Suppl. Table 11).
After 7 days from the start of the hyperosmotic challenge there were 423 DEGs compared to pre-transfer fish at time 0 ( Fig. 4; Suppl. Table 20). This observed high number of differentially expressed genes could be an effect of the small sample size for this time-point (n = 3). Numerous GO terms were enriched (Suppl . Table 21), and over 80 genes, mostly downregulated compared to time 0, pertained to mitotic cell cycle as well as cell proliferation and differentiation processes (Suppl. Table 11). Seven-day post-transfer fish gills were still showing significant changes in the expression of several ion channels, as well as genes implicated in energy production and tissue modifications, such as cell adhesion, cytoskeleton organization, mucopolysaccharide metabolism, and extracellular matrix degradation (Suppl. Table 11).
The last time point at 21 days post transfer exhibited 139 DEGs (Suppl. Table 22) in comparison with pretransfer fish. There was no functional enrichment for these genes, which were involved in a variety of processes (Suppl . Table 11). Osmoregulation was still under transition, with the upregulation of two typical seawater sodium-coupled transporters, the ATPase Na + /K + transporting subunit alpha 1 (ATP1A1) and the solute carrier family 12 member 2 (SLC12A2). While cell cycle and energy production genes were not differentially expressed at 21 days, gene expression related genes and tissue organization were upregulated compared to time 0 (Suppl .  Table 11). Moreover, two genes involved in the biosynthesis of mucopolysaccharides (GALT-1, ST3GAL1) showed upregulation, joining B3GNT7, CHST3, ST3GAL2 that were already upregulated earlier in the acclimation.
Expression signatures leading to adaptation. Additional clustering analyses focused on the identification of genes putatively critical to long-term acclimation to elevated salinity in the Arabian pupfish. Starting from a subset of genes differentially expressed between the two populations at time 0, the focus was placed on 69 DEGs which progressively resembled the expression levels of the Red Sea population along the experimental timeline and are agreed on by the clustering and ImpulseDE2 (Suppl . Table 23). Overall, some functions appeared earlier than others along the acclimation window (Fig. 6). For example, many genes implicated in osmotic stress response and osmoregulation (CA4, CFTR, CLCN2, IMPA1, KCNJ1, STK39, USP2, WNK2) were among the first to steadily change expression, between 6 and 24 h from the start of the exposure. Cornifelin (CNFN), a component of the cornified cell envelope, and glycoprotein M6B (GPM6B), a key upstream regulator of genes involved in actin cytoskeleton organization, were upregulated at every time point post-transfer. Similarly, genes involved in the metabolism of polyamines (ODC1, PAOX), important in hyposaline acclimation in teleosts, were downregulated from the first 24 h of exposure. As time progressed more genes changed their expression to match that of the seawater population. At 72 h post-transfer, the functions of these genes were related to changes in tissue organization and permeability, with the upregulation of calpain 2 (CAPN2) and microfibril associated protein 4 (MFAP4), both part of the extracellular matrix degradation pathway, and the downregulation of the water channel aquaporin 3 (AQP3). Moreover, some genes involved in mucopolysaccharide metabolism (B3GNT7, CHST3, GALT-1) started to be upregulated at 72 h, while others (ST3GAL1, ST3GAL2) showed upregulation after 7 days. From 7 days of exposure to sea water, the transferred fish started to change the expression of branchial genes involved in lipid metabolism, showing upregulation of neutral cholesterol ester hydrolase 1 (NCEH1) and downregulation of low-density lipoprotein receptor (LDLR), both involved in lipoprotein and cholesterol metabolism. The transient receptor potential cation channel subfamily V member 4 (TRPV4), a non-selective cation channel involved in osmotic pressure regulation, was also downregulated from 7 day post-transfer. Finally, at the end of the experimental acclimation period, 21 days post-transfer, there was a downregulation of glycine decarboxylase (GLDC), involved in osmotic regulation, and arginase 1 (ARG1), another component of the polyamine metabolism pathway.

Discussion
Arabian pupfish from isolated nearly-freshwater ponds (1.5 ppt) have been found to be genetically connected to individuals found in highly saline Red Sea (43 ppt) lagoons downstream of dry riverbeds (called "wadies") 39 . One of the possible explanations for these findings is for desert pond pupfish to be sporadically flushed out to the sea through the wadies by flash floods, where they would need to rapidly acclimate to the new environment www.nature.com/scientificreports/ conditions, salinity in particular. This study time course transplantation experiment tested and demonstrated the potential of Arabian pupfish to survive and acclimate to such an abrupt salinity increase. Moreover, comparing the gill gene expression profiles of the native populations, as well as the changes across a salinity challenge from hours to weeks post-transfer, enabled the separation between short-and longer-term osmotic stress responses, and the investigation of the processes that would allow pupfish to adapt to the high salinity typical of the Red Sea environment. The native nearly-freshwater pupfish were able to acclimate to the abrupt increase in water salinity by means of expression changes in a large number of branchial genes. A subset of genes whose expression changed with the salinity exposure to resemble those of native seawater individuals revealed the importance of related functions in long-term acclimation.
The two populations at time 0 displayed the largest difference in gene expression of the whole experiment. Although partially expected because of the different genetic origin and native environment, the magnitude of this result was still surprising given the months spent in the same acclimation conditions, salinity excepted. The main source of expression divergence between the two native populations was osmoregulation, as well as the primary function at play during pupfish acclimation to high salinity. When in seawater, teleosts need to extrude passively accumulated salts through a variety of branchial ion transporters 2 . Accordingly, the Arabian pupfish Red Sea population upregulated ion transporters and osmoregulation related genes typically found in marine teleosts 9,72 . The expression of many of these genes was quickly upregulated following high salinity exposure in desert pond fish and maintained along the entire acclimation timeline. For instance, genes involved in chloride ion secretion in marine-type ionocytes, such as cystic fibrosis transmembrane conductance regulator (CFTR) and Na + /K + /Cl − transporter (NKCC1 or SLC12A2) 9,73 , were upregulated in transferred pupfish already from 6 and 24 h, respectively. Accordingly, in another pupfish species, Cyprinodon nevadensis amargosae, CFTR and NKCC1 gill mRNA levels increased within a similar time frame post-transfer into seawater and remained elevated throughout the 2 week experiment 74 , and upregulation of these two genes was also found in the gills of mummichog (Fundulus heteroclitus) within 24 h post-transfer from brackish to seawater 5 . Involved in the activation of NKCC1 are two other genes, WNK lysine deficient protein kinase 2 (WNK2) and serine/threonine kinase 39 (STK39) 73,75,76 , also differentially expressed across the timeline in Arabian pupfish. WNK2 has indeed been found to regulate NKCC1 activation in Xenopus oocytes 77 , and its upregulation at every time point post-transfer to hyperosmotic water was also reported in the gills of the butterfish Pampus argenteus 78 . Likewise, Flemmer et al. 75 reported the upregulation of STK39 in seawater-acclimated mummichog gills, and additionally linked it to the increased expression of NKCC1 in the same samples. The prompt upregulation of marine-type ion transporters and osmoregulatory genes in desert pond pupfish exposed to seawater to quickly resemble the Red Sea population levels confirms the importance of these genes in osmoregulation to highly saline conditions and reveals the high degree of conservation in osmoregulatory mechanisms among different fish species, Arabian pupfish included.
Another fundamental component of teleost osmoregulation in salt water is the myo-inositol biosynthesis (MIB) pathway, through which the compatible organic osmolyte myo-inositol is synthetized and accumulated inside cells during osmotic stress for protection from salinity-induced damage 79 . The two enzymes constituting the MIB pathway, inositol monophosphatase 1 (IMPA1) and myo-inositol phosphate synthase (MIPS or ISYNA1), were both differentially expressed in high-salinity exposed pupfish. Accordingly, both genes were also upregulated following seawater transfer in other euryhaline fishes, such as eels, Anguilla anguilla 80 , and turbots, Scophthalmus maximus 81 , where the MIB pathway knockdown was directly implicated in causing weakened gill osmoregulation and reduced survival 82 . The MIB pathway is therefore an important osmoregulation mechanism in seawater across several different species. However, the transient upregulation of MIPS in the first 24 h only highlights the importance of this pathway in the short-term osmotic stress response of the Arabian Pupfish.
In hyposaline waters fish are susceptible to passive ion loss and need to compensate by active uptake of osmolytes from the surrounding water 1 . In Arabian pupfish, chloride uptake is likely accomplished by chloride ion channel protein 2 (CLCN2), typically found in freshwater ionocytes 83,84 , which was highly expressed in desert pond individuals pre-transfer, and decreased from 24 h of high salinity exposure onward. Rainbow trout (Oncorhynchus mykiss) ionocytes 84 and a Sacramento splittail, Pogonichthys macrolepidotus, population exposed to seawater 20 also show this pattern, indicating the importance of this chloride channel in a variety of fish species inhabiting freshwater environments, where chloride uptake from the surrounding water is a key aspect of osmoregulation. Conversely, desert pond Arabian pupfish lack most of the previously described mechanisms for the uptake of sodium 8,10 , and might possibly exclusively rely on specialized isoforms of Na + /K + -ATPase (NKA) subunits for its import. Although the NKA gene is usually identified in marine acclimated teleosts, where it is part of the ion secretion machinery, in Arabian pupfish different transcripts annotated to the α-1 subunit gene were upregulated in desert pond individuals, and only one transcript was upregulated in Red Sea pupfish. While this could be an indication of population-specific NKA subunit α-1 isoforms 20 , it is likely that salinity-dependent isoforms with opposite functions of ion uptake in hyposaline water and salt secretion in seawater exist in this species, as previously described for several euryhaline teleosts [85][86][87][88][89][90] . In the climbing perch (Anabas testudines) and in salmonids, NKA α-1 isoform a expression levels are highest in freshwater and decrease post-transfer to seawater, while other isoform mRNA expressions increase following exposure 87,91 . Accordingly, some of the NKA α-1 subunit coding transcripts showed a steady decreasing pattern along the exposure timeline in transferred pupfish, with one transcript resulting in a 4.5-fold downregulation at the end of the experimental timeline, which could be an indication of acclimation to the seawater habitat, where sodium uptake is not needed anymore. Revealed Arabian pupfish hyperosmoregulatory mechanisms could represent a new model for fish branchial ion absorption at low salinities and confirm the wide diversity of evolutionarily distinct branchial adaptations to hyposaline environments in teleosts. Overall, the main mechanisms differentiating Arabian pupfish natural populations pertain to osmoregulation, and rapid acclimation from near-freshwater to highly saline waters involves adjusting the gill gene expression to resemble the osmoregulatory processes typical to the long-term adapted seawater population. www.nature.com/scientificreports/ The second major difference between the two Arabian pupfish populations concerned the immune system, with several genes involved in the inflammatory and immune responses also differentially regulated post-transfer. In teleosts, salinity has been known to have intricate impacts on the immune system 34,92,93 . While osmotic stress has been found to increase the nonspecific immune response, a depression of the acquired immune response has also been reported owing to trade-offs in resource allocation 34 . Red Sea population and translocated desert pond fish showed however overexpression of acquired immune response components, such as immunoglobulins (Ig) and major histocompatibility complex class I-related (MR1). The immune response capacities are therefore likely not impacted during acclimation to seawater in Arabian pupfish, but rather they might actually be suppressed in near-freshwater, as opposed to what usually happens in salmonids and other fish species 34 . Another possible explanation for the elicitation of the Arabian pupfish adaptive immune response in seawater could lie in differences in pathogen attributes between salinities of the experimental setup. Gills indeed represent a major site for pathogen entry and consequent defence mechanism elicitation 94 , and salinity has been known to influence several pathogen characteristics, such as distribution 95,96 and virulence 97,98 . Regardless of the reason, other fishes have been shown to not suffer from immune depression in seawater, such as Nile tilapia, Oreochromis niloticus 99 , as well as Acanthopagrus latus and Lates calcarifer, where plasma Ig levels increase with water salinity 100 . Likewise, hyperosmotic immersion of Paralichthys olivaceus boosts branchial major histocompatibility complex expression and the overall mucosal immune response 24-48 h post-exposure 101 . Indeed, a crucial role in immunity and defence in teleosts is exerted by mucosal surfaces 102 . Fish gills, skin and gut are coated with a thin mucus layer which acts as a barrier from the surrounding environment and is characterized by physical and antimicrobial defensive functions 94,103 , but a role in osmoregulation has also been suggested 104,105 . Two mucin-like transcripts were upregulated in Red Sea Arabian pupfish, as well as at 24 h and 7 days post-transfer in desert pond fish. Moreover, several genes related to O-linked glycosylation of mucins, glycosaminoglycan metabolism and mucus production 106 were upregulated both in Red Sea individuals and in seawater-exposed desert pond fish. Accordingly, in Anguilla japonica mucosal tissues, seawater elicits an increase in mucus cell numbers and secretion, possibly to trap sodium ions 105 . Similar findings were reported for Salmo salar, in addition to salinity-driven modifications of mucin biochemistry 107 , that were also later described in other euryhaline fishes 108,109 . As the molecular results suggest, increased gill mucus production might be at play in seawater in Arabian pupfish and could be another aspect of their acclimation strategy to hyperosmotic environments.
While osmoregulation and immune response were revealed to be important in long-term adaptation to seawater, a series of short-term and transient response mechanisms were also elicited during the acclimation timeline. Abrupt increases in ion concentration can lead to macromolecular damage in exposed fish epithelia. An increase in sodium, for example, has been linked to cell membrane damage 23 through lipid peroxidation, catalysed by lipoxygenases in response to stress 28 . Two lipoxygenases (ALOXE3, ALOXE15B) were indeed upregulated in Arabian pupfish gills 6 h post-transfer, and were also identified in similar seawater transfer experiments on eels 110 . Such membrane disruption can lead to the activation of the so-called cellular stress response (CSR), which encompasses defence mechanisms to protect and repair damaged cellular components and restore homeostasis 23,28 . The CSR machinery responds to sodium-destabilized proteins by overexpressing molecular chaperones, such as heat shock proteins 23 , two of which were temporarily upregulated in Arabian pupfish at 24 h, as well as in a hyperosmotic challenged Sacramento splittail population 20 . Rising intracellular sodium levels can also cause nucleic acid structural disruptions 23 , and consequently, increased expressions of DNA damage related genes are often reported as part of the CSR following hyperosmotic stress in teleosts 26,27,70 . In Arabian pupfish several DNA damage related genes were transiently upregulated especially in the first 24 h of exposure, like serum/glucocorticoid regulated kinase 1 (SGK1), and similarly increased in other fish following acute seawater challenges 12,111 . Hence, in the first hours of high salinity challenge, gill acclimation in Arabian pupfish is dominated by the onset of macromolecular damage followed by the cellular stress response machinery initiating the repair of the compromised processes to restore homeostasis.
Another aspect of hyperosmotic-induced CSR, at least in cultured human cells, is the inhibition of transcriptional and translational activities 24 . Equivalent to other species, such as the climbing perch 31 , Arabian pupfish exhibited a downregulation of transcription related genes following seawater exposure, which might be a mechanism to prevent the replication of high salinity-damaged DNA. An inhibition in transcription and translation might also explain the onset of cell cycle arrest that was identified in both Arabian pupfish and climbing perch 31 via a downregulation of large sets of cell cycle and mitosis involved genes during the acclimation to seawater. In support of these findings, an immunocytochemistry study in tilapia (Oreochromis mossambicus) observed a G2 phase arrest in the mitotic cycle of gill cells over a period of 16-72 h post-seawater exposure 112 . Days to weeks after the start of the hyperosmotic challenge, an inhibition of transcriptional activity and simultaneous cell cycle arrest might represent a strategy for the fish to prevent the replication of damaged DNA, preserve energy and buy time to respond to macromolecular damage caused by the increase in ion concentration.
For longer-term acclimation to seawater, euryhaline fish gills must undergo profound remodelling events in order to switch from an ion absorbing epithelium to an ion secreting one. Arabian pupfish started displaying processes involved in tissue remodelling from the first hours of salinity exposure. Likely to allow a rapid reorganization of the gill epithelium to reverse the ion transport direction, a transient increase in cell proliferation and differentiation related genes was uncovered, as seen in the gills of euryhaline tilapia in the first 8 h of seawater exposure 112 . At the same time, genes involved in keratinization, a process by which keratin accumulate inside epithelial tissue cells to provide barrier-like functions, started to be upregulated. Keratinization gene expression has been found to be salinity dependent in tilapia 113 and to be upregulated following air exposure in the skin of the mangrove rivulus, Kryptolebias marmoratus 114 . Keratinization may represent a strategy to reduce the amount of water loss during dehydration, possibly also following increased environmental salinity, as seen in Arabian pupfish. Abrupt seawater transfer therefore elicited the rapid activation of gill tissue remodelling pathways, such www.nature.com/scientificreports/ as cell proliferation and keratinization, already in the first 6 h of exposure. Such a prompt response to salinity increases is one of the most striking evidence of the extensive branchial plasticity in this pupfish. Programmed cell death, or apoptosis, of high salinity-damaged cells and freshwater-type ionocytes represents another of the first steps in gill epithelium remodelling, essential for full acclimation to seawater 23 . A transient upregulation of cytochrome c (CYC) and other genes involved in the mitochondrial apoptotic pathway was seen in Arabian pupfish at 24 h. This has been previously recorded in mummichog in response to osmotic stresses 19 , and is supported by a microscopy study in Mozambique tilapia revealing increased branchial apoptotic freshwater ionocytes one day after transfer to seawater 115 . At 24 h post-transfer, genes related to cell adhesion began to be upregulated, and cytoskeleton and extracellular matrix organization functions showed increased expression from 72 h of exposure. Cell adhesion and extracellular matrix pathway upregulation was similarly reported in Sacramento splittails between one and seven days into the acclimation to elevated salinity 20,21 , and branchial cell cytoskeleton reorganization is largely recognized as a fundamental aspect of salinity acclimation in teleosts 12,13,32 . In the first days following seawater exposure, Arabian pupfish gill remodelling hence continued through apoptosis and modifications of cell adhesion, cytoskeleton and extracellular matrix, in order for the transferred fish branchial epithelium to better perform in the changed environment.
As a consequence of tissue remodelling events, an upregulation of mitochondrial respiratory chain and metabolism related genes is also expected to support the increased energy demand, as previously found in similar experiments of euryhaline fish translocation 19,22,31 . Analogously, in Arabian pupfish there was an overall upregulation of metabolism related genes up to 7 days post-transfer, while genes involved in mitochondrial respiration were overexpressed especially between days one and seven, which is consistent with the time frame for major gill remodelling processes in other species 16,17,20 . In a similar fashion to other euryhaline teleosts, seawater exposed pupfish are therefore affected by transient and longer-lasting gill tissue modifications occurring from the first hours to several days after the beginning of the exposure, and potentially resulting in durable modifications which allow longer-term acclimation to the highly saline environment.
Arabian pupfish inhabit profoundly divergent environments of the Arabian Peninsula, ranging from nearlyfreshwater ponds found in desert areas to highly saline Red Sea coastal lagoons. The plasticity of these fish under steep increases in water salinity is hypothesized to play a major role in the colonization potential of this species 39 . By simulating the exposure to high salinity from near-freshwater, this study tested Arabian pupfish potential to survive flash flood events and to move between highly different salinity environments. Not only key processes for a successful acclimation were identified, but also the importance of their timing was uncovered. Arabian pupfish branchial salinity-elicited pathways revealed osmoregulation, immune system and mucus production to be rapid but also long-term acclimation mechanisms to the new environment. In the short-term, cellular stress response processes were triggered, which prevented the fish from suffering permanent damage following acute hyperosmotic exposure. Later in the acclimation, pathways involved in gill epithelium modification and remodelling equipped the organism with lasting adaptations to the increased salinity. While some of the processes occurring during the acclimation timeline resembled mechanisms of seawater exposure previously reported in other euryhaline fish species, others, such as increased mucus production and keratinization, represent less common strategies for high salinity acclimation in teleosts. Overall, the branchial processes revealed in this nearly-freshwater Arabian pupfish population during high salinity acclimation sheds light into this nonmodel euryhaline species colonization potential of seawater habitats. A large set of differentially timed molecular mechanisms plays a role in the plastic reorganization of the gills in hyperosmotic environments that allows for the expansion of euryhaline teleosts into a wide variety of different habitats.

Data availability
Raw sequence data are available through the National Center for Biotechnology Information Sequence Read Archive under BioProject PRJNA722804.