Comparative analysis of the liver transcriptome in the red-eared slider Trachemys scripta elegans under chronic salinity stress.

The red-eared slider (Trachemys scripta elegans), identified as one of the 100 most invasive species in the world, is a freshwater turtle originally from the eastern United States and northeastern Mexico. Field investigations have shown that T. s. elegans can survive and lay eggs in saline habitats. In order to understand the molecular mechanisms of salinity adaptation, high-throughput RNA-Seq was utilized to identify the changes in gene expression profiles in the liver of T. s. elegans in response to elevated salinity. We exposed individuals to 0, 5, or 15 psu (practical salinity units) for 30 days. A total of 157.21 million reads were obtained and assembled into 205138 unigenes with an average length of 620 bp and N50 of 964 bp. Of these, 1019 DEGs (differentially expressed genes) were found in the comparison of 0 vs. 5 psu, 1194 DEGs in 0 vs. 15 psu and 1180 DEGs in 5 vs. 15 psu, which are mainly related to macromolecule metabolic process, ion transport, oxidoreductase activity and generation of precursor metabolites and energy by GO (Gene Ontology) enrichment analyses. T. s. elegans can adapt itself into salinity by balancing the entry of sodium and chloride ions via the up-regulation expression genes of ion transport (potassium voltage-gated channel subfamily H member 5, KCNH5; erine/threonine-protein kinase 32, STK32; salt-inducible kinase 1, SIK1; adiponectin, ACDC), and by accumulating plasma urea and free amino acid via the up-regulation expression genes of amino acid metabolism (ornithine decarboxylase antizyme 3, OAZ3; glutamine synthetase, GLUL; asparaginase-like protein 1b, ASRGL; L-amino-acid oxidase-like, LAAO; sodium-dependent neutral amino acid transporter B, SLC6A15s; amino acid permease, SLC7A9) in response to osmotic regulation. An investment of energy to maintain their homeostatic balance is required to salinity adaptation, therefore, the genes related to energy production and conversion (F-ATPase protein 6, ATP6; cytochrome c oxidase subunit I, COX1; cytochrome c oxidase subunit III, COX3; cytochrome b, CYTb; cytochrome P450 17A1, CYP17A1) were up-regulated with the increase of gene expression associated with lipid metabolism (apolipoprotein E precursor, APoE; coenzyme Q-binding protein, CoQ10; high-density lipoprotein particle, SAA) and carbohydrate metabolism (HK, MIP). These findings improve our understanding of the underlying molecular mechanisms involved in salinity adaptation and provide general guidance to illuminate the invasion potential of T. s. elegans into saline environments.

Introduction

The red-eared slider (Trachemys scripta elegans) has been introduced into diverse aquatic habitats worldwide (including many countries in Africa, Asia and Europe as well as Australia) via the pet-release pathway and, as a result, is classified as a highly invasive species by the International Union for Conservation of Nature (Luiselli et al., 1997; Martins, Assalim & Molina, 2014). It is native to freshwater habitats in 19 states of the eastern United States and two states of northeastern Mexico (Mittermeier et al., 2015). Recently these turtles have been found to lay eggs in the low salinity (0.1–26‰) estuary of the Nandujiang in Hainan Province, China (Liu et al., 2011; Yang & Shi, 2014). The extent of the saltwater adaptability of T. s. elegans is not fully understood. Studies of endocrine stress responses by T. s. elegans in the Lake Pontchartrain Basin of Louisiana suggest that these turtles may serve as a sentinel species for elevated salinity in environments where salinity is rising due to saltwater intrusion (Thompson, Franck & Valverde, 2011). These studies indicate that T. s. elegans can invade not only fresh water, but also saline water environments. Therefore, the invasion potential and the mechanism of response to salinity adaptation of T. s. elegans is of serious concern.

Changing levels of salinity are a crucial environmental stress factor for aquatic species that can disrupt electrolyte balance, cell energetics, and various other physiological responses, including activating stress hormones (Lushchak, 2011). Species show altered composition in osmolality of body fluids in response to changing salinity (Charmantier et al., 2011) and adaptation to salinity change typically involves the physiological solution: tolerance of elevated inorganic ion concentrations (mainly sodium and chloride) in plasma (Gordon & Tucker, 1965), and accumulation of organic osmolytes (e.g., urea) to counteract cell-volume changes. The most dramatic changes in urea concentration are seen in plasma and tissues such as skeletal muscle, resulting from the up-regulation of hepatic urea (Wright et al., 2004). Accumulation of intracellular free amino acids via hepatic protein degradation or de novo amino acid synthesis can also contribute to osmotic balance (McNamara et al., 2004; Yancey, 1985). When subjected to ambient salinity change, T. s. elegans increased serum glucose levels, and the activities of creatine kinase (CK), aspartate aminotransferase (AST), lactate dehydrogenase (LDH), and alkaline phosphatase (ALP) in liver (Shu et al., 2012). Our previous studies have shown that T. s. elegans can increase blood osmotic pressure by balancing the entry of sodium and chloride ions with a decrease in the secretion of aldosterone, and by accumulating plasma urea for osmoregulation when ambient salinity was lower than 15‰  (Hong et al., 2014). However the molecular basis of these adaptive responses has not been studied in T. s. elegans.

