Transcriptome Analysis Reveals Impaired Fertility and Immunity Under Salinity Exposure in Juvenile Grass Carp

Grass carp (Ctenopharyngodon idellus) is one of the most economically important aquaculture species and is widely cultured in China. However, its wild populations in many rivers are increasingly declining, and seawater intrusion is one of the most important threats to their survival. However, the mechanisms underlying the decline due to salinity pressure are still unknown. Here, we performed a comparative transcriptome analysis of C. idellus larvae in response to salinity exposures; a total of 481 differentially expressed genes (DEGs) were identified. These DEGs were significantly enriched in eight Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, among which steroid biosynthesis was the most important one, with the highest enrichment score. The pathway plays an important role in the development of the testes and ovary. Interestingly, all DEGs in steroid biosynthesis showed a down regulation, indicating that salinity exposure may pose damage to the fertility of C. idellus. Furthermore, three immunity-associated pathways (cytokine–cytokine receptor interaction, Toll-like receptor signaling pathway, and NOD-like receptor signaling pathway) were also significantly enriched, suggesting impaired immunity and a high risk of disease infection under salinity exposure. Overall, damage to both fertility and immunity would decrease the number of offspring and increase the risk of death due to disease infection. Our results provide a potential molecular mechanism underlying the decline of wild C. idellus populations in the Pearl River.


INTRODUCTION
Seawater intrusion is a natural phenomenon along shores. However, seawater intrusion has intensified and frequently occurs downstream of rivers in recent years mainly due to global climate warming, which is accelerating the rise of sea levels (Melloul and Collin, 2006;Werner and Simmons, 2010;Werner et al., 2013). Seawater intrusion to inland river regions can increase the salinity concentrations and thus pose a threat to freshwater species (Nguyen et al., 2009). High salinity can affect physiological and biochemical reactions, such as fish growth and survival, by decreasing the digestive enzyme activity, immunity, and metabolic activity (Tong et al., 2007;Kang et al., 2013). Hence, it is an urgent need to understand the mechanisms underlying the negative influence of seawater intrusion on freshwater species.
Grass carp (Ctenopharyngodon idellus) is one of the most economically important aquaculture species due to its high value as food and is widely cultured in China. In 2019, grass carp yields account for 21.72% of the freshwater aquaculture products in China (Fishery Statistical Yearbook of China, 2020), highlighting its importance in freshwater aquaculture. Grass carp is naturally distributed in lakes, streams, rivers, and other water systems in the wild. However, its wild resource in many systems, such as the Pearl River, has declined rapidly in the past years (Fishery Statistical Yearbook of China, 2017. The decline not only poses a threat to the protection of biodiversity in freshwater realms but also has potential negative effects on the aquaculture industry of grass carp, given that wild populations can provide excellent germplasm resources to avoid inbreeding depression (Berthelot et al., 2014;Lu and Luo, 2020). Seawater intrusion is considered as one of the major environmental factors responsible for the decline of its wild populations (Melloul and Collin, 2006;Nguyen et al., 2009;Werner and Simmons, 2010;Kang et al., 2013). However, the mechanisms underlying the decline due to salinity pressure are still unknown.
RNA sequencing (RNA-seq) is a powerful technology and has been widely used for the identification of crucial genes targeted by external pressures in many aquaculture species, such as catfish (Ictalurus spp.) (Li et al., 2012), yellow catfish (Pelteobagrus fulvidraco) (Zou et al., 2015), Mandarin fish (Siniperca chuatsi) , top mouth culter (Culter alburnus) (Zhao et al., 2016), Gymnocypris namensis , and rainbow trout (Zhou et al., 2021). For example, RNA-seq was used to identify the differentially expressed genes (DEGs) of grass carp liver imposed by ammonia exposure. The results showed that the DEGs were significantly enriched in the MAPK, FOXO, AMPK, and NF-kB signaling pathways, cytokine-cytokine receptor interaction, and the apoptosis signaling pathway. These pathways are mainly involved in apoptosis, indicating that grass carp liver apoptosis was induced by ammonia exposure through apoptosis pathways (Jin et al., 2017). Moreover, the draft genome of grass carp (Hu and Chen, 2015) can provide a reference for read mapping, promoting accuracy of subsequent analysis.
In this study, we performed comparative transcriptome expression analyses using RNA-seq technology in grass carp after salinity exposure for 48 h to identify DEGs and characterize the function categories of these DEGs, aiming to understand the mechanisms underlying death due to salinity pressure in the wild. The study would allow us to comprehensively understand the molecular mechanisms behind the population decline of grass carp in the freshwater realm.

MATERIALS AND METHODS
The procedures used in this study were approved by the Pearl River Fisheries Research Institute Animal Care Advisory Committee.

Sample Collection and Acclimation
Grass carp (2.0 ± 0.4 g) were purchased from an aquaculture farm (Huashan Yueqiangfeng Aquaculture Farm, Huadu District, Guangzhou, China) in November 2019. Fish were transported to our laboratory using plastic bags specially designed for fish transportation, and they were acclimatized in plastic tanks (4 × 200 L) with freshwater for 1 week. During the period of acclimation, grass carps were fed twice a day with commercial diet (Tongwei, Chengdu, China). Feeding was stopped 24 h before the experiment. After acclimation, a total of 90 fish were randomly selected and evenly divided into three groups: the first group was used for tolerance test to confirm the salinity concentration used in this experiment; the second and third groups were used as the experimental group and the control group, respectively.

Acute Salinity Stress
In order to confirm the salinity concentration used in the experimental group, we carried out a tolerance test of grass carp in different salinity gradients, which were prepared with sea crystal (25 kg/bag, purchased from Shanghai Yuepeng, Shanghai, China). The results showed that the 50% lethal concentration (LC 50 ) for grass carp was 0.5 psu at 72 h. Hence, a salinity concentration of 0.5 psu was used for salinity exposure of grass carp. For the control group, there was only fresh water without sea crystal.
The experiment and control groups were kept at 26 • C, pH 7.9-8.3. After 48 h, five individuals from each group were randomly chosen and dissected; the gills were collected and transferred into 2-ml tubes with RNA-EZ Reagents D RNA-Be-Locker A. The collected gill samples were kept at 4 • C overnight and then transferred into -80 • C until the extraction of total RNA.

RNA Extraction, Library Construction, and Transcriptome Sequencing
The total RNA of each sample was extracted according to the instruction manual of the TRlzol Reagent (Life Technologies, Carlsbad, CA, United States). RNA concentration and purity were measured using NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE, United States). Library construction was performed following the manufacturer's recommendation in the NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, United States). For the library fragments, purification was performed with the AMPure XP system (Beckman Coulter, Beverly, MA, United States) and size selection was carried out with USER Enzyme (NEB, Ipswich, MA, United States). Adaptor ligation cDNA was yielded at 37 • C for 15 min followed by 5 min at 95 • C. Then, PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. After purification of the PCR products with the AMPure XP system, library quality was assessed on the Agilent Bioanalyzer 2100 system. As a last step, clustering of the indexed samples was performed according to the instruction manual in the cBot Cluster Generation System using the TruSeq PE Cluster Kit v4-cBot-HS (Illumina). The libraries were then sequenced on an Illumina platform and paired-end reads were generated.

