Transcriptome and molecular regulatory mechanisms analysis of gills in the black tiger shrimp Penaeus monodon under chronic low-salinity stress

Background: Salinity is one of the main influencing factors in the culture environment and is extremely important for the survival, growth, development and reproduction of aquatic animals. Methods: In this study, a comparative transcriptome analysis (maintained for 45 days in three different salinities, 30 psu (HC group), 18 psu (MC group) and 3 psu (LC group)) was performed by high-throughput sequencing of economically cultured Penaeus monodon. P. monodon gill tissues from each treatment were collected for RNA-seq analysis to identify potential genes and pathways in response to low salinity stress. Results: A total of 64,475 unigenes were annotated in this study. There were 1,140 upregulated genes and 1,531 downregulated genes observed in the LC vs. HC group and 1,000 upregulated genes and 1,062 downregulated genes observed in the MC vs. HC group. In the LC vs. HC group, 583 DEGs significantly mapped to 37 signaling pathways, such as the NOD-like receptor signaling pathway, Toll-like receptor signaling pathway, and PI3K-Akt signaling pathway; in the MC vs. HC group, 444 DEGs significantly mapped to 28 signaling pathways, such as the MAPK signaling pathway, Hippo signaling pathway and calcium signaling pathway. These pathways were significantly associated mainly with signal transduction, immunity and metabolism. Conclusions: These results suggest that low salinity stress may affect regulatory mechanisms such as metabolism, immunity, and signal transduction in addition to osmolarity in P. monodon. The greater the difference in salinity, the more significant the difference in genes. This study provides some guidance for understanding the low-salt domestication culture of P. monodon.


Introduction
The stress response is a complex physiological mechanism activated by stimuli that an organism perceives as potentially threatening (Vallejos-Vidal et al., 2022). It involves the sum of continuous and coordinated activation mechanisms at different functional levels, including molecular processes that have a direct impact on gene expression and protein synthesis, and will ultimately activate responses from the cellular to the systemic and performance levels (Parra et al., 2015). Salinity is one of the key factors regulating the distribution, abundance and diversity of aquatic animals and is one of the main environmental factors exerting selection pressure on aquatic animals (Sun et al., 2020). Changes in environmental salinity can directly affect the composition of aquatic animal fluids and the homeostasis of the internal environment (Chen et al., 2015). It has been shown that salinity stress can alter the diffusion of salts between the hemolymph and the environment, leading to altered cellular states. Within the critical point salinity range, significant changes in protein function and structure as well as in the physicochemical properties of NaCl solutions can occur (Normant et al., 2005). Low salinity alters the cell membrane potential, the optical density of solubilized proteins and the activity of various enzymes (Esparza-Leal et al., 2019). Prolonged hyposalinity exposure affects the physiology of the organism, reducing the immune response of aquatic animals to attack by external factors and leading to reduced survival (Chen and He, 2019). Under moderate hypoosmotic stress, aquatic animals typically meet the increased cost of osmoregulation by increasing their regular metabolic rate or reallocating energy to osmoregulatory processes (Rivera-Ingraham et al., 2016). Under extreme hypoosmotic stress, many aquatic animals are limited to maintaining basic functions to conserve energy reserves until the environment improves (Podbielski et al., 2022). Therefore, elucidating the adaptation mechanisms of aquatic animals to different salinities is crucial for the continued health of aquaculture.
Penaeus monodon (also known as black tiger shrimp) is one of the most commercially valuable species of aquaculture shrimp Li et al., 2020). P. monodon is popular with farmers and consumers due to its high economic and nutritional value (FAO, 2018). P. monodon has a wide range of salinity adaptations and can grow in a wide range of salinities from one psu~57 psu (Joseph and Philip, 2020). P. monodon has an isosmotic point of 750 mOsm·kg -1 (equivalent to 25 psu). In recent years, its culture has advanced significantly from coastal marine waters to inland low-salinity waters due to mariculture disease outbreaks and near-coastal pollution (Vezzulli, 2022). Therefore, understanding the mechanisms of P. monodon adaptation in long-term low-salt stress will be beneficial to its healthy culture. The regulatory mechanisms of various genes in salinity adaptation of P. monodon have been reported, including cold shock domain containing protein E1 (Si et al., 2022a), ecdysoneinduced protein E74 (Si et al., 2022b), calreticulin, calnexin and endoplasmic reticulum protein 57 (Visudtiphole et al., 2017). However, systematic studies on the adaptation mechanisms of P. monodon in long-term low-salt stress environments are still lacking.
With the development of high-throughput sequencing technologies, RNA-seq offers the opportunity to study the mechanisms of transcriptional variation in organisms under different environmental conditions (de Jong et al., 2019). Transcriptome sequencing not only allows the detection of almost all validated genes expressed in a specific cell or organ but also allows a comprehensive resolution of the regulatory mechanisms involved in a specific biological process based on the structure and function of the differential genes (Stupnikov et al., 2021). In addition, transcriptome sequencing allows the simultaneous analysis of all physiological processes, including metabolism (Xiao et al., 2020), proteostasis , osmoregulation (Gao et al., 2021) and other cellular processes (Ge et al., 2022). Thus, transcriptome technologies address key gene expression patterns in non-model organisms in specific environments and enable the detection of unknown genes and the discovery of novel transcripts. In recent years, transcriptome technologies have been widely used for transcriptome assembly and annotation in many aquatic animals, such as Macrobrachium rosenbergii (Liu X. et al., 2021), Coilia nasus (Gao et al., 2021) and Penaeus vannamei (Yu et al., 2022). Thus, RNA-seq is expected to make an important contribution to the resolution of long-term low salinity adaptation mechanisms in P. monodon.
The gill is an important respiratory organ of P. monodon and is involved in its respiratory metabolism, acid-base balance, ammonia waste excretion, ion transport and osmotic pressure regulation (Henry et al., 2012). Therefore, in this study, the gene expression profiles of P. monodon gill tissues reared at three salinity levels were examined by RNA-seq. The results of this study will provide new insights into the mechanisms of osmotic pressure regulation in P. monodon under chronic low salinity stress.

