A key evolutionary step determining osmoregulatory ability for freshwater colonisation in early life stages of fish

The colonisation of novel environments significantly increases biodiversity (Schluter, 2000; Losos, 2010). A variety of marine animals have repeatedly colonised freshwater habitats, leading to the diversification of species and life histories of snails (Abdou et al., 2015), crustaceans (Anger, 2013), fishes (Betancur-R et al., 2015) and mammals (Hamilton et al., 2001). Diadromous species (e.g. fishes) regularly move between marine and freshwater habitats and are often regarded as evolutionary intermediates, potentially undergoing freshwater colonisation (Gross, 1987; Dodson et al., 2009). The derivation of strictly freshwater species from diadromous species has promoted species radiation in freshwater areas (McDowall, 2001; Waters and Wallis, 2001; Burridge and Waters, 2020). This has contributed to the conspicuous biodiversity in freshwater habitats (Corush, 2019), which contains approximately 40% of teleost fishes despite freshwater habitats covering only 1% of the Earth's surface (Lévêque et al., 2007). To understand the origin of freshwater biodiversity, it is necessary to elucidate the underlying environmental, ecological and physiological factors associated with the loss of diadromy (Burridge and Waters, 2020; Waters et al., 2020).

In teleost fishes, drastic changes in osmoregulatory mechanisms are presumably required for their habitat shift to novel environments with different salinities. Because teleosts maintain their body fluid at approximately one-quarter to one-third the osmotic concentration of seawater, they must discharge water and absorb salt ions (i.e. hyper-osmoregulation) to survive in freshwater conditions (Evans et al., 2005; Edwards and Marshall, 2013). To colonise freshwater habitats, marine fish need to acquire the ability to hyper-osmoregulate, often by developing osmoregulatory plasticity (Bamber and Henderson, 1988; McCairns and Bernatchez, 2010). In several marine-originated diadromous species, adult (or juvenile) fish have been reported to have developed osmoregulatory plasticity and actively enter freshwater areas (Shrimpton, 2013; Inokuchi et al., 2017; Oto, 2019). However, even in such diadromous species, larval fish often appear not to have the ability to inhabit freshwater habitats (Matsumiya et al., 1982; Tsukamoto et al., 2002; Dodson et al., 2009; Schmidt et al., 2014). Most of the osmoregulation-related organs, such as gills, kidneys and intestines, are underdeveloped in early larval fish (Varsamos et al., 2005), which possibly makes it more difficult for them to flexibly change the osmoregulatory process and acclimatise to novel salinity environments than for fish in later stages (McDowall, 1997; Lee and Bell, 1999). It is, therefore, likely that the acquisition of larval freshwater tolerance, which is potentially achieved by enhancing osmoregulatory plasticity, is a crucial step in the completion of freshwater colonisation by fish, even though the physiological and genetic mechanisms underlying the acquisition have not yet been clarified.

Among the three types of diadromy (i.e. anadromy, catadromy and amphidromy), amphidromy provides the most obvious example of how the difficulty in the evolution of larval physiology is a bottleneck to freshwater colonisation. Amphidromy can be found in various marine-originated lineages, including snails, shrimps and fishes, and is the most speciose (273 or more species) among the three types of diadromy (i.e. anadromy, 53 species; catadromy, 109 species; Augspurger et al., 2017). Marine-originated amphidromous animals only visit the sea during the short larval period because the larvae generally lack the physiological capabilities necessary to survive or grow under unusual freshwater conditions (Rome et al., 2009; Iida et al., 2010; Abdou et al., 2015). Notably, in some groups, amphidromous species have repeatedly evolved landlocked or resident forms in freshwater areas (hereinafter collectively referred to as freshwater forms) (Augspurger et al., 2017; Delgado and Ruzzante, 2020). Phylogenetic studies have provided evidence that species radiation has occurred after the loss of amphidromy (i.e. the loss of marine larval dispersal) in some fish families, such as Galaxiidae (Waters and Wallis, 2001; Burridge and Waters, 2020), Gobiidae (Yamasaki et al., 2015) and Eleotridae (Stevens and Hicks, 2009). Thus, focusing on the freshwater forms derived from amphidromous species would help us comprehend the role of physiological evolution in species diversification during freshwater colonisation. However, compared with other types of diadromy, amphidromy has received far less attention from evolutionary research, including from studies on physiological evolution (Shrimpton, 2013).

