PAHs and PCBs Affect Functionally Intercorrelated Genes in the Sea Urchin Paracentrotus lividus Embryos

Polycyclic aromatic hydrocarbons (PAHs) and polychlorinated biphenyls (PCBs) represent the most common pollutants in the marine sediments. Previous investigations demonstrated short-term sublethal effects of sediments polluted with both contaminants on the sea urchin Paracentrotus lividus after 2 months of exposure in mesocosms. In particular, morphological malformations observed in P. lividus embryos deriving from adults exposed to PAHs and PCBs were explained at molecular levels by de novo transcriptome assembly and real-time qPCR, leading to the identification of several differentially expressed genes involved in key physiological processes. Here, we extensively explored the genes involved in the response of the sea urchin P. lividus to PAHs and PCBs. Firstly, 25 new genes were identified and interactomic analysis revealed that they were functionally connected among them and to several genes previously defined as molecular targets of response to the two pollutants under analysis. The expression levels of these 25 genes were followed by Real Time qPCR, showing that almost all genes analyzed were affected by PAHs and PCBs. These findings represent an important further step in defining the impacts of slight concentrations of such contaminants on sea urchins and, more in general, on marine biota, increasing our knowledge of molecular targets involved in responses to environmental stressors.

Human activities represent a major source of stress for marine organisms. Anthropic activities can change, for example, the natural compositions, rate, deposition and transport of sediment, causing a radical shift of biological responses [18]. Since sediments produce the habitat for several benthic communities, their pollution could negatively influence the fitness and survival of marine flora and fauna [19]. Sediment drains and temporarily accumulates organic and/or inorganic compounds, including polycyclic aromatic hydrocarbon (PAHs), polychlorinated biphenyls (PCB) and organotin compounds (OTCs) derived from natural and anthropogenic sources, such as industrial, commercial, agricultural and urban activities [20][21][22][23]. Following several disruptive events (bioturbation and/or dredging), marine sediments can release these accumulated pollutants facilitating the re-allocation of contaminants within the water column [20,24,25]. This resuspension affects the marine environment and consequently human health, because pollutants can accumulate in bottom feeders and magnify through marine food webs, reaching humans through fishery products [26,27].
To our knowledge, there are several studies evaluating adaptability tolerance of sea urchins to several stress conditions, such as warming and ocean acidification [28][29][30]. However, concerning anthropogenic pollutants, the long-term exposure of sea urchins to heavy metals, PAHs and PCBs did not display any tolerance or adaptation [31][32][33][34].
PAHs and PCBs are among the most toxic pollutants commonly present in the sea. PAHs consist of a group of over 100 different organic compounds, constituted of two or more fused benzene rings and pentacyclic molecules, differing in their physical and chemical properties. PAHs can be released in the marine environment both through natural events and anthropogenic activities [35,36], with a noteworthy toxicity depending on the number of benzene rings [37]. PCBs are organic compounds, mainly derived from human activities, and considered among the most dangerous pollutants because of their toxic and carcinogenic properties, vast distribution and high persistence in the environment [37,38]. Several studies reported that PCBs were able to induce harmful effects on the embryo development of sea urchins [39][40][41]. Their direct involvement in the negative impact on both marine invertebrates' and vertebrates' embryo development was demonstrated by toxicity assays with pure compounds [42][43][44]. For example, the toxicity model proposed by Di Toro et al. [45], relied on assays using pure hydrocarbon compounds in dissolved phase, showing that dissolved PAHs were toxic in the absence of others compounds. However, few studies demonstrated possible negative influences of these compounds on the survival of adults and their progenies. Our recent study proved that PAHs and PCBs induced toxigenic effects on the sea urchin Paracentrotus lividus embryos derived from adults exposed to sub-chronic concentrations in contaminated mesocosm, using a multi-endpoint approach and de novo transcriptomic [33].
This investigation was aimed at further exploring the target genes involved in the response of the sea urchin P. lividus to PAHs and PCBs, and other environmental stressors. To perform this study, sea urchin plutei (48 h post-fertilization, hpf) were exposed to slight PAH-and PCB-contaminated mesocosms, also adopted in our previous experiments [33]. Firstly, an interactomic analysis was conducted on 62 genes, previously described in P. lividus [46][47][48][49] in order to understand the biological and functional relationships and the gene networks in which they were involved. Twenty-five new genes were identified, being functionally interconnected to 10 genes already described in previous studies [46][47][48][49]. Two networks were obtained, including most of the correlated genes, involved in stress responses, detoxification and developmental processes. Finally, the expression levels of these twenty-five genes were followed by Real Time qPCR to identify the possible gene targets affected by PAH and PCB.

