Gene Expression Analysis of the Stress Response to Lithium, Nickel, and Zinc in Paracentrotus lividus Embryos

Many anthropogenic pollutants such as metals are discharged into the marine environment through modern sources. Among these, lithium (Li), nickel (Ni), and zinc (Zn) can interfere with biological processes in many organisms when their concentration rises. These metals are toxic to sea urchin embryos, affecting their development. Indeed, animal/vegetal and dorso/ventral embryonic axes are differently perturbed: Li is a vegetalizing agent, Ni can disrupt dorso-ventral axis, Zn can be animalizing. To address the molecular response adopted by embryos to cope with these metals or involved in the gene networks regulating embryogenesis, and to detect new biomarkers for evaluating hazards in polluted environments in a well-known in vivo model, we applied a high-throughput screening approach to sea urchin embryos. After fertilization, Paracentrotus lividus embryos were exposed to Li, Ni, and Zn for 24/48 h. At both endpoints, RNAs were analyzed by NanoString nCounter technology. By in silico analyses, we selected a panel of 127 transcripts encoding for regulatory and structural proteins, ranked in categories: Apoptosis, Defense, Immune, Nervous, Development, and Biomineralization. The data analysis highlighted the dysregulation of many genes in a metal-dependent manner. A functional annotation analysis was performed by the KEEG Orthology database. This study provides a platform for research on metals biomarkers in sea urchins.


