Transcriptome Analysis of Gill Tissues from Neptunea cumingii in Different Seasons

1. Introduction

Neptunea cumingii (Gastropoda, Probranchia, Neptunea, Mothidae) is a carnivorous snail and the most economically important shellfish along the northern coast of China [31]. The meat of Neptunea cumingii is delicious, nutritious, and an ideal source of dietary protein and fat [11], and it is greatly favored by consumers. Previous studies of Neptunea cumingii mainly focuses on its distribution [9], growth and reproduction [8,59], genetic diversity, and the effects of salinity and dissolved oxygen on its physiological activities [60,61], whereas few have considered the mechanism associated with temperature adaptation by Neptunea cumingii.

Temperature is an important environmental factor that affects the metabolism, growth, and life activities of marine bivalves [17]. When shellfish perceive short-term fluctuations in temperature or seasonal changes, they will adapt to the change in the environmental temperature by actively adjusting their physiological and biochemical activities [41]. Previous studies have investigated the physical response mechanisms in various shellfish under changing environmental temperatures and shown that the expression of heat shock proteins (HSPs) is stimulated in shellfish to maintain cell homeostasis under high temperature conditions, where the response to damage also involves the activation of immune cells and accelerated metabolism. Under low temperature conditions, shellfish respond to oxidative stress by stimulating the production of glutathione S-transferases (GSTs), such as in Ruditapes philippinarum [17], Patinopecten yessoensis [16,28], and Chlamys farreri [28]. However, the regulatory mechanism that allows Neptunea cumingii to adapt to seasonal temperature changes has not been reported previously. Therefore, in this study, we conducted transcriptome sequencing using gill tissues collected from Neptunea cumingii in different seasons with natural temperature variability. Genes related to environmental adaptation were screened and the functions of these genes were verified in order to understand the mechanisms that allow Neptunea cumingii to adapt to seasonal temperature changes, and thus provide a reference to improve the production of Neptunea cumingii.

2. Materials and Methods

2.1. Experimental materials

In this study, samples of Neptunea cumingii were collected from Zhangzidao Island in Dalian City (ZZ; longitude: 121.66°N, latitude: 39.054°E) and Bajiao Port in Yantai City (YT, longitude: 121.15°N, latitude: 37.641°E), China. In particular, samples were collected during 2021 in the spring (March 28, ZZ temperature: 7°C, YT temperature: 8°C), summer (June 25, ZZ temperature: 21°C, YT temperature: 23°C), autumn (September 23, ZZ temperature: 18°C, YT temperature: 19°C), and winter (December 7, ZZ temperature: 5°C, YT temperature: 6°C). Twenty samples collected from each site were cleaned and dissected. The gill tissues from each individual were frozen in liquid nitrogen and stored in a freezer at –80°C for further analysis. During sample collection, five different samples from each sampling site were mixed in equal amounts.

2.2. Experimental methods

2.2.1. RNA extraction, library preparation, and Illumina sequencing

The total amounts and integrity of RNA were assessed using an RNA Nano 6000 Assay Kit with the Bioanalyzer 2100 system (Agilent Technologies, CA, USA) according to the manufacturer’s protocol. Briefly, mRNA was purified from total RNA by using poly-T oligo-attached magnetic beads. First strand cDNA was synthesized using random hexamer primers and M-MuLV reverse transcriptase, before applying RNaseH to degrade the RNA. Second strand cDNA synthesis was then performed using DNA polymerase I and dNTPs. The library fragments were purified with the AMPure XP system (Beckman Coulter, CA, USA). After PCR amplification, the PCR products were purified by using AMPure XP beads to obtain the final library. Then these libraries were sequenced on the illumina NovaSeq 6000（illumina，USA） and 150bp paired-end reads were generated. Transcriptome sequencing and analysis was performed by Novogene Biotechnology Co., Ltd. (Beijing, China).

2.2.2. Data filtering and de novo assembly

