O R Cooper

D D Parrish

A Stohl

M Trainer

P Nédélec

V Thouret

Jean-Pierre Cammas

Increasing springtime ozone mixing ratios in the free troposphere over western North America

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Increasing springtime ozone mixing ratios in the free troposphere over western North America In the lowermost layer of the atmosphere-the troposphereozone is an important source of the hydroxyl radical, an oxidant that breaks down most pollutants and some greenhouse gases [START_REF]Ground-level Ozone in the 21st century: Future Trends, Impacts and Policy Implications Royal Society policy document 15/08[END_REF] . High concentrations of tropospheric ozone are toxic, however, and have a detrimental effect on human health and ecosystem productivity [START_REF]Ground-level Ozone in the 21st century: Future Trends, Impacts and Policy Implications Royal Society policy document 15/08[END_REF] . Moreover, tropospheric ozone itself acts as an effective greenhouse gas [START_REF]Climate Change 2007-The Physical Science Basis, Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF] . Much of the present tropospheric ozone burden is a consequence of anthropogenic emissions of ozone precursors [START_REF] Horowitz | Past, present, and future concentrations of tropospheric ozone and aerosols: methodology, ozone evaluation, and sensitivity to aerosol wet removal[END_REF] resulting in widespread increases in ozone concentrations since the late 1800s [START_REF] Horowitz | Past, present, and future concentrations of tropospheric ozone and aerosols: methodology, ozone evaluation, and sensitivity to aerosol wet removal[END_REF][START_REF] Staehelin | Trends in surface ozone concentrations at Arosa (Switzerland)[END_REF][START_REF] Oltmans | Long-term changes in tropospheric ozone[END_REF][START_REF] Thouret | Tropopause referenced ozone climatology and inter-annual variability (1994-2003) from the MOZAIC programme[END_REF][START_REF] Zbinden | Mid-latitude tropospheric ozone columns from the MOZAIC program: climatology and interannual variability[END_REF] . At present, east Asia has the fastestgrowing ozone precursor emissions [START_REF] Van Der | Trends, seasonal variability and dominant NO x source derived from a ten year record of NO 2 measured from space[END_REF] . Much of the springtime east Asian pollution is exported eastwards towards western North America [START_REF] Stohl | On the pathways and timescales of intercontinental air pollution transport[END_REF] . Despite evidence that the exported Asian pollution produces ozone [START_REF] Hudman | Ozone production in transpacific Asian pollution plumes and implications for ozone air quality in California[END_REF] , no previous study has found a significant increase in free tropospheric ozone concentrations above the western USA since measurements began in the late 1970s [START_REF] Oltmans | Long-term changes in tropospheric ozone[END_REF][START_REF] Jaffe | Increase in surface ozone at rural sites in the western US[END_REF][START_REF] Schnadt Poberaj | Long-term changes in UT/LS ozone between the late 1970s and the 1990s deduced from the GASP and MOZAIC aircraft programs and from ozonesondes[END_REF] . Here we compile springtime ozone measurements from many different platforms across western North America. We show a strong increase in springtime ozone mixing ratios during 1995-2008 and we have some additional evidence that a similar rate of increase in ozone mixing ratio has occurred since 1984. We find that the rate of increase in ozone mixing ratio is greatest when measurements are more heavily influenced by direct transport from Asia. Our result agrees with previous modelling studies, which indicate that global ozone concentrations should be increasing during the early part of the twenty-first century as a result of increasing precursor emissions, especially at northern mid-latitudes [START_REF] Stevenson | Multimodel ensemble simulations of present-day and nearfuture tropospheric ozone[END_REF] , with western North America being particularly sensitive to rising Asian emissions [START_REF] Jacob | Effect of rising Asian emissions on surface ozone in the United States[END_REF] . We suggest that the observed increase in springtime background ozone mixing ratio may hinder the USA's compliance with its ozone air quality standard.

Inventoried global anthropogenic ozone precursor emissions of nitrogen oxides (NO 1 NO 2 5 NO x ) and volatile organic compounds increased rapidly during 1960-1990, but levelled off in the 1990s [START_REF]RETRO Emission Data Sets and Methodologies for Estimating Emissions, Work Package 1[END_REF] . Comprehensive inventories of ozone precursors have not yet been compiled for the years since 2000 and uncertainties remain whether total global emissions have increased since then. However, satellite measurements of column NO 2 indicate increased NO x emissions during 1996-2005 in China (up to 29% per year) and other locations in Asia, while emissions decreased in Europe and the USA [START_REF] Van Der | Trends, seasonal variability and dominant NO x source derived from a ten year record of NO 2 measured from space[END_REF] . The most recent bottom-up inventories indicate that south and east Asia's NO x emissions increased 44% during 2001-2006, with an increase of 55% within China [START_REF] Zhang | Asian emissions in 2006 for the NASA INTEX-B mission[END_REF] , while ozone precursor emissions decreased by more than a third across Europe 1 (1990-2005) and the USA 17 (1985-2008).

It is possible that increasing emissions from emerging economies like China, with relatively limited emissions controls, are outpacing reductions in the developed economies. This scenario would probably increase tropospheric ozone, which should be observable in the midtroposphere above western North America. This region is downwind of Asian emissions and also receives aged emissions from Europe and North America that have circled the northern hemisphere.

At the surface significant daytime ozone increases (0.19-0.50 parts per billion by volume per year, p.p.b.v. yr -1 ) have occurred at several rural sites in the western USA based on 1987-2004 yearly data [START_REF] Jaffe | Increase in surface ozone at rural sites in the western US[END_REF] , while springtime ozone increased by 0.33-0.59 p.p.b.v. yr -1 . However, the causes of the ozone increases at these sites were not clearly identified. Merged ozone records from several US west coast sites show that springtime ozone increased by 0.46 p.p.b.v. yr -1 (1985-2007) when winds from the North Pacific Ocean were strongest, indicating an increase in marine boundary layer background ozone [START_REF] Parrish | Increasing ozone in marine boundary layer inflow at the west coasts of North America and Europe[END_REF] . These rates are much higher than the surface ozone rate of increase (Northern Hemisphere annual average) attributed by chemical transport models to rising Asian emissions [START_REF] Fiore | Multimodel estimates of intercontinental source-receptor relationships for ozone pollution[END_REF] .

Regarding the free troposphere above the western USA, no previous study had found a significant ozone increase, despite frequent impact by Asian pollution plumes [START_REF] Stohl | On the pathways and timescales of intercontinental air pollution transport[END_REF] . Ozone measurements in the western North American free troposphere are relatively few, with the longest records being the weekly ozonesonde profiles above Edmonton, Alberta; Boulder, Colorado; and Trinidad Head, California. No significant ozone trend occurred above Edmonton during 1980-2001, but ozone did increase significantly during 1991-2001 [START_REF] Tarasick | Changes in the vertical distribution of ozone over Canada from ozonesondes: 1980-2001[END_REF] . No significant increase was found above Boulder (1985-2004) [START_REF] Oltmans | Long-term changes in tropospheric ozone[END_REF] or Trinidad Head (1997-2007) [START_REF] Oltmans | Background ozone levels of air entering the west coast of the US and assessment of longer-term changes[END_REF] . Upper tropospheric ozone measurements made by commercial aircraft in the late 1970s and late 1990s revealed significant increases above the eastern USA, among other regions, but not above western North America where small sample sizes in the late 1990s may have prevented the identification of significant changes [START_REF] Schnadt Poberaj | Long-term changes in UT/LS ozone between the late 1970s and the 1990s deduced from the GASP and MOZAIC aircraft programs and from ozonesondes[END_REF] . In contrast, mid-tropospheric ozone measurements made by research aircraft above coastal California during the springtime of 1984 and 2002 show an increase of 15 p.p.b.v. However, this apparent increase may simply have been due to interannual variability [START_REF] Parrish | Changes in the photochemical environment of the temperate North Pacific troposphere in response to increased Asian emissions[END_REF] .

Our goal is a reanalysis of all available free tropospheric ozone measurements above western North America during springtime of 1984-2008 to quantify the change in ozone using the largest data set possible. The analysis was restricted to April-May, when Asian emissions have their greatest impact on western North America 9 .E m p h a s i sw a s placed on the mid-troposphere (3.0-8.0 km) where the most data are available, but the upper troposphere (8.0-12.0 km) was also addressed. Ozone measurements were gathered from ozonesondes, lidar and commercial and research aircraft across western North America (25u-55u N, 130u-90u W). A particle dispersion model (PDM) was used to determine the 15-day transport history of every measurement (see Methods).

Continuous ozone data with adequate sample sizes (see Methods) are available in the mid-troposphere during 1995-2008 (Fig. 1a), where median ozone values have a significant rate of increase of 0.63 6 0.34 p.p.b.v. yr -1 . Similar values are found for the 33rd and 67th percentiles, whereas rates for the 95th and 5th percentiles are not significant. The only other year with a large sample size is 1984, with a median ozone value (46 p.p.b.v.) that is 9-20 p.p.b.v. less than any year in the period 1995-2008. When 1984 is included in the analysis, the ozone rate of increase is slightly greater for all percentiles (Fig. 1a), suggesting that the 1995-2008 ozone rate of increase can be extended to 1984. Transport patterns do vary interannually, but PDM analysis shows that no transport anomaly occurred in 1984, indicating that mid-tropospheric ozone was substantially lower in 1984 owing to reduced photochemical production rather than unusual transport. No other large April-May data sets are available from the 1980s to corroborate the 1984 measurements. However, mid-tropospheric ozone data measured above the western USA by a research aircraft in June 1985 [START_REF] Luke | Tropospheric chemistry over the lower Great Plains of the United States. 2: Trace gas profiles and distributions[END_REF] and May-June 1986 [START_REF] Luke | Reactive Nitrogen Compounds in the Troposphere: Observations, Transport, and Photochemistry[END_REF] have mean values of 47 and 53 p.p.b.v., respectively. Although these measurements occur later in the year, and have not undergone the same data averaging and processing as applied to the data in Fig. 1a, their low values are consistent with the 1984-2008 rate of increase found here.

Figure 1b-g shows the average transport history of air parcels in both the atmospheric column and the 300 m layer adjacent to the Earth's surface, associated with the lowest, middle and highest third of ozone measurements. The lowest ozone values have source regions to the west, stretching from the mid-latitudes to the tropics, with surface source regions mainly in the sub-tropical North Pacific. As ozone increases the source region shifts further west with greater influence from east Asia. Influence from the westernmost North Pacific is notable because it contains one of the world's major shipping routes, with global shipping producing about 13% of anthropogenic NO x emissions [START_REF]RETRO Emission Data Sets and Methodologies for Estimating Emissions, Work Package 1[END_REF] .