Quality Filtration and Read Mapping
Raw data were firstly processed through in-house perl scripts. In this step, clean reads were obtained by removing the reads containing an adapter, ploy-N, and low-quality sequences (lowquality base Q < 10, but up to 50% of a whole read). At the same time, the Q20 (Ewing et al., 1998), Q30, GC content, and the sequence duplication level of the clean data were calculated. Hisat2 (Kim et al., 2015) was used to map the clean reads to the C. idellus genome (Wang et al., 2015). Only reads with a perfect match or with one mismatch were further analyzed and annotated based on the reference genome.

Identification of DEGs
Gene reads were normalized using FPKM (fragments per kilobase of exon per million fragments mapped) values by Cufflinks (Trapnell et al., 2010). DESeq2 (Anders and Huber, 2010) was used to identify the DEGs. DESeq2 provides statistical routines using a model based on a negative binomial distribution. The resulting p-values were adjusted using the Benjamini and Hochberg approach for controlling the false discovery rate. Genes with an adjusted p-value < 0.01 and an absolute value of log2 ratio ≥ 2 were identified as DEGs.

Gene Ontology and KEGG Pathway Enrichment Analysis of DEGs
To achieve standardized gene functional classification, Gene Ontology (GO) analysis based on DEGs was performed by the GOseq R packages with Wallenius' non-central hypergeometric distribution (Young et al., 2010), which can adjust for gene length bias in DEGs. KOBAS (Mao et al., 2005) was used to test the statistical enrichment of the DEGs in Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