The image data measured by the high-throughput sequencer are converted into sequence data (reads) by CASAVA base recognition. Raw data (raw reads) of fastq format were firstly processed through in-house perl scripts. The reads containing ploy-N and the low-quality reads were removed to obtain the clean reads. After removing adaptor and low-quality sequences, theclean reads were assembled into expressed sequence tag clusters (contigs) and de novo assembled into transcript by using Trinity(version: 2.6.6) in paired-end method. Using the spliced transcriptome as a reference sequence (Ref), RSEM software [20]was used to compare the clean reads of each sample to the Ref. reads with a comparison mass value lower than 10 and those without paired comparison were filtered out, respectively. reads were compared to multiple regions of the genome, and then the readcount numbers were converted by Fragments Per Kilobase Million（FPKM）. Gene function was annotated based on the following databases: NCBI non-redundant protein sequences (Nr), NCBI non-redundant nucleotide sequences (Nt), Protein family (Pfam), Clusters of Orthologous Groups of proteins (KOG/COG), A manually annotated and reviewed protein sequence database (Swiss-Prot), KEGG Ortholog database (KO), Gene Ontology (GO).

2.2.3. DEGs screening and enrichment analysis

DEGs were identified using the DESeq2 R package (1.20.0) functions Benjamini and Hochberg. |log2(FoldChange)| > 1&P-value（padj） < 0.05 was set as the threshold for significantly differential expression. Hierarchical cluster analysis of DEGs was performed to explore transcripts expression pattern. Based on the hypergeometric distribution principle, GOseq (1.10.0) and KOBAS (v2.0.12) software were used for GO functional enrichment analysis and KEGG pathway enrichment analysis of DEGs.

2.2.4. Quantitative real time-PCR validation

To examine the reliability of the RNA-Seq results, we randomly selected 11 DEGs from 34 temperature-related DEGs for qRT-PCR validation，including caspase activity and apoptosis inhibitor 1（CAAP1）, toll-like receptor 8 （TLR8）, toll-like receptor 13（TLR13）, toll-like receptor 3（TLR3）, heat shock protein 70（HSP70）, heat shock protein 90（HSP90）, calreticulin （CRT） , phosphatidylinositol 3,4,5-trisphosphate 5-phosphatase 2（Ptp2）, glutathione S-transferase kappa 1-like （GSTK1）, interferon regulatory factor 5（ IRF5）, Mitogen-activated protein kinase kinase 6 (MAP2K6).The RNA was reverse transcribed to cDNA using TIANGEN^® FastKing gDNA Dispelling RT SuperMix (Tiangen Biotech, China). Suitable primers were designed using Primer 5.0 (Table 1) and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). qRT-PCR was performed with TIANGEN^® Tal-ent qPCR PreMix (SYBR Green) (Tiangen Biotech, China) on a Roche LightCycler^®96 (Roche Diagnostics, Switzerland) real-time PCR system according to the manufacturer's protocol. Samples from the same batch experiment of the RNA-seq were used for confirmation experiment. All reactions were performed in triplicate. The qRT-PCR conditions were as follows: 3 min at 95℃, followed by 45 cycles of 5 s at 95℃, and 15 s at 60℃. Dissociation curve analysis was carried out to determine the target specificity. The relative expression of tar-get genes was calculated using the 2^-ΔΔCT method [26].

3. Results

3.1. Transcriptome sequencing analysis

Gill tissue transcriptomes were sequenced by Illumina high-throughput sequencing to obtain 219.71G raw bases, and 206.28G clean bases (Table 2) were obtained after filtering, which comprised 93.89% of the original data. The average Q20 and Q30 scores were highest in the spring in ZZ, thereby indicating that the sequence quality was the greatest in this period. The average Q20 and Q30 scores for the other samples ranged among 95.67–97.55% and 90.11–93.41%, respectively. The GC contents of the samples were 37–46%, which were normal, and thus suitable for subsequent data analysis. The raw data have been submitted in the Sequence Read Archive (SRA) of NCBI with accession number PRJNA1023007.

3.2. Transcription splicing, comparison, and annotation

After removing adapters and low-quality reads using Trinity software, 187640 unigenes were obtained, with a longest length of 30,252 bp and shortest length of 301 bp. Among them, 80055 unigenes (42.66%) were less than 500 bp in length, 62261 (33.18%) were 500–1000 bp, 30404 (16.20%) were 1000–2000 bp, and 14,920 (7.95%) were greater than 2000 bp (Table 3). The clean data were compared with the unigenes as reference sequences, and the gene expression levels in different seasons were mapped, as shown in Figure 1.

