Reproductive toxicity induced by benzo[a]pyrene exposure: first exploration highlighting the multi-stage molecular mechanism in female scallop Chlamys farreri

Reproductive toxicity induced by benzo[a]pyrene (B[a]P) exposure has received great ecotoxicological concerns. However, huge gaps on the molecular mechanism still exist in bivalves. In this study, reproduction-related indicators were investigated in female scallops Chlamys farreri during life cycle of proliferative, growth, mature, and spawn stages, under gradient concentrations of B[a]P at 0, 0.04, 0.4, and 4 μg/L. Meanwhile, a multi-stage ovarian transcriptome analysis under 4 μg/L B[a]P exposure was also conducted to elucidate the potential molecular mechanisms. The results indicated that life-cycle exposure to 0.4 and 4 μg/L B[a]P significantly decreased GSI and sex steroid levels. Even 0.04 μg/L B[a]P could play the adverse role in DNA integrity at the mature and spawn stages. Ovarian histological sections showed that B[a]P inhibited the maturation and release of oocytes. Through the functional enrichment analysis of differentially expressed genes (DEGs) from transcriptome data, 18 genes involved in endocrine disruption effects, DNA damage and repair, and oogenesis were selected and further determined by qRT-PCR. The downregulation of genes involved in steroidogenic and estrogen signaling pathways indicated that B[a]P could cause endocrine disruption through both receptor-dependent and receptor-independent pathways. The variations of gene expressions involved in DNA single-strand break and repair implied the presence of toxic mechanisms similar with vertebrates. Additionally, the changes of gene expressions of cell cycle, apoptosis, and cell adhesion suggested that exposure to B[a]P possibly caused the reproductive toxicity effects by affecting oogenesis. Taken together, this study was a pioneer in combining genome-wide transcriptomic analysis with its corresponding reproductive indicators (GSI, sex steroid levels, DNA single-strand break, and histological sections) to explore the bivalves’ toxic mechanisms under B[a]P exposure. Meanwhile, some genes involved in estrogen signaling pathway and DNA damage were firstly analyzed in bivalves, and the expression data might be useful in establishing new hypotheses and discovering new biomarkers for marine biomonitoring.


Introduction
With the rapid pace of industrialization and frequent occurrence of marine oil spill accidents, persistent organic pollutants (POPs) are detected in terrestrial and marine environment globally (Heskett et al. 2012;Wurl and Obbard et al. 2004). Polycyclic aromatic hydrocarbons (PAHs) are one of the most widespread POPs in marine ecosystem, which could pose a potential threat to the health of aquatic organisms (Ankley et al. 2003;Kim et al. 2013). Sixteen PAHs have been listed as the priori contaminants by US Environmental Protection Agency (US EPA) due to their strong carcinogenic, teratogenic, and mutagenic properties (Office of the Federal Registration, 1982). During the last decades, numerous studies have observed that PAHs pollution could induce adverse effects on reproduction of marine organisms, including sex reversal, delayed gonadal development, and reduced spawning success (Chikae et al. 2004;Rochman et al. 2014;Vignet et al. 2016). For example, the success rate of reproduction and the survival rate of offspring have been observed to be reduced in maternal oyster Crassostrea gigas exposed to a representative PAH, benzo[a]pyrene (B[a]P), during ovarian development (Choy et al. 2007). Consequently, the phenomenon has attracted great concerns due to the strong evidence that PAHs possess high potential to induce reproductive toxicity in marine organisms. Currently, plenty of studies suggested that there is an association between PAHs exposure and reproductive toxicity in marine organisms, such as endocrine disruption (Collier et al. 2013;Yamamoto et al. 2017), DNA damage (Everaarts and Sarkar 1996;Meier et al. 2020), and histological alterations of gonads (Chukwuka et al. 2019;Louiz et al. 2009). However, information on the underlying toxic mechanisms of PAHs to marine invertebrates, such as mollusks (e.g., bivalves), is relatively scarce (Cuevas et al. 2015;Jing-jing et al. 2009;Sarker et al. 2018;Simão et al. 2020). Traditional studies, which investigate the toxic mechanism of POPs on reproduction of bivalves are based on the adverse effects of POPs as well as the transcriptional changes of well-known functional genes (Boulais et al. 2018;Yang et al. 2020aYang et al. , 2020b, leading to the remaining data gaps in this field . Hence, there is an urgent need to find a more efficient method to comprehensively elucidate the potential molecular mechanism of PAHs exposure on reproduction hazards in mollusks (Honda and Suzuki 2020).
The applications of omics in toxicological study have created opportunities for exploring the transcriptome of nonmodel species at unprecedented sensitivity and depth (Chen et al. 2019a;Zhang et al. 2010). Over the past decade, transcriptomic analysis has been well developed and widely introduced in reproductive toxicology studies in many aquatic species, such as fishes (Berg et al. 2016;Colli-Dula et al. 2018;Schiller et al. 2013), crustaceans (Liu et al. 2019;Yu et al. 2018), and bivalves Bachère et al. 2017). Recently, Meng et al. (2020) determined the toxic effects of benzophenones (BPs) on zebrafish larvae and reported detailed transcriptomic data to enhance the current understanding of the toxicity mechanisms. Overall, the integration of transcriptomic profiling and toxicological data can provide a deeper understanding of the mechanisms behind toxic processes.
Special attention was paid to the observation that animals at early stage were more sensitive when exposed to POPs for a given dose (Collier et al. 2013). Meanwhile, the whole life-cycle exposure of POPs has been proved to cause more severe impairment in reproductive functions of fish (Chen et al. 2019b;Horri et al. 2018;Timme-Laragy et al. 2006). However, majority of previous studies were short-term exposure experiments in bivalves, and only focused on the specific reproductive stages (mainly at mature stage) during gonadal development (Alonso et al., 2019;Smolarz et al., 2017). The limited researches have restricted our comprehensive knowledge on the reproduction-toxic effects of PAHs in bivalves.
Scallop Chlamys farreri is a typical bivalve for coastal monitoring programs, which is also frequently used for toxicological studies (Zhang et al. 2011). In the present study, female C. farreri was dosed with different concentrations of B[a]P (0, 0.04, 0.4, and 4 μg/L) during ovarian development stages (proliferative, growth, mature, and spawn stage), respectively. A multi-stage ovarian transcriptome profiling was performed to screen the key reproductive toxicology genes under 4 μg/L B[a]P exposure in bivalves. Overall, this study firstly revealed the multi-stage molecular toxic mechanisms by combining transcriptomic analysis and indicators of reproductive toxicity in bivalves. This would not only promote our understanding of the molecular toxicity of B[a] P to bivalves, but also provide vital data for future relevant ecotoxicological studies.