Ethical statement
Experimental shrimp were approved by the Animal Care and Use Committee of the Chinese Academy of Fisheries Sciences (CAFS) of the South China Sea Fisheries Research Institute of the Chinese Academy of Fisheries Science (NO: SCSFRI96-253). Our experiments were implemented in accordance with national and institutional guidelines for the protection and use of CAFS laboratory animals.

Experimental animals
Healthy P. monodon were obtained from the Shenzhen base of the South China Sea Fisheries Research Institute of the Chinese Academy of Fisheries Science. The experimental P. monodon were 7.8 ± 0.8 cm in length and 7.3 ± 1.0 g in body mass. The experiments were carried out in concrete ponds at a salinity of approximately 30 psu, 27-28°C and pH 7.6-7.8 for 1 week prior to the start of the experiment. During this period, the ponds were aerated and fed commercial bait four times a day. One-third of the water volume was replaced daily using multistage filtered and disinfected seawater.

Low salinity stress trials and sampling procedures
The experiment was divided into three groups: a pure seawater group with a salinity of 30 psu (high-salinity challenge, HC), a group with a salinity of 18 psu (moderatesalinity challenge, MC), and a group with a salinity of 3 psu (lowsalinity challenge, LC). The different salinities required for the experiments were adjusted by mixing bay water disinfected by multistage sand filtration with tap water in proportions, and the salinity was measured using a high-precision salinometer (ATAGO, Guangzhou). The formal experiments were carried out in special buckets for shrimp culture with a volume of 500 L and 400 L of culture water. Four parallels were set up in each group, using one bucket per parallel. Each bucket contained 60 juvenile shrimp. The salinities used were all 18 psu at the beginning and then adjusted by 3 psu each day until each group reached the desired salinity of 30 psu, 18 psu and 3 psu, respectively, before the official start of the culture trials at the three salinities. For 45 days of culture, water quality was tested daily, and shrimp growth and feeding were observed. Feeding was performed 3 times a day, weighing the same bait and feeding regularly. Water was changed completely every 3 days to maintain the water quality. At the end of the culture, gill tissue samples were taken from each of the three groups for transcriptome sequencing. Every three shrimp from the same group were mixed into one tube. All tissues were immediately frozen in liquid nitrogen and stored at −80°C until RNA extraction.