The unigene annotations are shown in Table 4, where 73085 unigenes were successfully annotated, which accounted for 38.94% of the total unigenes. The numbers of unigenes annotated using the Pfam and GO databases were 51,868 and 51,864, respectively, and they accounted for 27.64% of the total. The number of genes annotated using the Nr database was 38,761, which accounted for 20.65% of the total. After GO enrichment analysis, the unigenes were assigned to three categories comprising Biological Process, Cell Component, and Molecular Function, and 42 GO terms, including environmental adaptation and immune-related terms, such as "response to stimulus" and "immune system process." Most genes were enriched for cellular process and cellular anatomical entity terms (Figure 2). The enrichment analysis results obtained based on KEGG terms are shown in Figure 3. These unigenes were involved in 34 functional pathways and five categories comprising Cellular Processes, Environmental Information Processing, Genetic Information Processing, Metabolism, and Organismal Systems. The number of genes enriched for Organismal Systems was largest (1825), and most genes were involved in regulating the endocrine system and immune system.

3.3. Statistics for DEGs

DEGs were screened in gill tissue samples collected from ZZ and YT in the spring, autumn, winter and summer, and differential gene volcano maps (Figure 4) were drawn. The number of DEGs between ZZ_D vs ZZ_X was 10043, with 5461 upregulated genes and 4582 downregulated genes. Among the sampling points in YT, YT_D vs YT_X had the highest number of DEGs (20909), with 9503 upregulated genes and 11406 downregulated genes. The results showed that the number of upregulated genes was higher than that of downregulated genes in most comparisons, and the number of DEGs was highest in the comparison of winter and summer.

3.4. Enrichment analysis of DEGs

The enriched GO and KEGG terms determined for DEGs screened in ZZ_D vs ZZ_X and YT_D vs YT_X are shown in Figure 5 and Figure 6, respectively. GO analysis showed that 16 terms were significantly enriched in ZZ_D vs ZZ_X, and nine terms were significantly enriched in YT_D vs YT_X. In particular, many DEGs were related to the cellular nitrogen compound metabolic process and biosynthetic process terms, and only two terms were significantly enriched for cell components in YT_D vs YT_X. The significantly enriched terms in ZZ_D vs ZZ_X included intracellular, protein-containing complex. Among the molecular functions, the significantly enriched terms included oxidoreductase activity and structural constituent of ribosome. The results obtained by KEGG term analysis showed that the enriched pathways included ribosome, spliceosome, and proteasome pathways, and they were related to protein synthesis. Some pathways were related to immunity, such as the cytosolic DNA-sensing pathway and NF-kappa B signaling pathway. It should be noted that the energy metabolism-related pathway for oxidative phosphorylation was significantly enriched in both comparisons, and terms for many significant disease-related pathways were also enriched.

3.5. Seasonal changes in DEGs

We focused on DEGs related to immunity, oxidative stress, and metabolism. We screened 13 immune-related genes, which mainly included pattern recognition receptors, cytokines, phosphatases, and small molecule proteins (Figure 7A). The clustering results showed that the expression levels of most genes were higher in summer and autumn, and lower in spring and winter, such as pattern recognition receptors comprising toll-like receptor 2 (TLR2), toll-like receptor 3 (TLR3), interferon regulatory factor IRF5, small molecule protein interleukin 17-like protein (IL-17), and CRT related to apoptosis. We suggest that the high temperature season was more like to promote the expression of these genes, and thus immunity was stronger in summer and autumn. We then screened 12 genes related to oxidative stress, which were mainly assigned to three categories comprising HSPs, mitogen-activated protein kinases (MAPKs), and GSTs (Figure 7B). The results showed that the expression levels of most MAPKs and HSPs increased with the increase in seasonal temperature, such as p38 MAPK, mitogen-activated protein kinase kinase kinase 13-like (MAP3K13), HSP 68-like (HSP68), HSP70, and HSP90. Figure 8 shows that the expression levels of these genes were low in spring, high in summer, and lowest in winter. However, the expression levels of GSTs were relatively low in spring and summer, with gradual increases in autumn, and highest in winter. The nine metabolism-related genes were mainly involved in protein synthesis and energy metabolism (Figure 7C). The expression levels of these genes were highest in winter under low temperature conditions and relatively low in summer under high temperature conditions. We consider that Neptunea cumingii can adapt to low temperatures by speeding up its metabolism.

