1 Introduction

Salinity in coastal seawater and estuaries shifts fast, which is significantly influenced by multiple factors. For instance, the influx of freshwater from land runoff (Reguero et al., 2015) plays an instrumental role in shaping the dynamics of salinity, coupled with the tide and seasons (Schumann et al., 2006; Benham and Cheviron, 2020). The salinity is one of the main environmental factors shaping microbial community (Song et al., 2022). Salinity undulation shifts the community structure of both benthic and planktonic phytoplankton assemblages and alters their ecological and physiological performances (Velasco et al., 2019).

Protist is an essential component of the microbial food web, playing a central role there in biogeochemical cycling, and sustaining the healthy function of ecosystems (De Ruiter et al., 1995; Tsai et al., 2013). Deciphering the acclimation and adaptation of protists to saline stressing is important for evaluating both the flexibility and functional dynamics of ecosystems. Unicellular protists are vulnerable to salinity stress; they grapple with the ionic imbalance and osmotic stress induced by high salinity. Numerous studies have examined the saline tolerance and stress response of microalgae (e.g., green algae, diatoms), in which their ecological and physiological performances under saline stress and their capacity for salinity acclimation are determined from the perspectives of morphology, physiology and genomics (Neelam and Subramanyam, 2013; Fang et al., 2017; Li et al., 2019a). Some microalgae sensitive to salinity grow slowly and survive less in environments with high salinities (Talebi et al., 2013). Certain salt-tolerant algae are found to alter their morphology and osmolyte contents in shortterm salt stressing and accumulate favorable mutations in long-term salt stressing (Haris et al., 2022). To date, most studies focused on microalgae while very less attention was paid to heterotrophic protists represented mainly by the predators, e.g., ciliates, flagellates and amoeboid. These organisms play a vital role in nutrient transfer to high trophic levels within aquatic and soil food webs (De Ruiter et al., 1995).

At present, data of the saline tolerance for most phyla of heterotrophic protists are quite limiting (Ekelund, 2002; Booton et al., 2004; Hauer and Rogerson, 2005). Recent publications have dedicated mainly to the mechanisms underlying the osmotic adaptation of a limited number of halophilic/halotolerant flagellates and ciliates (Harding et al., 2016; Weinisch et al., 2018, 2019). These unicellular eukaryotes commonly employ the excretion strategy in coping with the persistent high salinity, and accumulate organic compatible solutes, e.g., glycine betaine, myo-inositol and ectoine, as the osmoprotectants (Harding et al., 2016; Weinisch et al., 2018). However, these mechanisms mainly function in the halo adaptations of organisms that have undergone aeons of evolutionary time. Therefore, these mechanisms are not applicable to explain the short-term response and the tolerance of freshwater protists to salinity stressing.

At the mixing areas of freshwater and seawater, e.g., the coasts and the estuaries, salinity alters frequently, and a salinity sustains shortly, hours for tides and months among seasons. The water in these areas is usually hyposaline (salinity 3 – 14) (Hou et al., 2020), not favorable for freshwater protozoans. Certain haloduric freshwater protozoans may survive salinity stressing; unfortunately, the studies on their physiological response and adaptation to such stress are scarce, and the studies on their genomic response and adaptation are even less. For these haloduric freshwater protozoans, our curiosities include their tolerance range to salinity, the balance between their growth and salinity tolerance, their acclimation capacity referring halophile species, and the genomic mechanism supporting their physiological response to salinity.

Diverse ciliates inhabit various salinities (Zhao and Xu, 2016) and represent the best studied phagocytic protists. Therefore, they can be used as a model for salinity acclimation in heterotrophic protozoa. Less previous studies focused on the salinity tolerance of two model ciliates, Paramecium and Tetrahymena (Smurov, 2000; Podlipaeva et al., 2008). For others, no information in the physical response and assimilation to saline is available, and the associating genomic mechanism is not clear. In this study, we tried to understand the response and acclimation to saline stressing of a non-model hypotrichid ciliate, Gastrostyla setifera, by comparing its transcriptome before and after a mid-term (> 30 d) acclimation to a low salinity. In our preliminary studies, we found that the ciliate survives and even grow when it is cultured at salinity 3. RNA-seq of 30 – 100 cells were used in order to overcome the less availability of ciliate individuals. From the cultured ciliates, the transcriptome of G. setifera was assembled de novo, on which we deciphered its response to low salinity. Our findings provided insights into the physiological mechanism underlying its response to saline stressing.