Effects of PAHs and PCBs on Gene Expression by Real Time qPCR
The pollutants concentrations (i.e., as mean value, n = 3) in the spiked mesocosms at the end of the experiments were 5 µg/kg, 5 µg/L and 29.1 µg/L for PAHs in sediment, water and sea urchins (i.e., gonads), respectively; while PCBs concentrations were 5 µg/kg, 5 µg/L and < 2 µg/L in sediment, water and sea urchins (i.e., gonads), respectively. PAHs and PCBs were below the relative LOD in mesocosms.
PAHs and PCBs affected almost all the genes under analysis (Figure 1; fold changes are reported in Supplementary Table S1). At the pluteus stage, PAHs and PCBs had several common targets. Real-time qPCR at the pluteus (48 hpf) stages of the sea urchin P. lividus embryos, deriving from adults exposed for 2 months to sediments contaminated with PAHs and PCBs. Data are reported as a fold difference compared with control embryos, deriving from adults reared in tanks with noncontaminated sediments (mean ± SD). Histograms show the fold-changes of 25 genes involved in three functional processes: stress response, development/differentiation and detoxification. Fold differences greater than ±1.5 (see red dotted horizontal guidelines at values of +1.5 and −1.5) were considered significant. Real-time qPCR reactions were carried out in triplicate. Statistical differences were evaluated by nonparametric Mann-Whitney test. p-values < 0.05 were considered significant.

Figure 2.
Gene networks obtained by STRING interactome analysis of ten genes involved in stress response and detoxification. Correlation confidence score cut-off of 400 was reported. Among functionally correlated genes, those with up (red) and down (green) expression affected by PAHs (a) and PCBs (b) were reported. Color shading depends on fold-change values. Gray spheres represent additional connections. The list of human orthologous genes used for Network analysis is reported in Supplementary Table S2. The heat shock protein hsp90 clearly showed a key role in several biological processes since a huge number of additional connections (grey spheres) were observed ( Figure 2). TNF was strongly correlated to hsp90, hsp75, GAPDH, GST, PKS and NADH, while less connections have been observed for SULT1 and NADH, particularly with hasp90, hsp75, GAPDH and PKS (Figure 2). ChE and CYP2-UI showed only one connection to hsp90 through the involvement other associated genes ( Figure 2).
Among development and differentiation genes, a huge functional correlation was appreciable ( Figure 3). Several genes resulted extremely important within this biological cascade, including HH, CREB, JAK, STAT1, M-Vg1, Ptc, Smo, FZ-7, NLK, NOTCHLESS and CM-K, since a large number of connections was observed among them. In addition, these latter genes reported an independent association to other effectors, with NOTCHLESS and CM-K genes displaying the majority of them ( Figure 3). Concerning the less correlated genes, PLC and EGF were found functionally interconnected between each other and with other genes, such as JAK, STAT1 and CAMP, and also Ptc gene, in the case of EGF. Moreover, PLAUF3 revealed only three functional relationships with Smo, Ptc and STAT1 ( Figure 3).