3.6. Validation analysis using qRT-PCR

To validate the sequencing data, 11 DEGs were selected for qPCR analysis (Figure 8). The results showed that IRF5, GSTK1, CRT, TLR8, and CAAP1 were upregulated, whereas HSP90, HSP70, TLR, TLR13, MAP2K6, and Ptp2 were downregulated. The expression patterns of the selected genes detected by qPCR were consistent with the RNA-seq expression profiles, thereby corroborating our results.

4. Discussion

Temperature is the main environmental factor that affects the immune function and metabolic processes in mollusks [16]. Studies have shown that changes in water temperature can cause significant changes in immune and metabolism related genes in shellfish [7,16,28,46]. Gills are the main organs involved in respiration and immune regulation in aquatic organisms [2,12], and they can rapidly sense changes in the water temperature because they are in direct contact with the water environment. Moreover, the large surface area of shellfish gills provides a greater contact area to absorb nutrients in water, thereby facilitating energy metabolism and material exchange. In the present study, transcriptome analysis was conducted using gill tissue samples from Neptunea cumingii collected in spring, summer, autumn, and winter to elucidate the regulatory mechanisms related to adaptation to seasonal temperature changes. The results showed that more DEGs were found in the comparison of winter and summer samples, with the highest number of DEGs (20909) for YT_D vs YT_X. Further screening of DEGs identified 34 temperature-related genes, where 13 played roles in innate immunity, 12 were related to oxidative stress, and nine were related to protein synthesis and energy metabolism. Similarly, Jiang et al. [16]conducted transcriptome analysis based on gill tissues from Patinopecten yessoensis under fluctuating water temperatures and found 696 DEGs, including 41 genes related to immunity and 16 genes related to protein metabolism. Furthermore, most of the immune pathways and apoptosis pathways in Chlamys farreri were activated under heat stress, where 239 DEGs identified by transcriptome analysis were mostly related to immune signaling pathways [28]. Similar results were obtained in the present study.

In addition, we found that most genes related to immune regulation were highly expressed in summer and autumn when the temperature was higher, with lower expression levels in spring and winter when the temperature was lower, including those of pattern recognition receptor TLRs, transcription factor TNF, and small molecule protein IL-17. Pattern recognition receptors are sensory receptors in organisms that allow them to resist invading pathogens and they initiate immune responses [39]. TLRs are the main pattern recognition receptors that recognize different pathogen-associated molecular patterns, before triggering signal transduction pathways and regulating the production of antimicrobial peptides [1]. In particular, TLR2 is the main mediator of macrophage activation, where it can recognize various bacterial products and activate the myeloid differentiation factor 88 (MyD88) signaling pathway [23]. In the present study, four types of TLRs (TLR2, TLR3, TLR10, and TLR13) were significantly expressed in summer and autumn, whereas their expression levels were low in spring and winter, thereby indicating that long-term low temperature conditions reduce the activity of immune factors and increase the susceptibility to pathogens. In addition, our results indicated the inhibition of TNF-family members (TNFSF10, TNFSF11, and TNFAIP3) as downstream signaling intermediates of the TLRs signaling pathway under low temperature conditions in the winter. Similarly, Jiang et al. [16] stimulated Patinopecten yessoensis with fluctuating water temperatures and found that TNF was downregulated in the gills and the TLRs pathway was inhibited. TNF can coordinate with the Imd, Toll, and JNK pathways to control the phenoloxidase activity and expression of antimicrobial peptide genes to regulate the innate immune response [42], and the NF-κB signaling pathway is crucial for mediating the cellular TNF response [43]. The NF-κB signaling pathway promotes cell survival by inducing anti-apoptotic gene expression to regulate the body's immune and inflammatory responses [29]. ILs are small molecular proteins secreted by immune cells that play immunomodulatory roles between cells. At present, only two types of ILs [36] comprising IL-12 [55] and IL-17 [14,22,34,54] have been found in shellfish. IL-17 plays an important role in regulating the inflammatory response and innate immune response in shellfish. For example, CgIL17-1 can bind to its receptor CgIL17R1 to promote the proliferation of blood lymphocytes in Crassostrea gigas [5]. In Patinopecten yessoensis [16], IL-17 was significantly downregulated after several temperature fluctuations, and the susceptibility to bacterial infection increased under these conditions. Yang et al. [58]observed the upregulation of IL-17 in Ctenopharyngodon idella after high temperature stimulation, thereby indicating that high temperature could promote the expression of this gene. In the present study, the expression level of IL-17 was higher in summer and lower in winter. We found that the high temperatures in summer promoted the expression of IL-17 in the gills of Neptunea cumingii, but the lower temperatures in winter decreased the expression of this gene. In summary, we suggest that metabolism was slow in Neptunea cumingii when the water temperature was low in the winter, and thus the immune cell activity was low and it was more susceptible to infection by pathogens. In shellfish, the immune function is mainly mediated by blood cells. Studies have shown that low temperatures can reduce the cellular immune activity mainly by decreasing cellular metabolism and protein inactivation [56]. In addition, the lower water temperature in winter will decrease the intake of food by shellfish and affect their immune function.