Real-Time PCR Validation of Genes
A total of five key genes including three immune-associated DEGs (CXCL11, CCL20, and CCR6) from cytokine-cytokine receptor interaction and two chaperones (HSP90 and DNAJB1) from the NOD-like receptor signaling pathway and protein processing in endoplasmic reticulum were selected for validation using quantitative real-time PCR (qPCR), with the same samples as those for transcriptome analysis. qPCR was performed on a Roche LightCycler 96 real-time PCR system (SN:10328) according to the manufacturer's instruction. Each PCR reaction was carried out in a total volume of 12.5 µl containing 6.25 µl of the Biomarker 2X SYBR Green Fast qPCR Mix (Biomarker Technologies, Beijing, China), 2.5 µM of each primer, 2 µl of cDNA template, and 3.75 µl of RNase-free ddH 2 O. The PCR program was composed of a denaturation step of 3 min at 95 • C, followed by 45 cycles of 95 • C for 5 s, 60 • C for 34 s, and 72 • C for 10 s. To rule out potential contamination, the negative controls in each run were set up with ddH 2 O as templates. Melting curve analysis of the PCR products was performed at the end of each PCR reaction to confirm the specificity of the primers. β-actin  was used as an internal control to normalize the relative expression of mRNA with the 2 − Ct method (Pfaffl, 2001). The forward and reverse primers of the genes for qPCR are shown in Table 1. Statistical analysis was performed using GraphPad Prism (version 8.0.0 for Windows, GraphPad Software, San Diego, CA, United States) 1 .

Transcriptome Sequencing and DEGs
Transcriptome sequencing produced a total of 61.16 Gb clean data with an average of 6.12 Gb clean data for each sample ( Table 2). The clean reads were firstly mapped to the reference genome of C. idellus (Wang et al., 2015). About 87.20% (± 0.62%) reads have only one mapping position, while 2.76% (± 0.13%) reads were mapped multiply ( Table 2). The transcriptome sequencing data from this study were deposited in the NCBI Sequence Read Archive (SRA) with the accession number PRJNA725283. To understand the molecular mechanisms underlying the damage to grass carp caused by sea intrusion, we performed comparative transcriptome expression analyses to identify the DEGs between the salinity stress group and the control group. Compared to the control group, we identified a total of 481 DEGs, among which 273 DEGs showed down regulation and 208 DEGs were upregulated.

GO Analysis of the DEGs
In order to understand the functions of the DEGs, we firstly performed GO analysis (Figure 1). The results showed that these DEGs were divided into three categories: biological process, cellular component, and molecular function. The category enrichment bar plot based on the DEGs showed that there were significant GO terms (adjusted p < 0.05) in both biological process and molecular function, while no significant GO terms were detected in cellular component, highlighting the importance of both biological process and molecular function in response to salinity exposure (Figure 2). In the biological process category, response to lipopolysaccharide, inflammatory response, immune response, and receptor-mediated endocytosis were the most represented. Additionally, most molecular function-related 1 www.graphpad.com   DEGs were associated with phosphoenolpyruvate carboxykinase activity and squalene monooxygenase activity.

KEGG Analysis of the DEGs
To identify the biological pathways that play an important role in response to salinity exposure, we performed KEGG enrichment analysis against the KEGG database. A total of eight pathways were significantly enriched (Figure 3), including one fertility-related pathway (steroid biosynthesis) with the highest enrichment score, three immunity-related pathways (cytokine-cytokine receptor interaction, Toll-like receptor signaling pathway, and NOD-like receptor signaling pathway), three metabolism-related pathways (fructose and mannose metabolism, purine metabolism, and pyruvate metabolism), and one ion exchange or signal transductionrelated pathway (gap junction). Specifically, all DEGs of steroid biosynthesis, cytokine-cytokine receptor interaction, and fructose and mannose metabolism had downregulated expressions (Supplementary Table 1).

Validation by Real-Time PCR
In order to validate the expression levels of the DEGs in transcriptome sequencing, we performed qPCR to measure the expression levels of five DEGs: heat shock protein 90 alpha (HSP90), C-X-C motif chemokine 11 (CXCL11), chemokine ligand 20 (CCL20), DnaJ homolog subfamily B member 1 (DNAJB1), and C-C chemokine receptor type 6 (CCR6). For all genes, the expression tendencies between transcriptome sequencing and qPCR showed good consistency (Figure 4).