The analysis so far contains ozone measurements that experienced varying degrees of transport from the North American boundary layer. To focus on the background air entering western North America, measurements with a strong and recent influence (,5 days) from the North American boundary layer were identified by the PDM and removed (Supplementary Information). This procedure resulted in a median ozone rate of increase of 0.71 6 0.45 p.p.b.v. yr -1 for 1995-2008 (compared to 0.63 6 0.34 p.p.b.v. yr [START_REF] Oltmans | Background ozone levels of air entering the west coast of the US and assessment of longer-term changes[END_REF] for the unfiltered data), and 0.76 6 0.29 p.p.b.v. yr -1 for 1984-2008 (compared to 0.70 6 0.22 p.p.b.v. yr -1 for the unfiltered data). The stronger rate of increase found in the background air demonstrates that recent North American influence is not responsible for the increasing ozone. Given the decreasing North American ozone precursor emissions since 1985 [START_REF] Van Der | Trends, seasonal variability and dominant NO x source derived from a ten year record of NO 2 measured from space[END_REF][START_REF] Jacob | Effect of rising Asian emissions on surface ozone in the United States[END_REF][START_REF]National Emissions Inventory (NEI) Air Pollutant Emissions Trends Data AEhttp[END_REF] , and the decrease in urban extreme ozone events across much of the USA [START_REF] Lefohn | Characterizing long-term changes in surface ozone levels in the United States (1980-2005)[END_REF] , it is unlikely that aged North American emissions are contributing to the increase of background ozone. PDM output can also be used to identify trends in the frequency of transport from the important ozone precursor emission region of south Asia and east Asia (referred to as SA/EA and including India, southeast Asia and China). Although the quantity of air transported from the SA/EA boundary layer to the North American free troposphere varied interannually, we found no significant trend during 1995-2008.

To remove the influence of interannual variability associated with transport from various emission regions we examined the ozone rate of increase by dividing the full 1995-2008 ozone data set (Fig. 1a) according to weaker or stronger transport from a particular emission region. Figure 2 shows the residence time of the air parcels associated with weaker or stronger transport from the anthropogenic NO x sources within SA/ EA. Ozone data with weaker transport from SA/EA have a strong association with transport from the central North Pacific and increase by 0.45 6 0.32 p.p.b.v. yr -1 (P 5 0.01), which is less than the overall rate of 0.63 6 0.34 p.p.b.v. yr -1 (Fig. 1a). In contrast, data with stronger transport from SA/EA have a strong surface influence from China and subtropical regions of the western North Pacific, with an additional influence from southeast Asia and northern India, and increase at a higher rate of 0.80 6 0.34 p.p.b.v. yr -1 (P 5 0.00). Weaker or stronger transport from China alone yields similar results of 0.42 6 0.39 p.p.b.v. yr -1 (P 5 0.04) and 0.72 6 0.36 p.p.b.v. yr -1 (P 5 0.00), respectively. Because the direct transport routes from India and southeastern Asia to western North America traverse China [START_REF] Stohl | On the pathways and timescales of intercontinental air pollution transport[END_REF] , this analysis cannot quantify the relative contributions of China, India and southeast Asia to the ozone increase. SA/ EA was the only region in the Northern Hemisphere associated with an ozone rate of increase greater than the overall rate of 0.63 p.p.b.v. yr -1 .For example, strong transport from Japan/Korea has an ozone rate of increase of 0.52 6 0.40 p.p.b.v. yr -1 (P 5 0.02), much lower than when transport from SA/EA is strong. The strong transport pattern from Japan/Korea has relatively long surface residence times above the western North Pacific, and shorter residence times above continental Asia (Supplementary Information). This small shift in transport away from continental Asia has a large impact on the ozone rate of increase above western North America.

Because this study shows a stronger ozone rate of increase above North America when air masses originate in the rapidly growing emission regions of SA/EA, recent ozone changes in Asia deserve attention (Supplementary Information). The limited long-term observations within Asia [START_REF] Chou | The trend of surface ozone in Taipei, Taiwan, and its causes: implications for ozone control strategies[END_REF][START_REF] Ding | Tropospheric ozone climatology over Beijing: analysis of aircraft data from the MOZAIC program[END_REF][START_REF] Tanimoto | Increase in springtime tropospheric ozone at a mountainous site in Japan for the period 1998-2006[END_REF] show ozone increases of 0.58-1.2 p.p.b.v. yr -1 at surface sites in China, Japan and Taiwan, comparable to our rate of 0.80 p.p.b.v. yr -1 above North America when transport from Asia is strong. Recent model analysis indicates that about half of the ozone transported from Asia to North America is produced in Asia and the rest is produced during trans-Pacific transport [START_REF] Zhang | Transpacific transport of ozone pollution and the effect of recent Asian emission increases on air quality in North America: an integrated analysis using satellite, aircraft, ozonesonde, and surface observations[END_REF] .

Lastly, the ozone rate of increase was calculated for the upper troposphere, using all available data during 1995-2008. Only the 50th and 33rd percentiles show significant increases in ozone-0.58 6 0.52 and 0.66 6 0.52, respectively-with values and transport history similar to the mid-tropospheric analysis. The strong influence of stratospheric intrusions in the upper troposphere probably explains the lack of significant rates of increase in the 67th and 95th percentiles. Further data analysis in the Supplementary Information shows that the mid-and upper tropospheric results are robust and not strongly influenced by air mass sampling biases, vertical ozone variation, or stratospheric impact.

We explored the discrepancy between this study, which shows a statistically significant ozone increase above western North America, and previous ozonesonde analyses which do not. Mid-tropospheric ozone rates of change above Trinidad Head (1998-2008), Boulder (1995-2008) and Edmonton (1995-2008), were 0.37 6 0.58 p.p.b.v. yr -1 , 0.05 6 0.46 p.p.b.v. yr -1 , and 0.12 6 0.52 p.p.b.v. yr -1 , respectively. Although these values are positive, they are much lower than the rate in Fig. 1a and statistically insignificant. The annual sample sizes for the ozonesonde sites are fairly low (32-564 data points), with larger sample sizes occurring in more recent years at Edmonton and Trinidad Head. The lack of significant rates of increase from ozonesondes could be due to small sample sizes in a region of high springtime ozone variability, requiring a much longer data record to detect any significant increase [START_REF] Jaffe | Increase in surface ozone at rural sites in the western US[END_REF][START_REF] Logan | Trends in the vertical distribution of ozone: a comparison of two analyses of ozonesonde data[END_REF] . Removing the ozonesondes from the full data set has little impact on this study, slightly increasing the 1995-2008 midtropospheric ozone rate of change from 0.63 6 0.34 p.p.b.v. yr -1 to 0.69 6 0.43 p.p.b.v. yr -1 .

This study has implications for North American air quality [START_REF] Jacob | Effect of rising Asian emissions on surface ozone in the United States[END_REF][START_REF] Fiore | Multimodel estimates of intercontinental source-receptor relationships for ozone pollution[END_REF] and quantification of the radiative forcing due to ozone 2 . Our data provide an important benchmark for model simulations of free tropospheric ozone downwind of Asia. Models that can reproduce the observed ozone rate of increase in this region will provide more accurate estimates of changes in air quality and ozone radiative forcing since pre-industrial times.

We suggest that the free tropospheric ozone increase could also affect the surface. Previous studies show that polluted midtropospheric air entering western North America can descend to the surface and modify the composition of the boundary layer (see 1984 1986 1988 1990 1992 1994 1996 1998 2000 2002 2004 review in the Supplementary Information). Our results support earlier work that indicates that rising Asian ozone precursor emissions would cause springtime surface ozone to increase in western North America since the 1980s, despite decreasing domestic emissions [START_REF] Jacob | Effect of rising Asian emissions on surface ozone in the United States[END_REF] .F i n a l l y , summertime extreme ozone events in many US urban areas have decreased [START_REF] Lefohn | Characterizing long-term changes in surface ozone levels in the United States (1980-2005)[END_REF] , while some rural and marine sites in the western US show increasing ozone [START_REF] Jaffe | Increase in surface ozone at rural sites in the western US[END_REF][START_REF] Parrish | Increasing ozone in marine boundary layer inflow at the west coasts of North America and Europe[END_REF] , possibly due to increasing background ozone. Future studies are required to quantify free tropospheric ozone trends above North America in other seasons.

METHODS SUMMARY