Oxidative stress is a common problem that causes huge economic losses in the animal production industry. Oxidative damage is caused by the increased generation of free radicals and/or a decrease in the ability to scavenge them due to disturbance of the oxidation system and antioxidant system, thereby resulting in the accumulation of free radicals in the body [21]. Temperature is one of the main causes of oxidative stress in shellfish. Studies of Patinopecten yessoensis [16], Haliotis discus hannai Ino [15], Sinonovacula constricta [10], Crassadoma gigantea [51], Scapharca broughtonii [57], and other shellfish have confirmed that temperature can directly affect their antioxidant systems. GST is a key antioxidant enzyme involved in the clearance of reactive oxygen species (ROS), ligand transport, and regulating stress-mediated cell signaling pathways [19]. In the present study, we found that GSTs were significantly expressed at higher levels in the gills of Neptunea cumingii under low temperature conditions in the winter, and similar results were also obtained previously for Pomacea canaliculata. Lin et al. [24] found that PcGST5 and PcGST11 were significantly upregulated in the hepatopancreas of Pomacea canaliculata after cold stress. The production of ROS significantly increased in Mytilus galloprovincialis after cold stress at 4°C (P < 0.01), where the peak value was 8.4 times that in the control group [47]. Studies have shown that low temperature can induce oxidative stress in aquatic animals [27], and they may produce more ROS under low temperature stress. Subsequently, the production of GST will be activated in the body to adapt to low temperature. The MAPK signaling pathway plays an important role in mollusk immunity and the response to external environmental stresses [4], and MAPKs are important components of this pathway. In the present study, MAPKs (MAP2K6, MAP3K5, MAP3K13, and p38 MAPK) were expressed at significantly higher levels in the summer but inhibited in the winter, and p38 MAPK is involved in regulating the P38 signaling pathway. We suggest that due to the low water temperature in winter, p38 MAPK was inhibited in the gills of Neptunea cumingii to affects its cold resistance. Previous studies showed that inhibiting the p38 MAPK activity in zebrafish (Danio rerio) could significantly reduce the cold resistance by larvae [32], and the apoptosis rate and ROS content of blood lymphocytes were reduced in Sinonovacula constricta by P38 MAPK expression [52]. Moreover, Sun et al. found that P38 MAPK could regulate the expression levels of IL-17 and TNF in Crassostrea gigas [37,38]. Therefore, we consider that the MAPK signaling pathway is involved in regulating the responses to external environment stresses by shellfish, as well as playing an important role in regulating the immune and inflammatory responses induced by pathogens. In addition, HSPs act as molecular chaperones and play important roles in stress resistance and environmental adaptation [33]. Studies have shown that HSPs are associated with heat resistance by mollusks, and their expression levels are related to the environmental conditions. Thus, HSPs are often used as stress biomarkers [35]. HSPs can also influence the immune response in mollusks by participating in macrophage activation and cytokine or chemokine production [45]. HSP70 and HSP90 are the two main molecular chaperones [6,25,35]. For example, treatment of Mercenaria mercenaria by heat stress and hypoxia stress increased the expression levels of HSP90 and HSP70 [13]. Moreover, cold treatment of Ruditapes philippinarum [30], Mytilus galloprovincialis [47], and Crassostrea gigas [62] decreased the expression levels of HSP70 in the gill tissues. In the present study, we found that the expression levels of HSP70 and HSP90 were higher in the gills of Neptunea cumingii during the summer. However, these genes were inhibited and their expression levels were lower under low temperatures in the spring and winter. We suggest that the cold environment negatively affected the transcription of HSPs in the gills of Neptunea cumingii by reducing the protein activity levels, disrupting the integrity of organelles, and inhibiting important processes such as transcription and mRNA translation.

