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.

1. Introduction

Marine organisms are permanently exposed to multiple stressors, such as climate changes [1,2,3] and the consequential ocean acidification, deoxygenation and sea-level rise [4,5,6], natural toxic metabolites [7,8,9,10,11] and compounds deriving from human activities [12,13,14,15]. The exposure to these stressors induces marine organisms to adopt strategies against external or internal environmental factors through metabolic and/or molecular changes to maintain cellular homeostasis [16,17].

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.

2. Results

2.1. 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.

Figure 1

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 non-contaminated 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.

Stress genes:

All nine genes were targeted by PAHs and PCBs. Both contaminants increased the expression levels of three genes (TNF, GST and CYP-2UI) and decreased GAPDH, PKS, ChE, SULT1, hsp90 and hsp75.

Genes involved in development/differentiation:

Out of 15 genes analyzed, only STAT1 and FZ-7 were not targeted by PAHs, whereas these two genes were up-regulated by PCBs. Common molecular targets switched on by PAHs and PCBs were CM-K, JAK, PLAUF3, M-Vg1, NOTCHLESS, PLC, EGF, NLK, HH and CREB which were up-regulated; whereas Lefty, Ptc and Smo were down-regulated.

Genes involved in detoxification:

The only gene analyzed, NADH, was targeted by both PAHs and PCBs, which were able to decrease its expression levels.

2.2. Network Analysis

Twenty-five genes analyzed by Real Time-qPCR were functionally intercorrelated among them and also correlated with several genes previously isolated from Paracentrotus lividus. Depending on gene functions, they were ascribed to two functional networks, reporting the highest correlations among stress/detoxification and developmental genes.

Among the 10 genes involved in stress responses, hsp90, hsp75 and GAPDH shared the largest number of connections (Figure 2).

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).

Figure 3

Gene network obtained by STRING interactome analysis of 15 genes involved in developmental processes. Correlation confidence score cut-off of 400 was reported. Among functionally correlated genes, those with up (red), down (green) and unchanged (blue) expression affected by PAHs (a) and PCBs (b) were reported. Color shading depends on fold-change values. Gray spheres represent additional connections. The whole list of human orthologous genes used for Network analysis is reported in Supplementary Table S2.

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).

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).

Figure 4

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.

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.

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.

4. Materials and Methods

4.1. Experimental Conditions

Adult Paracentrotus lividus (seven females and three males) were reared in each of experimental tanks of a mesocosm, previously spiked with PAHs (acenaphthene (ACE), acenaphthylene (ACY), anthracene (ANT), benzo(a)anthracene (BaA), chrysene (CHR), fluoranthene (FLT), fluorine (FLR), phenanthrene (PHE), pyrene (PYR)), and PCBs (standard solution) [33]. Each mesocosm was spiked with 192 µg/L and 0.15 µg/L of PAHs and PCBs (i.e., nominal concentrations), respectively, to investigate any sub-chronic effect at concentrations below sediment threshold limit values (TLVs) (TLVPAHs = 900 µg/L and TLVPCBs = 8 µg/L) according to the Italian regulation D.M. 173/2016. To evaluate the compartmentalization of PAHs and PCBs in sediment and seawater, their concentrations were evaluated at three times (before (t0) and at the end (tf) of the experiment) and were measured according to Trifuoggi et al. [92]. The limit of detection and limit of quantification for PAHs and PCBs were: LODPAHs = 0.004 µg/L and LOQPAHs = 0.01 µg/L, LODPCBs = 0.002 µg/L and LOQPCB = 0.05 µg/L for the seawater samples; LODPAHs = 0.016 µg/Kg and LOQPAHs = 0.01 µg/Kg, LODPCBs = 0.03 µg/Kg and LOQPCBs = 0.01 s µg/Kg for the sediment; LODPCBs = 0.4 µg/kg w.w. and LOQPCBs = 2 µg/kg w.w. for sea urchin tissues. All details about methods for PAHs and PCBs chemical assessment were summarized in Trifuoggi et al. [92].

4.2. RNA Extraction and cDNA Synthesis

Adult sea urchins were collected after 2 months of exposure in PAH- and PCB-contaminated 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.

4.3. Network Analysis

Sixty-two genes [46,47,48,49] (see also Supplementary Table S4) were analyzed by NetworkAnalyst 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).

Table 1

Gene Name, Acronym, Accession Numbers, Primer Names, Primer Sequences and Amplicon Lengths of PCR Fragments.

Gene Names	Acronym	Accession Number	Primer	Sequence (5′→3′)	Amplicon Lenght (bp)
Calcium/calmodulin-dependent protein kinase type 1D	CM-K	PRJNA495004	Pl_CM_F1	GTTATCCTCCATTTTACGATGAG	168
			Pl_CM_R1	GCAGATATACGTGTGAGGAAG
Camp-responsive element	CREB	PRJNA495004	Pl_CREB_F2	GTAACTAAAGCATCTGGGAGAC	158
			Pl_CREB_R2	GGTTCAGATATTAGTGGATGC