Introduction
The use of metals has accompanied civilization and is closely linked to anthropogenic activities due to their chemical characteristics [1]. Metals are present in many environments at variable concentrations. Some of them are essential, for example, zinc (Zn), copper, and manganese, which play key roles in many cellular metabolic processes in living organisms, while other metals are non-essential, not having a physiological role, as for example lead, cadmium (Cd), and lithium (Li) [1,2]. Others, such as nickel (Ni), are essential elements in plants and microorganisms, while in animals their essentiality is controversial [3]. Metals can be released into the environment by different types of human activities such as mining, urban runoff, sewage discharge, and chemical compounds applied to crops for insect or disease control [2]. Regardless of their essentiality or not, when metal concentrations raise, they become chronic and toxic with harmful effects for both marine and terrestrial organisms, including humans [2]. In this regard, several metals have been included in the 2019 Substance Priority List by the Agency for Toxic Substances and Disease Registry (ATSDR) (https://www.atsdr.cdc.gov/spl/ (accessed on 30 June 2021)), because they pose the most significant potential threat to human health due to their known or suspected toxicity. The potential hazards for human exposure, and the bioavailability of some metals, such as lead, cadmium, zinc, copper, manganese, iron, mercury, arsenic, and barium, have been described in marine environments, as well as metals bioaccumulation by numerous marine organisms [4,5]. Although there are many studies on the toxic effects and accumulation of metals in different marine organisms [5], relatively little is known about The aim of this study was to obtain an overview of the genes modulated in P. Lividus sea urchin embryos to cope with acute doses of LiCl, NiCl 2, and ZnSO 4 after 24 and 48 h of treatments, by the multiplex analysis of NanoString nCounter technology. In total, we selected 127 genes organized into six categories, i.e., Apoptosis (Ap), Defense (Def ), Immune (Imm), Nervous (Ner), Development (Dev), and Biomineralization (Bm), and took into consideration the territorial expression and functions of the analyzed genes. In addition, we tried to identify the potential biomarker genes particularly sensitive to the three metals, i.e., Li, Ni, and Zn, respectively, and greatly affected by them.

Embryo Culture
Gametes were collected from gonads of the P. lividus sea urchin harvested along the northwestern coast of Sicily. Eggs were fertilized and embryos were reared at 18 • C in Millipore (Billerica, MA, USA) filtered seawater (MFSW) containing antibiotics (50 mg/L streptomycin sulfate and 30 mg/L penicillin).

Metals Treatment of P. lividus Embryos
Caution: LiCl, NiCl 2 and ZnSO 4 are hazardous chemicals and should be handled carefully.
To study the embryonic response induced by Li, Ni, and Zn in P. lividus sea urchin embryos, we took advantage of a previous study [26]. Briefly, embryos were continuously cultured in the presence of 30 mM LiCl (Li) (Sigma-Aldrich, St. Louis, MO, USA), 0.5 mM NiCl 2 (Ni) (Sigma-Aldrich, St. Louis, MO, USA) and 0.1 mM ZnSO 4 (Zn) (Sigma-Aldrich, St. Louis, MO, USA), henceforth defined as Li-embryos, Ni-embryos, and Zn-embryos, respectively. Metals were added to embryo cultures 30 min. after fertilization. To collect embryos for subsequent analyses, we chose the endpoints corresponding to 24 and 48 h postfertilization (hpf), when controls (unexposed embryos) reached gastrula and pluteus stages, respectively. Embryo morphologies were monitored using a microscope Axioscop 2 plus (Zeiss, Jena, Germany) both on live and fixed embryos with 0.1% formaldehyde in seawater. Images were recorded by an Axiocam camera (Model 412-312, Zeiss, Jena, Germany).
The ciliary band was specifically labeled with a monoclonal antibody 295 (kindly provided by Prof. McClay, Department of Biology, Duke University, 124 Science Drive, Box 90338, Durham, NC 27708, USA) diluted 1:5 in TBST. Rabbit polyclonal anti-Sox2 (Sigma-Aldrich, St. Louis, MO, USA) was used at a dilution of 1:25 in 0.25% DS-TBST. Both primary antibodies were incubated overnight at 4 • C with specimens and, after washing in TBST, the fluorescent-labeled secondary and appropriate antibodies (Alexa Fluor 488conjugated goat anti-mouse or Alexa Fluor 555-conjugated donkey anti-rabbit) (Invitrogen Molecular Probes, Carlsbad, CA, USA) were diluted 1:200 in TBST and incubated for 1 h at room temperature.
For double-labeling of 295 and Sox2, embryos were incubated in blocking buffer DS/BSA-TBST for 1 h and incubated with both primary antibodies (diluted 1:5 for 295 and 1:25 for Sox2) in 0.25% DS-TBST overnight at 4 • C. Secondary antibodies were used as described in the single immunostaining.
All embryos were observed under a Leica DM4000 microscope (Leica Microsystems, Wetzlar, Germany), equipped for epifluorescence, and images were recorded using a digital camera system. Images were edited using Adobe Photoshop CS2 software. Negative controls were performed for every set of experiments by omitting the primary antibodies.

nCounter NanoString Gene Expression Assay
Embryos from the samples treated with Li, Ni, and Zn, respectively, were processed following the manufacturer's instructions of the RNeasy mini-Kit (Qiagen, Germantown, MD, USA) to extract total RNA, then quantified using the D30 bio-photometer (Eppendorf, Hamburg, Germany). 100 ng of total RNA was used to apply nCounter NanoString technology (https://www.nanostring.com/company/about-us (accessed on 1 January 2018)) provided by Diatech Labline SRL (Jesi, Italy), in order to analyze the expression of a panel of 127 P. lividus transcripts, selected and retrieved from the National Centre for Biotechnology Information (NCBI) nucleotide database (http://NCBI.nlm.nih.gov (accessed on 1 January 2018)). Briefly, unique custom-made probes developed for the selected P. lividus sequences were hybridized with total RNA samples. Then, each hybridized barcoded transcript was detected and counted by a digital analyzer for image acquisition and data recording. The expression level was measured after normalization of the resulting digital counts of each transcript with those of the Pl-Z12-1 reference gene (Accession Number: LT900344) [28]. Then, the normalized counts of each transcript were compared between control and treated embryos to obtain the fold-change values.
We organized the analyzed genes into six arbitrary categories, in alphabetical order: Apoptosis (Ap), Biomineralization (Bm), Defense (Def ), Development (Dev), Immune (Imm), and Nervous (Ner). The heat map was generated by applying the conditional formatting in Microsoft Excel, setting as rule three range of fold values, <−2 (violet), >+2 (green), and between −2 and +2 (yellow) corresponding to the colors indicated in brackets.

Real-Time Quantitative PCR (qPCR)
Total RNA from P. lividus control and treated embryos, the same utilized for the nCounter Nanostring analysis, was reverse transcribed to obtain the corresponding cDNA as described in Russo et al. [26]. In order to quantify gene expression, the cDNAs were amplified by using SYBR Green technology, based on a Comparative Threshold Cycle Method [29], according to the manufacturer's instructions (Applied Biosystems StepOne-Plus instrument, Grand Island, NY, USA) and as previously described [26]. Pl-Z12-1 mRNA was used as a reference gene. The primer sequences of the genes utilized for qPCR were previously reported: 14-3-3 [30], metallothionein (mt) [31], jun [32], alx1, carbonic anhydrase (can), galectin-8 (gal-8) [25]. The primer specificity and accuracy were confirmed by the "melting curve", carried out during the real-time PCR.

KEGG Enrichment Analysis
Enrichment analysis was performed for all the genes modulated by each metal and those modulated by all the three metals, using the integrated KEGG Orthology (KO) database (in particular, KEGG pathway map and BRITE hierarchy databases; https:// www.genome.jp/kegg/ko.html (accessed on 1 June 2022)), selecting S. purpuratus as an organism (entry code: spu; genome number: T01019). The KO database includes molecular functions that are represented in terms of functional orthologs. Each gene is given a KO identifier, which represents a functional ortholog defined on the basis of the similarity of its sequence with orthology genes from other organisms, to generate organism-specific versions of KEGG pathways (database, path) and BRITE hierarchies (database, br), to deduce high-level functions of the organism.

Statistical Analysis
QPCR values are reported as the mean of three independent qPCR analyses ± SD, using cDNAs obtained by the same experiment of P. lividus embryos exposed to Li, Ni, and Zn used for the nCounter NanoString analysis. QPCR decimal values were transformed into negative values to compare them to NanoString values. Results were compared using OneWay analysis of variance, ANOVA followed by Tukey's HSD test, used as a post-hoc test for mean comparison. All the analyses were performed using the OriginPro 8.1 statistical program (OriginLab Corp), and the level of significance was set to p ≤ 0.05.

Developmental Effects of Li, Ni and Zn on Sea Urchin Embryos
It is known that metals such as Li, Ni, and Zn perturb territorial differentiation during the development of diverse sea urchin species. We have previously shown that the sublethal doses of 30 mM LiCl, 0.5 mM NiCl 2, and 0.1 mM ZnSO 4 produced homogeneous populations of abnormal P. lividus embryos, each with a characteristic morphology [26,33]. Here, we selected the same doses to continuously treat the embryos for 24 and 48 hpf for subsequent molecular analyses. Figure 1A shows the experimental design and Figure 1B-I confirms the previously observed morphologies and therefore the reproducibility of the experiment. In particular, at 24 hpf, the untreated normal embryos were at the gastrula stage with an elongated archenteron and two triradiate spicules (corresponding to the initial formation of the skeleton) ( Figure 1B). The vegetalized Li-embryo (reduced ectoderm and increased endomesoderm, see Ghiglione et al. [20]) had an early gastrula shape with a very short invaginated archenteron and many cells spread inside the blastocoelic cavity ( Figure 1C). The radialized Ni-embryo also had a short archenteron (indicated by a blue ellipse in Figure 1D) with cells grouped inside the blastocoelic cavity, some of which were above the equatorial plane (black arrow in Figure 1D). The animalized Zn-embryo looked like a blastula with a group of cells arranged as a presumptive invaginating archenteron (black arrow in Figure 1E). At 48 hpf, untreated embryos reached the pluteus stage ( Figure 1F, lateral view), with a tripartite gut, elongated and patterned spicules, and pigment cells (red arrow), see also the drawing in Figure 1J. At this stage, all treatments greatly affected embryos development, as shown by images ( Figure 1G-I) and drawings ( Figure 1K-M) reporting the different embryonic morphologies and the direction of the developmental axes, i.e., animal-vegetal (A/V) and dorso-ventral (D/V). Li-embryos had a typical endoderm-derived everted gut and a sphere with many pigment cells (red arrow) ( Figure 1G,K). Ni-embryos had a straight gut, a thickened animal cap that provided a typical bell shape, and some cells grouped near the archenteron (blue ellipse in Figure 1H,L). Zn-embryos maintained a blastula-like shape ( Figure 1I,M). To highlight the morphological diversities, we performed an indirect immunofluorescence (IF) analysis using the anti-Sox2 antibody ( Figure 1N-Q), a marker of the nervous system [34], and 295 monoclonal antibody (Ab-295), which recognizes the epithelium of ciliated cells of the embryonic arms, surrounding the mouth, called ciliary band associated to the embryonic nervous system ( Figure 1R-U) [35]. The anti-Sox2 is a commercial antibody against a synthetic peptide corresponding to residues 32-47 of human Sox2, homologous to the Strongylocentrotus purpuratus SoxB1. In normal embryos, anti-Sox2 labeled the ciliary band (see white arrow in Figure 1N) and some cells in the pluteus apex (see 1.4× magnification within the white frame). In metals-treated embryos, anti-Sox2 labeled a belt of cells in Li and Ni embryos, and very few cells in the Zn-embryos (see white arrows in Figure 1O-Q). In normal embryos, in addition to recognizing the ciliary band ( Figure 1R, white arrow), the Ab-295 also labeled two small central ganglionic structures, indicated by white ellipses in Figure 1R. In Li-embryos, the Ab-295 clearly labeled the animal ectoderm separated from vegetal endomesoderm ( Figure 1S). In Ni-embryos, the Ab-295 labeled the ciliary band delimited from the vegetal region ( Figure 1T, white arrow), while in Zn-embryos, it labeled the embryo uniformly ( Figure 1U). In the merged images ( Figure 1V-Y), the yellow color highlights the co-localization of anti-Sox2 and Ab-295 in some cells of the ciliary band in normal embryos ( Figure 1V, arrow) as well as in Li-and Ni-embryos (see arrows in Figure 1W and 1X, respectively).

Profiling Transcriptional Outputs in Li, Ni, Zn Treated Embryos
We applied a high-throughput molecular approach to analyze the gene expression profile of sea urchin embryos exposed to acute doses of Li, Ni, and Zn, with the aim of focusing on those genes differentially expressed to cope with these agents singly applied. 127 transcripts of the sea urchin P. lividus were specifically selected after a careful study from the NCBI database and used to perform the multiplex RNA analysis by nCounter NanoString technology in embryos at 24 and 48 hpf.
We chose to organize the selected genes into six categories described in Section 2. For some of them, we also reported some subcategories and their territorial expression, taking into consideration their functions (see Table S1).
In Figure 2, the heat map shows an overview of the genes modulated by the treatments with the three different metals, where fold-change values higher or lower than +2 and −2 (compared to controls) were considered as significant increase/decrease, as suggested by the NanoString technique. The fold-change expression values of the 127 genes, calculated on previously normalized data (see M and M), are reported in Table S1. The cell fill color identifies the responsive genes, i.e., those upregulated in green (>+2), downregulated ones in violet (<−2), and not responsive ones in light yellow. The light blue identifies the no-quantifiable genes (digital counts lower than 50), assuming that they were very poorly expressed in the control and/or in the samples and thus to be excluded according to the manufacturer's instructions. Figure 3 summarizes the percentages of genes modulated or not by the metals at 24 and 48 hpf using the color code previously established for the heat map. In general, most of the genes were not affected (values between −2 and +2, light yellow in Figure 3) by metal treatments at both endpoints analyzed or were not quantified (light blue in Figure 3). Notably, the highest percentages of unaffected genes (60.6% and 61.4% at 24 and 48 hpf, respectively) were observed in Ni-embryos. As for the affected genes, although the three metals modulated the analyzed genes differently, these modulations did not substantially vary at the two endpoints of 24 and 48 hpf for each metal. In particular, Li triggered the modulation of 38.6% of the genes (49 in total) at both 24 and 48 hpf, with higher values for the downregulated genes (28.3% and 27.6%) than for the upregulated ones (10.2% and 11%) at both endpoints. As already highlighted, Ni treatment modulated a low percentage of genes, i.e., 26% (33 out 127 genes) at 24 hpf and 22.8% (29 out 127 genes) at 48 hpf, with similar percentages between downregulated (14.2% and 9.4%) and upregulated (11.8% and 13.4%) genes at both endpoints. Even in Zn-embryos, similar percentages of modulated genes were observed at both endpoints, i.e., 34.6% (44 genes) at 24 hpf and 35% (45 genes) at 48 hpf, with noteworthy differences between up-(24.4% vs. 17.3%) and downregulated (10.2% vs. 18.1%) genes ( Figure 3).
We further summarized nCounter NanoString data in the scatter plots of Figures 4-6, to analyze in detail the various genes divided into categories and sub-categories.
In Figure 4, the data for the genes belonging to Ap, Def, Imm, and Ner categories at 24 hpf ( Figure 4A) and 48 hpf ( Figure 4B) are reported. Focusing on the Ap category, Ni did not affect genes at both endpoints analyzed, differently from Li and Zn which modulated the expression of these genes. In particular, apaf1 was downregulated by Zn at 24 hpf (−3.16-fold) and Li at 48 hpf (−2.94-fold), while bax was only upregulated by Zn at 48 hpf (2,68-fold). Li upregulated caspase-8 (3.20 and 2.70-fold) at both times, while Zn upregulated it at 48 hpf (2.51-fold). p63 gene was upregulated at both 24 and 48 hpf by Li (2.89 and 2.44-fold) and Zn (2.52 and 6.36-fold). Among Def genes at both analyzed endpoints, 14-3-3epsilon gene was downregulated by Li, Ni, and Zn with the greatest effect due to Zn treatment at 48 hpf (−9.16-fold).     For greater clarity in the presentation of the data, the genes of the Dev category have been divided into sub-categories, taking into account their spatial expression, where known, and are shown in Figures 5 and 6. Figure 5 shows genes of the Dev category with ectodermic expression, which were modulated by the metal exposures at 24 hpf ( Figure 5A) and 48 hpf ( Figure 5B In Figure 6, genes belonging to the Dev and Bm categories, having endomesoderm and mesoderm expression respectively, are reported, at 24 ( Figure 6A) and 48 hpf ( Figure 6B Among the Bm genes, we considered a group of TFs expressed in the skeletogenic cells, the primary mesenchyme cells (PMCs). At 24 hpf, the most affected TFs were coquillette and ske-t, which were down-and upregulated by the three metal treatments, respectively ( Figure 6A), while at 48 hpf coquillette was downregulated by Zn (−2.98) and ske-t was upregulated by Li and Ni (4.26 and 2.84-fold) ( Figure 6B).
Among the 127 genes analyzed by nCounter NanoString methodology, we selected six target genes, representing the categories that most interested us, namely Pl-alx1 (TF within Bm), Pl-can (skeleton within Bm), Pl-gal-8 (Imm), Pl-jun (TF within Bm), Pl-mt (Def ), Pl-sox9 (Dev), and performed comparative qPCR analyses to compare the results obtained with NanoString analysis (Figure 7). Comparison of qPCR data of gene transcription levels in Li-, Ni-, and Zn-embryos with NanoString data. Pl-alx1, Pl-can, Pl-gal8, Pl-jun, Pl-mt, and Pl-sox9 mRNA levels were analyzed compared to control embryos using the endogenous gene Pl-Z12-1 for normalization. The gray boxes indicate the threshold for insignificant changes following Li, Ni, and Zn treatments, i.e., fold values between −2 and + 2, at 24 (A) and 48 (B) hours post-fertilization (hpf). Symbols: NanoS, NanoString data, blue bars; qPCR data, red bars; Li, Li-embryos; Ni, Ni-embryos, Zn, Zn-embryos. The light blue triangles indicate not quantifiable values. Each bar of qPCR data represents the mean of three independent experiments of qPCR ± SD. QPCR mean values were significantly different according to the one-way ANOVA (p < 0.05), followed by the Tukey's test. The asterisks (*) indicate statistically not significant variations to the relative control.
QPCR is an alternative sensible technology to measure changes in mRNA expression levels and, also in this case, we considered significant values of fold change those higher or lower than ±2, compared to controls. In this range, the qPCR data are in agreement with those of nCounter NanoString analysis, except for a slight difference at 48 hpf in the trend of the genes Pl-alx1 and Pl-mt in Zn-embryos, Pl-gal-8 in Ni-embryos, although such variations are to be considered not significant as they remain within the threshold ±2 ( Figure 7B). As for the other samples, the trend of fold changes was in good agreement between the two analyses, although the magnitude of values was sometimes quite different, as the great difference observed for Pl-sox9 in Li-and Zn-embryos at 48 hpf. However, these differences are probably due to the diverse technologies used for the analysis, in agreement with what was reported by Prokopec et al. [36].

Metal-Dependent Sensitive Genes as Potential Biomarkers
Through our deep analysis, we might identify specific responsive genes for each metal, grouping those genes that are particularly sensitive to only one of the three metals analyzed, or sensitive to more than one of the metals, to be used as analytical tools for future environmental toxicity studies. The Venn diagram, and the relative table below the diagram (Figure 8), allowed to represent the amount of overlapping modulated genes and those not shared by the three metals, the latter being 18 and 17 for Li, 12 and 9 for Ni, 14 and 16 for Zn, at 24, and 48 hpf, respectively. In the Supplementary Table S2, we reported the lists of these genes, i.e., those specifically modulated by each of the three metals, those in common between Li and Ni (Li/Ni), Li and Zn (Li/Zn) or Ni and Zn (Ni/Zn) and those in common among all the three metals (Li/Ni/Zn), at 24 and 48 hpf. In Li-embryos at least one gene for each category was sensitive to it at 24 hpf, and many were highly modulated (8 out 18 decreased/increased more than 3-fold). Therefore, in addition to apaf1 (Ap), mt (Def ), and pax2/5/8 (Ner), it would be possible to choose among apkc, fzd1/2/7, id, paxB, smo, tetraspanin, wnt16 and wnt6 of the Dev category and alk1/2, bmp5/8, fgfr1, jun, p16 and sm50 of the Bm category, all genes that were not modulated by the other two metals. A similar approach could be followed to identify Zn and Ni responsive genes.
In particular, genes belonging to five out six categories were specifically modulated by Zn (5 out 14 and 4 out 16 decreased/increased more than 3-fold at 24 and 48 hpf, respectively), whereas Ni modulated mainly genes belonging to Ner, Dev, and Bm categories (6 out 12 and 6 out 9 decreased/increased more than 3-fold at 24 and 48 hpf, respectively).

Metal-Dependent Sensitive Genes and Functional Enrichment
Metal-related impacts on gene expression were seen in all the treatments with some overlapping in the enrichment of various KEGG Orthology (KO) pathways.

Discussion
In this study, we report an in-depth characterization of the response at the molecular level to Li, Ni, and Zn of the P. lividus sea urchin embryo, through a high-throughput screening approach using the nCounter NanoString technology. We have analyzed 127 genes selected from the NCBI database and arbitrarily organized them into six different functional categories based on an in-depth study of the data in the literature. As a whole, the modulation of the expression of the various genes represents a synthesis of the molecular defense systems adopted by P. lividus embryos to counteract Li, Ni, and Zn effects, respectively. The activation of the defense systems involves an energetic cost for the embryos, which results in developmental alterations and the formation of alternative morphotypes. Furthermore, we need to take into consideration the dual role of some of the defensive genes and proteins, i.e., in the regulation of developmental processes and in the protection against environmental risks, so that their recruitment into protective mechanisms might affect embryonic development.
We focused our study on a short period of time, i.e., 0−48 h of sea urchin embryo development, to evaluate the effects of Li, Ni, and Zn treatments at sublethal concentrations, not environmentally relevant, but capable of producing homogeneous populations of abnormal embryos, each with its own characteristic morphology. Indeed, various morphological malformations have been already described in different sea urchin species treated with Li, Ni, and Zn [20,21,37,38], including our recent study on P. lividus [26]. The most striking morphological effects of Li, Zn, and Ni treatments are on the axes: the animal/vegetal (A/V) one in vegetalized Li-embryos, which have less ectoderm and more endomesoderm than animalized Zn-embryos that, on the contrary, have more ectoderm and less endomesoderm, while Ni treatment affects the establishment of dorso/ventral (D/V), also called oral/aboral (O/A), axis. While the A/V axis is already present in the egg and it is rigidly fixed, the D/V axis is established after fertilization through complex molecular events, recently reviewed by Molina and Lepage [39], and is well known to be easily disturbed by many physical or chemical treatments and agents.
In addition to perturbing the axes patterning, the embryonic nervous system appeared to be affected mainly by Zn, but also somewhat by Li and Ni, on the basis of the embryonic localization of two neural markers, Sox2, and ciliary band. Firstly, a slight difference in the localization of Sox2 in P. lividus compared to S. purpuratus embryos has to be highlighted. In particular, the soxB2 gene, homolog to the human sox2, is expressed in the animal pole domain (which contains neurogenic ectoderm) of early S. purpuratus gastrula at 30 hpf and in the oral ectoderm around the ciliary band in pluteus at 72 hpf, while no signal was detected in the aboral ectoderm of pluteus [40]. Partially in agreement, the anti-Sox2 antibody labeled mainly the ciliary band in P. lividus pluteus at 48 hpf, but the protein was also present in some cells of the aboral ectoderm of the pluteus apex, which can therefore be neurogenic cells too. The sox2-expressing cells were lost in nearly all metal-treated embryos, in addition to the two small central ganglionic structures labeled by Ab-295. These results reveal that Li, Ni, and Zn might alter the embryonic nervous system, in agreement with data obtained from the Nanostring analysis as discussed in the following.

Looking for Potential Biomarkers
The morphological abnormalities of P. lividus embryos treated with metals are obviously associated with changes in gene expression and, consequently, in signaling pathways. Indeed, several genes appeared modulated by Li, Ni, and Zn at 24 (49,33, and 44 out of 127, respectively) and 48 hpf (49, 29, and 45 out of 127, respectively), of which 13 and 6 were the same for the three treatments (shown in the Venn diagram of Figure 8). In Table 1, we listed the genes affected by Li, Ni, and Zn that can be used as potential biomarkers for all of them. In the following, we will provide some information on each of these genes, when available, to highlight the molecular mechanism of the stress response of sea urchins exposed to Li, Ni, and Zn. Three genes were affected by all three metals at both endpoints, 24 and 48 hpf, i.e., 14-3-3epsilon, gata1/2/3, and nectin. The 14-3-3epsilon belongs to a highly conserved family of eukaryotic adaptor proteins found in diverse organisms and involved in many cellular processes. For example, 14-3-3 proteins interact with bax preventing its translocation into mitochondria and playing an anti-apoptotic role [41]. In sea urchin embryos, increased Pl-14-3-3epsilon mRNA levels were found following UVB irradiation [30], suggesting its implication in the regulative cascade activated in the stress response.
The gata1/2/3 gene codes for a DNA binding zinc finger TF that is involved in vertebrate hematopoiesis. Regarding its sea urchin homologs, Spgatac regulates immune cell specification in embryos and larvae [42], while Spgatae might be a useful biomarker for assessing the sea urchin hypoxic response, as suggested by its downregulation in the adult immune cells of S. nudus maintained in hypoxic condition [43].
Nectin is one of the proteins constituting the extracellular matrix (ECM) of the sea urchin embryo, known to be involved in the ecto-mesoderm signaling regulating skeleton development [44]. The changes in nectin gene expression are in good agreement with the defects in the development of the skeleton observed in all metals-treated embryos, so much so that it can be taken into account among the potential biomarkers.
Among the genes affected by the three metals only at 24 hpf, chordin and pax2/5/8 [45] belong to the Ner category, interesting for new studies of the nervous responses to the toxic action of pollutants. A recent review well described the embryonic neurogenesis in echinoderms, which proved to be a very complex process, and the molecular toolkit involved [35], which can therefore be considered in the search for new biomarkers of the stress response.
In the Dev category, six genes were modulated by the three metals at 24 hpf. Although little is known about the blastula protease (bp) 10 gene [46] and its ectoderm expression, we consider it rather interesting as bp10 expression was modulated also by Cd treatment [25], in addition to Li, Ni, and Zn.
The paxB gene belongs to the PAX family of transcription factors, which have been characterized in sea urchin embryos and are reported to be involved in visual system development in a variety of species [47].
Four genes constitute the sea urchin repertoire of smad genes that are a family of TFs activated downstream the TGF-beta signaling [48] and, among these, smad6 is known to be expressed by the presumptive dorsal ectoderm [49].
Univin is a growth factor of the TGF-beta superfamily involved in sea urchin skeleton morphogenesis [50]. It could be an interesting biomarker of different pollutants, for example, its downregulation was observed in P. lividus embryos exposed to gadolinium (20 µM) for 48 h [51].
Wnt5 acts with a short-range signal that activates a narrow band of ectodermal cells from the endoderm, integrating different types of positional information to correctly locate the production site of signals needed to guide skeletogenesis [52].
Certainly, the biomineralization genes would be interesting molecular markers, to put beside the morphological observation of the absence of the skeleton in embryos exposed to the three metals. In this respect, p16, coding for an acidic phosphoprotein of the skeleton [28], as well as ske-t and coquillette, two TFs belonging to the T-box family, should thus be taken into consideration. Ske-T is expressed early in the S. purpuratus embryos [53]. Coquillette, a member of the Tbx2 subfamily, is expressed late in the S. purpuratus embryo, restricted to the PMCs of gastrula embryos, involved in the skeleton formation [54].
Only six were the genes affected by all three metals at 48 hpf. Within the Def category, in addition to the 14-3-3epsilon gene previously described, members of the Hsp family, such as hsp70-II, were expected to be modulated. The upregulation of hsp genes is a typical stress response, well known in many organisms, which can have an anti-apoptotic effect [55], as also suggested for sea urchins [56].
The exogastrula-inducing peptides (egip) gene codes for small epidermal growth factorrelated peptides, which probably are secreted and then associated with the ECM, although their role is not yet known [57].
Vegfr is involved in the VEGF signaling that guides skeleton morphogenesis from the ectoderm [58] and it is one of the genes grouped in the Bm category. Most of these genes are involved in embryonic skeleton development and are mainly expressed by the PMCs, clearly recognizable in Li, Ni, and Zn embryos, even if not properly organized [26]. Together with other genes of the Bm category, this gene seems an interesting potential molecular biomarker, since its modulation could be associated with morphological observations on the skeleton development.
Beyond the potential biomarkers proposed in Table S2, some other genes could be considered potential biomarkers, both because they are greatly affected by one or more of these metals and because they can have interesting functions in the sea urchin. In particular, we could suggest can, hox7, foxA, pks1, irxa, sox9, and wnt16.
Can codes for a carbonic anhydrase metalloenzyme, which has a key role in the biomineralization process in many metazoans, also including sea urchins [59]. Interestingly, its expression is often affected by negative external stimuli, such as it is strongly downregulated by Cd and Cd in combination with UVB radiation in P. lividus embryos [25].
Within the Imm category, the pks1 gene, found in many organisms, is greatly affected by Ni and Zn in agreement with studies reporting its reduced expression by Ni in L. variegatus [60] and in P. lividus embryos [33]. In sea urchins, the PKS enzyme is involved in the synthesis of the red pigment echinochrome A and is expressed by pigment cells, a type of embryonic immune cell [61]. The modulation of pks1 can be associated with the morphological observation, checking the presence/absence of pigment in exposed embryos. For example, the absence of pigment cells was observed in P. lividus Ni-and Zn-embryos (this study; [26]) and in L. variegatus Ni-embryos [21]. Growing evidence supports the importance and role of metal ions in the immune system, both regulating the innate immunity and host defense functions, in what is now called "metal-controlled immunity" [62]. In this regard, the studies on the immune response induced by metals in a simple model organism such as the sea urchin might be of great support.
Sox9 belongs to the Dev category and, more in detail, to the gene subgroup having endo-mesoderm territorial localization. Interestingly, other agents, such as UVB and Cd, caused its downregulation in P. lividus embryos [25]. Sox9 gene, also known as Sp-soxE in the S. purpuratus genome (www.echinobase.org (accessed on 1 January 2019)) [63], belongs to a family of minor groove DNA binders, widely expressed during sea urchin embryo development, and is one of the regulatory genes expressed during embryonic myogenesis [64].
Irxa and hox7 are two dorsal genes [49] grouped in the Dev category, together with many genes involved in axes specification and patterning, including A/V and D/V, which were differently altered by the three metals, as described above.
The foxa gene, a forkhead transcription factor, is an integral component of the endoderm specification subcircuit of the endomesoderm gene regulatory network, as characterized in the S. purpuratus embryo [65]. The Authors showed that blocking foxa expression, by the injection of morpholino antisense oligonucleotides at high concentrations, caused the failure of gut formation, a process that involves the regulation of structural genes and cell adhesion molecules [65]. The foxa gene, and the morphology of the embryonic gut, could be interesting to consider as pollutant biomarkers.
Most of the Wnt and their receptor frizzled (fzd) genes have expression domains associated with archenteron development, thus it was not surprising that wnt16 was greatly upregulated by Li. Specifically, at the late gastrula stages, different populations of endodermal cells expressed wnt16. At later stages, prism and early pluteus, the expression of wnt16, as many other Wnt genes, becomes even more restricted in a full ring around the anus [66].
Altogether these results provide a general framework for the possible responses at the molecular level of sea urchin embryos to different metals, from which to extrapolate useful information on possible new markers to be taken into account in ecotoxicological studies. Indeed, the reorganization of the embryonic territories induced by the different metals, clearly distinguishable at the morphological level, i.e., animalization (Li), vegetalization (Zn), and radialization (Ni), could now be associated with some genes modulated exclusively by each of the three metals. Nevertheless, among the various aspects that must be considered in the evaluation of metals toxicity, i.e., their concentration, the different sensitivities of the embryonic stages, and of the various species of sea urchins, it would be worthwhile to identify a panel of modulated genes for each metal.

Analysis of Pathways Affected by Metals
The data obtained by the analysis of pathways enrichment can provide further information on the molecular responses of sea urchin embryos exposed to Li, Ni, and Zn, useful for future exploitations in ecotoxicological studies.
Only 39 of the 82 genes modulated by each metal and those modulated by all the three metals were found in 13 pathways by KEGG Orthology (KO, https://www.genome.jp/ kegg/ko.html (accessed on 1 June 2022)) analysis (Table S3). The pathway of TFs was the most enriched, followed by the signaling pathways of Wnt, TGF-beta, FoxO, and notch, as well as the pathways of protein processing in the endoplasmic reticulum, endocytosis, membrane trafficking, cytoskeleton proteins, G protein-coupled receptors, DNA repair and recombinant proteins, peptidases and inhibitors, and ubiquitin-mediated proteolysis. All the three metals impacted the TFs pathway, modulating 11 genes involved in numerous and different cellular processes.
It seems evident that the main enriched pathways found by the KO analysis were the signaling ones, which control all aspects of embryo development. In particular, all the metal treatments here analyzed impacted the Wnt signaling pathway, up or downregulating some of the related genes (fzd1/2/7, fzd4, wnt16, wnt6, wnt4, wnt5, tcf/lef ). The Wnt signaling pathway is highly evolutionarily conserved in organisms, as it is involved in body axis patterning, cell proliferation and migration, and cell fate specification. In sea urchins, this signaling pathway has been extensively studied and is known to be a crucial regulator of embryo development [67], including the A/V axis formation. Furthermore, it is known to be involved in neurogenesis, restricting neuroectoderm to the anterior pole of the embryo with precise spatio-temporal control, involving also JNK signaling, in addition to wnt1 and wnt8 [68]. The main receptors for Wnts are Frizzled receptors that are involved in the formation and function of the neuronal circuits in mammalian central nervous system [69].
The third most enriched pathway for all metals was the TGF-beta signaling (bmp2/4, nodal, pitx2, chordin, smad6/7), involved in a wide range of cellular processes and highly conserved all over the animal world. In addition to ligands and receptors, the TGF-beta signaling pathway requires the involvement of receptor-regulated proteins (the intracellular mediators Smads) and various TFs to regulate the expression of target genes. Among the genes belonging to this KO pathway, chordin was included in our Ner category, as it is one of the marker genes for apical and neural territories, in addition to gfi1 [49], onecut/hnf6 [70], pax6, sox9, and pax2/5/8.
Most interesting for our study, are the three genes (foxG, foxO, p38MAPK) that enriched the FoxO signaling pathway. The Forkhead (Fox) family includes TFs regulating target genes involved in diverse biological processes, proliferation, differentiation, longevity, metabolism, and apoptosis. In the sea urchin, at least 22 genes belonging to the FOX family have been annotated [71], although many of them have not yet been characterized. The two fox genes analyzed are differently expressed in late sea urchin embryos, i.e., foxO is expressed in PMCs, suggesting its involvement in biomineralization, whereas foxG expression is restricted to the ciliated band, suggesting an involvement in nervous processes [71]. By our KO analysis, p38MAPK gene enriched the FoxO signaling pathway, even though this gene is known to have a dual role in sea urchin, i.e., an essential regulator of embryonic development, involved in D/V axis specification and skeletogenesis [72], and a stress marker [73].
Interestingly, the functional categories of protein processing in the endoplasmic reticulum and DNA repair and recombinant proteins were impacted by the three metals, i.e., the genes bax, hsp70-II, and 14-3-3epsilon that indeed are stress marker genes.
We believe that this study not only can greatly facilitate the identification of metalspecific biomarkers but can also be useful for the identification of functional developmental pathways altered by pollutants in P. lividus embryos. This is obviously a preliminary study, considering that the metals were administered individually, while in nature the effects are due to the combined action of several toxic agents.

Conclusions
As many industrial activities release metals that can be harmful to human health as well as to the environment, it is extremely important to develop new tools that can detect these pollutants quickly and easily, but at the same time with high levels of sensitivity. In this direction, our study contributes both to expanding the knowledge of the molecular biomarkers that could be used in the environmental monitoring, and to facilitating the identification of panels of modulated genes to be used in the evaluation of metals toxicity.
The sea urchin embryo is a well-known model used in ecotoxicological studies to detect a wide range of chemical, physical and natural agents, such as metals, radiations, ocean acidifications, and sediments as well as to study the effects of combined agents. In addition, a huge amount of data is now available concerning the defensive, nervous, and immune systems and the related GRN sub-circuits governing them in sea urchin embryos, which provides new attractive objects of study, such as the immune and nervous responses to the toxic action of pollutants.
All of this can be of great help to provide new approaches to address ecotoxicological studies, encouraging new possible monitoring and intervention strategies.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/toxics10060325/s1, Table S1: Fold-change values of the analyzed genes.  Institutional Review Board Statement: All applicable international, national, and/or institutional guidelines for the care and use of invertebrate animals were followed.
Data Availability Statement: Not applicable.