Oxidative phosphorylation is an important component of biological energy metabolism. In this process, ADP synthesizes ATP by coupling the energy released via the oxidation of substances in the body with inorganic phosphorus in mitochondria. In particular, 95% of the ATP in the body is produced by oxidative phosphorylation metabolism [47], and it has very important roles in biological growth, metabolism, and environmental adaptability. Oxidative phosphorylation occurs via the electron transfer chain, which is also known as the respiratory chain, in a continuous reaction system comprising a series of hydrogen and electron reactions, where the hydrogen atoms removed by metabolites are transferred to oxygen to generate water, as well as generating ATP [47]. The energy generated by oxidative phosphorylation is strongly related to environmental stress in marine organisms [44]. In the present study, we found that NADH dehydrogenase, V-type H+-transport ATPase subunits (ATPeV0E), and F-type H+-transport ATPase subunits (ATPeF0C, ATPeF1A, and ATPeFG) were significantly upregulated in the gills of Neptunea cumingii under low temperatures in the winter. These key enzymes are involved in ATP synthesis and transport, and thus our findings indicate that ATP synthesis and transport were enhanced under low temperatures. Previous studies have demonstrated significant increases in ATPase in the hepatopancreas and gill of Litopenaeus vannamei [49], and ATPase in the adductor muscle tissue of Pinctada fucata martensii [50] under low temperature conditions. Johnston et al. [18] measured the ATPase activities in the muscles of 19 fish species under long-term low temperature stress and found that the ATPase activities were higher compared with those in fish in warm water. In addition, low temperature can lead to reductions or the loss of enzyme activities, and thus we consider that increasing the transcription of enzymes during oxidative phosphorylation is effective for improving the production of ATPase and slowing decreases in enzyme activities under low temperature conditions. The KEGG enrichment analysis results obtained in this study showed that many DEGs were significantly enriched in the ribosome pathway, and most were ribosomal proteins (RPs). A previous study also showed that the expression of RP-L19 decreased in the gills of Patinopecten yessoensis after repeated heat stimulation [16], and similar results were obtained in our study. In particular, the expression levels of several RP genes (RP-S5, RP-L19, RP-L32, RP-L28e, RPL28, and RP-L19e) decreased under high temperatures in the summer but increased under low temperatures in the winter. Interestingly, the expression patterns of these RP genes differed during the autumn at the two study locations, where the expression levels of these RP genes were upregulated during the autumn in ZZ but downregulated during the autumn in YT, and this difference may have been related to the latitudes of the two locations, and thus further study is required. We suggest that Neptunea cumingii adapted to the low temperature environment during the winter by synthesizing RPs, which had stable functions in the subregions of ribosomes [3]. These RPs promoted or inhibited ribosomes in different seasons by replacing damaged subregions in order to ensure the structural and functional stability of ribosomes.

5. Conclusions

In the present study, the transcriptome results were significantly different for the gill tissues of Neptunea cumingii in different seasons. The expression levels of immune regulatory genes were promoted in Neptunea cumingii during the high temperature season, whereas the expression levels of genes related to energy metabolism and RP genes were inhibited. Genes related to immune regulation were inhibited in Neptunea cumingii during the low temperature season and its immunity was relatively low. Neptunea cumingii adapted to the low temperature environment by activating GSTs in the body and accelerating energy metabolism. These results provided a basis for understanding how the molecular defense mechanism of Neptunea cumingii might adapt to climate change, as well as a reference for improving the production of Neptunea cumingii.