B[a]P preparation
B[a]P (CAS No. 50-32-8, 99% purity) and dimethyl sulfoxide (DMSO, 99% purity) were purchased from Sigma. B[a]P was first dissolved in DMSO, and then added to seawater to achieve a final DMSO concentration of 0.001%. The exposure concentrations of B[a]P at 0.04, 0.4, and 4 μg/L were chosen according to the reported levels in the coastal seawater and surface sediments of China and abroad (Sarria-Villa et al. 2016;Wang et al. 2016;Zhang et al. 2016).

Scallop treatment and collection
The experiment was conducted in strict accordance with the national and institutional guidelines for the protection of animal welfare, and this study was approved by the Committee on the Ethics of Animal Experiments at Ocean University of China. Female scallops (C. farreri) with the shell length of 5.8 ± 0.3 cm, mainly at pre-proliferative stage, were obtained from Nanshan aquatic market in Qingdao, China. Scallops were placed in 60-L aquaria with 40 L sand-filtered seawater (1 L water per scallop), with a continuously aerated 14-h dark/10-h light cycle at 15 ± 1 °C. The seawater in every tank was renewed completely once a day. Scallops were mainly fed with dried powder of Spirulina platensis, and the temperature, salinity, pH, and dissolved oxygen were daily monitored during the acclimatization period.
After 1 week's acclimatization, scallops were exposed to the three gradient concentrations of B[a]P at 0.04, 0.4, and 4 μg/L, and a B[a]P-free solution group was set as the control group. Each group was set in three replicates. Meanwhile, a multi-stage ovarian transcriptome analysis of B[a] P-free and 4 μg/L B[a]P exposure groups was conducted to elucidate the potential molecular mechanisms. Based on the ovarian histology by Chung (2008) of C. farreri, ovarian samples at four reproductive stages during the exposure experiment of 29 days were continually collected: proliferative stage (7 days, named Pro.); growth stage (14 days, named Gro.); mature stage (28 days, named Mat.); and spawn stage (29 days, named Spa.). The culture management was the same as the acclimatized period during the whole exposure experiment.

Gonadosomatic index (GSI) determination
The wet weight of the whole soft tissue and ovaries, of the control and three B[a]P experimental groups, was registered for GSI determination. The GSI was calculated as the formula of GSI = [wet weight of ovary (g) /wet weight of soft tissue (g)] * 100 (Frantzen et al. 2016).

Steroid hormones extraction and assays
The three steroid hormones of progesterone (P), testosterone (T), and 17β-estradiol (E2) levels in ovaries were analyzed with the commercial enzyme-linked immunosorbent assays kits (Lengton Biological Technology Co. LTD, China) after steroids extraction. Ovaries samples (0.1 g) were homogenized in PBS (1 × , 0.01 M, pH = 7.4) on ice for 10 min, then centrifuged at 3000 rpm for 20 min at 4 °C, finally the supernatants were collected carefully for further measurement. Concentrations of P, T, and E2 were measured by ELISA kit, following the manufacturer's instructions, then quantified with a MULTISKAN GO (Thermo Scientific) at 450 nm.

