The mean survival rate between 12 and 16 dph was 29 ± 14%, which is consistent with the survival rates recorded in aquaculture at these early life stages (Villamizar et al., 2011). The mean growth percentages were +18% for solvent control (0.70 ± 0.03 cm) and E2 0.4 ng.L^-1 (0.70 ± 0.06 cm), +14% for E2 40 ng.L^-1 (0.67 ± 0.04 cm) and + 20% for BPA 1.6 μg.L^-1 (0.71 ± 0.04 cm) and BPA 160 μg.L^-1 (0.71 ± 0.04 cm) between 12 and 16 dph, without any significant difference in larval size between the different treatments at the end of the exposure. Abnormal swim bladder was observed in 33% of 12-dph individuals at the beginning of the exposure, while it was 37 ± 10% in 16-dph individuals with about ⅔ of over-inflation and ⅓ of under-developed swim bladders. Larvae with abnormal bladder were discarded for subsequent mRNA or Alizarin analyses.

Fixed larvae were gradually rehydrated in Ethanol/KOH 0.5% (v/v) then in KOH 0.5%. They were depigmented by incubation in a 0.3% hydrogen peroxide in KOH 0.5% solution for 30 minutes in the dark, washed in KOH 0.5%, incubated overnight in a solution of 0.01% Alizarin red in KOH 0.5%, and washed in KOH 0.5%. Stained larvae were transferred in graded series of KOH 0.5%/glycerol and then stored and imaged in 100% glycerol. Staining was carried out simultaneously for all treatments to avoid technical variability of the staining.

For each pictured individual after Alizarin staining (ventral and lateral views), head mineralized structures were identified and counted. Based on the mineralization level of each cranial unit studied (0: no mineralization, 0.25: 1 quarter mineralized, 0.5: half mineralized, 0.75: 3 quarter mineralized, 1: fully mineralized), an overall mineralization score was assigned to each larva corresponding to the sum of the cranial unit mineralization scores (see Supplementary data Figure S2 as an example of mineralization level quantification).

Statistical analyses were performed on GraphPad Prism (version 6, GraphPad Software Incorporated, La Jolla, CA 268, USA). Since data did not fit normality or homoscedasticity assumptions, differences between groups were tested using the nonparametric Kruskal-Wallis test followed by Dunn’s multiple comparison test. Multiplicity adjusted p-value are given, accounting for familywise error rate due to multiple pairwise comparisons. The strength and direction of the linear relationship between fish standard length and head mineralization level was tested using the nonparametric Spearman correlation test (p<0.05). To compare the effect of treatment on mineralization levels, nonparametric ANCOVA was performed to test differences between predicted non-linear regression curves using the r package fANCOVA (0.6-1). Three pairwise comparisons were made for each data set. Therefore, the significance threshold was adjusted to p <0.017 according to the Bonferroni correction.

Only the otoliths, cleithrum and jaws are mineralized at 12 dph in control conditions, similarly to the previous description of (Darias et al., 2010) ( Figure 1A ). In addition to these structures, the parasphenoid, branchial arches and basioccipital bones are mineralized at 16 dph ( Figure 1C ). At these stages, the vertebrae are not visible as they are not yet mineralized under our rearing conditions ( Figure 1 ).

Three larval developmental time points were chosen to compare the expression levels of selected genes involved in early skeletal initiation in a first cohort: (i) 6 dph, when only the cleithrum and maxillary have started their mineralization; (ii) 10 dph, when the mandibular has started its mineralization; (iii) 17 dph, when the mineralization of many other cephalic structures is undergoing ( Table 2 ).