Recently, high-throughput next-generation sequencing techniques have allowed researchers to broadly explore the extent and complexity of the transcriptomes of a wide range of eukaryotic species and to gain novel information about the gene responses of aquatic species to physiological and environmental stresses (Wang et al., 2009). RNA-seq is an efficient technique to probe the gene responses to physiological stress (Li et al., 2013; Smith, Bernatchez & Beheregaray, 2013; Xia et al., 2013), particularly when working with species that do not have a sequenced genome. The present study used RNA-seq to analyze the transcriptomic response of T. s. elegans to salinity stress and identify the genes (and their metabolic functions) that are involved in salinity adaptation of turtles challenged by a brackish water environment. These results provide insights into the molecular mechanisms underlying osmoregulation in T. s. elegans and address the potential for this species to further invade and spread through new aquatic and brackish territories.

Materials and Methods

Animals

Healthy T. s. elegans were obtained from a local turtle farm in Hainan Province, China and were acclimated in three cement pools half-filled with freshwater for two weeks. After acclimatization, nine healthy T. s. elegans (BW: 424–478 g, 2 years old) were divided into three groups in pools (190 cm ×65 cm ×32 cm) of differing salinity: one in freshwater serving as the control (0 practical salinity units, psu), and the other two challenged with 5‰  (5 psu) or 15‰ (15 psu) saltwater. Turtles were fed a commercial diet each Monday and Thursday and 24 h after feeding, unused feed was siphoned out followed by replacement of one-third of the water in each pool. Water salinity was measured every day and adjusted to the proper salinity as needed. Other water quality parameters were monitored 2–3 times a week with steady values of pH 7.5–7.9, total ammonia nitrogen of <0.02 mg L−1 and temperature 26–28 °C. Photoperiod was 12 h:12 h L:D throughout.

After 30 days of exposure to the three experimental conditions, the three turtles from each group were subjected to a 24 h fast, then anesthetized by cryo-anesthesia moving turtles to −20 °C for 0.5–1 h. Experimental animal procedures had the prior approval of the Animal Research Ethics Committee of Hainan Provincial Education Centre for Ecology and Environment, Hainan Normal University (permit no. HNECEE-2014-004). Following euthanasia, the liver of each individual was sampled and divided into two sections, flash frozen in liquid nitrogen, and stored at −80 °C until used for RNA extraction. One liver section was used for quantitative real-time PCR (qRT- PCR), and the other was mixed from each group of three turtles for RNA-seq analysis.

Total RNA extraction, library construction and sequencing

Extraction of total RNA from liver samples used TRIzol® Reagent, following manufacturer’s instructions. Total RNA purity and concentration were determined using a NanoDrop 2000. The sequencing library was then constructed from high-quality RNA (OD260∕280 = 1.8–2.2, OD260∕230 ≥ 1.5, RIN ≥ 8.0, 28S:18 ≥ 1.0, >10 µg). RNA-seq analysis was provided by Novel Bioinformatics Co., Ltd using the Sanger/Illumina method. Subsequently, cDNA libraries were made using the Hiseq4000 Truseq SBS Kit v3-HS using 5 µg total RNA, following manufacturer’s instructions. Poly(A) mRNA was isolated with Dyabeads (Life Technologies, USA), fragmented with RNaseIII and purified. The fragmented RNA was added and ligated with ion adaptor. Then double-stranded cDNA was synthesized and purified using magnetic beads. The molar concentration of the purified cDNA in each cDNA library was then quantified with a TBS-380 fluorometer using Picogreen. The paired-end RNA-seq library was sequenced with an Illumina HiSeq 4000. The RNA-Seq data were deposited in the NCBI with accession number GSE117354.

De novo assembly and annotation

