Tolerance at the genetic level of the brine shrimp Artemia salina to a wide range of salinity

Background: The brine shrimp Artemia salina can thrive in a variety of salinities and is commonly distributed in natural hypersaline lakes and solar salterns. The zooplankter A. salina proves to be a lter feeder, consuming the alga Dunaliella and prokaryotes and plays a critical role in the hypersaline food web. However, the high salinity adaptation mechanisms of A. salina remain poorly understood through transcriptome analysis. Here, we examined the gene expression patterns of A. salina adults that were salt-adapted for 2–4 weeks at ve salinities (35, 50, 100, 150, and 230 psu), and generated long-read isoform sequencing (IsoSeq) data to construct a high-quality transcriptome assembly of A. salina. The patterns of A. salina along the salinity gradient provide evidence for halotolerant and euryhaline adaptations at the genetic level. Results: We conrmed that the activity of sodium/potassium ATPase was up-regulated at the genetic level in high salinity waters. Interestingly, genes related to beta-mannosidase and mannose activities were also up-regulated, suggesting that mannose and mannose derivatives may be accumulated as organic osmolytes. Alternatively, considering that glucose and galactose-related activities were suppressed at high salinities, mannose may be the primary sugar involved in the glycolytic pathway under such conditions. This result further supports the theory that mannose is the main energy source used by A. salina in highly saline environments. The gene expression patterns of A. salina may also be affected by increased thickness of the cuticle, increased numbers of mitochondria, and low dissolved oxygen in high salinity waters. Furthermore, the cellular response of A. salina to acclimation to intermediate salinities depends on the number and type of genes expressed; differential expression patterns are likely to uctuate at the population level. Conclusions: Our results provide a high-quality transcriptome assembly of the cosmopolitan brine shrimp Artemia salina at ve different salinities (35, 50, 100, 150, and 230 psu) for the rst time. The gene expression patterns of salt-adapted A. salina display

The brine shrimp A. salina (Branchiopoda; Pancrustacea; Crustacea) is often observed in natural hypersaline lakes (e.g., the Great Salt Lake) as well as in arti cial solar salterns [18,19]. This zooplankter can survive in saline waters ranging from 5 psu (practical salinity unit) to > 300 psu, demonstrating that it is a halotolerant and euryhaline species [20][21][22]. A. salina is a lter feeder, consuming phytoplankton (e.g., Dunaliella) and prokaryotes in these high salinity waters. Thus, this extremophile species must be a strict osmoregulator, as it regularly ingests high salinity water along with prey [23][24][25][26]. Several studies report that an increase in salinity stimulates Na + /K + -ATPase activity in multiple organs (e.g., metepipodites, maxillary glands, and gut) of adult A. salina for the excretion of accumulated salts [23][24][25][26]. In addition, the hemolymph of A. salina is maintained at a constant ~ 300 mOsm/kg (9 psu salinity), which is similar to the osmolality of teleost shes, when exposed to a range of external high salinities [25,27]. However, the differential gene expression patterns in adult A. salina along salinity gradients have not been investigated [22]. Moreover, other adaptive osmoregulatory mechanisms of A. salina to differing salinity ranges remain poorly understood.
Here, we generated a de novo transcriptome assembly of adult A. salina using a hybrid assembly approach, with long-read isoform sequencing and high-throughput RNA sequencing data, to analyze the differential gene expression patterns in a variety of salinities (i.e., 35,50,100,150, and 230 psu). In the present study, A. salina signi cantly up-regulated genes (e.g., Na + /K + -ATPase gene) related to osmoregulation with increasing salinity from 35 psu to 230 psu. In addition, we note that low oxygen concentration and thickness of A. salina's body surface may greatly affect its metabolic processes at high salinities.

