Field application of de novo transcriptomic analysis to evaluate the effects of sublethal freshwater salinization on Gasterosteus aculeatus in urban streams

Freshwater salinization poses global challenges for aquatic organisms inhabiting urban streams, impacting their physiology and ecology. However, current salinization research predominantly focuses on mortality endpoints in limited model species, overlooking the sublethal effects on a broader spectrum of organisms and the exploration of adaptive mechanisms and pathways under natural field conditions. To address these gaps, we conducted high-throughput sequencing transcriptomic analysis on the gill tissue of the euryhaline fish Gasterosteus aculeatus, investigating its molecular response to salinity stress in the highly urbanized river Boye, Germany. We found that in stream sections with sublethal concentrations of chloride costly osmoregulatory systems were activated, evidenced by the differential expression of genes related to osmoregulation. Our enrichment analysis revealed differentially expressed genes (DEGs) related to transmembrane transport and regulation of transport and other osmoregulation pathways, which aligns with the crucial role of these pathways in maintaining biological homeostasis. Notably, we identified candidate genes involved in increased osmoregulatory activity under salinity stress, including those responsible for moving ions across membranes: ion channels, ion pumps, and ion transporters. Particularly, genes from the solute carrier family SLC, aquaporin AQP1, chloride channel CLC7, ATP-binding cassette transporter ABCE1, and ATPases member ATAD2 exhibited prominent differential expression. These findings provide insights into the potential molecular mechanisms underlying the adaptive response of euryhaline fish to salinity stress and have implications for their conservation and management in the face of freshwater salinization.