DISCUSSION
Seawater intrusion is one of the most important threats to freshwater ecosystems. High salinity due to seawater intrusion can induce a change of osmotic pressure, affecting physiology processes and even leading to large-scale death of freshwater species. Here, we identified a total of eight KEGG pathways with high enrichment, among which the steroid biosynthesis pathway  had the highest enrichment score. In addition, GO analysis based on the DEGs in this study highlighted the importance of squalene monooxygenase, which participates in synthesis of sterol (Yamamoto and Bloch, 1970;Han et al., 2010). These results suggested that the synthesis of steroid was definitely influenced by salinity exposure. There were four DEGs in this pathway that were annotated to three proteins: lanosterol synthase, methylsterol monooxygenase, and squalene monooxygenase. The first one is widely detected in various organisms (Suzuki et al., 2006) and converts (S)-2,3-epoxysqualene to lanosterol in the biosynthesis of cholesterol (Wada et al., 2020), which is a precursor of steroid hormones and plays a role in maintaining Na + balance in salt-sensitive hypertension (Cuka et al., 2019). The second one is localized on the endoplasmic reticulum membrane and is believed to function in cholesterol biosynthesis (Li and Kaplan, 1996). The last one is a sub-pathway of the lanosterol biosynthesis pathway (Yalinkaya and Lai, 2000). Steroid had been shown to play an important role in the early development of germ cells in Drosophila (Morris and Spradling, 2012). The downregulated expressions of all these DEGs in this pathway indicated that steroid biosynthesis in larval grass carp was inhibited by salinity exposure, further prohibiting the normal development of spermatocytogenesis and oogenesis and ultimately posing damage to the reproductive capacity of grass carp. The decreased fertility would reduce the number of offspring and may become one of the important mechanisms underlying the population decline of wild grass carp in rivers, such as the Pearl River. Interestingly, three immunity-related pathways-cytokinecytokine receptor interaction, Toll-like receptor signaling pathway, and NOD-like receptor signaling pathway-were also significantly enriched based on the DEGs between the salinity exposure group and the control group. Furthermore, four GO terms (response to lipopolysaccharide, inflammatory response, immune response, and receptor-mediated endocytosis) represented in biological process were also immunity-related processes (Raetz and Whitfield, 2002;Dauphinee and Karsan, 2006), providing further evidence of the affected immunity caused by salinity exposure. Cytokine-cytokine receptor interaction is critical in determining the effects of inflammation in the development of diseases (Reeth, 2000) and associates with innate and adaptive inflammatory host defenses (Chen et al., 2000). There are eight DEGs with annotation of seven proteins in the cytokine-cytokine pathway (Supplementary Table 1). Among which, interleukin-8 (IL-8) is a type of neutrophilactivating cytokine that plays a role in inflammation and host defense (Baggiolini and Clark-Lewis, 1992). CXCL11 is associated with inflammatory and angiogenic processes (Liu et al., 2011). CCL20 is the unique chemokine interacting with CCR6 and is essential in promoting the growth of immature dendritic cells and lymphocytes (Schutyser et al., 2003). Another crucial protein, Cytokine receptor family member b8 (preferred name: IL20Rα), is part of the IL20R complex and is associated with both proand anti-inflammatory functions (Ha et al., 2020). In this study, two crucial DEGs (IL-8 and CXCL11) were upregulated and three DEGs (CCL20, CCR6, and IL20Rα) were downregulated in the cytokine-cytokine receptor interaction pathway, suggesting that high salinity may change the activation of the DEGs, further damaging the function of cytokine-cytokine receptor interaction. Toll-like receptors (TLRs), the main components of the second pathway, are transmembrane molecules with a leucine-rich repeat motif in the external domain. The structure of Toll-like receptors is traditionally conserved, attributing to their important role in host defense (O'Neill, 2000;Kimbrell and Beutler, 2001;Kaisho and Akira, 2002;Janeway and Medzhitov, 2002). The Toll-like receptor signal pathway could induce inflammatory cytokines in macrophages and dendritic cells, triggering innate and adaptive immune responses (Krieg, 2002;Ozato et al., 2002;Tsujimura et al., 2002). TLR25 is a major member of the teleost TLR1 subfamily, belonging to the group "non-mammalian TLRs" (Temperley et al., 2008;Wang et al., 2016), and activates the transcriptomes of some proinflammatory cytokines and contributes to immune response (Lee et al., 2020). The downregulation of TLR25 highlights the dysfunction of the Toll-like receptor signaling pathway. As a main component of the third immunity-related pathway, NODlike receptors are evolutionarily conserved intracellular proteins with at least 57 members and play key roles in host innate immune response (Franchi et al., 2010). They can be induced by both biotic factors such as bacteria (Swain et al., 2012a,b) and abiotic factors such as thermal stress (Basu et al., 2015). The observation suggested that the immunity of larval grass carp was largely affected by salinity exposure. Similar findings were also reported in other aquaculture species, such as Mauremys mutica (Gao et al., 2021), Apostichopus japonicus (Xu et al., 2016), Catla catla (Basu et al., 2015), and Labeo rohita (Swain et al., 2012a,b). It is widely recognized that impaired immunity may decrease the ability to deal with exogenous pathogens, making species more susceptible to illness or even mass mortality when populations suffer from adverse environmental challenges (Bowden et al., 2007). Taken together, impaired immunity may be another important mechanism responsible for the population decline of wild grass carp in rivers.
The maintenance of metabolism balance is a fundamental demand for the normal performance of physiological activities such as dealing with adverse environmental stress. Here, we identified one metabolism-related GO term (phosphoenolpyruvate carboxykinase activity) and three metabolism-related pathways (fructose and mannose metabolism, purine metabolism, and pyruvate metabolism). Pyruvate in the human body is produced by the degradation of glucose and glycogen, or the degradation of amino acids, and is involved in energy metabolism. In this study, three DEGs were detected in the pyruvate pathway and annotated to two proteins (phosphoenolpyruvate carboxykinase and mitochondrial trifunctional enzyme subunit beta). The former acts as the rate-limiting enzyme in gluconeogenesis, which is a process that converts various non-sugar substances into glucose or glycogen (Latorre et al., 2016). The latter is a protein of mitochondrial trifunctional enzyme that catalyzes the last part reactions of the mitochondrial beta-oxidation pathway, a major energy-producing process in tissues (Kamijo et al., 1994;Liang et al., 2018). Interestingly, the three genes in the pyruvate pathway were all downregulated. In addition, the expressions of the DEGs in both fructose and mannose metabolism showed similar patterns, suggesting that salinity exposure may inhibit the energy metabolism; thus, not enough energy was available for resisting environmental stress. Moreover, purines are essential components in energy-requiring enzymatic reactions, such as the metabolism of the nervous system (Jinnah et al., 2013). In this pathway, seven DEGs were detected and annotated to seven proteins, including guanylyl cyclases, 3 -phosphoadenosine 5 -prime-phosphosulfate, TWISTNB protein, ectonucleoside triphosphate diphosphohydrolase, NTPDase 3, NTPDase 8, cGMP-specific-3 ,5 -cyclic phosphodiesterase, and guanine deaminase. They are mainly involved in the regulation of the nervous system (Snyder et al., 2000;Xu et al., 2000;Belcher et al., 2006). In this study, most of the DEGs in this pathway were downregulated, suggesting that the function of the nervous system is damaged, which may lead to uncoordinated movements in animals. This may be the reason for the unusual performance of the up or down swimming movements observed in the present study. Hence, the imbalance of metabolism may be another reason for the mortality and population decline of grass carp in the wild.
It is traditionally thought that salinity exposure can result in the imbalance of osmotic pressure and loss of integrant water in organisms. Gap junctions, dynamic hexameric protein channels, play important roles in ion exchange and balance of osmotic pressure. Inorganic ions or small molecules such as Na + are necessary components in maintaining osmotic balance through gap junctions (Kumar and Gilula, 1996). In this study, there were four DEGs showing significant enrichment in the gap junction pathway. Among which, tubulin is one major component of microtubules that plays important roles in maintaining cell shape and transcription of cell signaling (Jordan and Wilson, 2004). The down regulation of tubulin 5 under salinity exposure indicated that the function of the gap junctions may be impaired, thus has potential to break the balance of osmotic pressure. The imbalance of osmotic pressure may cause tissue edema. Expectedly, intestinal edemas in the experiment group were observed, while no similar phenomenon was detected in the control group.
In conclusion, our results revealed that salinity exposure has a significant influence on gene expression in larval grass carp, and many DEGs were identified. The DEGs were mainly involved in fertility, immunity, metabolism, and osmotic pressure-related categories, which indicated that salinity exposure posed a significant influence on fertility, immunity, metabolism, and osmotic pressure. These changes may be the major mechanisms underlying the population decline of grass carp in river systems. Our study provided molecular evidence in support of the population decline of wild aquaculture species, which are conducive to the management of wild aquaculture species and recovery of germplasm resource.

DATA AVAILABILITY STATEMENT
The data presented in the study are deposited in the NCBI Sequence Read Archive (SRA) repository, accession number is PRJNA725283.

ETHICS STATEMENT
The animal study was reviewed and approved by the Pearl River Fisheries Research Institute Animal Care Advisory Committee.

AUTHOR CONTRIBUTIONS
JZ, JL, and XL contributed to the conception of the study. JZ, ZW, and YH performed the experiment. JZ, ZW, and JL contributed significantly to the analysis and manuscript preparation. JZ performed the data analyses and wrote the manuscript. JL and XL helped perform the analysis with constructive discussions. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We thank Yangchun Gao, Ph.D, from the Institute of Zoology, Guangdong Academy of Sciences, for editing the English text of a draft of this manuscript.