The ozone data were measured by: ( of the true ozone median value, at the 99% confidence level. The data have an approximately normal distribution and the straight-line fits through the ozone percentiles (Fig. 1a) are calculated using the least-squares method of simple linear regression. The slopes of these lines indicate the increase of ozone per year for 1984-2008. The P value indicates the statistical significance of the linear relationship, determined by first calculating R, the correlation coefficient between ozone and time. We then test the null hypothesis that R 2 5 0 (no linear relationship) using the standard F-statistic (ratio of the mean square regression to the mean square residual). If the probability P associated with the F statistic is small (P # 0.05), the null hypothesis is rejected with a confidence level $95%.

The FLEXPART Lagrangian particle dispersion model produced the 15-day retroplumes for each ozone measurement with back-trajectory calculations based upon global ECMWF wind fields with a temporal resolution of 3 h, horizontal resolution of 1u 3 1u, and 60 (90) vertical levels in 1984-2005 (2006-2008). Back-trajectory particles are transported by the resolved winds, parameterized sub-grid motions, and a parameterized convection scheme. Further details on FLEXPART and the ozone data sets are given in the Supplementary Information.

O

  . R. Cooper 1,2 , D. D. Parrish 2 , A. Stohl 3 , M. Trainer 2 ,P.Ne ´de ´lec 4 , V. Thouret 4 , J. P. Cammas 4 , S. J. Oltmans 2 , B. J. Johnson 2 , D. Tarasick 5 , T. Leblanc 6 , I. S. McDermid 6 , D. Jaffe 7 , R. Gao 2 , J. Stith 8 , T. Ryerson 2 , K. Aikin 1,2 , T. Campos 9 , A. Weinheimer 9 & M. A. Avery 10

Figure 1 |

 1 Figure 1 | Springtime ozone distributions for 1984, 1995-2008 in the midtroposphere (3.0-8.0 km), and air mass source regions. a, Distributions of springtime ozone measurements madeinthetropospherebetween3.0and 8.0 km (stratospheric samples have been filtered out). The green line and data points are the median, and the yellow data points are means. The upper and lower blue lines (and data points) indicate the 95th and 5th percentiles. The upperandlowerredlines(anddatapoints)indicatethe67thand33rd percentiles. Ozone sample sizes range from 1,663 in 1984 to 8,587 in 2006 (see the Supplementary Information). Also shown are the ozone rates of increase for 1984-2008 and 1995-2008, as determined from the slope of the linear regression. The range on the slope indicates the 95% confidence limit that the slope lies within that range. Ozone data were gathered over mid-latitude western North America (25u-55u N, 130u-90u W), as shown in the Supplementary Information. The transport history of each ozone measurement was determined by calculating a retroplume with the FLEXPART PDM (see Methods and Supplementary Information). Every retroplume consisted of 40,000 back-trajectory particles released from the time and location of each measurement and advected backwards in time for 15 days. b-d,Theaverage 1984-2008 retroplume for three ranges of ozone measurements, expressed as column residence times. e-g, The corresponding retroplume residence times in the lowest 300 m of the atmosphere (the footprint layer). Ozone percentile ranges: b and e,0-33rd,c and f,34th-66th;d and g,67th-99th.Columnand footprint sample sizes are equal because every 15-day retroplume has some degree of transport through the lowest 300 m of the atmosphere.

Figure 2 |

 2 Figure 2 | Average 1995-2008 FLEXPART retroplume residence times.Average retroplume residence times in the atmospheric column (a and b) and footprint layer (c and d) for weaker (a and c) and stronger (b and d) transport from the SA/EA region. The retroplumes associated with every 1995-2008 ozone measurement in Fig.1awere ranked according to the amount of time spent in the boundary layer above the SA/EA emissions regions. The 50% of retroplumes with the longest residence times were categorized as having stronger transport from SA/EA, while the rest were categorized as having weaker transport from SA/EA. For each group of retroplumes the average column and footprint was calculated (1995-2008), and for each group the ozone rate of change was also calculated (reported in the main text).

Cooperative Institute for Research in Environmental Sciences, University of Colorado, Boulder, Colorado 80309, USA.

NOAA Earth System Research Laboratory,

Broadway, Boulder, Colorado 80305, USA.[START_REF] Horowitz | Past, present, and future concentrations of tropospheric ozone and aerosols: methodology, ozone evaluation, and sensitivity to aerosol wet removal[END_REF] Department of Regional and Global Pollution Issues, Norwegian Institute for Air Research, PO Box 100, N-2027 Kjeller, Norway.[START_REF] Staehelin | Trends in surface ozone concentrations at Arosa (Switzerland)[END_REF] Laboratoire d'Ae ´rologie, Centre National de la Recherche Scientifique, Observatoire Midi-Pyre ´ne ´es, 14 Avenue E. Belin, 31400 Toulouse, France.[START_REF] Oltmans | Long-term changes in tropospheric ozone[END_REF] Experimental Studies Research Division, Meteorological Service of Canada, Environment Canada, 4905 Dufferin Street, Downsview, Ontario M3H 5T4, Canada.[START_REF] Thouret | Tropopause referenced ozone climatology and inter-annual variability (1994-2003) from the MOZAIC programme[END_REF] Table Mountain Facility, Jet Propulsion Laboratory, California Institute of Technology, 24490 Table Mountain Road, Wrightwood, California 92397-0367, USA. 7 Department of Atmospheric and Environmental Chemistry, University of Washington-Bothell, 18115 Campus Way NE, Bothell, Washington 98011-8246, USA. 8 NCAR Research Aviation Facility, National Center for Atmospheric Research, 10802 Airport Court, Broomfield, Colorado 80021, USA. 9 Atmospheric Chemistry Division, National Center for Atmospheric Research, PO Box 3000, Boulder, Colorado 80305-3000, USA. 10 NASA Langley Research Center, Hampton, Virginia 23681-2199, USA.

Acknowledgements This work was supported in part by NOAA's Climate Goal Program. We acknowledge the support of MOZAIC by the European Communities, EADS, Airbus and the airlines (Lufthansa, Austrian, Air France) who have carried MOZAIC equipment free of charge since 1994. B. Ridley, NCAR (retired), measured ozone from the NCAR C-130 during TOPSE. G. L. Gregory, NASA (retired), measured ozone from the Convair CV-990 during CITE-1C. M. Proffitt, NOAA (retired) measured ozone from the NASA ER2 during STRAT and POLARIS, and measured ozone from the NASA WB57 during WAM. PACDEX ozone data were provided by NCAR/EOL under sponsorship of the National Science Foundation (http://data.eol.ucar.edu/). EDGAR (http://www.mnp.nl/edgar) is a product of the National Institute for Public Health and the Netherlands Organisation for Applied Scientific Research and is part of the Global Emissions Inventory Activity of IGBP/IGAC. Finally we thank R. Dickerson for providing the mean ozone values for June 1985 and May-June 1986. Author Contributions O.R.C., D.D.P., A.S. and M.T conceived the study and guided its development. O.R.C. merged the ozone datasets and ran the PDM with assistance from A.S. O.R.C. analysed the data and wrote the text. P.N., V.T., J.P.C., S.J.O., B.J.J., D.T., T.L., I.S.M., D.J., R.G., J.S., T.R., K.A., T.C., A.W. and M.A.A. made the ozone measurements and guided interpretation of the datasets.

Supplementary Information is linked to the online version of the paper at www.nature.com/nature.

Author Information Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests. Correspondence and requests for materials should be addressed to O.R.C. (owen.r.cooper@noaa.gov).


Shaohui Zhang
email: shaohui.zhang@uef.fi

Cédric Vega
email: cedric.vega@ign.fr

Olivier Bouriaud
email: olivier.bouriaud@ign.fr

Sylvie Durrieu
email: sylvie.durrieu@inrae.fr

Jean-Pierre Renaud
email: jean-pierre.renaud-02@onf.fr

Unit-level Small Area Estimation of Forest Inventory with GEDI Auxiliary Information

Introduction

National Forest Inventories (NFIs) play an important role in understanding the state of forests at the national and regional levels. Forest inventory for small territorial areas, such as municipalities, is also important for decision-makers. However, information is relatively limited at this level. As a result, developing small area estimation (SAE) approaches has gained increasing popularity in the field of forest inventory. It enables prediction of forest attributes for sub-populations using regression models based on auxiliary data commonly derived from remote sensing techniques over an area of interest (AOI). It has been reported that SAE can improve the precision of forest inventory without increasing costs [START_REF] Mandallaz | New regression estimators in forest inventories with two-phase sampling and partially exhaustive information: a design-based monte carlo approach with applications to small-area estimation[END_REF] and may produce reliable predictions of forest attributes locally, even when field plots are not available [START_REF] Rao | Small-area estimation[END_REF]. [START_REF] Tomppo | The Finnish multi-source national forest inventory-small area estimationand map production[END_REF] is a pioneer in the use of auxiliary data for multi-source forest inventory. Previously, common sources of auxiliary data often came from satellite-based imagery [START_REF] Mcroberts | Estimating areal means and variances of forest attributes using the k-Nearest Neighbors technique and satellite imagery[END_REF]), digital aerial photogrammetry [START_REF] Breidenbach | Unit-level and area-level small area estimation under heteroscedasticity using digital aerial photogrammetry data[END_REF], and airborne laser scanning [START_REF] Magnussen | National forest inventories in the service of small area estimation of stem volume[END_REF]). NASA's newly-launched Global Ecosystem Dynamics Investigation (GEDI) is a full waveform LiDAR instrument aboard the International Space Station (ISS). Its products consist of footprint measurements projected to cover 4% of the global land surface by the end of its mission [START_REF] Dubayah | The Global Ecosystem Dynamics Investigation: High-resolution laser ranging of the Earth's forests and topography[END_REF]. This will provide an unprecedented opportunity to systematically collect samples of forest information that can be used in SAE on a large scale.

The objective of this study is to explore the possibility of using GEDI auxiliary data to improve the accuracy of forest inventory for a large natural area in central France (Sologne), as well as for smaller sub-areas defined by French administrative boundaries (departments). The results will then be compared against estimates obtained from simple random sampling (SRS), to assess the efficiency of the auxiliary data.

Data and Methods

Study Area

Our study is based in Sologne, central France, which covers an area of approximately 6000km 2 . The topography is mostly flat, with most elevations within the range of 110-250m, except the south-eastern part, where undulating terrain reaches 400m. Forests cover approximately 48% of the area and are dominated by pure broadleaved species (75.3%). Conifer and mixed stands account for 15.5% and 9.2% of the forest areas respectively. The climate is temperate Atlantic, with mean annual temperature and precipitation of 11°C and 725mm.

NFI Data

For the study area, 902 permanent NFI plots, surveyed between 2015 and 2019, were available. This five-year timeframe is routinely used in official French NFI statistics. Each plot contains detailed inventory information, including density (trees/ha), quadratic mean diameter (cm), basal area (m2/ha), dominant height (m), and volume (m3/ha). However, this paper focuses solely on the aspect of forest volume estimations. Details of inventory schemes and methods can be found in [START_REF] Hervé | L'inventaire des ressources forestières en France: un nouveau regard sur de nouvelles forêts[END_REF].

Auxiliary Data Sources and Processing

GEDI Level 2A products, acquired between 2019-04-22 and 2020-04-14, were used as auxiliary data. Each footprint is of 25-meter diameter and has the information of beam types (i.e., full or half power), sensitivity, geo-located elevation and height metrics. Only the best quality footprint data were retained based on the following filtering criteria: Firstly, we only kept footprints whose quality flag is 1, degrade flag is 0, and sensitivity is ≥ 0.9. This is a standard filtering criterion applied to select suitable GEDI footprints, which also ensures that the beam power is strong enough to penetrate the canopy and reach the ground with 90% probability [START_REF] Hancock | The GEDI Simulator: A Large-Footprint Waveform Lidar Simulator for Calibration and Validation of Spaceborne Missions[END_REF]). Next, we followed the definition of forest from FAO, which states that "trees in a forest reach a minimum height of 2-5 m at maturity" (FAO 2021). Therefore, we used an average height of 3m as threshold and removed those footprints whose RH100 values were smaller than 3 and thus did not qualify as forest. As a result, a total number of 112,569 footprints were included for further analysis. In addition, the footprints were intersected with forest masks to determine in which types of forest they were located. Height metrics and forest types were then extracted from the footprints and formed the auxiliary data frame.

Lastly, seven nearest neighbouring GEDI footprints around individual NFI plots were identified based on their Euclidean distances. We set an additional distance threshold of 200m to filter out those NFI plots that had neighbouring footprints located farther than the threshold and thus may misrepresent the plot information. Based on the identified shot numbers, the remaining 105 NFI plots were then joined with one of the seven neighbouring footprints that shared the closest forest height, as defined by the smallest value of |NFI dominant height -RH100|. This formed the calibration data frame.

Unit-level Small Area Estimation

Unit-level SAE was performed using the two-phase estimation procedure provided in R package "forestinventory" and described in [START_REF] Hill | The R Package forestinventory: Design-Based Global and Small Area Estimations for Multiphase Forest Inventories[END_REF]. The first phase of auxiliary GEDI information was used to generate model predictions based on linear regression using the method of ordinary least square. The second phase contains the targeted NFI plot attributes, i.e., forest volume alone in this case, that is used to generate model coefficients and correct bias.

Results and Discussion

A total of 101 auxiliary GEDI variables (100 height metrics and one forest type) were available and tested to predict forest volumes. Variable selection was done using an exhaustive search with the help of "randomForestSRC" R-package [START_REF] Ishwaran | RandomForestSRC: Fast Unified Random Forests for Survival, Regression, and Classification (RF-SRC)[END_REF]. The most relevant variables were manually verified to yield the best model fit and a low variance inflation factor (< 5). The final linear model retained was:

Ŷ = 𝛽0 + 𝛽1 * RH100 + 𝛽2 * RH20 + 𝛽3 * Forest Type (1)
Loir-et-Cher Loiret Cher

Departments in Sologne

Where Ŷ is the predicted volume, 𝛽0 is the intercept, 𝛽1, 𝛽2 and 𝛽3 represent coefficients. RH100 and RH20 represent respectively the relative heights at which the maximum and the 20 th percentile of waveform energy above the ground peak are reached.

Based on this model formulation, SAE of forest volume was performed at the Sologne and sub-area levels. Results showed that GEDI auxiliary information significantly improved estimate accuracy compared with results obtained without this auxiliary information (Table 1). At the whole AOI level, the variance was significantly reduced and an increase in relative efficiency by a factor of 2.6 was obtained. In each sub-area, similar results were achieved, with reduced variance and an increase in relative efficiency by factors varying between 1.6 and 2.6. However, the mean forest volume was somehow underestimated using SAE estimation with the help of GEDI auxiliary data at both AOI and sub-area levels. We further calibrated the model and discovered that this was likely caused by model extrapolation, as a considerable number (16%) of predicted forest volumes fell outside the model calibration domain.

Conclusions

This paper performed unit-level small area estimation using GEDI Level 2A as auxiliary data. We associated GEDI auxiliary information with NFI plots based on Euclidean distance to assess forest volume estimation. It is shown that GEDI auxiliary information can help improve forest volume estimation significantly when compared to simple random sampling alone. The fact that GEDI data are open-access and cover the entire country makes it a particularly attractive tool for improving forest inventory at regional and local levels.

Figure 1 :

 1 Figure 1: Location of the Study Area (departments have different colours) and availability of GEDI footprints over the AOI.

Table 1 .

 1 Volume estimations of SAE and SRS at both AOI and department levels

	Area	Plot N	SRS Estimation	SRS Variance	SAE Estimation	SAE Variance	Relative Efficiency
	Overall AOI	105	192.0 ± 25.6	164.3	170.2 ± 15.8	62.4	2.6
	Cher	20	198.0 ± 51.1	653.0	185.3 ± 40.2	404.0	1.6
	Loiret	46	224.1 ± 43.6	475.5	181.7 ± 27.2	224.1	2.6
	Loir-et-Cher	39	151.1 ± 34.8	302.2	150.6 ± 22.0	120.4	2.5


Canesi

A Miglioli

Teresa Balbi
email: teresa.balbi@unige.it

M Montagna

R Dumollard B#

L Canesi

Tetrabromobisphenol A acts a neurodevelopmental disruptor in early larval stages of Mytilus galloprovincialis

Keywords: tetrabromobisphenol A, Mytilus, early larval stages, shell formation, neuroendocrine disruptor

  L'archive ouverte pluridisciplinaire

Introduction

Tetrabromobisphenol A (TBBPA) is one of the most widely used brominated flame retardants in the production of epoxy and polycarbonate resins employed in electronic devices, furniture, plastics, textiles, and other products [START_REF] Covaci | Analytical and environmental aspects of the flame retardant tetrabromobisphenol-A and its derivatives[END_REF]. TBBPA is highly produced worldwide, and its global market is expected to increase over the next five years (WorldAnalytics, 2019). EPA designated TBBPA as a high priority chemical currently undergoing risk evaluation (EPA, 2019).

TBBPA is frequently detected in different ecosystem compartments, including aquatic environments [START_REF] Liu | A review of status of tetrabromobisphenol A (TBBPA) in China[END_REF][START_REF] Zhou | Tetrabromobisphenol A (TBBPA): a controversial environmental pollutant[END_REF] and it can be also accumulated in the food chain [START_REF] Morris | Distribution and fate of HBCD and TBBPA brominated flame retardants in North Sea estuaries and aquatic food webs[END_REF][START_REF] Yang | Tetrabromobisphenol A: tissue distribution in fish, and seasonal variation in water and sediment of Chaohu Lake, China[END_REF][START_REF] Shi | A national survey of tetrabromobisphenol-A, hexabromocyclododecane and decabrominated diphenyl ether in human milk from China: Occurrence and exposure assessment[END_REF].

TBBPA is likely to reach the marine environment largely through industrial waste waters from landbased industrial activities, and coastal areas worldwide are subjected to increasing TBBPA contamination, thus raising great environmental concern [START_REF] Gu | Aerobic cometabolism of tetrabromobisphenol A by marine bacterial consortia[END_REF]. Levels of TBBPA in urbanized coastal areas range from ng to µg in both sediments (g -1 dw) and waters (L -1 ) [START_REF] Liu | A review of status of tetrabromobisphenol A (TBBPA) in China[END_REF][START_REF] Sühring | Fingerprint analysis of brominated flame retardants and Dechloranes in North Sea sediments[END_REF]Gong et al., 2017 and refs therein). High concentrations (>1 µg L -1 ) detected at polluted sites [START_REF] Liu | A review of status of tetrabromobisphenol A (TBBPA) in China[END_REF][START_REF] Gong | The occurrence and spatial-temporal distribution of tetrabromobisphenol A in the coastal intertidal zone of Qingdao in China, with a focus on toxicity assessment by biological monitoring[END_REF][START_REF] Li | A verification of correlation between chemical monitoring and multibiomarker approach using clam Ruditapes philippinarum and scallop Chlamys farreri to assess the impact of pollution in Shandong coastal area of China[END_REF] are comparable to TBBPA levels in urban soils and freshwater, suggesting that coastal seawater may receive discharge from industrial effluents or landfill leakages without effective treatment.

Due to its almost ubiquitous presence, the potential impact of TBBPA on human and environmental health has been widely investigated [START_REF] Wang | Contamination level, distribution characteristics, and ecotoxicity of tetrabromobisphenol A in water and sediment from Weihe river basin, China[END_REF][START_REF] Zhou | Tetrabromobisphenol A (TBBPA): a controversial environmental pollutant[END_REF]. Studies based on cell and animal models (mammals and fish) have demonstrated that TBBPA can induce pleiotropic effects on liver, kidney, heart, neural and reproductive functions; in particular, TBBPA has been shown to be a potential endocrine disruptor (ED), because of its reproductive and developmental effects including interference with thyroid hormone signalling [START_REF] Zhang | Tetrabromobisphenol A disrupts vertebrate development via thyroid hormone signaling pathway in a developmental stage-dependent manner[END_REF][START_REF] Zhou | Tetrabromobisphenol A (TBBPA): a controversial environmental pollutant[END_REF]. Significant effects of TBBPA on early embryo development have been observed in mammalian and zebrafish models in the nmol/L range, underlying that particular attention should be paid to the potential effects of TBBPA on early developmental stages [START_REF] Zhou | Tetrabromobisphenol A (TBBPA): a controversial environmental pollutant[END_REF]. Moreover, it has been shown that molluscan shell growth may be uniquely sensitive to TBBPA, more sensitive than chronic fish or crustacean toxicity endpoints [START_REF] Pittinger | Review of historical aquatic toxicity and bioconcentration data for the brominated flame retardant tetrabromobisphenol A (TBBPA): effects to fish, invertebrates, algae, and microbial communities[END_REF].

We have previously demonstrated that in the model marine bivalve, the mussel Mytilus galloprovincialis, environmental concentrations of TBBPA affected early larval development evaluated in the standard 48 h embryotoxicity assay, with a LOEC of 0.1 µg/L and an EC50 of 5.52 µg/L [START_REF] Fabbri | Adaptation of the bivalve embryotoxicity assay for the high throughput screening of emerging contaminants in Mytilus galloprovincialis[END_REF]. The effects were comparable with those observed with the parent compound Bisphenol A-BPA [START_REF] Fabbri | Adaptation of the bivalve embryotoxicity assay for the high throughput screening of emerging contaminants in Mytilus galloprovincialis[END_REF], that has been subsequently shown to interfere with the molecular mechanisms involved in shell biogenesis and development of the serotoninergic system in early mussel larvae [START_REF] Balbi | Impact of bisphenol A (BPA) on early embryo development in the marine mussel Mytilus galloprovincialis: Effects on gene transcription[END_REF]Miglioli et al., 2021a).

In this work, the effects and mechanisms of action of TBBPA were further investigated in mussel larvae at different times post fertilization, with special emphasis on early shell formation, the key step in the transition from the trocophora (24 hpf) to the first D-veliger stage (48 hpf), as well as on neurodevelopment. Concentrations were chosen among those encompassing the EC50 values obtained in the 48 h larval developmental assay (1 and 10 µg/L, depending on the endpoint measured) [START_REF] Fabbri | Adaptation of the bivalve embryotoxicity assay for the high throughput screening of emerging contaminants in Mytilus galloprovincialis[END_REF], and environmental levels detected in polluted coastal areas [START_REF] Liu | A review of status of tetrabromobisphenol A (TBBPA) in China[END_REF][START_REF] Gong | The occurrence and spatial-temporal distribution of tetrabromobisphenol A in the coastal intertidal zone of Qingdao in China, with a focus on toxicity assessment by biological monitoring[END_REF][START_REF] Li | A verification of correlation between chemical monitoring and multibiomarker approach using clam Ruditapes philippinarum and scallop Chlamys farreri to assess the impact of pollution in Shandong coastal area of China[END_REF]. Higher, embryotoxic concentrations (100 µg/L) [START_REF] Fabbri | Adaptation of the bivalve embryotoxicity assay for the high throughput screening of emerging contaminants in Mytilus galloprovincialis[END_REF] were also utilized to evaluate the effects on shell formation at earlier stages.

The effects of TBBPA on shell biogenesis were investigated by different techniques, following the experimental setup previously described [START_REF] Miglioli | Characterization of the main steps in first shell formation in Mytilus galloprovincialis: possible role of tyrosinase[END_REF](Miglioli et al., , 2021a)). Deposition of the organic matrix and of the calcified shell were evaluated at 24, 28, 32, and 48 hpf by Calcofluor/Calcein staining, and quantified by digital imaging. In addition, shell morphology at 48 hpf was investigated by Scanning Electron Microscopy-SEM. Across different developmental stages, expression of key genes involved in shell matrix deposition (chitin synthase-CS; tyrosinase-TYR) and calcification (carbonic anhydrase-CA; extrapallial protein-EP) was quantified by qPCR; the expression pattern of tyrosinase and chitinase was also evaluated by In Situ Hybridization-ISH.

Moreover, the possible effects of TBBPA on neurodevelopment were investigated focusing on the main neurotransmitters identified in bivalves [START_REF] Liu | D1 dopamine receptor is involved in shell formation in larvae of Pacific oyster Crassostrea gigas[END_REF][START_REF] Joyce | Molluscan bivalve settlement and metamorphosis: Neuroendocrine inducers and morphogenetic responses[END_REF]. Localization of components of the dopaminergic system (Tyrosine Hydroxylase-TH, Dopamine-β-Hydroxylase-DβH, Dopamine Receptor 1-DR1) was investigated by ISH; development of serotoninergic and GABAergic neurons was evaluated by immunohistochemistry.