While col1a2 expression appears constant over time, the chondrocyte-associated col2a1a showed a decrease in expression level between 6 and 10 dph which was maintained at 17dph, while the perichondrium/osteoblast-associated col2a1b increased expression between 6 and 17 dph. These patterns were consistent with the timing of mineralization initiation in the head of the larvae. The three genes sparc, bgp1a and alp exhibited an increase in their expression levels detected at 17 dph, with a start as early as 10 dph for bgp1a ( Table 2 ). Regarding ihha and ihhb, sox9a and sox9b, runx2 and sp7 genes, no significant transcriptional regulation of these signaling and transcription factors was detected according to developmental timing ( Table 2 ). Though a 1.9-fold transient decrease in ctsk mRNA expression was observed at 10 dph, when compared with its expression at 6 dph, the expression of rankl and acp5 genes did not significantly change in the earliest phase of skeletal mineralization ( Table 2 ).

Regarding estrogenic signaling related genes, we detected an increased expression of erb1 between 10 and 17 dph, but no variation of erb2 in this time window ( Table 2 ). Concerning the gper family, only gperb was detected and its mRNA expression increased between 6 and 10 dph. Among the five estrogen-related receptors, errb2 and errga expression is constant over time, erra and errb1 expression decreases at 10 dph and 17 dph, respectively.

Comparisons between gene expression levels show that extracellular matrix genes have the highest levels of expression, i.e. the three structural collagen genes (col1a2, col2a1a, col2a1b) and sparc are 10 to 100 times more represented than the others ( Figure 2 ), which highlights the intense secretory activity of skeletal cells (chondrocytes and osteoblasts) over this period. In addition, the estrogen-related receptors erra, errb1 and errga are about one order of magnitude more expressed than the nuclear estrogen receptors erb1 and erb2 ( Figure 2 ) suggesting their potential functioning in ER signaling modulation at these developmental stages, despite not identifying in which cell types they are expressed.

Variation of head mineralization levels in 16-dph larvae is represented in Figure 3 . Head mineralization score is positively correlated to the standard length for every treatment with Spearman correlation coefficient values ranging between 0.68 and 0.88 (p<0.001). Larvae exposed to E2 40 ng.L^-1 for 4 days exhibit significantly higher head mineralization levels compared to the solvent control, along the full range of observed sizes (from 0.55 to 0.75 cm; nonparametric ANCOVA, p<0.01, Figure 3A ). Conversely, no significant difference in head mineralization level was statistically detected between larvae exposed to E2 0.4 ng.L^-1, BPA 1.6 and 160 μg.L^-1 compared to the solvent control.

No significant effect of E2 0.4 ng.L^-1 treatment could be detected on the expression of any of the selected skeletogenesis-related genes. Conversely, the E2 40 ng.L^-1 induced a significant increase in expression of 9 out of 16 genes ( Figure 4 ), four of them are expressed by osteoblasts (col1a2, col2a1b, sparc and runx2), one (ihha) is expressed by hypertrophic chondrocytes to activate osteoblast differentiation, both sox9a and sox9b are expressed during chondrogenesis and two genes (rankl, ctsk) are expected to be expressed by osteoclasts. The expression of skeletal genes is significantly modified only with high concentrations of E2 and seems to act mostly on bone secreting (osteoblasts) and degrading (osteoclasts) cells.

In contrast, we observed a general tendency towards a dose-dependent up-regulation of most tested genes in both BPA exposure experiments ( Figure 4 ). Exposure to BPA 1.6 μg.L^-1 induced a significant over-expression of 6 of the 16 studied genes: three extracellular matrix markers (col1a2, bgp1b and sparc), the transcription factor sox9b and the osteoclast markers rankl and ctsk. This effect is quite similar to the observed results of E2 40 ng.L^-1 treatment. In contrast, the exposure of BPA 160 μg.L^-1 has a significant effect for 11 of the 16 genes studied, and often with a stronger overexpression than the values obtained after E2 40 ng.L^-1 exposure ( Figure 4 ; Table 3 ). The BPA treatment with high concentration therefore has a stronger impact on gene expression for all skeletal cells, both quantitatively and qualitatively.