2 Materials and Methods 2.1 Ciliate and Its Culture

The ciliate Gastrostyla setifera was collected in July, 2020, from a freshwater lake in Yantai University campus (37˚28΄N, 121˚27΄E). The water temperature was 23℃. Cells were isolated under a stereo microscope (Guiguang GL-6345BI, China) and cultured in freshwater at 25℃. Single individual cells were manually selected using a micropipette and transferred into the medium with different salinities (0 and 3). The 3-salinity medium was prepared by diluting sterile seawater with double-distilled water, whereas the freshwater medium was made using sterilized double distilled water. We found that the ciliate survives in 3 saline water for up to six months at 25℃. Three replicate cultures were set each salinity. Sterilized rice was added to the culture to amplify bacteria as the feed of the ciliate. Ahead of RNA extraction, the ciliate was cultured in freshwater at 25℃ until mid-exponential growth phase. For saline stressing, the ciliate was pre-assimilated at 25℃ for more than 30 days with the individuals inoculated into a new saline medium, and grew into mid-exponential growth phase, during which the cells were collected for RNA extraction.

2.2 RNA Extraction, cDNA Library Construction and Transcriptome Sequencing

In total, 50 – 100 ciliate individuals (cells) were manually isolated from the culture, washed three times with sterilized water, and starved for two hours. Approximately 30 cells and 100 cells were collected from freshwater and saline water culture, respectively, by centrifuging at 650 g for 4 min. The cells were lysed with 2 μL lysis buffer. The gene expression was immediately stopped by placing the lysate on ice. The total RNA was extracted and reversely transcribed into cDNA. The cDNA was amplified using QI-Aseq FX Single Cell RNA Library Kit (Qiagen, Hilden, Germany) following manufacturer's instructions. The quality and quantity of the cDNA were assessed on NanoDrop 2000C Spectrophotometer (Thermo, Wilmington, DE, USA) and through electrophoresis in 1% agarose gel. The qualified cDNA was sent for sequencing by Novogene Company (Tianjin, China). A library was constructed each replicate. Sequencing was performed on an Illumina HiSeq 2000 platform.

2.3 Data Processing and Transcriptome Analysis

The reads containing adapters were trimmed and those containing either > 10% N or 50% bases with ≤ 5 Q-value were filtered out using Trimmomatic v0.33 with LEADING: 5 TRAILING: 5 SLIDINGWINDOW: 5:16 MINLEN: 80 as parameters (Bolger et al., 2014). Clean reads were assembled into unigenes (genes) using Trinity v2.4.0 (Grabherr et al., 2011) with default parameters. To eliminate rRNA and potential contaminants, the assemblies were blasted against NCBI-nt and NCBI-nr databases using blastN and blastX using rRNA/18S rRNA/16S rRNA/5.8S rRNA/23S rRNA/ribosomal RNA/26S rRNA/ribosomal DNA/ITS/LSU/SSU as the keywords. A sequence was identified as the putative contaminant if it is over 50% in length and at least 90% in sequence similar with a sequence from nonprotists (Kolisko et al., 2014).

The ciliate specific unigenes were clustered using CD-HIT v4.7 with CD-HIT-EST, -c 0.9 as the parameters (Fu et al., 2012) to eliminate the redundancy of contigs with a sequence similarity threshold of 90%. The unigenes were annotated against databases including NCBI nonredundant protein sequence (nr), NCBI nucleotide sequence (nt/nr), Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) using BLAST with a cut off E-value of < 10⁻⁵. In GO annotation, Benjamin-Hochberg correction for multiple testing was used to control the P-values, and the GO term fusion option was applied to reduce the term complexity of GO. The KEGG Ortholog (KO) assignments and pathway maps were obtained using the bidirectional best hit method (BBH) on the KEGG Automatic Ontology Annotation Server (KAAS).