Raw reads were trimmed and quality controlled using SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) using default parameters. High-quality trimmed sequences were used for sequence assembly with Trinity (https://github.com/trinityrnaseq/trinityrnaseq/wiki) (Grabherr et al., 2011). Q20, Q30, GC-content and sequence duplication level of the clean data were all calculated. After that, All assembled transcripts were identified by using BLASTX against the databases of NR (NCBI non-redundant protein sequence), Swissprot (a manually annotated and reviewed protein sequence database), Pfam (Protein family), GO (Gene ontology), COG (Clutsters of Orthologous Groups of proteins) and KEGG (Kyoto Encyclopaedia of Genes and Genomes). The BLAST2GO (http://www.blast2go.com/b2ghome) program (Conesa et al., 2005) was used to obtain gene ontology (GO) annotations of unique assembled transcripts for describing biological processes, cellular components, and molecular functions.

Analysis of differential expression and functional enrichment

Expression levels of transcripts were calculated as fragments per kilobase of exon per million mapped reads (FPKM). RNA-Seq by Expectation-Maximization (RSEM; http://deweylab.biostat.wisc.edu/rsem/) was used to quantify gene transcripts, and DEGseq (http://www.bioconductor.org/packages/release/bioc/vignettes/DEGseq/inst/doc/DEGseq.pdf) was used to conduct differential expression analysis. 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 <0.05 found by DEGseq were assigned as differentially expressed. GO enrichment analysis of the DEGs was implemented by the GOseqR packages based on Wallenius non-central hyper-geometric distribution (Xie et al., 2011), which can adjust for gene length bias in DEGs.

Experimental validation by qRT-PCR

Eighteen genes identified as significantly expressed from the GO terms related to osmotic regulation were selected to understand the gene expression levels in different groups, and also used for validation by qRT-PCR. Table 1 shows the specific primers used. TRIzol® Reagent (Invitrogen, Carlsbad, USA) was used to extract total RNA from liver, followed by reverse-transcription using First-strand cDNA Synthesis Kit (Invitrogen, Carlsbad, CA, US). The qRT-PCR for gene expression was analyzed by an Applied Biosystems 7500 Fast Real-Time PCR System in 96-well plates with a 20 µl reaction volume containing 1 × SYBR Green qPCR Master Mix, gene-specific forward and reverse primers (0.4 µM) and cDNA (8 ng). The cycling conditions were 95 °C for 2 min followed by 40 cycles of 95 °C for 5 s and 60 °C for 30 s. β-actin was chosen as the reference gene, as it was expressed in the three groups and unaffected under salinity stress, and relative fold changes were determined using Relative Expression Software Tool v.2009 based on the cycle threshold (Ct) values generated by qRT-PCR.

Table 1

Sequences of primers for qRT-PCR validation.

Gene	Forward primer (5′ to 3′)	Reserse primer (5′ to 3′)	Product length (bp)
COX3	TCACTTGAGCCCACCATAGC	AGAGCCGTACACACCATCAG	155
ATP6	CTTGATGCCCTCTTCCCGTG	TTCCTCGTTCTCCACAGCCT	145
CYP17A1	GATCGGCTTCGAGAGACACC	GGATCAGCAGAGGGGAGACA	110
APoE	GTCTGGAGCGGGCTTAGTAG	CATTCCCAGGTCTCCCACAG	117
CoQ10	GCGAGTGCTGGGCTACT	TGAGCCCTTTGCGGTAGG	119
SAA	TCTAGGCGCTGGGGATATGT	CCACTAATGCCATCCTGCCA	180
FADS6	GCTGCCATACAACGAGGACT	AGCCCTATGTCTTGCTGTCG	146
HK	TAAAGGCGTAACCAGGCTGC	AATCGCACGTCAGAGTCAGG	126
GCK	CGGGAACTGCTGAAATGCTC	GAATGTGAAGCCCAGAGGCA	104
GLUL	GTTGCCACACCAACTTCAGC	AAGCGCGGATATGGTACTGG	111
ASRGL	TAGCACCTGTTCCAGTGAGC	GCTGTGTTTGATGCAGGTCA	193
TAT	CATCCACCAGCGACTCCAAG	CATCCTGGTGCCAAGACCTG	137
ASS1	CGGGCTGTACCAGAAACCAT	GGGACCATCCTGTACCATGC	129
STK32	TCCAGTGCTAATGCCAGCTC	TGGAACACCCCTTCCTGGTT	172
SIK1	TGGTGTGGTGCTGTATGTCC	TCCACAACTAGCATCCGTCG	156
SIK2	GCTGGTCCTAGACCCATCCA	GAAGGCTCGTTCTCCTGTCC	123
INSRR	CCGAGTACCGTGATCTGCTC	GGCAGCTCCACATCTACCAC	127
STK33	AGGCAGTTTTGGGGTGGTAA	TAAGATGCTCACCTCCCGTTC	130
β-actin	GCACCCTGTGCTGCTTACA	CACAGTGTGGGTGACACCAT	190

The mRNA expression levels were expressed as mean ± standard error. Statistical analyses were conducted with SPSS 19.0. After testing the homogeneity of variance, statistical difference between treatments and controls were determined by one-way analysis of variance (ANOVA). LSD multiple comparison tests were carried out when the variances were homogeneous. Significant differences were set at p < 0.05.

Results

Analysis of sequenced data quality

Because we haven’t got the genome of T. s. elegans yet, a transcriptome was used as a reference to identify the differentially expressed genes induced by salinity exposure. Therefore, a mixed RNA pool from the samples of the control and salinity treatment groups was sequenced as the reference transcriptome. The data qualities from each sample were shown in Table 2. A total of 157.21 million reads and 23.58 billion bases were obtained from the liver transcriptome of T. s. elegans, including 50.68 million in 0 psu, 57.91 million in 5 psu group and 48.62 million in 15 psu group. After filtering low quality sequences by, trimming sequencing adapters/poly-N and removing poor quality reads, there were 152.52 million clean reads (97.02% of raw reads) were retained including 48.92 million in 0 psu, 56.50 million in 5 psu and 47.10 million in 15 psu. Subsequently, 205,138 unigenes with an average length of 620 bp and N50 of 964 bp were obtained by de novo assembly. The largest and smallest unigenes were 22,866 bp and 201 bp, respectively (Table 3).

Table 2

Summary of Illumina expressed short reads production and filtering.

	Groups	Total_Reads	Total_Bases	Error%	Q20%	Q30%	GC%
Raw data	0 psu	50681398	7602209700	0.0126	96.31	92.31	49.34
	5 psu	57912616	8686892400	0.0121	96.72	92.87	49.18
	15 psu	48616356	7292453400	0.0129	96.21	91.98	49.93
Clean data	0 psu	48918412	7157032240	0.0104	98.29	95.12	49.19
	5 psu	56498176	8281257833	0.0103	98.34	95.18	49.09
	15 psu	47101028	6881004183	0.0107	98.17	94.77	49.83

Table 3

Summary of de novo assembly results of Illumina sequence data.

Type	Unigene	Transcripts
Total sequence number	205,138	244,815
Total sequence base	127,206,404	185,942,194
Percent GC	46.09	46.65
Largest	22,866	22,866
Smallest	201	201
Average	620.1	759.52
N50	964	1,535
N90	249	268

Annotation and differential expression of genes

Putative functions of proteins encoded by the 205,138 genes were predicted by NR, Pfam, COG, Swissprot, GO and KEGG database. The results showed that there were 20,458 (9.97%), 20,421 (9.95%), 13,362 (6.51%), 10,673 (5.2%), 25,438 (12.36%), and 38,651 (18.84%) in Pfam, KEGG, GO, COG, Swissprot and NR databases, respectively (Table 4).

Table 4

Functional annotation of the Trachemys scripta elegans transcriptome.

Annotated database	Number of unigenes	Precent(%)
Pfam	20,458	9.97
KEGG	20,421	9.95
GO	13,362	6.51
COG	10,673	5.2
Swissprot	25,438	12.36
NR	38,651	18.84
Total	205,138	1

By GO annotation, the genes up-regulated or down-regulated were divided into three categories including biological process, cellular component and molecular function (Fig. 1). Among the category of biological process, the number of differentially expressed gene was higher in the GO terms of cellular process, metabolic process, single-organism process and biological regulation. Among the category of cellular component, the number of differentially expressed gene was higher in the GO terms of cell, cell part, macromolecular complex, and membrane. In relation to molecular function, the number of differentially expressed gene was higher in the GO terms of binding, catalytic activity, and transporter activity. Moreover, there was almost the same trend in the comparisons of 0 vs. 5 psu, 0 vs. 15 psu, and 5 vs. 15 psu.

Figure 1

GO categories in the comparison of 0 vs. 5 psu (A), 0 vs. 15 psu (B), and 5 vs. 15 psu (C).

All genes were divided into several functional groups within three categories: cellular component, molecular function, and biological process. The below x-axis indicated the number of genes in each category, while the above x-axis indicated the percentage of total genes in that category.

On the basis of criteria of two-fold or greater change and Q of p < 0.05, 3,393 unigenes were identified as significant differentially expressed genes (DEGs), including 1019 (445 up-regulated and 574 down-regulated) DEGs in 0 vs. 5 psu, 1,194 (526 up-regulated and 668 down-regulated) DEGs in 0 vs. 15 psu, and1180 (548 up-regulated and 632 down-regulated) DEGs in 5 vs. 15 psu. Many DEGs were classified into some dominant categories by GO enrichment analyses, including the macromolecule metabolic process (52 DEGs), ion transport (35 DEGs), ion transmembrane transport (24 DEGs) in biological process, and nucleic acid binding (37 DEGs), oxidoreductase activity (35 DEGs), transporter activity (33 DEGs), and transmembrane transporter activity (27 DEGs) in molecular function (Table 5).

Table 5

Summary of GO term enrichment results on ion-regulation and macromolecular metabolism in T. s. elegans under salinity stress.

GO ID	Description	Ratio in study	Ratio in pop	p-value (FDR)
GO:0006091	generation of precursor metabolites and energy	17/357	119/13,362	4.33E − 05
GO:0009060	aerobic respiration	7/357	17/13,362	0.0003
GO:0045333	cellular respiration	7/357	19/13,362	0.0004
GO:0015980	energy derivation by oxidation of organic compounds	8/357	40/13,362	0.0065
GO:0034220	ion transmembrane transport	24/357	336/13,362	0.0081
GO:0098660	inorganic ion transmembrane transport	18/357	209/13,362	0.0081
GO:0098655	cation transmembrane transport	18/357	210/13,362	0.0081
GO:0006119	oxidative phosphorylation	4/357	7/13,362	0.0081
GO:1902991	regulation of amyloid precursor protein catabolic process	3/357	3/13,362	0.0081
GO:1902992	negative regulation of amyloid precursor protein catabolic process	3/357	3/13,362	0.0081
GO:1902430	negative regulation of beta-amyloid formation	3/357	3/13,362	0.0081
GO:1902003	regulation of beta-amyloid formation	3/357	3/13,362	0.0081
GO:0006811	ion transport	35/357	623/13,362	0.0109
GO:0006812	cation transport	23/357	335/13,362	0.0114
GO:1902600	hydrogen ion transmembrane transport	10/357	77/13,362	0.0117
GO:1900221	regulation of beta-amyloid clearance	3/357	4/13,362	0.0177
GO:0098662	inorganic cation transmembrane transport	15/357	176/13,362	0.0178
GO:0050773	regulation of dendrite development	6/357	28/13,362	0.0184
GO:0015672	monovalent inorganic cation transport	14/357	171/13,362	0.0410
GO:0005576	extracellular region	18/357	220/13,362	0.0097
GO:0044463	cell projection part	16/357	206/13,362	0.0281
GO:0015078	hydrogen ion transmembrane transporter activity	11/357	76/13,362	0.0047
GO:0015077	monovalent inorganic cation transmembrane transporter activity	16/357	176/13,362	0.0083
GO:0022890	inorganic cation transmembrane transporter activity	18/357	236/13,362	0.0167
GO:0008324	cation transmembrane transporter activity	20/357	293/13,362	0.0268
GO:1902991	regulation of amyloid precursor protein catabolic process	3/357	3/13,362	0.0081
GO:1902992	negative regulation of amyloid precursor protein catabolic process	3/357	3/13,362	0.0081
GO:0034364	high-density lipoprotein particle	5/357	13/13,362	0.0081
GO:0032994	protein-lipid complex	5/357	16/13,362	0.0129
GO:1990777	lipoprotein particle	5/357	16/13,362	0.0129
GO:0034358	plasma lipoprotein particle	5/357	16/13,362	0.0129

Notes.

GO names were retained only from GO terms of levels >2.

As for osmoregulation, there were 15 major GO terms related to osmotic function in turtle liver in the comparison of 0 vs. 5 psu (Table S1) and 12 in the comparison of 0 vs. 15 psu (Table S2). Based on GO enrichment analysis, a list of candidate genes involved in salinity acclimation in T. s. elegans were identified. Many of these were identified as encoding proteins involved in ion transport, energy production and conversion, and macromolecule metabolic process including lipid, protein and carbohydrate.

Based on the annotation of DEGs, we selected ten genes related to ion regulation, five genes related to energy production and conversion, five genes related lipid metabolism, eight genes related amino acid metabolism and six genes related carbohydrate metabolism (Table 6). Among the DEGs associated with ion regulation, five genes (adipocytokine, ACDC; insulin receptor-related protein, INSRR; serine/threonine-protein kinase 32, STK32; salt-inducible kinase 1, SIK1; potassium voltage-gated channel subfamily H member 5, KCNH5) were up-regulated in the comparison of 0 vs. 5 psu, and four genes (salt-inducible kinase 2, SIK2; SIK1; ACDC; STK32) were up-regulated in the comparison of 0 vs. 15 psu, and five genes (SIK1, SIK2, STK33, ACDC, and solute carrier family 26 member 9, SLC26A9) were up-regulated in the comparison of 5 vs. 15 psu. All of the five DEGs (cytochrome c oxidase subunit I, COX1; cytochrome c oxidase subunit III,COX3; cytochrome b, CYTb; F-ATPase protein 6, ATP6, and cytochrome P450 17A1, CYP17A1) associated with energy production and conversion showed up-regulation under salinity stress, which indicated that there was a higher energy demand in response to salinity exposure. In general, lipid and carbohydrate were the main sources of energy. As for five DEGs associated with lipid metabolism, three DEGs (apolipoprotein E precursor, ApoE; coenzyme Q-binding protein, CoQ10; high-density lipoprotein particle, SAA) increased with ambient salinity increased, while one DEGs (alcohol dehydrogenase 4, ADH4) showed down-regulation in the comparison of 0 vs. 5 psu and one DEGs (fatty acid desaturase 6, FADS6) down-regulation in the comparison of 0 vs. 15 psu. Among six DEGs related to carbohydrate metabolism, two DEGs (hexokinase, HK and lens fiber major intrinsic protein, MIP) increased with ambient salinity increased, four DEGs (glucokinase, GCK; tagatose 1,6-diphosphate aldolase, LacD; L-gulono-gamma-lactone oxidase, GLO and ribulose bisphosphate carboxylase small chain, RBCs) decreased under salinity stress. In addition, six DEGs (ornithine decarboxylase antizyme 3, OAZ3; glutamine synthetase, GLUL; asparaginase-like protein, ASRGL; L-amino-acid oxidase-like, LAAO; sodium-dependent neutral amino acid transporter, SLC6A15s; amino acid permease, SLC7A9) related to amino acid metabolism showed up-regulation and two DEGs (tyrosine aminotransferase, TAT and argininosuccinate synthase, ASS1) showed down-regulation in the comparison of 0 vs. 5 psu, while four DEGs (OAZ3, GLUL, ASRGL, SLC7A9) showed up-regulation and three DEGs (LAAO, TAT, ASS1) showed down-regulation in the comparison of 0 vs. 15 psu.

Table 6

DEGs related to energy production and conversion, macromolecule metabolic process, and ion transport in the liver of T. s. elegans.

Unigene ID	Description	Log2 (5/0 psu)	Log2 (15/0 psu)	Log2 (15/5 psu)
DEGs related to energy production and conversion
c198757_g1	cytochrome c oxidase subunit I (COX1)	7.03	7.14	–
c205807_g1	cytochrome c oxidase subunit III (COX3)	5.79	5.57	–
c184332_g1	cytochrome b (CYTb)	6.37	6.51	–
c168776_g1	F-ATPase protein 6 (ATP6)	7.93	8.41	0.67
c122818_g1	cytochrome P450 17A1 (CYP17A1)	0.81	1.63	0.82
DEGs related to lipid metabolic process
c183562_g1	apolipoprotein E precursor (APoE)	4.59	4.18	−0.40
c108401_g1	coenzyme Q-binding protein (CoQ10)	2.51	3.44	0.93
c106996_g1	high-density lipoprotein particle (SAA)	1.24	4.13	2.89
c101752_g1	alcohol dehydrogenase 4 (ADH4)	−3.8	–	3.81
c5275_g1	fatty acid desaturase 6 (FADS6)	−3.38	−1.54	1.85
DEGs related to carbohydrate metabolic process
c184502_g1	hexokinase-1 (HK)	4.49	5.22	0.73
c141564_g1	lens fiber major intrinsic protein (MIP)	3.37	3.75	–
c121000_g1	glucokinase (GCK)	−0.22	−1.94	−1.85
c198426_g1	tagatose 1,6-diphosphate aldolase (LacD)	−5.38	−5.12	–
c139209_g1	L-gulono-gamma-lactone oxidase (GLO)	−5.38	−5.12	–
c44501_g1	ribulose bisphosphate carboxylase small chain (RBCs)	−4.38	−4.12	–
DEGs related to amino acid metabolic process
c198734_g1	ornithine decarboxylase antizyme 3 (OAZ3)	6.28	7.00	0.72
c184545_g1	glutamine synthetase (GLUL)	4.03	4.3	0.26
c197142_g1	asparaginase-like protein 1b (ASRGL)	4.49	4.96	0.47
c117856_g3	L-amino-acid oxidase-like (LAAO)	1.62	−2.06	−3.68
c169209_g1	sodium-dependent neutral amino acid transporter B (SLC6A15s)	3.59	–	−1.54
c103361_g1	amino acid permease (SLC7A9)	1.78	5.72	3.93
c108456_g1	tyrosine aminotransferase (TAT)	−1.72	−1.41	0.31
c99414_g1	argininosuccinate synthase (ASS1)	−0.69	−1.1	−0.41
DEGs related to ion transport
c182997_g1	potassium voltage-gated channel subfamily H member 5 (KCNH5)	3.5	–	−3.5
c121057_g1	serine/threonine-protein kinase 32 (STK32)	1.95	1.80	–
c121806_g3	salt-inducible kinase 1 (SIK1)	1.23	2.20	0.97
c123793_g3	salt-inducible kinase 2 (SIK2)	–	0.27	0.40
c114797_g4	adiponectin (ACDC)	0.60	2.12	1.52
c119528_g1	insulin receptor-related protein (INSRR)	0.30	−1.70	−2.00
c95283_g1	serine/threonine-protein kinase 33 (STK33)	−1.5	−0.65	0.85
c58533_g1	sodium channel subunit beta-1 (SCN1B)	−5.31	−5.31	–
c94169_g1	natriuretic peptides A-like (NPPA)	−5.72	−5.72	–
c95302_g1	solute carrier family 26 member 9 (SLC26A9)	−3.55	−0.44	3.11

Notes.

The values above zero show up-regulation of gene expression, while the values below zero show down-regulation. “–” means that the level of gene expression is so low that it could not be detected.

Verification of gene expression by SYBR Green qRT-PCR

Some genes related to ion transport, energy production and conversion and macromolecule metabolic process mentioned above were selected for qRT-PCR analysis in order to validate the differentially expressed genes that were identified by RNA-Seq and gain detailed quantitative information on their differing expression patterns. As shown in Table 7, the mRNA expression levels of COX3, ATP6 and CYP17A1 related to energy production and conversion in the 5 psu group were almost 3.6-, 2.2- and 1.5-fold of the control respectively, while those in the 15 psu group were almost 2.4-, 3.4- and 2.1-fold of the control respectively. In relation to ion transport, the mRNA expression levels of STK32 and SIK1 in the 5 psu group were almost 3.6- and 2.1-fold of the control respectively, while those in the15 psu group were almost 3.0- and 3.5-fold of the control respectively. The mRNA expression levels of INSRR increased 2-fold in the 5 psu group and decreased 2.6-fold in the 15 psu group compared to the control, however, that of STK33 decreased 2.2- and 1.5-fold in the group of 5 psu and 15 psu compared to the control. As for macromolecular metabolism, the mRNA expression levels of APoE, CoQ10, SAA, HK, GLUL, and ASRGL increased with ambient salinity increased. Especially for the genes of CoQ10 and SAA, the mRNA expression levels in the 15 psu group were almost 8.4- and 15.3-fold of the control respectively. However, the mRNA expression levels of GCK and ASS1 decreased with ambient salinity increased, and those in the group of 5 psu were almost 1.5- and 2.0-fold of the control, and those in the group of 15 psu were 3.7- and 9.7-fold of the control. The mRNA expression levels of FADS6 and TAT in the 5 psu group were lowest among the three groups, and decreased almost 2.3- and 3.6-fold compared to the control.

Table 7

The expression levels of some genes in relation to osmotic adjustment in the liver of T. s. elegans by qRT-PCR.

Unigene ID	Description	Control	5 psu	15 psu
Genes related to energy production and conversion
c205807_g1	cytochrome c oxidase subunit III (COX3)	1.01 ±  0.06c	3.57 ± 0.27a	2.40 ± 0.29b
c168776_g1	F-ATPase protein 6 (ATP6)	1.09 ± 0.10c	2.42 ± 0.09b	3.67 ± 0.21a
c122818_g1	cytochrome P450 17A1 (CYP17A1)	2.13 ± 0.11c	3.15 ± 0.22b	4.46 ± 0.38a
Genes related to lipid metabolic process
c183562_g1	apolipoprotein E precursor (APoE)	0.86 ± 0.06c	2.77 ± 0.12a	1.80 ± 0.13b
c108401_g1	coenzyme Q-binding protein (CoQ10)	0.64 ± 0.08c	2.66 ± 0.17b	5.39 ± 0.23a
c106996_g1	high-density lipoprotein particle (SAA)	0.48 ± 0.06c	1.18 ± 0.08b	7.32 ± 0.27a
c5275_g1	fatty acid desaturase 6 (FADS6)	2.23 ± 0.27a	0.97 ± 0.23b	1.37 ± 0.13b
Genes related to carbohydrate metabolic process
c184502_g1	hexokinase-1 (HK)	1.37 ± 0.09c	2.35 ± 0.21b	3.27 ± 0.22a
c121000_g1	glucokinase (GCK)	1.79 ± 0.07a	1.19 ± 0.22b	0.48 ± 0.05c
Genes related to amino acid metabolic process
c184545_g1	glutamine synthetase (GLUL)	1.82 ± 0.14c	5.10 ± 0.23b	6.00 ± 0.19a
c197142_g1	asparaginase-like protein 1b (ASRGL)	1.01 ± 0.05b	1.64 ± 0.17a	1.94 ± 0.31a
c108456_g1	tyrosine aminotransferase (TAT)	3.37 ± 0.17a	0.93 ± 0.03c	1.38 ± 0.06b
c99414_g1	argininosuccinate synthase (ASS1)	5.16 ± 0.29a	2.59 ± 0.20b	0.53 ± 0.07c
Genes related to ion transport
c121806_g3	salt-inducible kinase 1 (SIK1)	1.19 ± 0.13c	2.46 ± 0.31b	4.16 ± 0.15a
c123793_g3	salt-inducible kinase 2 (SIK2)	1.06 ± 0.15a	0.85 ± 0.09a	1.13 ± 0.10a
c119528_g1	insulin receptor-related protein (INSRR)	2.04 ± 0.18b	4.13 ± 0.26a	0.78 ± 0.13c
c121057_g1	serine/threonine-protein kinase 32 (STK32)	1.48 ± 0.16c	5.39 ± 0.24a	4.49 ± 0.33b
c95283_g1	serine/threonine-protein kinase 33 (STK33)	3.93 ± 0.37a	1.79 ± 0.16b	2.57 ± 0.23b

Notes.

Different lowercase letters represent significance among different groups (p < 0.05).

The qRT-PCR results were significantly correlated with the RNA-seq results with correlation coefficients of 0.744 in 0 vs. 5 psu, 0.862 in 0 vs. 15 psu group and 0.748 in 5 vs. 15 psu group (p < 0.05) (Fig. S1). This provides strong evidence that RNA-Seq data can be a reliable indicator of the expression patterns of the hundreds of genes identified as differentially expressed in the current study.

Discussion

Osmoregulation in some aquatic animals can be a complex process because individuals must deal with fluctuating salinity levels in their natural habitats, often on a daily or seasonal basis. As a normally freshwater species, the red-eared slider T. s. elegans, does not possess salt glands, and requires osmoregulation to survive when entering environments of higher salinity (e.g., brackish or estuarine water). The physiological research has shown that T. s. elegans can increase blood osmotic pressure by balancing the entry of NaCl with the decreased secretion of aldosterone, and accumulating urea and free amino acids in blood (Hong et al., 2014). In our study, many DEGs were classified into some dominant categories including the macromolecule metabolic process, ion transport and ion transmembrane transport in biological process, which provide a strong evidence for the physiological mechanism.

Under salinity stress, loss of water from the cells can cause cells shrink and potentially die. In this sense, ion regulation is important for a cell to balance osmotic change. KCNH5 (http://www.ncbi.nlm.nih.gov/gene/27133) and SCN1B (Qin et al., 2003) are involved in potassium/sodium voltage-gated ion channels and fluid balance that controls arterial blood pressure by altering blood electrolyte composition and/or volume. Natriuretic peptide A (NPPA) is well known to regulate body fluid levels and electrolytic homeostasis and has natriuretic, diuretic, and vasodilatory actions (Espiner et al., 2014). NPPA is highly expressed and associated with H2O/Na+ absorption and protein Ser/Thr phosphatases (Espiner et al., 2014). SIK acts to modulate adrenocortical function particularly in response to high plasma Na+, K+, ACTH, or stress (Wang et al., 1999). SIK1 also has a role in steroidogenesis whereas SIK2 is implicated in gluconeogenesis regulation in liver; both belong to the AMPK (AMP-activated kinase) subfamily of serine/threonine kinases (Berggreen et al., 2012). The AMPK is a crucial regulator of cellular energy levels (Hardie & Ashford, 2014) and under stress conditions that deplete ATP, AMPK action promotes ATP-producing catabolic pathways while inhibiting ATP-consuming anabolic functions (Rider et al., 2009; Rider et al., 2006). Transcripts of adiponectin (ACDC) were also enriched in liver under both salinity stresses and this hormone participates in the pathway of fatty acid oxidation by regulating AMPK (Chong et al., 2013). So, in our study, five genes (ACDC, INSRR, STK32, SIK1, KCNH5) were up-regulated in the comparison of 0 vs. 5 psu, and three genes (ACDC, STK32, SIK2) were up-regulated in the comparison of 0 vs. 15 psu, which suggested that T. s. elegans can adapt itself into saline environment by increasing the expression levels of genes related to ion regulation.

Genes associated with transporting molecules related to metabolic processes were also modulated in association with an up-regulation of genes involved in ATP energy production. As salinity level increases, acclimation of T. s. elegans to elevated salinity conditions requires investment of energy to maintain their homeostatic balance (Hong et al., 2014). In our study, the DEGs related to energy production and conversion including ATP6, COX1, COXIII, CYTb, and CYP17A1 increased with ambient salinity increased, which suggested that a requirement for more energy by T. s. elegans during raised salinity conditions. The result is associated with the ATP needed for the synthesis and operation of transport-related proteins that drive ion- and osmoregulatory processes (Lee et al., 2003). It is also documented that changes in components related to the glycolysis, fatty acid metabolism, and ATP production are often associated with raised salinity conditions in freshwater fish (Lavado, Aparicio-Fabre & Schlenk, 2014; Tine et al., 2008). Our previous study has shown that T. s. elegans increased serum glucose and triglyceride levels when subjected to salinity stress (Shu et al., 2012). This study also indicated that salinity stress influences glycolysis/gluconeogenesis and fatty acid metabolism-related genes. The mRNA expressions of HK and MIP were increased, which suggested that utilization of glycogen might be increased and more glucose might be used for energy expenditure under salinity exposure. Moreover, the gene expressions of APoE, CoQ10, and SAA were up-regulated while ADH4 and FADS6 were down-regulated, which suggested that the lipolysis increased and lipogenesis decreased to produce more energy in response to salinity stress.

The ureogenesis is a strategy for diamondback terrapins (Malaclemys terrapin ) that are known to inhabit brackish water, and for the desert tortoise, Gopherus agassizii under dehydration stress (Dantzler & Schmidt-Nielsen, 1966). The red-eared slider T. s. elegans, can also increase urea content in plasma and tissues by synthesizing organic osmolytes (e.g., free amino acids and urea) in order to provide a colligative defense against water loss (Hong et al., 2014). Some DEGs associated with cellular amino acid metabolism including GLUL, ASRGL, LAAO, SLC6A15s and SLC7A9 were up-regulated in the comparison of 0 vs. 5 psu, which indicated that amino acid metabolism and transport were strengthened by salinity stress. The up-regulated expression of GLUL showed the increase of glutamine synthetase. Glutamine is the most prevalent amino acid in body fluids and muscle, is mainly transported by a Na+-dependent neutral amino acid system, and its turnover rate exceeds those of other amino acids (Zander et al., 2015). Glutamine is synthesized from the ATP-dependent conjugation of ammonia to glutamate, and is a well-known defense against ammonia (Cooper & Plum, 1987; Essexfraser et al., 2005). When salinity returns to normal, glutamine pools can be utilized as precursors for a variety of important cell molecules (e.g., purines, pyrimidines, mucopolysaccharides) or, in the presence of glutaminase, glutamine can be deaminated for direct excretion of ammonia in the kidney or used for urea synthesis in liver before excretion. The accumulation of urea in response to high salinity can be due to both urea retention and elevated rates of urea synthesis via the ornithine-urea cycle, uricolysis of uric acid, or hydrolysis of arginine (Dépêche & Schoffeniels, 1975; Gordon & Tucker, 1965). OAZ3, involved in the ornithine-urea cycle, showed up-regulation under salinity exposure in this study, which suggested that ureogenesis may also be activated in response to hyperosmotic conditions in this species. Our result correlated well with studies in crab-eating frog (Fejervarya cancrivora), where its resistance to hyperosmotic environmental conditions is generally linked to the accumulation of urea in plasma and tissues, urea resulting from an up-regulation of the hepatic urea synthesis machinery (Gordon & Tucker, 1965). Also, the striped catfish (Pangasianodon hypophthalmus) can counteract osmotic imbalance by triggering a regulatory volume increase, an internal process that initiates a net gain in osmolytes and/or water, increasing cell volume that re-establishes normal values and prevents further cell shrinkage (Nguyen et al., 2016).

The current study was designed not only to generate a catalogue of differentially expressed genes involved with salinity exposure, but also to allow the data to be integrated to identify the relationships on the adaptive response. Based on gene ontology information and from data in published studies, functional categories of differentially expressed genes were identified. Therefore, the potential interactions of the differentially expressed genes that responded to high salinity in T. s. elegans are shown in Fig. 2. The pathways and processes that are targeted provide us with numerous candidate genes for future investigations about the molecular mechanisms that underlie high salinity tolerance.

Figure 2

Interactions of positively selected genes and differentially expressed genes involved in the adaptation of T. s. elegans to high salinity.

(1) NPPA plays roles in the regulation of body fluid levels and electrolytic homeostasis pathway, while ACDC and SIK in the pathway of lipid (glucose) metabolism by the regulation of AMPK (STK32, STK33); (2) GLUL, TAT, LAAO, SCL6A9, SLC7A9 and OAZ3, ASRGL, ASS1 play roles in the accumulation of free amino acid and urea, while KCNH5, SCN1B, SLC26A9 mainly in the process of H2O/Na+ absorption; (3) ATP6, COX1, COX3, CYTb, CYP17A1 are associated with energy production and mediated by lipid and glucose metabolism (APoE, CoQ10, SAA, ADH4, FADS6, HK, MIP, GCK, LacD, GLO and RBCs).

Conclusions

In this study, we report the first transcriptome analysis of T. s. elegans under salinity stress. When T. s. elegans was subjected to salinity exposure, 3,393 unigenes in the liver were identified as DEGs, which were classified into some dominant categories including macromolecule metabolic process, ion transport, ion transmembrane transport in biological process. Also, a list of candidate DEGs potentially involved in salinity acclimation in T. s. elegans were identified into three kinds such as ion transport, energy production and conversion, and macromolecule metabolic process including lipid, protein and carbonhydrate.

The genes related to macromolecule metabolic process (OAZ3, GLUL, ASRGL, LAAO, SLC6A15s, SLC7A9, APoE, CoQ10, SAA, HK, and MIP), ion transport (KCNH5, STK32, SIK1, and ACDC) and ATP synthesis (ATP6, COX1, COX3, CYTb, and CYP17A1) were up-regulated by salinity stress, which indicated that T. s. elegans could adapt itself into salinity stress by balancing the entry of NaCl and accumulating urea and free amino acids in blood in response to osmotic pressure with higher ATP energy production. However, some genes related to macromolecule metabolic process (ADH4, FADS6, LacD, GLO, RBCs, TAT, and ASS1) and ion transport (SCN1B, NPPA, SLC26A9, and STK33) were down-regulated by salinity stress. Finally, we combined the data on functional salinity tolerance genes into a hypothetical schematic model that describe potential relationships and interactions among target genes to explain the molecular pathways related to salinity responses in T. s. elegans.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.