Disparate responses to salinity across species and organizational levels in anchialine shrimps

Crustaceans are among the most well studied group for physiological responses to changes in environmental salinity (Henry et al., 2012; McNamara and Faria, 2012). Models of crustacean osmoregulation suggest that most marine euryhaline species act as osmoconformers in seawater but transition to hyper-osmoregulation upon exposure to more dilute waters (e.g. below 27‰; Henry, 2005). While likely adaptive, osmoregulation necessitates an energetically expensive physiological transformation on molecular, cellular and ultrastructural scales (Henry et al., 2012; McNamara and Faria, 2012). For instance, ATP is expended to power Na⁺/K⁺-ATPase (NKA), which pumps sodium and potassium against their concentrations and establishes the electrochemical gradient allowing ion transport (Henry et al., 2012). Unsurprisingly, during salinity transfer the expression of NKA and accessory enzymes like carbonic anhydrase (CA) tends to increase, often dramatically (Havird et al., 2013; Henry, 2001; Luquet et al., 2005; Serrano et al., 2007; Serrano and Henry, 2008). Mitochondrion-rich cells (MRCs), which provide the ATP required for osmoregulation (Copeland and Fitzjarrell, 1968; Havird et al., 2014c; Lovett et al., 2006), also proliferate during salinity transfer (Genovese et al., 2000; Taylor and Taylor, 1992). These processes occur in the gills, the primary site of osmoregulation. However, only posterior gills typically contribute to ion transport in many taxa, while anterior ones are specialized for respiration (Freire et al., 2008; Henry, 1984, 2001; Henry and Cameron, 1982; Neufeld et al., 1980; Taylor and Taylor, 1992).

However, this paradigm of physiological response by crustaceans to varying salinity may be incomplete for at least two reasons. Firstly, prior work has predominately focused on the decapod brachyurans (i.e. crabs), but this represents just a limited slice of ecological and phylogenetic diversity. For example, crabs usually only invade estuaries and streams seasonally or for reproduction (Warner, 1987). Therefore, salinity changes in this group may proceed in a gradual and predictable manner, which might not be representative for many taxa or environments. Consistent with this, alternative osmoregulatory strategies have been identified among other crustaceans. In freshwater prawns, uncharacteristic patterns of ion transporter activity, expression and localization were found during seasonal migrations (Faleiros et al., 2010; Huong et al., 2010; Maraschi et al., 2015; McNamara and Faria, 2012; McNamara and Lima, 1997; Ordiano et al., 2005). Likewise, specific gills involved in ion transport differ in freshwater crayfish and anomurans (Dickson et al., 1991; Wheatly and Henry, 1987) and at least one freshwater crab shows stable expression of NKA during salinity change (Tsai and Lin, 2007). Whether the taxa are of marine or freshwater ancestry may therefore be important. Regardless, the mechanisms of osmoregulation at whole-organism, tissue and molecular levels may be more diverse than previously believed.

Secondly, most previous studies have focused on individual genes functioning directly in ion transport. Consequently, possible alternative or accessory enzymes may have been largely overlooked. Supporting this, transcriptome-based studies have identified: (1) novel ion transporters with similar salinity-mediated expression patterns to NKA, (2) possible hormonal regulators, and (3) accessory systems involved in the energetics of osmoregulation (Chen et al., 2015; Havird et al., 2016; Hu et al., 2015; Li et al., 2014; Lv et al., 2013; Towle et al., 2011). Unfortunately, many RNA-Seq (RNA-sequencing) studies in crustaceans suffer from a lack of biological replication, making results from this technique potentially unreliable in crustacean ecophysiology (Havird and Santos, 2016b).

Shrimp species from the anchialine ecosystem are ideal models for further characterizing osmoregulatory mechanisms. Anchialine habitats are defined as landlocked ponds, pools and caves fed from both marine and freshwaters of subterranean origin (Holthuis, 1973; Sket, 1996). By definition they are characterized by extensive spatial and temporal fluctuations in salinity due to strong haloclines and tidal influences, respectively (Sket, 1996). As a result, anchialine crustaceans may often transition from full-strength seawater (i.e. 32‰) to freshwater (i.e. 0–5‰) and even hypersaline waters (i.e. 50‰). Shrimps have also independently invaded the anchialine ecosystem at least 3 times, from both freshwater and marine origins (von Rintelen et al., 2012) (Fig. 1). Moreover, previous work on the anchialine endemic shrimp Halocaridina rubra (Holthuis, 1963) suggests it may employ unusual mechanisms of osmoregulation: while H. rubra is among the strongest hyper-osmoregulators described, expression of typical osmoregulatory genes such as NKA and CA remains unchanged during salinity transfer (Havird et al., 2014c). Furthermore, all gills of H. rubra participate in ion transport. Together, these physiological attributes imply that the gills of H. rubra may constitutively express elevated levels of ion transporters and are specialized for osmoregulation (Havird et al., 2014c), possibly at the expense of gas exchange (Havird et al., 2014a).

The utilization of multiple osmoregulatory mechanisms within crustaceans raises intriguing hypotheses. For example, despite inhabiting the same ecosystem, taxa from different aquatic origins (i.e. freshwater versus marine) might exhibit different mechanisms of osmoregulation. Closely related taxa may also exhibit more similar responses to salinity compared with more distantly related species. Finally, differing patterns might be observed among species regardless of phylogenetic or ecological ancestry. To more thoroughly explore these hypotheses, we quantified osmoregulation, gill physiology and salinity-induced changes in gene expression across the gill transcriptome from five phylogenetically diverse shrimp species representing at least three independent invasions of the anchialine ecosystem (Fig. 1).

Five shrimp species were collected from habitats in the Hawaiian Islands and Okinawa, Japan (Fig. 1) as described previously (Havird et al., 2014a,c). Specifically, Halocaridina rubra (Hawaii), along with both Antecaridina lauensis (Edmondson, 1935) and Halocaridinides trigonophthalma (Fujino and Shokita, 1975) (Japan), belong to the ‘anchialine clade’ of the Atyidae (Fig. 1; von Rintelen et al., 2012). Given the freshwater ancestry of the family (Fryer, 1977; Glaessner, 1969; Smith and Williams, 1981), the common ancestor for members of this clade likely represents a single invasion of the anchialine ecosystem from a freshwater habitat. In contrast, Caridina rubella (Fujino and Shokita, 1975) (Japan) is the only member of the ‘Caridina-like clade’ of the Atyidae (Fig. 1) endemic to the anchialine ecosystem, thus implying a second, independent invasion by the family from freshwater. Lastly, Metabetaeus minutus (Whitelegge, 1897) (Japan) represents a third invasion, but likely from a nearshore marine habitat given the ancestry of most other alpheids (Anker, 2010). Importantly, each species encompasses multiple, genetically distinct lineages that likely represent cryptic species because of their high levels of molecular divergence (Craft et al., 2008; Weese et al., 2013). Thus, our sampling strategy (Fig. 1) targeted only a single genetic lineage per species in order to avoid any sources of within-species variation (Havird et al., 2014a).