Among teleosts, the goby family Gobiidae includes more than a hundred amphidromous species (Augspurger et al., 2017; Delgado and Ruzzante, 2020). The marine-originated goby genus Gymnogobius, distributed in East Asia, includes at least three amphidromous, five freshwater, and eight marine and estuarine species (Ellingson et al., 2014; Chiba et al., 2020; Ito et al., 2022). This genus is appealing for studying the processes of freshwater colonisation because of its diversity in life histories and habitats (Fig. 1A). In particular, one clade in this genus contains three amphidromous species that differ in their dependence on freshwater environments: several landlocked populations of Gymnogobius urotaenia (Hilgendorf 1879) have been identified in freshwater lakes, while no landlocked populations have been reported in Gymnogobius petschiliensis (Rendahl 1924) and Gymnogobius opperiens Stevenson 2002 (Aizawa et al., 1994; Harada et al., 2002; Stevenson, 2002; Sota et al., 2005). The latter two species possibly need to migrate to the sea just after hatching. Here, we hypothesise that: (1) G. urotaenia, which flexibly establishes freshwater forms, has evolved a higher freshwater tolerance at the larval stage than the other two congeners, and (2) the freshwater tolerance was acquired through the enhancement of osmoregulatory plasticity.

To investigate these two hypotheses, we conducted a comparative transcriptomic analysis to ascertain whether larvae of G. urotaenia, the amphidromous species that has successfully established freshwater populations, can regulate the expression of osmoregulatory genes more flexibly than the other obligate amphidromous congeners, G. petschiliensis and G. opperiens (Fig. 1B). First, we confirmed the migratory history of the three Gymnogobius species by analysing the trace element concentrations in the otoliths, which vary with habitat salinity (Elsdon et al., 2008). We then quantified the larval survival rates in freshwater and seawater conditions in the laboratory to compare the extent of osmoregulatory ability among the three species. Finally, to identify the evolutionary change responsible for larval survival in freshwater conditions, we analysed the whole-body transcriptomes of larval fish and sought differentially expressed genes between freshwater- and seawater-reared individuals in each species. We investigated both the up-regulation of genes responsible for hyper-osmoregulatory capacity and down-regulation of genes responsible for water uptake and ion discharge (hypo-osmoregulation), which are similarly important for freshwater acclimation in teleosts (Velotta et al., 2017; Nakamura et al., 2020). When the larvae were reared in freshwater, we predicted that the expression of some genes important for hyper- and hypo-osmoregulation would be up- and down-regulated, respectively, in G. urotaenia but would show less/no regulation in the other amphidromous congeners. The results should enable us to confirm the importance of larval freshwater tolerance and osmoregulatory plasticity and clarify key evolutionary changes for enhancing larval freshwater tolerance, which possibly enabled the species to establish freshwater forms.

