Genome-Wide Transcriptional Responses of Marine Nematode Litoditis marina to Hyposaline and Hypersaline Stresses

Maintenance of osmotic homeostasis is essential for all organisms, especially for marine animals in the ocean with 3% salinity or higher. However, the underlying molecular mechanisms that how marine animals adapt to high salinity environment compared to their terrestrial relatives, remain elusive. Here, we investigated marine animal’s genome-wide transcriptional responses to salinity stresses using an emerging marine nematode model Litoditis marina. We found that the transthyretin-like family genes were significantly increased in both hyposaline and hypersaline conditions, while multiple neurotransmitter receptor and ion transporter genes were down-regulated in both conditions, suggesting the existence of conserved strategies for response to stressful salinity environments in L. marina. Unsaturated fatty acids biosynthesis related genes, neuronal related tubulins and intraflagellar transport genes were specifically up-regulated in hyposaline treated worms. By contrast, cuticle related collagen genes were enriched and up-regulated for hypersaline response. Given a wide range of salinity tolerance of the marine nematodes, this study and further genetic analysis of key gene(s) of osmoregulation in L. marina will likely provide important insights into biological evolution and environmental adaptation mechanisms in nematodes and other invertebrate animals in general.


INTRODUCTION
Salinity, as an important ecological factor, affects the physiology and behavior of marine and terrestrial animals. As a nutrient element in the diet, salt is of vital importance to the health of animal and human. In human, chronic high dietary salt intake gradually causes an increased risk for cardiovascular disease, particularly hypertension; as well as other disease such as stroke, gastric cancer, kidney disease and obesity (Rust and Ekmekcioglu, 2017;He and MacGregor, 2018). Therefore, studies on the underlying mechanisms of animals' sensation, response, and adaptation to environmental salinity have always been a hot topic.
The multicellular model organism, Caenorhabditis elegans senses most of the environmental disturbance by the terminal cilia of sensory neurons (Bargmann, 2006). It is known that ASH neurons are required for the perception of high salt, while ASEs are for the low salt. The sensation of salinity stress can trigger subsequent avoidance behavior to protect worms from harmful salinity conditions (Hilliard et al., 2005;Kunitomo et al., 2013). Once the stressed salinity environment is unavoidable, C. elegans will engage a sleeplike quiescent behavior and cease locomotion and feeding, which is dependent on ALA neuron (Hill et al., 2014). Due to the imbalance between internal and external osmotic pressure, the body volume of nematodes undergoes significant changes under salinity stresses, manifests as shrinkage rapidly when environmental salinity increased, however, hypertonicityacclimated worms swell and then return to their initial body volume when exposed to low-salt condition (Lamitina et al., 2004). Organic osmolytes play an important role in osmotic regulation and salinity stress adaptation for all organisms. In C. elegans, cellular osmotic homeostasis can be maintained by rapid accumulation of glycerol upon high salinity challenge (Lamitina et al., 2004;Lamitina et al., 2006). It is well accepted that the C. elegans' cuticle might act as a "sensor" in responding to salinity stress damage, which in turn triggers downstream physiological changes (Choe, 2013;Dodd et al., 2018). On the other hand, numerous genes involved in osmotic regulation have been identified in C. elegans, such as the osmolyte glycerol synthesis enzyme gene (Lamitina et al., 2004(Lamitina et al., , 2006Choe, 2013), transient receptor potential cation channel TRP subfamily genes (Choe, 2013), chloride channel genes (Choe, 2013), aquaporin water channel genes (Igual Gil et al., 2017), extracellular matrix component genes (Lamitina et al., 2006;Choe, 2013), as well as genes related to MAPK, WNK-1/GCK-3, Notch and insulin-like signaling pathways (Choe, 2013;Dresen et al., 2015;Burton et al., 2017). Many of the above osmotic regulation genes play evolutionarily conserved roles in systemic osmotic homeostasis in yeast, flies, plants and mammals Burg et al., 2007;Brewster and Gustin, 2014;Pasantes-Morales, 2016;Zhou et al., 2016;Yang and Guo, 2018), providing clues for treatment of human disease that accompany osmotic perturbation.
C. elegans is one of the typical free-living terrestrial nematode species, whereas many extremophilic nematodes have been isolated in multiple extreme environments, such as Plectus murrayi, which is capable of tolerating Antarctic environmental extremes (Adhikari et al., 2009), Halicephalobus mephisto inhabiting a fluid-filled aquifer accessed from the Beatrix Gold Mine in South Africa at 1.3 km below the surface (Weinstein et al., 2019), and Auanema species in the arsenicrich, alkaline, and hypersaline Mono Lake (CA, US) (Shih et al., 2019). Multiple stress related genes have also bene unveiled in those extremophilic nematodes, including heat-shock protein (HSP) genes, late embryogenesis abundant (LEA) protein genes, trehalose-6-phosphate synthase (TPS) genes and others, ensuring their survival in a harsh and variable environment (Adhikari et al., 2009;Shih et al., 2019;Weinstein et al., 2019). Of note, about 43% of the known nematode species are distributed in the ocean (Appeltans et al., 2012;Zhang et al., 2015). It is speculated that nematodes may have emerged from a marine habitat during the Cambrian explosion (Van Megen et al., 2009), and colonized land about 442 million years ago (Rota-Stabelli et al., 2013). Salinity is obviously one of the most significant factors that changed during this successful terrestrialization. However, the underlying mechanisms are largely unexplored.
Litoditis marina is a dioecious free-living marine nematode, which is widely distributed in the littoral zone of coasts and estuaries, and plays an important role in these marine ecosystems (Derycke et al., 2016;Xie et al., 2020). It possesses some promising characteristics similar as C. elegans, such as short generation time, clear genetic background and a sequenced genome , which facilitated its laboratory application for the in-depth study of molecular biology, cell biology, physiology and behavior regulation in this species. Generally, the habitat salinity for intertidal marine nematodes, including L. marina, is frequently changed due to the influence of many factors such as tides, sun exposure, rainfall, ocean currents and climate. The effective sensation and response to the dynamic salinity environments is of great significance for marine nematodes' survival. However, the underlying molecular mechanism is still unknown.
In this study, we challenged L. marina L1 larvae with hyposaline and hypersaline stresses, respectively, and further demonstrated their genome-wide transcriptional signatures via RNA sequencing (RNA-seq) analysis. Both common and specific responding genes were identified in hyposaline and hypersaline stressed worms. These results not only provide a basis for understanding the salinity response mechanism for L. marina, but also might provide new clues for in-depth exploration of osmoregulation and environmental adaptation mechanisms for other marine animals.