Discussion
Our results greatly expand the reach of previous investigations on genes involved in the response of the sea urchin Paracentrotus lividus to PAH-and PCB-contaminated sediments [33]. Previous results revealed morphological malformations in the sea urchin plutei, which were correlated to a significant alteration of the expression of several genes involved in different functional processes, such as stress response, development/differentiation, detoxification and skeletogenesis.
Here, 25 new genes involved in stress response, development/differentiation and detoxification processes were isolated in P. lividus. Our findings showed that PAHs and PCBs were able to switch on almost all genes under analysis ( Figure 4; Supplementary  Table S3).
Thus, as already reported by Albarano et al. [33], a possible correlation could be hypothesized between morphological malformations of plutei and the molecular responses exhibited. In fact, among the genes analyzed, twenty-three were altered by both pollutants, revealing comparable values of gene expression in a huge number of shared targets ( Figure 1). These results could indicate that organic compounds affected some common molecular pathways, so changing physiological mechanisms in sea urchin adults, which in turn may induce abnormal offspring.
Concerning gene expression results, no significant differences were observed between PAHs and PCBs, suggesting that similar biological pathways were activated in response to PAHs and PCBs exposure. The only exception was GAPDH, whose expression increased in the PAHs treatment. Heatmap showing the expression profiles and hierarchical clustering of 25 genes analyzed by Real Time qPCR in embryos deriving from adult P. lividus sea urchins exposed for 2 months to sediment contaminated with PAHs and PCBs. Color code: green, up-regulated genes with respect to the control; red, down-regulated genes with respect to the control; black, genes for which there is no variation in gene expression with respect to the control.
As for the stress response (Figure 2), both PAHs and PCBs targeted all genes analyzed in the present study. Over the last 10 years it has been shown that alteration of ChE, PKS and GST activity may cause morphological alterations of the spicules, similar to those shown in Figure 4 [50][51][52]. In fact, skeletal malformations and delayed larvae correlated to significant down-regulation of PKS gene was recently found after a 24-and 48-h exposure of P. lividus eggs to nickel [53]. Similarly, CYP-2AI, hsp70 and hsp90 are involved in defense against chemical stressors, and their alteration is able to impact the gastrulation process [54][55][56]. For instance, malformed embryos and up-regulation of hsp90 was found in sea urchins exposed to thermal stress and bacterial LPS, corroborating a fundamental role of this gene in the stress response [57][58][59][60]. Concerning GAPDH, this gene possesses key functions in the glycolytic pathway, and also roles in nuclear RNA export, membrane fusion, and DNA repair [61]. Although GAPDH was considered as a housekeeping gene in several organisms, including sea urchins [62][63][64], in the present study this gene was found differentially expressed after PAH and PCB exposure. The alteration of the expression level of GAPDH gene in sea urchin could be considered a suitable bioindicator for human health since its up-regulation was detected in cervical and ovarian tissue during cancer development [65].
As reported for the stress response network, the affected developmental and differentiation genes were almost exclusively the same for both organic compounds. A strong up-regulation of these genes was induced by sea urchin exposure to PAHs and PCBs ( Figure 3). Ten functionally connected genes, HH, CREB, NOTCHLESS, JAK, NLK, EGF, PLC, M-Vg1, PLAUF3 and CM-K, were up-regulated in the same biological pathway (Figure 3), thus corroborating previous results reporting several aberrations in sea urchin progenies [33]. FZ-7 and STAT1 were not affected by PAHs treatment (Figure 3a), indicating that the biological cascade altered by PAHs and PCBs might be slightly different. Specifically, these genes were up-regulated only by PCBs treatment. As previously mentioned for GAPDH gene the variation of mRNA levels might also be correlated to human diseases. In fact, the overexpression of FZD7 was detected in 90% of tumors (most in the human hepatocellular carcinoma cells and in breast cancer cells) where causes the accumulation of β-catenin protein affecting the Wnt/β-catenin signal transduction pathway [66][67][68]. Similarly, the up-regulation of STAT1 represents a serious risk for human health since its activation is able to down-regulate programmed cell death genes in several cancer cell lines (i.e., breast cancer and adenocarcinoma cells) [69][70][71]. As reported in Figure 3, CREB, HH, PLAUF3, NLK, EGF, JAK, M-Vg, PLC and CM-K were up-regulated in both treatments. Similar to our study, the expression of the transcription factor CREB was targeted in response to heat stress and starvation in sea urchin of the genus Strongylocentrotus [72,73]. The PLC gene is typically regulated by CM-K, and together they are involved in the events following fertilization and embryo development in sea urchins [74][75][76][77]. The M-Vg gene together with nodal is involved into left-right axis formation, regulating some down-stream effectors, including the BMP2 gene [78]. The up-regulation of these genes could be a possible consequence of PAHs and PCBs exposures by inducing malformations in the embryo structure. The NOTCHLESS gene is involved in the NOTCH-DELTA pathway which, in turn, regulates the expression of NLK and Lefty genes [46,79,80]. This pathway has a key role during the induction and differentiation of secondary mesenchyme cells [80][81][82] and could be involved in the formation of malformed gastrulae. Ptc and Smo genes showed a down-regulation after both PAHs and PCBs treatments (Figure 3). The Smo and Ptc proteins are co-receptors of HH ligand and are expressed by the skeletogenic and non-skeletogenic mesoderm [83][84][85]. Probably, the defects in embryo development were caused by the perturbation of this pathway.
As described above, almost all genes were targeted by both pollutants. The sole exceptions were represented by FZ-7 and STAT1, which were up-regulated by PCBs (Figure 3b). FZ-7 receptor binding to Wtn6 controls the β-catenin signal during the specification of the endomesoderm [86]. The STAT1 gene together with JAK gene constitute the JAK/STAT pathway, which plays a fundamental role in the regulation of the cellular complex [87][88][89].
On the whole, molecular analyses of the aforementioned 25 genes revealed interesting results, which might justify the clear malformations affecting the apex and arms of sea urchin plutei obtained from adults exposed to PAHs and PCBs. These pollutants are released into the sea by anthropogenic activities and accumulated in marine sediments where sea urchins, typical hosts of benthic environments, could suffer for their long-lasting impacts.
In addition, an increase of knowledge on P. lividus gene functions has a great significance for molecular analysis on this well-established model organism. In fact, despite of its importance for the scientific community, the complete genome of the sea urchin P. lividus is still not fully annotated and this represents a limit in its use for molecular studies.
In the sea, when evaporation exceeds the contribution of rivers, the concentration of environmental pollutants increases with a consequent effect on organisms, biodiversity and as long-term effect on human health [90]. In fact, embryos and larvae of marine invertebrates seem to be a suitable indicator in understanding the toxicological response induced by chemical pollutants marine invertebrates, accumulating high levels of them in their tissues. Furthermore, these invertebrates have a key role in the pelagic as well in benthic food webs, representing intermediate consumers. Then, the toxicological risk faced by marine organisms and even by humans through the ingestion of contaminated edible species is concrete. Consequently, studies on the status of chemical pollutants in marine ecosystems represent an important step in evaluating possible risks on human health. On this line is also our approach with a 2-month exposure of the sea urchins to PAHs and PCBs, aiming at detecting the long-term morphological and molecular effects of these contaminants on marine invertebrates. More in general, our findings are in agreement with the idea that variations in gene expression represent an early indicator of the presence of stressful conditions in various marine environments. In fact, the identification of molecular pathways in which the targeted genes were involved represents a key step in understanding how marine organisms attempt protection against toxicants. In the case of PCBs and PAHs we found the alteration of almost the same target genes, revealing that both pollutants could activate similar biological pathways in sea urchin. This result might suggest the hypothesis that common responses of PCBs and PAHs could be also triggered in other marine invertebrates. Furthermore, this work increased the pool of genes named "defensome", as described in Marrone et al. [91], used by marine organisms to avoid deleterious consequences and irreversible damages. In conclusion, target genes for PAHs and PCBs may be considered as possible universal biomarkers to detect the presence and the effects of key environmental pollutants impacting the physiology of marine invertebrates.

