TY - JOUR
AU - Greenway, Ryan
AB - Introduction Salinity variation imposes osmoregulatory challenges on aquatic organisms, and contact zones between freshwater and saltwater environments act as barriers that limit the ability of animals to move from one habitat to the other [1]. Many aquatic taxa have consequently failed to cross natural salinity gradients [2]. Those that achieve such habitat shifts overcome osmoregulatory challenges through plasticity or adaptation, and both responses have greatly shaped aquatic species distributions [3–5]. Due to the ecological expansion that accompanies colonization of novel habitats, invasions from freshwater to saltwater environments, or vice versa, that result in adaptation are of particular interest for elucidating the genetic mechanisms underlying salinity tolerance and the diversification of aquatic organisms [1, 2, 6]. Salinity transitions in fishes Among fishes, diversification has coincided with repeated transitions between freshwater and saltwater habitats [2]. Many saltwater-freshwater transitions in fishes have occurred over long evolutionary timescales, and deep evolutionary divergences have resulted in many species only tolerating a narrow range of salinities, restricting them to either freshwater or saltwater environments [2, 6–9]. A considerably smaller portion of fish species can survive in both freshwater and saltwater environments [10]. Among species that tolerate a broad range of salinities, movement along salinity clines is characteristic of diadromous lineages with life histories involving migration between freshwater and saltwater environments during an individual’s lifetime [10, 11]. While some species cannot cross the saltwater-freshwater boundary and some do so routinely in their lifetimes, there are also lineages between these two extremes that have made transitions between saltwater and freshwater environments at microevolutionary scales [6, 12, 13]. Few studies have investigated mechanisms of salinity adaptation in species where closely related populations experience different salinity regimes, and the evolutionary repeatability of these mechanisms across lineages remains to be explored. Transitions between freshwater and saltwater environments are challenging for animals that actively maintain internal solute homeostasis. The many physiological processes involved in stable state osmotic and ionic balance necessitate changes in multiple interdependent processes when crossing a salinity barrier [14]. To maintain homeostasis, fishes in freshwater must actively absorb salt and excrete water in the form of dilute urine to counteract their passive loss of salt and absorption of water [14]. In contrast, fishes in saltwater environments must remove salt and retain water [14]. As a result, crossing a salinity barrier requires a shift between absorption and excretion of ions and water in multiple organs, involving both active and passive processes [14, 15]. Remodeling of gill epithelia and regulation of ion transporters, aquaporins, and tight junctions in the gill are particularly central to this process [15–18]. Gene expression responses to salinity variation One way to quantify such complex physiological responses to variation in salinity is to compare patterns of gene expression across populations in different habitat types. Several studies have investigated the physiological and transcriptomic responses to changes in salinity between populations of the same species or among closely related species of fish [13, 18–20]. Differential gene expression associated with osmoregulation and ion transport is commonly found between populations inhabiting environments of different salinities [13, 18, 20]. In addition, differential expression has also been documented in genes associated with other biological functions, including immune processes [19, 20], cell communication [13], stress tolerance [13], and gill membrane permeability [19]. These studies highlight candidate pathways that play important roles in divergence between saltwater and freshwater ecotypes within species, but there have been few comparisons investigating the repeatability of gene expression responses to variation in salinity across phylogenetically disparate taxa. Investigating evidence of convergence in gene expression patterns across taxa that have undergone similar salinity transitions will provide insight into possible shared and unique pathways involved in adaptation to different salinity regimes. Salinity tolerance in Limia perugiae (Poeciliidae) In this study, we investigated patterns of gene expression in a freshwater and a hypersaline population of a livebearing fish species, Limia perugiae (Poeciliidae), and compared transcriptomic variation between these populations to those observed in other species that have undergone similar salinity transitions. Freshwater fishes in the genus Limia are endemic to the islands of the West Indies [21, 22]. Limia perugiae is a widespread species across the southern portion of Hispaniola, occurring in freshwater artesian springs and low-order creeks, as well as hypersaline inland lakes and coastal lagoons [23–25]. Exposure to high salinities in L. perugiae has been shown to decrease metabolic rate [25], increase the production of Na+/K+-ATPase and oxidative phosphorylation proteins in the gills [23], and reduce adult body size [23]. Though predominantly associated with freshwater habitats, many fishes in the family Poeciliidae are able to tolerate a broad range of salinities, a factor potentially responsible for facilitating their dispersal across a wide geographic range [26–29]. While the mechanisms and consequences of salinity tolerance at the biochemical and physiological levels have been a focus of research in poeciliids [23, 30–33], the genetic and regulatory underpinnings of salinity tolerance have yet to be investigated in this family. We used a natural system with conspecific populations occurring in both a freshwater and hypersaline habitat to characterize the potential molecular mechanisms underlying high salinity tolerance in poeciliid fishes. Objectives We used RNA-sequencing (RNA-seq) to study genome-wide patterns of gene expression between freshwater and hypersaline L. perugiae. From this analysis, we aimed to identify genes and physiological pathways associated with salinity tolerance in this species. We then compared the L. perugiae population pair to other population pairs of freshwater and saltwater ecotypes in disparate actinopterygian taxa to understand if mechanisms of osmoregulatory capability are shared across divergent lineages of teleost fishes. We utilized a comparative transcriptomics approach that leverages new and pre-existing gene expression datasets to address the following questions: 1) What genes are differentially expressed between freshwater and saltwater L. perugiae populations, and with what physiological processes are they associated? 2) Is there evidence for commonalities in gene expression among phylogenetically disparate teleosts with freshwater and saltwater populations? Salinity transitions in fishes Among fishes, diversification has coincided with repeated transitions between freshwater and saltwater habitats [2]. Many saltwater-freshwater transitions in fishes have occurred over long evolutionary timescales, and deep evolutionary divergences have resulted in many species only tolerating a narrow range of salinities, restricting them to either freshwater or saltwater environments [2, 6–9]. A considerably smaller portion of fish species can survive in both freshwater and saltwater environments [10]. Among species that tolerate a broad range of salinities, movement along salinity clines is characteristic of diadromous lineages with life histories involving migration between freshwater and saltwater environments during an individual’s lifetime [10, 11]. While some species cannot cross the saltwater-freshwater boundary and some do so routinely in their lifetimes, there are also lineages between these two extremes that have made transitions between saltwater and freshwater environments at microevolutionary scales [6, 12, 13]. Few studies have investigated mechanisms of salinity adaptation in species where closely related populations experience different salinity regimes, and the evolutionary repeatability of these mechanisms across lineages remains to be explored. Transitions between freshwater and saltwater environments are challenging for animals that actively maintain internal solute homeostasis. The many physiological processes involved in stable state osmotic and ionic balance necessitate changes in multiple interdependent processes when crossing a salinity barrier [14]. To maintain homeostasis, fishes in freshwater must actively absorb salt and excrete water in the form of dilute urine to counteract their passive loss of salt and absorption of water [14]. In contrast, fishes in saltwater environments must remove salt and retain water [14]. As a result, crossing a salinity barrier requires a shift between absorption and excretion of ions and water in multiple organs, involving both active and passive processes [14, 15]. Remodeling of gill epithelia and regulation of ion transporters, aquaporins, and tight junctions in the gill are particularly central to this process [15–18]. Gene expression responses to salinity variation One way to quantify such complex physiological responses to variation in salinity is to compare patterns of gene expression across populations in different habitat types. Several studies have investigated the physiological and transcriptomic responses to changes in salinity between populations of the same species or among closely related species of fish [13, 18–20]. Differential gene expression associated with osmoregulation and ion transport is commonly found between populations inhabiting environments of different salinities [13, 18, 20]. In addition, differential expression has also been documented in genes associated with other biological functions, including immune processes [19, 20], cell communication [13], stress tolerance [13], and gill membrane permeability [19]. These studies highlight candidate pathways that play important roles in divergence between saltwater and freshwater ecotypes within species, but there have been few comparisons investigating the repeatability of gene expression responses to variation in salinity across phylogenetically disparate taxa. Investigating evidence of convergence in gene expression patterns across taxa that have undergone similar salinity transitions will provide insight into possible shared and unique pathways involved in adaptation to different salinity regimes. Salinity tolerance in Limia perugiae (Poeciliidae) In this study, we investigated patterns of gene expression in a freshwater and a hypersaline population of a livebearing fish species, Limia perugiae (Poeciliidae), and compared transcriptomic variation between these populations to those observed in other species that have undergone similar salinity transitions. Freshwater fishes in the genus Limia are endemic to the islands of the West Indies [21, 22]. Limia perugiae is a widespread species across the southern portion of Hispaniola, occurring in freshwater artesian springs and low-order creeks, as well as hypersaline inland lakes and coastal lagoons [23–25]. Exposure to high salinities in L. perugiae has been shown to decrease metabolic rate [25], increase the production of Na+/K+-ATPase and oxidative phosphorylation proteins in the gills [23], and reduce adult body size [23]. Though predominantly associated with freshwater habitats, many fishes in the family Poeciliidae are able to tolerate a broad range of salinities, a factor potentially responsible for facilitating their dispersal across a wide geographic range [26–29]. While the mechanisms and consequences of salinity tolerance at the biochemical and physiological levels have been a focus of research in poeciliids [23, 30–33], the genetic and regulatory underpinnings of salinity tolerance have yet to be investigated in this family. We used a natural system with conspecific populations occurring in both a freshwater and hypersaline habitat to characterize the potential molecular mechanisms underlying high salinity tolerance in poeciliid fishes. Objectives We used RNA-sequencing (RNA-seq) to study genome-wide patterns of gene expression between freshwater and hypersaline L. perugiae. From this analysis, we aimed to identify genes and physiological pathways associated with salinity tolerance in this species. We then compared the L. perugiae population pair to other population pairs of freshwater and saltwater ecotypes in disparate actinopterygian taxa to understand if mechanisms of osmoregulatory capability are shared across divergent lineages of teleost fishes. We utilized a comparative transcriptomics approach that leverages new and pre-existing gene expression datasets to address the following questions: 1) What genes are differentially expressed between freshwater and saltwater L. perugiae populations, and with what physiological processes are they associated? 2) Is there evidence for commonalities in gene expression among phylogenetically disparate teleosts with freshwater and saltwater populations? Materials and methods Sample collection Limia perugiae were collected using a seine from a hypersaline lagoon (Laguna Oviedo: 17.801°N, 71.363°W) and a geographically proximate freshwater stream (Los Cocos: 17.905°N, 71.286°W) in the Dominican Republic. Laguna Oviedo had a salinity of 54–61 ppt, and Los Cocos had a salinity of 0–1 ppt. Following capture, adult females (N = 6 per site) were euthanized using MS222, and all efforts were taken to alleviate suffering during handling and euthanasia. Both sets of gills were extracted using sterilized scissors and forceps. Tissues were immediately preserved in RNAlater (Ambion Inc., Austin, TX, USA). All samples were collected with permits from the corresponding authorities in the Dominican Republic (permit number 0092–11). All procedures involving animals followed established best-practices and were approved by the Institutional Animal Care and Use Committee of Kansas State University (protocol 3473). RNA-seq library preparation RNA extraction, library preparation, and sequencing of samples followed procedures previously employed for related poeciliid species [34, 35]. Approximately 10–30 mg of tissue from each individual was sealed in a Covaris TT1 TissueTube (Covaris, Inc., Woburn, MA, USA), frozen in liquid nitrogen, and pulverized. Total RNA was extracted from the pulverized tissue using the NucleoSpin RNA kit (Machery-Nagel, Düren, Germany). mRNA isolation and cDNA library preparation were completed with the NEBNext Poly(A) mRNA Magnetic Isolation Module (New England Biolabs, Inc., Ipswich, MA, USA) and NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England Biolabs, Inc., Ipswich, MA, USA), with minor modifications to the manufacturers’ protocol [34–36]. cDNA libraries were individually barcoded, quantified with Qubit and an Agilent 2100 Bioanalyzer High Sensitivity DNA chip, and then pooled with cDNA samples from other projects in sets of 11–12 samples based on nanomolar concentrations. Samples were split across pools such that samples from each habitat type were not all sequenced together, and there was no evidence for lane effects. Libraries were sequenced on an Illumina HiSeq 2500 using paired-end 101-base-pair (bp) reads at the Washington State University Spokane Genomics Core. Mapping All raw reads were trimmed twice (quality 0 to remove Illumina adapters, followed by quality 24) using Trimgalore! v.0.4.0 [37]. Trimmed reads were mapped to the Poecilia mexicana reference genome (RefSeq accession number: GCF_001443325.1 [38]) with an appended mitochondrial genome (GenBank Accession Number: KC992998.1) using BWA-MEM v.0.7.12 [39]. We annotated genes from the P. mexicana reference genome by extracting the longest transcript for each gene (with the perl script gff2fast.pl: https://github.com/ISUgenomics/common_scripts/blob/master/gff2fasta.pl) and comparing them against entries in the human SWISS-PROT database (critical E-value 0.001; access date 04/15/2017) using BLASTx [40]. Each gene was annotated with the best BLAST hit from the human database based on the top high-scoring segment pair. Differential gene expression We used StringTie (v.1.3.3b) [41, 42] to quantify the number of reads mapped to each gene for each individual (measured in counts per million mapped reads) and then used the prepDE.py script (provided with StringTie) to generate a read count matrix [42]. We removed genes that did not have at least two counts per million in three or more individuals across both populations, resulting in 18,659 genes that were included in differential gene expression analysis. We identified differentially expressed genes using generalized linear models (GLMs) in R, as implemented in the Bioconductor package edgeR [43]. We fit a negative binomial GLM to the normalized read counts of each gene based on tagwise dispersion estimates and a design matrix describing the comparison between the saltwater and freshwater population using glmFit. The tagwise dispersion estimates were generated using the estimateDisp function in edgeR, which employs a weighted likelihood empirical Bayes approach [44]. We assessed statistical significance using the GLM likelihood-ratio test with a false discovery rate (FDR) of q-value < 0.05, calculated with the Benjamini-Hochberg correction [45]. After identifying the set of differentially expressed genes between the saltwater and freshwater population, we used a Gene Ontology (GO) enrichment analysis to explore putative biological functions of these genes. We annotated all differentially expressed genes that had a match in the human SWISS-PROT database with GO IDs [46] and tested for the enrichment of specific GO IDs separately in up and downregulated genes relative to the full background set of 18,659 genes using GOrilla (FDR q-value < 0.05, accessed May 26, 2022) [47]. A total of 10,935 genes in the background set were associated with a GO term in the database. Weighted gene co-expression network analysis We constructed weighted gene co-expression networks to identify clusters of genes that were co-expressed across our samples [48]. To prepare the gene expression data for this analysis, we applied a variance-stabilizing transformation to the filtered reads using the varianceStabilizingTransformation function from the DESeq2 package (v.1.36.0) [49] in R, which allowed us to normalize the read counts relative to library size and appropriately scale the data for clustering [50, 51]. After transforming and normalizing the read counts, we used the cpm function in the edgeR package (v.3.38.1) [43] to generate a gene matrix of scaled read counts (log2-cpm, counts per million mapped reads) from the transformed read counts [52]. As outlined in the Weighted Gene Co-expression Network Analysis (WGCNA) package documentation [50], we used hierarchical clustering to cluster the samples based on their gene expression profiles. We used the hclust function from the flashClust package (v.1.1.2) [51] to cluster the samples. We then created a weighted network adjacency matrix using the adjacency function from the WGCNA package (v.1.71) [50, 51]. The adjacency matrix was constructed by calculating pairwise co-expression similarities (Pearson correlation coefficients) and raising them to a power of β, a soft thresholding power [50]. We used the pickSoftThreshold function in the WGCNA package [50] to assist in selecting the lowest possible value of β that ensured our network fit the approximate scale free topology criterion while retaining the highest possible mean connectivity between the network genes. Based on the scale free topology model fit and mean connectivity of our network, and following the recommendations outlined by Zhang and Horvath [48], we selected β = 7. From our correlation network, we then generated a topological overlap dissimilarity matrix to identify modules of co-expressed genes. We calculated dissimilarity between the genes by converting the adjacencies into topological overlap similarities using the TOMsimilarity function in the WGCNA package [50] and then subtracting these topological overlap measures from 1. To identify modules of co-expressed genes, we used the hclust function from the flashClust package [51] for hierarchical clustering of the genes and then created a hierarchical clustering dendrogram. We used the cutreeDynamic function from the dynamicTreeCut package [50, 53] to extract the modules from the dendrogram. To summarize the gene expression variation in each module, the first principal component of each module in the expression matrix (the eigengene) was calculated using the moduleEigengenes function from the WGCNA package [50]. We used the function mergeCloseModules in the WGCNA package to merge eigengenes that were highly correlated. We included eigengenes with a correlation greater than 0.9 for merging. Finally, we identified modules that were significantly associated with the presence or absence of salinity by calculating correlation coefficients between the eigengenes and the habitat type. P-values of the correlation coefficients were calculated using the corPvalueStudent function from the WGCNA package, and we retained modules with P-values less than 0.01 for functional enrichment analyses. Similar to our differential gene expression analysis, we used GOrilla [47] for Gene Ontology enrichment analysis of genes contained in modules exhibiting significant correlations with habitat type. Comparisons of L. perugiae with other species We mined previously published datasets to identify gene expression patterns commonly associated with salinity tolerance in disparate taxa, including South American silversides (Odontesthes bonariensis and Odontesthes argentinensis) from Hughes et al. [20], three-spine sticklebacks (Gasterosteus aculeatus) from Gibbons et al. [19], Amur ide (Leuciscus waleckii) from Xu et al. [13], and L. perugiae (this study; Table 1). Each of these experiments generated paired-end RNA-seq raw reads from gill tissue in two ecotypes (one freshwater and the other saltwater) of the same lineage. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. Species included in the analysis, including environment (freshwater [FW] or saltwater [SW]), collection location, sample size (N), NCBI Sequence Read Archive (SRA) accession numbers, and study reference. https://doi.org/10.1371/journal.pone.0315014.t001 Raw RNA-seq reads from each transcriptomics project were downloaded in FASTQ format, and reference genomes or transcriptomes for each species were downloaded in FASTA format from Genbank (see Table 1 for accession numbers). Reads were trimmed and mapped to their respective reference genomes (Poecilia mexicana: GCF_001443325.1 [38]; Cyprinus carpio: GCF_000951615.1 [54]; Gasterosteus aculeatus: Broad S1 v. 93 [55]) or reference transcriptome (Menidia menidia: GEVY00000000.1 [56]) following the same methods described above for L. perugiae. We then quantified the number of reads mapped to each gene in the annotation file for each reference genome and created a read counts matrix for each species, which were used for subsequent expression analyses. Expression analyses were performed in R version 4.1.2. The 10,000 genes with the highest standard deviation between freshwater and saltwater samples were abstracted from each read counts matrix, and overall expression patterns were visualized with multi-dimensional scaling (MDS) plots. To make comparisons across species, we used OrthoFinder v2.2.6 to identify orthologous genes among the reference genomes [57, 58]. For OrthoFinder, we used the ‘dendroblast’ option for gene tree inference, ‘blast’ for the sequence search program, ‘mafft’ for the multiple sequence alignment, and ‘fasttree’ for the tree inference method. A total of 18,419 orthogroups were identified. Of the total orthogroups identified, 13,899 orthogroups had at least one copy in each species. Of those, 1,735 were 1:1 orthologs. To calculate counts per orthogroup, we used the gene counts matrix of each species to sum up the counts across all loci contained in an orthogroup. Based on this orthogroup counts matrix, we retained only the orthogroups that were present in all species and expressed in all individuals (cpm > 0 per individual), resulting in 12,743 retained orthogroups. To evaluate expression differences for each orthogroup, we made pairwise comparisons between ecotypes following the same methods described above for the L. perugiae comparisons. Briefly, we normalized reads, created and compared generalized linear models of the normalized read counts, generated a design matrix, estimated tagwise dispersion, and conducted GLM likelihood-ratio tests to test whether differences in expression were statistically different between the freshwater and saltwater population for each orthgroup. To identify orthogroups exhibiting convergent expression patterns across lineages, we intersected the significantly upregulated and downregulated orthogroups from all lineage-specific comparisons, identifying orthogroups that were differentially expressed in the same direction in pairwise, three-, and four-way comparisons among the lineages. After identifying the set of orthogroups with differential expression across three or more lineages, we used a GO enrichment analysis as described above to explore the putative biological functions of these candidate gene sets. Sample collection Limia perugiae were collected using a seine from a hypersaline lagoon (Laguna Oviedo: 17.801°N, 71.363°W) and a geographically proximate freshwater stream (Los Cocos: 17.905°N, 71.286°W) in the Dominican Republic. Laguna Oviedo had a salinity of 54–61 ppt, and Los Cocos had a salinity of 0–1 ppt. Following capture, adult females (N = 6 per site) were euthanized using MS222, and all efforts were taken to alleviate suffering during handling and euthanasia. Both sets of gills were extracted using sterilized scissors and forceps. Tissues were immediately preserved in RNAlater (Ambion Inc., Austin, TX, USA). All samples were collected with permits from the corresponding authorities in the Dominican Republic (permit number 0092–11). All procedures involving animals followed established best-practices and were approved by the Institutional Animal Care and Use Committee of Kansas State University (protocol 3473). RNA-seq library preparation RNA extraction, library preparation, and sequencing of samples followed procedures previously employed for related poeciliid species [34, 35]. Approximately 10–30 mg of tissue from each individual was sealed in a Covaris TT1 TissueTube (Covaris, Inc., Woburn, MA, USA), frozen in liquid nitrogen, and pulverized. Total RNA was extracted from the pulverized tissue using the NucleoSpin RNA kit (Machery-Nagel, Düren, Germany). mRNA isolation and cDNA library preparation were completed with the NEBNext Poly(A) mRNA Magnetic Isolation Module (New England Biolabs, Inc., Ipswich, MA, USA) and NEBNext Ultra Directional RNA Library Prep Kit for Illumina (New England Biolabs, Inc., Ipswich, MA, USA), with minor modifications to the manufacturers’ protocol [34–36]. cDNA libraries were individually barcoded, quantified with Qubit and an Agilent 2100 Bioanalyzer High Sensitivity DNA chip, and then pooled with cDNA samples from other projects in sets of 11–12 samples based on nanomolar concentrations. Samples were split across pools such that samples from each habitat type were not all sequenced together, and there was no evidence for lane effects. Libraries were sequenced on an Illumina HiSeq 2500 using paired-end 101-base-pair (bp) reads at the Washington State University Spokane Genomics Core. Mapping All raw reads were trimmed twice (quality 0 to remove Illumina adapters, followed by quality 24) using Trimgalore! v.0.4.0 [37]. Trimmed reads were mapped to the Poecilia mexicana reference genome (RefSeq accession number: GCF_001443325.1 [38]) with an appended mitochondrial genome (GenBank Accession Number: KC992998.1) using BWA-MEM v.0.7.12 [39]. We annotated genes from the P. mexicana reference genome by extracting the longest transcript for each gene (with the perl script gff2fast.pl: https://github.com/ISUgenomics/common_scripts/blob/master/gff2fasta.pl) and comparing them against entries in the human SWISS-PROT database (critical E-value 0.001; access date 04/15/2017) using BLASTx [40]. Each gene was annotated with the best BLAST hit from the human database based on the top high-scoring segment pair. Differential gene expression We used StringTie (v.1.3.3b) [41, 42] to quantify the number of reads mapped to each gene for each individual (measured in counts per million mapped reads) and then used the prepDE.py script (provided with StringTie) to generate a read count matrix [42]. We removed genes that did not have at least two counts per million in three or more individuals across both populations, resulting in 18,659 genes that were included in differential gene expression analysis. We identified differentially expressed genes using generalized linear models (GLMs) in R, as implemented in the Bioconductor package edgeR [43]. We fit a negative binomial GLM to the normalized read counts of each gene based on tagwise dispersion estimates and a design matrix describing the comparison between the saltwater and freshwater population using glmFit. The tagwise dispersion estimates were generated using the estimateDisp function in edgeR, which employs a weighted likelihood empirical Bayes approach [44]. We assessed statistical significance using the GLM likelihood-ratio test with a false discovery rate (FDR) of q-value < 0.05, calculated with the Benjamini-Hochberg correction [45]. After identifying the set of differentially expressed genes between the saltwater and freshwater population, we used a Gene Ontology (GO) enrichment analysis to explore putative biological functions of these genes. We annotated all differentially expressed genes that had a match in the human SWISS-PROT database with GO IDs [46] and tested for the enrichment of specific GO IDs separately in up and downregulated genes relative to the full background set of 18,659 genes using GOrilla (FDR q-value < 0.05, accessed May 26, 2022) [47]. A total of 10,935 genes in the background set were associated with a GO term in the database. Weighted gene co-expression network analysis We constructed weighted gene co-expression networks to identify clusters of genes that were co-expressed across our samples [48]. To prepare the gene expression data for this analysis, we applied a variance-stabilizing transformation to the filtered reads using the varianceStabilizingTransformation function from the DESeq2 package (v.1.36.0) [49] in R, which allowed us to normalize the read counts relative to library size and appropriately scale the data for clustering [50, 51]. After transforming and normalizing the read counts, we used the cpm function in the edgeR package (v.3.38.1) [43] to generate a gene matrix of scaled read counts (log2-cpm, counts per million mapped reads) from the transformed read counts [52]. As outlined in the Weighted Gene Co-expression Network Analysis (WGCNA) package documentation [50], we used hierarchical clustering to cluster the samples based on their gene expression profiles. We used the hclust function from the flashClust package (v.1.1.2) [51] to cluster the samples. We then created a weighted network adjacency matrix using the adjacency function from the WGCNA package (v.1.71) [50, 51]. The adjacency matrix was constructed by calculating pairwise co-expression similarities (Pearson correlation coefficients) and raising them to a power of β, a soft thresholding power [50]. We used the pickSoftThreshold function in the WGCNA package [50] to assist in selecting the lowest possible value of β that ensured our network fit the approximate scale free topology criterion while retaining the highest possible mean connectivity between the network genes. Based on the scale free topology model fit and mean connectivity of our network, and following the recommendations outlined by Zhang and Horvath [48], we selected β = 7. From our correlation network, we then generated a topological overlap dissimilarity matrix to identify modules of co-expressed genes. We calculated dissimilarity between the genes by converting the adjacencies into topological overlap similarities using the TOMsimilarity function in the WGCNA package [50] and then subtracting these topological overlap measures from 1. To identify modules of co-expressed genes, we used the hclust function from the flashClust package [51] for hierarchical clustering of the genes and then created a hierarchical clustering dendrogram. We used the cutreeDynamic function from the dynamicTreeCut package [50, 53] to extract the modules from the dendrogram. To summarize the gene expression variation in each module, the first principal component of each module in the expression matrix (the eigengene) was calculated using the moduleEigengenes function from the WGCNA package [50]. We used the function mergeCloseModules in the WGCNA package to merge eigengenes that were highly correlated. We included eigengenes with a correlation greater than 0.9 for merging. Finally, we identified modules that were significantly associated with the presence or absence of salinity by calculating correlation coefficients between the eigengenes and the habitat type. P-values of the correlation coefficients were calculated using the corPvalueStudent function from the WGCNA package, and we retained modules with P-values less than 0.01 for functional enrichment analyses. Similar to our differential gene expression analysis, we used GOrilla [47] for Gene Ontology enrichment analysis of genes contained in modules exhibiting significant correlations with habitat type. Comparisons of L. perugiae with other species We mined previously published datasets to identify gene expression patterns commonly associated with salinity tolerance in disparate taxa, including South American silversides (Odontesthes bonariensis and Odontesthes argentinensis) from Hughes et al. [20], three-spine sticklebacks (Gasterosteus aculeatus) from Gibbons et al. [19], Amur ide (Leuciscus waleckii) from Xu et al. [13], and L. perugiae (this study; Table 1). Each of these experiments generated paired-end RNA-seq raw reads from gill tissue in two ecotypes (one freshwater and the other saltwater) of the same lineage. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. Species included in the analysis, including environment (freshwater [FW] or saltwater [SW]), collection location, sample size (N), NCBI Sequence Read Archive (SRA) accession numbers, and study reference. https://doi.org/10.1371/journal.pone.0315014.t001 Raw RNA-seq reads from each transcriptomics project were downloaded in FASTQ format, and reference genomes or transcriptomes for each species were downloaded in FASTA format from Genbank (see Table 1 for accession numbers). Reads were trimmed and mapped to their respective reference genomes (Poecilia mexicana: GCF_001443325.1 [38]; Cyprinus carpio: GCF_000951615.1 [54]; Gasterosteus aculeatus: Broad S1 v. 93 [55]) or reference transcriptome (Menidia menidia: GEVY00000000.1 [56]) following the same methods described above for L. perugiae. We then quantified the number of reads mapped to each gene in the annotation file for each reference genome and created a read counts matrix for each species, which were used for subsequent expression analyses. Expression analyses were performed in R version 4.1.2. The 10,000 genes with the highest standard deviation between freshwater and saltwater samples were abstracted from each read counts matrix, and overall expression patterns were visualized with multi-dimensional scaling (MDS) plots. To make comparisons across species, we used OrthoFinder v2.2.6 to identify orthologous genes among the reference genomes [57, 58]. For OrthoFinder, we used the ‘dendroblast’ option for gene tree inference, ‘blast’ for the sequence search program, ‘mafft’ for the multiple sequence alignment, and ‘fasttree’ for the tree inference method. A total of 18,419 orthogroups were identified. Of the total orthogroups identified, 13,899 orthogroups had at least one copy in each species. Of those, 1,735 were 1:1 orthologs. To calculate counts per orthogroup, we used the gene counts matrix of each species to sum up the counts across all loci contained in an orthogroup. Based on this orthogroup counts matrix, we retained only the orthogroups that were present in all species and expressed in all individuals (cpm > 0 per individual), resulting in 12,743 retained orthogroups. To evaluate expression differences for each orthogroup, we made pairwise comparisons between ecotypes following the same methods described above for the L. perugiae comparisons. Briefly, we normalized reads, created and compared generalized linear models of the normalized read counts, generated a design matrix, estimated tagwise dispersion, and conducted GLM likelihood-ratio tests to test whether differences in expression were statistically different between the freshwater and saltwater population for each orthgroup. To identify orthogroups exhibiting convergent expression patterns across lineages, we intersected the significantly upregulated and downregulated orthogroups from all lineage-specific comparisons, identifying orthogroups that were differentially expressed in the same direction in pairwise, three-, and four-way comparisons among the lineages. After identifying the set of orthogroups with differential expression across three or more lineages, we used a GO enrichment analysis as described above to explore the putative biological functions of these candidate gene sets. Results Comparative analysis of freshwater and hypersaline L. perugiae We used RNA-seq to characterize the transcriptomes of L. perugiae from a freshwater (n = 6) and a hypersaline populations (Table 1). 73,590,325 total raw reads were obtained across all individuals: 34,998,157 from freshwater L. perugiae (n = 6) and 38,592,168 from saltwater L. perugiae (n = 6) before trimming (Table A in S1 Appendix). After trimming, 95.5% of reads from the freshwater individuals mapped to the Poecilia mexicana reference genome, and 95.2% mapped for the saltwater individuals (Table A in S1 Appendix). We identified 4,895 differentially expressed genes between saltwater and freshwater ecotypes of L. perugiae, 2,437 of which were upregulated and 2,458 of which were downregulated in the saltwater ecotype (Fig 1A). The genes upregulated in the saltwater population were largely associated with ion transport, maintaining chemical homeostasis, and cell signaling (FDR < 0.05) (Table 2 and S1 Fig). Processes relevant to chemical homeostasis included several solute carrier genes, such as SLC9A3, SLC8B1, SLC12A8, and SLC30A9. Na+/K+-ATPase and other ATPase genes, such as ATP1B1 and ATP6V1A, were also upregulated among the genes involved in chemical homeostasis and signal transduction. Genes downregulated in the saltwater population corresponded to GO process terms such as mitotic cell cycle process, protein folding, chromosome segregation, rRNA processing, and mitochondrial translational elongation (FDR < 0.05; Table 2 and S2 Fig). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. Differentially expressed genes and gene expression profiles of hypersaline and freshwater Limia perugiae. A. Volcano plot depicting differentially expressed genes between hypersaline and freshwater Limia perugiae. Genes that were significantly differentially expressed between hypersaline and freshwater populations (FDR < 0.05) are indicated by the blue and red points—blue points represent genes downregulated in the hypersaline populations, while red points represent upregulated genes. B. Multi-dimensional scaling (MDS) plot of hypersaline and freshwater L. perugiae gene expression profiles. MDS axis 1 separated samples by freshwater vs. hypersaline environments. https://doi.org/10.1371/journal.pone.0315014.g001 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. GO process terms with significant enrichment in genes upregulated and downregulated in the hypersaline ecotype of Limia perugiae (FDR < 0.05). https://doi.org/10.1371/journal.pone.0315014.t002 Weighted gene co-expression network analysis (WGCNA) revealed five modules of co-expressed genes that were significantly correlated with salinity (Fig 2). The turquoise module was positively correlated with salinity (P-value < 0.01), and the black, red, royalblue, and blue modules were all negatively correlated with salinity (Fig 2). Out of the 18,659 genes included in our analysis, the positively correlated module (turquoise module) contained 4,056 genes. The negatively correlated modules contained 2,166 genes (black module), 837 genes (red module), 173 genes (royalblue module), and 3,300 genes (blue module). The correlation coefficients and their associated P-values between each gene and the environmental condition (salinity), and between each gene and each module, are included in Table B in S1 Appendix. Each gene’s module assignment can also be found in Table B in S1 Appendix. From the functional enrichment analysis, we found that the turquoise, royalblue, and blue modules were significantly enriched for biological processes, and these modules of co-expressed genes largely corroborated the differential expression results. Like the biological processes that were enriched among the up-regulated genes in the saltwater population, the module that was positively correlated with salinity (turquoise) was functionally enriched for GO terms involved in ion transport and cell signaling (Table C in S1 Appendix). There were, however, GO terms associated with the turquoise module that were not identified from the differential expression analysis, including regulation of autophagy and lipid transport (Table C in S1 Appendix). The modules that were negatively correlated with salinity also reflected biological processes that were associated with down-regulated genes in the saltwater L. perugiae, including regulation of the cell cycle and protein folding (Tables D and E in S1 Appendix). The royalblue module only had one significantly enriched GO term, chemokine-mediated signaling pathway, which was a GO term also associated with down-regulation of genes in the saltwater environment (Table D in S1 Appendix). In contrast, the blue module contained several significant GO terms, including mitotic cell cycle process, protein folding, rRNA metabolic process, and chromosome organization (Table E in S1 Appendix). There were many GO terms associated with the blue module that were not represented in the differential expression analysis, including sarcomere organization, macromolecule biosynthetic process, regulation of cellular response to heat, and nucleic acid metabolic process (Table E in S1 Appendix). Several of the GO terms unique to the blue module were related to cellular metabolism and biosynthesis pathways (Table E in S1 Appendix). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Weighted gene co-expression network analysis. A. Average linkage clustering tree based on topological overlap distances in gene expression patterns of L. perugiae from freshwater and saltwater habitats. Branches of the dendrogram correspond to modules, as shown in the color bars below. DC is an abbreviation for Dynamic Tree Cut; MD for Merged Dynamic. B. Correlation between module eigenvalues and habitat type (freshwater vs. saltwater). Each row corresponds to a module of coexpressed genes, and values are Pearson correlation coefficients (left column) and P-values (right column in parentheses). Color coloration scales with the correlation coefficient according to the scale bar to the right. https://doi.org/10.1371/journal.pone.0315014.g002 Comparisons of Limia with phylogenetically disparate teleosts We compared transcriptomic differences between L. perugiae populations to previously published transcriptome data from teleosts with populations from freshwater and saline habitats (Table 1). Mapping statistics for all four population pairs can be found in Table A in S1 Appendix. As expected, MDS plots indicated that orthogroup expression variation was primarily driven by differences among taxonomic groups, with much smaller differences between ecotypes within species (Fig 3A). We then compared orthogroup expression profiles of all the freshwater and saltwater lineages based on mean expression values and found that the variation in orthogroup expression largely reflects phylogenetic divergence among lineages (Fig 3B). Closely related lineages exhibited more similar expression profiles, irrespective of environmental conditions (saltwater vs. freshwater; Fig 3B). Across all groups, 122 differentially expressed orthogroups were shared across at least three lineages, but only 10 shared orthogroups were differentially expressed across all freshwater and saltwater population pairs (Fig 3C). Of those 10 orthogroups, 9 had annotations in the SWISS-PROT database (Table 3 and Fig 4), including a Na+/H+-exchanger (SLC9A3) involved in osmoregulation. The directionality of differential expression varied among lineages, and none of the 10 shared differentially expressed orthogroups were up-regulated or down-regulated among all four population pairs (Fig 4). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. Gene expression profiles and shared differentially expressed genes across lineages. A. Multi-dimensional scaling plot (MDS) of the general expression patterns of all the populations included in our analysis. B. Similarity of gene expression profiles of saltwater (SW) and freshwater (FW) populations across different lineages. The majority of variation in gene expression reflects phylogenetic divergence among lineages. C. Shared differentially expressed genes across lineages. The large, central number in each section represents the total number of shared differentially expressed genes among the lineages in that intersection. The top number in each section represents the number of shared up-regulated genes, and the bottom number is the number of shared down-regulated genes in that intersection. Only 10 genes were consistently differentially expressed between all of the SW and FW populations. https://doi.org/10.1371/journal.pone.0315014.g003 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Examples of expression variation in shared differentially expressed genes. Nine of the ten shared differentially expressed genes have annotations, and those genes are included here. Magnitude and direction of differential expression is lineage specific. https://doi.org/10.1371/journal.pone.0315014.g004 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 3. Shared differentially expressed genes among all four population pairs and their associated proteins. https://doi.org/10.1371/journal.pone.0315014.t003 To investigate the biological processes reflected in shared differentially expressed orthogroups, we analyzed the differentially expressed orthogroups that were shared among three or more lineages (Fig 3C). GO analysis indicated enrichment in biological processes associated mostly with transmembrane transport (particularly ion transport) and some associated with immune function, which were all significant based on P-value (P-value < 0.001) but not after FDR correction (Table 4 and S3 Fig). Some of the GO processes specific to transmembrane transport included anion transport, inorganic cation import across plasma membrane, potassium ion import, urea transmembrane transport, and pyruvate transmembrane transport (P-value < 0.001; FDR > 0.05). The shared differentially expressed immune genes were reflective of negative regulation of T-helper cell differentiation, negative regulation of leukocyte differentiation, negative regulation of CD4-positive, alpha-beta T cell differentiation, and regulation of adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily (P-value < 0.001; FDR > 0.05). Download: PPT PowerPoint slide PNG larger image TIFF original image Table 4. GO process terms that had significant enrichment in orthogenes with significant differential expression in at least three lineages (P-value < 0.001), but not after FDR correction (FDR > 0.05). https://doi.org/10.1371/journal.pone.0315014.t004 Comparative analysis of freshwater and hypersaline L. perugiae We used RNA-seq to characterize the transcriptomes of L. perugiae from a freshwater (n = 6) and a hypersaline populations (Table 1). 73,590,325 total raw reads were obtained across all individuals: 34,998,157 from freshwater L. perugiae (n = 6) and 38,592,168 from saltwater L. perugiae (n = 6) before trimming (Table A in S1 Appendix). After trimming, 95.5% of reads from the freshwater individuals mapped to the Poecilia mexicana reference genome, and 95.2% mapped for the saltwater individuals (Table A in S1 Appendix). We identified 4,895 differentially expressed genes between saltwater and freshwater ecotypes of L. perugiae, 2,437 of which were upregulated and 2,458 of which were downregulated in the saltwater ecotype (Fig 1A). The genes upregulated in the saltwater population were largely associated with ion transport, maintaining chemical homeostasis, and cell signaling (FDR < 0.05) (Table 2 and S1 Fig). Processes relevant to chemical homeostasis included several solute carrier genes, such as SLC9A3, SLC8B1, SLC12A8, and SLC30A9. Na+/K+-ATPase and other ATPase genes, such as ATP1B1 and ATP6V1A, were also upregulated among the genes involved in chemical homeostasis and signal transduction. Genes downregulated in the saltwater population corresponded to GO process terms such as mitotic cell cycle process, protein folding, chromosome segregation, rRNA processing, and mitochondrial translational elongation (FDR < 0.05; Table 2 and S2 Fig). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. Differentially expressed genes and gene expression profiles of hypersaline and freshwater Limia perugiae. A. Volcano plot depicting differentially expressed genes between hypersaline and freshwater Limia perugiae. Genes that were significantly differentially expressed between hypersaline and freshwater populations (FDR < 0.05) are indicated by the blue and red points—blue points represent genes downregulated in the hypersaline populations, while red points represent upregulated genes. B. Multi-dimensional scaling (MDS) plot of hypersaline and freshwater L. perugiae gene expression profiles. MDS axis 1 separated samples by freshwater vs. hypersaline environments. https://doi.org/10.1371/journal.pone.0315014.g001 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. GO process terms with significant enrichment in genes upregulated and downregulated in the hypersaline ecotype of Limia perugiae (FDR < 0.05). https://doi.org/10.1371/journal.pone.0315014.t002 Weighted gene co-expression network analysis (WGCNA) revealed five modules of co-expressed genes that were significantly correlated with salinity (Fig 2). The turquoise module was positively correlated with salinity (P-value < 0.01), and the black, red, royalblue, and blue modules were all negatively correlated with salinity (Fig 2). Out of the 18,659 genes included in our analysis, the positively correlated module (turquoise module) contained 4,056 genes. The negatively correlated modules contained 2,166 genes (black module), 837 genes (red module), 173 genes (royalblue module), and 3,300 genes (blue module). The correlation coefficients and their associated P-values between each gene and the environmental condition (salinity), and between each gene and each module, are included in Table B in S1 Appendix. Each gene’s module assignment can also be found in Table B in S1 Appendix. From the functional enrichment analysis, we found that the turquoise, royalblue, and blue modules were significantly enriched for biological processes, and these modules of co-expressed genes largely corroborated the differential expression results. Like the biological processes that were enriched among the up-regulated genes in the saltwater population, the module that was positively correlated with salinity (turquoise) was functionally enriched for GO terms involved in ion transport and cell signaling (Table C in S1 Appendix). There were, however, GO terms associated with the turquoise module that were not identified from the differential expression analysis, including regulation of autophagy and lipid transport (Table C in S1 Appendix). The modules that were negatively correlated with salinity also reflected biological processes that were associated with down-regulated genes in the saltwater L. perugiae, including regulation of the cell cycle and protein folding (Tables D and E in S1 Appendix). The royalblue module only had one significantly enriched GO term, chemokine-mediated signaling pathway, which was a GO term also associated with down-regulation of genes in the saltwater environment (Table D in S1 Appendix). In contrast, the blue module contained several significant GO terms, including mitotic cell cycle process, protein folding, rRNA metabolic process, and chromosome organization (Table E in S1 Appendix). There were many GO terms associated with the blue module that were not represented in the differential expression analysis, including sarcomere organization, macromolecule biosynthetic process, regulation of cellular response to heat, and nucleic acid metabolic process (Table E in S1 Appendix). Several of the GO terms unique to the blue module were related to cellular metabolism and biosynthesis pathways (Table E in S1 Appendix). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Weighted gene co-expression network analysis. A. Average linkage clustering tree based on topological overlap distances in gene expression patterns of L. perugiae from freshwater and saltwater habitats. Branches of the dendrogram correspond to modules, as shown in the color bars below. DC is an abbreviation for Dynamic Tree Cut; MD for Merged Dynamic. B. Correlation between module eigenvalues and habitat type (freshwater vs. saltwater). Each row corresponds to a module of coexpressed genes, and values are Pearson correlation coefficients (left column) and P-values (right column in parentheses). Color coloration scales with the correlation coefficient according to the scale bar to the right. https://doi.org/10.1371/journal.pone.0315014.g002 Comparisons of Limia with phylogenetically disparate teleosts We compared transcriptomic differences between L. perugiae populations to previously published transcriptome data from teleosts with populations from freshwater and saline habitats (Table 1). Mapping statistics for all four population pairs can be found in Table A in S1 Appendix. As expected, MDS plots indicated that orthogroup expression variation was primarily driven by differences among taxonomic groups, with much smaller differences between ecotypes within species (Fig 3A). We then compared orthogroup expression profiles of all the freshwater and saltwater lineages based on mean expression values and found that the variation in orthogroup expression largely reflects phylogenetic divergence among lineages (Fig 3B). Closely related lineages exhibited more similar expression profiles, irrespective of environmental conditions (saltwater vs. freshwater; Fig 3B). Across all groups, 122 differentially expressed orthogroups were shared across at least three lineages, but only 10 shared orthogroups were differentially expressed across all freshwater and saltwater population pairs (Fig 3C). Of those 10 orthogroups, 9 had annotations in the SWISS-PROT database (Table 3 and Fig 4), including a Na+/H+-exchanger (SLC9A3) involved in osmoregulation. The directionality of differential expression varied among lineages, and none of the 10 shared differentially expressed orthogroups were up-regulated or down-regulated among all four population pairs (Fig 4). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. Gene expression profiles and shared differentially expressed genes across lineages. A. Multi-dimensional scaling plot (MDS) of the general expression patterns of all the populations included in our analysis. B. Similarity of gene expression profiles of saltwater (SW) and freshwater (FW) populations across different lineages. The majority of variation in gene expression reflects phylogenetic divergence among lineages. C. Shared differentially expressed genes across lineages. The large, central number in each section represents the total number of shared differentially expressed genes among the lineages in that intersection. The top number in each section represents the number of shared up-regulated genes, and the bottom number is the number of shared down-regulated genes in that intersection. Only 10 genes were consistently differentially expressed between all of the SW and FW populations. https://doi.org/10.1371/journal.pone.0315014.g003 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Examples of expression variation in shared differentially expressed genes. Nine of the ten shared differentially expressed genes have annotations, and those genes are included here. Magnitude and direction of differential expression is lineage specific. https://doi.org/10.1371/journal.pone.0315014.g004 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 3. Shared differentially expressed genes among all four population pairs and their associated proteins. https://doi.org/10.1371/journal.pone.0315014.t003 To investigate the biological processes reflected in shared differentially expressed orthogroups, we analyzed the differentially expressed orthogroups that were shared among three or more lineages (Fig 3C). GO analysis indicated enrichment in biological processes associated mostly with transmembrane transport (particularly ion transport) and some associated with immune function, which were all significant based on P-value (P-value < 0.001) but not after FDR correction (Table 4 and S3 Fig). Some of the GO processes specific to transmembrane transport included anion transport, inorganic cation import across plasma membrane, potassium ion import, urea transmembrane transport, and pyruvate transmembrane transport (P-value < 0.001; FDR > 0.05). The shared differentially expressed immune genes were reflective of negative regulation of T-helper cell differentiation, negative regulation of leukocyte differentiation, negative regulation of CD4-positive, alpha-beta T cell differentiation, and regulation of adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily (P-value < 0.001; FDR > 0.05). Download: PPT PowerPoint slide PNG larger image TIFF original image Table 4. GO process terms that had significant enrichment in orthogenes with significant differential expression in at least three lineages (P-value < 0.001), but not after FDR correction (FDR > 0.05). https://doi.org/10.1371/journal.pone.0315014.t004 Discussion In this study, we used an RNA-seq approach to better understand mechanisms of salinity tolerance in the livebearing fish Limia perugiae (Poeciliidae), as well as to investigate whether mechanisms of osmoregulation are shared across distantly related species of teleosts that have all made recent transitions between freshwater and saltwater environments. We compared patterns of gene expression between a freshwater and a hypersaline population of L. perugiae and found that the genes upregulated in the saltwater population were largely related to ion transport and maintaining chemical homeostasis, while downregulated genes were associated with processes involved in the cell cycle regulation and protein folding. These results provide insight into how L. perugiae has colonized a novel environment and maintains homeostasis under extreme salinity stress. Comparisons of our Limia results with pre-existing gene expression data collected from freshwater and saltwater ecotypes of South American silversides (Odontesthes spp.) [20], three-spine stickleback (Gasterosteus aculeatus) [19], and Amur ide (Leuciscus waleckii) [13] indicated that there were few shared differentially expressed genes among all four ecotype pairs. Variation in gene expression was largely shaped by phylogeny rather than environment, and the shared differentially expressed genes among all four pairs showed strong variation in the direction and magnitude of differential expression across lineages. Overall, these results suggest that disparate lineages utilize different mechanisms for overcoming salinity challenges—at least at the level of gene expression. We found that various patterns of gene expression can emerge from crossing saltwater-freshwater boundaries, providing evidence of diverse responses of teleost lineages to a similar environmental challenge. A major question remaining is to what degree variation in gene expression in L. perugiae and the other lineages were shaped by phenotypic plasticity and by genetic differences in gene regulation. A previous study in stickleback found evidence for heritable gene expression differences between freshwater and saltwater ecotypes, evidence for shared plastic responses between ecotypes, but only little evidence for ecotype-specific plasticity [19]. Similar studies that combine field studies with laboratory experimentation are highly warranted for other study systems. Responses to variation in salinity in L. perugiae Regulation of transmembrane transport and gill epithelial permeability. Overcoming the physiological challenges associated with transitioning between saltwater and freshwater environments requires modification of ion transport across the gill epithelia [16]. In our analysis of differentially expressed genes between a freshwater and hypersaline population of L. perugiae, we found evidence for differential regulation of ion transport. Among the genes that were upregulated in the hypersaline population, GO enrichment analysis revealed that several terms were associated with anion transport, sodium ion transport, metal ion transport, and maintaining chemical homeostasis. Genes corresponding to solute carrier families (e.g., SLC9A3, SLC8B1, SLC12A8, SLC30A9, SLC7A1) and ATPases (e.g., ATP1B1, ATP6V1A, ATP13A3) were among the upregulated genes that are important in maintaining ion and chemical homeostasis [59]. Increasing the rate of inorganic ion, amino acid, and nucleotide transport via upregulation of solute carriers allows aquatic organisms to maintain osmotic balance in saline environments, potentially facilitating acclimation or adaptation to hypersalinity [59]. Specifically, ATPases—such as Na+/K+-ATPase—have been well-studied for their role in salinity tolerance during both acclimation and adaptation [23, 31, 60–64]. Consistent with our study, previous Western blot analyses of freshwater and hypersaline L. perugiae populations have revealed an increase in Na+/K+-ATPase expression in hypersaline L. perugiae when compared to their freshwater conspecifics, which is essential for excreting Na+ and Cl- out of the body to maintain homeostasis [23]. At very high salinity levels, it is especially difficult to pump Na+ against the chemical gradient and out of the gill [65]. In a hypersaline environment, it therefore may be beneficial to maintain a lower epithelial salt permeability to avoid a back-flux of salts across the tight junctions of the gill epithelia [65, 66]. We found up-regulation of genes involved in cell-cell junction organization, including claudins involved in the formation of tight junctions (CLDN3, CLDN4, CLDN5, CLDN8, and CLDN10). Expression of gill claudin genes has been associated with salinity acclimation in fish [64, 66], and up-regulation of claudin-8 (CLDN8) has been shown to reduce the paracellular barrier permeability to Na+ [67], suggesting a role for regulation of gill epithelial permeability via tight junctions during osmoregulatory challenges. Our findings support the hypothesis that modifications to transmembrane transport and gill epithelial permeability jointly contribute to the salt exclusion necessary for survival in a hypersaline environment. Cell cycle regulation and cellular metabolic and biosynthesis pathways. Populations living in hypersaline environments must cope with the adverse effects of high salinity levels, often resulting in strategies for damage repair and energy reallocation [23, 62]. Functional analysis of the downregulated genes in the hypersaline population of L. perugiae showed signatures of downregulation of cell cycle and protein folding processes. Downregulation of genes involved in the cell cycle is expected to occur under high salinity stress to stop the replication of damaged cells and allow time for DNA repair [68, 69]. Similarly, locally adapted killifish living in freshwater and brackish environments exhibit evidence for divergence in the expression of genes underlying control of the cell cycle, supporting a central role in cell cycle regulation in adaptation to salinity stress [68, 70]. However, our findings may also suggest a general re-allocation of energy resources [68]. In hypersaline L. perugiae, the downregulation of genes involved in the cell cycle and protein folding could be reflective of an energetic trade-off between increased demand for processes involved in osmoregulation and the energetic cost of growth under stressful conditions [23]. In support of such energetic trade-offs, hypersaline L. perugiae have morphological differences from their freshwater counterparts, including a significantly smaller body size and reduced secondary sex features in males [23, 25]. Furthermore, energetic trade-offs associated with increased energy investment in osmoregulation in a hypersaline environment may also result in reduced investment in cellular metabolic and biosynthetic processes. In addition to genes involved in cell cycle processes, the weighted gene co-expression network analysis revealed that several cellular metabolic and biosynthetic processes were negatively associated with salinity (associated with genes in the blue module), including cellular aromatic compound metabolic process, organic cyclic compound metabolic process, nucleic acid metabolic process, cellular metabolic process, small molecule biosynthetic process, cellular macromolecule biosynthetic process, and carbohydrate biosynthetic process. Cellular metabolism and biosynthesis of compounds are energetically costly pathways for organisms [71], so these functions may receive less energy investment due to the high amounts of ATP required for ion transport and other osmoregulatory functions [68]. Future research will be important for identifying how osmoregulation in hypersaline environments may carry costs that influence cellular processes and how these costs impact organismal performance under different environmental conditions. Regulation of cell signaling. To respond to osmoregulatory challenges, fishes must perceive their environmental salinity and maintain intracellular signals that modulate ion transport and other processes [72, 73]. We found differentially expressed genes associated with a variety of cell signaling processes. Several upregulated genes in the hypersaline population were enriched in GO terms associated with cell signaling, including signal transduction, regulation of signaling, G-protein-coupled receptor signaling pathway, and second-messenger-mediated signaling. G-protein-coupled receptor signaling is among the pathways involved in allowing aquatic organisms to sense their environmental salinity [74, 75]. Some of these signaling pathway components may also play a role in regulating ion transport, such as mitogen-activated protein kinases (MAPKs) and other serine/threonine protein kinases [76–79], which were among the upregulated genes. MAPK signaling pathways are also implicated in regulation of the cell cycle, some of which trigger cellular growth arrest and DNA damage repair in response to osmotic stress [79, 80]. The upregulation of signaling pathways, especially MAPK genes, may consequently play a role in the downregulation of cell-cycle genes that we found in the hypersaline L. perugiae population. Comparisons among replicated lineages Lineage-specific responses to variation in salinity. Convergence in gene expression among independently evolved populations can occur in response to shared environmental stressors [81, 82]. Among fishes, cases of convergence in gene expression have been identified in response to selective pressures such as pollution [81], hydrogen sulfide [82], and the absence of light in caves [83]. In contrast, our comparative transcriptomics analysis indicated there is little evidence for convergence in gene expression patterns among fish lineages in response to variation in salinity. Even among the shared differentially expressed genes, there is considerable variation in the magnitude and direction of expression differences between ecotypes across lineages. This finding mirrors genomic analyses that looked for convergent signatures associated with adaptation to salinity variation, which found osmoregulation genes to be common targets of selection but also variation in the exact genes that were involved across different lineages [84]. Our finding of lineage-specific responses to variation in salinity, and therefore low levels of convergence in gene expression, could be explained by several non-mutually exclusive hypotheses [85]. First, changes in both protein-coding DNA sequences and gene expression may occur during adaptation to a novel environment, or they may occur independently [55, 86, 87]. Under certain conditions, selection may favor changes in protein structure and function via modifications at the sequence level, regardless of whether or not gene expression is affected [86]. Alternatively, changes in gene expression may be favored without selection acting at the sequence level [86]. Variation in how selection acts on gene expression could contribute to the lack of convergence among the lineages in our study, and there also may be stronger signals of convergence at levels other than gene expression. Secondly, differences in genetic architecture among lineages may cause different responses to selection, leading to diverse evolutionary outcomes [85]. Convergence is therefore less likely to occur with increasing genetic divergence between lineages [88], so it is not necessarily surprising that we did not find convergence among distantly related fishes. Even if the genetic architectures are similar among lineages, idiosyncratic responses to selection can also arise as a consequence of functional redundancy [89]. Specifically, modification of different genes and pathways may have equivalent functional and fitness consequences [78, 84]. Functional redundancy is particularly common in complex traits like salinity tolerance that involve many genes and physiological pathways [89]. While salinity tolerance is a shared outcome among the lineages in our analysis, the molecular mechanisms underlying osmoregulation may be unique to each lineage. Third, the directionality of habitat transition may impact gene expression responses, especially if there are genetic adaptations to the habitat of origin. For example, a high-activity version of an ion transporter may be downregulated during transitions to freshwater, while a low-activity version of the same enzyme may be upregulated in the opposite direction. Among the four lineages we included in our study, two of them transitioned from saltwater to freshwater environments (Gasterosteus and Odontesthes), and two transitioned from freshwater to saltwater environments (Limia and Leuciscus). If the direction of transition elicits similar responses, then we would expect populations that transitioned in the same direction to share more differentially expressed genes than those that did not transition in the same direction. However, we did not find this to be the case, as the pairs transitioning in the same direction share fewer differentially expressed genes than those that made opposite transitions (Fig 3C). Finally, covariation with other sources of selection—both abiotic and biotic—may cause idiosyncratic gene expression patterns across lineages. Beyond the challenges directly imposed by variation in salinity, such habitat transitions are often accompanied by other environmental challenges, such as variation in temperature, exposure to novel parasites, and restructuring of host-associated microbial communities [2, 20, 90, 91]. For example, hypersaline Amur ide must cope with high alkalinity stress in addition to salinity, and hypersaline L. perugiae have a warmer environment than freshwater L. perugiae [13, 23]. Additionally, there are differences in the gill microbial communities of saltwater and freshwater populations of South American silversides [20]. Such correlated environmental factors can contribute variation in gene expression responses to salinity transitions we observed. Shared responses involving ion transport and immune function. Although there is little evidence for convergence among all four lineage pairs included in our analysis, there was some functional overlap in the differentially expressed genes shared among three or more lineages. Most of the shared differentially expressed genes were associated with ion transport and immune system processes. Genes implicated in ion transport and immune system processes were also identified among the ten differentially expressed genes that were shared among all four lineages. For example, SLC9A3, a solute carrier gene involved in osmoregulation, was among the ten shared genes, suggesting its role in salinity transitions among fishes [59]. Another shared differentially expressed gene, CEACAM1 (Carcinoembryonic antigen-related cell adhesion molecule 1), has been implicated in immune processes in other systems and is thought to be involved in a variety of pathways, but its function is not well known [92, 93]. As previously discussed, regulation of ion transport is expected to be crucial when crossing a saltwater-freshwater boundary [16], and osmoregulation genes also exhibit evidence of convergent evolution during salinity transitions [84].Variation in expression of immune genes, particularly those related to inflammation and adaptive immunity, has also been documented in fish during salinity acclimation [94, 95]. Under a variety of selection pressures, locally adapted populations also frequently show divergence in immune genes due to other factors such as life history differences, different parasite exposure, and shifts in the microbiome [20, 95–98]. Immune loci are consequently evolutionary hotspots in diversification [99–101]. It needs to be tested whether changes in the expression of immune-related genes are directly linked to variation in salinity or whether these genes generally respond to changes in correlated biotic sources of selection. Overall, our analysis of gene expression patterns between locally adapted freshwater and hypersaline populations of L. perugiae provide insight into how this livebearing fish maintains homeostasis in a hypersaline environment. In addition, comparisons between four population pairs of freshwater and saltwater ecotypes in disparate teleost lineages showed little evidence for convergence, as there were only ten differentially expressed genes that were shared among them all. Despite this, we found that the differentially expressed genes shared in three or more of the lineages reflected biological processes related to ion transport and immune functioning. Our findings provide insight into shared and unique gene expression responses to salinity variation, broadly informing our understanding of salinity tolerance in aquatic organisms. Furthermore, comparisons in gene expression across species that have made habitat transitions are important for understanding mechanisms of adaptation to novel environments. Our results shed light on the repeatability of transcriptomic responses to salinity variation among fishes, which remains a relatively underexplored area of research despite its relevance for aquatic ecology and evolutionary biology. Future research comparing gene expression in fishes from freshwater and saltwater environments in both laboratory and field settings will be important for identifying the degree in which shared or unique gene expression responses are due to plasticity or adaptation, and additional studies that take a comparative transcriptomic approach with more lineages—representing transitions from freshwater to saltwater and vice versa as well as spanning a range of transition time and divergence time between populations—will provide further insight into shared mechanisms of osmoregulation among fishes. Responses to variation in salinity in L. perugiae Regulation of transmembrane transport and gill epithelial permeability. Overcoming the physiological challenges associated with transitioning between saltwater and freshwater environments requires modification of ion transport across the gill epithelia [16]. In our analysis of differentially expressed genes between a freshwater and hypersaline population of L. perugiae, we found evidence for differential regulation of ion transport. Among the genes that were upregulated in the hypersaline population, GO enrichment analysis revealed that several terms were associated with anion transport, sodium ion transport, metal ion transport, and maintaining chemical homeostasis. Genes corresponding to solute carrier families (e.g., SLC9A3, SLC8B1, SLC12A8, SLC30A9, SLC7A1) and ATPases (e.g., ATP1B1, ATP6V1A, ATP13A3) were among the upregulated genes that are important in maintaining ion and chemical homeostasis [59]. Increasing the rate of inorganic ion, amino acid, and nucleotide transport via upregulation of solute carriers allows aquatic organisms to maintain osmotic balance in saline environments, potentially facilitating acclimation or adaptation to hypersalinity [59]. Specifically, ATPases—such as Na+/K+-ATPase—have been well-studied for their role in salinity tolerance during both acclimation and adaptation [23, 31, 60–64]. Consistent with our study, previous Western blot analyses of freshwater and hypersaline L. perugiae populations have revealed an increase in Na+/K+-ATPase expression in hypersaline L. perugiae when compared to their freshwater conspecifics, which is essential for excreting Na+ and Cl- out of the body to maintain homeostasis [23]. At very high salinity levels, it is especially difficult to pump Na+ against the chemical gradient and out of the gill [65]. In a hypersaline environment, it therefore may be beneficial to maintain a lower epithelial salt permeability to avoid a back-flux of salts across the tight junctions of the gill epithelia [65, 66]. We found up-regulation of genes involved in cell-cell junction organization, including claudins involved in the formation of tight junctions (CLDN3, CLDN4, CLDN5, CLDN8, and CLDN10). Expression of gill claudin genes has been associated with salinity acclimation in fish [64, 66], and up-regulation of claudin-8 (CLDN8) has been shown to reduce the paracellular barrier permeability to Na+ [67], suggesting a role for regulation of gill epithelial permeability via tight junctions during osmoregulatory challenges. Our findings support the hypothesis that modifications to transmembrane transport and gill epithelial permeability jointly contribute to the salt exclusion necessary for survival in a hypersaline environment. Cell cycle regulation and cellular metabolic and biosynthesis pathways. Populations living in hypersaline environments must cope with the adverse effects of high salinity levels, often resulting in strategies for damage repair and energy reallocation [23, 62]. Functional analysis of the downregulated genes in the hypersaline population of L. perugiae showed signatures of downregulation of cell cycle and protein folding processes. Downregulation of genes involved in the cell cycle is expected to occur under high salinity stress to stop the replication of damaged cells and allow time for DNA repair [68, 69]. Similarly, locally adapted killifish living in freshwater and brackish environments exhibit evidence for divergence in the expression of genes underlying control of the cell cycle, supporting a central role in cell cycle regulation in adaptation to salinity stress [68, 70]. However, our findings may also suggest a general re-allocation of energy resources [68]. In hypersaline L. perugiae, the downregulation of genes involved in the cell cycle and protein folding could be reflective of an energetic trade-off between increased demand for processes involved in osmoregulation and the energetic cost of growth under stressful conditions [23]. In support of such energetic trade-offs, hypersaline L. perugiae have morphological differences from their freshwater counterparts, including a significantly smaller body size and reduced secondary sex features in males [23, 25]. Furthermore, energetic trade-offs associated with increased energy investment in osmoregulation in a hypersaline environment may also result in reduced investment in cellular metabolic and biosynthetic processes. In addition to genes involved in cell cycle processes, the weighted gene co-expression network analysis revealed that several cellular metabolic and biosynthetic processes were negatively associated with salinity (associated with genes in the blue module), including cellular aromatic compound metabolic process, organic cyclic compound metabolic process, nucleic acid metabolic process, cellular metabolic process, small molecule biosynthetic process, cellular macromolecule biosynthetic process, and carbohydrate biosynthetic process. Cellular metabolism and biosynthesis of compounds are energetically costly pathways for organisms [71], so these functions may receive less energy investment due to the high amounts of ATP required for ion transport and other osmoregulatory functions [68]. Future research will be important for identifying how osmoregulation in hypersaline environments may carry costs that influence cellular processes and how these costs impact organismal performance under different environmental conditions. Regulation of cell signaling. To respond to osmoregulatory challenges, fishes must perceive their environmental salinity and maintain intracellular signals that modulate ion transport and other processes [72, 73]. We found differentially expressed genes associated with a variety of cell signaling processes. Several upregulated genes in the hypersaline population were enriched in GO terms associated with cell signaling, including signal transduction, regulation of signaling, G-protein-coupled receptor signaling pathway, and second-messenger-mediated signaling. G-protein-coupled receptor signaling is among the pathways involved in allowing aquatic organisms to sense their environmental salinity [74, 75]. Some of these signaling pathway components may also play a role in regulating ion transport, such as mitogen-activated protein kinases (MAPKs) and other serine/threonine protein kinases [76–79], which were among the upregulated genes. MAPK signaling pathways are also implicated in regulation of the cell cycle, some of which trigger cellular growth arrest and DNA damage repair in response to osmotic stress [79, 80]. The upregulation of signaling pathways, especially MAPK genes, may consequently play a role in the downregulation of cell-cycle genes that we found in the hypersaline L. perugiae population. Regulation of transmembrane transport and gill epithelial permeability. Overcoming the physiological challenges associated with transitioning between saltwater and freshwater environments requires modification of ion transport across the gill epithelia [16]. In our analysis of differentially expressed genes between a freshwater and hypersaline population of L. perugiae, we found evidence for differential regulation of ion transport. Among the genes that were upregulated in the hypersaline population, GO enrichment analysis revealed that several terms were associated with anion transport, sodium ion transport, metal ion transport, and maintaining chemical homeostasis. Genes corresponding to solute carrier families (e.g., SLC9A3, SLC8B1, SLC12A8, SLC30A9, SLC7A1) and ATPases (e.g., ATP1B1, ATP6V1A, ATP13A3) were among the upregulated genes that are important in maintaining ion and chemical homeostasis [59]. Increasing the rate of inorganic ion, amino acid, and nucleotide transport via upregulation of solute carriers allows aquatic organisms to maintain osmotic balance in saline environments, potentially facilitating acclimation or adaptation to hypersalinity [59]. Specifically, ATPases—such as Na+/K+-ATPase—have been well-studied for their role in salinity tolerance during both acclimation and adaptation [23, 31, 60–64]. Consistent with our study, previous Western blot analyses of freshwater and hypersaline L. perugiae populations have revealed an increase in Na+/K+-ATPase expression in hypersaline L. perugiae when compared to their freshwater conspecifics, which is essential for excreting Na+ and Cl- out of the body to maintain homeostasis [23]. At very high salinity levels, it is especially difficult to pump Na+ against the chemical gradient and out of the gill [65]. In a hypersaline environment, it therefore may be beneficial to maintain a lower epithelial salt permeability to avoid a back-flux of salts across the tight junctions of the gill epithelia [65, 66]. We found up-regulation of genes involved in cell-cell junction organization, including claudins involved in the formation of tight junctions (CLDN3, CLDN4, CLDN5, CLDN8, and CLDN10). Expression of gill claudin genes has been associated with salinity acclimation in fish [64, 66], and up-regulation of claudin-8 (CLDN8) has been shown to reduce the paracellular barrier permeability to Na+ [67], suggesting a role for regulation of gill epithelial permeability via tight junctions during osmoregulatory challenges. Our findings support the hypothesis that modifications to transmembrane transport and gill epithelial permeability jointly contribute to the salt exclusion necessary for survival in a hypersaline environment. Cell cycle regulation and cellular metabolic and biosynthesis pathways. Populations living in hypersaline environments must cope with the adverse effects of high salinity levels, often resulting in strategies for damage repair and energy reallocation [23, 62]. Functional analysis of the downregulated genes in the hypersaline population of L. perugiae showed signatures of downregulation of cell cycle and protein folding processes. Downregulation of genes involved in the cell cycle is expected to occur under high salinity stress to stop the replication of damaged cells and allow time for DNA repair [68, 69]. Similarly, locally adapted killifish living in freshwater and brackish environments exhibit evidence for divergence in the expression of genes underlying control of the cell cycle, supporting a central role in cell cycle regulation in adaptation to salinity stress [68, 70]. However, our findings may also suggest a general re-allocation of energy resources [68]. In hypersaline L. perugiae, the downregulation of genes involved in the cell cycle and protein folding could be reflective of an energetic trade-off between increased demand for processes involved in osmoregulation and the energetic cost of growth under stressful conditions [23]. In support of such energetic trade-offs, hypersaline L. perugiae have morphological differences from their freshwater counterparts, including a significantly smaller body size and reduced secondary sex features in males [23, 25]. Furthermore, energetic trade-offs associated with increased energy investment in osmoregulation in a hypersaline environment may also result in reduced investment in cellular metabolic and biosynthetic processes. In addition to genes involved in cell cycle processes, the weighted gene co-expression network analysis revealed that several cellular metabolic and biosynthetic processes were negatively associated with salinity (associated with genes in the blue module), including cellular aromatic compound metabolic process, organic cyclic compound metabolic process, nucleic acid metabolic process, cellular metabolic process, small molecule biosynthetic process, cellular macromolecule biosynthetic process, and carbohydrate biosynthetic process. Cellular metabolism and biosynthesis of compounds are energetically costly pathways for organisms [71], so these functions may receive less energy investment due to the high amounts of ATP required for ion transport and other osmoregulatory functions [68]. Future research will be important for identifying how osmoregulation in hypersaline environments may carry costs that influence cellular processes and how these costs impact organismal performance under different environmental conditions. Regulation of cell signaling. To respond to osmoregulatory challenges, fishes must perceive their environmental salinity and maintain intracellular signals that modulate ion transport and other processes [72, 73]. We found differentially expressed genes associated with a variety of cell signaling processes. Several upregulated genes in the hypersaline population were enriched in GO terms associated with cell signaling, including signal transduction, regulation of signaling, G-protein-coupled receptor signaling pathway, and second-messenger-mediated signaling. G-protein-coupled receptor signaling is among the pathways involved in allowing aquatic organisms to sense their environmental salinity [74, 75]. Some of these signaling pathway components may also play a role in regulating ion transport, such as mitogen-activated protein kinases (MAPKs) and other serine/threonine protein kinases [76–79], which were among the upregulated genes. MAPK signaling pathways are also implicated in regulation of the cell cycle, some of which trigger cellular growth arrest and DNA damage repair in response to osmotic stress [79, 80]. The upregulation of signaling pathways, especially MAPK genes, may consequently play a role in the downregulation of cell-cycle genes that we found in the hypersaline L. perugiae population. Comparisons among replicated lineages Lineage-specific responses to variation in salinity. Convergence in gene expression among independently evolved populations can occur in response to shared environmental stressors [81, 82]. Among fishes, cases of convergence in gene expression have been identified in response to selective pressures such as pollution [81], hydrogen sulfide [82], and the absence of light in caves [83]. In contrast, our comparative transcriptomics analysis indicated there is little evidence for convergence in gene expression patterns among fish lineages in response to variation in salinity. Even among the shared differentially expressed genes, there is considerable variation in the magnitude and direction of expression differences between ecotypes across lineages. This finding mirrors genomic analyses that looked for convergent signatures associated with adaptation to salinity variation, which found osmoregulation genes to be common targets of selection but also variation in the exact genes that were involved across different lineages [84]. Our finding of lineage-specific responses to variation in salinity, and therefore low levels of convergence in gene expression, could be explained by several non-mutually exclusive hypotheses [85]. First, changes in both protein-coding DNA sequences and gene expression may occur during adaptation to a novel environment, or they may occur independently [55, 86, 87]. Under certain conditions, selection may favor changes in protein structure and function via modifications at the sequence level, regardless of whether or not gene expression is affected [86]. Alternatively, changes in gene expression may be favored without selection acting at the sequence level [86]. Variation in how selection acts on gene expression could contribute to the lack of convergence among the lineages in our study, and there also may be stronger signals of convergence at levels other than gene expression. Secondly, differences in genetic architecture among lineages may cause different responses to selection, leading to diverse evolutionary outcomes [85]. Convergence is therefore less likely to occur with increasing genetic divergence between lineages [88], so it is not necessarily surprising that we did not find convergence among distantly related fishes. Even if the genetic architectures are similar among lineages, idiosyncratic responses to selection can also arise as a consequence of functional redundancy [89]. Specifically, modification of different genes and pathways may have equivalent functional and fitness consequences [78, 84]. Functional redundancy is particularly common in complex traits like salinity tolerance that involve many genes and physiological pathways [89]. While salinity tolerance is a shared outcome among the lineages in our analysis, the molecular mechanisms underlying osmoregulation may be unique to each lineage. Third, the directionality of habitat transition may impact gene expression responses, especially if there are genetic adaptations to the habitat of origin. For example, a high-activity version of an ion transporter may be downregulated during transitions to freshwater, while a low-activity version of the same enzyme may be upregulated in the opposite direction. Among the four lineages we included in our study, two of them transitioned from saltwater to freshwater environments (Gasterosteus and Odontesthes), and two transitioned from freshwater to saltwater environments (Limia and Leuciscus). If the direction of transition elicits similar responses, then we would expect populations that transitioned in the same direction to share more differentially expressed genes than those that did not transition in the same direction. However, we did not find this to be the case, as the pairs transitioning in the same direction share fewer differentially expressed genes than those that made opposite transitions (Fig 3C). Finally, covariation with other sources of selection—both abiotic and biotic—may cause idiosyncratic gene expression patterns across lineages. Beyond the challenges directly imposed by variation in salinity, such habitat transitions are often accompanied by other environmental challenges, such as variation in temperature, exposure to novel parasites, and restructuring of host-associated microbial communities [2, 20, 90, 91]. For example, hypersaline Amur ide must cope with high alkalinity stress in addition to salinity, and hypersaline L. perugiae have a warmer environment than freshwater L. perugiae [13, 23]. Additionally, there are differences in the gill microbial communities of saltwater and freshwater populations of South American silversides [20]. Such correlated environmental factors can contribute variation in gene expression responses to salinity transitions we observed. Shared responses involving ion transport and immune function. Although there is little evidence for convergence among all four lineage pairs included in our analysis, there was some functional overlap in the differentially expressed genes shared among three or more lineages. Most of the shared differentially expressed genes were associated with ion transport and immune system processes. Genes implicated in ion transport and immune system processes were also identified among the ten differentially expressed genes that were shared among all four lineages. For example, SLC9A3, a solute carrier gene involved in osmoregulation, was among the ten shared genes, suggesting its role in salinity transitions among fishes [59]. Another shared differentially expressed gene, CEACAM1 (Carcinoembryonic antigen-related cell adhesion molecule 1), has been implicated in immune processes in other systems and is thought to be involved in a variety of pathways, but its function is not well known [92, 93]. As previously discussed, regulation of ion transport is expected to be crucial when crossing a saltwater-freshwater boundary [16], and osmoregulation genes also exhibit evidence of convergent evolution during salinity transitions [84].Variation in expression of immune genes, particularly those related to inflammation and adaptive immunity, has also been documented in fish during salinity acclimation [94, 95]. Under a variety of selection pressures, locally adapted populations also frequently show divergence in immune genes due to other factors such as life history differences, different parasite exposure, and shifts in the microbiome [20, 95–98]. Immune loci are consequently evolutionary hotspots in diversification [99–101]. It needs to be tested whether changes in the expression of immune-related genes are directly linked to variation in salinity or whether these genes generally respond to changes in correlated biotic sources of selection. Overall, our analysis of gene expression patterns between locally adapted freshwater and hypersaline populations of L. perugiae provide insight into how this livebearing fish maintains homeostasis in a hypersaline environment. In addition, comparisons between four population pairs of freshwater and saltwater ecotypes in disparate teleost lineages showed little evidence for convergence, as there were only ten differentially expressed genes that were shared among them all. Despite this, we found that the differentially expressed genes shared in three or more of the lineages reflected biological processes related to ion transport and immune functioning. Our findings provide insight into shared and unique gene expression responses to salinity variation, broadly informing our understanding of salinity tolerance in aquatic organisms. Furthermore, comparisons in gene expression across species that have made habitat transitions are important for understanding mechanisms of adaptation to novel environments. Our results shed light on the repeatability of transcriptomic responses to salinity variation among fishes, which remains a relatively underexplored area of research despite its relevance for aquatic ecology and evolutionary biology. Future research comparing gene expression in fishes from freshwater and saltwater environments in both laboratory and field settings will be important for identifying the degree in which shared or unique gene expression responses are due to plasticity or adaptation, and additional studies that take a comparative transcriptomic approach with more lineages—representing transitions from freshwater to saltwater and vice versa as well as spanning a range of transition time and divergence time between populations—will provide further insight into shared mechanisms of osmoregulation among fishes. Lineage-specific responses to variation in salinity. Convergence in gene expression among independently evolved populations can occur in response to shared environmental stressors [81, 82]. Among fishes, cases of convergence in gene expression have been identified in response to selective pressures such as pollution [81], hydrogen sulfide [82], and the absence of light in caves [83]. In contrast, our comparative transcriptomics analysis indicated there is little evidence for convergence in gene expression patterns among fish lineages in response to variation in salinity. Even among the shared differentially expressed genes, there is considerable variation in the magnitude and direction of expression differences between ecotypes across lineages. This finding mirrors genomic analyses that looked for convergent signatures associated with adaptation to salinity variation, which found osmoregulation genes to be common targets of selection but also variation in the exact genes that were involved across different lineages [84]. Our finding of lineage-specific responses to variation in salinity, and therefore low levels of convergence in gene expression, could be explained by several non-mutually exclusive hypotheses [85]. First, changes in both protein-coding DNA sequences and gene expression may occur during adaptation to a novel environment, or they may occur independently [55, 86, 87]. Under certain conditions, selection may favor changes in protein structure and function via modifications at the sequence level, regardless of whether or not gene expression is affected [86]. Alternatively, changes in gene expression may be favored without selection acting at the sequence level [86]. Variation in how selection acts on gene expression could contribute to the lack of convergence among the lineages in our study, and there also may be stronger signals of convergence at levels other than gene expression. Secondly, differences in genetic architecture among lineages may cause different responses to selection, leading to diverse evolutionary outcomes [85]. Convergence is therefore less likely to occur with increasing genetic divergence between lineages [88], so it is not necessarily surprising that we did not find convergence among distantly related fishes. Even if the genetic architectures are similar among lineages, idiosyncratic responses to selection can also arise as a consequence of functional redundancy [89]. Specifically, modification of different genes and pathways may have equivalent functional and fitness consequences [78, 84]. Functional redundancy is particularly common in complex traits like salinity tolerance that involve many genes and physiological pathways [89]. While salinity tolerance is a shared outcome among the lineages in our analysis, the molecular mechanisms underlying osmoregulation may be unique to each lineage. Third, the directionality of habitat transition may impact gene expression responses, especially if there are genetic adaptations to the habitat of origin. For example, a high-activity version of an ion transporter may be downregulated during transitions to freshwater, while a low-activity version of the same enzyme may be upregulated in the opposite direction. Among the four lineages we included in our study, two of them transitioned from saltwater to freshwater environments (Gasterosteus and Odontesthes), and two transitioned from freshwater to saltwater environments (Limia and Leuciscus). If the direction of transition elicits similar responses, then we would expect populations that transitioned in the same direction to share more differentially expressed genes than those that did not transition in the same direction. However, we did not find this to be the case, as the pairs transitioning in the same direction share fewer differentially expressed genes than those that made opposite transitions (Fig 3C). Finally, covariation with other sources of selection—both abiotic and biotic—may cause idiosyncratic gene expression patterns across lineages. Beyond the challenges directly imposed by variation in salinity, such habitat transitions are often accompanied by other environmental challenges, such as variation in temperature, exposure to novel parasites, and restructuring of host-associated microbial communities [2, 20, 90, 91]. For example, hypersaline Amur ide must cope with high alkalinity stress in addition to salinity, and hypersaline L. perugiae have a warmer environment than freshwater L. perugiae [13, 23]. Additionally, there are differences in the gill microbial communities of saltwater and freshwater populations of South American silversides [20]. Such correlated environmental factors can contribute variation in gene expression responses to salinity transitions we observed. Shared responses involving ion transport and immune function. Although there is little evidence for convergence among all four lineage pairs included in our analysis, there was some functional overlap in the differentially expressed genes shared among three or more lineages. Most of the shared differentially expressed genes were associated with ion transport and immune system processes. Genes implicated in ion transport and immune system processes were also identified among the ten differentially expressed genes that were shared among all four lineages. For example, SLC9A3, a solute carrier gene involved in osmoregulation, was among the ten shared genes, suggesting its role in salinity transitions among fishes [59]. Another shared differentially expressed gene, CEACAM1 (Carcinoembryonic antigen-related cell adhesion molecule 1), has been implicated in immune processes in other systems and is thought to be involved in a variety of pathways, but its function is not well known [92, 93]. As previously discussed, regulation of ion transport is expected to be crucial when crossing a saltwater-freshwater boundary [16], and osmoregulation genes also exhibit evidence of convergent evolution during salinity transitions [84].Variation in expression of immune genes, particularly those related to inflammation and adaptive immunity, has also been documented in fish during salinity acclimation [94, 95]. Under a variety of selection pressures, locally adapted populations also frequently show divergence in immune genes due to other factors such as life history differences, different parasite exposure, and shifts in the microbiome [20, 95–98]. Immune loci are consequently evolutionary hotspots in diversification [99–101]. It needs to be tested whether changes in the expression of immune-related genes are directly linked to variation in salinity or whether these genes generally respond to changes in correlated biotic sources of selection. Overall, our analysis of gene expression patterns between locally adapted freshwater and hypersaline populations of L. perugiae provide insight into how this livebearing fish maintains homeostasis in a hypersaline environment. In addition, comparisons between four population pairs of freshwater and saltwater ecotypes in disparate teleost lineages showed little evidence for convergence, as there were only ten differentially expressed genes that were shared among them all. Despite this, we found that the differentially expressed genes shared in three or more of the lineages reflected biological processes related to ion transport and immune functioning. Our findings provide insight into shared and unique gene expression responses to salinity variation, broadly informing our understanding of salinity tolerance in aquatic organisms. Furthermore, comparisons in gene expression across species that have made habitat transitions are important for understanding mechanisms of adaptation to novel environments. Our results shed light on the repeatability of transcriptomic responses to salinity variation among fishes, which remains a relatively underexplored area of research despite its relevance for aquatic ecology and evolutionary biology. Future research comparing gene expression in fishes from freshwater and saltwater environments in both laboratory and field settings will be important for identifying the degree in which shared or unique gene expression responses are due to plasticity or adaptation, and additional studies that take a comparative transcriptomic approach with more lineages—representing transitions from freshwater to saltwater and vice versa as well as spanning a range of transition time and divergence time between populations—will provide further insight into shared mechanisms of osmoregulation among fishes. Supporting information S1 Fig. Biological processes from the GO enrichment analysis that were associated with upregulated genes in the saltwater population of L. perugiae. https://doi.org/10.1371/journal.pone.0315014.s001 (PNG) S2 Fig. Biological processes from the GO enrichment analysis that were associated with downregulated genes in the saltwater population of L. perugiae. https://doi.org/10.1371/journal.pone.0315014.s002 (PNG) S3 Fig. Biological processes from the GO enrichment analysis that were associated with differentially expressed genes shared among three or more lineages. https://doi.org/10.1371/journal.pone.0315014.s003 (PNG) S1 Appendix. Contains Tables A-E, which include the mapping statistics for all four population pairs, the weighted gene co-expression network analysis (WGCNA) results, and the significantly enriched GO terms from the WGCNA modules. https://doi.org/10.1371/journal.pone.0315014.s004 (XLSX) Acknowledgments We thank Ingo Schlupp for help with the field collections and the constructive feedback on the manuscript. We thank the authors of past salinity transcriptomics papers [13, 19, 20] who made their data publicly available and enabled the comparative analyses presented here.
TI - Gene expression signatures between Limia perugiae (Poeciliidae) populations from freshwater and hypersaline habitats, with comparisons to other teleosts
JF - PLoS ONE
DO - 10.1371/journal.pone.0315014
DA - 2024-12-05
UR - https://www.deepdyve.com/lp/public-library-of-science-plos-journal/gene-expression-signatures-between-limia-perugiae-poeciliidae-0s0vjkL9cE
SP - e0315014
VL - 19
IS - 12
DP - DeepDyve
ER -