Changes in transcriptomic response to salinity stress induce the brackish water adaptation in a freshwater snail

Studying the mechanisms of the establishment of a population in a novel environment allows us to examine the process of local adaptations and subsequent range expansion. In a river system, detecting genetic or phenotypic differences between a freshwater and brackish water population could contribute to our understanding of the initial process of brackish water adaptation. Here, we investigated behavioural and gene expression responses to salt water in a freshwater and brackish water population of the freshwater snail, Semisulcospira reiniana. Although the individuals in brackish water exhibited significantly higher activity in saltwater than freshwater individuals just after sampling, the activity of freshwater individuals had increased in the second observation after rearing, suggesting that their salinity tolerance was plastic rather than genetic. We found 476 and 1002 differentially expressed genes across salinity conditions in the freshwater and brackish water populations, respectively. The major biological process involved in the salinity response of the freshwater population was the biosynthesis and metabolic processing of nitrogen-containing compounds, but that of the brackish water population was influenced by the chitin metabolic process. These results suggest that phenotypic plasticity induces adaptation to brackish water in the freshwater snail by modifying its physiological response to salinity.

Studies concerning distribution patterns of species or populations along the environmental continuum play a significant role in discovering evidence of environmental adaptations. When a species encounters a novel environment during range expansion, adaptation to the environment could aid in the establishment of a population. Adaptations are achieved via two different processes: adaptive evolution and phenotypic plasticity 1,2 . For instance, invasive species, which are introduced into non-native regions by human activity, have been studied extensively with respect to the ongoing process of their adaptive evolution, as they have acquired various traits that are suitable for a novel environment through rapid evolution [3][4][5] . Evolutionary adaptations, of course, contribute to colonization and expansion into a novel environment in native species [6][7][8] . Phenotypic plasticity could enhance adaptation to a novel environment, as well as range expansion, because it should enable a species to cope with unfamiliar conditions [9][10][11] . Both theoretical and empirical studies suggest that individuals with high levels of phenotypic plasticity dominate the edges of distribution ranges 12,13 . For example, plasticity in thermal tolerance is an important trait facilitating survival at the range margins for the fruit fly Ceratitis capitata 14 . In general, plasticity makes a significant contribution in adaptation to heterogeneous environments at small spatial scales, or to rapid environmental changes 15,16 . Adaptive plasticity is hypothesized to provide a rapid response to the environment, and also to precede and facilitate adaptive evolution [17][18][19][20][21] . Therefore, understanding the process of both adaptive evolution and phenotypic plasticity is important for revealing the mechanism of range expansion.
In systems in which different environmental conditions are adjacent to each other at relatively small spatial scales, individuals have opportunities to colonize new habitats outside of their current distribution, and thus range expansion can be achieved if a species or population adapts to the environment. Therefore, we can study the ongoing process of adaptation by examining the evolutionary or plastic change in traits along an environmental gradient. A river is one such system, because a steep gradient of salt concentration, which is a significant environmental factor influencing the distribution range of many aquatic species [22][23][24] , is observed near the estuary