Methods
Sample preparation, isoform, and RNA sequencing of Artemia salina Artemia salina was isolated from 140 psu saline water collected from a solar saltern in Uiseong (36°60′20.86″N, 126°29′71.16″E), Republic of Korea, in April 2018 ( Fig. 1; an additional movie le shows this in more detail [see Additional File 1]). A population of male and female A. salina in a 20-L aquarium, also containing sediment from the saltern and 15 L arti cial saline water made by dilution of Medium V (300 psu) [28], was maintained at 150 psu. Approximately 20 adults with a half-and-half sex ratio were transferred and maintained in arti cial saline water at 35, 50, 100, 150, and 230 psu for 2-4 weeks. After this time, we collected 10 adults ( ve males and ve females) from each salinity medium to perform transcriptome-level analyses of isoforms. Total RNA was extracted from A. salina using the RNeasy Plus Mini Kit (Qiagen, Hilden, Germany). The RNA sequencing libraries were constructed using the Illumina Truseq Standard mRNA Prep Kit (Illumina Inc., San Diego, CA, USA), and the qualities of the libraries were checked using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The RNA sequencing libraries were sequenced using 100-bp paired-end reagents with an Illumina Novaseq6000 (Table S1 in Additional File 2). The isoform sequencing (IsoSeq) library was constructed using a SMARTer PCR cDNA Synthesis Kit and DNA Template Prep Kit 1.0; the IsoSeq data were generated (1.0 Gbp; 780,667 reads) using the PacBio Sequel platform (Paci c Biosciences, Menlo Park, CA, USA). The high-quality (HQ) isoform consensus (16.6 Mbp; 16,555 reads) of the IsoSeq data was constructed using a SMRT Link 5.1.0 with the Iso-Seq2 application platform. Raw reads of RNA sequencing and high-quality Isoform consensus data were uploaded to the NCBI SRA (Sequence Read Archive) database with accession numbers (SRR12358297 -SRR12358302; BioProject PRJNA649341).
De novo transcriptome assembly, protein predictions, and functional annotations To construct high-quality isoform sequences, we generated a hybrid de novo transcriptome assembly of Artemia salina using a long-read HQ isoform consensus (16.6 Mbp; 16,555 reads; PacBio Sequel) and high-throughput RNA sequencing (a total of 34.6 Gbp; 100 × 100 bp Illumina Novaseq6000 reads). To improve sequence accuracy, the error correction step of the IsoSeq consensus was performed using the high-accuracy Illumina sequencing reads, which were aligned by Bowtie2, and variant calling in the alignment was processed using the Samtools program ( Figure S1 in Additional File 3) [29,30]. The corrected IsoSeq reads were assembled with the short-read Illumina sequencing data using the Trinity assembler (v2.8.4; default option of de novo transcriptome assembly with '--long_reads') [31,32]. To reduce redundancy, the assembled transcripts were clustered using the cd-hit-EST program (v4.6; -c 0.95 -aS 0.95 -n 10; https://github.com/weizhongli/cdhit). The RNA-seq raw reads were mapped into the assembled transcripts using the Salmon program (default options) [33], and the potential false-positive transcripts including zero-tpm (transcripts per million) values were excluded. Potential ribosomal RNA fragments in the assembled transcripts were excluded based on the results of a local BLASTn search (evalue cutoff = 1.e-10) using Artemia salina ribosomal RNA sequences (X01723.1 and AF169697.1).
All possible open reading frames (ORFs) in the assembled transcriptome data were predicted by a 6frame translation by Python (minimum length of protein-coding sequences = 150 bp), so that all three forward and reverse transcript frames were translated with the standard code (NCBI genetic code table 1).

Analyses of conserved eukaryotic gene sets and phylogeny
Conserved eukaryotic gene set completeness analyses were conducted using BUSCO (Benchmarking Universal Single-Copy Orthologs; v3.0.2) [47,48]. Sixteen taxa from Pancrustacea, including A. salina, were used to reconstruct the molecular phylogenetic tree. The phylogenetic position of A. salina was inferred using 5,604 concatenated protein sequences, which were related to Pancrustacea (local BLASTp e-value cutoff = 1.e-05). All protein sequences were aligned using MAFFT v7.313 (default options: --auto) [49]. Each alignment was trimmed when an aligned locus included more than 70% gap sequences. The phylogenetic analysis was inferred using maximum likelihood (ML) analysis. The ML tree was estimated using IQ-tree v.1.6.12 [50] with the best-t evolutionary model, and statistical support was estimated using bootstrapping with 1,000 replicates.

Differentially expressed gene analysis of Artemia salina
Mapping of RNA-seq raw reads was performed using Salmon (default option) [33]  Based on the z-score dynamics of each gene candidate, we analyzed up-and down-regulated gene expression patterns using customized Python scripts, and sorted the genes by salinity-dependent gene expression patterns including only those with following criteria: maximum fold-change ≥ 2-fold between the highest and the lowest tpm value, and the lowest tpm ≥ 5 in a target gene. Gene ontology (GO) terms were extracted from the functional annotations of the target genes described by the eggNOG-mapper results, and the enrichment of GO terms was performed using the topGO package in R (Fisher test, p-value < 0.05).

Results
De novo transcriptome assembly and phylogenomic analysis A total of 18,301 protein-coding sequences from Artemia salina were obtained through ltration of the assembled transcripts (i.e., removing too-short, redundant, and false-positive transcripts). Reference data were used for the homologous pancrustacean sequences (Figure S1 in Additional File 3; see above). Based on the BUSCO analysis, the protein pro les of A. salina included 92.3% conserved eukaryotic gene sets, a reliability quality comparable to those of other crustaceans (85-94% of BUSCOs; Table S2 in Additional File 2). The class Branchiopoda, which includes A. salina, formed a robust clade, and A. salina was closest to A. franciscana with maximal bootstrap support (ML: 100%; Fig. 2). Branchiopoda was the sister group to the class Hexapoda (ML: 100%), which also formed a monophyletic group (Fig. 2).
Differentially expressed gene patterns along the salinity gradient Due to the wide salinity tolerance range of A. salina, we evaluated variations in the gene expression patterns of salt-adapted cells along a salinity gradient (35 psu to 230 psu) based on z-scores using custom Python scripts (Python v2.7.16; Fig. 3). The enriched GO terms of the salinity-dependent genes (topGO package in R, Fisher test, p-value < 0.05) revealed that A. salina signi cantly up-regulated genes related to Na + /K + -ATPases, mannose and carotenoid metabolism, and monocarboxylate transmembrane transporters with an increase in salinity from 35 psu to 230 psu ( Fig. 3; Table S3 and S4 in Additional File 2; Figure S2 in Additional File 3). Conversely, A. salina signi cantly down-regulated genes associated with galactose, and glucose metabolism, N-acetylgalactosamine-4-sulfatases, amino acid symporters, and amino acid transmembrane transporters in a salinity-dependent trend ( Fig. 3; Table S3, and S4 in Additional File 2; Figure S2 in Additional File 3). Interestingly, the frequency of acidic, basic, hydrophilic, and hydrophobic amino acids in A. salina was identical to that of other pancrustacean species (Table S5 in Additional File 2). This result suggests that the osmoregulation system of the extremely euryhaline A. salina may operate by excreting excess salts instead of involving conformational changes to the proteomes within the animal's cells.
In addition to the those involving osmoregulatory processes of A. salina (e.g., Na + /K + -ATPase), transcripts related to the visual cycle and mitochondrial morphogenesis were also up-regulated in hypertonic waters ( Fig. 3; Table S4 in Additional File 2). Furthermore, genes involved in oxidative stress regulation and glycoprotein catabolic processes were up-regulated at high salinities ( Fig. 3; Table S4 in Additional File 2).
In contrast, the transcript encoding N-acetylgalactosamine-4-sulfatase, which hydrolyzes sulfates, was down-regulated in high salinity waters ( Fig. 3; Table S4 in Additional File 2). Furthermore, A. salina noticeably down-regulated genes related to diverse transporter activities in response to high salinity ( Fig. 3; Table S4 in Additional File 2). This indicates that the transporter systems in adult A. salina may be repressed in highly saline environments. U-shaped or inverted U-shaped expression patterns Artemia spp. grow best in a salinity range of 100-150 psu and can tolerate < 5 psu or > 300 psu [20-22, 51, 52]. Therefore, some genes may be expressed in a quadratic pro le, such as a U-shaped (i.e., upregulated at the two salinity extremes) or inverted U-shaped (i.e., down-regulated at the two salinity extremes) curve. A total of 326 A. salina genes were expressed in U-shaped curves with the lowest point at 50 psu; 222 and 411 genes were expressed in U-shaped curves with the lowest points at 100 and 150 psu, respectively ( Fig. 4; Table S6 in Additional File 2). Furthermore, inverted U-shaped expression patterns were observed in 379 genes with the highest point at 50 psu; 289 and 63 genes were detected with the highest points at 100 and 150 psu, respectively ( Fig. 4; Table S6 in Additional File 2).

Discussion
Halophilic and halotolerant organisms are capable of adapting to high salinity waters and maintaining their internal osmotic balance in order to thrive in harsh environments. In general, the salt-in and salt-out processes are regarded as two distinct osmoregulation strategies. Salt-in organisms typically utilize potassium as the main intracellular cation in high salinity waters, while salt-out organisms accumulate organic solutes (e.g., betaines, glycerol, ectoine, sucrose, and mannitol) as the main osmolytes [17,59]. In addition, salt-in organisms characteristically contain acidic proteins with a negative charge, which are not observed in the salt-out organisms, to avoid protein aggregation [17,60]. In the present study, Artemia (class Branchiopoda) may be considered a salt-in organism due to increased Na + /K + -ATPase activity in high salinity water (230 psu). However, the notable absence of acidic protein pro les in the brine shrimp Artemia was identical to other crustaceans in seawater. Thus, our transcriptomic analysis con rms that A. salina is a strict hyporegulator [24][25][26]61].
In crustaceans, the enzyme Na + /K + -ATPase serves to maintain sodium and potassium ion homeostasis across cell membranes [62][63][64][65]. Several studies have used the silver staining method to demonstrate that the metepipodites, digestive gut, and maxillary glands in A. salina exhibit high Na + /K + -ATPase activity under highly saline conditions [24,26,61]. This con rms that the enzyme Na + /K + -ATPase plays an essential role in A. salina osmoregulation at the genetic level. Interestingly, genes related to betamannosidase and other enzymes involved with mannose metabolism were up-regulated in high salinity media. This result implies that A. salina may accumulate mannose or mannose derivatives as organic osmolytes. Mannose and mannitol are widely recognized as organic osmolytes in diatoms, brown algae, green algae, and terrestrial plants [66,67]. However, most marine animals accumulate betaine, taurine, trimethylamine oxide, glycine, alanine, proline, homarine, or arginine as organic osmolytes [68]. Thus, the accumulation of mannose and mannose derivatives at high salinities in A. salina is particularly interesting. Alternatively, considering that glucose and galactose-related activities were suppressed at high salinities, mannose may be the primary sugar involved in the glycolytic pathway under such conditions. Furthermore, A. salina down-regulated UDP-glucose:hexose-1-phosphate uridylyltransferase (synonym: galactose-l-phosphate uridylyltransferase), which converts galactose-1-phosphate into glucose-1-phosphate, at high salinities. This result further supports the theory that mannose is the main energy source used by A. salina in highly saline environments. Bunn & Higgins [69] reported that organisms that accumulate high concentrations of aldohexoses with unstable ring structures (i.e., mannose and galactose) may experience enzymatic malfunctions as a result of covalent modi cations to individual proteins. Therefore, the covalent modi cation of A. salina proteins may also be detected in high salinity waters. Furthermore, Horst [70] noted that mannose and mannose derivatives in A. salina are substantially involved in glycoprotein catabolic processes, suggesting that exposure to increased salinity may result in an increased activity of mannose-related processes. Investigation of the primary role of mannose and mannose derivatives in A. salina under various salinities is warranted.
While Artemia can maintain ion homeostasis in its cells as the external salinity increases, more energy is required for growth and reproduction [71]. Theoretically, the dissolved oxygen (DO) concentration in 35 psu media at 25 °C is 3.2 times lower than that in 220 psu media [72]. In this study, we present several pieces of evidence that suggest that low oxygen may substantially impact the growth and reproduction of A. salina when salinity increases from 35 psu to 230 psu. Additionally, marine invertebrate vision is one of the most energetically demanding functions and is highly susceptible to dissolved oxygen uctuations [73]. In high-salinity waters, transcripts related to retinal metabolic processes were signi cantly up-regulated; A. salina signi cantly up-regulated ninaB homologous genes in high-salinity media. Many diverse metazoan species contain ninaB homolog genes ( Figure S3 in Additional File 3), which are associated with the synthesis of visual pigment and oxidative stress [74][75][76]. Moreover, the upregulation of genes related to oxidative stress suggests that A. salina can adapt to waters simultaneously high in salinity and low in oxygen. The activity of N-acetylgalactosamine-4-sulfatase, which was one of the down-regulated genes in A. salina, can decline when exposed to low oxygen [77]. Therefore, it is likely that the gene expression patterns of A. salina are greatly affected by low dissolved oxygen levels in its media.
Few mitochondria and a thin body surface cuticle are characteristic of Artemia in low salinity waters, whereas many mitochondria and a thick cuticle have been observed in high salinity waters [23,27,54]. In this study, transcripts related to mitochondrion morphogenesis and glycoprotein catabolic processes were up-regulated at high salinities, which is consistent with previous investigations. It is reasonable to infer that the rate of mitochondrion morphogenesis is inversely related to the dissolved oxygen concentration in saline waters. Meanwhile, the increased thickness of the A. salina cuticle associated with glycoprotein catabolic processes lowers the water permeability of the body surface, resulting in a relatively low rate of oxygen diffusion [23,54]. Thus, A. salina may require more mitochondria in high salinity than low salinity environments. Most sugar (i.e., glucose and galactose) transport systems, as well as those for nucleotides and amino acids, were repressed in A. salina at high salinities. These results imply that Artemia may be unable to properly use these essential macromolecules for growth and reproduction under highly saline conditions. Even though Artemia spp. are extremely halotolerant and euryhaline, signi cant energy expenditures are likely required to accommodate these adaptations.
However, the increased thickness of A. salina body surface layers at high salinities indicates that the reduction in permeability to water, oxygen, and essential nutrients occurs passively. Moreover, several studies have reported decreased protein, carbohydrate, and glycogen contents in some crustacean species with increasing salinity [78][79][80].
The cellular response of A. salina required for acclimation to intermediate salinities depends on the number and type of genes expressed. A U-shaped gene expression pattern observed along the salinity gradient in A. salina implies that it can adapt well to intermediate salinities [81,82], whereas an inverted U-shaped pattern indicates the presence of salinity response at intermediate salinities [83]. At 150 psu, 411 genes were expressed at low points in U-shaped patterns and only 63 genes were expressed at high points in inverted U-shaped patterns, implying that A. salina is well-adapted to 150 psu; most genes were actively expressed to tolerate this salinity. Furthermore, the U-shaped expression patterns observed in A. salina represent cellular responses related to cellular signaling, catabolic processes, morphogenesis, and development. Meanwhile, the inverted U-shaped expression patterns were related to diverse environmental response factors (e.g., salt aversion, sensory perception of salty taste, cellular response to light intensity, UV-B, ozone, response to oxidative stress, and starvation). Based on enriched GO terms and KEGG metabolic pathways, the functional annotations of the U-shaped expression pattern genes were usually different from those of the inverted U-shaped expression pattern genes ( Fig. 4; Table S6 in Additional File 2; Figure S4 in Additional File 3). This result indicates that the genes related to U-shaped and inverted U-shaped patterns are unequally expressed, and the expression pattern observed at a speci c salinity is related to the physiological characteristics of A. salina.
Within the genus Artemia, the gene expression patterns we observed with increasing salinity were not always consistent with previously identi ed patterns ( Fig. 5; Table S7 in Additional File 2). Intriguingly, the types of genes displaying U-shaped and inverted U-shaped expression patterns were usually differentiated at each intermediate salinity (i.e., 50, 100, and 150 psu). These results demonstrate that A. salina could be a potential model organism to study locally adapted populations called 'local adaptation' [19,84,85], and differential expression patterns of Artemia are likely to uctuate at the population level. Campillo et al. [85] reported that rotifer populations (Brachionus plicatilis) might present different cellular responses depending on the salinity of their medium. Thus, the gene expression of A. salina could also uctuate depending on the occurrence of ecological specialization in each population at speci c salinities. Further study is needed to con rm that A. salina can provide a model for local adaptation.

Conclusions
The gene expression patterns of salt-adapted A. salina at ve different salinities (35,50,100,150, and 230 psu) displayed greater osmoregulation process complexity than previously thought. A. salina along the salinity gradient could differentially express the genes related to the environmental variables (e.g., high salinity and low oxygen) and morphological transformations (e.g., the thickness of the cuticle and the numbers of mitochondria). Moreover, the types of genes displaying U-shaped and inverted U-shaped expression patterns suggested that A. salina appeared to have substantially differential adaptive osmoregulatory mechanisms to intermediate salinities.