Worms
The wild strain of marine nematode L. marina, HQ1, was isolated from intertidal sediments (Huiquan Bay, Qingdao). Healthy worms were cultured on SW-NGM agar plates (prepared with seawater with a salinity of 3%) seeded with Escherichia coli OP50 as a food source, as reported previously . Worms were maintained and propagated at 20 • C in the laboratory for about 3 years till this study.

Behavioral and Developmental Analysis Under Salinity Stresses
Artificial seawater-NGM agar plates were prepared by Sea Salt (Instant Ocean) in 0.3, 3, and 6% salinity, respectively. Two sets of salinity conditions, 0.3% (hyposaline) and 6% (hypersaline), were salinity stress treatment groups, while 3% was the control.
For behavioral and developmental analysis, 30 newly hatched L. marina L1 larvae were transferred onto each indicated 3 cmdimeter agar plates seeded with 15 µl OP50. Worms were scored as active if response was detected after prodding with a platinum wire 24 h post-treatment. The number of adult worms was scored 120 h (5 days) post-treatment.
Statistical analysis was performed using IBM SPSS Statistics 19. Data represent the average of at least three replicates unless specified otherwise. The comparisons between two groups were performed using the two-tailed Student's t-test. P < 0.05 was considered statistically significant. Particularly, we also treated C. elegans on 0.3 and 3% salinity plates as described above, 30 L1 larvae per plate in triplicates for each condition. We found that 100% L1 worms were active on 0.3% salinity plates, while none of tested L1 worms could survive on 3% salinity plates.