No significant effect of the E2 0.4 ng.L^-1 treatment could be detected on the expression of any of the ER or ERR encoding genes ( Figure 5 ). Conversely, E2 40 ng.L^-1 induced a significant overexpression of gperb, erra, errb2 and errga. This exposure condition also induced a weak (about 2-fold) but non-significant up-regulation of erb1 and erb2 ( Table 3 ). Three estrogen-related receptors were induced in response to BPA: errb2 and errga at BPA 1.6 μg.L^-1, and erra at BPA 160 μg.L^-1. In contrast to the results obtained with E2 treatments, gperb mRNA expression was not significantly modified by BPA exposure but the expression of erb2 was significantly overexpressed compared to solvent control in response to both tested concentrations of BPA. As a consequence, in both cases, the effects of exposure might be amplified by a putative better availability of estrogen receptors (Gperb with E2, or ERβ with BPA), and of estrogen-related receptors.

The estrogen-ER genomic signaling pathway involves the binding of ERs to ERE in the promoter of the target genes. ERRs can also bind the ERE and they may act as modulators of this genomic pathway (Acconcia et al., 2015; Tanida, 2022). We tested the presence of putative ERE consensus sequences in the non-coding sequences of all 23 genes, with 1/Kd>0.08 threshold: nine genes displayed one putative ERE (col1a2, col2a1a, col2a1b, sox9a and sox9b, rankl, erb2, gperb, errb2) while acp5, erra and errga had two putative EREs and sparc and ctsk had three in their vicinity ( Table 3 ). In these conditions, no ERE was identified in the 15kb surrounding the promoter region for bgp1a and bgp1b, ihha and ihhb, runx2, sp7, alp, erb1 or errb1. Putative EREs were therefore identified in the vicinity of 14 of the tested genes, of which 13 showed modified expression after E2 and/or BPA exposure. Only acp5 showed no significant change in its expression after E2 or BPA exposure despite the detection of two putative EREs. The expression of col2a1a and erb2 was not altered in the E2 treatments but increased in the BPA treatments. Of the 9 genes without any detected ERE, 5 genes (ihhb, sp7, alp, erb1 and errb1) display no response to E2 and BPA treatments, supporting independence of their transcriptional regulation from a genomic estrogen signaling. Despite the absence of putative ERE, two genes (bgp1a and bgp1b) have their expression induced by the BPA treatments only, the ihha gene had its expression induced by E2 40 ng.L^-1, and runx2 had its expression increased at E2 40 ng.L^-1 and BPA 160 μg.L^-1.

Since early developmental stages are often more sensitive to xenobiotic exposure, studying the potential precocious deleterious effects of xenobiotics is of great interest in ecotoxicology and aquaculture. The potential physiological disturbance due to estrogenic compounds exposure occurring at early life stages, when blood estrogen concentrations are low or undetectable, have to be addressed, because low concentration of estrogens may potentially lead to significant endocrine disruption or other developmental failures (Chin et al., 2018). In non-model species of ecological interest, early life stages are generally less considered due to rearing or other technical difficulties, yet in the European seabass Dicentrarchus labrax, as in other fish species, critical developmental steps for survival occur (Chatain and Dewavrin, 1989): mouth opening and exotrophy start around 4-5 dph (Barnabé, 1974; Cucchi et al., 2012) and the correct establishment of the swim bladder occurs between 3 and 7 dph (Chatain, 1986). In our study, the rates of swim bladder anomalies observed were in the range of those commonly observed in aquaculture in Dicentrarchus labrax (Chatain and Dewavrin, 1989; Miller, 2009). Its early skeletal development is cephalic and is linked with the function of collecting food and the development of the respiratory function. We detected an acceleration of head mineralization as a result of specifically the E2 40 ng.L^-1 treatment (not detected in any BPA treatment, or lower E2 treatment), showing that exposure to waterborne estrogens has an impact on fish skeletal development even at early larval stages.