DNA alkaline unwinding assay
Under the pre-defined alkaline denaturing conditions, the conversion rate of double-stranded DNA (dsDNA) to singlestranded DNA (ssDNA) was proportional to the number of breaks in the phosphodiester backbone (Daniel et al. 1985), which was employed as an indicator of DNA integrity.
The method performed in this study was adapted from Ching et al. (2001) regarding DNA extraction and alkaline unwinding assay. Briefly, three equal portions for fluorescence determination of dsDNA, ssDNA, and alkaline unwound DNA (auDNA) were diluted from DNA samples, and the ratio between dsDNA and total DNA (F value) was determined using the following formula: F value = (auDNA-ssDNA) / (dsDNA-ssDNA).

Ovarian histological evaluation
Ovaries were excised and fixed in Bouin's solution (Solarbio, China). After 18 h of fixation, samples were dehydrated through a series of alcohol and xylene and then embedded in paraffin (m.p. 56-58 °C). Sections (4 μm) of the ovary were cut and stained by Ehrlich's hematoxylin and eosin (H&E) (Solarbio, China). Then, the ovarian sections were observed and photographed by the Nikon ECLIPSE Ci-L fluorescence microscope (Nikon, Japan) to determine the reproductive stages and analyze the structural damages caused by B[a]P.

Reproduction-toxicology genes exploration by RNA-Seq analysis
Ovarian samples of C. farreri from the control and high-dose B[a]P exposure of 4 μg/L, at the Pro., Gro., Mat., and Spa., were used for high-throughput RNA-seq, respectively. Each group consisted of three biological replicates. The RNA extraction, library construction, and RNA-Seq were performed as we described previously (Xu et al. 2020). And the clean data are available at the NCBI SRA under the accession numbers of SRP238487 and SRP268573.
The previous paper aimed to explore interference molecular mechanisms on lipid metabolism under B[a]P exposure in bivalves (Xu et al. 2020). At the present transcriptome analysis, genes with the threshold of p-value < 0.05 and |log2(FoldChange)|≥ 1 between control and its corresponding exposure group at each stage were identified as differentially expressed genes (DEGs). The identified DEGs were carried out into Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways functional enrichment analysis. GO enrichment analysis was implemented with the GOseq R package. The KOBAS (2.0) software package was used to test the statistical enrichment of DEGs in KEGG pathways, and the enrichment of KEGG Level 2 functions was to reflect the classification of KEGG functions. In addition, the top 20 of KEGG pathways (sorting by p-value for each entry) were further selected to create a bubble diagram. Consequently, based on the above functional enrichment analysis, key genes related to reproductive toxicology were selected for the following analysis.

Gene expression analysis
Total RNA for qRT-PCR was extracted using reagent RNAiso Plus (TaKaRa Corp., Dalian, China) with genomic DNA removed by RNAse-free DNase (TaKaRa Corp., Dalian, China). Reverse transcription of each RNA sample was performed to get the first-strand cDNA using Prime Script Reverse Transcriptase kit (TaKaRa, Dalian, China). Then, qRT-PCR analysis was used the SYBR Premix Ex TaqTM kit (TaKaRa, Dalian, China) on a PikoReal 96 Real-Time PCR System (Thermo Scientific., USA). The thermal profile for SYBR Green RT-PCR was 95 °C for 30 s, followed by 40 cycles of 95 °C for 10 s, then 51-55 °C for 30 s, 72 °C for 30 s, 65 to 95 °C in increments of 0.5 °C per 5 s.
Gene-specific primers were designed according to Illumina sequencing data with Primer Premier 5.0 software. The primer sequences are listed in Table 1. Three commonly used reference genes, β-actin, GAPDH, and Ribosomal 18 s, were tested to evaluate their stability as endogenous control gene by using GeNorm v3.5 software, and the M value of the test gene was 0.233, 0.286, and 0.277, respectively. Then, the relative quantification of each target gene was determined using β-actin as the internal reference gene. The expression levels of the target genes were calculated with the relative Ct method (Schmittgen and Livak 2008). Additionally, TBtools software was used to draw the heatmap of genes expression (Chen et al. 2020). And the relationship between RNA-Seq and qRT-PCR was calculated by R ggplot2 package of R software.