Introduction
Urbanization, as reported by the United Nations Population Division (UNPD), is poised to expand urban land cover by 1.2-1.8 million Km 2 by 2030 [1].The negative impacts of this global urban expansion on freshwater habitats and biodiversity have led to the coining of the term "Urban Stream Syndrome" to describe the myriad impacts on streams, including altered channel morphology, highly variable hydrographs, diminished biotic richness, and elevated concentrations of nutrients and contaminants [2][3][4].One particularly striking consequence of urbanization is the accumulation of various ions, such as Na + , Ca 2+ , Mg 2+ , K+, NO 3 -, SO 4 2 -, Br, and Cl -, contributing to freshwater salinization-an increasingly pressing issue drawing attention from researchers and policymakers [5][6][7][8].A variety of factors can lead to heightened levels of ion accumulations, including background concentrations in watersheds, as well as human-driven changes in land use, sea level rise, agriculture, road salts, wastewater, and mining [8].Given this background, the Ruhrgebiet region in Western Germany offers an excellent case study for examining freshwater salinization in urban streams.This area is one of the most densely populated in Europe [9] and encompasses a wide variety of activities that are known to contribute to a higher vulnerability of freshwater systems to salinization.Mining is one of the major threats in this region, with coal extraction requiring the infiltration of groundwater that is subsequently pumped into nearby rivers and streams [10].As a result, mining water discharge has been recorded to reach chloride concentrations as high as 3500 mg/L in the River Lippe which is 175 times higher than the typical natural background chloride concentrations in the world surface freshwater ecosystems of <20 mg/L [11].However, the implementation of new management measures and the closure of many mines has led to a significant reduction in concentration, with current levels now below 400 mg/L [10].Despite this progress, the ongoing inputs from abandoned coal mines in Western Germany continue to contribute to freshwater salinization in the region [12].Another significant factor contributing to freshwater salinization in the Ruhrgebiet region is the changing land cover.As urbanization has intensified, rivers and wetlands have become increasingly surrounded by impervious man-made surfaces.This has led to a greater runoff of road salts and sewage waters, particularly in urban surface waters which turns them particularly vulnerable to salinization [13].
Most freshwater animals are strictly adapted to low ion concentrations in their surrounding water compared to their internal ion concentrations (stenohaline), and spend a lot of energy for osmoregulation by actively pumping ions inside their bodies and excreting water [14].Thus, hyperosmotic stress due to freshwater salinization poses a significant physiological challenge to many freshwater animals, including fish which must maintain their ion homeostasis in the face of ion loss through diffusion and water absorption across their permeable membranes [15].Active membrane transport of ions against concentration gradients is required to maintain internal salt concentrations, and increased environmental salinity above the internal ion concentrations leads to greater energy expenditure for osmoregulation [16].As a result of the advancement in fish osmoregulation physiology research, and their importance in bioindication, fisheries, and the economy, fish have become essential models for studying the effects of environmental change [17].Chloride is widely recognized as the most important anion for osmoregulation in freshwater, making it one of the most commonly used indicators to measure salinity effects on fish [18].The annual upper chloride mean threshold concentrations in German surface waters, as per the European water directive, range between 40-90 mg/ L, with Canada and the U.S. having chronic thresholds of 160 mg/L and 230 mg/L, respectively [11,19].However, these thresholds have been based on ecotoxicological tests on a limited number of species that often rely on measuring the salinity effect on endpoint survival and may overlook sublethal effects on non-model species.In ecotoxicology, sublethal effects are defined as those that do not directly cause the death of an individual but rather have an effect on individual fitness (i.e.behavior, cognition, physiology) [20,21].This is critical, as it is recognized that one of the major challenges of understanding the impact of environmental change is determining the point at which sublethal responses to stressors begin to adversely impact the organism [22].Focusing solely on tissue, organ, or endpoint mortality effects and using a limited number of model species may overlook critical information regarding the physiological impacts of sublethal salinity stress on organisms.Moreover, until recently, studies on fish physiology have relied on technologies with reduced resolution and a limited range of species models [18].Thus, there is a pressing need to apply novel technologies and widen the model species pool to gain a better understanding of the effects of salinity stress on freshwater species.
In this context, the application of molecular techniques has greatly improved our understanding of the toxicity mechanisms of major anthropogenic ions in various fish species at high resolution.This development aligns with the recent scientific consensus on the analysis of the response of fish to environmental stressors at different levels, including molecular gene expression and potential adaptation [23].Furthermore, the molecular adaptation to salinity by freshwater species has been identified as one of the major gaps in the research agenda of global freshwater salinization, and transcriptomic technologies have been proposed as a promising approach to tackle it [6].Transcriptomics allow capturing a snapshot in time of the total gene expression (mRNA) present in a cell under a given condition [24].And as noted by Connon et al., (2018) [25], can reveal sublethal stress thresholds that extend beyond the negative impacts on individual and population fitness.This provides a more comprehensive understanding of species' habitat requirements, which in turn facilitates their management.Recently, transcriptomics has been applied to identify the function of novel and conserved genes and the underlying pathways of salinity stress in fish species such as Lateolabrax maculatus, Poecilia reticulate, Pseudopleuronectes yokohamae, Lates calcarifer, Acipenser baeri, Oreochromis mossambicus female × O. urolepis hornorum male, Oncorhynchus keta, and Gasterosteus aculeatus [26][27][28][29][30][31][32][33].While the studies previously cited have provided insight into the mechanisms underlying tolerance thresholds and adaptations to osmotic stress, they have been limited by their reliance on laboratory manipulations and exposure to salinity concentrations that exceed sublethal thresholds.As highlighted by Whitehead & Crawford (2006) [34], it is crucial to investigate expression variation in ecological contexts to fully comprehend the role of gene expression in the adaptive response to environmental stressors.Under natural conditions, organisms are exposed to a variety of fluctuating environmental signals, providing a valuable opportunity to discern gene expression patterns that are not discernible under laboratory conditions [35].
The threespine stickleback (Gasterosteus aculeatus) is an opportunistic species that has demonstrated an exceptional capacity for adaptation to diverse aquatic habitats and is widely distributed in the Ruhrgebiet region in Western Germany.One key trait that has enabled this species to thrive in various environments and colonize systems throughout the Northern Hemisphere is its ability to regulate its osmotic balance in response to changes in salinity (Euryhalinity) [36].Even landlocked freshwater populations of stickleback, which have been separated from the marine environment for thousands of years, possess the osmoregulatory machinery needed to handle abrupt changes in salinity [37].This suggests that sticklebacks have maintained their physiological adaptation to saltwater despite repeated colonization of freshwater habitats.Their wide distribution and adaptability have made them an important model in the study of adaptive evolution and have been subject to the development of largescale genetic and genomic resources [38].Furthermore, experimental evidence has shown that their gill transcriptome responds to subtle changes in environmental salinities within the freshwater range [28].Although their euryhalinity suggests they can tolerate the salinity changes, it has been noted that the activation of the hyperosmotic osmoregulatory mechanisms comes at a high energetical cost [39].Cellular remodeling, enzyme expression, and transport protein synthesis linked to salinity changes have been estimated to account for 20-68% of the total energy expenditure in certain fish species [40].These costs may drive the need for physiological and behavioral optimization to minimize energy expenditure, resulting in narrow physiological tolerance windows and behavioral avoidance of salinities outside the optimal range [39,41].However, there is still a lack of knowledge on the molecular mechanisms underlying this species' responses to salinity stress.Therefore, the threespine stickleback provides an excellent model for studying the mechanisms underlying the osmoregulation of a euryhaline opportunistic species under sublethal salinity concentrations.
Taking into consideration the research gaps and opportunities introduced above, we designed a study with the objective of understanding the physiological effects on G. aculeatus under subtle changes in concentrations of chloride from anthropogenic sources in the field.We are shifting our focus from controlled manipulative studies to the application of transcriptomics in environmental field research.Specifically, we aim to understand how the salinity gradient in the urban stream affects the gene expression profile in the gills of G. aculeatus, with a particular focus on pathways related to osmoregulation.We hypothesize that the salinity gradient in the urban stream will affect the gene expression profile of G. aculeatus gills, with an enrichment of osmoregulation-related pathways.To test this hypothesis, we sampled fish from an urban stream in North Western Germany with a historical salinization influence at three stations along a salinity gradient.We then evaluated the transcriptomic fingerprint of G. aculeatus gills, performed a comparative transcriptomic study, evaluated the gene ontology of the differentially expressed genes and identified candidate genes related to osmoregulation.Ultimately, our goal was to utilize de novo transcriptomic analysis to gain insights into the impact of sublethal freshwater salinization on this euryhaline species in urban streams.This research not only aims to uncover molecular adaptation mechanisms but also to contribute in broadening the range of model species used in the study of sublethal effects of salinity stress on freshwater organisms.