Salinity was quantified at each sampling location as previously described (Havird et al., 2014a,c). Animals were transported either to Auburn University (Auburn, AL, USA – for H. rubra) or the Sesoko Marine Station – Tropical Biosphere Research Center (Okinawa, Japan – for the other four species) following collection and held in aquaria at ∼32‰ (with the exception of C. rubella – see Results) for 1 month prior to salinity transfers to ensure chronic acclimation (Jillette et al., 2011).

Quantification of hemolymph osmolality during salinity transfer was conducted for the four species from Japan as previously described for H. rubra (Havird et al., 2014c). Briefly, shrimp were anesthetized on ice, excess water wicked away and hemolymph extracted via dissection and subsequent centrifugation. Hemolymph samples were then frozen at −80°C, transported on ice to Auburn University, and osmolality measured on a vapor pressure osmometer (Wescor 5100C, Logan, UT, USA) from 10 μl samples prepared from either single (i.e. C. rubella) or the pooling of 1–3 (i.e. A. lauensis, H. trigonophthalma and M. minutus) individuals per sample. Hemolymph osmolality was quantified at a minimum when shrimp were chronically acclimated to high (32‰) versus low (2‰) salinity. Furthermore, H. rubra is known to alter hemolymph osmolality within 48 h of transfer (Havird et al., 2014c). Therefore, a minimum chronic acclimation duration of 48 h was utilized (though some individuals were acclimated for up to 12 days prior to hemolymph extraction). Because of limitations associated with specimen availability, the details of each experiment (sample size, the number of acclimation salinities examined and the number of time points utilized) varied among species (see below).