Cholinesterase	ChE	SUB6701449	Pl_ChE_F2	CGAGATGGCGTATGTTTTGAG	160
			Pl_ChE_R2	GACTATGTTCCCGCTGACTG
Citochrome P450 2UI isoform X2	CYP-2UI	SUB6701449	Pl_CYP-2UI_F1	GCGCCTCTTCGTTCTATTCC	174
			Pl_CUP-2UI_R1	CGGCATAGTAGTAGACTAGC
Epidermal growth factor	EGF	PRJNA495004	Pl_EGF_F1	CGGCGGTGTGTGTATCGATG	189
			Pl_EGF_R1	CAGTAGCCATCCTAGTGTTCC
Frizzled7	FZ-7	PRJNA495004	Pl_FZ_F1	GATCGTGAGCGTAGCATATAC	175
			Pl_FZ_R1	CATGGTCTTTTTGGGCACTA
Glutathione-S-transferase	GST	SUB6701449	Pl_GST_F4	GCCCGACTTACCTACTTTGC	165
			Pl_GST_R4	CTTGCAGCTCATCACTGATGG
Glyceraldehyde-3-phosphate dehydrogenase	GAPDH	PRJNA495004	Pl_GAPDH_F1	GTACTACTTCTCATTCACCTTG	213
			Pl_GAPDH_R1	CATAGCTCTGACACCGCCAC
heat shock protein 75	hsp75	PRJNA495004	Pl_hsp75_F2	GGACTGGTGGAACAACTATATC	173
			Pl_hsp75_R2	CGATCACCACTCTCTGTCAC
heat shock protein 90	hsp90	SUB6701449	Pl_hsp90_F1	GGGTGTGGTAGATTCTGATG	148
			Pl_hsp90_R1	GCTCTCCATGTATTCATCAG
Hedgehog	HH	PRJNA495004	Pl_HH_F1	GGTACATGAGGCACAAGCTAG	193
			Pl_HH_R1	CCACTTCACATCACTTGACC
Janus kinase	JAK	XM_030985987	Pl_JAK_F2	CACCTTACCCATACTAGACAG	192
			Pl_JAK_R2	CTTGCCAGAGCCTCCGCTGAC
Lefty	Lefty	SUB6701449	Pl_Lefty_F2	CAGTCCAGACATGGGTGGCAG	182
			Pl_Lefty_R2	CATTTCGTCGACCACCTGCTG
maternal Vg1	M-Vg1	SUB6701449	Pl_M-Vg1_F1	GCACCTGCACCTAGAGACTC	145
			Pl_M-Vg1_R1	GCATGACCTTTTCCGGCCTG
NADH dehydrogenase	NADH	SUB2817153	Pl_NADH_F1	GTCTCCGTCGGATAAATCAAAG	194
			Pl_NADH_R1	CCGAAAAGGAAATAACGAAGC
Nemo-like kinase protein	NLK	AY442297	Pl_NLK_F1	CCTCTACCAGATTCTCAGAG	192
			Pl_NLK_R1	GTGACACAGTACTACCGCGC
Notchless protein	NOTCHLESS	PRJNA495004	Pl_Notchless_F1	GGGAAGCTAAGGCATCAGAC	145
			Pl_Notchless_R1	CGATCCTCTCAAGCACTTTAG
Patched	Ptc	SUB6701449	Pl_Ptc_F1	CGGTCGTCAGTATCATCATG	135
			Pl_Ptc_R1	GCAACCACGACTCCGTAAGC
Phospholipase C	PLC	AJ012336	Pl_PLC_F1	GAGACATTCACAGTGCCCAC	139
			Pl_PLC_R1	CTGACCGATACCAAGCTGTAC
PLAUF 3 RNA-binding protein AUF1 mRNA	PLAUF 3	AY682309.1	Pl_PLAUF3_F2	GGAGGATACGGCGGTGGCGG	182
			Pl_PLAUF3_R2	GTGTTGACTCCACAGGAGTG
Polyketide synthase	PKS	SUB6701449	Pl_PKS_F1	GCTTCCTCGACCAGTCTGTC	142
			Pl_PKS_R1	CCTCCGAAGACAGTCATCTG
signal transducer and activator of transcription	STAT1	PRJNA495004	Pl_STAT_F1	GTGTGTCAATCAGCCAGTGC	196
			Pl_STAT_R1	GTACATCATGAGCTTACCATTTC
Smoothened	Smo	SUB6701449	Pl_Smo_F1	CACGATCCATTACGGCGTTG	217
			Pl_Smo_R1	GCCCAACTCACACCCATGAC
sulfotransferase 1C2-like	SULT1	SUB6701449	Pl_SULT1_F2	CAGGCACTCACTGGCTCATG	140
			Pl_SULT1_R2	CTCTTCAGCTCTCGTCTTCG
Tumor necrosis factor alpha	TNF	SUB2817153	Pl_TNF_F1	CCTGATGTGTATGCCTCTATC	144
			Pl_TNF_R1	CAAGATCCTCATGTCAGGAAG

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.

4.4. 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.

4.5. 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.