Statistical analysis
The results were presented as mean ± standard deviation (SD). SPSS (version 25.0, SPSS, Chicago, IL, USA) software was used for statistical analysis. The data were tested for normality (Shapiro-Wilk) and homogeneity of variance (Levene's), and the statistical differences between the three B[a]P exposure groups and their corresponding control group were evaluated by Dunnett's post hoc test of one-way analyses of variance (one-way ANOVA) (*, p-value < 0.05; **, p-value < 0.01).

Effects of B[a]P exposure on GSI
The GSI of female C. farreri at 0.4 and 4 μg/L B[a]P exposure groups decreased significantly at growth and mature stages. However, no significant difference occurred in 0.04 μg/L B[a]P exposure group during ovarian development process. It was noteworthy that the 4 μg/L B[a]P treatment significantly induced GSI at spawn stage compared with its control group (Fig. 1A). Figure 1B shows the effects of B[a]P exposure on P, T, and E2 levels in the ovary. The trends of the three sex steroids were similar throughout the experiment, that all of them continued to increase during ovarian maturation process but decreased at spawn stage. The 0.4 and 4 μg/L B[a]P significantly reduced P, T, and E2 levels at growth, mature, and spawn stages. However, 0.04 μg/L B[a]P only significantly inhibited P concentrations at mature stage and T levels at spawn stage, respectively. In addition, no significant difference in E2 levels was observed during the exposure period between control and 0.04 μg/L B[a]P exposure group. Overall, some obvious dose effects were found from the interference effects of B[a]P exposure on sex steroids, and the 0.4 and 4 μg/L B[a]P significantly interfered the synthesis of steroid hormones.

Effects of B[a]P exposure on DNA damage
The effects of B[a]P exposure on DNA alkaline unwinding are shown in Fig. 1C. The results reflected that the 0.4 and 4 μg/L B[a]P exposure induced serious DNA alkaline unwinding during the whole experimental period, while the 0.04 μg/L B[a]P played its negative role in the latter experiment (including mature and spawn stages).

Effects of B[a]P exposure on histological alterations
Ovarian histological sections of female scallops at four control groups displayed the normal characteristics of oogenesis during development process. At proliferative stage, the oogonia and previtellogenic oocytes were found in the ovary and scattered on the follicular wall ( Fig. 2A). Figure 2E shows that both the number and size of oocytes increased and filled the follicles gradually, indicating the coming of growth stage. C. farreri developed into mature stage when ovarian follicles were filled with polygonal-shaped mature oocytes dense packing (Fig. 2I). In this study, C. farreri developed into the first spawn stage, and the ovarian follicles showed obvious cavity with the discharge of mature oocytes, while the immature oocytes and a small number of mature oocytes still remained in the follicles (Fig. 2M). Ovarian development of female scallops was severely affected by B[a]P exposure. At growth stage, exposure of 0.4 and 4 μg/L B[a]P induced the disorganization of follicular structure, while oocytes showed signs of degeneration with cell lysis, and disruption of the nuclear membrane ( Fig. 2F, G, H). At mature stage, exposure of the three B[a]P groups altered the normally polygonal mature oocytes, disordered oocytes, and delayed ovarian development. The severe vacuolization of the cytoplasm occurred at 4 μg/L B[a]P exposure group, with cell lysis and infoldings of the nuclear membrane ( Fig. 2J, K, L). After spawning, interspace among follicles occurred and some degenerated oocytes were observed in three B[a]P exposure groups (Fig. 2N, O, P). In summary, the ovarian development was influenced by the passage of B[a]P exposure time and the increase of exposure dose, respectively.

Identification of DEGs
The overview of transcriptome sequencing results was performed before (Xu et al., 2020). In general, the RNA-seq results indicated that the sequence data were of high quality and appropriate for following bioinformatic analysis. The Venn diagram of Fig. 3A shows that there were totals of 2986 DEGs generated from the four comparisons in this study. Specifically, 754 genes at proliferative stage, 960 genes at growth stage, 1708 genes at mature stage, and 495 genes at spawn stage were differentially expressed at 4 μg/L B[a]P exposure condition compared with the control groups, respectively.

Functional annotation and enrichment analysis of DEGs
Functional classifications of the total 2986 DEGs were performed via GO and KEGG pathway annotation analysis. Figure 3B shows the top 30 GO terms at the second level of GO annotation analysis. For molecular function (MF) category, catalytic activity (with 894 DEGs) and binding (with 843 DEGs) terms were the most dominated DEGs. In addition, cellular process and metabolic process with 553 DEGs and 388 DEGs, respectively, had the most DEGs in biological process (BP) category. What's more, membrane part with 1039 DEGs in cellular component (CC) was the most enriched GO term, which indicated the intense toxic effect to membrane under B[a]P exposure. Figure 3C reveals the histogram of KEGG annotation results of the total 2986 DEGs. Based on the first classifications of KEGG, six types of KEGG pathways were annotated, which included "Metabolism" (with 428 DEGs), "Genetic Information Processing" (with 65 DEGs), "Environmental Information Processing" (with 245 DEGs), "Cellular Processes" (with 248 DEGs), "Organismal System" (with 521 DEGs), and "Human Diseases" (with 662 DEGs). Overall, "Human Diseases," "Organismal System," and "Metabolism" dominated the top 3 KEGG classifications with enriched DEGs. Additionally, "Signal transduction" in "Environmental Information Processing" and "Cancers: Overview" in "Human Diseases" had the most DEGs. The results of KEGG annotation indicated the destructions of organism homeostasis under B[a]P exposure, and the carcinogenesis of B[a]P was confirmed in the ovaries of C. farreri.
To further elucidate the detailed pathways and genes linked to reproductive toxicity at each specific stage caused by B[a]P, the integrated top 20 KEGG level 2 classification was shown in a bubble diagram (Fig. 4). Figures A, B, C, and D represent the KEGG enrichment diagrams of DEGs generated from the four reproductive stages (proliferative, growth, mature, and spawn stages) under B[a]P exposure, respectively. For "Metabolism" category, it was notable that "Steroid hormone biosynthesis" pathway was enriched at the four comparisons. For "Environmental Information Processing" category, "Cell adhesion molecules (CAMs)" pathway was enriched in the early stage of ovarian development. For "Cellular Processes" category, "Lysosome" pathway was significantly enriched along with numerous DEGs, and "Oocyte meiosis" pathway of cell growth was significantly enriched at spawn stage. Furthermore, signaling pathways that were closely related to oogenesis and ovarian development, including "Estrogen signaling pathway," "Thyroid hormone signaling pathway," "Insulin signaling pathway," "GnRH signaling pathway," "Oxytocin signaling pathway," and "PPAR signaling pathway," were prominently enriched under B[a]P exposure. In addition, DEGs were enriched in "Single strand break repair (SSBR)" pathway.

Screening of reproduction-toxicology genes under B[a]P exposure
In this report, the DEGs related to endocrine disruption effects, DNA damage, and oogenesis induced by B[a]P were mainly focused on. According to the above functional enrichment The asterisks indicated statistically significant differences when comparing to the corresponding control group (* p-value < 0.05; ** p-value < 0.01). Data are given as mean ± SD of three replicates ◂ analysis, a total of 30 DEGs were screened with their transcriptional changes and are listed in Table 2. These 30 DEGs were grouped into three parts, involving (1) endocrine disruption effects, which includes receptor-independence pathway as well as receptor-dependence pathway, (2) DNA damage and repair, and (3) oogenesis, which includes ell cycle, apoptosis, and cell adhesion.
Additionally, correlations of 18 candidate DEGs between RNA-Seq and qRT-PCR were assessed with Pearson's coefficient, which indicated that the RNA-Seq results were highly correlated with qRT-PCR and the RNA-Seq results were reliable (regression p-value < 0.001, Pearson's correlation coefficient r = 0.853, Fig. 5).