RNA extraction, library construction and sequencing
Total RNA was extracted from each sample of P. monodon gill tissue according to the instructions of the mirVana ™ miRNA ISOlation Kit (Ambion, America). The product (concentration ≥50 ng mL −1 ) was examined by 1.5% agarose gel electrophoresis for completeness, which showed clear 28S and 18S bands with no obvious dragging. RNA purity was measured by a NanoDrop 2000 (Thermo, America) with A260/280 ≈ 2.0. A total of 4 μg of total RNA was used to construct RNA-seq libraries according to the instructions of the TruSeq Stranded mRNA LTSample Prep Kit (Illumina, America) (Loureiro et al., 2019). The RNA libraries were tested by an Agilent 2100 microarray for length and quality. After passing quality control, all RNA libraries were sequenced using an Illumina HiSeqTM 2500 instrument (Illumina, United States). The raw reads were archived in the sequence read archive of the National Center for Biotechnology Information under accession number SRP262105.

De novo assembly and annotation
Transcript sequences were obtained using the paired-end splicing method with Trinity trinityrnaseq_r20131110 software (Grabherr et al., 2011), and the final unigene was obtained using TGICL 2.1 software (Pertea et al., 2003) for clustering and deredundancy extension to obtain the final unigene. Unigene sequences were aligned by blastx to non-redundant (NR), Swiss-Prot protein (Swiss-Prot), Clusters of orthologous groups for eukaryotic complete genomes (KOG), Kyoto encyclopedia of genes and genomes (KEGG) and gene ontology (GO) libraries for comparison. The annotations with e < 1e −5 were taken to obtain the protein with the highest sequence similarity to the given Unigene. The KEGG annotation information of UniGene was obtained by using KAAS (http://www.genome.jp/kaas-bin/kaas_main), and based on the annotation results of SWISSPROT, the GO term was mapped to obtain the GO protein function annotation information of the UniGene. The annotation summary table information was obtained by merging Unigene with the best one from each database.

Differentially expressed gene analysis and functional enrichment analysis
The expression levels of all DEGs were quantified based on the negative binomial distribution test in DESeq software (http:// bioconductor.org/packages/release/bioc/html/DESeq.html).
The NB (negative binomial distribution test) was used to test the significance of differences in the number of reads, and baseMean values were used to estimate the amount of gene expression.
The number of differential genes included in each GO entry was counted, and the significance of differential gene enrichment in each GO entry was calculated using the hypergeometric distribution test. A p-value <0.05 and log 2 -fold change value >2 indicated that differential genes were enriched in that GO entry.

Functional annotation and classification of P. monodon transcripts
A total of 64,475 unigenes were annotated in five commonly used databases, including NR, Swiss-Prot, KOG, KEGG and GO. E-values less than 1e −60 were present in 55.02% of the 20,990 unigenes (32.56%) annotated to the NR database. When mapped to other species, 10,780 unigenes (27.75%) showed a preference for sequence similarity to Hyalella azteca species. Sequences from 472 unigenes (2.25%) showed similarities to P. monodon species. Due to the shortage of decapoda genomic resources, up to 51.36% and 26.65% of unigenes could not be annotated in the NR and Swiss-Prot databases, respectively, and may represent the genetic specificity of P. monodon.

Simple sequence repeat (SSR) analysis
The 64,475 unigenes (total sequence length 82,251,807bp) were examined using MISA and Primmer 3 software, and their SSR loci were analyzed. We identified 51,544 SSR loci (31.70% of the total unigenes), of which 12,634 sequences contained more than one SSR. Among all SSRs, 5,387 SSRs were present in compound formation. Among them, the most represented repeat type was mononucleotide (46.04%), and the least represented repeat types were pentanucleotide and hexanucleotide (0.17%) ( Table 2).

FIGURE 1
The number of differentially expressed genes in the three groups.

Differentially expressed gene screening
The genes responsive to low salinity stress and the preliminary mechanism of P. monodon adaptation to low salinity stress were screened. The expression levels were calculated and compared between the three groups (HC, MC and LC groups), and the DEGs of the MC and LC groups compared to the HC group were determined (the conditions for significant differences were p-value <0.05 and difference multiplicity >2). The DEG results showed that compared to the HC group, 1,000 upregulated genes and 1,062 downregulated genes were observed in the MC group, and 1,140 upregulated genes and 1,531 downregulated genes were observed in the LC group compared to the HC group (Figure 1). Figure 2 shows the grouped clustering analysis of differentially expressed genes in three salinity groups, HC, MC and LC. The results showed that the differentially expressed genes basically had a clear upand downregulation relationship in the two samples. The genes thus classified are good references for conducting gene enrichment analysis and the resolution of salinity adaptation mechanisms.

GO and KEGG enrichment analysis of DEGs
Functional annotation analysis of DEGs was used to further compare and explore the transcriptome regulatory patterns under different salt concentration stresses. Figure 3 shows the GO enrichment analysis of all differentially expressed genes. The GO enrichment results showed the highest enrichment (p-value <0.05) in the three clusters of biological processes, molecular functions and cellular components. Figure 4 shows the integration of GO Level 2 classification of upregulated genes and downregulated genes in the differential grouping. The results showed that the differentially expressed genes were more annotated in the biological process category in cellular process, single-organism process, metabolic process, biological regulation, regulation of biological process and other functional localizations. The enriched cellular components were mainly functions such as cell, cell part, organelle, organelle part and membrane. The enrichment in molecular function was mainly in the functional classification of binding, catalytic activity, structural molecule activity, and transporter activity.
KEGG pathway analysis can provide detailed information on the pathways in which DEGs are involved. In this study, 1,223 DEGs were annotated into 36 pathways. These pathways were divided into six main groups at Level 1: Metabolism, genetic information processing, environmental information processing, cellular processes, biological systems, and human diseases ( Figure 5). There were significant differences in the enriched pathways between the experimental groups. In the HC vs. LC group, 583 DEGs significantly mapped to 37 signaling pathways. In the HC vs. MC group, 444 DEGs were significantly mapped to 28 signaling pathways. These pathways were mainly significantly associated with infectious diseases, signal transduction, immunity and metabolism. The top twenty pathways enriched in each group are shown in Figure 6. Frontiers in Physiology frontiersin.org

Metabolic and immune-related differentially expressed genes
To investigate the effects of low salinity on the metabolic regulation of P. monodon, gill tissues were screened for metabolism-related DEGs after 45 days of stress with different salinities. Twelve metabolic pathways were identified, involving processes in amino acid metabolism, glucose metabolism and lipid metabolism. Five significantly enriched metabolic pathways were detected in the LC vs. HC group, including pentose and  Frontiers in Physiology frontiersin.org 06 glucuronate interconversions, ascorbate and aldarate metabolism, cutin, suberin and wax biosynthesis, and sulfur metabolism. Three significantly enriched metabolic pathways were detected in the MC vs. HC group, including cutin, suberin and wax biosynthesis, and butanoate. Some representative differentially expressed genes are shown in Table 3.
To investigate the effect of chronic low salinity on immune regulation in P. monodon, the study analyzed immune system pathways. A total of 17 immune system regulatory signaling pathways were detected, including the NOD-like receptor signaling pathway, Toll-like receptor signaling pathway, and Toll and IMD signaling pathways. The study detected more differentially expressed genes in the NOD-like receptor signaling pathway and Toll and IMD signaling pathways. In addition, antigen processing and presentation were significantly different only in the LC vs. HC group, and cell differentiation was only different in the MC vs. HC group. Some representative differentially expressed genes are shown in Table 4.

Signal transduction-related differentially expressed genes
Understanding signal transduction pathways is an important guide for resolving the adaptation mechanisms of organisms to environmental changes. Therefore, the signal transduction pathways of P. monodon under chronic low salinity stress were analyzed during this study. A total of 31 signal transduction pathways were identified, including the PI3K-Akt signaling pathway, MAPK signaling pathway, AMPK signaling pathway, Hippo signaling pathway, calcium signaling pathway, sphingolipid signaling pathway, cGMP-PKG signaling pathway, HIF-1 signaling pathway, cAMP signaling pathway and Ras signaling pathway. Among them, the MAPK signaling pathway was significantly different in both the LC vs. HC and MC vs. HC groups, while the Hippo signaling pathway was only significantly different in the MC vs. HC groups. Some representative differentially expressed genes are shown in Table 5.

Discussion
Although previous studies have shown that P. monodon can tolerate a wide range of salinities, studies on its osmotic adaptation mechanisms under chronic low salinity stress are still lacking. The mechanisms of adaptation of P. monodon to long-term varying salinity were assessed by full-length transcriptome amplification in our study. This study identified 2,062 DEGs from the HC vs. MC group, including 1,000 upregulated and 1,062 downregulated genes, and identified 2,671 DEGs from the HC vs. LC group, including 1,140 upregulated and 1,531 downregulated genes. This difference in gene regulation suggests that although P. monodon can cope with this stress, 18 psu and 3 psu may cause subtle stress. The GO analysis results are consistent with previous studies on aquatic animals that classified DEGs into three categories of biological processes, molecular functions and cellular components, such as P. vannamei (Yu et al., 2022), Pelteobagrus fulvidraco  and Oratosquilla oratoria (Lou et al., 2019). Furthermore, KEGG analysis showed that changes in salinity significantly affected metabolic, immune and signaling pathways in P. monodon, consistent with the results of transcriptome analyses of P. vannamei (Chen et al., 2015), Acanthogobius ommaturus (Sun et al., 2020) and Eriocheir sinensis (Yang et al., 2019) under low salinity stress. Therefore, the present study focused on the adaptation of metabolic, immune and signal transduction mechanisms in relation to salinity. Salt stress is largely responsible for altering the metabolic state of organisms (Sun et al., 2020;Liu Z. et al., 2021). Our study found that the expression of metabolism-related genes was significantly affected and showed metabolic imbalance in P. monodon during prolonged low salinity stress. Many studies have reported similar results, but most of the studies on metabolism in osmoregulation have focused on osmotic physiological responses (Guo et al., 2020). In the present study, three significantly enriched metabolic pathways were detected in the MC vs. HC group, including Cutin, suberin and wax biosynthesis, Butanoate metabolism and Valine, leucine and isoleucine biosynthesis; five significantly enriched metabolic pathways were detected in the LC vs. HC group, including Pentose and glucuronate interconversions, Ascorbate and aldarate metabolism, Cutin, suberin and wax biosynthesis, Sulfur metabolism, Arginine biosynthesis and Phenylalanine metabolism. These pathways are mainly related to the regulation of glucose, lipid and amino acid metabolism. It has been found that crustaceans exposed to salinity fluctuations require additional energy to ensure osmotic pressure regulation and ion exchange (Chen et al., 2014). Lipids may provide sufficient energy to maintain ion homeostasis and regulate the permeability of membrane structures and therefore may play an important role in osmotic pressure regulation (Palacios et al., 2004;Tseng and Hwang, 2008). Free amino acids are osmotic effectors in crustaceans, and the metabolism of free amino acids may provide energy and eventually lead to cell growth in crustaceans (Lemos et al., 2001). Thus, metabolic pathways play an important role in the regulation of osmotic pressure in crustaceans.
Chronic low salinity stress can damage the immune system of aquatic animals. Antigen processing and presentation and the NODlike receptor signaling pathway are the immune pathways that are more affected by chronic low salinity stress in P. monodon. Antigen processing and presentation are the cornerstones of adaptive immunity, through which organisms immunologically recognize some bacteria, viruses and parasites (Pishesha et al., 2022). Therefore, changes in genes associated with this pathway suggest Frontiers in Physiology frontiersin.org that chronic low salinity stress may disrupt the innate immunity of shrimp and increase their susceptibility to bacteria, viruses, and parasites, among others (Kelly and Trowsdale, 2019). Complement and coagulation cascades are evolutionarily relevant blood enzyme cascade reactions in the blood. The complement system has immunoprotective and regulatory functions, while the coagulation cascades are responsible for ensuring hemostasis, and together, they limit pathogen infection through innate immune mechanisms. Thus, both protein C (PROC) and plasminogen (PLG), important regulators in the complement and coagulation cascades in this study, were significantly upregulated, suggesting that prolonged hyposalinity stress may have some effect on coagulation function in P. monodon (Ayón-Núñez et al., 2018). In addition, the NOD-like receptor signaling pathway (NLR) immune recognition factor was significantly altered in both groups. This regulation has important implications for NLR and complement components, which are pattern recognition receptors that recognize molecular patterns associated with viral, bacterial and eukaryotic pathogens (Ward and Rosenthal, 2014). For crustaceans, the innate immune system plays a crucial role in defense against microbial infections and environmental stress (Tran et al., 2022). These results suggest that P. monodon may have an increased susceptibility to pathogen attack in low salinity environments. The adaptive and stimulatory responses of aquatic animals to salt stress depend on effective mechanisms of osmotic sensing and osmotic stress signaling (Nie et al., 2017). Furthermore, osmotic sensing signals are not only rapidly transduced intracellularly but may also be amplified and distributed to multiple downstream osmotic effects (Sun et al., 2020). cAMP, as a signaling factor, is involved in the regulation of osmotic pressure by mediating the increased uptake of Na + by cells in the internal environment . At the same time, cAMP stimulates the production of calmodulin (CaM) and ATP, thereby regulating stress, apoptosis and calcium homeostasis .  Thus, as important ion transport regulatory proteins, CaM and NKA are significantly upregulated at low salinity. The PI3K-AKT pathway is one of the most actively studied kinase pathways, mediating cell growth, survival, cell cycle progression, differentiation, transcription, translation, and glucose metabolism in animals (Liu et al., 2018). Therefore, downregulation of this pathway may be responsible for the slow growth of P. monodon under prolonged hyposaline conditions, as shrimp require more energy for osmoregulation and ion homeostasis during normal growth (Cui et al., 2020). The MAPK pathway is important for osmotic pressure regulation, and chilling and hypoxia are important. In the present study, DEGs in the gill tissue of P. monodon under salt stress were enriched in the MAPK signaling pathway, indicating that the pathway was induced. Similar results were found in P. vannamei hemocytes under salt stress (Zhao et al., 2015). Furthermore, AMPK is suitable as an intracellular energy and pressure sensor to reveal the molecular mechanisms underlying temperature, hypoxia and salinity fluctuations in animals (Hardie, 2011). Salinity stress can temporarily activate AMPK, while glucose depletion leads to sustained AMPK activation. Moreover, in both cases, AMPK activation levels depend on the degree of stress (Schutt and Moseley, 2017). The activation of the AMPK signaling pathway and the expression of related genes in P. monodon under low salinity stress may be due to its continuous activation by glucose depletion. The differences in the expression of genes related to the AMPK signaling pathway in gill tissues of P. monodon at different salinities may be due to different salt stress stresses, resulting in different energy regulation mechanisms.

Conclusion
In this study, RNA-seq analysis of P. monodon exposed to three salinity gradients (namely, HC, MC and LC) was carried out, and the high-quality data obtained are the best guarantee of experimental reliability. These changes were verified by this experiment to be related to metabolic, immune system and signal transduction properties. Thus, P. monodon may adapt to changes in environmental salinity by regulating mechanisms such as metabolism, immunity and signal transduction. These data will provide basic data for future selection and breeding for resistance in the P. monodon industry, which will contribute to the sustainable and healthy development of shrimp farming.

Ethics statement
The animal study was reviewed and approved by Animal Care and Use Committee of South China Sea Fisheries Research Institute, Chinese Academy of Fishery Sciences (No. SCSFRI96-253).