Sampling design
We choose three stations in a gradient of chloride concentration in two tributaries and in the main channel of the river Boye (Fig 1).The station with the highest salinity of the stations is in the Kirchschemmsbach, a stream that has historically affectation by salinity intrusion (Station C).Followed by the main channel of the Boye in its confluence with the Kirchschemmsbach (Station A).Station B is a station in a stream close to the other stations but with no historical or present influence of secondary salinization.It was therefore used ascontrol site.The chloride concentrations in the period of sampling for station B were 12.9 mg/L, Station A 26.4 mg/ L, and station C 32.0 mg/L.In the sites that we sampled, the concentration was below the annual chloride mean upper threshold from German streams in the central highlands which is in the range of 40-50 mg/L [12].And below the reported chronic or acute toxicity level for species such as trout and fathead minnow [42].However, both stations A and C are above the typical Cl − concentrations from natural sources in surface freshwater ecosystems that are in general below 20 mg/L [11].
Gasterosteus aculeatus (Linnaeus) individuals were sampled using electrofishing equipment.Fishing and electrofishing permits were given by the Untere Fischereibeho ¨rde Kreis Recklinghausen (Genehmigung_Recklinghausen_20210824) and the Emschergenossenschaft (Genehmigung_EGLV_Entnahme_20210901114640736).Fish were then quickly euthanized using an overdose of tricaine methanesulfonate (MS222) before harvesting the gill tissue, immediately fixing it in RNAprotect (QIAGEN), and storing it at -20˚C.In total nine samples were obtained, three from the control (Station B), three from the mid-salinity station (Station A), and three from the high-salinity station (Station C).

RNA isolation and Illumina sequencing
The fixed tissue of each sample was homogenized in a solution of buffer 700μl RLT with a 10μl:1ml concentration of β-mercaptoethanol in a FastPrep-24™ (MP Biomedicals) bead beater during 30s at 5 m/s in a screw cap vial with 1mm diameter glass beads.Afterwards, the RNA was extracted using a RNeasy Mini Kit (QIAGEN) following the manufacturer's instructions.The quality of the extractions was checked using a Nanodrop 1000 Spectrophotometer (Peqlab Biotechnologie, Germany), ensuring that the concentration was >50ng/μl, and that OD 260/ 280 = 1.8-2.1, and OD 260/280>1.5.Next, the RIN e value quality was assessed using the RNA ScreenTape system in an Agilent 2200 TapeStation, ensuring that the samples had a RIN e score over 7.0.An Illumina Tru-Seq™ RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) was applied in generating sequencing libraries following the manufacturer's protocol.After purifying, Illumina sequencing was carried out on an Illumina/Hiseq-2500 platform (Illumina, San Diego, CA, USA) to generate 100 bp paired-end reads with 50 million reads depth by the Cologne Center for Genomics (CCG), Germany.The choice of read depth followed CCG technical advice, which estimated that given a published genome size of approximately 200,000 genes, this read depth would be adequate for constructing a de novo transcriptome with nine biological replicates.

Sequence quality control, de novo assembly, and annotation
Raw sequences of all nine samples were checked in FASTQC [43] for quality, followed by a quality improvement via an adaptor and low-quality reads filtering using TRIMMOMATIC, with a threshold of phred + 33 quality score of at least 25 for all bases, and a minimum length of 50 bp for all reads kept for downstream analyses [44].TRINITY v2.9.1 [45] was used to create a de novo assembly using the nine pairs of clean sequences.The reads of the raw nine paired sequences were aligned to the de novo transcriptome using SALMON [46] as the abundance estimation method.BUSCO v5.2.2 [47] was used to assess transcriptome assembly completeness by searching against the actinopterygii_odb10 (Creation date: 2021-02-19) dataset.The expression matrix was generated and the low expressions transcripts (minimum expression of 1.0) were filtered using TRINITY v2.9.1, the expression levels were normalized as transcripts per million transcripts (TPM).Identification of likely protein-coding regions in transcripts was performed with TRANSDECODER v5.5.0 [48].The filtered transcriptome sequencing reads were aligned to protein, signal peptides, and transmembrane domains using the tools, DIAMOND v2.0.8 [49], SIGNALP 6.0 [50], TMHMM v2.0 [51], and, HMMER v3.3.2 [52].The de novo transcriptome was then functionally annotated using the tool TRINO-TATE v3.2.2 [45].

