Transcriptomic and Proteomic Analysis of Marine Nematode Litoditis marina Acclimated to Different Salinities

Salinity is a critical abiotic factor for all living organisms. The ability to adapt to different salinity environments determines an organism’s survival and ecological niches. Litoditis marina is a euryhaline marine nematode widely distributed in coastal ecosystems all over the world, although numerous genes involved in its salinity response have been reported, the adaptive mechanisms underlying its euryhalinity remain unexplored. Here, we utilized worms which have been acclimated to either low-salinity or high-salinity conditions and evaluated their basal gene expression at both transcriptomic and proteomic levels. We found that several conserved regulators, including osmolytes biosynthesis genes, transthyretin-like family genes, V-type H+-transporting ATPase and potassium channel genes, were involved in both short-term salinity stress response and long-term acclimation processes. In addition, we identified genes related to cell volume regulation, such as actin regulatory genes, Rho family small GTPases and diverse ion transporters, which might contribute to hyposaline acclimation, while the glycerol biosynthesis genes gpdh-1 and gpdh-2 accompanied hypersaline acclimation in L. marina. This study paves the way for further in-depth exploration of the adaptive mechanisms underlying euryhalinity and may also contribute to the study of healthy ecosystems in the context of global climate change.


Introduction
Salinity is an important abiotic environmental factor, which affects survival, growth, reproduction and ecological distribution of living organisms. Efficient sensation, response and adaption to their external salinity environment is vital for all living individuals. The imbalance of salt intake also affects human health, which is associated with a variety of cardiovascular diseases and other physiological pathologies [1][2][3]. Therefore, osmoregulation mechanisms have always been an important part of biology.
In the process of response and adaptation to salinity changes in the surrounding environment, certain universal strategies are applied by different organisms, yet with diversity in the details of regulation. Studies using brewer's yeast as a model, have demonstrated that osmotic sensation and transduction within a single eukaryotic cell can be highly complex, acting in parallel pathways, and often cross-communicate with other signaling processes [4,5]. Yeast cells accumulate glycerol as a compatible osmolyte under hyperosmotic stress, and the high osmolarity glycerol (HOG) response pathway controls glycerol accumulation at various levels including the activation of mitogen-activated protein kinase (MAPK) pathway genes [5][6][7]. Upon acute osmotic shocks, cell volume changes rapidly. Hyperosmotic shock causes cell shrinkage, while hypoosmotic stress leads to cell swelling. Along with volume perturbation, the biophysical properties of the cell membrane, physiological state in the cytosol, and DNA structure as well as gene expression in the nucleus will

Developmental Analysis under Different Salinity Conditions
For developmental analysis, 100 newly hatched L. marina L1 larvae were transferred onto each indicated 35 mm diameter agar plate, seeded with 30 µL OP50. The number of adult worms was scored every 24 h.
Data represent the average of three replicates. Statistical significance was determined using the two-tailed Student's t-test between two groups. p value < 0.05 was considered statistically significant.
For egg laying time analysis, the earliest observed egg laying time on each assay plate was recorded, 100 L1 larvae per plate in six replicates for each condition.

Lifespan Assay
Worms were cultured under optimal growth conditions for at least 3 generations before lifespan assay. The lifespan assay was performed starting at day 1 of adulthood as described previously [38,56], with minor modification. In brief, 35 mm diameter assay plates seeded with 30 µL OP50 were prepared every day. Forty L4 females were transferred to each assay plate, and incubated at 20 • C. The number of live and dead worms was determined using a dissecting microscope every 24 h. Alive worms were transferred to fresh OP50-seeded agar plates daily. Three replicates were analyzed. Worms were scored as dead if no response was detected after prodding with a platinum wire. Dead worms on the wall of the plate were not counted.
Statistical analysis of the average lifespan for worms acclimated to each salinity condition was performed. Data represent the average of three replicates. The comparisons between two groups were performed using the two-tailed Student's t-test. p value < 0.05 was considered statistically significant.

Synchronized L1 Larvae Collection for Each Salinity Condition
The synchronized L1 larvae under each salinity condition were collected as previously reported by Xie et al. [39]. Worms, acclimated to S3, S30 and S50 conditions, were allowed to lay eggs overnight at 20 • C. Eggs were washed off and collected using corresponding suitable solutions, specifically, M9 buffer for the S3 group, filtered sterile seawater for the S30 group, and filtered sterile artificial seawater containing 5% sea salt (Instant Ocean) for the S50 group. The eggs were then 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 the corresponding suitable solution. The clean eggs hatched overnight and developed into L1 larvae in the corresponding suitable solution at 20 • C. The synchronized L1 larvae were collected by filtration using a 500-grid nylon filter with 25 µm mesh size. The samples were frozen immediately in liquid nitrogen.