Methods

Mussels, larval development and exposure conditions

All procedures were carried out as previously described [START_REF] Miglioli | Characterization of the main steps in first shell formation in Mytilus galloprovincialis: possible role of tyrosinase[END_REF](Miglioli et al., , 2021a)). Sexually mature specimens of M. galloprovincialis were collected from a natural population in the Bay of Villefranche-sur-mer (43.682°N, 7.319°E -France) during the spawning season of 2019/20 (December-March). Animals were maintained at the Institut de la Mer de Villefranche (IMEV) by the Centre de Ressources Biologiques Marines of the institute (CRBM), transferred to the laboratory and acclimatized in flow-through vessels containing filtered natural seawater MFSW (Millipore filtered seawater) (pH 7.9-8.0, 38 ppt salinity, 15°C).

Gametes were obtained by spontaneous spawning; after checking egg quality (shape, size) and sperm motility, fertilization was carried out with an egg:sperm ratio of 1:10. After 30 min, fertilization success (n. fertilized eggs/n. total eggs x 100) was verified by microscopic observation (usually > 90%). Fertilized eggs were transferred to 24-well plates (200 larvae mL -1 ) and exposed to TBBPA 3,3′,5,5′-Tetrabromobisphenol A (Sigma-Aldrich, Milan), from stock solutions (0.5 M in Dimethylsulfoxide/DMSO) suitably diluted in filtered marine sea water-MFSW to obtain the selected final concentrations. Parallel control samples were run in MFSW added with 0.01% DMSO. In all experimental conditions, samples were represented by larvae obtained from 5 independent parental pairs (N=5), unless otherwise indicated.

TBBPA final concentrations (1 and 10 µg/L, corresponding to 1.8 and 18 nM, respectively), depending on the endpoint measured, were chosen among those encompassing higher levels detected in polluted coastal areas [START_REF] Gong | The occurrence and spatial-temporal distribution of tetrabromobisphenol A in the coastal intertidal zone of Qingdao in China, with a focus on toxicity assessment by biological monitoring[END_REF][START_REF] Li | A verification of correlation between chemical monitoring and multibiomarker approach using clam Ruditapes philippinarum and scallop Chlamys farreri to assess the impact of pollution in Shandong coastal area of China[END_REF], as well as the EC50 values previously obtained for TBBPA in the 48 h embryotoxicity assay [START_REF] Fabbri | Adaptation of the bivalve embryotoxicity assay for the high throughput screening of emerging contaminants in Mytilus galloprovincialis[END_REF]. Higher, embryotoxic concentrations (100 µg/L) were also utilized to confirm the effects on shell formation at early stages.

Concentrations in test solutions were checked by LC/MS as previously described [START_REF] Fabbri | Adaptation of the bivalve embryotoxicity assay for the high throughput screening of emerging contaminants in Mytilus galloprovincialis[END_REF].

Effect of TBBPA on shell formation

At 48 hpf, the effect of TBBPA (1, 10 and 100 µg/L) was evaluated by scoring larval phenotypes in at least 50 larvae/sample. A larva was considered normal when the shell was D-shaped (straight hinge) and the mantle did not protrude out of the shell, and malformed if the larva had not reached the stage typical for 48 h (trocophore or earlier stages) or when developmental defects were observed (concave, malformed or damaged shell, protruding mantle) [START_REF] Fabbri | Adaptation of the bivalve embryotoxicity assay for the high throughput screening of emerging contaminants in Mytilus galloprovincialis[END_REF]. Data are reported as % values of normal D-veligers. The acceptability of test results was based on control samples with a percentage of normal D-larvae >75% (ASTM, 2012).

The time-course of shell formation was also investigated at different times pf between the trochophora and the D-veliger stage (24, 28, 32 and 48 hpf). Larval shell components were visualized by Calcofluor White Fluorescent Brightener 28 (Sigma Aldrich, Lyon, France) for organic matrix, and Calcein (Sigma Aldrich, Lyon, France) for CaCO3 deposition, respectively [START_REF] Miglioli | Characterization of the main steps in first shell formation in Mytilus galloprovincialis: possible role of tyrosinase[END_REF](Miglioli et al., , 2021a)). Larvae were imaged with a Leica SP8 Confocal Laser Scanning Microscope (CLSM -Leica, France) scanning sequentially Brightfield, Calcofluor, and Calcein signals with a 0.5 µm Z-stack interval. Channels were merged, 3D rendered and rotated to measure the area (in µm 2 ) of each shell component in a single valva of each larva by manual drawing using IMAGEJ software [START_REF] Miglioli | Characterization of the main steps in first shell formation in Mytilus galloprovincialis: possible role of tyrosinase[END_REF](Miglioli et al., , 2021a)). Analyses were performed on at least 50 larvae for each sample. Data were normalized with respect to controls for each parental pair and experimental condition.

Role of thyroid hormone and dopamine signaling in early larval development

The role of thyroid hormone-TH and dopamine signaling were investigated by exposure of fertilized eggs to different TH agonists and antagonist: Triiodotyronine (T3) and Tetraiodotyronine (T4), the selective and high-affinity THR antagonist 1-850 [START_REF] Schapira | Discovery of diverse thyroid hormone receptor antagonists by high-throughput docking[END_REF], and the antithyroid drug thiourea (TU), a TH depleting compound [START_REF] Davidson | Thiourea and cyanamide as inhibitors of thyroid peroxidase: the role of iodide[END_REF] and to the selective antagonist of dopamine receptor DR1, SCH 23390 [START_REF] Liu | D1 dopamine receptor is involved in shell formation in larvae of Pacific oyster Crassostrea gigas[END_REF][START_REF] Liu | Ocean acidification inhibits initial shell formation of oyster larvae by suppressing the biosynthesis of serotonin and dopamine[END_REF]) (See SI).

Scanning electron microscopy (SEM) of larval shell

Control and TBBPA exposed larvae (1 and 10 µg/L) at 48 hpf were fixed in 3% glutaraldehyde in ASW. After fixation, samples were placed onto Whatman 22 µm filters, dehydrated in an ascending ethanol series (50% -80% -90% -100%) and air-dried. Samples were sputter-coated with gold and observed at 20 kV with a Vega3 -Tescan scanning electron microscope. All procedures were carried out as previously described [START_REF] Balbi | Impact of bisphenol A (BPA) on early embryo development in the marine mussel Mytilus galloprovincialis: Effects on gene transcription[END_REF][START_REF] Balbi | Impact of cationic polystyrene nanoparticles (PS-NH2) on early embryo development of Mytilus galloprovincialis: Effects on shell formation[END_REF][START_REF] Balbi | Diclofenac affects early embryo development in the marine bivalve Mytilus galloprovincialis[END_REF].

qPCR

Fertilized eggs were exposed to TBBPA (10 μg/L) in polystyrene 6-well plates [START_REF] Balbi | Impact of bisphenol A (BPA) on early embryo development in the marine mussel Mytilus galloprovincialis: Effects on gene transcription[END_REF]Miglioli et al., 2021a). Control samples were run in parallel. Larvae were collected at 24, 28, 32 and 48 hpf by a nylon mesh (20 μm pore-filter) and washed with filter sterilized artificial seawater (ASW) [START_REF] Balbi | Impact of bisphenol A (BPA) on early embryo development in the marine mussel Mytilus galloprovincialis: Effects on gene transcription[END_REF]. Three wells for each sample were pooled in order to obtain approximately 7000 larvae/replicate. All procedures (RNA extraction, retro-transcription and qPCR) were carried out as previously described [START_REF] Balbi | Impact of bisphenol A (BPA) on early embryo development in the marine mussel Mytilus galloprovincialis: Effects on gene transcription[END_REF][START_REF] Miglioli | Characterization of the main steps in first shell formation in Mytilus galloprovincialis: possible role of tyrosinase[END_REF]Miglioli et al., , 2021a)). Primers utilized for qPCR are reported in Tab. S1.

In Situ Hybridization

In situ hybridization (ISH) was used to evaluate the expression pattern of enzymes involved in matrix deposition (tyrosinase and chitinase) and of components of the dopaminergic system (Tyrosine Hydroxylase, Dopamine-β-Hydroxylase, Dopamine Receptor 1-DR1) following the protocol developed for M. galloprovincialis larvae [START_REF] Miglioli | Characterization of the main steps in first shell formation in Mytilus galloprovincialis: possible role of tyrosinase[END_REF](Miglioli et al., , 2021a

5-HT and GABA immunocytochemistry

Antibodies against the neurotransmitter serotonin (5-HT) and γ-Aminobutyric acid (GABA) were utilized to visualize the development of the serotonergic and GABAergic systems in both control and TBBPA-exposed larvae (10 µg/L) at different times pf as previously described (Miglioli et al., 2021a;

Vaasjo et al., 2018).

Samples were incubated with primary antibodies either against Rabbit 5-HT (1:10.000; Immunostar, Hudson, WI, USA) or GABA (1:4.000; A2052, Sigma Aldrich) for 5 days at 4°C. Subsequently, samples were incubated with Anti-Rabbit Rhodamine Red™-X (RRX) IgG secondary antibody

(1:500, Ex/Em: 570/590 nm -Jackson laboratory, Bar Harbor, Maine, USA) for 3 days at 4°C; nonspecific binding was evaluated in subsamples incubated only in the presence of the secondary antibody (not shown) as previously described (Miglioli et al., 2021a). Larvae were then stained with 1 µg/mL Hoechst in PBS (UV, Ex/Em: 352⁄461 nm-Hoechst 33342, Invitrogen) and mounted with glycerol phosphate buffered solution containing antifading agent (CitiFluorTM AF1, pH 10, Agar scientific). Samples were imaged with a Leica Sp8 CLSM. Hoechst, 5-HT immunoreactivity and Brightfield signals were sequentially scanned with a 0.3 µm Z-stack interval. Using ImageJ software, channels were merged and subsequent Z-stacks were assembled in order to assess the number of 5-HT-ir (immuno-reactive) cells as previously described (Miglioli et al., 2021a). Larvae stained with GABA Abs were entirely imaged in 0.6 µm Z-stacks for a total of 100 z planes, with Z-projections made with 10 sequential z planes (from z10 to z100). At least 50 individuals were analyzed for each sample and experimental condition (N=5).

Statistics

Data obtained from 4 or 5 independent parental pairs, depending on the experiment (mean ± SD) were analysed by one-way nonparametric ANOVA (Kruskal-Wallis test) followed by the Tukey's test (p<0.05) using GraphPad Prism 5 software (GraphPad Inc.).

Results

Phenotypic changes induced by TBBPA in mussel larvae at 48 hpf

The phenotypic effects of TBBPA were first evaluated in 48 hpf larvae and the results are reported in Fig. 1. TBBPA induced a concentration-dependent decrease in the percentage of normal D-veligers as previously described (Fabbri et al., 2014) (Fig. 1A). At all concentrations, immature larvae (preveligers), as well as larvae withheld at the trocophora stage were observed, together with shell malformations. These effects could be better appreciated by Calcofluor/Calcein double staining, that shows deposition of each shell component, the organic matrix (blue) and CaCO3 (green) (Fig. 1B).

The shell of normal D-veligers showed extensive calcification, accretion rings and straight hinge. In TBBPA-exposed samples (10 µg/L), malformed larvae were characterized by irregular calcification and hinge malformations; immature D-veligers were identified by incomplete calcification and protrusion of the velum, the larval feeding organ; other larvae showed the arrested trocophora phenotype (Fig. 1B).

The effects of TBBPA (10 and 1 µg/L) on larval shells were investigated in more detail by SEM, and representative images are reported in Fig. 2 and3, respectively. As previously described [START_REF] Balbi | Impact of bisphenol A (BPA) on early embryo development in the marine mussel Mytilus galloprovincialis: Effects on gene transcription[END_REF][START_REF] Balbi | Impact of cationic polystyrene nanoparticles (PS-NH2) on early embryo development of Mytilus galloprovincialis: Effects on shell formation[END_REF][START_REF] Balbi | Diclofenac affects early embryo development in the marine bivalve Mytilus galloprovincialis[END_REF]) normal D-veliger shells were characterized by a straight hinge, symmetric valvae and a uniform, smooth surface (Fig. 2A, Fig. 3A). In TBBPA-exposed larvae (10 μg/L) shells showed irregular surfaces, holes, convex hinges, and fractures (Fig. 2B-E) often perpendicular to the hinge axis (Fig. 2F). At lower concentrations (1 μg/L) in addition to shell fractures (Fig. 3B), immature shells, characterized by polygonal shapes, asymmetric valvae (Fig. 3C,D) and protruding velum (Fig. 3E) were observed. Irregular calcification at the center of shell was also detected (Fig. 3F).

