The sexual maturity S. esculenta were collected in Qingdao sea area, China, and temporarily raised until laying eggs. The eggs were collected and placed in a breeding pool with flowing seawater and continuously oxygenated during hatching. The eggs hatched after a month and were divided into two groups of 50 individuals in each group. Larvae of control group (C) grew in normal seawater, and exposed group (NPs) larvae grew in seawater with NPs (100 mg/L). Then, the above larvae were collected at 0 h (C_0 h) and 4 h (C_4 h and NPs_4 h) respectively. The collected larvae were stored in liquid nitrogen for future use.

The activity of SOD were measured by Total Superoxide Dismutase (T-SOD) assay kit (Hydroxylamine method) purchased from Nanjing Jiancheng Bioengineering Institute. In addition, the MDA levels were detected using Malondialdehyde (MDA) assay kit (TBA method) according to the instructions. Eight replicates were set in each group (C_4h and NPs_4h), and the measured tissue fluid was obtained by grinding 8 randomly selected whole larvae.

The calculation formula for SOD activity is as follows: A = O D   v a l u e   c o n t r a s t   t u b e − t e s t i n g   t u b e O D   v a l u e   o f   c o n t r a s t   t u b e ÷ 50 % × V 1 V 2 ÷ C Notes: A, total SOD activity (U/mg prot); V1, total volume of reaction solution (mL); V2, volume of sampling amount (mL); C, protein concentration of testing sample (mg prot/mL).

The calculation formula for MDA level is as follows: L = O D   v a l u e   t e s t i n g   t u b e − c o n t r a s t   t u b e O D   v a l u e   s t a n d a r d   t u b e − b l a n k   t u b e × C 1 ÷ C 2 Notes: L, MDA level (nmol/mg prot); C1, concentration of standard substance (10 nmol/mL); C2, protein concentration of testing sample (mg prot/mL).

We used the TRI Reagent method (Rio et al., 2010) with the manufacturer’s protocol to extract total RNA and identified the integrity using Agilent 2100 bioanalyzer (Sodowich et al., 2007). Nine larvae were randomly selected for RNA extraction from groups C_0h, C_4h, and NPs_4h, respectively. Then, the RNA of 9 larvae with equal molar masses in each group was mixed into 3 replicates for subsequent sequencing. Using NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® to construct the transcriptome library (Parkhomchuk et al., 2009). Raw reads were sequenced by Illumina NovaSeq 6000 (Illumina, United States), whose SRA accession number were SRR23936172, SRR23936173, SRR23936174, SRR23936175, SRR23936181, SRR23936182, SRR25114243, SRR25114244 and SRR25114245. Removing low quality reads from raw reads to obtain clean reads. The obtained clean reads were mapped to the reference genome (unpublished) using HISAT2.

In this study, the DESeq2 software of R was used as a model to screen differentially expressed genes (DEGs). First, the data were involved in constructing the ddsmodel, after which the dispersion of the samples was estimated using the DESeq function, and finally the differences in gene expression were analyzed. DEGs with p-value ≤0.05 to compare groups C_4h and NPs_4h were screened out (Love et al., 2014).

