All E. leucolegnote collected in Qidiao Village (110°38′26″ E, 19°56′31″ N, Haikou, Hainan, China) were cultured in laboratory (117°16′72″ E, 31°91′38″ N, Hefei, Anhui, China) for 7 days before experiments, to allow the sea slugs to acclimate to the experimental conditions. E. leucolegnote were placed on 1.1 T neodymium iron boron (NdFeB) N38 permanent magnets (length× width× height: 5.8 × 4.8 × 3.8 cm), namely 1.1 T static magnetic field (1.1 T SMF). To measure the distributions of the magnetic fields at different positions, a magnet analyzer (FE-2100RD, Forever Elegance, China) was used to scan the SMF distribution above the magnets. Pilot tests were carried out to determine the suitable SMF treatment time. The preliminary data suggested that the mortality rate of sea slugs in 1.1 T SMF group reached 70% at 14 days, and decreased to 40% at 13 days. And all sea slugs survived after 10 days with 1.1 T magnetic field treatment. Thus, 10 days of SMF treatment was chosen to explore the cause of death and the biological effects of magnetic field treatment for this study. To minimize the experimental variations, the control group (here named geomagnetic field treatment group, GMF group, with inclination of 47.6765°, declination of -5.5585° and total field intensity of 49,892.9 nT, which corresponding to 0.49 Gauss) was placed on aluminium block of the same size as 1.1 T SMF, and keep distant from the 1.1 T SMF group to eliminate the magnetic disturbance from artificial magnetic field. A total of 20 E. leucolegnote with an average length 4 mm were randomly assigned to two groups and placed in a plastic Petri dish (length× width× height: 4.3 × 1.1 × 1.1 cm) filled with self-configured seawater with a salinity of 25, pH 8.1 and sealed with a sealer. E. leucolegnote were placed in each group in an air-conditioned room maintained at 25°C, water temperature 18°C, and follow with the changes of 12 h light: 12 h dark for consecutive 10 days. After experiments, E. leucolegnote from each group were quickly frozen at −80°C ultra-low temperature refrigerator individually until subsequent enzyme activity and transcriptomics analysis.

Glucose (Glu), Triglyceride (TG), Total cholesterol (TCH), enzyme activity of Catalase (CAT), Superoxide dismutase (SOD), Glutathione (GSH), Glutathione peroxidase (GSH-PX), Lysozyme (LZM), Aspartate aminotransferase (AST), Alanine transaminase (ALT), Amylase (AMS), Pepsin (PEP), Lipase (LPS) and Trypsin (TRY) were measured using kits purchased from Nanjing Jiancheng Bioengineering Institute (Nanjing, China).

Total RNA was extracted from the E. leucolegnote using TRIzol® Reagent according the manufacturer’s instructions (Invitrogen) and genomic DNA was removed using DNase I (TaKara). The integrity and purity of the total RNA quality was determined by 2100 Bioanalyser (Agilent Technologies) and quantified using the ND-2000 (NanoDrop Technologies). Only high-quality RNA sample (OD 260/280 = 1.8–2.2, OD 260/230 ≥ 2.0, RIN≥ 6.5, 28 S: 18 S ≥ 1.0, >1 μg) was used to construct sequencing library.

Isolation of mRNA by oligo (dT) beads following the polyA selection method, and sequencing libraries were constructed by the TruSeqTM RNA sample of preparation Kit (Illumina, San Diego, CA). The first strand cDNA was synthesized using reverse transcriptase, random primers, and the short fragments, followed by second strand cDNA synthesis. Secondly double-stranded cDNA was synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen, CA) with random hexamer primers (Illumina). Then, the synthesized cDNA was subjected to end-repair, phosphorylation and “A” base addition. Libraries were size selected for cDNA target fragments of 300 bp on 2% Low Range Ultra Agarose followed by PCR amplified using Phusion DNA polymerase (NEB) for 15 PCR cycles. After quantified by TBS 380, paired-end RNA-seq sequencing library was sequenced with the Illumina NovaSeq 6000 sequencer (2 × 150 bp read length).

The clean reads data were obtained from the raw reads trimmed. All clean data from the samples were used to do de-novo assembly with Trinity (Grabherr et al., 2011). Then the assembled transcripts were assessed and optimized with BUSCO (Benchmarking Universal Single-Copy Orthologs) (Manni et al., 2021), TransRate (Smith-Unna et al., 2016) and CD-HIT (Fu et al., 2012). All the assembled transcripts were searched against the NCBI protein non-redundant, Swiss-Prot12, Pfam, Clusters of Orthologous Groups of proteins, GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) databases using BLASTX to identify the proteins that had the highest sequence similarity with the given transcripts to retrieve their function annotations and a typical cut-off E-values less than 1.0 × 10⁻⁵ was established.