Results
Genetic differentiation. We obtained 983,061 transcripts with a mean length of 556 bp via de novo assembly. Among the 978 metazoan core gene orthologues, 969 genes (99.1%) were identified completely. We ran the CD-HIT program using 164,656 putative coding sequences (CDS), which provided 94,676 clustered CDS transcripts. Using these CDS transcripts, we identified 5871 SNPs on the contigs annotated against all protein sequences of Crassostrea gigas by BLAST search. We used these SNPs for the following population genetic analysis.
The BayeScan program did not identify any outlier loci, indicating that all SNPs were selectively neutral between the two populations. Using all 5871 SNPs, the F ST value between the freshwater and brackish water populations was estimated at 0.011.
Behavioural response to saltwater. The locomotive activity in 0% saltwater was not significantly different between populations, either in the first observation immediately after collecting, or in the second observation after rearing in experimental freshwater for one week following the first observation (population [P]: χ 2 = 0.02, P = 0.89; the number of days after collecting [D]: χ 2 = 0.03, P = 0.85; P × D: χ 2 = 0.03, P = 0.87, Fig. 1a). In contrast, we found interaction effects for the locomotive activity in 0.4% saltwater (P: χ 2 = 1.06, P = 0.30; D: χ 2 = 0.64, P = 0.42; P × D: χ 2 = 11.3, P < 0.001, Fig. 1b), indicating a differential response to the experience of being reared in an experimental freshwater setup for one week between populations. Post hoc comparison found significant differences in locomotive activity between populations on the first observation, immediately after sampling. The locomotive activity of the freshwater population in 0.4% saltwater was higher during the second observation, after rearing in freshwater. Although we were not able to identify a statistically significant differ- www.nature.com/scientificreports/ ence, the locomotive activity of the brackish water population in 0.4% saltwater tended to decrease between the first observation and the second observation. Consequently, no significant difference between populations in locomotive activity in 0.4% saltwater was found at the second behavioural observation.
Gene expression response to salinity. All freshwater and brackish water individuals were active in 0% saltwater. In 0.6% saltwater, all brackish water individuals were active, while three of the four freshwater individuals were not. In 0.8% saltwater, two of the four brackish water individuals were active, but all the freshwater individuals withdrew into their shells to avoid the salinity. These results supported the hypothesis that brackish water individuals have a higher tolerance to salinity, as shown previously. Among the164,656 putative CDS, 81,495 (49.5%) of them were annotated against all protein sequences of C. gigas, using BLAST search. These annotated sequences comprised 34,816 genes. We used these genes for gene expression analysis.
The number of differentially expressed genes (DEGs) among the three salt concentrations was higher in the freshwater population (1002) than in the brackish population (476), indicating that many more genes were differentially expressed in the freshwater population than in the brackish water population ( Fig. 2a; Supplementary Tables S1 and S2). Among these DEGs, 48 genes were shared between the two populations, but the expression patterns of the shared DEGs differed from each other ( Fig. 2b; Supplementary Table S3). The number of genes highly expressed in 0.8% saltwater was higher than that expressed in other saltwater concentrations, irrespective of the population, probably due to strong salinity stress (Fig. 2c).
The GO terms that were significantly enriched within the set of DEGs in the freshwater and brackish water populations are summarized in Supplementary Tables S4 and S5, respectively. We found 62 terms in the freshwater population and 53 terms in the brackish one. In the freshwater population, 45 terms were unique to the www.nature.com/scientificreports/ population (Table 1), including those involved in the biosynthesis and metabolic processes of nitrogen-containing compounds. In the brackish water population, 36 terms were unique to the population (Table 2), including those involved in the chitin metabolic process and urea transport.