Analysis of expression patterns of endocrine disruption effects genes
Analysis of expression patterns of steroidogenesis genes: receptor-independence pathway The effects of B[a]P exposure on the mRNA expression of steroidogenesisrelated genes were determined (Fig. 6A). At proliferative stage, CYP17 was significantly downregulated at three B[a] P exposure groups compared with the control group, and  The top 30 GO terms of the total 2986 DEGs at the second level of GO annotation analysis among molecular function, cellular component, and biological process. (C) The histogram of KEGG annotation of the total 2986 DEGs among its first classifications, which include "Metabolism," "Genetic Information Processing," "Environmental Information Processing," "Cellular Processes," "Organismal System," and "Human Diseases". Abbreviation: Con, control group; Exp, exposure group; Pro, proliferative stage; Gro, growth stage; Mat, mature stage; and Spa, spawn stage the differential expression showed the dose effects, while the expressions of 3β-HSD and 17β-HSD only significantly responded to 4 μg/L B[a]P exposure group. At growth stage, the mRNA expression of GnRHR showed significant downregulation at three B[a]P exposure groups, which had the dose effect relations with the exposure concentration. However, the expression levels of CYP17 and Foxl2 were significantly downregulated with their control groups at 0.4 and 4 μg/L B[a]P groups, while 3β-HSD and MAPKAPK5 only showed significant differences at 4 μg/L B[a]P group. The most differential expression of genes occurred at mature stage, that is, the expression levels of CYP17, 3β-HSD, and 17β-HSD significantly responded to the three B[a]P exposure groups, while the expressions of GnRHR, MAP-KAPK5, and Foxl2 were significantly downregulated at 0.4 and 4 μg/L B[a]P exposure groups. At spawn stage, except for CYP17 and 17β-HSD that showed no significant difference with the control groups, the other four genes all significantly responded to B[a]P exposure. The expression levels of GnRHR, 3β-HSD, and Foxl2 only showed significant downregulation at 4 μg/L B[a]P group, and MAPKAPK5 mRNA expression was significantly downregulated at 0.4 and 4 μg/L B[a]P groups. Notably, the expression levels of , Con_Mat_vs_Exp_Mat (C), and Con_Spa_vs_ Exp_Spa (D) DEGs, respectively. KEGG diagram description: on the left of every chart, the red-filled square with red text (a) represented the pathways of "Metabolism"; the green-filled square with green text (b) represented the pathways of "Environmental Information Processing"; the yellow-filled square with yellow text (c) represented the pathways of "Cellular Processes"; the purple-filled square with purple text (d) represented the pathways of "Organismal Systems"; the pinkfilled square with pink text (e) represented the pathways of "Human Diseases." Additionally, the rectangular at right represented the saliency, and the colors from blue to red represent the decrease of saliency with p-value, which indicated from 0 to 0.3 (A), from 0 to 0.2 (B), from 0 to 0.04 (C), and from 0 to 0.2, respectively. Larger bubbles indicate a higher number of DEGs Foxl2 were significantly lower compared with the control group among all three B[a]P exposure groups at spawn stage.