Differential gene expression analysis
Before the differential gene expression analysis, the SALMON tool was used to align the clean reads of the three sites to the filtered transcriptome and the mapping rate was calculated.The mapping tables were merged according to the stations and using TMM [53] normalization.Differential expression analysis with three biological replicates of the two stations with higher chloride (A and C) were compared to the reference condition (Station B) using the DESeq2 [54] package, the fold change cut-off was set at >2, and the FDR�0.05.This analysis assesses thousands of gene expression differences between the two stations and the control.Therefore, the incorporation of the false discovery rate (FDR) within DESeq2 is essential for mitigating type I errors, making it a widely accepted standard in genomics for correcting false positive results in differential expression analysis [55,56].A heatmap plot using GGPLOT2 [57] was used to represent the differentially expressed genes (DEGs) between each treatment and the control site.The significant upregulated and downregulated DEGs ontology terms were annotated and a gene ontology (GO) pathway enrichment analysis was performed using ShinyGo [58] (version 0.61).An enrichment analysis finds the GO terms with unbalanced distribution between two groups of genes or probe sets.To account for multiple testing issues, a FDR correction was performed on the P-values of this comparison to control falsely rejected hypotheses [59].The genes identified in the GO analysis to be related to the pathways involved in osmoregulation were extracted and listed.The assembly was performed in the High-performance computing system of the university of Cologne (CHEOPS) and the rest of the analysis was performed using the Galaxy project platform [60].

RNA-seq qPCR validation
To validate the results from the RNA-seq analysis, a subset of eight osmoregulation-related upregulated genes was chosen and assessed via real-time quantitative PCR analysis (qPCR).The primer design was performed using the Primer-Blast software [61], and all primer sequences are listed in Table 1.Firstly, reverse transcription of the total RNA samples was performed using the ReverAid first-strand cDNA synthesis kit (Thermo Scientific), and the resulting cDNA was diluted 10-fold and used as a template for qPCR.The qPCR was performed using PowerUp™ SYBR™ Green Master Mix in a StepOne Real-Time PCR system following the manufacturer's indications.The qPCR was carried out in triplicate for each sample with a final well volume of 20μl.An initial 10-minute hot start at 95˚C, followed by 45 cycles of 30 seconds of denaturation at 95˚C and 1-minute annealing at 59˚C.Followed by a melting curve ranging from 60˚C to 95˚C with data acquisition every 0.3˚C.Each gene was assessed using the 2 -44Ct method, with GADPH as the housekeeping gene for normalization, the optimized sequence for the gene was obtained from [37].Finally, the sample with the lower expression of the control station was used as the reference to compare the expression.

Sequence quality control, de novo assembly, and annotation
The RNAseq was performed on the library products from the gill of G. aculeatus, resulting in an average of 47,945,147, 64,032,899.33, and 68,881,883.67 total reads from stations A, B and C respectively.The quality score %Q30, which states the sum of the reads with a base call accuracy of 99, 9% was 99.32%, 99.36% and 99.37% respectively.After quality control, an average of 44,571,717, 59,441,980.67 and 64,633,936 were retained for stations A, B and C respectively.All the retained reads had >Q30.The assembly with all sequences produced a raw transcriptome of 559.3mb, with a N50 of 3797bp, with 3319 complete BUSCO's and a completeness of 91.2%.The transcriptome was filtered for low-expression transcripts, and from a total 340418 in the raw transcriptome, 98643 were retained (28.98%).The filtered transcriptome outcome was 182.4mb of sequences, with a N50 of 3735bp, 2334 complete BUSCO's and a completeness of 64.1%.In total 50155 transcripts were annotated which represents 50.8% of the de novo transcriptome.