To identify DEGs (differential expression genes) between two different treatments, the expression level of each gene was calculated according to the transcripts per million reads (TPM) method. RSEM (Li and Dewey, 2011) was used to quantify gene abundances. Essentially, differential expression analysis was performed using the DESeq2 (Love et al., 2014) with |log₂ (foldchange)| ≥ 1 and P-adjust ≤.05. Functional-enrichment analysis including GO and KEGG were performed to identify which DEGs were significantly enriched in GO terms and metabolic pathways at P-adjust ≤.05 compared with the whole-transcriptome background. GO functional enrichment and KEGG pathway analysis were performed using Goatools and KOBAS (Xie et al., 2011), with adjusted p < .05 using the Benjamini–Hochberg method.

Eight DEGs from transcriptomic studies were randomly selected, and one gene (SSR3) with no expressional differences was chosen to validate the RNA-seq results by RT-qPCR. A housekeeping gene β-actin was used as an internal control to normalize the expression of the target genes. Primers for RT-qPCR were designed using Primer Premier 5.0 software, and primers’ details were listed in Supplementary Table S1. Total mRNA of the samples was extracted using above-mentioned method. Then, cDNA was synthesized using the PrimeScript^® RT Reagent Kit with gDNA Eraser (Takara, China), according to the manufacturer’s protocol. The SYBR Green RT-PCR assay was conducted to determine mRNA expressions of the genes. The temperature programming conditions were: denaturation at 95°C for 30 s, followed by 35 cycles of 95°C for 10 s, 60°C for 30 s, 95°C for 15 s, 60°C for 60 s, and 95°C for 1 s the melting curve was analyzed to confirm the presence of a single product in these reactions. Gene expression results were obtained using the 2^−ΔΔCT method (Livak and Schmittgen, 2001).

There are at least three biological replicates, excluding RNA-seq, for each sample. GraphPad Prism five was used for histogram and statistical analysis. Student’s t-test was used to examine the raw data. Differences were considered significant at *p < .05.

E. leucolegnote in SMF group were exposed to static magnetic field generated by permanent magnets as shown in Figure 1A. The surface of SMF distribution is uneven and ranges from 0.1 T to 1.1 T. The mollusc were placed in the middle area (black rectangle in Figure 1B), where the intensity of magnetic field could reach the highest 1.1 T. In this study, we firstly measured the glucose and lipid metabolism of E. leucolegnote and compared them in both GMF and SMF groups. Significant increase of blood glucose and lipid levels were observed after 10 days of SMF exposure (Figure 2, p < .05), represented by the increased glucose (Glu, Figure 2A), triglyceride (TG, Figure 2B) and total cholesterol (TCH, Figure 2C).

Static magnetic field treatment has significant effects on oxidative stress in E. leucolegnote. Obvious changes of antioxidant enzyme activities upon exposure to 1.1 T SMF were observed (Figure 3). The activity of several key antioxidant enzymes including CAT (Figure 3A), SOD (Figure 3B), GSH (Figure 3C) and GSH-PX (Figure 3D) in the 1.1 T SMF group were significantly higher than those in the GMF group (p < .05), while the LZM activity of 1.1 T SMF group was lower than those in GMF group (p < .05) (Figure 3E). And the activities of AST (Figure 3F) and ALT (Figure 3G) in the 1.1 T SMF group were significantly higher than those in GMF group (p < .05).

Magnetic field has a significant impact on digestive enzyme activity as well (Figure 4). Several enzymes related to digestive performance of digestive glands including AMS (Figure 4A), PEP (Figure 4B) and LPS (Figure 4C) were down-regulated in the 1.1 T SMF group, compared with those in GMF group (p < .05), but no significant differences in TRY activity were detected among two groups (p > .05, Figure 4D). Since AMS, PEP and LPS are associated with the primary macronutrients in our diet, such as carbohydrates, proteins, and fats respectively, the activity of these enzymes represents a foundational aspect of gastrointestinal health. The decreased activity of these enzymes suggested 1.1 T SMF exposure led to decreased digestive performance in E. leucolegnote. The same finding was seen for the downregulated genes in the transcriptome as shown below.

To further investigate the possible molecular mechanism of the magnetic effects on E. leucolegnote, a comparative analysis of the gene expression profiles was performed. We sequenced and compared the transcriptome of E. leucolegnote exposed to either the 1.1 T static magnetic field or the geomagnetic field. More than 300 million raw were obtained, and an average of 49.7 and 47.6 million clean reads were localized to the E. leucolegnote genome from the different magnetic field treatment, respectively (Table 1). The GC contents of two groups showed the values with 43.96% and 43.36%. And Q20 (those with a base quality greater than 20) contents were in the average of 97.78% and 97.68%, Q30 (those with a base quality greater than 30) contents were in the average of 93.63% and 93.40% (Table 1). These data suggested that the quality of sequencing results was sufficiently high and reliable for the subsequent transcriptome analysis.