Discussion
Either adaptive evolution or phenotypic plasticity, or both, can contribute to the establishment of a population in a novel environment 1,2 . In a river, salinity gradually increases toward the estuary, creating a steep environmental gradient. Freshwater organisms typically do not possess the physiological mechanisms to cope with high salinity. www.nature.com/scientificreports/ The distribution of freshwater species is limited to a specific range between brackish and freshwater areas. Therefore, range expansion to brackish areas requires tolerance to salinity, and thus brackish water populations display different salinity responses than freshwater populations. In the present study, we confirmed that individuals from a population living in brackish water displayed different salinity responses from those of a freshwater population, in terms of both activity and physiology. Our findings could clarify the understanding of the mechanisms of brackish water adaptations via adaptive evolution and/or phenotypic plasticity. The F ST value between two populations, which we observed, is very low (0.011) compared with previous studies in snails 28,29 . This value of F ST is considered to indicate that the two adjacent populations exhibit little genetic differentiation 30 , supporting our assumption. The brackish water population is suggested to have colonized their novel environment recently or regularly receive immigrants from the freshwater population. The brackish water population appears to be in the initial stages of adaptation to brackish water, and thus few mutations have accumulated in the population. Therefore, in our study system, we could observe a phenomenon that occurs in the early stages of brackish water adaptation.
In 0.4% saltwater, we found significantly higher activity levels in individuals in brackish water than in individuals in fresh water immediately after transfer from one environment to the other, indicating that individuals in brackish water can be active even in conditions of salinity. The lack of significant differences between the activity levels of populations at the second observation is probably due to acclimation to the same rearing environment for one week. The increased activity of freshwater individuals at the second observation may be due to their acclimation to unexpected changes in water conditions, because the 0% saltwater used for keeping individuals www.nature.com/scientificreports/ may be slightly hypertonic due to the effect of their faeces and a crushed oyster shell. The lack of difference between the activity levels of populations in 0.4% saltwater after one week suggests that individuals have the ability to change their activity level according to water conditions. If higher activity in saltwater is derived from high tolerance to salinity, our results suggest that phenotypic plasticity may contribute to tolerance to salinity in individuals in brackish water. The hypothesis that phenotypic plasticity may contribute to salinity tolerance is supported by the absence of outlier loci in the two populations, which showed distinct responses to salinity. However, we cannot rule out adaptive evolution caused by genetic differences, especially in non-coding regions. The relationship between phenotypic plasticity and adaptation to novel environments has been long studied. Baldwin 17 suggested that adaptive plasticity plays an important role in the establishment and persistence of populations in new environments. Also, an adaptive phenotype initially accomplished by plasticity is hypothesized to sometimes become genetically-encoded, so that establishment in a new environment is achieved 18,19 . Our results indicate that S. reiniana can acclimate to salinity even in individuals grown in fresh water. Individuals in brackish water may cope with salinity via phenotypic plasticity, and may be in the process of adaptation to a brackish environment as suggested by Baldwin's hypothesis.
In our transcriptome analysis, we found that the expression level of 48 genes changed with salinity in the two populations, suggesting that these genes could be involved in the basis of saltwater response in S. reiniana. We detected four genes encoding proteins with von Willebrand factor domains within the shared DEG set. Von Willebrand factor, which plays a key role in normal hemostasis, is suggested to be expressed when extracellular NaCl levels are elevated 31 . High salinity would have a deleterious effect on freshwater snails, by promoting excess thrombus formation. In the freshwater population we found several genes involved in the biosynthesis and metabolic processing of nitrogen-containing compounds, such as glutamate synthase and glutamine synthetase. Gene ontology analysis also revealed DEGs enriched in nitrogen compound biosynthetic and metabolic processes. Nitrogen-containing compounds accumulate in plant species subjected to salinity stress and are involved in osmoregulation 32 . Nitrogen metabolism is also involved in the salinity responses of many molluscan species 33 , indicating that individuals in freshwater have common processes in response to salinity. A different pattern of gene expression was observed in the brackish water population, in which we noted that several genes associated with chitin metabolism were differentially expressed, including chitin synthase C and chitotriosidase-1. While chitin is an essential component of mollusc shells, its function in osmoregulation has not been reported. However, Lv et al. 34 suggested that chitinase is associated with the response to salinity in crustacea, leading to the hypothesis that a similar process occurs in S. reiniana. DEGs in the brackish water population were enriched not only in chitin metabolism but also in urea transport. Urea functions as an osmolyte, and holds special importance for cell volume preservation in an aquatic snail under hyperosmotic stress 35 . However, it is not clear why genes associated with urea were enriched only in the brackish water population and those associated with nitrogen-containing compounds only in the freshwater population; although both of them are osmolytes which act to prevent water loss in hyperosmotic environments. In summary, more than half of the GO terms that were enriched within each DEG set did not overlap, indicating that the biological processes in response to salinity were different in the two populations. These differences in gene expression patterns may be caused by plasticity in brackish water adaptations. Freshwater individuals would acclimate to saltwater in a few days, and individuals in both populations would exhibit similar gene expression patterns.
Our results indicate that individuals in fresh water and brackish water show different responses to salinity with respect to activity levels and gene expression patterns. The establishment of a population in a brackish water area requires the activation of biological processes that can cope with high salinity. Given the results of our behavioural and transcriptome analyses, metabolic processing of chitin and urea may lead to higher activity in saltwater in individuals from brackish water. However, our behavioural observations suggest that individuals from fresh water may have enhanced activity under salinity stress. In summary, acclimation to hypertonic conditions and alteration of salinity responses via physiological processes can contribute to the early stage of brackish water adaptation in freshwater snails.
Note that our results do not explain why range expansion to a brackish water environment has not been achieved in all rivers, despite their potential for acclimation to salinity. Evolutionary differences in the strength of plasticity and long-term salinity tolerance may be associated with the success of brackish water adaptations. Further investigation into the effects of long-term exposure to salinity on the survival, growth, and reproduction of individuals in fresh and brackish water is required.