Effects of TBBPA on shell formation at earlier larval stages

The effects of TBBPA were further investigated by Calcofluor/Calcein staining at shorter times postfertilization (24, 28, and 32 hpf), and the results are shown in Fig. 4. In Fig. 5, quantification of areas of occupied by the organic matrix and calcified shell in different experimental conditions are reported.

In control larvae, at 24 hpf (the trocophora stage) the initial shell was only made of organic matrix (blue) and characterized by a saddle shape (the shell field); calcification (green) was detectable from 28 hpf, starting from the center of each valva, and increased at 32 hpf, with evident shell accretion rings, whereas the hinge region was not calcified yet. Complete calcification, including the hinge, was observed at 48 hpf (Fig. 4). Exposure to TBBPA (10 µg/L) induced a significant decrease in the calcified shell area from 28 hpf, and in matrix deposition at 32 hpf (Fig. 4 and Fig. 5B,C).

Interestingly, in TBBPA-exposed larvae, absence of calcification was observed in the center of the valvae from 28 to 48 hpf (Fig. 4), corresponding to the irregular shell calcification observed by SEM in some individuals at 48 hpf. Overall, TBBPA (10 µg/L) induced a significant increase in the ratio organic matrix/calcified shell at 28 hpf (Fig. 5D). Stronger effects were observed at higher concentrations (100 µg/L) (Fig. 4 and Fig. 5A-D). All the following experiments were performed at 10 µg/L TBBPA.

Effects of TBBPA on transcription of genes involved in early shell formation

Transcription of selected genes involved in first shell formation [START_REF] Miglioli | Characterization of the main steps in first shell formation in Mytilus galloprovincialis: possible role of tyrosinase[END_REF]) was evaluated by qPCR in larvae exposed to TBBPA (10 µg/L) at different times pf. In Fig. 6 data are reported on changes in mRNA level for genes involved in matrix deposition and remodeling (chitin synthase-CS; tyrosinase-TYR) and calcification (carbonic anhydrase-CA; extrapallial protein-EP) in control and TBBPA-exposed samples across different stages with respect to eggs. Transcription of CS was upregulated at 32 and 48 hpf in control samples, and TBBPA induced a further small but significant increase at 48 hpf (+20% with respect to controls) (Fig. 6A). Although basal expression of TYR in control larvae is very high with respect to eggs at all times pf (Miglioli et al., 2019, 2021a), TBBPA induced a further general upregulation at all stages, that was highest at both 28 and 48 hpf (+85% vs controls, p<0.05) (Fig. 6B). Transcription of CA was unaffected by TBBPA (Fig. 6C); however, TBBPA induced a general increase in mRNA levels for EP, with highest upregulation at 28 and 48 hpf (+105 and +135% vs controls, respectively, p<0.05) (Fig. 6D).

The expression pattern of tyrosinase and of chitinase, the enzyme involved in degradation of chitin, the main component of the organic matrix, was also investigated by ISH, and the results are reported in Fig. 7. As previously demonstrated (Miglioli et al., 2021a), tyrosinase was strongly expressed at all time pf in the areas corresponding to the organic matrix of the growing shell. TBBPA exposure resulted in a general upregulation of tyrosinase, in line with the results obtained by PCR, without significantly affecting its localization. Chitinase showed expression in the central part of the shell in larvae from 24 to 32 hpf, followed by a sharp decrease in expression in D-veligers. TBBPA did not affect chitinase localization at 24 and 28 hpf, whereas it induced some downregulation at 32 hpf and some upregulation at 48 hpf (Fig. 7). In analogy with data obtained in oysters [START_REF] Liu | D1 dopamine receptor is involved in shell formation in larvae of Pacific oyster Crassostrea gigas[END_REF][START_REF] Liu | Ocean acidification inhibits initial shell formation of oyster larvae by suppressing the biosynthesis of serotonin and dopamine[END_REF], the possible involvement of the dopaminergic system in mussel larval development was investigated. Exposure to the DR1 antagonist SCH 23390 resulted in a complete disruption of shell formation. In particular, no calcification was observed up to 32 hpf (Fig. S2) and, from 24 to 32 hpf, the area of the shell field occupied by the organic matrix was reduced to less than 20% of controls (Fig. S2B). At 48 hpf, a severe impairment of calcification was observed (Fig. S2B), the overall effects resulting in about an 80% decrease in the percentage of normal D-larvae (Fig. S2C).

TBBPA affects the expression pattern of dopaminergic components

The effects of TBBPA (10 µg/L) on key components of the dopaminergic system were then investigated. To this aim, expression pattern of key enzymes involved in DOPA synthesis (Tyrosine Hydroxylase) and metabolism (Dopamine-β-Hydroxylase) and of the DOPA receptor DR1 were evaluated by ISH, and the results are shown in Fig. 8. Exposure to TBBPA induced changes in the expression pattern of all genes from 24 hpf. In control larvae, the dopaminergic marker TH, responsible for the synthesis of DOPA from the amino acid tyrosine, showed a diffuse localization at 24 hpf, followed by progressive concentration in two clusters dopaminergic neurons at 28 hpf and 32 hpf; at 48 hpf, a clear localization of TH in a single cluster was observed. TBBPA-exposed larvae showed a stronger and more diffuse signal at earlier stages, with two observable clusters remaining at 48 hpf (suggesting a delay in the differentiation of dopaminergic cell ganglion). Expression of DβH (the enzyme that catalyzes the conversion of DOPA to noradrenaline) was detectable in control samples from 28 hpf at the center of the larval body and showed a distinct localization only at 48 hpf.

TBBPA resulted in a general upregulation of DβH at all times pf, except for at 28 hpf. DR1 in control larvae showed a weak expression in the trocophorae at 24 hpf and subsequently increased during transition to the D-veliger, showing a sharp localization along the margins of the growing shell. In TBBPA-exposed larvae, DR1 was not expressed at 24 hpf; a weaker signal and less distinct localization were observed from 28 to 32 hpf, followed by strong upregulation at 48 hpf.

TBBPA affects serotonin and GABA immunoreactivity

Development of serotonin immunopositive cells (5-HT-ir), a key step in early mussel larvae neurodevelopment, was evaluated by immunocytochemistry as previously described (Miglioli et al., 2021a), and the results are reported in Fig. 9 and Fig. 10. In control larvae, the number of 5-HT-ir cells progressively increased from 2 to 7 from 24 to 48 hpf (Fig. 9A and Fig. 10A) as previously reported (Miglioli et al. 2021a). In samples exposed to TBBPA (10 µg/L), the 5-HT positive ganglion was still present: however, the number of 5-HT-ir neurons was generally lower at all times pf (Fig. 9B and Fig. 10B). In particular, at 48 hpf, such a decrease, was clearly evident in immature veligers and arrested trocophorae (Fig. 10B), and statistically significant in all treated larvae with respect to controls (Fig. 11). However, TBBPA did not affect mRNA levels for the serotonin receptor 5-HTR at any time pf (not shown).

Development of GABA-ir cells was also evaluated from 24 to 48 hpf (Fig. 12 and Fig. 13). Due to strong GABA signal observed at all times pf, larvae were entirely imaged in 100 z-stacks, with Zprojections made every 10 sequential stacks (from z10 to z100) to better identify the localization of GABA-ir cells (Fig S3 andS4). In control larvae, several GABA-ir cells were detected at 24 and 28 hpf, along the shell field and at the center of the body (Fig. 12A and Fig S3). From 32 hpf, immunoreactivity was concentrated along the external margins of the larvae. In TBBPA-exposed samples (10 µg/L), a clear decrease in the GABA signal was observed from 24 to 32 hpf (Fig. 12B).

At 48 hpf, GABA immunoreactivity was concentrated along the external margin of the body in control samples (Fig. 13A and Fig. S4). At this stage, in samples exposed to 10 µg/L TBBPA, GABAir cells were localized in the protruding velum of immature veligers (pre-veligers), and arrested trocophores showed the strong immunoreactivity observed in control 24-28 hpf trocophores (Fig. 13B).

Discussion