For the otolith microchemical analysis and larval rearing experiment, male parents and egg clutches of three amphidromous Gymnogobius species were sampled from two or three populations per species inhabiting five rivers and one lake located on the Pacific and Sea of Japan coasts in Japan (population ID: G. petschiliensis PI and PA; G. opperiens OO and ON; G. urotaenia UE, UA and UB; 3–9 clutches, i.e. biological replicates, for each population; Fig. 2; Table S1). We tried to collect 10 clutches for each population but, in some cases, the number was much lower than that because of sampling failures or microbial egg spoilage. We failed to capture male parents protecting eggs for the UE population; therefore, other adult individuals captured at the same sampling site were used in the otolith analysis. The UE clutches were molecularly confirmed to be G. urotaenia using Gymno-cytbL and H15915 as forward and reverse primers, respectively (Aoyama et al., 2000; Tabata et al., 2016). All fish were handled according to animal welfare guidelines and policies of the Kyoto University animal welfare committee and Japanese Ministry of Education, Culture, Sports, Science and Technology (https://www.kyoto-u.ac.jp/en/research/research-compliance-ethics/animal-experiments).

Otolith microchemical analysis was conducted to examine the migratory history of male parents from the seven populations of the three Gymnogobius species (using 2–9 individuals from each population; Table S1). We measured the strontium (Sr) and calcium (Ca) concentrations accumulated in the otolith specimens. The Sr/Ca ratio is influenced by environmental salinity, and is often used to reconstruct migratory patterns of various fishes, including amphidromous gobies (Tsunagawa and Arai, 2008; Iida et al., 2017). By measuring the ratio from the otolith core (formed when hatching) to the edge (formed just before sampling) along the transverse section of the otolith, the history of salinity habitat use can be traced. The specimens were read with an electron probe microanalyser (EPMA: JXA-8230, JEOL Ltd, Tokyo, Japan). Detailed specimen preparation and analysis protocols are available in the Supplementary Materials and Methods.

We conducted a rearing experiment under laboratory conditions using newly hatched larvae from 8–17 clutches for each species to compare the freshwater tolerance of the three Gymnogobius species. When more than half of the eggs had hatched from each clutch, exactly 18 individual larval fish were carefully transferred into 2 l experimental tanks filled with aerated freshwater (salinity 0.1; freshwater-reared group) or seawater (salinity 34.5–35.5; seawater-reared group) using pipettes. During the experiment, the larvae were fed twice a day. Surviving and dead individuals were counted daily until all individuals died. A detailed experimental protocol is available in the Supplementary Materials and Methods.

To evaluate larval tolerance to freshwater in the three Gymnogobius species, we conducted a survival time analysis using the Cox proportional hazard model (Therneau and Grambsch, 2000), which accepts censored data caused by the sampling performed for other experiments. To compare survivorship between the freshwater and seawater groups, Cox proportional hazard models were established for each species using the coxph function in the survival R package version 3.2-13 (https://CRAN.R-project.org/package=survival). In these models, salinity and egg clutch ID were included as the variables explaining the survival (days) of individual larvae. Here, the positive/negative and significance of the coefficient of the salinity group were assessed.

In addition, to directly illustrate interspecific variation in freshwater tolerance and test whether freshwater-tolerant species also exhibit high survival in seawater (i.e. exhibit physiological plasticity), the survival in each salinity group was compared among the species by establishing the Cox proportional hazard models with the explanatory variable being species.

The Cox proportional hazard model assumes that the effects of the explanatory variables on mortality are constant regardless of time (Therneau and Grambsch, 2000). Accordingly, if the survival curves estimated by the Kaplan–Meier method (Kaplan and Meier, 1958) intersect between different experimental groups (i.e. inversion of mortality occurred), careful interpretation of results is necessary with respect to the survival time analysis. R statistical software version 4.1.0 (http://www.R-project.org/) was used for all analyses.

For the RNA-Seq analysis to compare transcriptional responses to different salinity conditions among the three species, basically, two individuals of 7 day old larvae were sampled from one rearing tank of the PA (G. petschiliensis), OO (G. opperiens), UE, UA (G. urotaenia in the rivers) and UB (G. urotaenia in the lake) populations in the above experimental conditions. This manipulation allowed us to gain 6–18 samples per species/population per salinity group in total. However, the actual number of successfully sampled individuals ranged from 3 to 14 due to mortality. After total RNA was isolated from these larval whole bodies, three samples with the highest RNA integrity number (4.0–8.5) for total RNA were selected for each salinity group in each population (30 samples in total). mRNA was isolated from the selected total RNA samples, and was converted into a cDNA library. Finally, libraries were sequenced by Macrogen Japan (Tokyo, Japan) using the Illumina HiSeq X sequencing platforms. This experiment generated a total of 14–38 million pair-end RNA-Seq reads per sample (Table S2). For a detailed protocol on the larval sampling and library preparation, see Supplementary Materials and Methods.

To identify key candidate genes involved in larval adaptation to freshwater, we investigated differentially expressed genes (DEGs) depending on salinity conditions in larvae of the three species. Prior to the differential expression analysis, the trimmed 4.2–10.9 million pair-end RNA-Seq reads from each of the 30 samples were mapped to a reference draft genome of Gymnogobius isaza (closely related to the target species; R.I., T.M., Y.H. and K.W., registration in progress with DNA Data Bank of Japan) using STAR version 2.7.4a (Dobin et al., 2013; mapping results are listed in Table S3). Subsequently, significant DEGs between freshwater- and seawater-reared groups were listed for each species with the edgeR package in R version 4.0.2 (Robinson et al., 2010) using TCC-GUI (Su et al., 2019). All the DEGs were functionally annotated by referring to the annotation information for the genomes of the zebrafish (Danio rerio) and threespine stickleback (Gasterosteus aculeatus) in the online translated nucleotide Basic Local Alignment Search Tool (TBLASTN; Altschul et al., 1990). Finally, the osmoregulatory importance of those DEGs was determined by analysing whether the list of Gene Ontology (GO) terms describing biological processes for each DEG included ion or water transport/homeostasis using the human genome database GeneCards version 5.7.0 (Stelzer et al., 2016). For a detailed protocol on read mapping and DEG screening, see Supplementary Materials and Methods.

To extract candidate genes involved in freshwater adaptation, we first focused on the significant DEGs shared by all three populations of G. urotaenia but not included in the DEG lists of the other species. We predicted that the functions of those DEGs were important for hyper-osmoregulation because of the high freshwater tolerance of G. urotaenia larvae (see Results).

Second, to identify candidate genes from other perspectives, we expanded our view to the entire list of significant DEGs (not limited to those found only in G. urotaenia). After identifying osmoregulatory DEGs using the GO term lists, we further extracted those which have been demonstrated to be responsible for osmoregulation in teleosts by exhaustively reviewing published research (see Results). For those genes of apparent importance in teleost osmoregulation, the magnitude of expression changes was compared between the species with and without freshwater forms. The differences in the magnitude were tested by investigating the significance of the interaction effect between salinity and the presence of freshwater forms in the species in linear mixed models (LMMs) using the lmer function (Bates et al., 2015). The variables in the models were as follows: the response variable: normalised read counts of each gene; the explanatory variables: salinity groups, presence/absence of freshwater forms within species and interaction term between these two explanatory variables; and the random effects: population ID and parent ID. We predicted that hyper- or hypo-osmoregulatory DEGs that were extracted by referring to the lists of GO terms and literature were more significantly up- or down-regulated in the freshwater-reared group of G. urotaenia than in the other species. The DEGs found only in any one of G. urotaenia populations were eliminated from this analysis to exclude the effects of local adaptation within the species.

We also examined the expression of the gene of FADS2, which decisively influences larval growth and survival in freshwater habitats in stickleback and other fishes through docosahexaenoic acid (DHA) synthesis (Ishikawa et al., 2019) as another candidate gene involved in freshwater adaptation.

The otolith microchemical analysis of adult fish from the targeted populations of the three Gymnogobius species revealed that all specimens of G. petschiliensis (n=7) and G. opperiens (n=12) showed high Sr/Ca ratios [5–8(×10^–3)] from near the core to the intermediate regions and low ratios (about 2×10^–3) in the outer regions, which aligns with the values in other amphidromous gobies [4–10(×10^–3) in the sea versus <4×10^–3 in freshwater areas (Tsunagawa and Arai, 2008; Iida et al., 2017)]. The results indicate that all specimens of these two species experienced amphidromous life histories (Fig. 3). Similarly, five of six individuals in the UE population of G. urotaenia had otoliths showing amphidromous Sr/Ca patterns (near the core: 5–7×10^–3; outer region: 2×10^–3; Fig. 3). As amphidromous Gymnogobius species are observed to migrate from the sea to freshwater areas while juveniles (Oto, 2021), the increase in Sr/Ca ratios is considered to reflect this juvenile migration. In contrast, otolith Sr/Ca ratios of one remaining UE individual and all the specimens from the UA (n=5) and UB (n=9) populations were constantly low (about 2×10^–3), suggesting that these individuals were freshwater forms.

The intraspecific comparison of larval survivorship between different salinity treatments revealed that the survival rate of G. petschiliensis was significantly higher in seawater than in freshwater (Cox proportional hazard model, z=−9.4, P<0.001; Fig. 4A). In contrast, G. opperiens and G. urotaenia exhibited significantly higher survival rates in freshwater than in seawater (z=10.0–17.3, P<0.001).

The interspecific comparison of survivorship in each salinity condition suggested that G. urotaenia larvae can survive well in the freshwater condition. Among the freshwater-reared groups, the survival rate of G. urotaenia larvae was significantly higher than that of G. petschiliensis and G. opperiens (Cox proportional hazard model, z=13.4–18.5, P<0.001; Fig. 4B). The magnitude of the difference in survival between the freshwater-reared G. petschiliensis and G. opperiens larvae could not be determined because of the intersection of the estimated survival curves (i.e. the occurrence of survival rate reversal). For the seawater groups, the survival rate was significantly lower in G. opperiens than in the two other species (z=−15.5–14.4, P<0.001; Fig. 4B), although this species spends the larval period in the sea, as indicated by the otolith microchemistry. The mortality of G. opperiens was especially high immediately after the transition from freshwater to seawater, perhaps due to their vulnerability to salinity changes. The magnitude of the difference in survival between seawater-reared G. petschiliensis and G. urotaenia larvae could not be determined because the estimated survival curves intersected. Intraspecific variation in larval survivorship in each salinity condition was sufficiently small compared with the variation between salinity groups or among species (Fig. 5).

Body size and starvation did not seem to critically affect survival rate because the notochord length and yolk sac size of G. opperiens (which showed poor survival in both salinity groups) and G. urotaenia (which displayed high survival in both salinity groups) were similar (Fig. S1; Table 1). Although G. petschiliensis larvae were slightly smaller than G. urotaenia larvae (Table 1), an unclear difference in clutch size between the two species is unlikely to affect the survival rate (Fig. S2).

The gene expression patterns of the three Gymnogobius species clearly differed depending on species rather than salinity conditions. The principal components analysis (PCA) for the RNA-Seq read count data of the 100 genes with the largest variance in read counts among the 30 samples represented three clusters which corresponded to each of the three species (for analysis procedures, see Supplementary Materials and Methods; Fig. 6A). The proportion of the variance explained by PC1 and PC2 was 44.0% and 24.6%, respectively. In G. urotaenia, the UE population tended to cluster separately from the other two geographically isolated populations (UA and UB) along PC2.

A total of 319 genes were detected as significant DEGs between the different salinity groups for at least one population (false discovery rate<0.1). Within these, 86 and 233 genes were up- and down-regulated in the freshwater group, respectively (Fig. 6B; Table S4). The number of the DEGs tended to be larger in G. urotaenia (246 and 30 in UE and UA stream populations, respectively, and 46 in UB landlocked lake population) compared with that in G. petschiliensis (9) and G. opperiens (33) (Fig. 6B). Of all the 319 significant DEGs, 265 were found in any one of the three G. urotaenia populations alone. In particular, the number of DEGs detected exclusively in the UE population amounted to 225, and 19 of the 225 DEGs contained osmoregulatory GO terms (e.g. slc5a1, kcna1a and atp1b1b; Table S4). However, those 265 DEGs were excluded from the detailed expression change analysis among species. The biological coefficient of variation (BCV), an index that can affect DEG detection (McCarthy et al., 2012), ranged from approximately 0.3 to 0.6 for each population. As the BCV values did not correlate with the number of DEGs (e.g. the greatest BCV in the UE with the most DEGs), the variation in the expression levels was not considered to significantly interfere with DEG detection.

The detected DEGs between the salinity groups included several osmoregulatory genes that were dramatically up-regulated in freshwater-reared G. urotaenia larvae. Most notable were two significant DEGs which were shared exclusively by the three populations of G. urotaenia. These two were candidate genes for freshwater adaptation. However, the expression change of one of the two (ca15b) was not considered to be significant only in G. urotaenia because the significant DEG list of G. opperiens also included another transcript similarly derived from ca15b (Table S4). The only significant DEG shared exclusively by G. urotaenia populations was the transcript of aqp3, encoding aquaporin-3 (AQP3), which is responsible for hyper-osmoregulation through the enhancement of basolateral hydraulic conductivity for cell volume regulation (Cutler et al., 2007; Cerdà and Finn, 2010; Whitehead et al., 2011).

Exhaustive screening of all the significant DEGs detected several genes responsible for hyper- or hypo-osmoregulation that were more significantly up- or down-regulated in the freshwater-reared larvae of G. urotaenia than in the other two species. After excluding the significant DEGs detected only from one G. urotaenia population, a total of 54 DEGs detected in any one (or more) species/population remained (intraspecific and interspecific variations in those expression levels can be simultaneously checked in Fig. S3), 49 of them being functionally annotated. These DEGs included the genes which have been demonstrated to be responsible for hyper-osmoregulation in teleosts encoding AQP3, solute carrier 13 (SLC13; Nakada et al., 2005) and solute carrier family 22 member 16 (SLC22A16; Lin et al., 2020) (Fig. 7; Table S4). Expression levels of these two genes were higher in the freshwater-reared G. urotaenia than in the other two species, but were similarly low in the seawater groups of all species. This was supported by the significant interaction effect between salinity condition and species' migratory pattern (presence of freshwater forms) on the expression levels of those genes (LMM, P<0.006; Table 2), justifying that those two were candidate genes for freshwater adaptation.

The 49 DEGs also included the genes which have been demonstrated to be responsible for hypo-osmoregulation in teleosts. The genes encode serum and glucocorticoid-induced protein kinase-1 (SGK-1; Shaw et al., 2008) and Na⁺/K⁺/Cl⁻ cotransporter (NKCC; Flemmer et al., 2010). In these genes, the gene expression of NKCC was significantly more down-regulated in freshwater-reared G. urotaenia than in the other species (LMM, P=0.023). The DEG list also included the gene for the claudin (CLDN) family, which is involved in hyper- or hypo-osmoregulation through tight junction pathways (Tipsmark et al., 2008; Guo et al., 2018). Its expression was more up-regulated in freshwater groups of G. urotaenia than in the other species (LMM, P=0.029). The genes encoding NKCC and CLDN were also candidate genes for freshwater adaptation. None of these DEGs, which are essential for teleost osmoregulation, exhibited greater expression changes in G. petschiliensis and G. opperiens than in G. urotaenia.

The gene expression of FADS2 was much lower in G. opperiens than in the other species, but did not significantly change between salinity groups in any species (Fig. 8).

The main objective of this study was to test hypothesis that G. urotaenia, which has established freshwater forms, has evolved a higher freshwater tolerance at the larval stage than the congeneric amphidromous species. As expected, G. urotaenia was confirmed to flexibly establish freshwater-resident forms in both lake and streams and exhibit high larval survival in the freshwater condition. Furthermore, the species showed greater changes in the expression levels of some important osmoregulatory genes than the obligate amphidromous G. petschiliensis and G. opperiens. The present results obtained by multifaceted experimental approaches provide the first empirical evidence supporting the importance of the evolution of larval freshwater tolerance for freshwater colonisation in diadromous fish. Moreover, the flexibility of osmoregulatory gene expression in G. urotaenia suggests that the enhancement of osmoregulatory plasticity is involved in the acquisition of freshwater tolerance, as discussed in more detail below.

The otolith microchemical analysis confirmed variation in the occurrence of freshwater forms (which remain in freshwater areas throughout life) among the three amphidromous Gymnogobius species. All G. petschiliensis and G. opperiens specimens exhibited amphidromous characteristics, with the hatchlings migrating from streams to the sea and returning to freshwater in the juvenile period. This result is consistent with the fact that no landlocked populations have previously been reported for these two species (Stevenson, 2002; Ishino, 2005). Similarly, the majority of G. urotaenia specimens from the Eiheiji River (UE) exhibited the amphidromous pattern. However, either all or some of the G. urotaenia specimens from all three populations were estimated to have remained in freshwater areas (both streams and lakes) throughout their life. This result is consistent with previous findings that several putatively freshwater populations exist in this species (Stevenson, 2002; Ishino, 2005). Notably, the individual variation in the migratory history of the UE population suggests that G. urotaenia can stay in freshwater habitats without being trapped in freshwater areas by physical barriers (e.g. naturally or artificially dammed lakes). This species may maintain the intra-populational polymorphism of migration as the alternative life history strategy. The lack of female samples in this study necessitates the further examination of migration patterns in both sexes, although there have been no previous reports confirming sexual differences in the migration patterns in gobies.

The rearing experiment supported the hypothesis that larval fish of G. urotaenia exhibit higher freshwater tolerance than the other two species. This physiological trait would partly explain the frequent occurrence of freshwater forms of this species.

As all individual larval fish in the experiment had died by 16 days after hatching, this study could not observe survival throughout the entire larval period, which lasts more than a month in their natural habitat. In addition, the higher mortality of G. opperiens in both freshwater and seawater conditions than the other species suggests that this species was particularly vulnerable to the local rearing conditions. It cannot be entirely ruled out that this species-specific larval debility may have affected the evaluation of physiological plasticity and gene expression changes. One possible reason for the mortality is the nutritional problem of the frozen rotifer, although its use was unavoidable to ensure uniform feeding conditions between different salinity groups. Nevertheless, the experiments under common conditions were judged to have effectively detected the interspecific differences in osmoregulatory ability in freshwater environments for the following two reasons: (1) each clutch was divided into freshwater and seawater groups; thus, most of the larval debility irrelevant to osmotic stress would have been cancelled out by this intra-clutch comparison; and (2) the analysis of 8–16 different clutches from 2–3 populations per species would have allowed for accurate interspecific comparisons.

The gene expression analysis provided evidence supporting that the evolution of gene expression associated with hyper-osmoregulation improved freshwater tolerance of G. urotaenia larvae. Most of the DEGs responsible for hyper-osmoregulation were more significantly up-regulated in freshwater-reared G. urotaenia than in the other species. One of the up-regulated genes was that encoding AQP3 (aqp3), a membrane water channel responsible for water discharge from the epithelial cells under hypotonic conditions (Cerdà and Finn, 2010). In adult fish, it is well evidenced that the enhanced expression of aqp3 is important for acclimation to freshwater in marine-originated species (for example, the killifish Fundulus heteroclitus: Whitehead et al., 2011; the alewife Alosa pseudoharengus: Velotta et al., 2017). Remarkably, aqp3 is practically the only significant DEG shared by all three G. urotaenia populations but not by the other species, which makes aqp3 a particularly strong candidate gene for larval freshwater tolerance. Another was the gene encoding SLC13, which contributes to the maintenance of sulphate homeostasis in freshwater-reared eels (Nakada et al., 2005). Similarly, G. urotaenia had a more pronounced up-regulation of the gene encoding CLDN, which can contribute to hyper-osmoregulation (Tipsmark et al., 2008; Guo et al., 2018). The enhanced expression of these hyper-osmoregulatory genes in larval fish would strengthen the freshwater tolerance of G. urotaenia, possibly facilitating landlocking or colonisation in freshwater areas. The importance of the evolution of hyper-osmoregulatory gene expression for freshwater residence has also been indicated in certain adult diadromous fish species (Velotta et al., 2015, 2017; Kusakabe et al., 2017). The present study provides the first evidence of gene expression evolution of larval fish in a lineage colonising freshwater habitats.

The down-regulation of gene expression for water absorption and ion discharge (hypo-osmoregulation) can also be important for freshwater acclimation in teleosts (Velotta et al., 2015, 2017; Nakamura et al., 2020). In our study, only one hypo-osmoregulatory gene (the gene encoding NKCC) was more significantly down-regulated in freshwater-reared G. urotaenia than in the other species. Although we cannot rule out the possibility that some important hypo-osmoregulatory genes may have been missed, the up-regulation of hyper-osmoregulatory genes may be particularly important for freshwater acclimation in Gymnogobius larvae.

The results of the gene expression analysis further suggest the importance of physiological plasticity for establishing freshwater forms at both individual and population levels. High physiological plasticity is widely considered to be essential for colonising novel environments (Ghalambor et al., 2007; Lande, 2015). For instance, high osmoregulatory plasticity has possibly facilitated freshwater colonisation in the threespine stickleback (McCairns and Bernatchez, 2010). In the present study, G. urotaenia larvae appeared to have higher osmoregulatory plasticity than larvae of the other species, as they could largely down-regulate the expression of hyper-osmoregulatory genes (e.g. aqp3) in seawater. The high osmoregulatory plasticity was also suggested by the greater number of significant DEGs between salinity conditions in G. urotaenia (30, 46 and 246) than in G. petschiliensis (8) and G. opperiens (33).

The survival analysis did not reject the possibility that the osmoregulatory plasticity was higher in G. urotaenia than in the other species, as evidenced by the high survival rates of G. urotaenia in the freshwater condition, which was harsh for larval fish of the marine-originated amphidromous species. In particular, the UE population of G. urotaenia, which had the greatest number of DEGs (246), exhibited the highest survival rate in both freshwater and seawater. In addition, some of the DEGs detected only in the UE population were related to osmoregulation. This population included both amphidromous and freshwater-resident individuals, suggesting that physiological plasticity may play an important role in maintaining such a life-history polymorphism. In contrast, the survival of G. urotaenia slightly decreased in the seawater condition, suggesting a trade-off between hyper- and hypo-osmoregulation, as observed in other fish species (Velotta et al., 2015). Thus, struggling with the osmoregulatory trade-off, larval fish of G. urotaenia may have evolutionarily enhanced physiological plasticity, so that they can flexibly choose a freshwater-resident life as an alternative strategy or establish freshwater populations when landlocked by geographic barriers.

The reason the larval habitats of marine-originated amphidromous species are usually limited to the sea (i.e. ancestral environment) might be explained by their much lower osmoregulatory ability than that of the adult fish, as osmoregulatory organs (gills, kidneys and intestines) are poorly developed at this early life stage (Kaneko et al., 2002; Varsamos et al., 2005). Interspecific variation in the gene expression patterns among the Gymnogobius larvae was much greater than intraspecific variation between the salinity groups, implying that the number of molecular mechanisms involved in larval freshwater adaptation is not large. As in the case of Gymnogobius species, aqp3 is already expressed at the early larval stage in some fish species (Deane and Woo, 2006; Cerdà and Finn, 2010). The evolution of expression levels of such a few specific genes may be crucial to enhancing the physiological plasticity of early-stage larval fish.

In addition to our findings, there could be other molecular mechanisms which enhance larval osmoregulatory ability in freshwater conditions. Other studies analysing the transcriptomes in osmoregulatory organs (e.g. gill and intestine) detected >100 DEGs between freshwater- and seawater-reared adults of anadromous species of alewife and puffer fish (Velotta et al., 2017; Nakamura et al., 2020). The present study detected smaller numbers of DEGs between the salinity groups of larval fish from each of the Gymnogobius species (9–246, mostly <50). This could be partly because the osmoregulatory organs of early larval fish are undeveloped as discussed above. Another possibility is that the efficiency of detecting DEGs was decreased by the highly expressed transcripts of non-osmoregulatory tissues in the whole-body transcriptome analysis and the small number of biological replicates (n=3). Further studies should therefore focus on direct quantification of the expression levels of commonly investigated genes controlling hyper-osmoregulation (e.g. Na⁺, Cl⁻ cotransporter and Na⁺/H⁺ exchanger; Edwards and Marshall, 2013).

Factors other than the modification of osmoregulatory systems could also be essential for freshwater colonisation by larval fish. One known example is the duplication and enhanced expression level of the gene encoding FADS2 (fads2) in marine-originated fishes, which increase their ability to synthesise DHA, dramatically improving larval survival in freshwater (Ishikawa et al., 2019). In the case of Gymnogobius, G. petschiliensis larvae with low freshwater tolerance showed a constant, moderate fads2 expression. Therefore, the DHA synthesis ability is unlikely to be a common restriction factor for freshwater colonisation. However, this might partly explain the absence of landlocked populations of G. opperiens within their natural environments, as the larvae exhibited fairly low fads2 expression.

In this study, otolith microchemical analysis, physiological experimentation and transcriptome analyses consistently supported the prediction that enhanced freshwater tolerance of early larvae is crucial in establishing freshwater forms at both individual and population levels in a marine-originated amphidromous group (Gymnogobius gobies). Moreover, interspecific comparison of the expression changes of osmoregulatory genes suggested the importance of osmoregulatory plasticity in acquiring freshwater tolerance. As the genus Gymnogobius includes several obligate freshwater species that invade various habitats (Tabata and Watanabe, 2013; Chiba et al., 2020), the evolution of larval freshwater tolerance through the enhancement of osmoregulatory plasticity may have played an important role in the ecological diversification of this genus. This study also found a strong candidate gene, aqp3, responsible for freshwater tolerance of Gymnogobius larvae, which is one of the few osmoregulatory genes already expressed during early developmental stages in teleost species (Deane and Woo, 2006; Cerdà and Finn, 2010). Further experiments are needed to test whether the evolution of expression of aqp3 and other similar candidate genes is responsible for larval freshwater colonisation, by manipulating their expression.

We thank K. Shirai, Y. Kusaba, K. Kido and A. Goto for their help with the otolith microchemical analysis, S. Kunimatsu for his help with the molecular species identification, M. Nakamura and Y. Morii for their advice on the manuscript, and members of Laboratory of Animal Ecology of Kyoto University for their helpful comments on the study design. The Atmosphere and Ocean Research Institute of the University of Tokyo provided an electron probe microanalyser for the otolith microchemical analysis. The larval rearing room was provided by the Graduate School of Science, Kyoto University.