In the present study, despite the lack of phenotypical effects in most of our 4-day treatments, results obtained after exposure to E2 (40 ng.L^-1) and BPA (1.6 and 160 µg.L^-1) suggest that both molecules exert a stimulatory effect on skeleton-related genes, in a dose dependent manner. The identification of putative ERE in the promoters of most of the up-regulated skeleton genes supports the hypothesis of genomic regulatory pathways as a result of E2 and BPA exposure. The concentration of E2 40 ng.L^-1 induced a significant increase in expression of genes known to be expressed by osteoblasts and hypertrophic chondrocytes, involved in osteoblast differentiation. This transcriptional induction of osteoblast-linked gene expression suggests a stimulation of osteoblast activity, but also of osteoblast differentiation (runx2 and ihha are early cell differentiation markers), having as a physiological consequence more head mineralization ( Figures 3A , 4 , 6 ). In addition, the concentration of E2 40 ng.L^-1 also induced a significant increase in expression of genes known to be expressed by osteoclasts involved in bone remodeling. A lower concentration of E2 did not induce such changes in gene expression, suggesting that estrogenic signal has to be strong to have an effect at early stages where the level of estrogen receptors, and therefore estrogen sensitivity, is low. Interestingly, ihha and runx2 do not have putative ERE while both their expression was increased upon treatment with E2 40 ng.L^-1. This up-regulation could be indirect and mediated by other transcription factors, or could be due to non-genomic regulation pathways (Björnström and Sjöberg, 2005) including, at least in part, Gperb (Pinto et al., 2014b).

Compared to E2, BPA treatment exerted an overlapping but wider positive effect on the transcriptional expression of skeletal genes. Notably, the matrix proteins encoding genes col2a1a, bgp1a and bgp1b were upregulated only by BPA and not by E2 ( Figures 4 , 6 ). These genes are considered to be chondrocytes and osteoblasts derived (Dale and Topczewski, 2011; Eames et al., 2012; reviewed in Leurs et al., 2021) though in situ hybridization techniques would help to localize mRNA expression to specific cell types in the seabass. However, the significant stimulatory effect of BPA was not correlated to any significant change in head mineralization ( Figure 3B ). The lack of significant induction of major precocious signaling and transcription factors (namely sox9a, ihha, ihhb, runx2) may explain that the head of larvae exposed to BPA 1.6 μg.L^-1 was not more mineralized than in controls. In contrast, BPA 160 μg.L^-1 exposure resulted in the same or stronger transcriptional activation of mineralization-related genes, as compared to E2 40 ng.L^-1 ( Table 3 ). However, a stronger activation of osteoclast-related genes was detected in BPA 160 μg.L^-1 exposed larvae, as shown by higher levels of transcription of rankl (10.3-fold changes versus 4.4-fold in E2 40 ng.L^-1) and ctsk (6.7-fold changes versus 4.0-fold in E2 40 ng.L^-1) ( Table 3 ). This suggests that the balance between bone mineralization (due to chondrocytes and osteoblast markers activation) and bone degradation (due to osteoclast markers activation) may be differentially regulated by E2 and BPA at elevated concentrations ( Figure 6 ). To confirm this hypothesis, osteoclast activity should be validated by TRAP/ACP5 staining in larval tissues. Overall, our results suggest that E2 exerts a positive physiological effect on head mineralization and that BPA does not strictly mimic this positive effect at skeletal level, though BPA can transcriptionally induce the same genes as E2 ( Figure 6 ).