The results here obtained extend previous observations indicating that TBBPA interferes with early larval development of M. galloprovincialis (Fabbri et The results here presented provide information on the mechanisms of action of TBBPA in developing larvae across different stages, focusing on shell morphogenesis and neurogenesis. For a mechanistic study, all experiments were carried out at 10 µg/L, a concentration about two-fold higher than both the EC50 obtained in the 48 h standard larval assay and TBBPA levels detected in polluted coastal Finally, TBBPA showed dramatic effects on the development of GABA-ir cells in mussel larvae at early stages (24 hpf). In mammalian cells, TBBPA and other brominated flame retardants inhibited uptake of GABA and 5-HT [START_REF] Mariussen | The effect of brominated flame retardants on neurotransmitter uptake into rat brain synaptosomes and vesicles[END_REF]. GABA is a conserved neurotransmitter in different phyla (Miller, this further emphasize the sensitivity of mussel early larval stages towards TBBPA.

Conclusions

The results demonstrate that TBBPA affected components of DOPA, 5-HT and GABA systems in mussel larvae from 24 hpf; from 28 hpf, these effects were accompanied by changes in deposition of organic matrix and calcification of the growing shell, and altered transcription pattern of genes crucial for shell formation. The overall effects of TBBPA resulted in larval malformations and delayed/arrested development at 48 hpf. Although for a mechanistic study all the endpoints were evaluated at a concentration two-fold higher than those detected in polluted coastal areas, to our knowledge, these results represent the demonstration of the effects of TBBPA at the lowest concentrations tested so far in bivalves. Data, reported as fold changes with respect to eggs, represent the mean ± SD (N=4). Data were analysed by one way ANOVA Kruskall-Wallis followed by Tukey's test and statistical differences between control and TBBPA-exposed samples at each developmental stage are reported (* = p<0.01).

As previously described [START_REF] Miglioli | Characterization of the main steps in first shell formation in Mytilus galloprovincialis: possible role of tyrosinase[END_REF], expression of all transcripts was significantly different at all times pf with respect to eggs (p<0.01, not shown). b) GABA and Hoechst merged channels; the z-stack acquisition is indicated in white (see Fig. S3).

Red arrowheads indicate the shell field region, and the perimeter of the larvae is indicated by a white dotted line; scale bar: 20 µm.

3. 4

 4 Role of thyroid and dopamine signaling in mussel shell morphogenesis as possible targets of TBBPA Since in vertebrate developmental models, including the zebrafish, TBBPA mainly affects thyroid hormone signalling (Zhu et al., 2018 and refs. therein; Yu et al., 2021) the role of this pathway was first investigated by exposure of fertilized eggs to the thyroid Receptor-THR agonists Triiodotyronine (T3) and Tetraiodotyronine (T4), the THR antagonist 1-850, and the antithyroid drug thiourea (TU).None of these compounds significantly affected larval development at 48 hpf, in terms of changes in proportion of larval phenotypes (Fig.S1). Only the THR antagonist 1-850 induced a slight decrease in the percentage of normal D-veligers; however, the values were above the acceptability of the test results (>75% normal larvae) (ASTM, 2012). The results indicate that thyroid hormone signalling does not play a key role in shell morphogenesis in mussel early larval stages.

  areas. The utilization of calcein/calcofluor staining to visualize and quantify the deposition of each shell component, clearly demonstrated that TBBPA significantly increased the ratio organic matrix/calcified shell at 28 hpf, when the first calcification of the growing shell occurs[START_REF] Miglioli | Characterization of the main steps in first shell formation in Mytilus galloprovincialis: possible role of tyrosinase[END_REF], indicating a perturbation of the firsts steps of shell biogenesis. Moreover, TBBPA affected the pattern of calcification; in particular, no CaCO3 deposition at the center of the valvae was observed, resulting in a 'key hole' phenotype similar to that previously observed in mussel larvae exposed to acidification[START_REF] Kapsenberg | Ocean pH fluctuations affect mussel larvae at key developmental transitions[END_REF]. These data indicate that TBBPA affects the mechanisms involved in deposition of the organic matrix, which represents the blueprint for biomineralization and normal shell biogenesis. The effects of TBBPA on shell formation could be appreciated by SEM imaging at 48 hpf, at concentrations as low as 1 µg/L, confirming delayed development and irregular calcification. Data on gene expression offer a first glimpse into the underlying mechanisms of TBBPA on shell formation, indicating a general dysregulation of the expression pattern of genes involved in matrix deposition and calcification across different developmental stages, with tyrosinase representing the most upregulated gene from 28 hpf. Tyrosinase has been identified as a key enzyme in organic matrix remodeling in early mussel larvae grown in physiological conditions[START_REF] Miglioli | Characterization of the main steps in first shell formation in Mytilus galloprovincialis: possible role of tyrosinase[END_REF], as well as a target of BPA(Miglioli et al., 2021a). The results here obtained indicate that TBBPA induced a stronger upregulation of tyrosinase with respect to BPA at similar concentrations (10 vs 11.4 µg/L) at all stages pf(Miglioli et al., 2021a); the effects were also stronger when compared on a molar basis (18 nM TBBPA vs 50 nM BPA). TBBPA also resulted in upregulation of EP, that plays a key role in calcification, with a similar trend across different stages. Only at 48 hpf, increased expression of chitin synthase and decreased that of chitinase were observed, as shown by both qPCR and ISH.These results indicate that increased synthesis of chitin, the main component of the shell organic matrix, typical of the initial stages of shell deposition, is upregulated by TBBPA at later stages, further supporting a delay in development. Overall, these data partly contribute to explain the phenotypic effects of TBBPA on shell formation.BPA has been previously shown to affect expression of genes involved in detoxification response in mussel larvae[START_REF] Balbi | Impact of bisphenol A (BPA) on early embryo development in the marine mussel Mytilus galloprovincialis: Effects on gene transcription[END_REF]. However, TBBPA (up to 10 µg/L) did not induce changes in mRNA levels for glutathione-s-transferase (GSTπ) and the ABC transporter p-glycoprotein (ABCB) with respect to controls at any developmental stage (not shown), suggesting that mussel early larvae have little capacity for TBBPA biotransformation/detoxification. The possible endocrine mechanisms of action of TBBPA were then investigated. In vertebrate systems TBBPA, structurally resembling thyroid hormones, can act as an endocrine disruptor affecting TH signaling and metabolism (Zhang et al., 2014; Zhou et al., 2020). These mechanisms have been shown to be associated with the neurodevelopmental toxicity of TBBPA in the zebrafish embryo (Zhu et al., 2018; Pang et al., 2020; Yu et al., 2021). The presence of THs and related signaling components has been demonstrated in developing and adult bivalves (Huang et al., 2015; Taylor and Heyland, 2017; Miglioli et al., 2021b). In clam juveniles, prolonged exposure to TBBPA caused alterations in tissue content of THs and in expression of TH synthesis-related enzymes (Jiang et al., 2019). In the oyster C. gigas, TH signaling was shown to be involved in embryogenesis and metamorphosis; however, from the trocophora to the first D-veliger stage, lowest values of T3 and T4 and absence of THR expression were observed (Huang et al., 2015). Accordingly, the results here obtained with M. galloprovincialis show that exposure to different TH agonists and antagonists did not affect larval development at 48 hpf. These observations support the hypothesis that TH signaling, which is considered a common target for the endocrine effects of TBBPA also in aquatic species, does not represent a key pathway in the transition from the trocophora and the first D-veliger of bivalve larvae. Therefore, other mechanisms may be involved in the action of TBBPA as a potential developmental disruptor in early developmental stages of mussel larvae. Neuroendocrine mechanisms are known to play a vital role in larval development of invertebrates (Cann-Moisan et al., 2002). Molecular components of the neuroendocrine system have been characterized in marine bivalves (catecholaminergic, cholinergic, enkephalinergic, serotoninergic, gamma-aminobutyric acid-ergic, and neuropeptide systems), and they share a similar molecular basis with their vertebrate counterparts (Wang et al., 2018; Liu et al., 2018; Joyce and Vogeler, 2018). In larvae of the oyster C. gigas, expression of a DOPA receptor Cg-DR-1 peaked during the transition from the trocophora to the D-veliger stage, indicating a role in shell biogenesis[START_REF] Liu | D1 dopamine receptor is involved in shell formation in larvae of Pacific oyster Crassostrea gigas[END_REF]. In a subsequent study, exposure of trocophorae to DR1 and 5-HTR inhibitors was found to prevent shell formation[START_REF] Liu | Ocean acidification inhibits initial shell formation of oyster larvae by suppressing the biosynthesis of serotonin and dopamine[END_REF]. Metabolomic and transcriptomic analyses suggested that the synthesis of DOPA and 5-HT may modulate initial shell formation through the TGF-β smad pathway, by activating the expression of tyrosinase and chitinase, that are essential for initial shell biogenesis[START_REF] Liu | Ocean acidification inhibits initial shell formation of oyster larvae by suppressing the biosynthesis of serotonin and dopamine[END_REF].The results here presented extend the findings in oysters by showing that in M. galloprovincialis, exposure of fertilized eggs to the selective DR1 inhibitor SCH 23390 disrupted shell formation at a concentration 20 times lower than those utilized in oyster trocophorae[START_REF] Liu | Ocean acidification inhibits initial shell formation of oyster larvae by suppressing the biosynthesis of serotonin and dopamine[END_REF]. Therefore, DOPA signaling has a key role in mussel shell morphogenesis. TBBPA induced alterations in the expression pattern of enzymes involved in DOPA synthesis and metabolism and of the DOPA receptor DR1 from 24 hpf, indicating a delay in the differentiation of dopaminergic neurons and DR1 expression, as well as an increased metabolism of DOPA. The effects on the dopaminergic system were detectable before the changes in expression of shell related genes observed from 28 hpf. We hypothesize that impairment of dopamine system interferes with the initial processes of shell deposition, resulting in altered shell phenotypes and delayed development. Only at later stages (48 hpf), TBBPA-induced upregulation of DR1 was accompanied by increased expression of tyrosinase, chitin synthase, EP, whereas chitinase was downregulated. Another link between dopaminergic signaling and shell formation is suggested by the results obtained on expression of mussel EP protein, that is involved in controlling the deposition of CaCO3 polymorphs, including vaterite, the most unstable crystalline phase formed during initial biomineralization[START_REF] Ji | Control of crystal polymorph in microfluidics using molluscan 28 kDa Ca 2+ -binding protein[END_REF]. In the pearl mussel Hyriopsis cumingii DOPA has been involved in stability of vaterite[START_REF] Kim | Dopamine-induced mineralization of calcium carbonate vaterite microspheres[END_REF].Serotonin represents a key navigational cue during bivalve neurodevelopment[START_REF] Yurchenko | Peripheral sensory neurons govern development of the nervous system in bivalve larvae[END_REF] and shell formation[START_REF] Liu | Ocean acidification inhibits initial shell formation of oyster larvae by suppressing the biosynthesis of serotonin and dopamine[END_REF]. In M. galloprovincialis, BPA affected expression of tyrosinase, larval patterning and 5-HT signalling(Miglioli et al., 2021a). The results here presented show that also exposure to TBBPA impairs the development of serotoninergic neurons from 24 hpf, whose number was decreased in 48 hpf larvae, characterized by delayed development and arrested phenotypes. Overall, the results indicate that tyrosinase, that represents a link between monoamine neurotransmitter signaling and shell biogenesis in both oysters and mussels[START_REF] Liu | Ocean acidification inhibits initial shell formation of oyster larvae by suppressing the biosynthesis of serotonin and dopamine[END_REF][START_REF] Miglioli | Characterization of the main steps in first shell formation in Mytilus galloprovincialis: possible role of tyrosinase[END_REF] Miglioli et al., , 2021a)), is a potential target for endocrine disruptors in early bivalve development.

Figure 1 -

 1 Figure 1 -Effect of TBBPA (1, 10 and 100 µg/L) on larval phenotypes at 48 hpf.