The functional enrichment analysis was performed on DEGs. To ascertain the GO terms and the distribution of DEGs, GO analyses were deployed on the union set distinguished at two distinct time points. Additionally, Gene Set Enrichment Analysis was employed to identify immune-related pathways and genes through the KEGG pathway analysis, thus elucidating the functions of DEGs. Enrichment analyses of GO and KEGG were executed using the DAVID database (https://david.ncifcrf.gov/) 2021 (Jiao et al., 2012). The construction of a protein-protein interaction (PPI) network can offer insights into the correlations amongst oxidative stress pathways, thereby simplifying the identification of pivotal genes. In this study, we leveraged the STRING database (https://cn.string-db.org) to construct a robust PPI network (Szklarczyk et al., 2019).

The accuracy of RNA-Seq was verified via qRT-PCR. In this study, 30 hub genes were pinpointed for validation via qRT-PCR. Utilizing Primer Premier 5.0 software, gene primer sequences were formulated based on the spliced transcriptome. The primers related information used in the qRT-PCR is shown in Supplementary Table S1. The gene of β-actin is used as housekeeper gene for qRT-PCR due to its evident stability within this experiment. The fluorescence quantification methods implemented were adapted from Liu et al.'s work (Liu et al., 2017).

As shown in Figure 1A, the activity of SOD was higher in the S. esculenta larvae exposed to PS-NPs than that in control group (p-value = 3.1E-5). The results of RNA-Seq data analysis showed that the transcript levels of SOD were also increased in S. esculenta larvae exposed to PS-NPs (Supplementary Table S2). Besides, the level of MDA in S. esculenta larvae whole body was increased after PS-NPs exposure (p-value = 6.0E-6) (Figure 1B).

The changes of physiological biomarkers pointed out that the acute exposure of PS-NPs induced oxidative stress in S. esculenta larvae. To explore the molecular mechanisms involved, transcriptome sequencing projects were performed. The average of 45,989,845 raw reads per sample were sequenced. Subsequently, after filtering, an average of 44,771,357 clean reads were generated for each sample. The average of Q20 and Q30 were 96.91% and 92.00%, respectively. And the average of GC content in clean reads was 40.24% (Supplementary Table S3). These results suggested a high quality of sequencing.

As the results of differential expression analysis, there were a total of 2570 DEGs (1,166 up- and 1,404 downregulated) at 4 h after PS-NPs exposure (Figure 2A). DEGs expression distribution of all groups was shown in the heatmap, which displayed an obvious difference in the expression pattern of DEGs between the PS-NPs exposure group and non-exposed groups (Figure 2B).

To investigate the function of DEGs, GO and KEGG enrichment analysis were conducted. As shown in Figure 2C, 145 significant GO terms were enriched, and translation, structural constituent of ribosome, and other terms are important for mediating oxidative stress. Based on level-2 KEGG enrichment analysis results, these DEGs were participating in multiple signaling pathways (Figure 2D). The enrichment of 14 KEGG signaling pathways (level-3), such as Ribosome, Proteasome, and MAPK signaling pathway, suggested that PS-NPs exposure have affected multiple biological processes in S. esculenta larvae (Table 1).

Subsequently, 139 DEGs involved in KEGG signaling pathways of Table 1 were used to construct the PPI network to scan hub genes regulating multiple biological processes affected (Figure 2E). The relevant parameters of PPI network are displayed in Supplementary Table S4. Considering the numbers of involved in KEGG signaling pathway and protein interaction, 30 key DEGs were obtained, for instance, RPS5, RPS9 and MRPL4 (Table 2).

The key DEGs were chosen for quantitative RT-PCR. The results showed the trend of qRT-PCR was consistent with that of RNA-Seq, which suggested that results of transcriptome profile were reliable and accurate (Figure 3).

Oxidative stress and even oxidative damage are common in various organisms exposed to NPs, such as mammalian cells, Danio rerio, D. pulex, and M. galloprovincialis (Sarasamma et al., 2020; Liu et al., 2021a; Banerjee and Shelver, 2021; Gonçalves et al., 2022). The occurrence of oxidative stress is triggered by high doses of reactive oxygen species (ROS, like superoxide anion). The removal of ROS in aerobic organisms relies on antioxidant systems. SOD, an important member of the antioxidant enzyme system, has the ability of catalyzing the disproportionation of superoxide anion to hydrogen peroxide, and is used as a biomarker for oxidative stress (Fukai and Ushio-Fukai, 2011). Studies have shown that NPs can enter cells through internalization pathways, causing raise of ROS production and SOD activity, and even directly binding to SOD to alter its activity (Schröter and Ventura, 2022; Wang et al., 2022). In addition, the increase of ROS production could lead to lipid peroxidation. The metabolomics analysis results confirmed that the metabolism of membrane lipids in Sinonovacula constricta was affected by acute exposure to PS-NPs (Jiang and Zhang, 2021). MDA, a biomarker of oxidative stress, is an important membrane lipid peroxidation product and can also intensify membrane damage (Sillero-Ríos et al., 2018). After exposure of PS-NPs, SOD activity significantly enhanced, indicating an increase in ROS production. This result is consistent with the improvement of SOD activity in the Ietalurus punetaus larvae and in the gills and digestive glands of Mytilus spp. caused by acute exposure to NPs (Cole et al., 2020; Jiang et al., 2022). Moreover, the improvement of MDA level suggested the occurrence of membrane lipid oxidation. Research has found that SOD activity and MDA levels also increased in intestine of A. fangsiao after exposure to high concentrations of MPs (Zheng et al., 2022). These results indicated that acute exposure to high concentration of PS-NPs can disrupt the redox balance and cause oxidative stress response, which reminded us the hazards of acute exposure to high concentrations of NPs.

Transcriptome analysis is widely used in the study of molecular level in biology along with the development of sequencing technologies. In recent years, transcriptome analysis has gradually been applied in the toxicology research of NPs on aquatic organisms, such as B. plicatilis, D. pulex, Isognomon alatus, Mytilus coruscus, Cherax quadricarinatus, C. robusta, D. rerio, Oreochromis mossambicus, I. punetaus (Liu et al., 2021a; Pang et al., 2021; Arini et al., 2022; Yu et al., 2022b; Cheng et al., 2022; Jiang et al., 2022; Shin and Jeong, 2022; Eliso et al., 2023; Qi et al., 2023). To research the effect of NPs on the molecular level of S. esculenta larvae, we applied RNA-Seq technology in this research.