Analysis of expression patterns of genomic and non-genomic pathways genes: receptor-dependence pathway
The effects of B[a]P exposure on mRNA expression of genes responsible for estrogen signaling pathways in female scallop C. farreri were determined (Fig. 6B). Expressions of ER and Vtg shared the same pattern: when exposed to 4 μg/L B[a]P, the expressions were significantly reduced compared to their corresponding controls during the whole experiment period. Additionally, 0.4 μg/L B[a]P exposure group significantly downregulated the expression of ER at growth and mature stages (Fig. 6B1).
The changes of Cav-1 and Cav-3 mRNA expressions indicated the non-genomic interference effect caused by B[a]P (Fig. 6B2). For 4 μg/L B[a]P exposure group, the trends of Cav-1 and Cav-3 mRNA expressions were different at proliferative stage, while they were all upregulated at growth stage. Cav-1 was sensitive to B[a]P exposure at spawn stage

Analysis of expression patterns of DNA single-strand break repair genes
Topo1 and XRCC1 were selected to represent the DNA alkaline unwinding and its repairment, respectively (Fig. 6C). The 0.4 and 4 μg/L B[a]P treatments significantly downregulated Topo1 mRNA expression at growth and mature stages. Conversely, the expression of XRCC1 was significantly upregulated at proliferative stage. However, at growth stage, the three B[a]P exposure groups all significantly regulated the mRNA expression of XRCC1, which indicated the interference of DNA repairment caused by B[a]P. Fig. 6 Heatmap showed the transcriptional levels of genes related to reproductive toxicology effects caused by B[a]P exposure in female C. farreri. I: genes involved in endocrine disruption effects, including steroidogenesis (A) and estrogen signaling pathway (B), further divided into genomic pathway (B1) and non-genomic pathway (B2). II: genes involved in DNA single-strand break and repair (C). III: genes involved in oogenesis, which include cell cycle (D1), apoptosis (D2), and cell adhesion (D3). The rectangular at left represented the transcriptional levels of genes, while the red indicated the upregulation and the blue indicated the downregulation. The asterisks indicated statistically significant differences when comparing to the corresponding control group (* p-value < 0.05; ** p-value < 0.01). Data are given as mean ± SD of three replicates  Figure 6D shows the effects of B[a]P exposure on the mRNA expressions related to oogenesis. The results indicated that B[a]P significantly influenced the expressions of CDK7 and CCNC, which were related to cell cycle, except at growth stage (Fig. 6D1), whereas the expressions of FasL and Caspase-8 that involved in apoptosis were upregulated by B[a]P exposure during the latter experiment (at mature and spawn stages) (Fig. 6D2). For cell adhesion, profilin-2 and COL6A3 mRNA expressions showed cooperatively significant downregulation under B[a]P exposure during the whole experiment (Fig. 6D3).