Figure 2 -

 2 Figure 2 -Effects of TBBPA (10 μg/L) on shell morphology of Mytilus larvae at 48 hpf, evaluated

Figure 3 -

 3 Figure 3 -Effects of TBBPA (1 μg/L) on shell morphology of Mytilus larvae at 48 hpf, evaluated

Figure 4 -

 4 Figure 4 -Effect of TBBPA (10 and 100 µg/L) on the time course of early shell formation in M.

Figure 5 -

 5 Figure 5 -Effect of TBBPA (10 and 100 µg/L) on the growth of organic matrix (blue) and

Figure 6 -

 6 Figure 6 -Expression of genes related to shell formation in control and TBBPA-exposed larvae

Figure 7 -

 7 Figure 7 -Expression pattern of Tyrosinase and Chitinase in control and TBBPA-exposed

Figure 8 -

 8 Figure 8 -Expression pattern of dopaminergic components control and TBBPA-exposed larvae

Figure 9 -M. galloprovincialis from 24 to 32 hpf. 5 -

 95 Figure 9 -Confocal images of serotonin immunoreactive cells (5-HT-ir) in early larval stages of

Figure 10 -

 10 Figure 10 -Confocal images of serotonin immunoreactive cells (5-HT-ir) in early larval stages

Figure 11 -

 11 Figure 11 -Effect of TBBPA (10 µg/L) on development of serotonin immunoreactive cells (5-

Figure 12 -

 12 Figure 12 -Confocal images of GABA immunoreactive cells (GABA-ir) in early larval stages of

Figure 13 -

 13 Figure 13 -Confocal images of GABA immunoreactive cells (GABA-ir) in early larval stages of

Figure 1

 1 Figure 1

Figure 2

 2 Figure 2

Figure 3

 3 Figure 3

Figure 4

 4 Figure 4

Figure 5

 5 Figure 5

Figure 6

 6 Figure 6

Figure 7

 7 Figure 7

Figure 8

 8 Figure 8

Figure 9

 9 Figure 9

Figure 10 Click

 10 Figure 10

Figure 11 Click

 11 Figure 11

Figure 12

 12 Figure 12

Figure 13 Click

 13 Figure 13

  Springborn Life Sciences, 1989). In Mytilus edulis, chronic exposure (70 days) to TBBPA affected growth, with NOEC and LOEC values of 17 and 32 µg/L, respectively (AstraZeneca, 2005a,b). In juveniles of the clam Ruditapes phillipinarum exposed to TBBPA (62.5 -1000 µg/L, 28 days) inhibition of shell growth and filtration rate were observed[START_REF] Jiang | Inhibition of growth in juvenile manila clam Ruditapes philippinarum: Potential adverse outcome pathway of TBBPA[END_REF].

	Overall, these data indicates that shell growth represents a target for TBBPA in adult and juvenile
	bivalves, although at higher concentrations than those detected in the environment.
	In M. galloprovincialis, TBBPA in a wide concentration range, encompassing environmental
	exposure levels (from ng to µg/L), affected early larval development, increasing the percentage of
	TBBPA has been reported as a potential thyroid endocrine disruptor and neurotoxicant in mammalian
	and fish (i.e. the zebrafish) models through multiple mechanisms, indicating that particular attention
	should be paid to the effects of this compound on early developmental stages (reviewed in Zhou et
	al., 2020). In marine organisms, the effects of TBBPA have been investigated in a range of species

al., 2014) at concentrations encompassing those detectable in contaminated coastal areas (Gong et al., 2017 and refs. therein). Moreover, we provide first information on the mechanisms of action of TBBPA on mussel shell biogenesis and neurodevelopment from earlier larval stages. The effects of TBBPA on all endpoints measured from molecular to individual level were observed at 10 µg/L, a concentration of the same order of magnitude as those detected in polluted environments (low µg/L) (Liu et al., 2016; Gong et al., 2017 and refs. therein; Li et al., 2020). across different trophic levels, indicating that some fish and molluscs are particularly sensitive to TBBPA (Pittinger and Pecquet, 2018; Nos et al., 2020). With regards to marine bivalves, acute exposure of the Eastern oyster Crassostrea virginica to TBBPA (18 µg/L, 96 h) significantly reduced shell growth (trocophorae/immature D-veligers

[START_REF] Fabbri | Adaptation of the bivalve embryotoxicity assay for the high throughput screening of emerging contaminants in Mytilus galloprovincialis[END_REF]

. However, the possible mechanisms of action were not investigated, due to little information available on the molecular pathways regulating shell biogenesis in mussels. Subsequent studies shed some light on the expression pattern of genes involved in shell formation, and in neuroendocrine pathways in early mussel larvae grown both in physiological conditions and in the presence of different emerging contaminants

[START_REF] Balbi | Impact of bisphenol A (BPA) on early embryo development in the marine mussel Mytilus galloprovincialis: Effects on gene transcription[END_REF][START_REF] Balbi | Impact of cationic polystyrene nanoparticles (PS-NH2) on early embryo development of Mytilus galloprovincialis: Effects on shell formation[END_REF][START_REF] Balbi | Diclofenac affects early embryo development in the marine bivalve Mytilus galloprovincialis[END_REF][START_REF] Franzellitti | Phenotypical and molecular changes induced by carbamazepine and propranolol on larval stages of Mytilus galloprovincialis[END_REF][START_REF] Miglioli | Characterization of the main steps in first shell formation in Mytilus galloprovincialis: possible role of tyrosinase[END_REF] Miglioli et al., , 2021a))

.

  2019;[START_REF] Smart | A half century of γ-aminobutyric acid[END_REF]: in molluscs, GABA has been shown to affect larval settlement and metamorphosis, with possible downstream effects on serotonin pathways[START_REF] Biscocho | GABA is an inhibitory neurotransmitter in the neural circuit regulating metamorphosis in a marine snail[END_REF] Joyce and Vogeler, 2018, and refs therein). The high abundance of GABA-ir cells observed in control mussel trocophorae (at 24 hpf) is in line with the role of this neurotransmitter in the early developing mammalian brain, where the GABAergic system plays an important trophic role, even prior to synaptogenesis, in the proliferation of neuronal progenitor cells (reviewed in[START_REF] Wu | GABA receptors in brain development, function, and injury[END_REF]. Overall, on the basis of our first observations, the dramatic decrease in GABA immunoreactivity induced by TBBPA further suggests that, in mussel larvae, TBBPA may act upstream of the first stages of neurodevelopment, to affect subsequent differentiation of serotoninergic or other neurons. Interference with GABA signalling has been involved in TBBPA neurotoxicity in Xenopus and zebrafish embryos[START_REF] Hendriks | Multiple novel modes of action involved in the in vitro neurotoxic effects of tetrabromobisphenol-A[END_REF][START_REF] Ye | Metabolomics approach reveals metabolic disorders and potential biomarkers associated with the developmental toxicity of tetrabromobisphenol A and tetrachlorobisphenol A[END_REF]. However, irrespective of the mechanisms of action, in vertebrate developmental models significant effects of TBBPA were observed from concentrations generally higher than those utilized in the present study;

  The signaling pathways linking development of neurotransmitters with shell formation still need to be elucidated; however, the data obtained so far (this work,Miglioli et al., 2021a) indicate that EDs such as TBBPA and BPA can affect the complex processes involved in early neurodevelopment and shell biogenesis in the model marine bivalve Mytilus.The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. 54. Zhu, B., Zhao, G., Yang, L., Zhou, B., 2018. Tetrabromobisphenol A caused neurodevelopmental toxicity via disrupting thyroid hormones in zebrafish larvae.

	Chemosphere 197, 353-361. https://doi.org/10.1016/j.chemosphere.2018.01.080
	55. Yu, Y., Hou, Y., Dang, Y., Zhu, X., Li, Z., Chen, H., Xiang, M., Li, Z., Hu, G., 2021.
	The identification of the mode of action (MOA) of EDs in marine invertebrates can contribute Exposure of adult zebrafish (Danio rerio) to Tetrabromobisphenol A causes neurotoxicity
	building Adverse outcome pathways (AOPs) (Ankley et al., 2010; Ford and Le Blanc, 2020), in larval offspring, an adverse transgenerational effect. J. Hazard Mater. 414, 125408.
	encompassing relationships from molecular initiating events (MIE) caused by the initial exposure, to https://doi.org/10.1016/j.jhazmat.2021.125408
	a series of key events (KE), that lead to identification of effects on individuals (i.e. developmental
	effects) and potential loss of population sustainability. Although the possible consequences of
	TBBPA at the population level cannot at present be demonstrated, the results of this work provide
	data on possible progression of the effects of this compound from the molecular to the individual
	level.
	A better understanding of the neuroendocrine function during development in model marine
	invertebrates will greatly help understanding the impact of emerging pollutants, in particular EDs, as
	well as of other environmental stressors, on more susceptible life stages of ecologically important
	species in coastal environments.
	CRediT authorship contribution statement
	Angelica Miglioli: Investigation, Methodology, Data curation, Writing. Teresa Balbi: Investigation,
	Data curation, methodology, Writing, review & editing, funding acquisition. Michele Montagna: data
	curation, methodology. Remi Dumollard: Conceptualization, Supervision, Writing, funding
	acquisition. Laura Canesi: Conceptualization, Supervision, Writing.
	Declaration of Competing Interest

Acknowledgments

The authors would like to thank Laurent Gilletta, Alexandre Jean, Régis Lasbleiz and the Centre de Ressources Biologiques Marines of the CRBM-IMEV that is supported by EMBRC-France, whose French state funds are managed by the Agence Nationale de la Recherche (ANR) within the 'Investissement d'Avenir' program (ANR-10-INBS-02); Sameh Benaicha; the Ascidian BioCell group members. This work was partially supported by an ANR grant (Marine-EmbryoTox project, ANR-14-OHRI-0009-01-1, R. Dumollard), ad by Fondi di Ricerca di Ateneo, FRA2020, University of Genoa (T. Balbi).

Supplementary material for on-line publication only

Click here to access/download Supplementary material for on-line publication only Miglioli et

Effect of Thyroid receptor agonists/antagonists on early larval development

Fertilized eggs were exposed to thyroid hormones [3,3,5,5-triiodothyronine-T4 (T2376, Sigma Aldrich) and 3,3,5-triiodothyronine-T3 (T6397, Sigma Aldrich)], the selective antagonist of human

Results

Figure S1 -Effect of THRs agonists and antagonists on larval development at 48 hpf.

Effects of T3 (3,3,5-triiodothyronine), T4 (3,3,5,5-triiodothyronine), 1-850 the selective antagonist of human THR, and Thiourea-TU on larval phenotypes: normally developed D-Veligers (white), with hinge malformations (light grey), pre-veligers (dark grey) and arrested trocophorae stage (black).

Measurements were made on at least 50 larvae for each experimental condition obtained from 5 parental pairs (N=5).