Differential gene expression analysis
The mapping rate of the individual sequences to the filtered transcriptome was an average of 60.68%, 59.78% and 60.18% on stations A, B and C respectively.When comparing the station with mid-salinity (A) with the reference site (B), a total of 627 genes (306 down-and 321-regulated) (Fig 2).For the comparison of the station with high-salinity (C) with the reference site (B), a total of 611 genes (270 down-and 341regulated) (Fig 3).In both heatmaps the pattern of expression was similar in all three sample replicates per station.The significant DEGs were assigned to three major Gene Ontology (GO) categories: biological process, cellular component, and molecular function.The DEGs were classified in a total of 200 and 193 higher-level GO terms for A and C, respectively, relative to the reference station (S1 File).Our comparative analysis of the high-level GO terms within the significant DEG dataset, identified potential osmoregulation-related terms, including transporter activity (GO:0005215), transmembrane transporter activity (GO:0022857), and channel regulator activity (GO:0016247).The GO ontology analysis was used to find which terms were enriched in the two higher salinity stations compared with the reference site for the upregulated and downregulated genes separately.Furthermore, to identify the candidate DEGs related to these particular osmoregulation-related GO terms the gene codes were extracted and assigned to the correspondent functional category term (S2 File).For the upregulated genes we found a set of enriched GO terms for the three categories, molecular function, cellular component and biological process (   found enriched GO terms related to osmoregulation only in the molecular function GO category: Anion binding (GO:0043168) and ATP binding (GO:0005524).The main enriched GO terms related to osmoregulation in the downregulated DEGs are in the molecular function GO category for station C vs B and are mostly related to the Ras, Rho and Rab small GTPases terms.These related terms are: Ras guanyl-nucleotide exchange factor activity (GO:0005085), Regulation of Ras protein signal transduction (GO:0046578), Rho guanyl-nucleotide exchange factor activity (GO:0005085), Rho GTPase binding (GO:0017048), GTPase regulator activity (GO:0030695), GTPase activator activity (GO:0005096) and Rho GTPase binding (GO:0031267).
In summary we found enriched GO terms specifically associated with osmoregulation in both upregulated and downregulated genes within the higher salinity stations compared to the reference site.
From the enriched GO terms that were identified to be related to osmoregulation for both higher salinity stations, we extracted all the DEGs and listed them in S3 File.Then we performed a literature search of all the extracted genes and listed the DEGs with referenced osmoregulation involvement (Table 2).For station C vs B, the genes with the higher upregulation were the protein transporters SEC31A, followed by SEC24B.Then, the regulators of the ion channel and ion transporters PTPN3, PTPN6, AHCYL1 and NEDD4 were also found to be upregulated.In the last position were other genes involved in ion transport, the cation channel protein anoctamin 9 ANO9 and the solute carrier SLC9B2.The strongest downregulated genes in C vs B, were the Ras, Rho and Rab small GTPases related genes ITPKB and AKAP13.And finally, the genes involved in ion transport, the ion channel TRPM7 and the solute carrier SLC9B2.For station A vs B, the genes with the higher upregulation were the regulators of ion channel and ion transporters PTPN3, NEDD4 and NEDD4L.Followed closely by the channel protein-encoding gene AQP1 and the solute carrier SLC35B2.In the last position with lower expression were the channel regulators PRF1, CLCN7, the solute carrier SLC6A15, and the channel protein Anoctamin 10 ANO10.The strongest downregulation was recorded for the membrane tension regulator SCARB2, the transporter ABCE1, the channel regulators VPS35 and VP26C and the protein complex gene PIKFYVE.In the last position was the anion exchanger ATAD2.

RNA-seq data validation by qPCR
The analysis of the melting curve of the qPCR resulted in a single product for all genes.The visual comparison of the relative fold changes of the subset of eight DEGs showed a consistently positive direction in both qPCR and RNA-Seq analysis (Fig 6), which suggests that the RNA-seq analysis is reliable.

Discussion
Our de novo transcriptome assembly yielded over 50 million reads with 100% of reads surpassing Q30, and a N50 length of 3735bp, indicating high-quality raw data.The mapping rate of the individual sequences to the transcriptome, which was consistent with previous studies [29,62], was close to 60% for all stations, demonstrating the effectiveness and reliability of the high-throughput sequencing transcriptomic analysis.Additionally, the maintained trend in gene expression between qRT-PCR and RNA-seq validated our results and confirms their suitability for further analysis.The results of the quality control mechanisms validate our sampling design, demonstrating that a depth coverage of 50 million reads and three biological replicates per treatment suffice for the differential gene expression analysis using a de novo transcriptome approach.According to our prediction, we found a differential expression of the genes in the gill tissue of G. aculeatus under a gradient of salinity in the field.Notably, even when chloride concentrations were well below reported toxic thresholds, our comparative transcriptomics approach was able to evidence the activation of the energetically costly osmoregulatory system.Furthermore, we were able to identify candidate genes related to increased osmoregulatory activity even for subtle changes of chloride concentration.In the ontology analysis, we found an enrichment of DEGs genes predominantly related to transmembrane transport and regulation of transport.It is to highlight, that in neither of the treatments we found DEGs with functional categories that could be related to the handling or euthanasia of the fish.The pathways and genes that we found fall in one or more of the big three classes of proteins in charge of moving ions across membranes, ion channels, ion pumps and ion transporters which are essential for all living organisms [63].These key groups of proteins maintain the osmotic pressure in exchange for considerable quantities of energy, and together receive the name "Transportome" [64].The transportome is in charge of maintaining biological homeostasis through a balance between the intracellular and extracellular environments and currently accounts for 10% of the protein-coding human genome [65].For instance, studies with Chinese mitten crab, Guppy Poecilia reticulata, Acipenser baeri, hybrid tilapia, and G. aculeatus showed a significant number of DEGs involved in ion transport and transmembrane transport in varying saline habitats from estuaries to marine environment [26,28,29,32,66].These results confirm the importance of the transportome in the adaptive response to salinity changes and align with the enrichment of genes in related categories in our study.
The enrichment of proteins related to gill tissue transportome in response to salinity changes in our study aligns with consistent findings across various transcriptomic investigations, spanning from invertebrates to fish, as elaborated below.Malik and Kim (2021) identified that about 40% (162/405) of the DEGs were classified as belonging to the transportome in a meta-analysis of four different salinity challenge studies in Chinese mitten crab.In a study with the Guppy Poecilia reticulata it was also shown that a large number of DEGs were involved in ion transport and transmembrane transport [26].Similarly, when hybrid tilapia underwent various laboratory salinity stress treatments,, numerous ion transport enzymes and ion transporters associated with osmoregulation were identified [32].Another example, involves the Siberian sturgeon (Acipenser baeri), where gill transcriptome was sequenced after a high salinity stress exposure.Several pathways were enriched under salinity stress, and as in our case there was a set of DEGs related to osmoregulation like genes of the CLC family and Small GTPases, together with other ion transport mechanisms when challenged by salinity [29].In additional comparative gill transcriptomic studies, the euryhaline fishes Acanthogobius ommaturus, Lates calcarifer and Lateolabrax maculatus were exposed to a salinity gradient manipulation.Strikingly, their GO analysis of the gill transcriptome responded similar to ours, with an enrichment of ion transport, cell volume regulation and signaling pathways [33,62,67].While we found only one example of comparative gill transcriptomics involving G. aculeatus, the study revealed differential expression of several genes related to 'transmembrane transport' when exposed to varying salinities [28].Even when all the previously cited studies converge with our findings in the enrichment of osmoregulatory, cellular signaling and transmembrane ion transport in the gill transcriptomic of fish when facing salinity stress.It is noteworthy that that our findings add to the previous knowledge by demonstrating this osmoregulatory pathway enrichment in wild-caught fish under a natural gradient of salinity stress.
In the following sections, we discuss the relevance of the enriched pathways for osmoregulation, the specific function of the candidate genes, and the ecological implications of our findings for euryhaline fish such as G. aculeatus in the face of global freshwater salinization.

Membrane transporters DEGs
Membrane transporters were the most abundant group of candidate genes that were differentially expressed in the stations with higher chloride concentrations.The candidate genes related to these pathways are VPS35, VPS26, SLC35B2, SLC9B2, SLCO5A1, SLC6A1, SEC31A, and SEC24.The Vacuolar protein sorting-associated proteins VPS35 and VPS26 are components of the cargo recognition sub-complex of the retromer, and together they recognize and bind to specific sorting signals on cargo proteins to be transported back to the Golgi apparatus [68,69].Both proteins are crucial in cell homeostasis, by regulating the transport of particles and by controlling the water reabsorption protein Aquaporin-3 (Aquaporin-3) [70].According to the review by Verri et al., (2012) [71], the SLCs transporters in teleost fish, are conformed by over 50 families with 380 known members and facilitate the transportation of nearly all soluble molecules across cellular membranes.SLC35B2 is part of the nucleosidesugar transporter family (SLC35) that is involved in the endoplasmic reticulum and Golgi of eukaryotic cells [72].Similar to our findings, differential expression of SLC35 family members was reported in the gills of the Chinese mitten crab in response to changes in salinity [66].A remarkable finding was the upregulation of SLC9B2, as the members of the SLC9 family play a critical role in regulating pH in cells and organelles, as well as maintaining acid-base and volume homeostasis [73].In previous studies a high upregulation of SLC9 genes has been recorded under salinity stress, in particular, SLC9B2 was found to function as a Na + /H + exchanger in an experiment using Xenopus oocytes [74].SLCO5A1 from the SLC5 Sodiumglucose cotransporter family and SLC6A15 from the SLC6 Sodium-and chloride-dependent neurotransmitter transporter family were also differentially expressed in our study and are also important transporters found in fish [71].Furthermore, The SLCs have shown to hold high importance in generating counter-ion fluxes to maintain ion homeostasis while avoiding overtly high membrane tension [75].In addition to observing an enrichment of SLC proteins under salinity stress, we also found an increase in proteins related to the vesicular traffic of the SEC family, such as SEC31A as it was also found in Malik and Kim (2021), and SEC24B.The impairment of the expression by glucocorticoids of members of the SEC24 family has been related to the aberrant ion and macromolecular transport in Zebrafish [76].Overall, the findings suggest that these candidate genes play an important role in maintaining ion and macromolecular transport under salinity stress.

Ion channel regulation and transport DEGs
Another important group of DEGs enriched terms was the one related to ion channel regulation and transport.Examples of this that we pinpointed as candidate genes are PTPN3, PTPN6, NEDD4, NEDD4L, AHCYL1, AQP1, TRPM7, PRF1, CLCN7, ANO9 and ANO10.This holds particular importance as ion channels such as voltage-gated, ligand-gated, and second messenger-gated channels, nonselective cation channels, and epithelial Na + and Cl -channels are key to maintaining cellular homeostasis [77].And as described in Davis et al., (2001) [78] protein tyrosine phosphatases as PTPN3 and PTPN6 are involved in tyrosine phosphorylation, which is the primary means of regulating the majority of voltage-gated, ligand-gated, and second messenger-gated channels.PTPN3 has been found to regulate the high-osmolarity glycerol (HOG) MAPK pathway, which is essential for yeast survival in a hyperosmotic environment [79].Other important regulators that we propose as candidate genes are NEDD4 and NEDD4L, which have been found to mediate ubiquitination and degradation of AQP2 which is a key protein in water homeostasis [80].Also, the AHCYL1 protein has been identified as a regulator of the intracellular Ca 2+ channel inositol 1,4,5-trisphosphate (IP3), and multiple ion channel and ion transporters [63].AQP1 is a member of the aquaporin gene family, a group of channel proteins that are involved in the transport and other solutes in the presence of an osmotic gradient.AQP1 differential expression has been previously detected in fish gill, kidney and gut tissue in response to salinity stress for various species [81,82].TRPM7 is a plasma-membrane protein expressed ubiquitously, which belongs to the melastatin-related transient receptor-potential ion channel TRPM subfamily and possesses both ion channel and α-kinase domains.It has been observed that TRPM7 is involved in channel regulation via changes in intracellular Ca 2+ concentration when subject to osmolality gradients [83].The perforin (PRF1) was also differentially expressed, which is relevant as it is known to be able to polymerize and form channels in cell membranes of target cells [84].Also relevant was the differential expression of the chloride channel (CLC7), which was also found to be expressed under osmotic stress in the gills of tilapia under osmotic stress [32].Finally, we found upregulation of both transmembrane proteins anoctamin 9 and 10 (ANO9 and ANO10).Anoctamins are a family of Ca 2+ -activated Cl − channels and phospholipid scramblases that have been shown to support cell volume regulation [85,86].The differential expression of these genes highlights their potential role in cellular homeostasis regulation under osmoregulatory stress conditions.Further studies on these candidate genes could lead to a better understanding of the adaptive mechanisms used by fish to cope with environmental stressors.

Small GTPases, anion and ATP binding DEGs
Another set of pathways that were found to be differentially expressed are those related to the Ras, Rho and Rab small GTPases, anion and ATP binding genes.The Ras, Rho and Rab small GTPases proteins are key regulator and signaling mediators of a wide variety of process that occur in eukaryotic cells [87].The main importance of these proteins in the context of our study lies in their direct interaction with ion channels to regulate ion homeostasis [88].Furthermore, Ras GTPases are involved in the process of membrane thinning and curvature with effects on osmoregulation [89].One important candidate gene found in our study was the Rho GTPase AKAP13, which is a key regulator of ion homeostasis under osmotic stress.In detail, AKAP13 attracts JIP4 and activates NFAT5 through the Rho-type small G-proteins and p38 MAPK signalling pathway, regulating intracellular osmolarity [90].ITPKB regulates the subcellular distribution of Rasa3, a Ras GTPase-activating protein and has been observed at is differentially expressed when hyper-osmotically stressed [91,92].In the gene ontology categories, anion and ATP binding we highlight the candidate genes ABCE1, SCARB2, PIKFYVE and ATAD2.ABCE1 is a member of the ABC gene family, which constitutes the largest group of transmembrane transporter proteins encoded in the human genome.These proteins utilize ATP as an energy source to facilitate the transportation of various molecules across cellular membranes keeping cell homeostasis [93].The lysosomal integral membrane protein 2 (SCARB2) is a transporter of cholesterol that is suspected to control membrane tension and ion trafficking via the direct insertion of lipids [94].PIKFYVE is the sole source of all the Phosphoinositide lipids PI(5)P pool in yeast membranes and is a key regulator of ion transport, and it is upregulated under hyperosmotic shock in yeast [95].Finally, another interesting candidate is ATAD2, a gene from the ATPase family.Its relevance resides in the fact that ATPases are enzymes that utilize ATP to power many cellular processes including transmembrane transport and osmoregulation in fish [96,97].Note that the differential expression of ABCE1 and ATAD2 provides direct evidence of the energetic cost of osmoregulation, even at low salinity concentrations in our study.The proposed candidate genes play important roles in transmembrane transport, ion trafficking, and osmoregulation.

Ecological implications
Due to the inherent complexities of our field study design, where several uncontrolled factors may influence the gene expression of G. aculeatus, we are cautious in establishing a direct causal relationship between the differential gene expression and the salinity gradient of the stations.Nonetheless, our results strongly suggest that the gill osmoregulatory mechanism of G. aculeatus is activated at sublethal concentrations of chloride, and we have identified candidate genes that could serve as indicators of salinity stress.Using a transcriptomic approach, we detected activation of osmoregulatory pathways in both study sites, even below the thresholds established by German authorities for running waters.This finding is significant, as transcriptomics has been successful in determining the sublethal effects of various stressors and giving new insights into the management of species and ecosystems.For example, Komoroske et al., (2016) [39], were able to determine that the endangered delta smelt (Hypomesus transpacificus) osmoregulatory system was activated at a high energetically cost under sublethal salinity thresholds.And for the same species, again using transcriptomic coupled with physiological methods, it was demonstrated that sublethal critical effects are observed 4-6˚C below the previously established acute tolerance limits [98].The relevance of this transcriptomic study was further demonstrated when Brown et al., (2016) [99], combined the information on these sublethal thresholds to model the future habitat suitability of this species.Based on these precedents, we predict that our results have the potential to inform models on the expanding distribution beyond its native range of G. aculeatus in response to freshwater salinization.G. acualeatus populations have been increasing rapidly in several ecosystems around the world in response to altered conditions, with important negative impacts on native species that have led to consider it an invasive species in places such as lake Constance (Germany) [100,101].An effect that has been previously described for other invasive species that have demonstrated adaptive advantage in increasing salinity scenarios, increasing their populations and their distributions [102,103].
It is worth noting that the candidate genes related to osmoregulation discussed in our study, which showed differential expression, have primarily been identified in manipulative experiments involving significant changes in salinity.This highlights the relevance of our findings, as we have identified the same candidate genes within the subtle salinity gradient of our field experiment, suggesting their potential as biomarkers for rapidly assessing sublethal salinity stress in field applications.Furthermore, we were able to measure a physiological response to a stressor in a wild non-model organisms in field conditions where many other stressors may be at interplay.Therefore, our study outlines the potential of transcriptomic approaches for assessing physiological responses of wild organisms in conditions with known or unknown stressors and be tools for informing management and policy decisions [22,25,104].As our understanding of the mechanisms of osmoregulation advances with the use of transcriptomics, more biomarkers or candidate genes are identified.The expression of these selected biomarkers can then be used in transcriptomic arrays to evaluate the responses of individual fish to several important biological and environmental stressors such as water-quality stressors, aquatic invasive species, and climate change [105].Thus, further studies into the mechanisms and candidate genes that we present here for sublethal salinity stress have the potential to have an impact on the conservation and management of freshwater ecosystems in the prospect of increasing freshwater salinization.In the light of our findings, we suggest to perform further studies into other non-model wild-caught fish physiological stress response using a combination of techniques that include transcriptomics.We consider that by including these, the current thresholds for salinity and other pollutants in freshwater could be revisited as sublethal effects would be accounted for.Our study has shown that even small changes in salinity significantly change gene expression patterns in fish.Enhanced salinity could pose energetic costs to fish thereby reducing long term survival and reproductive output.This emphasizes the critical need to preserve the highest water quality in streams, recognizing that sublethal stressors can exert profound impacts on ecosystems.

Conclusion
In conclusion, our transcriptomics analysis of the gill tissue of G. aculeatus under a gradient of salinity in the field strongly suggests a significant differential expression of genes related to membrane transport and regulation of transport.Our results support the notion that the transportome plays a crucial role in the adaptive response to salinity changes in euryhaline fish and highlight the activation of energetically costly osmoregulatory systems even when chloride concentrations are below the established toxic thresholds.We would like to highlight our finding of a set of candidate genes related to the sublethal concentration of chloride that is involved in transmembrane transport and ion homeostasis and can be key to the further development of transcriptomic tools for evaluating stress response.Overall, our study underscores the importance of understanding the molecular mechanisms underlying osmoregulation in the face of global freshwater salinization and provides a basis for further investigation into the impact of environmental stressors on fish populations.

Fig 1 .
Fig 1. Location of the sampling station in the Boye catchment.The Chloride concentrations in the period of sampling for the upstream station B (Control) are 12.912 mg/L.And the two downstream treatment stations: Station A 26.378 mg/L, and station C 32 mg/L.https://doi.org/10.1371/journal.pone.0298213.g001 Fig 4).The main enriched GO terms related to osmoregulation in the upregulated DEGs in the biological process GO category for station A vs B are: Regulation of sodium ion transmembrane transport (GO:1902305 and GO:1902306), regulation of sodium ion transmembrane transporter activity (GO:2000649 and GO:2000650), and regulation of sodium ion transport (GO:0002028 and GO:0010766).For the molecular function the GO categories main enriched terms were related with the activation and regulation of sodium channels (GO:0017080,

Fig 2 .
Fig 2. Heatmap displaying 627 differentially expressed genes between the mid-salinity station (A) and the reference station (B).The X-axis represents sample replicates per station, and the Y-axis represents individual gene expression.Upregulated genes are shown in yellow, with brighter colors indicating higher expression values.In contrast, purple shades indicate downregulated genes, with the brightest shade indicating the strongest downregulation.https://doi.org/10.1371/journal.pone.0298213.g002

Fig 3 .
Fig 3. Heatmap displaying 611 differentially expressed genes between the high-salinity station (A) and the reference station (B).The X-axis represents sample replicates per station, and the Y-axis represents individual gene expression.Upregulated genes are shown in yellow, with brighter colors indicating higher expression values.In contrast, purple shades indicate downregulated genes, with the brightest shade indicating the strongest downregulation.https://doi.org/10.1371/journal.pone.0298213.g003

Fig 4 .
Fig 4. Gene ontology analysis of enriched upregulated genes for stations A and C compared to the reference site B. The number of genes is displayed on the X axis, while the Y axis represents GO categories colour-coded based on molecular function (red), cellular component (green), and biological process (blue).Could you unify the X axis to make a comparison easier?(max value AB = 16, max value CB = 18).https://doi.org/10.1371/journal.pone.0298213.g004

Fig 5 .
Fig 5. Gene ontology analysis of enriched downregulated genes for stations A and C compared to the reference site B. The number of genes is displayed on the X axis, while the Y axis represents GO categories colour-coded based on molecular function (red), cellular component (green), and biological process (blue).X axes-see above.https://doi.org/10.1371/journal.pone.0298213.g005

Fig 6 .
Fig 6. qPCR validation plot of RNA-seq data.It shows the expression of eight genes, including NEDD4 and PTPN3, in two treatment groups.NEDD4 and PTPN3 refers to group C vs. B, while NEDD4-A and PTPN3-A refer to group A vs. C.The gene names are on the X axis, and the Y axis displays the relative fold change in gene expression (Log2FCs).Mean ± standard deviation Log2FC values are presented.https://doi.org/10.1371/journal.pone.0298213.g006

Table 2 . List of the selected candidate genes involved in osmoregulation generated from the Gene Ontology analysis.
Upregulated genes are highlighted in red, while downregulated genes are highlighted in blue, with intensity indicating expression level (Log2 FC).The false discovery rate (FDR) value is also included, indicating that an FDR p<0.05 suggests a non-random differential expression of the gene. https://doi.org/10.1371/journal.pone.0298213.t002