RNA Extraction and cDNA Synthesis
Adult sea urchins were collected after 2 months of exposure in PAH-and PCBcontaminated mesocosms (192 µg/L and 0.15 µg/L, respectively), and their gametes were used for in vitro fertilization [33]. Collection of embryos at the pluteus stage (about 5000 sea urchin plutei) for RNA extraction was performed according to Ruocco et al. [93]. For each sample, 1000 ng of total RNA was retrotranscribed with an iScript™ cDNA Synthesis kit (Bio-Rad, Milan, Italy), following the manufacturer's instructions.

Network Analysis
Sixty-two genes [46][47][48][49] (see also Supplementary Table S4) were analyzed by Network-Analyst 3.0 software (https://www.networkanalyst.ca/; accessed on 8 February 2021; [94]), using STRING interactome of protein-protein interactions [95]. Human orthologous genes were used to compute the network analysis (Supplementary Table S2). The relations among genes (confidence score cut-off of 400) displaying experimental evidence were highlighted. Twenty-five mostly connected genes were then chosen and analyzed. The sequences were retrieved from the transcriptome of the sea urchin P. lividus deposited in the SRA database (Sequence Read Archive, available at https://www.ncbi.nlm.nih.gov/sra (accessed on 5 February 2021), accession number PRJNA495004, [96]; accession number SUB2817153, [97]; accession number SUB6701449, [33] and from NCBI (available at https://www.ncbi.nlm.nih.gov; accessed on 8 February 2021). For each gene, specific primers were designed on the nucleotide sequences (see Table 1). PCR fragments were purified from agarose gel using the QIAquick Gel Extraction kit (Qiagen, Milan, Italy), and the specificity of the PCR product was checked by DNA sequencing. PCR products were aligned to the original sequences by MultAlin Software [98] (see Supplementary Figure S1). Five serial dilutions were set up to determine Ct values and PCR efficiencies for all primer pairs (for RT-qPCR conditions see below). Standard curves were generated for each oligonucleotide pair using Ct values versus the logarithm of each dilution factor (Supplementary Figure S2). The biological functions for 25 genes are reported in Supplementary Table S3.

RT-qPCR Analysis
Undiluted cDNA was used as a template in a reaction containing a 0.3 mM final concentration for each primer and 1× FastStart SYBR Green master mix (total volume of 10 µL) (Applied Biosystems, Monza, Italy). PCR amplifications were performed in a ViiATM7 Real Time PCR System (Applied Biosystems, Monza, Italy) thermal cycler using the following thermal profile: 95 • C for 10 min, one cycle for cDNA denaturation; 95 • C for 15 s and 60 • C for 1 min, 40 cycles for amplification; 72 • C for 5 min, one cycle for final elongation; one cycle for melting curve analysis (from 60 • C to 95 • C) to verify the presence of a single product. Each assay included a no-template control for each primer pair. To capture intra-assay variability, all RT-qPCR reactions were carried out in triplicate. Fluorescence was measured using the ViiATM7 software (Applied Biosystems, Monza, Italy). The relative expression ratios were calculated from quantification cycles (Cq) through an efficiency (E) corrected calculation method (Etarget ∆Cq target (Mean Control-Mean Sample) /Ereference ∆Cq reference (Mean Control-Mean Sample) ) [99,100] using REST software (Version updated 2009, Relative Expression Software Tool, Weihenstephan, Germany). The expression of each gene was analyzed and internally normalized against ubiquitin [101] and 18S rRNA [102,103]. Relative expression ratios above 1.5 were considered significant. Nonparametric Mann-Whitney test was applied to ∆Cq (Cq gene of interest-Cq reference) values between treated and control samples (n = 3). p-Values < 0.05 were considered significant. Statistical analyses were performed using GraphPad Prism Software (version 9.00 for Windows, GraphPad Software, La Jolla, CA, USA, www.graphpad.com, accessed on 1 February 2021).
Fold change values were represented through a Heatmap generated by Heatmapper Software (http://www.heatmapper.ca/ (accessed on 5 February 2021); [104]). Values were scaled up by column and hierarchical clustering was performed on rows by a routine adopting the Average Linkage method considering the Euclidean Distances.

Gene Networks
An interactomic analysis was performed by NetworkAnalyst 3.0 software [93] available at https://www.networkanalyst.ca/ (accessed on 8 March 2021), using STRING interactome of protein-protein interactions [94]. Human orthologs of selected genes were used to compute the network analysis. The most significant relations among genes (confidence score cut-off = 900) displaying experimental evidence were highlighted.  Table S1: Fold changes in embryos at the pluteus stage, deriving from sea urchins exposed to PAHs and PCBs, (in red down-expressed genes; in green up-expressed genes) in comparison with the control (represented by embryos deriving from adults of sea urchins reared non-contaminated mesoscosm) at 48 hpf., Table S2: Gene name and acronym of genes from Paracentrotus lividus and their human orthologs., Table S3: Gene name, acronym, function and reference of 25 new genes., Table S4: Gene name, acronym and reference of 62 genes.