RNA-Seq Analysis
HQ1 strain worms cultured on SW-NGM plates were allowed to lay eggs overnight at 20 • C. Eggs were washed off and collected using filtered sterile seawater, then the eggs were treated with Worm Bleaching Solution (Sodium hypochlorite solution: 10 M NaOH: H 2 O = 4: 1: 10, prepared in terms of volume ratio) at room temperature for 1.5 min. Eggs were then washed twice with sterile seawater. The above clean eggs hatched overnight and developed into L1 larvae in sterile seawater at 20 • C. The synchronized L1 larvae were collected by filtration using 500 grid nylon filter with 25 µm mesh size, and then the L1 larvae were transferred to each 9 cm-dimeter agar plates prepared by Sea Salt mentioned above, which were seeded with 100 µl OP50 per plate (covering the entire plate evenly with a coating stick), respectively. Worms were collected after L1 larvae were incubated for 3 h at 20 • C under each salinity condition. Worms were washed with M9 for three times to remove the bulk of the residual bacteria. Excess supernatants were removed carefully via centrifugation. The samples were frozen immediately in liquid nitrogen. Total RNA was then extracted using Trizol (Invitrogen).
With three biological replicates for each treatment, a total of nine RNA libraries were prepared with 3 µg RNA using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, United States) following manufacturer's recommendations. Then, RNA libraries were sequenced on an Illumina NovaSeq 6000 platform and 150 bp paired-end reads were generated.
Firstly, clean data were obtained by removing reads containing sequencing adaptors, reads having poly-N and low-quality ones from raw data. The minimum of base score Q20 was over 97.5% and Q30 was over 93.13%. Then, the clean data were aligned to the L. marina reference genome  by Hisat2 (v2.0.5, with the default parameters) , with mapping ratio from 66.91 to 69.72% (Supplementary Table 1). New transcripts for novel genes were predicted and assembled by StringTie (v1.3.3b, with the default parameters) (Pertea et al., 2015), then annotated with Pfam, SUPERFAMILY, Gene Ontology (GO) and KEGG databases. Briefly, the functional annotation was performed using InterProScan (v5.17-56.0, Jones et al., 2014) by searching against publicly available databases Pfam 1 , SUPERFAMILY 2 , and GO 3 , with an E-value cutoff of 1e-5. KEGG function (Kanehisa and Goto, 2000) was assigned using KOBAS 3.0 (Xie et al., 2011) by best hit (with an E-value cutoff of 1e-5) to KEGG database 4 . Further, the reads numbers mapped to each gene were analyzed using featureCounts (v1.5.0-p3, with parameter -Q 10 -B -C) (Liao et al., 2014), and FPKM (expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced of each gene) was calculated based on the length of the gene and reads count mapped to this gene, which was used for estimating gene expression levels. The correlation coefficients (Pearson R) of 9 samples for three conditions were calculated using the cor.test() function (method = "pearson") in R 5 , based on the FPKM value for all genes. The high correlation coefficient within three replicates for each group (>0.92, Supplementary  Figure 1), indicated the reliable sample preparation. Hierarchical clustering for the sample can be visualized in Supplementary  Figure 2, which also supported the reliable sample preparation. Subsequently, differential expression analysis of two conditions was performed using the DESeq2 R package (v1.16.1) (Love et al., 2014). The resulting P-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value (padj) < 0.05 found by DESeq2 were assigned as differentially expressed. Hierarchical clustering for differentially expressed genes (DEGs) was performed using the pheatmap package in R 6 . Moreover, GO enrichment analysis and KEGG pathway enrichment analysis of DEGs were achieved by clusterProfiler R package (v3.4.4), an adjusted P-value (padj) < 0.05 was considered significantly enriched. GeneRatio was defined as the ratio of the number of differential genes annotated to the GO term or on the KEGG pathway to the total number of differential genes, respectively.
Synchronized L1 worms were separately treated under each salinity condition (0.3, 3, and 6%) using artificial seawater-NGM plates (prepared by Sea Salt, Instant Ocean) at 20 • C for 3 h. Each treatment was performed for three biological repeats. Total RNA was extracted using Trizol (Invitrogen), reserve transcribed to cDNA using the ReverTra Ace R qPCR RT Master Mix with gDNA Remover kit (TOYOBO, Code No. FSQ-301), and the cDNA was used for qPCR analysis using the QuantStudioTM 6 Flex Real-Time PCR System (Applied Biosystems) and SYBR Green detection system (TOYOBO,. The primers information of totally 15 salinity related genes, listed above, was shown in Supplementary File 1. Each experiment was performed in triplicates for each biological replication. Values were normalized against the reference gene EVM0013809, which is ortholog of C. elegans gene cdc-42 (Hoogewijs et al., 2008). No adults was observed under 6% salinity condition. Four replicates were performed for each experimental condition. (B) For behavioral analysis, 30 newly hatched L1s were transferred onto each indicated 3 cm-dimeter agar plates seeded with 15 µl OP50. Worms were scored as active if response was detected after prodding with a platinum wire 24 h post-treatment. Specifically, 100% of tested worms are active on the control plates (3%). Six replicates were performed for each experimental condition. Error bars represent the standard error of the mean from replicated experiments. Differences between groups were analyzed statistically employing the two-tailed Student's t-test using IBM SPSS Statistics 19 software. P < 0.05 was considered statistically significant. ***P < 0.001.
Gene expression was presented as a fold change using the delta delta Ct method (Livak and Schmittgen, 2001). The comparisons between the treatment and control groups were analyzed statistically employing the two-tailed Student's t-test using IBM SPSS Statistics 19 software, values were considered to be significant at P < 0.05.
Further, correlation analysis of the results of qPCR and RNA-seq for interest genes was performed using GraphPad Prism 5 software.

L. marina Behavioral and Developmental Defects Under Salinity Stresses
L. marina is maintained under 3% salinity condition in the laboratory, around 91% newly hatched L1 larvae developed into adulthood after 5 days at 20 • C ( Figure 1A). To test its salinity tolerance, we first treated L1 worms under two conditions: hyposaline with a 0.3% salinity and hypersaline with a 6% salinity. We observed that L1 worms were paralyzed immediately on both salinity plates, with obvious body volume change in a manner similar to that reported in C. elegans (Lamitina et al., 2004). Worms can recover their motility afterward. Compared to the control group (3% salinity), approximately 88% worms under hyposaline could move normally after 24 h, whereas only 29.4% worms could recover motility under hypersaline ( Figure 1B). Next, we did the same test applied to even higher salinity conditions such as 7 and 8%, and observed that worms cannot survive under those conditions, indicating 6% is probably the extreme high salinity for L. marina to tolerate.
Taken together, worms exhibited both significantly behavioral and developmental defects when stressed with either low salinity or high salinity. Notably, enhanced defects were demonstrated under high salinity.

RNA-Seq Analysis in L. marina Under Hyposaline and Hypersaline Environments
To investigate genome-wide responses in L. marina to salinity stress, we used RNA-seq analysis. Newly hatched L1s were treated for 3 h on low salinity (0.3%), normal salinity (3%, control) and high salinity (6%) plates, respectively (Figure 2A). A total of 1209 DEGs were identified under low salinity, and 1330 DEGs under high salinity (Figure 2A). Interestingly, there were 108 up-regulated DEGs and 93 down-regulated DEGs shared in both conditions (Figure 2A), indicating common response patterns under hyposaline and hypersaline stresses. On the other hand, condition-specific DEGs exhibited salinitydependent responsive and regulatory mechanisms in L. marina. Details of significantly up-regulated and down-regulated DEGs were listed in Supplementary File 2.

Shared Transcriptomic Signature of L. marina Under Both Low and High Salinity Stress Conditions
As both hyposaline and hypersaline stresses lead to behavioral and developmental defects in L. marina, we found common transcriptomic signature between these two conditions based on GO enrichment analysis.
As shown in Figure 2B, extracellular region related genes were significantly enriched in DEGs in both examined salinity Frontiers in Physiology | www.frontiersin.org FIGURE 2 | RNA-seq identifies the transcriptomic signature of short-time salinity stressed marine nematodes. (A) Experimental design of this study and the resulting transcriptomic signature of salinity stressed worms. Differentially expressed genes (DEGs, | log2foldchange| > 1; DESeq2 padj < 0.05) were determined for each condition. (B) GO enrichment analysis for DEGs. | log2foldchange| > 1; DESeq2 padj < 0.05 was set as the differential gene screening threshold. GO enrichment analysis of DEGs were achieved by clusterProfiler R package (v3.4.4), an adjusted P-value (padj) < 0.05 was considered significantly enriched. The color from red to purple represents the significance of the enrichment. GeneRatio was defined as the ratio of the number of differential genes annotated to the GO term to the total number of differential genes.
Frontiers in Physiology | www.frontiersin.org conditions. We found that a series of TTL family genes, such as EVM0003584/ttr-48 and EVM0004159/ttr-59, were up-regulated under both conditions (Figure 3A), indicating that extracellular region related genes could be induced by either low or high salinity stress. In addition, the TPS gene (EVM0007411/tps-2, Figure 3B), which is crucial for trehalose biosynthesis, was significantly induced upon both salinity stresses.
Overall, these shared features indicate the existence of conserved strategies for response to stressful salinity environments in L. marina.
However, three of the above UFAs biosynthesis related genes as well as all the above tubulin and IFT genes were significantly down-regulated in response to high salinity stress (Figures 4B,C), indicating their critical roles in salinity stress response.

Quantitative Real-Time PCR Validation
We applied qPCR to validate the expression patterns of interest genes identified from our RNA-seq results ( Figure 6A). Consistent trends were demonstrated in Figure 6B and Supplementary File 3. The expression (G) Expression level of ion transporter genes. Fold change indicates the ratio of FPKM value of the treatment group (0.3%, 6%, as indicated) to that of the control group (3%). The error bars represent standard error of the mean of three biological replicates per condition. The comparisons between the treatment and control groups were analyzed statistically using the Benjamini and Hochberg's methods calculated by DESeq2. The value of padj < 0.05 was considered statistically significant. *padj < 0.05, **padj < 0.01, ***padj < 0.001, ns-not significant. levels of TTL family gene EVM0003534, TPS gene EVM0007411/tps-2, were significantly increased under both 0.3 and 6% salinity conditions ( Figure 6A and  Supplementary File 3). While, the expression levels of dopamine receptor gene EVM0000190/dop-1, glutamate receptor gene EVM0013383/glc-4, acetylcholine receptor gene EVM0009741/eat-2, serotonin receptor gene EVM0012843/ser-1, neuropeptide Y receptor genes EVM0015448/npr-6 and FIGURE 4 | Genes that up-regulated in expression with decreasing salinity. (A) Hierarchical clustering of 144 DEGs that up-regulated in expression with decreasing salinity, performed using the pheatmap package in R. | log2foldchange| > 0; DESeq2 padj < 0.05 was set as the differential gene screening threshold. Relative expression levels (z-score, log-transformed gene expression level based on FPKM value) were indicated for each gene (row) in each sample (column). Red indicates upregulation; blue indicates downregulation. The scale bar shows the z-score for a differentially expressed gene. (B) Expression level of genes involved in biosynthesis of unsaturated fatty acids (UFAs). (C) Expression level of cytoskeleton tubulin and related intraflagellar transport (IFT) genes. Fold change indicates the ratio of FPKM value of the treatment group (0.3%, 6%, as indicated) to that of the control group (3%). The error bars represent standard error of the mean of three biological replicates per condition. The comparisons between the treatment and control groups were analyzed statistically using the Benjamini and Hochberg's methods calculated by DESeq2. The value of padj < 0.05 was considered statistically significant. *padj < 0.05, **padj < 0.01, ***padj < 0.001, ns-not significant.
EVM0010018/npr-4, ion transporter genes EVM0004010/kcc-2 and EVM0012374/twk-24, were significantly decreased under both conditions (Figure 6A and Supplementary  File 3). In addition, we confirmed that fatty acid elongase gene EVM0013022/elo-2, fatty acid desaturase gene EVM0001302/fat-4, tubulin gene EVM0007116/tba-5 were significantly up-regulated under hyposaline environment ( Figure 6A and Supplementary File 3). By contrast, the expression of cuticle collagen genes EVM0002243/col-156 and EVM0005554/col-107 were validated to be up-regulated under high salinity stress (Figure 6A and Supplementary File 3). DEGs that up-regulated in expression with increasing salinity, performed using the pheatmap package in R. | log2foldchange| > 0; DESeq2 padj < 0.05 was set as the differential gene screening threshold. Relative expression levels (z-score, log-transformed gene expression level based on FPKM value) were indicated for each gene (row) in each sample (column). Red indicates upregulation; blue indicates downregulation. The scale bar shows the z-score for a differentially expressed gene. (B) Expression level of cuticle collagen genes. Fold change indicates the ratio of FPKM value of the treatment group (0.3%, 6%, as indicated) to that of the control group (3%). The error bars represent standard error of the mean of three biological replicates per condition. The comparisons between the treatment and control groups were analyzed statistically using the Benjamini and Hochberg's methods calculated by DESeq2. The value of padj < 0.05 was considered statistically significant. *padj < 0.05, **padj < 0.01, ***padj < 0.001, ns-not significant.

Synchronization in L. marina for Large-Scaled Analysis
According to the synchronization methods used for hermaphroditic worms such as C. elegans and P. pacificus (Lewis and Fleming, 1995;Pires da Silva, 2005), we could not obtain enough L1 larvae samples for large-scale analysis in terms of the dioecious marine nematode L. marina. Therefore, we instead allowed L. marina adults to lay eggs on plates overnight, and then collected eggs from these plates. After a short-time treatment with bleaching solution to kill E. coli OP50 in the washing solution, eggs were subsequently incubated in filtered sterile seawater to hatch. Finally, synchronized newly hatched L. marina L1 larvae were effectively obtained by filtration using a grid nylon filter with mesh size of 25 µm. The establishment of large-scale L1 larvae collecting protocol developed in this study will facilitate further L. marina multi-omics studies, which requires large-scale synchronized worms.

Marine Nematode L. marina Has a Wider Range of Salinity Tolerance
In the intertidal areas, habitat salinity of L. marina is subject to either sudden or gradual changes in response to tides, rainfall, ocean currents, seawater evaporation and climate. In the laboratory, we found that L. marina could survive a wide range of salinity from 0.3 to 6%. Of note, for marine nematodes, 6% salinity is obviously an extreme condition, which is almost twice of that of sea water. By contrast, we noticed that C. elegans couldn't survive at 3% salinity (30 L1s per plate in triplicates). Thus, L. marina is a euryhaline marine nematode and has a wider range of salinity tolerance than its terrestrial relative C. elegans. Further studies using marine nematode L. marina as a model, will provide universal mechanisms underlying marine invertebrates' euryhaline adaptation.

Transthyretin-Like Family Genes Are Presumably Involved in the Damage Control Mechanisms in Response to Salinity Stresses
Under diverse environmental and physiological stresses, organisms usually demonstrate various degrees of cell damage by stress-induced protein misfolding, denaturation or aggregation, thereby disrupting proteostasis and cell homeostasis (Lamitina et al., 2006;Galluzzi et al., 2018). In the process of stress response, common stress-inducible genes, such as HSP genes, are induced to protect cells (Spees et al., 2002;Lamitina et al., 2006;Adhikari et al., 2009;Shih et al., 2019). Such stress-inducible genes were also found in our transcriptome results, for example, a FIGURE 6 | Validation of the RNA-seq results using qPCR. (A) qPCR analysis of interest genes identified from the RNA-seq results. Gene expression was presented as a fold change, indicating the ratio of the treatment group (0.3%, 6%, as indicated) to the control group (3%), which was calculated by the delta delta Ct method. The mean fold changes and standard error of the mean of three biological replicates are graphed. The comparisons between the treatment and control groups were analyzed statistically employing the two-tailed Student's t-test using IBM SPSS Statistics 19 software. P < 0.05 was considered statistically significant. *P < 0.05, **P < 0.01, ***P < 0.001, ns-not significant. (B) Correlation analysis of the results of RNA-seq and qPCR for interest genes. Each dot represents a gene, with detailed information shown in Supplementary File 3. Correlation analysis was performed and graphed using GraphPad Prism 5 software, Pearson R = 0.9423, with a P < 0.0001. series of HSP20, HSP70 family chaperone genes and dozens of proteasome related genes were significantly up-regulated under both low and high salinity stresses (Supplementary File 4).
In terms of shared common DEGs between both salinity stresses, one prominent type of significantly up-regulated were the TTL family genes. In L. marina, at least 38 TTL family genes have been annotated by database mining. In the present study, a total of 19 genes encoding TTLs were up-regulated under both hyposaline and hypersaline environments (Figure 3A), suggesting that they might play important roles in responding to salinity stresses in L. marina. TTLs represent one of the largest nematode-specific protein families, sharing sequence similarity to vertebrate transthyretins (Parkinson et al., 2004). In vertebrates, transthyretins are present in extracellular fluids to transport thyroid hormones as well as vitamin A (Vieira and Saraiva, 2014). In terms of nematode TTLs, they were presumed to participate in disposal of toxic lipophilic moieties and hormonal signaling (Parkinson et al., 2004;Jacob et al., 2007). However, these elusive genes have not been implicated in the response to salinity stress up to now and their functions are largely unknown. TTR-52 was reported as a bridging factor involved in cell corps engulfment and apoptosis (Wang et al., 2010;Mapes et al., 2012), indicating that up-regulation of TTL genes might be part of the damage control mechanisms in response to either low or high salinity stresses in L. marina.

Trehalose-6-Phosphate Synthase Gene Is Up-Regulated in L. marina Upon Both Low and High Salinity Stresses
Almost all organisms accumulate organic osmolytes when exposed to hyperosmolarity, and more than one type of osmolytes might be utilized for a particular organism (Burg and Ferraris, 2008). Unlike most marine invertebrates, which mainly use free amino acids and methylamines as organic osmolytes (Niu et al., 2020), it has well demonstrated that hyperosmotic stress in C. elegans activates rapid accumulation of organic osmolyte glycerol via the rapid up-regulation of the glycerol-3-phosphate dehydrogenate enzyme gene gpdh-1, which is a key gene for de novo glycerol synthesis (Lamitina et al., 2004(Lamitina et al., , 2006. Based on the annotation information of L. marina genome, the predicted glycerol-3-phosphate dehydrogenase gene, EVM0001663/gpdh-1, was significantly up-regulated under high salinity condition (Supplementary Figure 5), indicating that L. marina might utilize glycerol as an osmolyte in response to high salinity stress similar to C. elegans.
Trehalose, a disaccharide of glucose, is present in a wide variety of organisms including nematodes, and is known to act as stress protectant to against effects of dehydration, desiccation, heat, freezing as well as high osmotic stress (Wharton, 2003;Adhikari et al., 2009;Erkut et al., 2011;Hibshman et al., 2020). It not only supports survival by stabilizing lipid membranes and improving proteostasis during water loss, but also serves as an energy source. It is known that TPS gene encodes the enzyme catalyzing the first step of trehalose biosynthesis (Watts and Ristow, 2017). Previously, it was reported that trehalose accumulates in several anhydrobiotic nematodes in response to desiccation and/or low temperatures (Wharton, 2003;Hibshman et al., 2020). Moreover, trehalose levels were elevated in response to high salinity environment in C. elegans age-1 mutants, and RNAi knockdown of tps-1 and tps-2 reduced trehalose levels by 90% and the hypertonic stress resistance of age-1 mutant was dramatically decreased, indicating an important functional role of TPS in hypertonic stress resistance in C. elegans (Lamitina and Strange, 2005). In this study, we found that L. marina TPS gene, EVM0007411/tps-2, was significantly up-regulated under both hyposaline and hypersaline stresses (Figure 3B), which could possibly cause accumulation of trehalose to facilitate its adaptation to both stress environments.

Certain Neuronal Signaling Are Transcriptionally Repressed by Salinity Stresses
In C. elegans, acetylcholine, serotonin, dopamine, glutamate, and neuropeptide are known important neurotransmitters, which play essential roles in a broad repertoire of behaviors, including locomotion, feeding, reproduction, social behavior, mechanosensation, chemosensation, learning, memory, behavioral plasticity, and adaptation (Chase and Koelle, 2007;Rand, 2007;Hukema et al., 2008;Frooninckx et al., 2012;Vidal-Gadea and Pierce-Shimomura, 2012). Under harsh conditions, the nervous system plays critical roles in worm stress response, facilitating worm survival and adaptation (Kim and Jin, 2015). Although, specific neurons that function in osmotic stress response have been studied extensively, for example, it is well known that ASH neurons are required for hyperosmotic sensation and ASEs are involved in hypoosmolarity response (Hilliard et al., 2005;Kunitomo et al., 2013). When osmotic stress is unavoidable, worms engage a sleep-like quiescent behavior and cease locomotion and feeding, which is dependent on ALA neuron (Hill et al., 2014). However, specific neurotransmitters, including serotonin, dopamine, glutamate and neuropeptide are only reported affecting avoidance behavior to NaCl in C. elegans (Hukema et al., 2008;Watteyne et al., 2020). To our knowledge, this is the first to describe a systematic repressed cholinergic, serotonergic, dopaminergic, glutamatergic and neuropeptide signaling in response to both low and high salinity stresses in L. marina.
In the present study, these reduced neural signaling are positively correlated with the reduced worm mobility in salinity stress environments. Furthermore, we speculate that these reduced signaling might further activate osmoregulation pathways in L. marina to promote its adaptation to the stressed salinity conditions. Unsaturated Fatty Acids Are Involved in Hyposaline Stress Response in L. marina In our study, genes involved in UFAs biosynthesis were significantly up-regulated under low salinity condition. Defects in UFAs biosynthesis have been reported to cause deficiencies in worm growth, development and neurological function (Watts and Browse, 2002;Kniazeva et al., 2003). Notably, we found that several genes encoding fatty acid elongase (EVM0013022/elo-2) and fatty acid desaturases (EVM0008235/fat-2 and EVM0011847/fat-3) were significantly decreased under hypersaline stress (Figure 4B), which might account for the developmental defects in 6% salinity stressed worms.
UFAs have a profound effect on the fluidity, flexibility and permeability of cell membranes, as well as play important roles in energy storage and signaling process (Zhu and Han, 2014). The synthesis of UFAs is regulated during changing environmental conditions, and the UFAs play a crucial role in environmental adaptation in a broad range of organisms from microorganisms to animals and plants (Wada et al., 1990;Miquel et al., 1993;Svensk et al., 2013;Yancey, 2020;Zhang Z. et al., 2020). It has been widely reported that there were more UFAs in salttolerant plants, yeast as well as bacterial cells (de Carvalho et al., 2014;Guo et al., 2019;Li et al., 2019). Moreover, Lucu et al. (2008) reported that UFAs might be important during the acclimation of the shore crab Carcinus aestuarii to hypoosmotic condition. However, nematode UFAs have not been implicated in the response to salinity stress as far as we know. In the present study, body swelling is initially observed for L. marina L1 larvae under hyposaline stress condition, while worms present a better ability to tolerate and acclimate to this condition, suggesting that the regulatory mechanism of body volume recovery is robust under hyposaline stress. This may be related to the cell membrane properties and UFAs compositions in L. marina. Thus, we propose that the effective induction of UFAs biosynthesis genes might act as part of the protective and adaptative strategies of marine nematodes upon low salinity stress.
Cuticle Collagen Genes Are Involved in Hypersaline Stress Response in L. marina Collagens are major structural proteins for the nematode exoskeleton, the cuticle. The identified C. elegans cuticle collagen mutants usually show either body morphology or locomotion defects, which is consistent with the cuticle's essential role for maintenance of body shape, as well as movement via attachments to muscles (Page and Johnstone, 2007). Moreover, several cuticle collagen mutants, such as dpy-2, dpy-7, and dpy-10, were reported to exhibit constitutive activation of gpdh-1 expression and glycerol accumulation, and show osmotic resistance phenotype (Lamitina et al., 2006;Wheeler and Thomas, 2006;Dodd et al., 2018). It is believed that the cuticle may serve as a sensor for osmotic stress and play important roles in C. elegans osmotic regulation (Lamitina et al., 2006;Dodd et al., 2018).
The cuticle in C. elegans is synthesized and secreted by underlying hypodermis, this process occurs five times during development, first at the end of embryogenesis before hatching and then again at the end of each larval stage before molting. C. elegans cuticle collagens are encoded by a multi-gene family, consisting over 170 genes (Page and Johnstone, 2007). In fact, these genes are not all expressed at the same time during cuticle synthesis. Clear spatial and temporal differences can be observed for individual genes (Johnstone, 2000;Page and Johnstone, 2007), indicating each cuticle collagen gene has specific roles in cuticle formation and function. The cuticle functions as the primary barrier between worm and its environment, and acts as the first line of defense against environmental stresses (Page and Johnstone, 2007). Previously, Dodd et al. (2018) has shown that multiple cuticle collagen genes can be induced when C. elegans exposed to high NaCl. Recently, we reported that a battery of collagen genes in C. elegans increased their expression to deal with acidic pH stress environments (Cong et al., 2020). Together, these results indicated a protective role of collagens in response to various stresses. Similarly, in the present study, an abundant group of over 20 cuticle collagen genes were significantly up-regulated upon hypersaline stress (Figure 5B), while almost half of them were remarkably down-regulated upon hyposaline stress, indicating certain collagens could specifically function to detect or transform salinity stress-induced signals, or just change the chemical and physical composition of the cuticle to provide the primary barrier to defend the dynamic osmotic variation.

Cytoskeleton Related Genes Are Differentially Regulated Under Different Salinity Stresses in L. marina
Tubulin is the basic component of cytoskeleton microtubules. It plays an indispensable role in structure maintaining, neuronal sensation as well as many other cell processes including intraflagellar transport. In our results, four tubulin genes, such as EVM0017317/tba-4, EVM0007116/tba-5, EVM0015250, and EVM0004244/ben-1, exhibited specific up-regulation when salinity is decreasing (Figure 4C). The tba-4 gene encodes αtubulin in C. elegans, both TBA-5 and BEN-1 are neuronal tubulins (Hurd, 2018). TBA-5 is an axonemal α-tubulin expressed in amphid and phasmid sensory neurons where it is localized to cilia. BEN-1 is a neuronal β-tubulin in C. elegans, which is also broadly expressed in the nervous system. Moreover, it is reported that IFT is essential for assembly, maintenance and function of sensory cilia in C. elegans (Hao et al., 2011). Here, we found that two IFT genes (EVM0001368/osm-3 and EVM0002085/daf-10) were induced specifically under low salinity condition ( Figure 4C). The osm-3 gene in C. elegans encodes a kinesin-2 family member of IFT motors, mediating IFT particles transport within sensory cilia (Prevo et al., 2015). daf-10 is required for IFT and for proper development of a number of sensory neurons (Bell et al., 2006). By contrast, both tba-5 and daf-10 genes were found to be upregulated in the osmo-resistant dpy-7 worms (Dodd et al., 2018). Moreover, the above tubulin and IFT genes exhibited significantly opposing changes between low and high salinity stresses (Figure 4C), indicating their critical roles in salinity stress response.
Recently, we found strong positive selections on EVM 0011198.1/col-109, EVM0009368.1/lec-1, EVM0011012.1/T13H 5.6, EVM0014994.2/osm-12, EVM0011797.1/gem-1, and EVM00 08924.1/imp-2 genes in L. marina genome , which might indicate an adaptation to the high salinity environment in the sea. Thus, we checked the expression level of these putative salinity related positive selected genes (PSGs) in our transcriptome data. We found that four PSGs were significantly regulated in response to salinity stresses (Supplementary Figure 6). The underlying mechanisms that how those PSGs function in osmoregulation in L. marina deserve to be further explored.

CONCLUSION
In conclusion, we have described for the first time the genome-wide transcriptional responses to both hyposaline and hypersaline stresses in the marine nematode L. marina. The present study will provide an essential foundation for identifying the key genes and genetic pathways required for osmoregulation in the marine nematodes. Given a wide range of salinity tolerance of the marine nematodes, our results and further genetic analysis of key gene(s) of osmoregulation in L. marina will likely provide important insights into biological evolution and physiological adaptation mechanisms in nematodes and other organisms in general.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI; PRJNA694479.

AUTHOR CONTRIBUTIONS
YX and LZ conceived and designed the experiments. YX carried out most of the experiments, analyzed the data, and wrote the manuscript. PZ contributed to the RNA-seq sampling and qPCR validation. LZ edited the manuscript and supervised the project. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
We are grateful to all members of the LZ laboratory for their helpful discussions.