Materials and methods
Species and sampling. Semisulcospira reiniana is a common freshwater snail in Japan and mainly inhabits the freshwater areas of rivers. We collected adult individuals of S. reiniana from freshwater (35° 15′ 15″ N, 136° 41′ 09″ E) and brackish water (35° 07′ 22″ N, 136° 41′ 27″ E) areas of Kiso river in Japan, in May 2018 and June 2019. While water level fluctuations caused by the tidal cycles occur at both sampling sites, seawater does not reach the freshwater site. The salt concentration of the freshwater site was almost 0%, while the salt concentration at the brackish water site fluctuates from 0 to 0.5% according to the tidal cycle. Since these two populations are located relatively close to each other (15 km apart), we assume that the two populations have not diverged genetically. The snail specimens were preserved in a container (13 × 13 × 20 cm) maintained at 23 °C for one day until the onset of experiments.
Behavioural response to salinity. Individuals collected in freshwater and brackish water populations were examined for activity under concentrations of 0% and 0.4% saltwater. Ten individuals were tested in each salinity condition for each of these two populations (total sample size was 40). Snails were put separately in individual bowls (φ13 × 3.5 cm) filled with 0% or 0.4% saltwater, prepared from decalcified tap water. The bowls Scientific RepoRtS | (2020) 10:16049 | https://doi.org/10.1038/s41598-020-73000-8 www.nature.com/scientificreports/ were kept under laboratory conditions at a temperature of approximately 23 °C. Thirty minutes after putting the individuals into the bowls, they were monitored by a camera (400-CAM061, Sanwa Direct) for two hours to quantify locomotive activity. Images were taken every 30 s and used to create time-lapse movies. We tracked the position of each individual using the tracking software UMATracker 36 . Using the snail trajectories, the total locomotion distance for two hours was calculated for each individual. After the first behavioural observation, we put individuals back into the containers filled with fresh water (decalcified tap water). We added a crushed oyster shell to the containers to keep the water clear. One week later, we conducted the second behavioural observation following the same procedure, and using the same sample set.
Gene expression response to salinity. Three freshwater and three brackish water individuals, reared under the freshwater condition for one week after collection, were exposed to 0%, 0.6%, or 0.8% saltwater for three hours in the laboratory. While 0% saltwater was assumed to be a typical optimum environment for S. reiniana in the natural population, 0.6% and 0.8% were expected to be high and extremely high salinity conditions in the tidal cycle, respectively. Soon after exposure, we checked to see whether the snails were active. Individuals were dissected and their epidermis was preserved in 750 μl of RNAlater Stabilization Solution (Invitrogen). The samples were kept at − 80 °C until the extraction of total RNA. Total RNA was extracted using Maxwell 16 LEV Plant RNA Kit with the Maxwell 16 Research Instrument (Promega) according to the manufacturer's instructions. Electrophoresis on 1% agarose gels was performed to check for RNA degradation. RNA concentrations were estimated using a Qubit 2.0 Fluorometer (Invitrogen). RNA purity was estimated using a BioSpec-nano (Shimadzu). The cDNA library was constructed using TruSeq RNA Sample Prep Kits. Paired-end (150 bp) RNA sequencing (RNA-Seq) was performed on the Illumina NovaSeq6000 platform. After the removal of adaptor sequences and low-quality reads using Trimmomatic 37 , we used FastQC (https ://www.bioin forma tics.babra ham.ac.uk/proje cts/fastq c/) for quality control. The remaining high-quality reads were used for de novo assembly using Trinity 38 .
To estimate gene expression levels, all reads of each sample were mapped to the reference transcripts using RSEM 39 . The read count data was used for gene expression analysis. We searched for homologues of every S. reiniana gene using BLAST searches for all protein sequences of C. gigas. Genes with the best hit and with an e-value < 0.0001 were used for gene expression analysis. DEGs among the three salt concentrations were detected using the TCC package 40,41 . We considered genes with q < 0.05 as DEGs. Gene ontology (GO) enrichment analysis of the DEGs was performed using Blast2GO software 42 . F St calculation. Putative coding regions were extracted from the reference transcripts using TransDecoder (https ://githu b.com/Trans Decod er/Trans Decod er/wiki), providing the reference transcripts only contained CDS. They were clustered based on sequence identities of 90% to remove redundancy, using the CD-HIT program 43 . We used STAR for mapping all reads of each sample to the clustered CDS reference. We then used GATK 44 to identify SNPs. We used the homologues of C. gigas, which were detected in the same manner as described in a previous section. Genes with the best hit and with an e-value < 0.0001 were used to estimate the genetic structure of the two populations. To identify outlier loci which were not selectively neutral between the two populations, we ran the BayeScan program 45 with default parameters. Putatively neutral loci were used to estimate F ST between the two populations using Arlequin 46 .
Statistical analyses. Locomotive distances in 0% and 0.4% saltwater were analyzed by the generalized linear model (GLM) with a Gamma distribution. Population (i.e., freshwater, and brackish water population) and the number of days after collecting (i.e., first and second observation) were included as explanatory variables. Post hoc test for all four pairwise comparisons was conducted using the GLM with a Gamma distribution. Since the analysis was performed four times, we applied the Bonferroni correction to account for multiple comparisons (Bonferroni-adjusted α was 0.0125). All statistical analyses were performed in R 3.4.2.

Data availability
All raw transcriptome read data were deposited in the DDBJ Sequenced Read Archive under accession numbers SAMD00218387-SAMD00218404. Behavioural and transcriptome data are available in the Dryad Data Archive at https ://doi.org/10.5061/dryad .jdfn2 z37w.