RNA-Seq Analysis
Total RNA was extracted using Trizol (Invitrogen, Carlsbad, CA, USA). With three biological replicates for each group, a total of nine RNA libraries were prepared with 3 µg RNA using NEBNext ® UltraTM RNA Library Prep Kit for Illumina ® (New England Biolabs, Ipswich, MA, USA) following the manufacturer's instructions. Then, RNA libraries were sequenced on an Illumina NovaSeq 6000 platform and 150 bp paired-end reads were generated.
First, 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.79% and Q30 was over 93.64%. Then, the clean data were aligned to the L. marina reference genome [38] by Hisat2 (v2.0.5, with the default parameters) [57], with mapping ratio from 68.57% to 70.70% (Supplementary Table S2). New transcripts for novel genes were predicted and assembled by StringTie (v1.3.3b, with the default parameters) [58], then annotated with Pfam, SUPERFAMILY, Gene Ontology (GO) and KEGG databases. Briefly, the functional annotation was performed using InterProScan v5. .0 [59] by searching against publicly available databases Pfam (http://pfam.xfam.org/, accessed on 27 July 2020), SUPERFAMILY (http://supfam.org, accessed on 27 July 2020), and GO (http: //www.geneontology.org/, accessed on 27 July 2020), with an E value cutoff of 1 × 10 −5 . KEGG function [60] was assigned using KOBAS 3.0 [61] by best hit (with an E value cutoff of 1 × 10 −5 ) to KEGG database (http://www.genome.jp/kegg/, accessed on 27 July 2020). Further, the read numbers mapped to each gene were analyzed using featureCounts (v1.5.0-p3, with parameter -Q 10 -B -C) [62], and FPKM (expected number of fragments per kilobase of transcript sequence per million base pairs sequenced of each gene) was calculated based on the length of the gene and read counts mapped to this gene, which was used for estimating gene expression levels. Hierarchical clustering for 9 samples performed using the pheatmap package in R and shown in Supplementary Figure S1 indicated the sample preparation was reliable. Subsequently, differential expression analysis of two conditions was performed using the DESeq2 R package (v1.16.1) [63]. The resulting p values were adjusted using the Benjamini and Hochberg approach for controlling the false discovery rate. Genes with an adjusted p value (padj) < 0.05 found by DESeq2 were assigned as differentially expressed. Moreover, GO enrichment analysis and KEGG pathway enrichment analysis for differentially expressed genes (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.

Proteomic Analysis
Worm samples were sonicated three times on ice using a high intensity ultrasonic processor (pulsed at 25% power for 3 s on and 5 s off for 3 min, Scientz) in lysis buffer (8 M urea, 1% protease inhibitor cocktail), with three biological replicates for each group. The remaining debris was removed by centrifugation at 12,000× g at 4 • C for 10 min. Then, the supernatant was collected, and the protein concentration was determined with a BCA kit (Beyotime Biotechnology, Shanghai, China) according to the manufacturer's instructions. Protein samples were digested with trypsin (Promega, Madison, WI, USA) at 1:50 trypsin-to-protein mass ratio overnight. After being desalted by Strata X C18 SPE column (Phenomenex, Torrance, CA, USA) and vacuum-dried, peptides were labelled with a tandem mass tags (TMT) kit (ThermoFisher Scientific, Waltham, MA, USA) according to the manufacturer's protocol.
Next, the TMT-labeled peptides were fractionated by high-pH reverse-phase HPLC using Agilent 300 Extend-C18 column (5 µm particles, 4.6 mm ID, 250 mm length; Agilent, Santa Clara, CA, USA). Briefly, peptides were first separated with a gradient of 8% to 32% acetonitrile (pH 9.0) over 60 min into 60 fractions. Then, the peptides were combined into 9 fractions and dried by vacuum centrifugation.
Further, peptides were identified and quantified by liquid chromatography-tandem mass spectrometry (LC-MS/MS). Briefly, the tryptic peptides were dissolved in solvent A (0.1% formic acid in 2% acetonitrile), directly loaded onto a home-made reversed-phase analytical column (75 µm ID, 15 cm length). The gradient comprised an increase from 9% to 22% solvent B (0.1% formic acid in 90% acetonitrile) over 40 min, 22% to 32% in 14 min and climbing to 80% in 3 min then holding at 80% for the last 3 min, all at a constant flow rate of 450 nL/min on an EASY-nLC 1200 UPLC system (ThermoFisher Scientific). The peptides were subjected to a nanospray ionization (NSI) source followed by tandem mass spectrometry (MS/MS) in Q Exactive HF-X (ThermoFisher Scientific) coupled online to the UPLC. The electrospray voltage was set as 2.2 kV. The m/z scan range was 350 to 1400 for full scan, and intact peptides were detected in the Orbitrap at a resolution of 120,000. Peptides were then selected for MS/MS using the normalized collisional energy (NCE) setting at 28 and the ion fragments were detected in the Orbitrap at a resolution of 30,000. A data-dependent procedure that alternated between one MS scan followed by 20 MS/MS scans with a 15 s dynamic exclusion was used. Automatic gain control (AGC) was set at 1 × 10 5 . The fixed first mass was set as 100 m/z.
The resulting MS/MS data were processed using the Maxquant search engine (v1.5.2.8). Tandem mass spectra were searched against the 17,661 protein database of L. marina [38] concatenated with the reverse decoy database. Trypsin/p was specified as the cleavage enzyme allowing up to 2 missing cleavages and 5 modifications per peptide. The mass tolerance for precursor ions was set as 10 ppm in the first search and 5 ppm in the main search, and the mass tolerance for fragment ions was set as 0.02 Da. Carbamidomethyl on Cys was specified as fixed modification, and acetylation on protein N-terminal, oxidation on Met, deamidation on Asn and Gln were specified as variable modifications. Minimum peptide length was set at 7. The quantitative method was set to TMT 10plex, and the false discovery rate (FDR) for protein identification was adjusted to < 1%. All the other parameters in MaxQuant were set to default values. A p value < 0.05 from t-tests and a fold change > 1.3 or < 1/1.3 were set as the thresholds for significantly differentially expressed proteins (DEPs).
A total of 306,324 spectrums were acquired, of which 78,015 unique spectrums were obtained, and a total of 45,669 peptides were identified (Supplementary Table S3). The length distribution analysis of peptides showed that most of them consisted of 7-20 amino acids, which is in accordance with the quality control requirements (Supplementary Figure S2). Subsequently, hierarchical clustering for 9 samples are shown in Supplementary Figure S3, indicated the reliable sample preparation.
Moreover, different databases were selected for protein functional annotation. GO annotation proteome was derived from the UniProt-GOA database (http://www.ebi.ac. uk/GOA/, accessed on 30 July 2020). If the proteins were not annotated by UniProt-GOA database, the InterProScan (https://www.ebi.ac.uk/interpro, accessed on 30 July 2020) would be applied to annotate by protein sequence alignment method, which was also used for protein domain functional description. KEGG online service tools KAAS was utilized to predict the pathways in which the identified proteins were involved. Then the annotation results were mapped on the KEGG pathway database using KEGG mapper. Subcellular localizations of the protein were predicted by wolfpsort (http://wolfpsort.seq.cbrc.jp/, accessed on 30 July 2020). COG annotation of the protein was achieved by eggnog-mapper software (v2.0) with the default parameters.
Additionally, enrichment analysis of GO and KEGG pathway was conducted for DEPs by Python, using an in-house script. Two-tailed Fisher's exact test was employed to test the enrichment of DEPs against the background of all identified proteins and a corrected p value < 0.05 was considered significantly enriched. Ratio was defined as the ratio of the number of differential proteins annotated to the GO term or on the KEGG pathway to the total number of differential proteins, respectively.

L. marina Is a Typical Euryhaline Marine Nematode
In the laboratory, marine nematode L. marina HQ1 is normally cultured on SW-NGM agar plates prepared with seawater with a salinity of 3%, and this condition is referred to as "S30" in this paper. Basic developmental characteristics were tested in the laboratory, including the adulthood percentage, earliest egg laying time and lifespan. Under the S30 condition, although only 3% of newly hatched L1 larvae had developed into the adult stage within 3 days at 20 • C, the adulthood percentage was as high as 92% within a week ( Figure 1A). The earliest observed egg laying time was around 84 h post L1 stage ( Figure 1B). Moreover, L4 worms lived as long as 29 days ( Figure 1C), with an average lifespan of about 16 days ( Figure 1D).
Interestingly, HQ1 worms acclimated quite well to the low-salinity condition (0.3% NaCl, referred to as "S3") which is the standard culture condition for terrestrial C. elegans. We observed that, under S3 condition, 2% of L1s developed into adulthood within 3 days, the adulthood percentage was about 92% on the 7th day, and the earliest egg laying time was around 83 h post L1 stage, which showed no significant differences compared with the S30 group ( Figure 1A,B). However, worms' lifespan under S3 were slightly shorter: L4 worms lived as long as 27 days ( Figure 1C), with an average lifespan of about 13 days ( Figure 1D). HQ1 worms were able to propagate and acclimate to an even higher salinity environment (5% artificial sea salt, referred to as "S50"). Notably, L1 worms could not develop into the adult stage within 3 days under the S50 condition. The percentage of adulthood on the 4th and the 5th days also showed significant differences between S50 and S30 groups, while no differences were observed afterwards ( Figure 1A). The earliest egg laying time for S50 worms was around 98 h post L1 stage, which showed a severe egg laying delay comparing with both S30 and S3 worms ( Figure 1B). L4 worms acclimated to the S50 condition could live as long as 29 days ( Figure 1C), with an average lifespan of about 15 days ( Figure 1D), which was similar to that of S30 worms.
Taken together, HQ1 worms could acclimate to a wide range of salinity and are one typical euryhaline marine nematodes. How the worms regulate gene expression to acclimate to different salinity environments, remains unknown. Next, we applied newly hatched L1s to investigate the basal differences at both transcriptomic and proteomic levels, respectively.

Analysis of the Basal Transcriptome for L. marina Acclimated to Different Salinity Environments
To investigate the basal transcriptomic differences of L. marina growing under different salinity environments, we applied newly hatched L1s for RNA-seq analysis. Significant DEGs were identified from different comparison groups (Figure 2A), with |log 2 foldchange| > 1 and DESeq2 padj < 0.05 setting as the differential gene screening thresholds. Details of these DEGs were listed in Supplementary File S1.
In the low-salinity S3 group, a total of 1191 genes were significantly up-regulated while 773 genes were down-regulated when compared with S30 group (Figure 2A). Based on GO enrichment analysis, significant enrichment was only observed within down-regulated DEGs. Such genes were mainly annotated to "chromosome", "cellular macromolecular complex assembly", "cellular component biogenesis", "intracellular non-membranebounded organelle", "plasma membrane" and "extracellular matrix" ( Figure 2B). Moreover, "beta-alanine metabolism" pathway-related genes were significantly enriched among upregulated DEGs via KEGG pathway enrichment analysis (Supplementary Figure S4).
including the adulthood percentage, earliest egg laying time and lifespan. Under the S30 condition, although only 3% of newly hatched L1 larvae had developed into the adult stage within 3 days at 20 °C, the adulthood percentage was as high as 92% within a week ( Figure 1A). The earliest observed egg laying time was around 84 h post L1 stage ( Figure  1B). Moreover, L4 worms lived as long as 29 days ( Figure 1C), with an average lifespan of about 16 days ( Figure 1D). to different salinity conditions. One hundred newly hatched L1s were transferred onto each indicated 35 mm dimeter agar plate seeded with 30 μL OP50. The earliest observed egg laying time on each plate was recorded. Six replicates were performed for each experimental condition. (C) For lifespan assay, 40 L4 females were transferred to each assay plate, incubated at 20 °C. The number of live and dead worms was determined using a dissecting microscope every 24 h. Live worms were transferred to fresh OP50-seeded plates daily. Three replicates were analyzed. Worms were scored as dead if no response was detected after prodding with a platinum wire. Dead worms on the wall of the plate were not counted. (D) Average lifespan of L. marina acclimated to different salinity conditions. Error bars represent the standard error of the mean from replicated experiments. Differences between groups were analyzed to different salinity conditions. One hundred newly hatched L1s were transferred onto each indicated 35 mm dimeter agar plate seeded with 30 µL OP50. The earliest observed egg laying time on each plate was recorded. Six replicates were performed for each experimental condition. (C) For lifespan assay, 40 L4 females were transferred to each assay plate, incubated at 20 • C. The number of live and dead worms was determined using a dissecting microscope every 24 h. Live worms were transferred to fresh OP50-seeded plates daily. Three replicates were analyzed. Worms were scored as dead if no response was detected after prodding with a platinum wire. Dead worms on the wall of the plate were not counted. (D) Average lifespan of L. marina acclimated to different salinity conditions. 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. p < 0.05 was considered statistically significant. * p < 0.05, ** p < 0.01, *** p < 0.001, ns-not significant.
There were relatively fewer DEGs identified in the S50 versus S30 comparison group compared with the S3 versus S30 comparison group (Figure 2A). Compared with the S30 group, 114 genes were significantly up-regulated in S50 group, of which "purine metabolism" pathway related genes were significantly enriched (Supplementary Figure S4). In addition, there were 185 DEGs exhibiting significant down-regulation. We found that genes annotated to "chromosome", "cellular component biogenesis", "intracellular organelle part", "intracellular non-membrane-bounded organelle", "serine-type endopeptidase inhibitor activity", "structural constituent of cuticle" and other GO terms were significantly enriched within these down-regulated DEGs ( Figure 2B). Volcano plot of identified DEGs in different comparison groups via RNA-seq analysis (|log2fold-change| > 1; DESeq2 padj < 0.05 was set as the differential gene screening threshold). Red dots (Up) represent significantly up-regulated genes, green dots (Down) represent significantly down-regulated genes, blue dots (NoSig) represent insignificantly differentially expressed genes. (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 was achieved using the 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. (C) Volcano plot of identified DEPs in different comparison groups via proteomic analysis (Ratio of fold change >1.3 or <1/1.3; p value < 0.05 was set as the differential protein screening threshold). Red dots (Up) represent significantly up-regulated proteins, blue dots (Down) represent significantly down-regulated proteins, grey dots (NoSig) represent insignificantly differentially expressed proteins. (D) GO enrichment analysis for DEPs. Ratio of fold change >1.3 or <1/1.3; p value < 0.05 was set as the differential protein screening threshold. A corrected p value < 0.05 was Figure 2. Differentially expressed genes (DEGs) and proteins (DEPs) identified via transcriptomic and proteomic analysis of L1 marine nematodes acclimated to low and high salinity conditions. (A) Volcano plot of identified DEGs in different comparison groups via RNA-seq analysis (|log 2 foldchange| > 1; DESeq2 padj < 0.05 was set as the differential gene screening threshold). Red dots (Up) represent significantly up-regulated genes, green dots (Down) represent significantly down-regulated genes, blue dots (NoSig) represent insignificantly differentially expressed genes. (B) GO enrichment analysis for DEGs. |log 2 foldchange| > 1; DESeq2 padj < 0.05 was set as the differential gene screening threshold. GO enrichment analysis of DEGs was achieved using the 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. (C) Volcano plot of identified DEPs in different comparison groups via proteomic analysis (Ratio of fold change > 1.3 or < 1/1.3; p value < 0.05 was set as the differential protein screening threshold). Red dots (Up) represent significantly up-regulated proteins, blue dots (Down) represent significantly down-regulated proteins, grey dots (NoSig) represent insignificantly differentially expressed proteins. (D) GO enrichment analysis for DEPs. Ratio of fold change > 1.3 or < 1/1.3; p value < 0.05 was set as the differential protein screening threshold. A corrected p value < 0.05 was considered significantly enriched. The color represents the significance of the enrichment. Ratio was defined as the ratio of the number of differential proteins annotated to the GO term to the number of background proteins.

Analysis of the Basal Proteome for L. marina Acclimated to Different Salinity Environments
In parallel, we also applied quantitative proteomic analysis to investigate the basal protein differences among worms growing under different salinity environments. Newly hatched L1s were used, the same developmental staged worms for above transcriptomic analysis. A total of 6068 proteins were identified (Supplementary Table S3) and significant DEPs were selected from different comparison groups ( Figure 2C), with ratio of fold change > 1.3 or < 1/1.3 and p value < 0.05 as the differential protein screening thresholds. Details of these DEPs are listed in Supplementary File S2.
In the low-salinity S3 group, 144 up-regulated proteins and 168 down-regulated proteins were selected, compared with the S30 group ( Figure 2C). Among these up-regulated DEPs, ribosome-related proteins were the most significantly enriched ones revealed by both GO enrichment and KEGG pathway enrichment analysis ( Figure 2D and Supplementary Figure S5), others such as "vacuolar membrane" proteins, "dauer larval development" proteins, "isomerase activity" proteins, "negative regulation of mitotic cell cycle" proteins, "nematode larval development" proteins, "positive regulation of catabolic process" proteins and "phosphatidylinositol phosphate binding" proteins were also significantly enriched ( Figure 2D). However, "extracellular region" proteins were the most significantly enriched ones in down-regulated proteins. Besides, cytoskeleton related proteins (Supplementary  Table S4), defense response related proteins, "DNA packaging" proteins, "modified amino acid binding" proteins, "metal ion binding" proteins and "divalent inorganic cation transport" proteins were also significantly enriched ( Figure 2D).
Compared with the S30 group, there were 56 proteins significantly up-regulated in the S50 group ( Figure 2C), of which ribosome-related proteins were also significantly enriched via both GO enrichment and KEGG pathway enrichment analysis ( Figure 2D and Supplementary Figure S5). Moreover, "chromosomal region" proteins, "midbody" proteins (Supplementary Table S4) and "phosphoric ester hydrolase activity" proteins were also significantly enriched among these up-regulated DEPs ( Figure 2D). Additionally, 105 down-regulated DEPs were found in the S50 group ( Figure 2C). Proteins related with "extracellular region" exhibited the most significant enrichment similar to the low-salinity S3 group. Microtubule-related proteins (Supplementary Table S4), defense response related proteins and "hydrolase activity, acting on ester bonds" proteins were also significantly enriched among down-regulated DEPs ( Figure 2D).
In addition, we also performed a comparison between the S3 and S50 groups, and found 154 up-regulated DEPs in the S3 group and 143 in the S50 group ( Figure 2C). Based on KEGG pathway enrichment analysis, "lysosome" pathway-related proteins and "DNA replication" pathway-related proteins were the most significantly enriched ones in the S3 and S50 groups, respectively (Supplementary Figure S5). Besides, based on GO enrichment analysis, "response to nicotine" proteins, "proton-transporting V-type ATPase complex" proteins (including VHA-5 and VHA-6), "whole membrane" proteins, "regulation of cellular catabolic process" proteins, "spindle microtubule" proteins (Supplementary Table S4), "negative regulation of mitotic cell cycle" proteins, "regulation of proteolysis" proteins, "regulation of synapse organization" proteins exhibited significant enrichment in the S3 group ( Figure 2D). In contrast, "MCM complex" proteins, "condensed chromosome" proteins, "extracellular region" proteins, "cation transport" proteins, "nucleotide metabolic process" proteins, "cell division site" proteins, "hydrolase activity, acting on glycosyl bonds" proteins and others were significantly enriched in the S50 group ( Figure 2D).

Identification of Genes Expressed Consistently at Both mRNA and Protein Levels in Different Salinity Environments
Interestingly, not only was the salinity difference between the S50 and S30 conditions smaller than that between the S3 and S30 conditions, but we also noticed a similar tendency in the number of DEGs and DEPs identified within different comparisons (Figure 2A,C), i.e., the difference between the S50 and S30 groups was smaller than that between the S3 versus S30 comparison pair, especially for DEGs. Herein, we will refer to these worms which were acclimated to S3 environment as the "low-salinity group", while worms grown under S30 and S50 environments as the "high-salinity groups". In order to identify crucial genes related to euryhalinity, we tried to screen genes that were expressed consistently at both mRNA and protein levels.
As shown in Figure 3A, 1638 genes were up-regulated specifically in the low-salinity group (S3), compared with the high-salinity groups (S30 and S50), with the screening thresholds set to fold change > 1.0, padj < 0.05 applied to the transcriptomic data. Similarly, with the screening thresholds set to fold change > 1.0, p value < 0.05 applied to the proteomic data, a total of 354 proteins were specifically up-regulated in the low-salinity group. Furthermore, we combined the results from transcriptomic and proteomic analysis and found that 78 genes exhibited consistent expression at both mRNA and protein levels. Interestingly, we also found that 29 out of the above 78 genes, including the trehalose biosynthesis gene tps-2, the transthyretin-like family genes ttr-15 and ttr-30, the ion transport genes mca-1, twk-33 and vha-5, demonstrated specific induction upon hyposaline stress in L. marina based on our previous data (Supplementary Table S5). Thus, these 78 genes were considered low-salinity-specific genes, and their detailed information is summarized in Table 1.
Genes 2022, 13, x FOR PEER REVIEW 11 of 26 Figure 3. Identification of genes expressed consistently at both mRNA and protein levels under lowand high-salinity conditions. (A) Identification of 78 genes with up-regulated expression at both mRNA and protein levels under the low-salinity condition (S3 group) compared with high-salinity conditions (S30 and S50 groups). DESeq2 padj < 0.05; fold change >1 was set as the differential gene screening threshold. p value < 0.05; Ratio > 1 was set as the differential protein screening threshold.
(B) Identification of 69 genes with up-regulated expression at both mRNA and protein levels under high-salinity conditions (S30 and S50 groups) compared with the low-salinity condition (S3 group). DESeq2 padj < 0.05; fold change >1 was set as the differential gene screening threshold. p value < 0.05; Ratio > 1 was set as the differential protein screening threshold. (C) Classification of differential genes expressed under low-and high-salinity conditions based on COG annotation. . DESeq2 padj < 0.05; fold change > 1 was set as the differential gene screening threshold. p value < 0.05; Ratio > 1 was set as the differential protein screening threshold.
(B) Identification of 69 genes with up-regulated expression at both mRNA and protein levels under high-salinity conditions (S30 and S50 groups) compared with the low-salinity condition (S3 group). DESeq2 padj < 0.05; fold change > 1 was set as the differential gene screening threshold. p value < 0.05; Ratio > 1 was set as the differential protein screening threshold. (C) Classification of differential genes expressed under low-and high-salinity conditions based on COG annotation.

Identification of Genes and Their Corresponding Proteins, the Abundance of Which Was Proportional to Environmental Salinity
In the present study, we found that DEGs enriched markedly based on RNA-seq analysis, hardly exhibited enrichment at the protein level ( Figure 2B,D), in fact, only limited genes were expressed consistently at both mRNA and protein levels ( Figure 3A,B). Next, we focused on those genes whose abundance was directly or inversely proportional to environmental salinity at both mRNA and protein levels, which were probably key regulators or effectors associated with the effective osmoregulation of euryhaline L. marina.
We therefore combined transcriptomic and proteomic profiles and further found that 66 genes and their corresponding proteins demonstrated environmental salinity-dependent patterns in expression (Supplementary File S3). Specifically, 38 genes were down-regulated when salinity increased, while 28 genes were up-regulated when salinity increased (Figure 4 and Supplementary File S3).  Note: Genes marked with a red asterisk can be specifically induced by high-salinity stress, as shown in Supplementary Table S6 analyzed based on data from our previous study [39].
Together, the above genes could be key regulators or effectors in L. marina, involved in its acclimation process to different salinity environments, and worthy of further in-depth functional study in the future. s 2022, 13, x FOR PEER REVIEW 18 of 26

Discussion
Effective osmotic regulation not only directly affects animals' survival, but also shapes their behaviors and distribution. There are various euryhaline fish in nature, such as eels, which spending most of their lives in freshwater until they return to their spawning grounds in the sea, whereas salmon migrate from ocean through their natal river for spawning [64][65][66]. During freshwater to seawater transition or vice versa, euryhaline fish although cope with external salinity changes in a species-specific way, evolutionary conserved strategies do exist among them. The gill, intestine and kidney are the major

Discussion
Effective osmotic regulation not only directly affects animals' survival, but also shapes their behaviors and distribution. There are various euryhaline fish in nature, such as eels, which spending most of their lives in freshwater until they return to their spawning grounds in the sea, whereas salmon migrate from ocean through their natal river for spawning [64][65][66]. During freshwater to seawater transition or vice versa, euryhaline fish although cope with external salinity changes in a species-specific way, evolutionary conserved strategies do exist among them. The gill, intestine and kidney are the major osmoregulatory tissues in fish and successful acclimation in both freshwater and seawater environments depends on proper physiological, metabolic and structural adjustments within these tissues, which also involves neuroendocrine regulation [29,67]. Euryhaline fish have the ability to perceive salinity changes through multiple osmo-sensors, including transmembrane proteins and cytoskeletal proteins, which results in early osmotic response and regulation processes and the following regulatory expression and activities of diverse genes and corresponding proteins involved in water-and ion transport, macromolecular damage control, accumulation and transport of organic osmolytes and other physiological processes together contribute to cellular stress response, cell volume regulation and tissue remodeling [64,[67][68][69][70]. In addition to the final physiological acclimation to external salinity, there are also phenotypic differences in behavior and body morphology or size between marine and freshwater populations, as reported in Fundulus and three-spined sticklebacks [71,72]. Similar regulatory mechanisms have also been reported in crustaceans, such as crabs and shrimps [30,31,73].
When animals encounter a stress-inducing habitat, they will first trigger early stress responses and then might adapt to the new environment through subsequent adaptive regulation. Our previous report showed that L. marina L1 larvae were paralyzed immediately upon both low-and high-salinity conditions, and the stressed worms exhibited developmental defects afterwards [39]. In the present study, we observed that L. marina could move, grow and propagate normally after acclimation to either hyposaline or hypersaline environments, except for a relatively shortened lifespan for hyposaline-acclimated worms and relatively delayed development of hypersaline-acclimated worms. We then analyzed the basal transcriptomic and proteomic differences among L. marina worms acclimated to different salinities, aiming to provide insight into invertebrate euryhalinity.

Cellular Stress Response Genes Might Play Essential Roles in L. marina Euryhalinity
Several basic aspects of the osmotic-induced cellular stress response are well documented, including regulation of the cell cycle and reallocation of metabolic energy [70]. Here, diverse genes related to "cell cycle control, cell division, chromosome partitioning" and "replication, recombination and repair", were identified in both low-and high-salinityacclimated groups, as summarized in Figure 3C. Similarly, metabolism-related genes associated with "amino acid transport and metabolism", "carbohydrate transport and metabolism", "coenzyme transport and metabolism", "lipid transport and metabolism", "nucleotide transport and metabolism" and "energy production and conversion" also exhibited differential expression patterns in worms growing under hyposaline and/or hypersaline-acclimated conditions ( Figure 3C). Interestingly, we found that many of these genes, especially the metabolism-related ones, could also be induced by early short-term salinity stresses [39], for instance, bus-2, gdh-1, pnp-1, tatn-1, C01B10.3, F08F3.4, Y43F8C.13 and Y71F9B.9 were responsive to hyposaline stress (Table 1 and Supplementary Table S5), while gpdh-1, F37C4.6 and Y71G12B.10 were responsive to hypersaline stress (Table 2 and  Supplementary Table S6), suggesting the involvement of these genes in both osmotic stress response and acclimation processes in euryhaline L. marina.
Another aspect of the osmotic-induced cellular stress response is programmed cell death [70]. Previously, we reported a group of TTL genes presumably involved in apoptosis and damage control regulation in response to both low-and high-salinity stresses in L. marina [39]. Here, we found that diverse TTL genes were differentially expressed in acclimation to different salinity environments. A total of 31 out of about 47 annotated TTL genes in L. marina's genome exhibited expression differences with significance in at least one comparison group at either mRNA or protein level or both (Supplementary  Table S7). For instance, 10 TTL genes were significantly increased under the low-salinity condition, 4 of which were increased at both mRNA and protein levels (EVM0016470/ttr-59, EVM0008626/ttr-46, EVM0003972/ttr-15 and EVM0004638/ttr-30). In contrast, 17 TTL genes were up-regulated under high-salinity condition(s), 4 of which were significantly up-regulated at both mRNA and protein levels (EVM0002004/ttr-51, EVM0005297/ttr-44, EVM0007550/ttr-30 and EVM0011170/ttr-27). Together, this distinctively differential regulation of diverse TTL genes in either hyposaline-or hypersaline-acclimated nematodes, suggests that certain TTL genes might play essential roles in hyposaline acclimation, while others might play essential roles in hypersaline acclimation. Moreover, 16 out of the above 31 salinity acclimation-related TTL genes were identified in our previous study [39], and could be induced by both hypo-and hyper-osmotic stresses, such as ttr-30, ttr-44, ttr-46 and ttr-59 (Supplementary Table S7). As one of the largest nematode protein families, most TTLs were predicted to be secreted [74,75]. Some TTL genes in C. elegans were responsive to diverse environmental challenges including oxidative stress, pathogen exposure and osmotic stress [52,[76][77][78][79]. According to studies on ttr-52, which was reported as a bridging factor involved in cell corpse engulfment, apoptosis and axon repair [75,80,81], and another TTL gene, ttr-33, was reported as a potential secreted sensor or scavenger of oxidative stress involved in neuroprotective mechanism [82]. We thus speculated that TTL genes presumably function in apoptosis and damage-control machinery involved in the cellular stress response to cope with salinity stresses and play important roles in both salinity response and acclimation processes in euryhaline L. marina. How each TTL gene functions in salinity acclimation deserves to be further explored.

Certain Cell Volume-Regulation Genes Might Contribute to Hyposaline Acclimation in L. marina
Cell volume is known to be affected by osmotic exposure [8,42]. Volume-dependent regulation will be elicited to restore near-normal cell volume to maintain homeostasis after volume perturbation [8]. The cytoskeleton was implicated as a potential volume sensor and a mediator of volume-dependent regulation of various ion transporters and channels [8,83,84]. It has been well documented that cytoskeletal remodeling, especially actin-related, is involved in volume changes upon osmotic stress [11,85]. Previously, we reported that several tubulin genes show significantly opposing transcriptional changes between low-and high-salinity-stressed worms, indicating their potential roles in salinity response in L. marina [39]. In the present study, although some tubulin-related proteins were enriched in different salinity comparison groups based on the proteomic data, they were hardly expressed consistently at the mRNA level ( Figure 2D and Supplementary Table  S4). However, five cytoskeleton-related genes, ben-1, gsnl-1, tni-3, tnt-3 and EVM0008515, exhibited remarkable elevation at both mRNA and protein levels in low-salinity-acclimated worms ( Table 1). Expression of all these five genes was inversely proportional to the environmental salinity, increasing their abundance as salinity decreases ( Figure 4B), suggesting their important roles in hyposaline acclimation in L. marina. Three of these are actin-regulatory genes, for example, tni-3 and tnt-3 are orthologs of human troponin I and troponin T, respectively [86,87], gsnl-1 is an ortholog of the human gelsolin-family gene [88,89]. Both troponin and gelsolin are involved in regulation of actin dynamics in a calcium-dependent manner, playing important roles in various actin-related processes, including cell structure, cell growth, cell motility, intracellular transport and muscle contraction [90][91][92]. Notably, to our knowledge, this is the first time a potential role for troponinand gelsolin-related proteins in hyposaline acclimation has been exhibited in L. marina, or even in marine invertebrates.
Interestingly, corresponding to the expression profiles of the above actin-regulatory genes, one plasma membrane calcium pump gene mca-1 [93,94] also exhibited elevation in both mRNA and protein abundances when salinity decreased ( Figure 4B). As an important component for the maintenance of calcium homeostasis in cells, mca-1 probably has a role in actin cytoskeleton regulation in L. marina. Besides, Rho family small GTPases have been suspected to participate in the osmotically induced remodeling of the cytoskeleton [11,84]. Consistent with these findings, we also found a set of small GTPase-related genes via GO enrichment analysis, which only increased their transcriptional levels under the low-salinity condition (Supplementary Table S8), indicating their potential involvement in cytoskeleton regulation associated with hyposaline acclimation in L. marina.
In marine invertebrates and fish, ion transporters and channels play important roles in osmoregulation [30,31,33,67,95]. One example is the upregulation of V-type H + -transporting ATPase genes in response to low-salinity stress, which have been documented not only in crustaceans, such as the shrimp Litopenaeus vannamei [96] and the mud crab Scylla paramamosain [31], but also in the marine nematode L. marina [39], indicating their conserved function among marine invertebrates. Besides, Bonzi et al. [66] reported that ion transport genes showed significant expression differences in gill tissues between two natural populations of Arabian pupfish (Aphanius dispar), which inhabit nearly fresh water and seawater areas. In our present study, in addition to the V-type H + -transporting ATPase gene vha-5 and a calcium pump gene mca-1 discussed above, expression of a potassium channel gene twk-33 [97] and a potassium/chloride cotransporter gene kcc-1 [98], was also found inversely proportional to the level of salinity, elevating both mRNA and protein abundances when salinity was decreasing ( Figure 4B). The increased abundance of potassium channel and potassium/chloride cotransporter genes was consistent with their general requirement mediating K + and Cl − efflux to defend against cell swelling under hypoosmotic conditions [8,99]. Of note, mca-1, twk-33 and vha-5 also demonstrated specific induction upon hyposaline stress in our previous report [39] (Table 1 and Supplementary Table S5), we therefore proposed that these above-mentioned ion transport-related genes play crucial roles in both osmotic stress response and long-term acclimation processes specifically under low-salinity conditions in L. marina.

The Glycerol Biosynthesis Genes gpdh-1 and gpdh-2 Accompanied Hypersaline Acclimation in L. marina
It is known that the accumulation of organic osmolytes is a ubiquitous mechanism in cellular osmoregulation [26,27]. There are a number of organic osmolytes such as glycerol, trehalose, inositol, betaine and taurine, which allow cells to counteract the effects of hyperosmolarity and to adapt to hyperosmotic conditions. For example, free amino acids and methylamines are mainly utilized by most marine invertebrates [30,31], whereas glycerol is the most important osmolyte for the terrestrial nematode C. elegans [42,46]. It is well documented that an increased level of glycerol is essential for C. elegans' survival in hypertonic environments mediated by upregulation of the glycerol biosynthetic enzyme gene gpdh-1 [43,46]. Burton et al. reported that if C. elegans was exposed to mild highsalinity stress, its progeny could be protected from strong osmotic stress via increasing the expression of gpdh-2, and this glycerol synthesis gene is essentially required for the intergenerational osmotic protection [45]. In line with its terrestrial relative C. elegans, we previously demonstrated that gpdh-1 was significantly induced by high-salinity stress in L. marina, suggesting that glycerol as an essential osmolyte in both nematodes' response to hyperosmotic stress [39]. In the current study, two glycerol synthesis genes, gpdh-1 and gpdh-2, both showed significant up-regulation in high-salinity groups at both mRNA and protein levels ( Table 2 and Figure 4B), whose expression was directly proportional to the level of salinity. Given gpdh-2 was not notably induced by salinity stress, our results indicate that gpdh-2 might play a role in intergenerational osmotic protection inheritance and hyperosmotic acclimation in the marine nematode L. marina. Together, our data suggest that both gpdh-1 and gpdh-2 are essential to hypersaline acclimation in L. marina.

Conclusions
In conclusion, we have described for the first time the genome-wide transcriptional and proteomic analysis of the marine nematode L. marina acclimated to either low-or highsalinity conditions. We found that various cellular stress response genes may function as the conserved regulators in both short-term salinity stress response and long-term acclimation processes. Additionally, we identified diverse genes involved in cell volume regulation which might contribute to hyposaline acclimation, and we also demonstrated that genes related to glycerol biosynthesis might accompany hypersaline acclimation in L. marina. Thus, our data might lay the foundation to identify the key gene(s) for further in-depth exploration on environmental adaptation mechanisms in euryhaline organisms, especially in the context of global climate change and the corresponding marine salinity stresses.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13040651/s1, Figure S1: Heatmap of expression levels for all DEGs among nine samples; Figure S2: Length distribution of peptides identified by mass spectrometry; Figure S3: Heatmap of expression levels for all DEPs among nine samples; Figure S4: KEGG pathway enrichment analysis for DEGs identified via RNA-seq; Figure S5: KEGG pathway enrichment analysis for DEPs identified via proteomic analysis; Table S1: Recipes for agar plates used for different salinity culture conditions; Table S2: Statistics of mapping ratio; Table S3: Overview of the LC-MS/MS analysis results; Table S4: Expression profiles of enriched tubulin related proteins based on proteomic data via GO enrichment analysis; Table S5: 29 out of 78 low-salinity-specific genes showed specific induction upon hyposaline stress in L. marina based on our previous data; Table S6: 11 out of 69 high-salinity-specific genes showed specific induction upon hypersaline stress in L. marina based on our previous data; Table S7: Expression profiles of 47 annotated transthyretinlike (TTL) family genes in L. marina; Table S8: Up-regulation of small GTPase-related genes at transcriptional level under low-salinity condition; File S1: Identified DEGs via RNA-seq analysis; File S2: Identified DEPs via proteomic analysis; File S3: Sixty-six genes expressed proportionally to environmental salinity.

Data Availability Statement:
The RNA-seq raw data were submitted to the NCBI under the accession code PRJNA778902. All raw proteomics data are available via ProteomeXchange under identifier PXD029671.