Metabolic rates
MNP exposure led to a significantly higher ingestion rate in oysters exposed to 1 µg L–1 (12.4 ± 0.6 × 107 cells h–1 g–1) than in oysters exposed to the control (10.7 ± 0.5 × 107 cells h–1 g–1) and 0.025 µg L–1 (10.5 ± 0.6 × 107 cells h–1 g–1) conditions (Tukey’s HSD, P < 0.001) (Fig. 1b). Furthermore, a significantly lower assimilation efficiency was observed in oysters exposed to 1 µg L–1 (26.1 ± 5.2%) than in control oysters (37.5 ± 3.4%; Tukey’s HSD, P < 0.001), while only a downward trend was observed in oysters in the 0.025 µg L–1 condition (31.4 ± 3.3%; Tukey’s HSD, P = 0.100) (Fig. 1d). No significant difference was observed regarding oxygen consumption (ANOVA, P = 0.420), with a mean of 0.79 ± 0.08, 0.81 ± 0.07 and 0.85 ± 0.06 mg O2 h–1 g–1 under control, 0.025 µg L–1 and 1 µg L–1 conditions, respectively (Fig. 1c).
Scope for growth
Metabolic rate analysis revealed a mean energy balance (scope for growth) of 31.2 ± 3.9, 22.2 ± 4.7 and 22.6 ± 6.8 J h–1 g–1 for the control, 0.025 µg L–1 MNP and 1 µg L–1 MNP conditions, respectively (Fig. 1e). The overall Kruskal‒Wallis (KW) test was at the limit of significance (P = 0.050), and multiple comparisons revealed a slight significant difference in means between the control and 0.025 µg L–1 conditions (Dunn’s test, P = 0.048) (Fig. 1e).
Energy reserve
The glycogen content analysis on pearl oyster muscle revealed a significantly lower energy reserve in the 0.025 µg L–1 (4.6 ± 0.8 µg mg–1) condition compared with the control (8.1 ± 1.2 µg mg–1) and 1 µg L–1 (7.9 ± 1.5 µg mg–1) conditions (Tukey’s HSD test, P < 0.001 and P = 0.001, respectively) (Fig. 1f).
Gametogenesis
All pearl oysters were of the male sex, except one individual of the female sex in the control condition. No sign of abnormal gametogenesis was observed by gonad cross-section histology analysis, regardless of condition.
Shell growth
No significant difference in shell growth was observed among conditions (ANOVA, P = 0.321), as suggested by the similar shell deposit rate between conditions, with values of 5.1 ± 0.5 µm day–1 in the control, 5.1 ± 0.6 µm day–1 in the 0.025 µg L–1 condition and 5.6 ± 0.5 µm day–1 in the 1 µg L–1 condition) (Fig. 1g).
Pearl rotation
The mean angular speed of pearl rotation in the pearl sac of P. margaritifera was similar among conditions (ANOVA, P = 0.355), with values of 5.9 ± 2.6, 4.0 ± 1.5 and 6.5 ± 1.3 deg min–1 in the control, 0.025 and 1 µg L–1 conditions, respectively.
Harvest
The number of pearls collected at the end of the experiment was statistically similar among conditions in terms of frequency distribution (χ2 = 0.164, P = 0.921), with a mean of 62.4 ± 12.8% (N = 23) for individuals that produced pearls in the control condition, 59.7 ± 17.4% (N = 27) for individuals that produced pearls in the 0.025 µg L–1 condition and 59.9 ± 11.4% (N = 31) for individuals that produced pearls in the 1 µg L–1 condition (Fig. 2a and b). However, a significant difference in the frequency distribution of the production of keshi pearls, i.e. pearl oysters that have rejected their nucleus (Fig. 2a and c), was observed between the control and 1 µg L–1 conditions (Fisher’s exact test, P = 0.039); the 0.025 µg L–1 condition had an intermediate value that was not significantly different (Fisher’s exact test, P = 0.144). Indeed, a mean of 9.9 ± 8.3% (N = 4) and 17.6 ± 9.1% (N = 9) of pearl oysters produced keshi pearls in the 0.025 and 1 µg L–1 conditions, respectively, while no keshi pearls were collected in the control condition (N = 0; Fig. 2a).
Pearl quality traits
Pearl nacre deposition measured on the nucleus and composed of periostracum, calcite and aragonite crystals (Fig. 2d) reached similar mean values of 80.2 ± 12.2 µg in the control condition, 85.4 ± 17.8 µg in the 0.025 µg L–1 condition, and 78.0 ± 21.1 µg in the 1 µg L–1 condition (KW test, P = 0.717). Similarly, no significant difference in biomineralization rate (thickness) was measured among conditions (ANOVA, P = 0.430), with values of 3.0 ± 0.5, 3.4 ± 0.7 and 2.9 ± 0.6 µm d–1 in the control, 0.025 and 1 µg L–1 conditions, respectively (Fig. 2e). However, a significantly higher proportion of organic material in the coating on the nucleus (periostracum) was measured in the 0.025 µg L–1 condition (7.3 ± 6.5%) compared to the control (1.2 ± 1.9%) (Dunn’s test, P = 0.025), while no difference was observed regarding calcite and aragonite crystals (KW test, P = 0.258 and P = 0.148, respectively) (Fig. 2f). Focusing on aragonite crystals at the microstructure level, no significant difference in aragonite front spacing (Fig. 2g and h) was measured on the pearl surface (KW test, P = 0.950), with a mean of 16.3 ± 1.9 µm in the control condition, 15.9 ± 1.7 µm in the 0.025 µg L–1 condition and 18.2 ± 3.9 µm in the 1 µg L–1 condition (Fig. 2g). However, a significant difference in aragonite platelet thickness (Fig. 2i and j) was measured (KW test, P = 0.010), with values of 491.9 ± 45.3, 442.4 ± 36.1 and 395.1 ± 18.5 nm in the control, 0.025 and 1 µg L–1 conditions, respectively, revealing a significantly lower value in the 1 µg L–1 condition compared to the control condition (Dunn’s test, P = 0.002; Fig. 2i).
Keshi pearl origin
First, keshi pearl observation with a stereomicroscope led to the visualization of a purple particle (~ 9 µm), the typical color of the synthetic rope, embedded in the mineral surface microlayer of a keshi pearl produced by a pearl oyster exposed to 1 µg L–1 MNPs (Fig. 2k and l). After a soft sanding step, the particle was tentatively identified using µ-Raman but remained embedded in a thin mineral layer and was lost during additional soft sanding. This approach thus failed to identify a PE signature corresponding to the MNP produced from weathered purple synthetic rope (Fig. 2m and n). To overcome losses due to sanding, the particles were first tentatively analyzed as a whole by Py-GC‒MS, but the MS signals were not workable. In a second attempt, a keshi pearl was dissolved with 98% sulfuric acid after confirming the innocuousness of this chemical for both PE and PP fragments. With this approach, some black fragments were observed and tentatively identified using µFT-IR, but none of them matched the signature of spat collectors (PP).
RNA-Seq data
RNA sequencing of mantle (N = 29), hemocyte (N = 29) and pearl sac (N = 29) samples yielded means of 36.6 ± 4.5, 33.2 ± 3.3 and 34.4 ± 2.3 M raw reads per individual, respectively. After trimming, 95.2% of reads were recovered in mantle, hemocyte and pearl sac samples and used for downstream analyses. The mapping rate reached 46.0 ± 1.1%, 36.1 ± 0.8% and 42.2 ± 1.0% in mantle, hemocyte and pearl sac samples, respectively, with no significant differences among conditions (ANOVA, P = 0.609, P = 0.149 and P = 0.191, respectively). Sequencing results and read survival after trimming and mapping are shown in Supplementary Table 1.
Global gene expression patterns across tissues and conditions
Patterns of up- and downregulation within euKaryotic Orthologous Groups (KOG) classes were used to compare MNP conditions with global gene expression data generated from mantle, hemocytes and pearl sac of P. margaritifera in response to 0 (control), 0.025 and 1 µg MNPs L–1. A hierarchical clustering analysis comparing the change in magnitude and direction of gene expression within each KOG class among datasets revealed that the patterns of the 0.025 and 1 µg L–1 conditions were mostly similar (Fig. 3a-c) but that largely differences occurred across tissues. Indeed, P. margaritifera KOG functional enrichment correlated across conditions and tissues, with Pearson’s r values of 0.74 (P < 0.001), 0.90 (P < 0.001) and 0.56 (P = 0.004) for the mantle, hemocytes and pearl sac, respectively (Fig. 3a-c). No significant correlations in gene expression were detected among tissue samples from each MNP condition, except between mantle and pearl sac datasets for the 1 µg L–1 condition (r = 0.59, P = 0.002) (Supplementary Fig. 1).
In the mantle, individuals exposed to MNPs exhibited upregulation of genes involved in “energy production and conversion” (0.025 and 1 µg L–1, Padj < 0.001), which was the most significant enrichment, followed by “replication, recombination and repair” (0.025 µg L–1 MNPs, Padj = 0.049; 1 µg L–1 MNPs, Padj = 0.016). Regarding downregulated genes, the most significant enrichment common to both MNP conditions was associated with “cytoskeleton” (0.025 and 1 µg L–1, Padj < 0.001), although this KOG class was upregulated in hemocyte samples, followed by “signal transduction mechanisms” (0.025 and 1 µg L–1 MNPs, Padj < 0.001, respectively) (Fig. 3d). Hemocytes of exposed oysters also exhibited downregulation of genes associated with “nuclear structure” (0.025 µg L–1 MNPs, Padj = 0.047; 1 µg L–1 MNPs, Padj < 0.001), “RNA processing and modification” and “translation, ribosomal structure and biogenesis” (0.025 and 1 µg L–1 MNPs, Padj < 0.001, respectively) (Fig. 3d). Finally, in the pearl sac, individuals exposed to MNP conditions exhibited common enrichment of upregulated genes involved in “chromatin structure and dynamics” (0.025 µg L–1 MNPs, Padj = 0.009; 1 µg L–1 MNPs, Padj = 0.031) (Fig. 3d). Some unique patterns specific to MNP conditions were also observed, as shown in Fig. 3d.
Gene coexpression modules associated with MNPs and physiological traits
Mantle module functional enrichment
The expression values of 23,610 genes in 29 mantle samples were used to construct the coexpression module by weighted gene coexpression network analysis (WGCNA). Overall, 26 module genes were detected and labeled with different colors (Supplementary Fig. 2). After clustering of module eigengenes (MEs) based on dissimilarity, a total of 10 modules were selected according to module-trait relationships (r ≥ 0.45, P ≤ 0.01) regarding correlations of ME expression with experimental conditions and physiological traits (Fig. 4a and Supplementary Fig. 3). Among the modules of interest, 5 modules (cluster 1) were characterized by reduced expression under MNP conditions compared to control conditions, with the darkgrey (N = 103 genes) and turquoise (N = 5,346 genes) modules showing significantly lower ME expression at 0.025 (Dunn’s test, P = 0.043) and 1 (Dunn’s test, P = 0.023) µg L–1 MNPs, respectively (Fig. 4b). In the darkgrey module, gene ontology (GO) enrichment analysis of biological process (BP) highlighted “negative regulation of calcium ion-dependent exocytosis” (GO:0045955, Padj < 0.001). The turquoise module exhibited several representative enrichments (N = 15 GO) associated with cell organization and biogenesis, metabolism, development and signal transduction (Fig. 4c). Four modules were also identified in cluster 2, with the brown module (N = 1,606 genes) showing significantly higher ME expression at 1 µg L–1 MNPs than under the control conditions (Tukey’s HSD test, P = 0.006). The most significant enrichments identified in this module were “chitin metabolic process” (GO:0006030, Padj < 0.001), “positive regulation of lipoprotein transport” (GO:0140077, Padj < 0.001), “ribonucleoside catabolic process” (GO:0042454, Padj = 0.006) and “biological adhesion” (GO:0022610, Padj = 0.005). Finally, the tan module (N = 718 genes; cluster 3) was identified as a key module showing a strong correlation with condition (r = 0.67, P < 0.001), specifically with 1 µg L–1 MNPs (r = 0.69, P < 0.001), and strong correlations with physiological traits such as ingestion (r = 0.54, P = 0.002) and assimilation (r = − 0.45, P = 0.01) (Fig. 4a). This module exhibited significantly higher ME expression in the 1 µg L–1 group than in both the control and 0.025 µg L–1 groups (Dunn’s test, P = 0.007 and P = 0.014, respectively) (Fig. 4b). The tan module was significantly enriched for “MAPK cascade” (GO:0000165, Padj < 0.001), “omega-hydroxylase P450 pathway” (GO:0097267, Padj = 0.008), “response to potassium ion” (GO:0035864, Padj < 0.001) and “inorganic anion transport” (GO:0015698, Padj = 0.006) (Fig. 4c). Details for all GO enrichments in BP and molecular function (MF) are provided for modules of interest in Supplementary Table 2.
Regarding the differentially expressed genes (DEGs) between the MNP conditions and the control conditions, a total of 438 DEGs were identified, with 88 DEGs vs. 0.025 µg L–1 MNPs and 405 DEGs vs. 1 µg L–1 MNPs; of these DEGs, 55 were common to both MNP conditions (Fig. 4d). A total of 24 and 64 up- and downregulated genes were found, respectively, from individuals treated under the 0.025 µg L–1 condition compared with the control and 246 and 159 up- and downregulated genes, respectively, from individuals treated under the 1 µg L–1 condition (Supplementary Fig. 7). Focusing on the tan module, a total of 82 DEGs (7 and 82 DEGs from 0.025 and 1 µg L–1 MNPs, respectively) were identified among the module genes, with 53 genes annotated according to UniProt entries (Supplementary Table 4). All DEGs were upregulated, and 7 DEGs were common to both MNP conditions, of which 6 were annotated (Fig. 4e and Supplementary Table 4). Furthermore, by focusing on shared DEGs among sample types (N = 5 DEGs; Fig. 5), a total of 4 DEGs were identified in the tan module (Fig. 4e and Supplementary Table 4), namely, CYP2C8, CYP2J2 (cytochrome P450 transcripts), HR4 (hormone receptor 4), and SULT1B1 (sulfotransferase family 1B member 1), as being specific to the 1 µg L–1 condition (Fig. 5). A complete list of DEG distribution among WGCNA modules in mantle samples is available in Supplementary Table 4.
Hemocyte module functional enrichment
A coexpression module was constructed by WGCNA of the expression values of 24,910 genes in 29 hemocyte samples. Overall, 12 modules were detected (Supplementary Fig. 4), and a total of 3 modules of interest were selected according to module-trait relationships (r ≥ 0.45, P ≤ 0.01), with 727, 1,020 and 2,092 genes in the magenta, pink and red modules, respectively (Supplementary Fig. 5). Among these 3 modules, none were related to condition, but the magenta module exhibited a strong correlation with 0.025 µg L–1 MNPs (r = − 0.67, P < 0.001) and 1 µg L–1 MNPs (r = − 0.52, P = 0.004), and there was a significant difference in ME expression between the control condition and both 0.025 and 1 µg L–1 conditions (Tukey’s HSD test, P = 0.022 and P = 0.017, respectively) (Supplementary Fig. 5). This module was significantly enriched for “maturation of LSU-rRNA from tricistronic rRNA transcript” (GO:0002108, Padj < 0.001), “cytoplasmic translation” (GO:0002181, Padj < 0.001), “ribosomal small subunit assembly” (GO:0000028, Padj < 0.001) and “establishment of protein localization to endoplasmic reticulum” (GO:0072599, Padj = 0.006). Details for all GO enrichments in BP and MF are provided for the magenta, pink and red modules in Supplementary Table 3.
A total of 204 DEGs were identified under MNP conditions compared to the control conditions, with 132 DEGs vs. 0.025 µg L–1 MNPs and 110 DEGs vs. 1 µg L–1 MNPs, of which 38 were common DEGs (Supplementary Fig. 5). A total of 73 and 59 up- and downregulated genes, respectively, were found in individuals under the 0.025 µg L–1 condition compared with the control condition, and 55 and 55 up- and downregulated genes, respectively, in the 1 µg L–1 condition compared to the control condition (Supplementary Fig. 7). A complete list of DEG distribution among WGCNA modules in hemocyte samples is available in Supplementary Table 5.
Pearl sac module functional enrichment
A coexpression module was constructed by WGCNA of the expression values of 28,116 genes in 29 pearl sac samples. Overall, 11 modules were detected (Supplementary Fig. 6), but no modules were selected according to module-trait relationships (r ≥ 0.45, P ≤ 0.01).
A total of 146 DEGs were identified in MNP conditions compared to the control condition, with 80 DEGs vs. 0.025 µg L–1 MNPs and 73 DEGs vs. 1 µg L–1 MNPs, with 7 common DEGs. We found a total of 22 and 58 up- and downregulated genes, respectively, in individuals treated under the 0.025 µg L–1 condition compared with individuals treated under the control condition and 44 and 29 up- and downregulated genes, respectively, in the 1 µg L–1 condition (Supplementary Fig. 6). A complete list of DEGs in the pearl sac is available in Supplementary Table 6.