In total, 836 differentially expressed genes (DEGs) were identified in 1.1 T SMF treatment compared with GMF group, among which 352 genes are upregulated, and 484 genes are downregulated (Figure 5A and Supplementary Table S2). To identify the processes enriched in significant DEGs, we subjected significant DEGs to gene ontology (GO) functional annotation and term enrichment analysis, a tool developed to represent common and basic biological information. Most DEGs were involved in the biological processes such as cellular process and metabolic processes, catalytic activity and binding in the molecular function, and cell part and membrane part in the cellular component. In our study, GO enrichment analysis corresponding to 836 significant DEGs were produced and assigned to 69 functional groups and three categories (Figure 5B). Most of DEGs in the comparison of 1.1 T SMF vs. GMF were enriched in cellular component such as “extracellular region,” “integral component of membrane,” “intrinsic component of membrane” and “cellular anatomical entity,” and three additional in molecular function such as “chitin binding,” “sulfuric ester hydrolase activity” and “carbohydrate binding” (Figures 5B–D). Mapping these DEGs to the pathways from databases KEGG suggested that these genes are significantly clustered into several key signaling pathways, namely phagosome lysosome, apoptosis and endocytosis in cellular processes, tuberculosis in human diseases, glycosaminoglycan degradation and other glycan degradation in metabolism, vitamin digestion and absorption in organismal systems and NF-kappa B signaling pathway in environmental information (Figures 5E, F).

Significantly upregulated genes and downregulated genes were selected based on the annotation of DEGs. The upregulated genes were significantly clustered into biological processes including “vacuolar transport,” “sulfuric ester hydrolase activity,” “cysteine-type peptidase activity” and “hydrolase activity” in molecular function, and “lysosomal membrane” and “lytic vacuole membrane” in cellular component (Figures 6A–C). Mapping these upregulated genes to the pathways from databases KEGG suggested that these upregulated genes are significantly enriched in the signaling pathways of lysosome and apoptosis in cellular processes (Figures 6D, E). It is worth pointing out that both the caspase and lysosome change of E. leucolegnote after 10 day of SMF treated are related and involved in apoptosis pathway based on the apoptosis map, this further validates the involvement of lysosomal proteases (Turk et al., 2002) and caspase family (Bearoff and Fuller-Espie, 2011) in the apoptotic pathways.

Taking together, our data suggested that the upregulated genes upon SMF treatment identified in this study were mostly leading to apoptosis and increased lysosomal activity, which is consistent with our biochemical studies as shown in Figures 3, 4, and in an agreement with previous reports as well. For example, SMF exposure alone or in combination with extremely low frequency MF (ELF) of more than 1 mT have a selective impact on cell signaling related to apoptosis, which may through magnetic field effect on motion of charged molecules (Tofani et al., 2001). In another study, after 12 days of exposure to the 1 mT electromagnetic field generated by the submarine cable, the number of apoptotic in the Limecola balthica was significantly elevated (Stankeviciute et al., 2019). It has been suggested that the underlying mechanism of how magnetic field affect apoptosis may through magnetic field effect on electron spin, which lead to increased reactive oxygen species (ROS) concentration (Tofani, 2022).

We also analyzed the downregulated genes by GO annotations and KEGG enrichment analysis (Figure 7). GO functional enrichment showed that most downregulated genes were enriched in molecular function such as “chitin binding” and “transmembrane transporter activity,” and in cellular component such as “extracellular region,” “cellular anatomical entity,” and “microtubule-based process” and “microtubule-based movement” as well in biological process (Figures 7A–C). Mapping these DEGs to the pathways from databases KEGG suggested that these downregulated genes are significantly clustered into the signaling pathways of phagosome, digestive systems and infectious diseases (Figures 7D, E).

Taking together, our data suggested that among the total of 484 downregulated genes, many of them might lead to the reduced digestive performance and immune systems, which could increase the risk of infectious diseases. The analysis of downregulated DEGs is largely consistent with our metabolism studies as shown in Figures 3, 4.

To further validate the accuracy and reliability of transcriptomic data, eight DEGs from transcriptomic studies including significantly up/down regulated genes, were chosen, and one gene with no expressional differences were used as control. Real-time quantitative PCR (RT-qPCR) was carried out on these nine genes to measure the expression levels (Figure 8). Notably, five DEGs involved in apoptosis and lysosome pathways such as cathepsin Z (CTSZ), cathepsin-L (CTSL), vacuolar-processing enzyme (LGMN), arylsulfatase B (ARSB) and N-sulphoglucosamine sulphohydrolase (SGSH) showed significant up-regulation. In contrast, two DEGs involved in phagosome pathway such as actin G1 (ACTB_G1) and macrophage mannose receptor 1 (MRC) and one DEG involved in Glycerolipid metabolism pathway, pancreatic lipase-related protein (PLRP1), appeared to be down-regulated. And one DEG involved in Protein processing in endoplasmic reticulum pathway, translocon-associated protein subunit gamma (SSR3), only showed mild change in RT-qPCR validation. The data suggested the several signaling pathways including apoptosis, lysosome and phagosome pathway were clearly involved in the SMF response, which shed light on the poorly understand molecular mechanism of magnetic field effects on biological organisms.