To visualize the proportion of gills participating in ion transport, silver nitrate (AgNO₃) was used, which stains ion-transporting epithelia black (Croghan, 1958; Holliday et al., 1990; Kikuchi and Shiraishi, 1997). Prior to staining, shrimp were chronically acclimated for at least 10 days to either high (32‰, except 25‰ for C. rubella) or low (2‰, except 4‰ for M. minutus) salinity prior to staining. While AgNO₃ staining specifically indicates chloride transport across the gill, it is possible changes in pH or ammonia concentration could influence ion transport (Henry et al., 2012). However, as these parameters were kept constant in the different salinity treatments, it is unlikely they influenced ion transport. With the limited availability of A. lauensis, ion transport staining at multiple salinities was not possible, although staining at 2‰ for three individuals was performed to visualize gill morphology. Stain concentrations, incubation times and washes followed the protocol developed for H. rubra (Havird et al., 2014c). Photographs of stained gills were captured with a dissecting microscope and analyzed using ImageJ version 1.45s (National Institutes of Health, Bethesda, MD, USA) in order to quantify the relative area stained (area fraction) of each gill. Statistical effects of salinity, gill number and their interaction on the area fraction stained of each gill for each species were assessed with a linear model using the lm function in R version 3.3.1 (http://www.R-project.org/).

Assembly of reference transcriptomes was previously described for H. rubra (Havird and Santos, 2016a) and the four species from Japan (Havird et al., 2014b). Briefly, total RNA was extracted from whole animals <1 month after preservation in RNAlater via bead-beating in Trizol and using an RNeasy Kit (Qiagen). Synthesis of cDNA libraries utilized a SMART cDNA construction kit (ClonTech) and libraries were sent to the Genomic Services Lab at the HudsonAlpha Institute for Biotechnology (Huntsville, AL, USA) for indexing and subsequent sequencing on an Illumina HiSeq 2000, with 100 bp paired-end reads delivered in FASTQ format. Raw reads were digitally normalized and assembled using Trinity version Trinityrnaseq_r20131110 (Grabherr et al., 2011) with default parameters. The ‘completeness’ of the five reference transcriptomes was assessed with BUSCO version 3 to search for 978 single-copy orthologs conserved across metazoans (Simão et al., 2015). Transcriptomes are available from the Transcriptome Shotgun Assembly Sequence Database (NCBI TSA; accession numbers in Table S1 of figshare dataset – see ‘Data availability’ section, below) and Dryad (doi:10.5061/dryad.jn8t1) while raw reads are available from the Sequence Read Archive (NCBI SRA; accession numbers in Table S2 of figshare dataset – see ‘Data availability’ section, below).

In general, environmental salinity differed among and within the anchialine habitats where these species were collected. However, decapod crustaceans make the molecular transition from osmoconformity to osmoregulation at a critical salinity, typically 26 ppt (Henry, 2005). Thus, shrimp were chronically acclimated for at least 1 month to a high ‘reference’ salinity and salinity changes were all designed to cross the critical salinity and induce the physiological and molecular transition between conformity and regulation. Because of differences in specimen availability and other logistical constraints, exact details of the transfers and sampling differed somewhat between species (see Table S3 of figshare dataset – ‘Data availability’ section, below). Briefly, a single transfer from 32‰ to 15‰ was used in H. rubra with sampling occurring at multiple time points after transfer, while in the Japanese species, multiple transfers were investigated but only at a single time point (24 h after transfer). Only intermolt individuals were used and shrimp were not fed prior to sampling as four of the five species (exception being M. minutus) are bacterial/algal grazers that forage substrate in their immediate environment. Water chemistry was monitored regularly prior to transfers to ensure pH and ammonia levels were similar among species. At each sampling point after salinity transfer, all gills were dissected and pooled as described previously (Havird et al., 2014c), although different pooling schemes were used because of differences in body size among species. Control individuals transferred to the same salinity as the reference were also sampled 24 h after transfer.

We utilized a cost-effective 3′-tag-based RNA-Seq approach, in which the terminal ∼50 bp of transcript are sequenced, to identify differentially expressed genes (DEGs) during salinity transfer (Meyer et al., 2011). For the four species from Japan, gills for RNA-Seq were preserved in RNAlater and stored at −80°C for ∼1 month, while for H. rubra, processing occurred immediately after dissection. RNA extraction from pooled gills followed a protocol similar to that for the reference transcriptomes, except that an E.Z.N.A. Total RNA Kit I (Omega Biotek) was employed. Quality and quantity of extracted RNA were assessed via a Tapestation (Agilent) and cDNA libraries were prepared as described previously (Meyer et al., 2011). Briefly, RNA was fragmented at 95°C for 20 min and first-strand cDNA synthesized using Superscript II reverse transcriptase (Thermo Fisher Scientific) and an oligo-dT primer to target 3′ ends. PCR was then performed with an initial denaturing step of 95°C for 5 min and 19 cycles of 95°C for 40 s, 65°C for 1 min and 72°C for 1 min. Amplified products were purified with a QIAquick PCR purification kit (Qiagen) followed by adapter ligation and gel electrophoresis for size selection (200–300 bp). Resulting cDNA libraries were sequenced as single-end 50 bp reads (Table S4 of figshare dataset – see ‘Data availability’ section, below) at the HudsonAlpha Genomic Services Lab on an Illumina HiSeq 2000 and are available from NCBI's SRA (accession numbers in Table S2 of figshare dataset – see ‘Data availability’ section, below).

Reads were screened to remove short and low-quality sequences, adapter contamination and homopolymer repeats using custom scripts (available from https://github.com/Eli-Meyer). Two different protocols were then utilized in identifying DEGs. The first followed Meyer et al. (2011), where the gmapper script from the SHRiMP package (Rumble et al., 2009) was used to initially map filtered reads to each reference transcriptome and the subcomponent_filter_and_count.pl script (found at https://github.com/Eli-Meyer) was used to generate a counts file per pooled gill sample. These were then combined into a single file containing read numbers mapping per treatment to each contig in the reference. Importantly, all ‘subcomponents’ for a given contig in the reference transcriptomes were collapsed to their ‘components’ so as to exclude isoforms, or alternative transcripts for the same gene, from the differential expression analyses. Finally, BioConductorR version 3.7 in R was then utilized to normalize read counts, generate FPKM values (fragments per kilobase of transcript per million mapped reads), and identify DEGs via DESeq (https://bioconductor.org/packages/release/bioc/html/DESeq.html) for pairwise comparisons of each experimental salinity compared with the reference salinity (minimum log₂ fold change=1, adjusted P<0.05, false discovery rate FDR<0.001) (Anders and Huber, 2010).

For some pairwise comparisons, a small number of DEGs (i.e. 12 or fewer) were recovered from these initial analyses, likely because default parameters in DESeq are relatively conservative (Love et al., 2014), as has been previously demonstrated in crustacean RNA-Seq experiments (Havird and Santos, 2016b). To supplement these DESeq analyses, DEGs were also identified with the ‘prepackaged’ analytical pipeline provided with Trinity and outlined in Haas et al. (2013). In this case, filtered reads were mapped to reference transcriptomes using Bowtie version 1.1 (Langmead, 2010), TMM (trimmed mean of M-values)-normalized counts and FPKM values generated with RSEM version 1.2.12 (Li and Dewey, 2011), and DEGs identified using DESeq2 (Love et al., 2014) and edgeR (Robinson et al., 2010) under the same statistical parameters as above.

Statistically significant DEGs were submitted to the Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.7 (Huang et al., 2007a,b) as described previously (Havird and Santos, 2016a) to identify enriched biological themes. Briefly, Uniprot accession numbers extracted from annotation reports were submitted to DAVID and significant enrichment for gene ontology (GO) terms assessed via a Fisher's exact test with a FDR of 0.05 (Hosack et al., 2003). Results were visualized as networks by inputting the ‘Functional Annotation Chart’ from DAVID into Cytoscape version 3.2 (Shannon et al., 2003) and utilizing the Enrichment Map version 2.1.0 (Merico et al., 2010) and WordCloud version 3.0 plugins with P- and FDR-value cutoffs of 0.05 and 0.1, respectively.

To further characterize significantly upregulated or downregulated genes and enriched biological themes, the expression of genes predicted a priori to be involved in crustacean osmoregulation was specifically examined, including ion transporters and accessory enzymes that are generally upregulated during salinity transfer: (1) NKAα; (2) NKAβ; (3) the Na⁺/K⁺/2Cl⁻ cotransporter (NKCC); (4) the V-type H⁺-ATPase (HAT); (5) CA; (6) Na⁺/H⁺ exchangers; (7) arginine kinase (ArgK); and (8) Na⁺/HCO₃⁻ cotransporters (Havird et al., 2013). Genes often upregulated as a result of general cellular stress, such as heat shock proteins (HSPs), superoxide dismutases (SODs), and glutathione peroxidases (GPXs), were also investigated. Finally, mitochondrial- and nuclear-encoded oxidative phosphorylation (OXPHOS) subunits that comprise the mitochondrial electron transport chain were examined given their upregulation in previous transcriptome-wide analyses of crustaceans (Havird et al., 2016; Towle et al., 2011): NADH dehydrogenase (complex I), cytochrome bc₁ (complex III), cytochrome c oxidase (complex IV) and ATP synthase (complex V).

Contigs corresponding to specific genes of interest were identified from the transcriptomes via a combination of keyword screening of annotation reports and similarity searches with BLAST version 2.2.29 utilizing known homolog sequences from closely related species as queries. Expression values, as TMM-normalized FPKM, were then extracted for each contig. Contigs with low expression values (i.e. an average FPKM<10 across all treatments for a given species) were excluded, with an exception made for subunit 2 of cytochrome bc₁ from M. minutus in order to have representation for this gene from all five species (see below). When multiple contigs were recovered for the same gene from a single species, only the one with the highest average expression value was retained.

To investigate whether taxa exhibited overall similarities or differences in salinity-induced gene expression among shared genes, a common set of orthologs present in all five species was examined. To supplement a curated set of orthologs previously generated for the four species from Japan (Havird et al., 2014b), orthologs from H. rubra were identified. Specifically, the corresponding H. trigonophthalma sequence for each ortholog was used as a query in tblastn version 2.2.29 searches against the H. rubra transcriptome with an e-value cutoff of 1e⁻²⁰. Each best BLAST hit was considered an ortholog from H. rubra if during reciprocal tblastx searches the original query sequence was identified as the best hit against the H. trigonophthalma transcriptome.

Similarities and differences in expression among shared orthologs were visualized using between-group analysis (BGA) (Culhane et al., 2002). Specifically, BGA is a type of multivariate discriminant analysis employable when the number of cases (i.e. orthologs) exceeds that of variables (i.e. treatments) (Ghalambor et al., 2015). BGA was implemented for comparisons across treatments from the five species using the made4 version 3.9 package (Culhane et al., 2005) in R.

The four species from Japan are strong hyper-osmoregulators, maintaining similar osmotic gradients in freshwater (∼750–1200 mOsm kg⁻¹ H₂O) as H. rubra (870 mOsm kg⁻¹ H₂O) (Fig. S1 of figshare dataset – see ‘Data availability’ section, below). Furthermore, these species also act as osmoconformers at salinities approaching seawater, as in H. rubra and other euryhaline crustaceans. However, the four species from Japan showed much larger variances in hemolymph osmolalities than H. rubra (Fig. S1 of figshare dataset – see ‘Data availability’ section, below).

While H. trigonophthalma, A. lauensis and M. minutus survived salinity transfer from freshwater to seawater and 45‰ (as in H. rubra; Havird et al., 2014c; Holthuis, 1973), C. rubella (collected at 0‰) exhibited total mortality <24 h after transfer from 0‰ to 45‰ or 32‰. Given this, specimens were chronically acclimated to 25‰ instead of 32‰. Following acclimation to 25‰, however, 100% survivorship occurred for at least 48 h after transfer to 32‰ (but not 45‰), thus allowing sampling for RNA-Seq experiments as with the other species.

Variation among which gills participated in ion transport and were salinity responsive was found among the species studied here (Fig. 2). All gills in H. rubra undergo ion transport regardless of salinity (Fig. 2A; Havird et al., 2014c). In H. trigonophthalma, gills showed the most staining (60–85%) under high salinity, with the anterior-most pair also having a large fraction stained (80%) regardless of salinity, resulting in a significant interaction between gill number and salinity (Fig. 2B; P=0.011). In contrast, the highest staining in C. rubella gills was at low salinity (50–75% versus 20–50%, P<0.001), and although posterior rather than anterior gills were more affected by salinity, the interaction between gill number and salinity was not significant (Fig. 2C; P=0.225). In M. minutus, staining was also most prominent in low salinity (70% versus 0–50%, P<0.001, linear regression; Fig. 2D), and the posterior gills responded the most to salinity change (interaction P<0.001). However, the maximal fraction stained under any treatment was 70–80% in all species.

Gill morphology also differed among the examined species (Fig. S2 of figshare dataset – see ‘Data availability’ section, below). In H. rubra, the phyllobranchiate gills possess thick, lobe-like lamellae, and the gills themselves are quite simple compared with previously studied euryhaline decapods (i.e. 10–16 lamellae per each of four gill pairs in H. rubra versus ∼300 lamellae per each of eight gills in crabs; Havird et al., 2014c; Fig. S2 of figshare dataset – see ‘Data availability’ section, below). Gill morphology of H. trigonophthalma was similar to that of H. rubra, with each of four thicker, lobe-like gills having 22–26 lamellae (Fig. S2 of figshare dataset – see ‘Data availability’ section, below). Antecaridina lauensis had similar gills, with 9–31 lobe-like lamellae per each of six gills. In contrast, gills of C. rubella were morphologically complex, with 31–53 lamellae per each of six gills. Gills of M. minutus also possess 19–25 thinner, sheet-like lamellae packed tightly together per each of five gills.

Transcriptomes constructed from an adult plus pooled larval sample of H. rubra (i.e. the ‘adult+larvae transcriptome’ in Havird and Santos, 2016a) or whole specimens of the four species from Japan (Havird et al., 2014b) contained ∼20–138 (×10³) contigs assembled from ∼41–59 (×10⁶) 100 bp paired-end reads (Table 1). The transcriptomes for C. rubella and M. minutus had considerably fewer contigs than the others and BUSCO analyses suggested these transcriptomes were less complete (33%–51%) compared with those of the other species (>90%, Table 1).

On average, 8.9±0.4 (×10⁶) (mean±s.e.m.) 3′-tag reads were generated for each of 90 gill samples generated for salinity transfer experiments (Table S4 of figshare dataset – see ‘Data availability’ section, below). Filtering and trimming resulted in discarding ∼8% of reads on average and subsequent mapping rates to the reference transcriptomes averaged 66% per sample. Rates of mapping were notably lower in C. rubella and M. minutus.

The number of significant DEGs identified using DESeq varied substantially among species and treatments (Table 2). For example, in C. rubella 559 DEGs were identified after transfer from 25‰ to 32‰ while no DEGs were identified after transfer to 15‰. Furthermore, low numbers of DEGs (≤12) were identified for six comparisons (Table 2). Additional analyses using DESeq2 and edgeR (data not shown) boosted the number of DEGs for some of these comparisons (Table 2).

For H. rubra, functional categories could not be assigned for many DEGs because of an overall low annotation rate, which is typical for crustacean transcriptomes (Das et al., 2016; Havird and Santos, 2016b). For instance, only 150 (38%) of the 395 unique DEGs identified were annotated. Visualized networks generated from the DAVID analysis revealed similarly enriched biological themes at 3, 8 and 48 h post-transfer (Fig. 3). These included upregulated calcium-binding and calcium-transport motifs as well as genes involved in muscle function and cell differentiation at each time point. The number of nodes increased in the corresponding subnetworks with time following transfer. A subnetwork corresponding to hemocyanin, oxygen transport and other extracellular functions was downregulated across the same time points, with the number of nodes peaking 8 h post-transfer. Other subnetworks were ephemeral at various time points. For example, a large subnetwork characterized by chitin-related processes was downregulated at 48 h post-transfer, while single nodes corresponding to the major facilitator superfamily and mitochondrial genome maintenance were upregulated at 3 and 8 h post-transfer, respectively.

Additional functional patterns emerged when specific annotated DEGs were examined (Fig. S3 of figshare dataset – see ‘Data availability’ section, below). Unexpectedly, transcripts related to vision and light perception were among those most differentially expressed in the gills of H. rubra, including two rhodopsin and four arrestin homologs. These encompassed at least five of the 10 highest upregulated, annotated DEGs at each time point (in many cases over 100-fold). Related DEGs, including retinol dehydrogenase and retinoid isomerohydrolase, were significantly downregulated (∼10-fold) at later time points.

Other highly upregulated transcripts in H. rubra were related to calcium transport, consistent with the DAVID analyses. These included genes possessing the major facilitator superfamily domain, which were not detected at 32‰ in some comparisons (i.e. thus considered ‘inestimably’ upregulated). Other upregulated transcripts attributable to ion transport were: (1) a putative inorganic phosphate cotransporter (38-fold); (2) a member of the solute carrier organic anion transporter family (71-fold); (3) an organic cation transporter protein (53-fold); (4) a Na⁺/Ca²⁺ exchanger (7-fold); (5) an aquaporin (4-fold); and (6) various calcium channels. An innexin was also downregulated 7-fold 48 h post-transfer. Importantly, ion transporters and accessory genes typically implicated in crustacean osmoregulation, including NKA, CA, NKCC and HAT, were not among the identified DEGs for H. rubra.

Additional DEGs of note in H. rubra included upregulation of multiple structural genes like actins, myosins, titin and troponin (up to 16-fold) as well as two HSPs (HSP22 and HSP90α, 8- and 3-fold, respectively). Again consistent with the DAVID analyses, hemocyanin and oxygen-transport transcripts were downregulated at most time points, including multiple hemocyanin chains (up to 26-fold), a hemolymph clottable protein (5-fold) and a heme-binding protein (13-fold). Also downregulated at all time points were vitellogenesis-related transcripts, including a vitelline membrane outer layer protein (5-fold), von Willebrand factors (3-fold) and several vitellogenin homologs (i.e. considered ‘inestimably’ downregulated). Genes generally involved in lipid and fatty acid metabolism were also downregulated. Mitochondrial OXPHOS transcripts, including ATP and NADH subunits, were also downregulated 8 h post-transfer (3-fold), as well as several chitin-related transcripts at later time points (7-fold). Finally, potential regulatory genes such as histone-related transcripts, prostaglandin synthases, an octopamine receptor and a diuretic hormone receptor were also identified as DEGs.

For H. trigonophthalma, 253 significant DEGs were recovered following transfer from 32‰ to 15‰ and 45‰, although only 69 (27%) were annotated. Multiple mitochondrial OXPHOS subunits (up to 16-fold) and ribosomal proteins (up to 15-fold) were upregulated during transfer to 15‰. mTOR, a sulfotransferase, and USMG5 (another OXPHOS subunit) were the most highly upregulated (up to 48-fold). Downregulated DEGs (from 4-fold to ‘inestimably’ downregulated) for the transfer to 15‰ included: (1) ribosomal subunits; (2) HSPs; (3) a SOD; (4) a hemocytin; (5) a subunit of HAT; (6) malate dehydrogenase; (7) an ammonium transporter; (8) a calcium transporter; and (9) apolipophorins. Similar DEGs were identified from the transfer to 45‰, with additional downregulation of structural genes (i.e. myosin, actin and tubulin – up to ‘inestimably’ downregulated), NKCC (9-fold) and a mitochondrial peroxide reductase (11-fold). Only two annotated DEGs were identified from the transfer to 2‰ (USMG5 and the immune gene NLRC3, upregulated and downregulated 5- and 2027-fold, respectively). Additional DEGs recovered for the transfer to 2‰ using DESeq2 generally fell into the categories mentioned above.

For A. lauensis, only 12 DEGs were recovered for the transfer from 32‰ to 2‰ (Table 2). Only two were annotated: an atrial natriuretic peptide-converting enzyme (‘inestimably’ upregulated) and a collagen protein (5-fold upregulated). No DEGs at all were recovered from A. lauensis when DESeq2 was utilized.

For C. rubella, hundreds of DEGs were recovered during salinity transfer from 25‰ to both 2‰ and 32‰ (Fig. S4 of figshare dataset – see ‘Data availability’ section, below). Notably, NKAα was one of the most highly upregulated (13-fold) transcripts after the transfer to 2‰. Additional upregulated genes were several cytoskeleton components (i.e. titin, troponin, gelsolin and nesprin – up to 9-fold), two subunits of the 40S ribosome (5-fold) and a heme-binding protein (5-fold). Downregulated genes following the transfer to 2‰ included several mitochondrial ribosome subunits (up to 17-fold), OXPHOS complex subunits (up to 9-fold), hemocyanin chain B (51-fold), clotting factor B (233-fold) and a HSP (9-fold). For the transfer to 32‰, similar transcripts were either upregulated (i.e. NKAα, cytoskeleton components, ribosomal subunits, a voltage-dependent calcium channel and the heme-binding protein) or downregulated (i.e. mitochondrial ribosome subunits, OXPHOS complex subunits, hemocyanin chains and the same heat shock protein). Additional downregulated genes included cytoskeleton proteins, a retinol dehydrogenase and NKCC (6-fold). A different transcript also annotated as NKAα was downregulated 10-fold during the transfer to 32‰, in contrast to the one that was upregulated 8-fold. For the 25‰ to 15‰ transfer, significant DEGs were only identified when DESeq2 was employed. These DEGs included upregulated ribosomal subunits (up to 31-fold) and an OXPHOS subunit (32-fold), along with downregulated cytoskeletal proteins (up to 122-fold) and a ribosomal subunit (24-fold).

For M. minutus, very few DEGs were initially identified using DESeq (Table 2). The only annotated ones were a solute carrier protein, a serine protease and a mitochondrial sarcosine dehydrogenase, which were all upregulated (up to 11-fold) 24 h following transfer from 32‰ to 2‰. Analyses with DESeq2 recovered dozens of DEGs for each comparison (Table 2), including upregulation (up to 40-fold) of tropomyosin, 60S ribosomal protein L13 and a small EDRK-rich factor 2 following transfer to lower salinities. Transcripts involved in mitochondrial function (e.g. OXPHOS subunits, subunit beta of the mitochondrial-processing peptidase and a mitochondrial ribosomal subunit) were downregulated (up to 28-fold) after transfer to lower salinities. Several of these same genes were also differentially expressed after transfer to 45‰. Subunit G of HAT was upregulated and downregulated 16- and 9-fold during transfer to low and high salinity, respectively.

When specific transcripts of interest were examined, dissimilar patterns of salinity-induced expression were apparent among species (Fig. 4; Fig. S5 of figshare dataset – see ‘Data availability’ section, below). For example, NKAα expression increased during transfer of A. lauensis, H. trigonophthalma and M. minutus from higher to lower salinities. However, NKAα expression decreased for C. rubella and remained essentially unchanged in H. rubra. Differences in salinity-induced expression among the species were also observed for transcripts of other ion transporters and accessory enzymes, general stress-response genes and OXPHOS subunits (Fig. 4; Fig. S5 of figshare dataset – see ‘Data availability’ section, below).

A total of 631 shared orthologs were identified for the five species. In the BGA, principal component 1 (PC1) explained 45.9% of the observed variance in ortholog gene expression, while PC2 explained 21.9% (Fig. 5). Groupings tended to occur by species and not shared salinity treatments, with no overlap among species (Fig. 5A). Moreover, salinity transfer affected the expression of these orthologs differentially among species. For H. rubra and H. trigonophthalma, treatments were generally separated on PC2, while for C. rubella and M. minutus, little separation occurred by differing treatment. The trajectories of change in ortholog expression were also variable among species (Fig. 5B). For example, during transfer from higher to lower salinities, values generally decreased on PC2 and then decreased on PC1 at later time points for H. rubra. For the other atyid species, however, these same values generally increased on PC2 and/or PC1. Thus, each species exhibited a distinct response, both in magnitude (Fig. 5A) and direction (Fig. 5B), in expression of shared orthologs, despite similar salinity transfers.

Here we show that different shrimp species from the anchialine ecosystem exhibit diverse responses to changing salinity at whole-animal, tissue and molecular scales. Anchialine crustaceans are ideal models to study diversity in the phenotypic response to salinity. First, anchialine crustaceans must contend with habitats that experience unpredictable temporal and spatial shifts in salinity. This is in contrast to well-studied species like blue crabs, which often undergo predictable, gradual changes in salinity. Second, previous studies in H. rubra showed unusual patterns of salinity-induced gene expression relative to previously examined crustaceans (Havird et al., 2014c). Finally, phylogenetic analyses suggest the anchialine ecosystem has been independently invaded multiple times by lineages from the Atyidae and Alpheidae representing freshwater and marine ancestry, respectively (Fig. 1; Fig. S7 of figshare dataset – see ‘Data availability’ section, below; von Rintelen et al., 2012), allowing us to test the possibility that differing ecological ancestries may manifest in differing strategies for osmoregulation.

All five species of anchialine shrimps examined here were strong osmoregulators in freshwater (Fig. S1 of figshare dataset – see ‘Data availability’ section, below). For context, the euryhaline crab Callinectes sapidus is generally considered a strong hyper-osmoregulator and maintains an osmotic gradient of ∼600 mOsm kg⁻¹ H₂O in freshwater (Cameron, 1978; Henry, 2001). Osmotic gradients of up to ∼1250 mOsm kg⁻¹ H₂O were measured in the anchialine species examined here. Unfortunately, relatively small sample sizes and potential evaporation during shipping likely resulted in the large variances in hemolymph osmolality for the four species from Japan (variances were small in H. rubra). Thus, while anchialine shrimps fit the general patterns of osmoregulation seen from other euryhaline crustaceans, these results should not be viewed as quantitative certainties of osmotic gradients in these species.

Unexpectedly, not all anchialine shrimps examined here survived the investigated salinity transfers. Specifically, C. rubella was not able to survive in seawater except with gradual acclimation. Unlike the others, this species was also unable to tolerate hypersaline waters (45‰) even with gradual acclimation. Both C. rubella and members of the anchialine clade likely invaded the anchialine ecosystem from freshwater (Fryer, 1977; Glaessner, 1969; Smith and Williams, 1981). However, C. rubella is the only member of the Caridina-like clade to have done so, likely representing a more recent invasion relative to the anchialine clade (von Rintelen et al., 2012). Caridina rubella has an apparently limited distribution, with adults being reported from just six islands in the Ryukyus and larvae only recently reported from Okinawa-jima (Fujino and Shokita, 1975; Fujita et al., 2019; Weese et al., 2012, 2013). Together, this suggests C. rubella adults may not be physiologically specialized to contend with variable salinities, but instead may be restricted to the freshwater portions of anchialine habitats. Traditionally, freshwater species are limited in their invasion of higher salinities by restrictions on their ability to increase the intracellular pool of osmolytes that would keep cell volume constant (e.g. Pierce and Amende, 1981). Genes responsible for amino acid synthesis and transamination were not upregulated during salinity transfer in C. rubella. However, it is expected that these processes would play out in the internal soft tissues (e.g. muscle and nerve), which need to retain normal cell volume in order to maintain function, rather than in the gills. In fact, higher survivorship during the gradual transfer in C. rubella suggests that the longer time periods allowed for the synthesis of such intracellular osmolytes. Examining differences in salinity-induced gene expression across other tissues might address the processes behind such species-specific patterns of salinity tolerance. Furthermore, if C. rubella is not physiologically capable of rapidly invading seawater, behavioral rather than physiological adaptations may have facilitated its invasion of the anchialine ecosystem. Unfortunately, given the predicted rises in salinity for anchialine habitats in general under current climate change scenarios (Marrack, 2016), this species may face possible expatriation or extinction from this ecosystem in the future.

Gill morphology and patterns of salinity-induced ion transport varied widely among species (Fig. 2; Fig. S2 of figshare dataset – see ‘Data availability’ section, below). In general, gill morphology was similar among species from the anchialine clade, but differed substantially from that of the more distantly related C. rubella and M. minutus, which had morphologies reminiscent of previously studied euryhaline crabs (Freire et al., 2008; Henry, 1984; Henry and Cameron, 1982; Neufeld et al., 1980; Taylor and Taylor, 1992). Salinity transfer tended to affect ion transport in posterior gills more than anterior ones across all the species (Fig. 2), although in H. rubra this pattern was fairly subtle relative to that in the other species. As predicted, lower salinity resulted in higher levels of ion transport in most species. However, the reverse pattern was observed in H. trigonophthalma. In all species, at least some treatments resulted in a very high proportion of the gills being used for ion transport (e.g. >75% compared with 52% in the euryhaline crab Chasmagnathus granulatus; Genovese et al., 2000). This suggests that despite differences among species, gills in most anchialine shrimps may be specialized for ion transport over other functions such as respiration, as appears to be the case in H. rubra (Havird et al., 2014d).

Salinity-induced gene expression also varied extensively among species. Similar salinity transfers resulted in anywhere from zero to hundreds of identified DEGs depending on the species (Table 2). The types of DEGs also varied. For instance, opsins and other vision-related genes were among the most upregulated in H. rubra (Fig. 3) but were not upregulated in any other species. The few categories of genes that were shared across species often showed differing patterns of expression. For example, OXPHOS genes were upregulated in H. rubra but downregulated in H. trigonophthalma, C. rubella and M. minutus after similar salinity transfers. Ribosomal components, structural genes such as actins, and genes related to oxygen transport were also differentially expressed across multiple species, but inconsistently. Similarly, genes that were predicted a priori to be important in euryhalinity showed species-specific expression patterns (Fig. 4). Ion transporters that are generally upregulated with salinity transfer (Havird et al., 2013) were upregulated in some species, downregulated in others, or remained unchanged – all in response to the same salinity transfer (e.g. NKAα after transfer to 15‰). Echoing previous qPCR results (Havird et al., 2014c), these genes had stable expression in H. rubra and were never identified as DEGs. Shared orthologs among the species also showed species-specific responses to salinity transfer, in both the magnitude and direction of change (Fig. 5). Although only a small minority of orthologs are likely to be salinity responsive, the lack of overlap in expression of these genes among species adds to the general finding that molecular responses to salinity transfer are highly species specific among these anchialine shrimps.

Why do anchialine shrimps show such varied responses to changes in salinity? Firstly, differences in methodological details among species may have played a role. For example, RNA was extracted from H. rubra immediately following treatment, while the other species were preserved in RNAlater for pre-processing transport. This may have contributed to artifacts like the incomplete reference transcriptome assemblies for M. minutus and C. rubella (i.e. RNAlater might not have penetrated the entire animal in these larger species), thus hindering our ability to detect some DEGs. However, this does not explain differences among species in salinity-induced ion transport in the gills (Fig. 2), expression of genes of interest shared across all transcriptomes and species (Fig. 4), or expression across all shared orthologs (Fig. 5). Similarly, environmental salinity differed among the habitats sampled here (see below) and technical constraints prevented utilizing identical salinity transfers or examining equivalent time points for every species. However, it is the critical low salinity at which the physiological transition from conformity to regulation takes place that governs the upregulation of osmoregulatory genes and not the absolute magnitude of the salinity transfer (Henry, 2005; Mitchell and Henry, 2014a). Moreover, salinity-induced changes in gene expression appear to be an all-or-none response (Mitchell and Henry, 2014a). Because transfers were generally from salinities where species act as osmoconformers to those where they act as osmoregulators, similar patterns of gene expression were expected. Finally, relatively small sample sizes in A. lauensis (n=3) may have resulted in the small number of DEGs compared with that in other species, although small numbers of DEGs were also found in some treatments with large sample sizes, suggesting variation in sample size does not fully explain differences in the number of DEGs.

Putting aside methodological discrepancies, several biological hypotheses may explain species-specific responses to salinity in anchialine shrimps. Although all species examined here are anchialine endemics, individual anchialine habitats can show remarkable heterogeneity. Habitats may range in size from small cracks and fissures accessing the groundwater, to lakes and caves. Although detailed salinity regimes are not available for most anchialine habitats, salinity challenges may also vary among habitats. For example, haloclines in Hawaiian anchialine habitats can range from being absent (constant salinity throughout the habitat) to strong haloclines showing 15‰ at the surface and over 30‰ a few meters below (Havird et al., 2014c). Because each species was sampled from a single geographic habitat to minimize genetic differences within species (Craft et al., 2008; Weese et al., 2013), habitat-specific environmental differences may have been magnified. However, by definition, all anchialine habitats show variable salinities. Atyid species also generally export larvae to the open ocean, and larvae of at least H. rubra do not tolerate freshwater (Havird et al., 2015). Taken together, while habitat-specific fluctuations in salinity may be important, all anchialine organisms must likely tolerate at least some fluctuations in salinity.

Another hypothesis is that the species examined here have adapted to the varying salinities of the anchialine ecosystem via differing evolutionary trajectories. In other words, differing responses to salinity in anchialine shrimps may represent a case of divergent evolution. At a broad scale, C. rubella may rely on behavioral adaptation while the other species employ physiological mechanisms, characterized by different patterns of ion transport and salinity-induced gene expression in the gills. Although divergent evolution is commonly invoked (Gervasi and Schiestl, 2017; Rundle et al., 2005), providing evidence for a shared phenotype arising in response to a common environment via divergent mechanisms is relatively difficult. One example comes from experimental evolution in the bacterium Pseudomonas, where different mutations can be responsible for adaptation to the same environment (MacLean and Bell, 2003; Vogwill et al., 2016). More research is needed to further explore this hypothesis in anchialine shrimps, and intra-specific variation among populations within a species should also be considered.

Another hypothesis is that ancestral salinity may dictate the specifics of osmoregulation. Ancestrally freshwater taxa previously showed stable expression of ion transporters regardless of salinity (Faleiros et al., 2010; Huong et al., 2010; Maraschi et al., 2015; McNamara and Faria, 2012; McNamara and Lima, 1997; Ordiano et al., 2005; Tsai and Lin, 2007), while in those with a marine ancestry these genes were upregulated during salinity transfers (Havird et al., 2013). While H. rubra and M. minutus followed this pattern (at least somewhat), other atyids of freshwater ancestry examined here also showed upregulation of these genes (Fig. 4). Exploring more diverse taxa with transcriptome-wide approaches is needed to elucidate this and previous hypotheses.

In euryhaline marine crabs, where osmoregulation has been most heavily studied, salinity transfers trigger a suite of physiological changes. Salinity-mediated gene expression is characterized by the repression of osmoregulatory genes at high salinity via a yet-unidentified peptide hormone localized in the major endocrine complex, the X-organ/sinus gland complex (Henry, 2006; Henry and Campoverde, 2006; Mitchell and Henry, 2014b). At high salinity, expression is kept at minimal, baseline levels. In low salinity, repression is removed and gene expression is induced. This includes upregulation of several key ion transporters and accessory enzymes (Havird et al., 2013; Henry, 2001; Luquet et al., 2005), resulting in ions being pumped against their concentration gradients and leading to an internal osmolality higher than that of external freshwater.

While such upregulation is a hallmark of euryhalinity in most osmoregulating animals (Evans et al., 2005; Henry et al., 2012; Sakamoto et al., 2015), qPCR in H. rubra revealed stable expression of these genes regardless of salinity (Havird et al., 2014c), implying alternative molecular mechanisms of osmoregulation in at least one anchialine crustacean. Moreover, examination of salinity-induced gene expression at the whole-transcriptome level in classic crustacean models has identified novel salinity-responsive genes and pathways (Havird et al., 2016; Towle et al., 2011), suggesting unacknowledged nuances underlying the process of osmoregulation. Here, we detected novel patterns of gene expression that further implicate underexplored mechanisms of osmoregulation in crustaceans.

We previously hypothesized that ion transporters in H. rubra were expressed constitutively at high levels in order to anticipate the constant salinity fluctuations occurring within anchialine habitats of the Hawaiian Islands (Havird et al., 2014c). The RNA-Seq analyses reported here for H. rubra support this hypothesis: such transporters were not identified as differentially expressed during salinity transfer and possessed relatively constant expression levels regardless of treatment (Fig. 4; Fig. S5 of figshare dataset – see ‘Data availability’ section, below). Moreover, these ion transporters were among the top 5% of transcripts (1% for NKAα) based on average expression levels across all treatments. Preliminary proteomics studies also identify NKA as being among the most abundant proteins in the H. rubra gill (J.C.H. and S.R.S., unpublished data). In contrast, expression of these same transporters varied predictably with salinity in RNA-Seq analyses of the euryhaline crab C. sapidus (Havird et al., 2016) as well for H. trigonophthalma (NKCC) and C. rubella (NKA and NKCC) as shown here (Fig. 4; Fig. S5 of figshare dataset – see ‘Data availability’ section, below). This suggests that while some anchialine shrimps may utilize mechanisms more similar to that of the well-recognized example of euryhaline crabs, others have constitutively high expression of the same ion transporters.

If traditional ion transporters are not upregulated with salinity in H. rubra, what genes, if any, might respond to salinity? Interestingly, some salinity-responsive genes were ion transporters not traditionally implicated in osmoregulation. For example, transporters associated with calcium (Fig. 3; Fig. S4 of figshare dataset – see ‘Data availability’ section, below) have not been identified as differentially expressed in C. sapidus (Havird et al., 2016), but were found to be so in H. rubra and a few of the other species examined here. In general, calcium and/or its transport is not predicted to be explicitly involved in osmoregulation, although for largely unknown reasons, such transporters do tend to be upregulated during seawater transfer in rainbow trout (Flik et al., 1997; Marshall, 2002). One possibility is that calcium transporters are upregulated in low salinity in anticipation of a molt, which often occurs in low salinity (Wheatly, 1996). Alternatively, calcium transport during salinity transfer could be related to the altered expression of chitin-related genes potentially underlying general modifications in cellular structure (see below). Because calcium transport likely occurs via MRCs in the gills (Flik et al., 1997; Marshall, 2002; Mccormick et al., 1992), increased transport of this ion may be a byproduct of post-transfer MRC proliferation. Moreover, the major changes in cytoskeletal and ion transport processes along with the downregulation of catabolic processes might indicate an overall decline in physiological condition of the experimental animals. However, it is unclear whether these patterns extend beyond the gill and instead may just be indicative of the remodeling, or transient cell swelling, that takes place in gills during low-salinity acclimation.

Genes that play a role in light perception or vision were among the most differentially expressed transcripts during salinity transfer in H. rubra. Notably, these patterns were unique to H. rubra as expression of these transcripts was not detected in gills from the other species. Furthermore, our dissection method excluded eye tissue, suggesting the unusual expression patterns of these genes exclusive originate from the gills of H. rubra. Intriguingly, co-option of light/stimulus detection genes to other physiological functions has been previously documented, including thermosensor functionality in Drosophila (Shen et al., 2011) and known structural homologs of opsins acting as light-sensitive chloride ion pumps in bacteria (Yoshizawa et al., 2014). This raises the possibility that opsins and associated proteins may have underappreciated functions outside of light/stimulus detection by playing direct and/or indirect roles in salinity sensing or ion transport in the gills of H. rubra. Such a situation deserves further consideration and exploration.

In many of the examined species, ribosomal proteins were also identified as DEGs. Unfortunately, the relative timing between transcriptional and translational responses to salinity is unknown for all of the study species and most of the genes discussed here. However, estimates from crabs based on CA suggest relatively rapid induction of gene expression (i.e. within 6 h of salinity transfer) but slower induction of enzyme synthesis (i.e. by 48 h) (Henry et al., 2006; Mitchell and Henry, 2014a; Serrano et al., 2007). Given this, other salinity-induced genes may follow a similar pattern of post-transcriptional rate limitation. As such, altered transcription of ribosomal subunits might represent a response to increased translation pressures for critical salinity genes. It should also be noted that translational responses to salinity differed substantially between cellular compartments: mitochondrial ribosomal subunits were significantly downregulated during transfer for several of the species examined here, similar to the crab Carcinus maenas (Towle et al., 2011), while cytosolic ribosomal subunits were generally upregulated.

Additional DEGs related to the cytoskeleton were also common, including actins, myosins and various cuticle and chitin-related genes. This parallels previous findings in C. sapidus (Havird et al., 2016) and likely relates to the structural remodeling occurring in gills experienced by many taxa following salinity transfer (Dickson et al., 1991; Evans et al., 2005; Freire et al., 2008; Genovese et al., 2000; Lovett et al., 2006). Similarly, upregulation of OXPHOS subunits and energy-related genes has been documented in other crustaceans (Havird et al., 2016; Towle et al., 2011), likely reflecting the increased energetic demands of osmoregulation during transfer to dilute waters. Perplexingly, these genes underwent downregulation in some species examined here. Furthermore and unlike in crabs (Havird et al., 2016; Towle et al., 2011), general stress genes (e.g. HSPs) were recovered as upregulated DEGs in certain transfers/species. When stress-response genes were examined as a whole, however, most exhibited little response to salinity, generally supporting the hypothesis that fluctuating salinity is typical for anchialine crustaceans (Maciolek, 1986; Sket, 1996). Oxygen-transporting genes were also often identified as DEGs, and the intimate inter-relationships between respiration and ion transport in the gills may explain some of these findings (Ern et al., 2014; Evans et al., 2005; Havird et al., 2014a; Mangum, 1986). Finally, many of the DEGs identified here are fodder for additional hypothesis-driven research, including possible utilization of lipid energy stores during salinity transfer in H. rubra (Chen et al., 2014; Chen et al., 2019), hormonal regulators of ion transport (e.g. prostaglandins, atrial natriuretic peptides and an octopamine receptor; Towle et al., 2011), as well as the many other annotated and unannotated contigs not specifically discussed. These data are provided via figshare (https://figshare.com/articles/Disparate_responses_to_salinity_across_species_and_organizational_levels_in_anchialine_shrimps/9976883) to encourage further exploration.

We thank M. Hidaka and the Sesoko Tropical Biosphere Research Center, Okinawa, Japan for graciously hosting J.C.H. Dwi Haryanti provided assistance in the Hidaka Laboratory while in Okinawa. We thank J. Snelling and N. Kirk for assisting with RNA-Seq experiments. This represents contributions #190 and #92 to the Auburn University (AU) Marine Biology Program and Molette Biology Laboratory for Environmental and Climate Change Studies, respectively.