The genomic effect of both E2 and BPA may act through the ERβ1 endogenously up-regulated in larvae during the studied developmental stages ( Table 2 ). In both BPA treatments, this genomic effect may be also amplified through ERβ2 whose expression was induced ( Figure 5 ). However, E2 40 ng.L^-1 exposure may activate or amplify a non-genomic activity as the expression level of gperb is increased ( Figure 5 ). In mammals, GPER is expressed in osteocytes, osteoclasts, osteoblasts and chondrocytes and is involved in estrogen-induced bone growth and development (Urushitani et al., 2002; Heino et al., 2008). GPERs are known to be involved in scale mineralization in seabass (Pinto et al., 2016). It may be hypothesized that the difference between E2 and BPA effects is an over-activation of the non-genomic pathway by E2, as opposed to an over-activation of the genomic pathway by BPA. This may explain the different outcome in mineralization. Despite this difference, both molecules similarly impact erra, errb2 and errga expression levels. During bone development, ERRs were previously shown to be involved in the differentiation and functional activation of osteoblasts and osteoclasts (Bonnelye and Aubin, 2013). The fact that EREs were found in the erb2, gperb, erra, errb2 and errga promoters supports the hypothesis for a potential crosstalk between nuclear ERs and orphan ERRs in the European seabass. Indeed, the ERRs have been postulated to modulate E2 signaling by either synergizing or competing with ERs in regulating multiple shared transcriptional targets (Horard and Vanacker, 2003; Cheung et al., 2013). In contrast to ERs, ERRs assume an active state without a ligand being bound to the ligand binding domain (LBD) (Kallen et al., 2004) which makes them constitutive transcription factors when linked to their ERR-response Element (ERRE). However, BPA was suggested to modify ERRs genomic activity through its binding, making BPA also a modulator of ERR activity (Tohmé et al., 2014). In addition, the presence of ERREs was identified in human bgp (Wang and Wang, 2013). As a consequence, some observed variations between E2 and BPA exposure effects on the transcriptional gene activation, such as those of bgp1a and bgp1b, might be coming from this effect of BPA on ERRs. In mice, in silico analysis found putative ERRE sites in the promoter regions of a number of skeletal development genes, including ihh and runx2 (Cardelli et al., 2013). This could explain observed upregulations of genes without an identified ERE (ihha and runx2) in our E2 and/or BPA treatments. Our data highlight the importance of including ERRs in future research on the effects of (xeno)estrogens on physiology, given their putative interaction with estrogen signaling pathways, especially in teleost where their role is probably very underestimated.

The action of E2 and BPA also do not seem to be strictly similar on ERRs: in our experiment, the expression of erra was activated by high concentrations of both E2 and BPA, whereas errb2 and errga were activated only following E2 40 ng.L^-1 and BPA 1.6 μg.L^-1 exposure ( Figure 5 ). In mouse, ERRα is highly expressed throughout the osteoblast developmental sequence and plays a physiological role in bone formation in vitro (Bonnelye et al., 2001) and has been shown to regulate the expression of sox9a/b and the downstream genes necessary for cartilage development during zebrafish embryogenesis (Kim et al., 2015). As a consequence, our variation of sox9a and sox9b expression might also result from disruption of the endogenous ERR effects. In mice, ERRγ negatively regulates chondrocyte proliferation and positively regulates matrix synthesis to coordinate growth plate organization (Cardelli et al., 2013), and its activity is thought to be dependent on Runx2 (Cardelli and Aubin, 2014). ERRγ was also previously shown to negatively regulate osteoblast differentiation and bone formation (Jeong et al., 2009). These variations effects of ERRγ on skeletal cells may also explain why E2 treatment result in more calcium deposition, while BPA treatment does not lead to this morphological effect.

This study brings new insights into the regulatory mechanisms of skeletogenesis by E2 and into the effects of waterborne exposure to BPA on early skeletal development of a marine teleost fish. More detailed insights into their effects on cell activity would require functional and cell-centered studies, in particular in situ hybridization of genes associated with mineralization for their specific cell type localization, and activity assays such as TRAP/ACP5 or ALP. Regarding the regulation of estrogenic signals, complex ligand-dependent and ligand-independent pathways could be at stake, involving the different subfamilies of estrogen receptors, with potential cross-talks. Our results suggest that the role of ERR and the interaction between estrogen-linked signaling pathways are probably very underestimated especially in teleost. Last, beside the membrane-bound and estrogen(-related) receptors, emerging data from in vitro and in silico models show that BPA binds with a significant number of hormone receptors, including androgen receptors, as well as the thyroid hormone receptor, glucocorticoid receptor, and PPARγ (Delfosse et al., 2012; MacKay and Abizaid, 2018; Cimmino et al., 2020), all of them playing a role in the development and maintenance of cartilage and bone (Zuo and Wan, 2017; Gouveia et al., 2018; Chou et al., 2021). This wider range of signalization for BPA might be another aspect explaining the differences observed between E2 and BPA treated larvae.