Discussions
The publishment of C. farreri genome represents the coming of post-genomic era, while there remains poor exploration from the genomic perspective (Li et al. 2017;. Additionally, the application of genome-wide transcriptomic analysis combined with relevant toxic data would greatly improve the understanding of the underlying toxic mechanisms (Meng et al. 2020;Suman et al. 2020). To our best knowledge, this is the first study to synthetically highlight endocrine-disrupting effect, DNA single-strand break repair, and oogenesis induction in bivalves under B[a] P exposure during the whole life cycle. Importantly, this study is a pioneer to take advantage of transcriptomic analysis and its corresponding toxic data to explore the underlying toxic mechanisms under B[a]P exposure in bivalves. The schematic representation of the possible mechanisms underlying B[a]P-induced reproductive toxicology is displayed at Fig. 7, and the specific functions would be discussed in the following part.

Exploration of the toxic mechanism of endocrine-disrupting effect under B[a]P exposure
Strong evidences have indicated that B[a]P could act as an endocrine-disrupting chemical (EDC) by interfering with the synthesis, secretion, and metabolism of steroid hormones in aquatic animals (Hoffmann and Oris 2006;Inyang et al. 2003). Numerous studies have proved that the steroidogenic and metabolic pathways in mollusks are similar with vertebrates, and it has been proposed that the vertebrate-like steroids P, T, and E2 also dominated the sexual maturation and spawn of mollusks (Fernandes et al. 2011;Köhler et al. 2007). In this study, KEGG pathway of "steroid biosynthesis" was enriched throughout the four reproductive stages, revealing the endocrine-disrupting effects at transcriptome level. In addition, the significant decrease of concentrations of P, T, and E2 directly confirmed the results from the biochemical level. This was consistent with the previous study by Gao et al. (2018), which found that both E2 and T levels were reduced in the ovary of adult female zebrafish under embryonic exposure to B[a]P. Steroid hormones play a crucial role in the regulation of ovarian processes, including follicle growth, oocyte maturation, and spawning (Craig et al. 2011). Thus, the decrease of level of sex steroid hormones could be a crucial underlying reason for postponed gonadal development in C. farreri under B[a]P exposure. It has been well characterized that steroidogenesis was regulated by the upstream protein kinases (Ras/Raf; MAPK/ ERK), Adcy-PKA pathway, and its downstream enzymes (Manna et al. 2009). According to the transcriptome analysis in this study, the collaborative downregulation of steroidogenesis-related genes (upstream: PKA,EGFR,GnRHR,and MAPKAPK5;downstream: CYP17, probably demonstrated the interference mechanism of B[a]P on steroidogenesis in bivalves (Fig. 7A1). It was noteworthy that B[a]P reduced the expression of Dax1 during ovarian development, which has been identified to repress the transcriptional activity of steroidogenic factor (SF-1) (Crawford et al. 1998). Deng et al. (2010) reported that chronic TBP exposure reduced the transcription of steroidogenic genes (3β-HSD, 17β-HSD, CYP17, CYP19A, CYP19B) in the ovary of zebrafish. However, CYP19 has not been sequenced in molluscan genome yet (Scott 2013). It is interesting that CYP19 has been identified under the temporal and spatial regulation of gene Foxl in fish, and both Foxl and CYP19 could participate in the regulation of sex steroid synthesis (Hudson et al. 2005;Nagahama 2005). Therefore, besides the widely discussed enzymes involved in steroidogenesis, some new genes such as Dax1 and Foxl2 should be paid attention to when the B[a]P-induced endocrine-disrupting effects were evaluated in mollusks (Fig. 7A2).
By mimicking endogenous estrogen, EDCs could bind to ERs and disturb estrogen signaling through genomic and non-genomic pathways. The consistent downregulation of nucleus-ER and Vtg mRNA expression under 4 μg/L B[a]P exposure condition indicated the genomic pathway caused by EDCs (Roepke et al. 2005). However, recent evidence indicated that the new non-genomic pathways also played crucial roles in estrogen signaling, where EDCs could bind to and activate some non-canonical estrogen receptors, thus leading to rapid stimulation in various protein kinase pathways (Lee et al. 2013;Rosenfeld et al. 2019). Contradictorily, the verified estrogen receptors in non-genomic pathway, such as membrane-ERs and G protein-coupled receptors, have not been sequenced in molluscan genome yet (Gomesdos-Santos et al. 2020). However, it has been demonstrated that membrane-ERs could locate in the caveolae, which was the invagination of plasma membrane vesicles (Razandi et al. 2002;Stan 2005). Blalock et al. (2018) reported that EDCs could bind to membrane-ERs, leading to the cell The downstream genes of estrogen signaling pathway involved in oogenesis, which include cell cycle, apoptosis, and cell adhesion signaling cascade, and subsequently have impact on cell cycle in exposed mussels Mytilus edulis. Taken together, the identified non-genomic estrogen signaling pathway-related genes (Cav-1, Cav-3, and CALM) revealed by transcriptomic analysis provided powerful evidence for the further study on B[a]P-induced endocrine-disrupting effects in bivalves (Fig. 7B).

Exploration of the toxic mechanism of DNA single-strand break repair under B[a]P exposure
DNA damage has been used as the biomarker under POPs exposure in marine bivalve monitoring program (Martins et al. 2013;Mitchelmore and Chipman 1998). It has been shown that exposure to 10 μg/L B[a]P could cause serious DNA single-strand break (SSB) in the ovary of C. farreri , and the method of DNA alkaline unwinding assay has been widely used to reflect DNA damage levels in bivalves (Slobodskova et al. 2012). The significant DNA SSB was observed in all three B[a]P exposure groups during the development process in the present study. Banni et al. (2010) reported that B[a]P-induced DNA damage degree in mussel Mytilus galloprovincialis was positively correlated with the exposure time. When exposed to benzo[k]fluoranthene, the DNA damage in the liver of common carp (Cyprinus carpio) increased with the exposure concentration (Kim et al. 2008). This study further confirmed that DNA damage could be a sensitive biomarker to evaluate marine pollutant stress during the life-cycle exposure in mollusks.
DNA damage is a major threat to genetic disorders, which may result in cancers, and DNA single-strand break repair (SSBR) is critical for the survival and genetic stability of mammalian cells (Tubbs and Nussenzweig 2018). However, data on the toxic mechanism of DNA damage in bivalves are limited ). It has been reported that SSB could be arisen either directly (e.g., excessive ROS production) or indirectly (e.g., inactivation of cellular enzymes Topo1) (Whitehouse et al. 2001). Topo1 is important in promoting the formation of DNA circular double-stranded molecules (Plo et al. 2003). Srinivasan et al. (2001) reported that PCBs increased DNA SSB in human HL-60 cells and furtherly found that PCBs could inhibit the Topo activity (Srinivasan et al. 2002). The significant decrease of Topo1 mRNA expression found in this study implied the potential disturbing mechanism of DNA damage under B[a]P exposure.
The knowledge of DNA repairment in aquatic species is of great importance as it is the primary line of defense against the gene-toxicants (Kienzler et al. 2013). SSBR can be divided into four basic steps: (1) DNA damage binding, which is achieved by the ADP-ribose polymerases, such as PARP-1 or PARP-2 (Langelier et al. 2013); (2) DNA end processing, which is stimulated by the interaction with XRCC1 (Ali et al. 2009); (3) DNA gap filling with the DNA polymerase Polβ (Sugo et al. 2000); and (4) DNA ligation with the involvement of DNA ligase (LIG) (Cappelli et al. 1997). Zuo et al. (2012) found that the expressions of both LIG1 and LIG3 were affected in the liver of fish Kryptolebias marmoratus after treatment with tributyltin (TBT). Surprisingly, the SSBR has been poorly investigated in aquatic species despite its key role in repair of oxidative damage (Kienzler et al. 2013). The results in this study suggested that there might exist similar toxic mechanism of DNA single-strand break and repair under B[a]P exposure in mollusks, while more research should be conducted to verify this hypothesis.

Exploration of the toxic mechanism of oogenesis under B[a]P exposure
Bivalve mollusks are descendants of an early-Cambrian lineage and have successfully evolved unique strategies for reproduction (Geist 2010;Regan et al. 2021). Normal ovarian growth contributes to fertilization and improvement of zygote quality substantially, while manipulating the development, maturation, and release of oocytes are of great significance to ensure the success of molluscan reproduction (Sühnel et al., 2014). However, oocyte is especially susceptible to the adverse effects of environmental toxicants, as well as EDCs (Gregoraszczuk and Ptak A 2013). The retardation of oogenesis and oocytes resorption was observed in the ovary of scallop (Mizuhopecten yessoensis) from the bay with heavy organic pollutions (Vaschenko et al. 1997). Tian et al. (2015) and Yang et al. (2020a) found that B[a] P exposure resulted in a significant decrease in the number of mature sperms and oocytes, respectively. Loughery et al. (2018) reported that the female GSI significantly decreased in the fathead minnow (Pimephales promelas) exposed to 202 μg/L phenanthrene (Phe) for 7 weeks. Interestingly, higher GSI and more oocytes occurred at 4 μg/L B[a]P group at spawn stage in this study (Fig. 1A, Fig. 2D). And the results indicated that B[a]P might induce anti-estrogenic effects in female scallops (Gao et al. 2014).
Numerous studies have reported the reproductive disorders under B[a]P exposure, particularly emphasizing damage on reproductive cell development and induction of serious apoptosis (Asweto et al. 2017;Kim et al. 2017). What's more, apoptotic rate of germ cell has been used as a reproductive biomarker of contaminant exposure in male fish (Zoarces viviparus) (Lyons et al. 2004). In this study, delayed and abnormal oogenesis was observed during ovarian mature process, with an increasing malignity order depending on dose effect. Additionally, the KEGG functional analysis indicated that "Oocyte meiosis" pathway was significantly enriched at spawn stage under B[a]P exposure. Overall, it could be concluded that B[a]P exposure seriously affected the normal oogenesis and ovarian functions in C. farreri. Meanwhile, genes involved in cell cycle (CDK7 and CCNC), apoptosis (FasL and Caspase-8), and cell adhesion (Profilin-2 and COL6A3) were screened and validated. The mRNA expressions of oogenesis-related genes exhibited the possibly disturbing mechanisms under B[a]P exposure, and the data might be useful in discovering new biomarkers to detect reproductive toxicity caused by B[a]P in bivalves (Fig. 7C).

Conclusion
This study explored the molecular mechanisms of reproductive toxicity induced by B[a]P in female scallop C. farreri during ovarian developmental process. The changes of reproduction-related indicators including GSI, sex steroid levels, DNA damage, and ovarian histology indicated the serious toxic effects caused by B[a]P. The expressions data of genes screened by transcriptome analysis elucidated the potential molecular mechanisms of reproductive toxicology under B[a]P exposure. The downregulation of steroidogenic and estrogen signaling pathways genes emphasized the functions of receptor-independent and -dependent pathways under B[a]P exposure. The variation of DNA single-strand break and repair gene expressions implied that there might exist the similar toxic mechanism with vertebrates. In addition, genes expression data involved in cell cycle, apoptosis, and cell adhesion exhibited the possibly toxic effects of oogenesis caused by B[a]P. In summary, this study promotes the current understanding of molecular toxicity of B[a]P in bivalves and provides vital data for future relevant ecotoxicological studies.
Funding This work was supported by the Key Research Project of Shandong (2018) (2018GHY115007).
Data availability All data generated or analyzed during this study are included in this published article.

Declarations
Ethics approval and consent to participate All experimental procedures were conducted in conformity with institutional guidelines for the care and use of laboratory animals, and protocols were approved by the Institutional Animal Care and Use Committee in Ocean University of China, Qingdao, China.
Consent for publication Additional informed consent was obtained from all individual participants for whom identifying information is included in this article.

Competing interests
The authors declare no competing interests.