To assess the completeness of the newly assembled G. setifera transcriptome, we searched for 20 tRNA synthetases and a set of 248 eukaryotic core proteins (http://korflab.ucdavis.edu/Datasets/genome_completeness/) (Kolisko et al., 2014). Tblastn was used to identify these genes in the transcriptome. The ortholog-grouping method was also applied to assess transcriptome completeness using the 'Orthovenn2' program integrated in the server https://orthovenn2.bioinfotoolkits.net (Xu et al., 2019). Genomes of two model species, Tetrahymena thermophila and Paramecium tetraurelia and two hypotrich species, Oxytricha trifallax and Stylonychia lemnae, retrieved from public databases were used to group genes by orthology (Table 1). Three transcriptomes of Pseudokeronopsis sp. and Euplotes focardi downloaded from MMETSP database were also added into the analysis as Table 1. The reference datasets were BLAST-ed against the predicted peptide dataset of G. setifera with E-value ≤ 1×10⁻² and similarity ≥ 30% and length ≥ 50 aa as the thresholds. Pairwise mutual best hits were identified as the putative orthologs.

The transcript abundance each gene was calculated each replicate using RSEM v1.3.1 (Li and Dewey, 2011) and PERL scripts in the Trinity package based on Bowtie v2.4.2 (Langmead and Salzberg, 2012) following manuals. Differentially expressed genes (DEG) were identified using the edgeR package of R based on gene count data (Squair et al., 2021). To calculate normalization factors among the replicates, we used the 'TMMwsp' method which is more suitable for zero-inflated RNA-seq data. The generalized linear model-likelihood ratio test (GLM-LRT) model was applied in the ANOVA test to assess the significant difference between two salinities (Squair et al., 2021). The false discovery rate (FDR), namely, q-values, was applied to determine the probability of differences between the two salinities. A Q value < 0.05 and |log₂ (Fold change) | > 2 were defined as the thresholds for significant differential expressions (DEGs). To account for gene expression profiles, we normalized the gene count data using the FPKM (Fragments Per Kilobase of exon model per Million mapped reads) method. We visualized the gene expression profiles of DEGs in a heatmap using the 'heatmap' package in R. To test the difference in gene expression profiles between salinities, principal component analysis (PCA) was conducted utilizing the 'PCAFactoMineR' package in R v4.2.2. Additionally, we performed GO and KEGG enrichment analyses using the 'clusterProfiler' package in R v4.2.2. We also utilized several other packages, such as 'ggplot2', 'LSD', 'ggsci' and 'Biostrings', to visualize the processing results of DEGs and select eligible sequences for subsequent analysis.

3 Results 3.1 Growth Performance of G. setifera at Salinity 3

The cells of G. setifera exhibited slower growth at salinity 3 than at 0 (freshwater) (Fig.1). The cell density in saline water was significantly lower than that in freshwater, and the final density was only approximately half of that in freshwater (Fig.1). The ciliate grew to plateau phase in a longer period, about 10 days, in saline water than in freshwater. The cell size and morphology were normal in saline water during the exponential growth phase, but abnormal and even showed the death characteristics after the plateau phase.

3.2 De novo Assembly and Evaluation of G. setifera Transcriptome

In total, 41.9 G raw reads were obtained for five replicates, 3 of freshwater culture, and 2 of saline water culture. One replicate of saline water culture failed in RNA-seq. After filtration, 25.82 G clean reads were retained (Table 2), which were then assembled into 635147 unigenes (Table 3). All data have been deposited in NCBI BioProject with accession no. PRJNA548847. Of the unigenes, 70161 were identified as rRNA, and 79514 matched with genes of prokaryotes and multicellular eukaryotes, which were eliminated as the potential contaminants (Table 3). Twenty aatRNA synthetases and 216 core eukaryotic proteins genes were identified, thus the recovery rate was 100% and 94% for aa-tRNA synthetase genes and 94% for 248 core-genes (Table 3). The clean unigenes were grouped into 458619 clusters, among them, 179227 and 189417 hit the deposited in nt and nr databases, 32793 were annotated against GO database, and 46590 were mapped onto KEGG database (Table 4).

Orthologous genes were identified among G. setifera, two model ciliates, Tetrahymena thermophila and Paramecium tetraurelia, and four phylogenic relatives, Oxytricha trifallax, Stylonychia lemnae, Pseudokeronopsis sp. and Euplotes focardi (Fig.2). The number of orthologous genes (4149) shared by G. setifera, two models and Euplotes focardi was lower (Fig.2A) than that (16248) identified between G. setifera and four hypotrich ciliates (Fig.2B). A proportion of 87% and 80% of the orthologs in Stylonychia lemnae and Oxytricha trifallax were identified in G. setifera (Fig.2C), suggesting the gene discovery rate of the newly assembled G. setifera transcriptome is high.

3.3 DEG Profile of G. setifera Cultured at Salinity 0 and 3 3.3.1 PCA and FPKM pattern

G. setifera cultures at the same salinity clustered together while those at different salinities did not (Fig.3A), implying the expression pattern is different between salinities. The FPKM value of > 45% of the mapped genes of the freshwater samples ranged from 1 to 100, while that of > 60% of the mapped genes of the saline water samples ranged from 0 to 1 (Fig.3B).

3.3.2 Overview of DEGs

Compared to the control (freshwater samples), a total of 4028 significant Up and 7662 Down DEGs were identified in the salinity 3 treatment, with the thresholds of |log (FC)| > 2 and FDR < 0.05 (Fig.3C), which accounted for 22% of the executed unigenes (52384). The FPKM value of each DEG was consistent among replicates within a salinity but significantly different between salinities (Fig.3D).

3.3.3 Functional annotation and GO Enrichment of DEGs

More Down DEGs were identified than Up ones (Fig. 4A). Out of 11690 DEGs, 31.6% were successfully annotated against at least one of the nr, KEGG and GO databases, of them the majority were derived from Down DEGs (Fig.4A). The 1439 DEGs annotated against GO database were significantly enriched in 79 GO terms (P < 0.05). The enriched GO terms in Up and Down DEGs showed distinct profile, and only 1 GO term – phosphorelay sensor kinase, was shared between the two (Fig.4B). A large number of genes involving in binding and phosphorylation activities were over expressed, including calcium ion binding, phosphorelay sensor kinase activity, and phosphorylation (Fig.4B). Significant number of genes related to phosphorelay sensor kinase activity (42 genes) were down regulated. Other Down DEGs were mostly enriched in catalytic activity of peptidase and proteins, comprising metallo- and serine-type endopeptidase, metallopeptidase, and proteolysis (Fig.4B).

3.3.4 KEGG enrichment of DEGs

A total of 1321 up- and 1087 down-regulated DEGs were annotated against KEGG database. The DEGs were significantly enriched in 36 KEGG pathways (P < 0.05), of them, no pathway was shared by Up and Down DEGs (Fig.5). More Up DEGs were enriched than Down DEGs (1321 vs. 1087 genes, Fig.5). A large proportion of genes encoding protein kinases/phosphatase, membrane and cytoskeleton proteins are up-regulated, including protein kinase, membrane trafficking proteins, cytoskeleton proteins, calcium signaling proteins, and protein phosphatases (Fig.5A). High portion of genes involved in carbon metabolism are down-regulated, including citrate cycle (TCA) and glycolysis/gluconeogenesis (Fig.5B). In addition, genes involved in protein synthesis and amino acid catabolism are also down-regulated significantly. These include those relating to ribosomes, amino acid-related enzymes, aminoacyl-tRNA biosynthesis, valine/leucine/isoleucine degradation (Fig.5B).

We further investigated the relative expression levels of annotated DEGs that were enclosed in the enriched KEGG pathways (Figs.6 and 7). The up-regulation of the expression of Up DEGs was slight (Figs.6 and 7); however, the down-regulation of the expression of Down SDEGs was dramatical, and most of them silenced even when the ciliate was cultured at high salinity. The genes with their expression significantly down regulated included those functioning in amino acid metabolism, aminoacyl-tRNA synthesis, citrate cycling, and valine/leucine/isoleucine degradation (Figs.6 and 7).

4 Discussion 4.1 G. setifera Acclimates High Salinity by Reprogramming Its Transcription

Salinity change affects the physiology and performance of planktonic heterotrophic protists. We have found that the growth of G. setifera under saline stress significantly declines (Fig.1), which is in accordance with the performance of other heterotrophic protists, e.g., Paramecium (Smurov, 2000), soil flagellates (Ekelund, 2002), and freshwater amoeba (Anderson et al., 2009).

Twenty two percent of the transcripts in the transcriptomes of G. setifera were found to alter their abundance after salt acclimation (Fig.3C). These DEGs were enriched in different GO terms and KEGG pathways (Figs.3D, 4 and 5). Such alteration reflected the acclimation of G. setifera in responding saline stress. The higher number of Down-DEGs than that of Up-DEGs and the more dramatically regulated expression levels of the Down-DEGs, implied that the ciliate tends to reduce its transcription or to attenuate its metabolism in responding saline stress. Similar adjustment has been documented in marine Daphnia and diatom Thalassiosira weissflogii during their two years saline acclimation (Latta et al., 2012; Bussard et al., 2017). Decreasing transcription cost thus saving energy have been proposed as the physiological mechanism underlining such adjustment of Daphnia and diatom, which should also applicable for the ciliate we studied.

4.2 Up-Regulated Genes May Function in Ion Balancing and Cellular Resistance

We found that a larger portion of Up DEGs of the ciliate at high salinity were enriched in cytoskeleton proteins and membrane trafficking pathways (Figs.5, 6), which associated with microtubule motor and ion channel activities, e.g., potassium ion transport and channel activity (Fig.4). These observations implied that ciliate may have increased ion trafficking through cellular membrane in responding saline stress. Such mechanism is partially similar to that of halophilic protists (Weinisch et al., 2019). The halophilic ciliates and flagellates exclude salts, and acquires and synthesizes rich organic compatible solutes (Harding et al., 2016; Weinisch et al., 2018, 2019).

Mitogen-activated protein kinases (MAPKs) are involved in the most ancient signal transduction pathways of eukaryotes, which can be activated by diverse stimuli including stressing (Jiao et al., 2019; Li et al., 2019b). Studies in plant Pisum sativum have shown that MAPKs phosphorylate various stress-relating proteins, inducing cellular resistance by producing the secondary metabolites. Microalgae highly express MAPK genes under NaCl stressing, indicating MAPKs involve in responding saline stress (Qiao et al., 2021). Protein kinases of bacteria Escherichia coli phosphorylate and thereby activate some proteins which maintain the integrity of cell membrane under salt stress (Joshi et al., 2010). In this work, the protozoa G. setifera cultured in saline water may also enhance its resistance to saline stress by increasing the expression of the MAPK genes. Seventeen Up-DEGs were annotated as MAPKs in our data (Fig.7). The up-regulation of protein kinases in G. setifera may associate with the activation of important proteins, maintaining cell membrane integrity under saline stress, as recorded in plant, microalgae and bacteria (Joshi et al., 2010; Qiao et al., 2021).

Protein phosphatase could regulate biological processes. For example, the serine/threonine-protein phosphatase 2B catalytic subunit (PP2B) plays an essential role in regulating bacterial stress response (Janczarek et al., 2018). Protein phosphatase could interact with signaling proteins, this process has been proved to influence the growth of microbes at high salinity (Ye et al., 2020). In this work, expression of a number of protein phosphatase genes were enhanced in the ciliate G. setifera acclimating salinity, which may regulate its biological processes to resist saline stress (Fig.7).

Over expression of heat shock protein genes (HSPs) stimulated by saline stress has been reported in bacterium Lactococcus lactis (Kilstrup et al., 1997) and antarctic bacterium Psychrobacter sp. (Che et al., 2013), which is believed to associate with bacterial resistance to salinity change (Che et al., 2013). In this work, G. setifera could have also resisted salinity change by higher expression of HSPs. HSPs among UP SDEGs were enriched into both GO term and KEGG pathways (Fig.4B, Fig.7). This result aligns with the findings in the model ciliate Paramecium nephridiatum, in which a high basal level of HSP70 protein accounts for the salt tolerance of the euryhaline strains (Smurov et al., 2007). These findings suggest a potential involvement of HSP as a universal mechanism for salt tolerance in euryhaline ciliates.

4.3 Down-DEGs Involved in Cellular carbon Metabolism and Amino Acid Catabolism

A large portion of Down-DEGs at saline water samples were enriched in carbon metabolism including citrate cycle and glycolysis/gluconeogenesis (Fig.5), and some of them even silenced completely (Figs.6, 7). These genes encode a variety of enzymes relating to energy generation, biosynthesis, and the conversion of organic carbon, suggesting that salinity stress has weakened the major physiological metabolism of the ciliates. As mentioned above, such weakening might function together to save energies and materials to maintain ciliate growth and development under salt stress. The weakening is also coordinated with the decreased growth (Fig.1). Specifically, depression of glycolytic enzyme genes, such as those encoding fructose-bisphosphate aldolase, glucose-6-phosphate isomerase and phosphoglycerate mutase, has been observed in Thalassiosira weissflogii after two years of salinity acclimation (Bussard et al., 2017). This adjustment has been hypothesized to alternate pathways to supply both inorganic and organic carbon to diatom (Bussard et al., 2017). The mechanism might be also applicable for the ciliate G. setifera.

Vast Down-DEGs of the ciliate culture in saline water function in amino acid catabolism, including ribosome, tRNA biosynthesis and valine/leucine/isoleucine degradation (Figs.6, 7). The down-regulation of these processes may contribute to amino acid accumulation in the ciliate cells. Amino acids could serve as the compatible solutes (Harding et al., 2016; Weinisch et al., 2018, 2019). The amino acid accumulation should reduce the osmotic potential between intra- and extra-cellular environments, thus assist the tolerance to high salinity. This mechanism has been documented in halophilic ciliates to minimize their water loss and assist them to maintain their cell turgor pressure (Weinisch et al., 2019). Specifically, two Down-SDEGs encode tyrosyl-tRNA and arginyl-tRNA synthetases (Fig.7), which have been recorded to cause amino acid accumulation in hybrid fish under light stressing (Zhong et al., 2022). Down-regulation of the SDEGs ornithine aminotransferase could lead to accumulation of ornithine (Lopez-Delacalle et al., 2021), which is the precursor of proline, a typical compatible solute found in halophilic ciliates (Weinisch et al., 2018). Moreover, it is known that eukaryotic cells contain about 4 × 10⁶ cytoplasmic ribosomes, accounting for 80% of the total cellular RNA and 5% – 10% of cellular proteins (Yoshihama et al., 2002). The down regulation of ribosomal protein in G. setifera at salinity 3 should decrease ribosome biosynthesis, and indirectly reduction of proteins and RNAs, otherwise lead to reservation of amino acids.

Collectively, the mechanisms of G. setifera accumulating to salinity could be speculated as follow: the ciliate may weaken their processes of protein and RNA synthesis to reduce energy consumption. Salts may be excreted by the higher regulation of ion binding and trafficking proteins. Amino acid may be acquired or accumulated in the cytoplasm to balance the osmotic inter- and intro-cell membrane, via reducing the cells' amino acid catabolism. Meanwhile, the cell also initiated a series signal proteins for protection and to respond to stress, e.g., the protein phosphatases and protein kinase.

5 Conclusions

Protists inhabiting coastal and estuarine environments always confront threats from salinity alterations. However, our understanding of the physiological response of fresh-water heterotrophic protists to salinity stressing remains severely limited. Our work was a pilot study to investigate the physiological mechanism underpinning the response of a freshwater heterotrophic ciliate, Gastrostyla setifera, in a 30 days acclimation to salinity 3. Dramatic alterations in gene expression indicated that the ciliate decreased its transcription, saving energy consumption under salinity stressing. Up SDEGs identified encode protein kinases, phosphatases and cytoskeleton proteins, and Down SDEGs function in amino acid metabolism. All these processes associate with energy saving and osmotic maintaining, which are partly in accordance with the documented in other saline acclimated protist and the halophilic ciliates. Our findings offered essential evidences for the plasticity of micro-eukaryotes and their functioning in coastal and estuarine ecosystems.

Acknowledgements

This study was supported by the National Natural Science Foundation of China (Nos. 32370488, 42176163, 319 70398 and 31672251), and the Youth Innovation Promotion Association of CAS (Nos. 2019216 and 2022211).