Laurène Mérillet

Sandrine Pavoine

Dorothée Kopp

Marianne Robert

Maud Mouchet

Biomass of slow life history species increases as local bottom trawl effort decreases in the Celtic sea

Keywords: Community weighted mean, Assembly processes, Biological traits, Elasmobranchs, Area closure

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Marine ecosystems experience multiple pressures, the most severe of which include climate change and fishing [START_REF] Halpern | Spatial and temporal changes in cumulative human impacts on the world's ocean[END_REF]. Climate change already affects the structure, functioning and adaptive capacity of ecosystems by modifying pH, temperature, oxygen concentration and food availability [START_REF] Henson | Rapid emergence of climate change in environmental drivers of marine ecosystems[END_REF], which notably results in displacement of species toward higher latitudes and deeper waters [START_REF] Cheung | Signature of ocean warming in global fisheries catch[END_REF][START_REF] Jorda | Ocean warming compresses the three-dimensional habitat of marine life[END_REF]. Its effects occur along with those of fishing [START_REF] Rogers | Shifting habitats expose fishing communities to risk under climate change[END_REF], which has caused a decrease in biomass and even local extinctions of species worldwide [START_REF] Worm | Impacts of Biodiversity Loss on Ocean Ecosystem Services[END_REF]. It is widely acknowledged that a population's ability to recover from disturbance depends on its demographic dynamics, which are related to life history strategies [START_REF] Le Quesne | Predicting species vulnerability with minimal data to support rapid risk assessment of fishing impacts on biodiversity[END_REF][START_REF] Wiedmann | Life history variation in Barents Sea fish: Implications for sensitivity to fishing in a changing environment[END_REF].

Historically, life history strategies were characterised by their position along the continuum of r strategies (energy allocated to producing many offspring) vs. K strategies (energy allocated to producing a few extremely fit offspring) [START_REF] Pianka | On r-and K-Selection[END_REF]. More recently, life history strategies of fish were described as a triangle: (i) an opportunistic strategy for small, short-lived and rapidly maturing fish; (ii) a periodic strategy for large, long-lived and highly fecund fish and (iii) an equilibrium strategy for intermediate-sized fish that produce a few large offspring for which they provide parental care [START_REF] Pecuchet | Spatial distribution of life-history traits and their response to environmental gradients across multiple marine taxa[END_REF][START_REF] Winemiller | Patterns of life-history diversification in North American fishes: implications for population regulation[END_REF]. These life history strategies can be summarised along a fast-slow continuum (Beukhof et al., 2019b;[START_REF] Juan-Jorda | Life in 3-D : life history strategies in tunas , mackerels and bonitos[END_REF][START_REF] Promislow | Living fast and dying young : A comparative analysis of life-history variation among mammals[END_REF][START_REF] Rochet | Comparative analysis of phylogenetic and fishing effects in life history patterns of teleost fishes[END_REF] that balances percentage of mortality and optimal size. Large and long-lived species that mature late and give birth to large offspring lie at the "slow" end of the continuum [START_REF] Wiedmann | Life history variation in Barents Sea fish: Implications for sensitivity to fishing in a changing environment[END_REF]. These species escape predation and have low natural mortality rate, while the fast end is characterised by the opposite properties. Species with a slow life history strategy are thus particularly sensitive to additional mortality, such as that caused by anthropogenic pressures (e.g. fishing) [START_REF] Brown | Toward a metabolic theory of ecology[END_REF][START_REF] Juan-Jorda | Population declines of tuna and relatives depend on their speed of life[END_REF][START_REF] Kozlowski | Why life history are diverse. State-of-the-art review[END_REF][START_REF] Promislow | Living fast and dying young : A comparative analysis of life-history variation among mammals[END_REF].

Understanding the spatio-temporal distribution of these slow life history species is thus particularly useful for ecosystem management [START_REF] Le Quesne | Predicting species vulnerability with minimal data to support rapid risk assessment of fishing impacts on biodiversity[END_REF].

Identifying and explaining how species are distributed has long been a core challenge in ecology. Evolutionary history, environmental variables and species interactions drive the spatio-temporal distribution of species [START_REF] Mouchet | Invariant scaling relationship between functional dissimilarity and co-occurrence in fish assemblages of the Patos Lagoon estuary (Brazil): Environmental filtering consistently overshadows competitive exclusion[END_REF][START_REF] Webb | Phylogenies and Community Ecology[END_REF]. However, the distribution of the taxa could also arise from random processes and the question of the relative importance of deterministic versus stochastic assembly rules remains central [START_REF] Vellend | Assessing the relative importance of neutral stochasticity in ecological communities[END_REF]. Three major nonexclusive assembly rules explain biodiversity patterns [START_REF] Kraft | Trait evolution, community assembly, and the phylogenetic structure of ecological communities[END_REF]: stochasticity [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography (MPB-32)[END_REF] and deterministic processes, namely environmental filtering [START_REF] Keddy | Assembly and response rules: two goals for predictive community ecology[END_REF][START_REF] Zobel | The relative role of species pools in determining plant species richness: An alternative explanation of species coexistence?[END_REF] and competitive exclusion [START_REF] Hardin | The Competitive Exclusion Principle[END_REF].

Environmental filtering implies strong abiotic control, which results in the survival of species that have a narrow range of traits that enable them to endure environmental pressures. The principle of competitive exclusion assumes that species can coexist if they have different niche-related biological traits [START_REF] Hardin | The Competitive Exclusion Principle[END_REF]. Thus, traits of species in an assemblage may be similar if they are selected via environmental filtering or different if they are selected via competitive exclusion. A careful interpretation is essential and should be based on the biological understanding of the processes at play. Indeed, competition (and more generally biotic interactions) can lead to the same patterns as those driven by environmental filtering [START_REF] Cadotte | Should Environmental Filtering be Abandoned ?[END_REF][START_REF] Kraft | Community assembly, coexistence and the environmental filtering metaphor[END_REF] and community patterns result in both processes interacting dynamically [START_REF] Callaway | Positive interactions among alpine plants increase with stress[END_REF][START_REF] Kraft | Community assembly, coexistence and the environmental filtering metaphor[END_REF]. Conversely to these two deterministic assembly rules, neutral theory hypothesises that species coexist regardless of their biological traits due to individual demographic stochastic events of dispersal, birth and death [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography (MPB-32)[END_REF][START_REF] Vellend | Assessing the relative importance of neutral stochasticity in ecological communities[END_REF]. The relative dominance of stochastic and deterministic processes can be disentangled notably by the comparison of the observed pattern with a null model [START_REF] Vellend | Assessing the relative importance of neutral stochasticity in ecological communities[END_REF]. The prevalence of one assembly rule over the others could depends on the scale considered. At the local scale, stochastic, biotic and abiotic parameters can act on species simultaneously, and determining their relative importance remains challenging [START_REF] Mouillot | Limiting similarity, niche filtering and functional diversity in coastal lagoon fish communities[END_REF]. At larger scale, environmental filtering was considered as having a larger effect than the other processes [START_REF] Freschet | Global to community scale differences in the prevalence of convergent over divergent leaf trait distributions in plant assemblages[END_REF].

Assembly processes in exploited marine ecosystems have attracted attention only recently [START_REF] Dencker | Temporal and spatial differences between taxonomic and trait biodiversity in a large marine ecosystem : Causes and consequences[END_REF][START_REF] Ford | Functional traits reveal the presence and nature of multiple processes in the assembly of marine fish communities[END_REF][START_REF] Pecuchet | Patterns and drivers of fish community assembly in a large marine ecosystem[END_REF], as studies have historically focussed on freshwater and estuarine communities [START_REF] Mouillot | Limiting similarity, niche filtering and functional diversity in coastal lagoon fish communities[END_REF][START_REF] Peres-Neto | Patterns in the co-occurrence of fish species in streams: The role of site suitability, morphology and phylogeny versus species interactions[END_REF][START_REF] Schmera | Habitat filtering determines spatial variation of macroinvertebrate community traits in northern headwater streams[END_REF]. Nevertheless the understanding of the processes responsible for the structuration of exploited marine communities can provide valuable inputs for their sustainable management [START_REF] Dencker | Temporal and spatial differences between taxonomic and trait biodiversity in a large marine ecosystem : Causes and consequences[END_REF]. The North East Atlantic fishing grounds have shown signs of overexploitation since the 1970s [START_REF] Gascuel | Fishing impact and environmental status in European seas: A diagnosis from stock assessments and ecosystem indicators[END_REF], which caused species distribution to be driven by fishing in addition to environmental and habitat characteristics, such as temperature, depth and sediment [START_REF] Foveau | Process-driven and biological characterisation and mapping of seabed habitats sensitive to trawling[END_REF]. We hypothesise that gradual enforcement of European Union (EU) fishing policies over the past 40 years that has decreased fishing pressure in the North East Atlantic [START_REF] Fernandes | Reversal of fish stock decline in the northeast atlantic[END_REF] may have helped the biomass of species with slow life history strategies to increase, reducing the strength of the environmental filter imposed by fishing on taxa.

However, temporal patterns of life history strategies of commercial and noncommercial taxa are rarely studied as a function of environmental or anthropogenic variables, and should be taken into account in management measures more frequently [START_REF] Fromentin | Fishing effects and life history traits: A case study comparing tropical versus temperate tunas[END_REF][START_REF] Matson | Resolving associative patterns in life history parameters among marine fish stocks in the Northeast Pacific Ocean[END_REF].

The Celtic Sea has been extensively fished for decades with a peak of the number of species significantly exploited in the 1990's [START_REF] Gascuel | Fishing impact and environmental status in European seas: A diagnosis from stock assessments and ecosystem indicators[END_REF][START_REF] Guénette | Shifting baselines in European fisheries: The case of the Celtic Sea and Bay of Biscay[END_REF]. Fishing thus likely imposed a major filtering on the distribution of the taxa. Fishing mortality has then decreased at the ecosystem scale since 2010 [START_REF] Moullec | Trophic models: What do we learn about Celtic Sea and Bay of Biscay ecosystems?[END_REF] with the implementation of more restrictive TAC and the creation of a closure area in the North of the Celtic Sea [START_REF] Bergen | ICES Advices 2007, EU request on Trevose closure[END_REF].

Environmental variables in the Celtic Sea were relatively stable over the recent period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016). Since a clear impact of climate change, such as an increase in bottom temperature, is not visible yet [START_REF] Mérillet | Environment outweighs the effects of fishing in regulating demersal community structure in an exploited marine ecosystem[END_REF], the decrease in fishing effort most likely had the largest influence on the recent distribution of life history strategies in the Celtic Sea. Based on a 17-year time series, we explore the influence of fishing on the spatio-temporal distribution of life history strategies in the benthodemersal community by considering 101 commercial and non-commercial taxa. We assessed (i) which assembly processes and (ii) which environmental and anthropogenic variables drive the spatio-temporal distribution of sensitivity to fishing.

Finally, (iii) we focussed on the trend of the biomass of sensitive taxa and bottom trawl fishing at the local scale.

Material and methods

Taxa biomass and life history traits

For the 2000-2016 period, data on species biomass were extracted from the EVHOE (Evaluation des ressources Halieutiques de l'Ouest de l'Europe) surveys, an International Bottom Trawl Survey (ICES, 2015) performed every November in the Celtic Sea. These surveys provide reliable and consistent data for evaluating the impact of fishing as they are independent from fishing effort and follow the same standardised protocol each year. EVHOE surveys use a 36/47 GOV (Grande Ouverture Verticale) net fitted with a 20 mm cod end mesh, with a horizontal opening of ca. 20 m and a vertical opening of 4 m. Sampling hauls are performed at 4 knots for 30 min and follow a random stratified design based on sediment-depth strata within a bank of possible hauls. They are not at the same location from one year to another, but due to random sampling a same haul can be visited several times across the time series. From 2000-2016, the surveys performed 53-84 valid hauls per year, for a total of 1175 hauls. Due to variations in the accuracy of species identification, some species were grouped into higher taxonomic levels (hereafter, 'taxon'; genus, for species of Alloteuthis, Argentina, Arnoglossus, Gaidrosparus, Loligo, Munida, Mustelus, Octopus, Pomatoschistus and Sepiola; family, for Ammodytidae and Gobiidae). To identify temporal changes in taxa biomass, only those taxa that occurred at more than 5% of the sampling sites in at least one EVHOE depth strata over the time series were kept for analysis and considered correctly sampled. This process identified a total of 101 taxa from six classes (Actinopterygii, Bivalvia, Cephalopoda, Elasmobranchii, Holocephali and Malacostraca). Biomass was standardised by the area sampled during each haul and then log-transformed to reduce variability. Sampling sites were located at depths of 57-340 m (Fig. 1), mainly on sand, mud, mixed sediments, coarse sediments and boulders (Fig S1). We selected the following biological traits for their known response to fishing pressure: maximum length, longevity, age at maturity, offspring size, reproductive guild and fecundity [START_REF] Barlow | Patterns of parental investment, dispersal and size among coralreef fishes[END_REF][START_REF] Brown | Toward a metabolic theory of ecology[END_REF][START_REF] Jennings | The Effects of Fishing on Marine Ecosystems[END_REF][START_REF] Jeschke | The roles of body size and phylogeny in fast and slow life histories[END_REF][START_REF] Le Quesne | Predicting species vulnerability with minimal data to support rapid risk assessment of fishing impacts on biodiversity[END_REF][START_REF] Pecuchet | Patterns and drivers of fish community assembly in a large marine ecosystem[END_REF][START_REF] Shephard | Spatial Heterogeneity in Fishing Creates de facto Refugia for Endangered Celtic Sea Elasmobranchs[END_REF][START_REF] Tillin | Chronic bottom trawling alters the functional composition of benthic invertebrate communities on a seabasin scale[END_REF][START_REF] Tyler-Walters | A method to assess the sensitivity of sedimentary communities to fishing activities[END_REF][START_REF] Wiedmann | Life history variation in Barents Sea fish: Implications for sensitivity to fishing in a changing environment[END_REF] (Table 1). Values of these traits were extracted from the literature and the PANGAEA data library (Beukhof et al., 2019a) (see Table S1 for more details).

Table 1: Life history traits used and their relevance, when used all together, to characterise life history strategy and sensitivity to fishing.

Traits Description Relevance

Maximum length in cm Relates to many characteristics of an organism [START_REF] Brown | Toward a metabolic theory of ecology[END_REF]. Fishing is known to have led to the decrease in abundance of large body size taxa [START_REF] Shephard | Spatial Heterogeneity in Fishing Creates de facto Refugia for Endangered Celtic Sea Elasmobranchs[END_REF]) Longevity in year (maximum reported age)

Fishing is known to have led to the decrease in abundance of long-lived taxa [START_REF] Shephard | Spatial Heterogeneity in Fishing Creates de facto Refugia for Endangered Celtic Sea Elasmobranchs[END_REF][START_REF] Tillin | Chronic bottom trawling alters the functional composition of benthic invertebrate communities on a seabasin scale[END_REF]) Offspring size in mm; Diameter of the released eggs, length of the egg case or length of the young for the bearer taxa Informs on the parental investment and the energy needed for the production of one offspring [START_REF] Winkler | Offspring Size and Number : A Life History Model Linking Effort Per Offspring and Total Effort[END_REF]. Perturbations, and thus fishing, select for taxa with small offspring [START_REF] Jeschke | The roles of body size and phylogeny in fast and slow life histories[END_REF][START_REF] Wiedmann | Life history variation in Barents Sea fish: Implications for sensitivity to fishing in a changing environment[END_REF]. Linked to the resilience of the taxa after fishing event.

Age at maturity in year

Informs on generation time [START_REF] Pecuchet | Patterns and drivers of fish community assembly in a large marine ecosystem[END_REF][START_REF] Tillin | Chronic bottom trawling alters the functional composition of benthic invertebrate communities on a seabasin scale[END_REF]. Fishing favors early maturing taxa [START_REF] Claireaux | Evolutionary effects of fishing gear on foraging behavior and life-history traits[END_REF][START_REF] Shephard | Spatial Heterogeneity in Fishing Creates de facto Refugia for Endangered Celtic Sea Elasmobranchs[END_REF]. Linked to the resilience of the taxa after fishing event.

Reproductive guild

Bearer -giving birth to a free living progeny.

Informs on the parental investment and the energy needed for the production of one offspring [START_REF] Barlow | Patterns of parental investment, dispersal and size among coralreef fishes[END_REF][START_REF] Tillin | Chronic bottom trawling alters the functional composition of benthic invertebrate communities on a seabasin scale[END_REF]. Linked to the resilience of the taxa after fishing event.

Guarder -eggs protected (carried under the abdomen for crustaceans) or guarded. Non guarder -pelagic eggs.

Fecundity

Segment of offspring number per female 1-1000: low 1000-100,000: medium 100,000-1,000,000: high >1,000,000: very high Capacity of the population to recover; productivity of the taxa [START_REF] Pecuchet | Patterns and drivers of fish community assembly in a large marine ecosystem[END_REF].

Fishing depleted low fecundity taxa [START_REF] Shephard | Spatial Heterogeneity in Fishing Creates de facto Refugia for Endangered Celtic Sea Elasmobranchs[END_REF]. Linked to the resilience of the taxa after fishing event.

We used Hill-Smith analysis [START_REF] Hill | Principal component analysis of taxonomic data with multi-state discrete characters[END_REF], which is a principal component analysis that considers both quantitative and qualitative variables, to assess the sensitivity of each taxon based on the six life history traits selected. The first axis of the Hill-Smith analysis explained the most variability in the distribution of life history traits across taxa. It was used as a sensitivity score (the more positive the score, the more sensitive the taxa are) and also reflects the life history strategy of the taxon.

Community assembly processes

To identify community assembly processes related to life history traits over the 2000-2016 period, we calculated the community-weighted variance (CWV) according

to [START_REF] Gaüzère | A framework for estimating species-specific contributions to community indicators[END_REF] for the sensitivity score (and each of the six life history traits, see Fig. S2). Assembly rules were assessed using a standard effect size (SES) by comparing the CWV to a null model based on random community assembly (i.e.

1000 random permutations of trait values among taxa) [START_REF] Gotelli | Species co-occurrence: A meta-analysis of J. M. Diamond's assembly rules model[END_REF].

The null model consists in shuffling the traits of occurring taxa at each site so that taxa richness remains constant at each site. If the observed CWV was not significantly different from the random CWV, stochasticity was at play. Significant deviation of the CWV from the random mean CWV and negative SES values indicated underdispersion, meaning that the observed trait variability was lower than expected by chance. Conversely, significant deviation and positive SES values indicated overdispersion and an observed trait variability higher than expected by chance.

The SES was calculated for each sampling site over the 2000-2016 period as ( 1).

(

) 1
with the observed CWV at site i, the mean of the theoretical values simulated by the null models at site i (hereafter referred as "mean random CWV") and the standard deviation of these theoretical values simulated by the null models at site i.

The significance of the deviation of the observed CWV from the random mean CWV was calculated with a permutation test, performed with the krandtest function of the "ade4" R package [START_REF] Thioulouse | Multivariate Analysis of Ecological Data with ade4[END_REF]. P-values of these tests were adjusted for multiple testing [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF]. Inverse-distance-weighting interpolation was used only for graphical representation.

Life history trait distribution as a function of environmental and fishing variables

To relate the spatio-temporal distribution of the sensitivity score and the life history traits to environmental and fishing variables over 2000-2016 at the scale of the Celtic Sea, community-weighted means (CWM) were calculated using the "FD" package [START_REF] Laliberté | Package 'FD'. Measuring functional diversity from multiple traits, and other tools for functional ecology[END_REF] for a community of R taxa as followed.

with p i the log(x+1) transformed biomass of taxon i, and t i the value of the sensitivity score of taxon i (or each life history trait, see Fig. S3). The CWM of the sensitivity score was modelled with a Gaussian distribution, so it was first transformed to make it positive before applying the log function, ie the CWM of the sensitivity score was log(x+c+1)-transformed, with c=-min(x) to get a positive or null variable. Like for CWV, inverse-distance-weighting interpolation was used only for graphical representation.

We selected environmental and fishing variables known to influence Celtic Sea ecosystems, and more generally North Atlantic shelf ecosystems [START_REF] Foveau | Process-driven and biological characterisation and mapping of seabed habitats sensitive to trawling[END_REF][START_REF] Mérillet | Environment outweighs the effects of fishing in regulating demersal community structure in an exploited marine ecosystem[END_REF] [START_REF] Mérillet | Environment outweighs the effects of fishing in regulating demersal community structure in an exploited marine ecosystem[END_REF].

We studied the distribution of sensitivity score CWM as a function of environmental variables (depth, sediment, chlorophyll a and temperature), fishing variables (bottom trawl, pelagic trawl) and time (year). In a preliminary step, a permutation test for Moran's I statistic [START_REF] Moran | Note on continous and stochastic phenomena[END_REF] indicated that spatial autocorrelation was significant for sensitivity score (Moran'I= 0.29; p-value = 0.001).

Collinearity between variables was assessed with a variance inflation factor (VIF) and a threshold of five.

The log(x+c+1)-transformed CWM of the sensitivity score followed a Gaussian distribution and were modelled using simultaneous autoregressive (SAR) models (after testing their advantage compared to GLM models with spatial terms, table S2, and looking at their relationship with environmental variables Fig. S5). This family of models assumes that the response at each location is a function of the explanatory variables but also of neighbouring locations [START_REF] Kissling | Spatial autocorrelation and the selection of simultaneous autoregressive models[END_REF]. The spatialerror model, which is the most reliable SAR model, captures spatial autocorrelation that is not completely explained by the explanatory variables or that is an inherent property of the response variable. To implement this model, we used the errorsarlm function in the "SpatialReg" package [START_REF] Bivand | Comparing implementations of estimation methods for spatial econometrics[END_REF]. A full model was built with all the environmental, fishing and time variables, as well as spatial variables (latitude, longitude, latitude², longitude² and latitude x longitude). The model with the lowest AIC was selected using the dredge function ("MuMIn" R package). Temporal autocorrelation in the residuals of the model was tested with a Durbin-Watson test [START_REF] Durbin | Testing for Serial Correlation in Least Squares Regression[END_REF]). S4). Then a Spearman correlation was computed between the significant slopes of the regressions of the biomass over the time and the sensitivity scores (Fig. 4). Based on a histogram of the sensitivity scores (Fig 4), the most sensitive taxa, which clustered apart from the other taxa (i.e.

before the stall in sensitivity score) were examined at a local scale to investigate more precisely the link between temporal evolution of the biomass and fishing effort.

In that purpose, the temporal trend over 2000-2016 of the biomass of the most sensitive taxa was investigated in each International Council for the Exploration of the Seas (ICES) statistical rectangle.

Since the number of sampling sites is not homogeneously distributed across rectangles, Y i the elevated biomass of a taxa i was calculated in each rectangle r with the following formula [START_REF] Mahe | Manuel des protocoles de campagne halieutique[END_REF]:

With Y r the mean biomass per haul in a rectangle r, for taxa i, A r the area of a given statistical rectangle (all similar in the Celtic Sea), A the sum of the areas of all the rectangles that are included or partially included within the limits of the Celtic Sea considered here. The log(x+1)-transformed elevated biomass of each taxa was modelled, in each statistical rectangle, as a function of the year with a linear model using a Gaussian distribution. P-values associated with the coefficient of the year variable were adjusted for multiple tests [START_REF] Benjamini | The control of the false discovery rate in multiple testing under dependency[END_REF]. For each taxon, the coefficient of the year variable was extracted when significant (p<0.05).

We graphically examined the spatial covariation of biomass and trends in fishing effort over 2000-2016 to assess whether the sampling sites that had an increased biomass of the most sensitive taxa also showed a decrease in bottom trawl effort.

All analysis were conducted R 4.2.1 software (R Core Team, 2020).

Results

Assembly rules based on sensitivity score

Hill-Smith analysis of the selected life history traits enabled us to characterise the extent to which each taxon's combination of life history traits made it sensitive to fishing. The first axis (used as the sensitivity score) explained 34.3% of the total variability, while the second axis explained 15.2%. The most sensitive taxa (i.e. large, long-lived taxa that mature late, bear their offspring and have low fecundity) had the highest sensitivity scores.

To distinguish the dominant assembly process between environmental filtering, competition and stochasticity at each site, the CWV of the sensitivity score was compared to a random CWV over the 2000-2016 period. At the majority of the sampling sites (1142 over the 1175 sampling sites), the difference between the observed CWV and the mean random CWV of the sensitivity score was not significant, meaning that stochasticity was the dominant assembly process at these sites. No sampling sites had significant negative SES values, indicating that variability in sensitivity score was not significantly lower than that expected by chance (i.e. underdispersed traits). Few significant positive SES values (i.e. wider range of traits values than that expected by chance), i.e. 33 out of 1175, were visible at sampling sites in the centre of the Celtic Sea and indicated overdispersion (Fig. 2).

Environmental drivers of the distribution of sensitivity score and traits

Overall, most of the Celtic Sea showed low values of sensitivity to fishing, which was indicated by the negative values of CWM of the sensitivity score (Fig. 3).

The biomass of sensitive taxa to fishing was highest in the centre-east of the Celtic Sea (Fig. 3). This area also had the widest range of trait values expressed (i.e. overdispersion), indicated by the significant positive SES values of sensitivity score (Fig. 2), and the significant positive SES values associated with maximum length, age at maturity and offspring size (Fig. S2). Thus, sensitive taxa, with large maximum length, age at maturity and offspring size were also present in the centre-east in addition of the low sensitivity taxa present in all the Celtic Sea. Having a closer look at each trait separately, we showed that only longevity showed a slightly different pattern, with the narrowest range of trait values in the north and the widest in the south (Fig. S2). For reproductive guild, two of the three levels (i.e. bearer and guarder) had their highest proportions of biomass in the centre and extreme northeast of the Celtic Sea and were nearly absent elsewhere. Thus, all three reproductive guilds occurred in these two areas (Fig. S2), which agreed with the overdispersion observed there (Fig. S2). Similarly, three levels of the four fecundity levels (i.e. low, medium and very high) had some of their highest proportions of biomass in the centre and were nearly absent elsewhere, which led to overdispersion of the fecundity trait in the centre (Fig. S2). The VIF of the model with environmental, fishing and time variable was below 5, so all variables were kept for the full model. Spatial variables were also kept to account for spatial autocorrelation. Moreover, Moran's I indicated no residual autocorrelation (p >0.05) in the selected SAR model of the CWM of the sensitivity score, which was not the case for the GLM (Table S2). In addition, no temporal autocorrelation could be found in the residuals of the selected SAR model for sensitivity score (Durbin-Watson statistic = 1.40, p-value= 1.96.10 -1 ), neither in the residuals of each traits modelled separately (Fig. S2). The spatial variables (latitude, longitude, their quadratic effects and interaction) as well as bottom trawl effort, sediment and year were significant drivers of the CWM of the sensitivity score (Table 2). Traits taken separately brought more details on which trait forming the sensitivity score responded to which variables. The models of the sensitivity score and each of the life history traits highlights the importance of the habitat, described by sediments and depth. The biomass of sensitive taxa tends to be higher on mud and sand sediments than on mix and coarse (boulders taken as a reference) (Table 2).

Depth was not a significant driver of the CWM of the sensitivity score, but was a significant structuring variable of some of the CWMs made on the traits taken separately: the relative biomass of taxa with large maximum length and high age at maturity increased in shallower areas (Table S3). Conversely, the relative biomass of taxa that are long-lived, bear their offspring (i.e. bearer) and have low fecundity, as well as that of non-guarder taxa and high fecundity taxa, tended to increase in deeper areas (Table S3).

All CWMs of the sensitivity score and each trait taken separately covaried negatively with bottom trawl effort (Table 2 and Table S3). Areas with little bottom trawl effort tended to have the highest relative biomass of large, long-lived taxa that have low fecundity, mature late and give birth to large offspring. The relative biomass of taxa with very high fecundity also decreased as bottom trawl effort increased.

Conversely, areas with large bottom trawl effort had the highest relative biomass of egg-guarding taxa and high fecundity taxa. In comparison, the spatial covariation of CWMs with pelagic trawl effort was significant only for the longevity trait and the nonguarder level. Pelagic trawl effort had a negative influence on the relative biomass of long-lived taxa but a positive influence on that of non-guarder taxa (Table S3).

Regarding temporal trends in CWMs in the Celtic Sea, the relative biomass of sensitive taxa increased over 2000-2016 (Table 2). More precisely, the relative biomass of large taxa that mature late, have low fecundity and large offspring increased from 2000-2016. Conversely, the relative biomass of taxa that guard their eggs and have medium or high fecundity decreased from 2000-2016 (Table S3). As bottom temperature increased, the relative biomass of long-lived taxa increased, while that of non-guarder taxa and taxa with very high fecundity decreased. Taxa with very high and medium fecundity tended to have lower relative biomass at sites with higher chlorophyll a concentrations (Table S3).

3.3. Local variation in the biomass of the most sensitive taxa in relation to local fishing patterns and closure area Among all taxa, the temporal trend in biomass significantly increased for 19 taxa, significantly decreased for 10 taxa and had no significant trend for 72 taxa (Table S4). At the scale of the Celtic Sea, biomass appeared to increase for the most sensitive slow life history taxa (Fig. 4). This is confirmed by the positive correlation between significant slope of the regression (p<0.05) of the biomass of each taxa over time and sensitivity scores (Spearman's ρ = 0.52; p = 0.004) (Fig S6 ).

Seven taxa had a distinctly higher sensitivity score than the other taxa ( Areas where the increase happened correspond to the ones that experienced a significant decrease in bottom trawling effort: the seasonal area closure Trevose box and its adjacent area. Nevertheless, the relationship was difficult to assess for S.

stellaris because of its small range in the Celtic Sea (Fig. 5 E

Discussion

To evaluate impacts of fishing on community structure, the distribution of life history traits that are sensitive to fishing was studied to identify assembly rules, as well as their environmental and anthropogenic drivers. At the scale of the Celtic Sea, sensitivity to fishing was mostly distributed according to stochastic assembly process.

Significant overdispersion was mainly detected in the centre-east of the Celtic Sea. Values of life history traits are more diverse than expected by chance there, which suggests more diversified life history strategies. This overdispersion in traits is driven by the co-occurrence of taxa with a slow life history strategy, such as the common skate Dipturus batis, and a fast life history strategy, such as the Atlantic horse mackerel Trachurus trachurus. Three non-exclusive deterministic processes could explain the distribution of life history traits in the centre-east of the Celtic Sea:

(i) Competition among coexisting taxa fosters trait overdispersion.

Competition might lead to either competitive exclusion or to spatial and/or temporal avoidance among co-occurring taxa [START_REF] Cahill | Does phylogenetic relatedness influence the strength of competition among vascular plants?[END_REF][START_REF] Mayfield | Opposing effects of competitive exclusion on the phylogenetic structure of communities[END_REF]. However, in the centre of the Celtic Sea were overdispersion occurs, competition between fast and slow life history taxa does not appear as a reasonable explanation of the large range of sensitivity score expressed, as for instance skate and Atlantic horse mackerels do not compete for food.

(ii) Small-scale seascape heterogeneity provides a greater diversity of niches, in agreement with the habitat heterogeneity hypothesis [START_REF] Tews | Animal species diversity driven by habitat heterogeneity/diversity: The importance of keystone structures[END_REF].

We observed that sediments, which are generally homogenous throughout the Celtic Sea (mud and sand), tend to be patchier where overdispersion occurs (i.e. centre-east). Life history traits define one of at least five dimensions of the species ecological niche (i.e. habitat, life history, trophic, defence and metabolic) [START_REF] Winemiller | Functional traits, convergent evolution, and periodic tables of niches[END_REF]).

The area where overdispersion was observed could thus provide more available niches than the rest of the Celtic Sea and allow more life history strategies to coexist in accordance with the limiting similarity process [START_REF] Macarthur | The Limiting Similarity , Convergence , and Divergence of Coexisting Species Author ( s )[END_REF].

(iii) Decrease in fishing pressure alleviated the filtering imposed on sensitive taxa. Overdispersion might occur when new taxa with different trait values arrive or due to changes in taxa biomass among trait values. For the Celtic Sea, we observed stability in the taxonomic structure (no large changes in biomasses of taxa) of these communities despite spatio-temporal changes in anthropogenic pressures [START_REF] Mérillet | Environment outweighs the effects of fishing in regulating demersal community structure in an exploited marine ecosystem[END_REF]. Focusing on traits, we highlight that the biomass of slow life history taxa increased over time in the centre-east and extreme north-east, thereby increasing the frequency of these original trait values in the community and flattening the distribution of trait values. In the extreme north-east, the biomass of sensitive taxa could have increased due to the decrease in fishing pressure. It is likely that when fishing was still high, it acted as a major filter on biodiversity as demonstrated in [START_REF] Guénette | Shifting baselines in European fisheries: The case of the Celtic Sea and Bay of Biscay[END_REF]. Sensitive taxa may have been filtered from the areas with intensive fishing, whereas the decrease in fishing allowed more diverse taxa to coexist again.

Stochasticity was the dominant assembly process over the majority of the Celtic Sea. Assembly processes thus appeared to be dominated by demographic processes (birth, death, immigration, emigration) occurring at random with respect to the taxon's identity [START_REF] Vellend | Assessing the relative importance of neutral stochasticity in ecological communities[END_REF]. This seems in line with the high number of mobile fish in the taxa studied and dominance of dispersion of egg as reproductive strategy. However, this might also result from the failure to detect deterministic processes [START_REF] Vellend | Assessing the relative importance of neutral stochasticity in ecological communities[END_REF] which can be due to the large number of taxa in this study interacting together and reacting to abiotic variables in many different ways.

Underdispersion (even if non-significant) observed in the most of the Celtic Sea could be due to environmental filtering that selects taxa that can endure specific environmental conditions. Depth, sediments, year and bottom trawl effort are the variables that significantly drive the largest number of traits and have thus a particularly structuring effect on their distribution. The dependence of the CWMs on variables that do not vary over time at a given site (i.e. depth and sediments)

illustrates the strong structuring effects of habitat. For variable that vary over time and space, bottom trawl effort thus appears to be a severe condition that restricts the distribution of life history traits in the Celtic Sea, with a globally negative relationship with the biomass of slow life history strategy taxa demonstrated in all CWMs models.

We observed temporal variation in the biomass of slow life history strategy taxa that increased from 2000-2016 at the regional scale. At the scale of the Celtic Sea, bottom trawling was stable over the period studied [START_REF] Mérillet | Environment outweighs the effects of fishing in regulating demersal community structure in an exploited marine ecosystem[END_REF], but this stability masks local differences that are visible at the scale of the statistical rectangle, as showed here. Almost all significant increases in the biomass of the most sensitive taxa occurred in statistical rectangles in which bottom trawling decreased or showed no significant trend. This is particularly clear for Mustelus sp., S. stellaris, D. Batis and S. acanthias, whose populations increase in areas where fishing pressure decreased. The biomass of taxa with a slow life history, especially elasmobranchs, is higher in the centre-east of the Celtic Sea, and the highest part of the increase in biomass of the five elasmobranchs occurred in adjacent areas: the extreme north-east and south-west. This pattern could be due to a local increase in abundance after a decrease in fishing pressure (in the north-east), but also to a replenishment from the larger number of individuals present in the centre-east of the Celtic Sea [START_REF] Green | Designing Marine Reserves for Fisheries Management, Biodiversity Conservation, and Climate Change Adaptation[END_REF]. As an open sea with relatively homogenous sediment, the Celtic Sea appears to be a seascape that large-bodied and mobile taxa can cross easily. We hypothesize that the centre of the Celtic Sea could have functioned as a refuge and then a source for recolonization of the north-east when fishing pressure began to decrease. A similar recovery of a depleted population from adjacent areas in a fishing ground was observed for predatory fish on the north-west Atlantic shelf [START_REF] Shackell | Spatial scale of similarity as an indicator of metacommunity stability in exploited marine systems[END_REF].

The observed increase in sensitive slow life history taxa could be related to a local decrease in fishing pressure due to management practices. In 2009, the EU began to enforce its Common Fisheries Policy, which decreased overall fishing effort in EU seas [START_REF] Walter | The EU's Common Fisheries Policy: A Review and Assessment[END_REF]. In this context, fishing effort was displaced from the northeast to the west and south-west of the Celtic Sea and decreased bottom trawl effort in certain statistical rectangles. It may also have benefited the prey of slow life history taxa, as demonstrated in Georges Bank, a formerly overexploited ecosystem east of Massachusetts, USA [START_REF] Mayo | Aggregate Fish Biomass and Yield on Georges Bank, 1960-87[END_REF], although we did not observe an overall increase in biomass of small pelagic fish. This lack of evidence may have been due

to incomplete sampling of pelagic communities. Moreover, the seasonal fishing closures of the which confirms the hypothesis of the benefits of this closure [START_REF] Horwood | Evaluation of Closed Areas for Fish Stock Conservation[END_REF].

Indeed, we observed the largest significant decrease in bottom trawl effort and increase in biomass for four of the five aforementioned elasmobranch taxa in this area. This area is also a known spawning ground and nursery for many species, including elasmobranchs [START_REF] Potter | Seasonal, annual and regional variations in ichthyofaunal composition in the inner Severn Estuary and inner Bristol Channel[END_REF], and could, in addition to serving as a refuge for adults from fishing pressure, favour the replenishment of the population by increasing reproductive success. The increase in biomass of slow life history taxa observed in the Celtic Sea is in line with the increase in equilibrium strategy taxa (e.g. elasmobranchs) observed in the North Sea [START_REF] Pecuchet | From traits to lifehistory strategies: Deconstructing fish community composition across European seas[END_REF] following a decrease in fishing effort and an increase in temperature.

A negative relationship was observed between bottom temperature and the biomass of short lifespan taxa as well as very high fecundity taxa (characteristic of the fast life history strategy). In addition, we also showed that taxa with a significant decrease in biomass tended to be least sensitive to fishing and had fast life history traits (i.e. small maximum length, short-lived, early-maturing, small offspring, nonguarder and high or very high fecundity). This decrease in the biomass of fast life history taxa could thus arise from the significant increase in bottom temperature that occurred in the north of the Celtic Sea over 2000-2016 [START_REF] Mérillet | Environment outweighs the effects of fishing in regulating demersal community structure in an exploited marine ecosystem[END_REF]. This decrease in fast life history taxa was also observed in the English Channel after a warm phase of the Atlantic Multidecadal Oscillation [START_REF] Mclean | A Climate-Driven Functional Inversion of Connected Marine Ecosystems[END_REF]. Fast life history taxa are highly responsive to an increase in temperature due to their short generation time [START_REF] Mclean | A Climate-Driven Functional Inversion of Connected Marine Ecosystems[END_REF]. In the North Sea, mixed trends were observed since the biomass of periodic strategy taxa have decreased over the past 30 years, but opportunistic taxa (at the fasted end of the continuum of life history strategy) have increased [START_REF] Pecuchet | From traits to lifehistory strategies: Deconstructing fish community composition across European seas[END_REF]. Overall, as in adjacent area, the observed life history patterns in the Celtic Sea most likely result from an interplay between temperature and fishing.

Conclusion

Analysis at various spatial and temporal scales enabled us to draw some conclusions about the recent evolutions of sensitivity patterns to fishing in the Celtic Sea. Assembly processes are mostly explained by stochasticity in the demographic events. Nevertheless, environmental filtering applied by fishing pressure on the taxa also appeared to have had an importance. In this historical fishing ground, fishing probably acted as a filter in the north-east and centre-east of the area before the beginning of the time series studied here. This filter got slowly removed allowing significant limiting similarity to be revealed (in a heterogeneous habitat or from the competition between taxa).

We showed the negative relationship between fishing effort and the biomass of taxa exhibiting sensitive traits, but also the increase in biomass of the most sensitive taxa from 2000-2016, in link with the creation of the Trevose box area closure. Rapid changes in the functional structure of the communities are expected in the future as climate change effects become more tangible and superimpose to fishing ones. Management measures leading to short-term improvement of the biomass of sensitive taxa to fishing are thus urgently needed [START_REF] O'leary | The resilience of marine ecosystems to climatic disturbances[END_REF].

This study highlights potential positive dynamics for a group of sensitive species due to a local decrease in fishing pressure. Within two decades, mitigation measures, such as area closures that restrict fishing effort spatially or temporally on essential habitats such as spawning and nursery grounds could significantly increase the biomass of the most sensitive taxa, which can then act as a source for recolonization.

This positive response to recent EU policies could be a useful feedback for managers. However, a continual caution during ecosystem management of fishing grounds is required, as commercially fished small pelagic taxa could experience a decrease in biomass in the future due to climate change and could eventually be replaced by warm water pelagic taxa.

Figure 1 .

 1 Figure 1. Map of sampling sites (grey dots) in the study area in the Celtic Sea. The red line delineates the limits of the study area, as defined in the EVHOE surveys.

2. 4 .

 4 Local temporal trends in the biomass of the taxa most sensitive to fishing Based on a histogram, the sensitivity score distribution was examined. To inform a global trend over the Celtic Sea between sensitivity and dynamics of the biomass, temporal trends of the log(x+1)-transformed biomass of each taxa in the whole Celtic Sea was modelled (biomasses were elevated to the whole Celtic Sea following (Mahe and Poulard, 2005)), as a function of the year with a linear model using a Gaussian distribution (Fig 4, Fig. S4, Table

Figure 2 .

 2 Figure 2. Maps of standard effect size of the community-weighted variance (CWV) of the sensitivity to fishing computed from six life history traits of 101 taxa in the Celtic Sea. Sampling sites surrounded by circles show significant departure from the stochastic assembly process. Positive values (in red) indicate overdispersion in the traits and negative values (in blue) indicate underdispersion.

Figure 3 .

 3 Figure 3. Maps of community-weighted means (CWMs) of the sensitivity to fishing computed from six life history traits of 101 taxa in the Celtic Sea.

  Fig. 4), namely spurdog Squalus acanthias, school shark Galeorhinus galeus, blue skate Dipturus batis, smooth-hound Mustelus sp., velvet belly Etmopetrus spinax, rabbit fish Chimaera monstrosa and nursehound Scyliorhinus stellaris (i.e. sharks, rays and chimaeras).

Figure 4 .

 4 Figure 4. Sensitivity scores from Hill-Smith analysis of the six life history traits of 101 taxa in the Celtic Sea, ranked from most to least sensitive to fishing. Sensitivity scores of individual taxa coloured according to the quartile of the slope of the regression of biomass (g.m -2 ) from 2000-2016. Non-significant slopes (p >0.05) are coloured in white. The most sensitive taxa before the stall in sensitivity values are framed in red. An asterisk indicates significant slope after adjusting p-values for multiple testing.

  .). Along the southern border of the Celtic Sea, biomass of Mustelus sp. and to a lesser extent G. galeus increased significantly in statistical rectangles with no significant trends in bottom trawling (Fig. 5 B. and D.). In addition, biomass of G. galeus significantly decreased in a statistical rectangle where bottom trawl effort significantly increased, in the centre of the Celtic Sea.

Figure 5 .

 5 Figure 5.(A to E) Map of the slope of the regression of biomass (g/m 2 , black circles) in each statistical rectangle (111.12 km × 55.56 km) of the Celtic Sea from 2000-2016. This figure focussed on the taxa the most sensitive to fishing that have a significant variation of their biomass in at least one statistical rectangle. Small open circles represent the distribution of the taxa and filled circle the significant (p<0.05) slope of the regression of the biomass over time. Statistical rectangles are coloured according the sign of temporal trend of bottom trawl effort (hours per year within the rectangle) on time (only slopes with p <0.05 are coloured). (F) Map of the seasonal area closure Trevose Box.

  . Depth was recorded at each sampling site. Sediment data according to Folk's 5-level classification system came from the EMODnet Geology Portal (https://www.emodnet-geology.eu). Data for the few (95) sampling sites that were not included in the EMODnet dataset were extrapolated from the SHOM database (https://diffusion.shom.fr/loisirs/sedim-mondiale.html). As our community data were based on annual sampling, and the 101 taxa studied display a wide variety

of life cycles, we assumed that annual means would be the most relevant values for the covariables. Bottom temperature data were downloaded from Copernicus (http://marine.copernicus.eu), from the Atlantic-Iberian Biscay Irish-Ocean Physics Reanalysis data product at 8 km × 8 km resolution, while chlorophyll a data were extracted from satellite observations at 1 km × 1 km resolution. Fishing effort data were derived from the New Fisheries Dependent Information from the EU's Scientific, Technical and Economic Committee for Fisheries (https://stecf.jrc.ec.europa.eu), which provided data on fishing effort for all EU countries in hours fished per statistical rectangle (111.12 km × 55.56 km). Bottom and pelagic trawl were the two dominant métiers over the Celtic Sea

Table 2

 2

	Trait	Depth	Sediment		Chlorophyll a	Bottom temp..	Bottom trawl	Pelagic trawl	Year	Lat.	Long.	Lat. 2	Long. 2	Lat. × Long.
	Sensitivity score		Coarse	-7.08×10	-3		-4.05×10	-6		5.69×10-3	4.20 ×10-2	1.36 ×10-1	-7.53 ×10-2	-4.96 ×10-2	7.52 ×10-3
			Mix	-7.50×10	-2
			Mud	7.54×10	-2
			Sand	6.56×10	-2

. Estimates of explanatory variables in spatial autoregressive model of the sensitivity score (log(x+c+1)transformed). Only significant estimates are shown (p<0.01).

Acknowledgement

We thank all those who made data sampling possible during the EVHOE surveys. This work is part of a PhD thesis jointly funded by the French Institute for the Exploitation of the Sea (IFREMER) and the National Museum of Natural History (MNHN).


Statistical Biases and Very-Long-Term Time Stability Analysis

François Vernotte and Éric Lantz

Abstract-The prediction of very-long-term time stability is a key issue in various fields, such as time keeping, obviously, but also navigation and spatial applications. This is usually performed by extrapolating the measurement data obtained by estimators such as the Allan variance, modified Allan variance, Hadamard variance, etc. This extrapolation may be assessed from a fit over the variance estimates. However, this fit should be performed on the log-log graph of the estimates, which corresponds to a leastsquares minimization of the relative difference between the variance estimates and the fitting curve. However, a bias exists between the average of the log of the estimates and the log of the true value of the estimated variance.

This paper presents the theoretical calculation of this loglog bias based on the number of equivalent degrees of freedom of the estimates, shows simulations over a large number of realizations, and provides a reliable method of unbiased logarithmic fit. Extrapolating this fit yields a more confident assessment of the very-long-term time stability.

I. I

T  stability of clocks and oscillators is generally as- sessed by estimators such as the Allan deviation, modified Allan deviation, Hadamard deviation, etc. However, unlike the variance estimates, the deviation estimates are severely biased and should not be used for further calculations, such as fit.

On the other hand, the prediction of very-long-term time stability is usually performed by fitting and extrapolating an Allan variance (AVAR) curve. Such a curve may be modeled as σ ττττ

y i ii i i C 2 =0 4 2
()=( )( )= .

∑ -ΦΦ ,w ith

An extrapolation may be achieved by extending the asymptotes of such a curve beyond the largest τ value (see Fig. 1).

The C i coefficient are related to the noise levels h α of the power law noise model by Table I, [START_REF] Rutman | Characterization of phase and frequency instabilities in precision frequency sources: Fifteen years of progress[END_REF]:

S fh f y ()= =2 2 α α α - + ∑ , (2)
and to the linear frequency drift D 1 :

yt DD tt ()=( ) 01 ++ η , (3)
where η(t) represents the pure random component of the frequency deviation y(t).

However such an extrapolation is very sensitive to estimation errors. It is then of importance to perform a reliable estimation of the C i coefficients. The weighted relative least squares are an effective method [START_REF] Vernotte | Oscillator noise analysis: Multivariance measurement[END_REF] which is weighted, to take into account the uncertainty ∆ j over each estimate ˆ, σ τ yj 2 () and relative, to be equivalent to a classical least square fit on the log-log plot by minimizing the relative distance between the curve and the estimate [(estimate -curve)/estimate],

j N j yj i ii j yj C =1 2 =0 4 2 2 1 () () () ∑ ∑ × -               ∆ Φ ˆστ τ στ minim mum. (4)
It is important to emphasize that for the relative fit, i.e., the fit of ∆ ˆσσ 22 / and not of ∆ σ 2 , that 1) it ensures that the fit will minimize the relative differences between the curve and the estimates (e.g., 10% for an estimate at the 10 -14 level as well as 10% for an estimate at the 10 9 level), whereas a classical fit would minimize the absolute differences, overwhelming the small estimates; and 2) because dx/x = d ln x, it is exactly equivalent to a classical fit on the log of the estimates.

As a consequence, this extrapolation method is based on the logarithm of the estimates and not the estimates themselves. We know that the estimate σ τ

y 2 () is an unbi- ased estimator of σ τ y 2 () , but it is not known if ln(ˆ) στ y 2 () is an unbiased estimator of ln() στ y 2 ().

II. S A

A. Statistics of the Allan Variance and the Allan Deviation

The definition of the Allan variance is [START_REF] Allan | Statistics of atomic frequency standards[END_REF], [START_REF] Barnes | Characterization of frequency stability[END_REF]

σ τ y yy 2 21 2 ()= 1 2 . 〈〉 - () (5)
The estimate of the Allan variance may be calculated as F. Vernotte is with UTINAM, Observatory Terre-Homme-Environnement-Temps-Astronomie (THETA) of Franche-Comté, University of Franche-Comté/Centre National de la Recherche Scientifique (CNRS), Besançon, France (e-mail: francois.vernotte@obs-besancon.fr).

E. Lantz is with Department of Optics, FEMTO-ST, University of Franche-Comté/CNRS, Unité de Formation et de Recherche (UFR) Sciences, Besançon, France. ∑ -

() (6)
Because y 2 and y 1 are Gaussian values

2 χ ν = 2ν (s td ⇒ () 2 χ ν = 2ν
). As shown by Fig. 2, the number of EDF reflects the AVAR estimates dispersion. A method for assessing the EDF may be found in [START_REF] Greenhall | Uncertainty of stability variances based on finite differences[END_REF].

B. Model World and Measurement World

Let us define [see Fig. 3(a)] θ as the model parameter and ξ as the measured data of the parameter. For example, we could consider the true value σ τ Two problems are generally addressed: the direct problem, which attempts to forecast the measurement data knowing the parameter and the inverse problem, which aims to estimate the parameter knowing the measurement data. In the same vein, Tarantola distinguishes the model space, i.e., the space in which the parameter is given, from the data space, i.e., the space in which the measurement data are given [START_REF] Tarantola | Inverse Problem Theory[END_REF]. In the following, we will use the terms model world and measurement world.

1) Model World (Direct Problem):

In the direct problem, the task is to determine how the measurement data ξ are distributed when the parameter θ 0 is known [see Fig. 3(b)]. However, the parameter θ 0 is precisely the unknown quantity that we want to estimate. This parameter is only known in theory and simulations.

2) Measurement World (Inverse Problem):

In the inverse problem, the task is to estimate a confidence interval over θ from the known measurement data ξ 0 [see Fig. 3(c)]. This is the usual task of the metrologist.

C. Study of a χ 2 Distribution With ν = 1 Degree of Freedom 1) Properties:

The properties of a χ 1 2 distribution are discussed in this section [START_REF] Evans | Statistical Distributions[END_REF]. The probability density function (PDF) is

S y ( f ) h +2 f +2 h +1 f +1 h 0 f 0 h -1 f -1 h -2 f -2 y(t) = D 1 t σ τ y 2 () 3 4 2 22 hf + h πτ 1.038 3( 2) 4 1 22 + [] + ln πτ πτ fh h h 0 2τ 2( 2) 1 ln h - 2 3 2 2 πτ h - D 1 22 2 τ
For the high frequency noises ( f +2 FM and f +1 FM), it is necessary to introduce a high cut-off frequency f h [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF].

Fig. 1. Extrapolation of a least square fit. The graph shows the Allan deviation (i.e., the square root of the Allan variance) but the fit was performed over the Allan variance estimates according to the criterion (4) using SigmaTheta-1.2 [START_REF] Vernotte | Sig-maTheta software package[END_REF]. The equation of the curve is [(4.7

• 10 -13 ) 2 τ -2 + (3.1 • 10 -13 ) 2 τ -1 + (4.8 • 10 -14 ) 2 + (1.2 • 10 -15 ) 2 τ +1 + (5.2 • 10 -18 ) 2 τ +2 ] 1/2 . p X e X X () = 2 2 -/ π . ( 7
)
The PDF is normalized:

0 () =1 ∞ ∫ pX X d. (8)
The mathematical expectation is

E XX pX X () =( )= =1 0 ∞ ∫ ⋅ d ν . (9)
We can define the logarithmic mean as the exponential of the mathematical expectation of ln (X )

e Xp XX e EX C (ln) expl n , () 0 =( )( ) = 2 0.2807 ∞ - ∫ ⋅       ≈ d ( 10
)
where C is the Euler constant; C ≈ 0.57722. The cumulative distribution function (CDF) is

P Xp yy X X () =( )= 1 (1 2, 2) 0 ∫ - d // inc Γ π , (11)
and the inverse CDF is

P -- - 11 ()=2 12 ,( 1) απ α Γ inc / () . ( 12
)
This relationship is useful for getting the bounds of a confidence interval. How is the reduced variable X of the relationships ( 7) to [START_REF] Benford | The law of anomalous numbers[END_REF] related to the measurement data ξ and the model parameter θ ?

2) Parameter, Measurement, and Reduced Variable: Let us consider the standard χ 1 2 variable X = x 2 where x is a Gaussian centered standard random variable:

E XE X () == 1= 1. ν ν ⇒ ( ) (13)
It is well known that σ τ

y 2 () = ξ is χ 1 2 distributed; ξ is an unbiased estimator of the parameter σ τ y 2 () = θ. Therefore, E ξ θ ( ) =1. ( 14
)
We can then define the reduced variable X as

X == . ν ξ θ ξ θ ( 15
)
The differential dX is then

dd d X XX =. ∂ ∂ + ∂ ∂ ξ ξ θ θ (16)
3) Conditional Probabilities in the Model World (Direct Problem): The parameter θ is a constant. Let θ = θ 0 ; the differential dX is then: From ( 17) and the relationships ( 7) to ( 9), we find the conditional probability density function:

p e (| )= 2 . 0 (2 ) 0 0 ξθ πθ ξ ξθ -/ (18)
ξ is unbiased:

E p ()=( |) =. 0 00 ξξ ξθ ξθ ∞ ∫ ⋅ d ( 19
)
The logarithmic mean is

e p e E C ln expl n () 0 0 00 =( )( |) = 2 0.2807 . ξ ξξ θξ θθ [] ∞ - ∫ ⋅       ≈⋅ d (20)
Thus, the mathematical expectation of the log of the measurement data does not converge to the log of the parameter. Therefore, there is a significant bias for log-log representation (and fit); let us call it the log-log bias.

4) Conditional Probabilities in the Measurement World (Inverse Problem):

The measurement data ξ is a constant. Let ξ = ξ 0 ; the differential dX is then

dd X =. 0 2 - ξ θ θ (21)
From (21) and the relationships ( 7)-( 9) comes the conditional probability density function:

p e (| )= 2 . 0 0 3 (2 ) 0 θξ ξ πθ ξθ -/ (22)
There is no mean value,

E p ()=( |) 0 0 θθ θξ θ ∞ ∫ ⋅ d ( 23
0 0 00 =( )( |) =2 3.562 θ θθ ξθ ξξ [] ∞ ∫ ⋅       ≈⋅ d )
In the measurement world also, the mathematical expectation of the log of the parameter does not converge to the log of the measurement data. However, this bias is exactly the inverse of the log-log bias defined previously. (This property will be demonstrated later in a more general context, see Section II-D.) This implies that we just have to renormalize the measurement data by the log-log bias to get both an unbiased estimator in the model world and a mathematical expectation of the parameter equal to the measurement data in the measurement world.

D. Generalization to Any Positive-Valued Random Process

Let ξ be a positive-valued estimator constructed with M measurement data resulting from a random measurement process. ξ is itself a random variable in the model world, which we assume to have a mean. The most obvious example is Allan variance. Another example is the estimation of the standard deviation. It is quite generally possible to define a new reduced variable,

Y k = 0 ξ θ , (25)
which obeys ( 13), unlike X:

E YY pY Y () =( )= 1 0 ∞ ∫ ⋅ d. (26)
Let Z = f (Y) be a monotonic function of Y. A fundamental property of density probabilities is

p YY pZ Z YZ () =( ) dd , (27)
where the subscripts emphasize the fact that the different density probabilities are different functions (these subscripts are, as usual, omitted in the following). An immediate consequence is, for Z = ξ,

0 0 () = ∞ ∫ ⋅ ξξ ξθ pk d/ . ( 28
)
If X is a standard χ M 2 variable with ν = M degrees of freedom, then Y = X/ν ; with ξ a variance measurement data, we have, in the model world, k = 1 and (28) means that ξ is an unbiased estimator of the true variance θ 0 . In the case of standard deviation estimation, the square root of the Allan variance is a biased estimator, but the value of k can be easily calculated as a function of M (the estimator is convergent: k → 1 if M → ∞) and (28) allows also the definition of an unbiased estimator.

The situation is worse in the measurement world. The variable is, in this case, θ :

Y k = 0 ξ θ , (29)
with, once more,

p YY p () =(|) 0 dd θξ θ. ( 30
)
Eqs. ( 29) and (30) give, for k = 1,

0 00 0 0 (| )= 1 () ∞∞ ∫∫ ⋅⋅ ⋅≠ θθ ξθ ξξ p Y pY Y dd . (31)
We have seen previously that this latter expression may not be convergent.

Let Z = ln (Y ). We can easily calculate

Eq. ( 33) means that ln (ξ ) -B is an unbiased estimator of ln (θ 0 ). Hence, the log-unbiased estimator of θ 0 appears to be, in the model world, ξ ′, defined as

′ ⋅ - ξ ξ =. e B ( 34
)
In the measurement world, we have, from (29

) with k = 1, Y e B = 0 ′ ⋅ ξ θ , (35)
and then

ln ln ln ()=( )( ). 0 θξ ′ -+ YB (36)
Unlike linear estimators, this estimator ensures that its mean is equal to the actual measurement data: (37)

To conclude this section, a log-estimator can easily be rendered unbiased by multiplying the linear estimator by the right constant, with the important further advantage that, in the measurement world, the unknown true value has, for its mean, the actual measurement data. This equivalence between the model and measurement worlds for a scale parameter was already noticed 50 years ago [START_REF] Lindley | Fiducial distributions and Bayes theorem[END_REF].

E. Generalization to a χ ν 2 Distribution

The results of Section II-C can be extended to any χ ν 2 distribution, regardless of the number of degrees of freedom ν.

1) Model World: The reduced variable is .

X = 0 ν ξ θ , (38)

Γ (39)

The mean is

E()= 0 ξθ . ( 40
)
The logarithmic mean is

e E ln exp, () (0) 0 = 2 2 ξ ν ψ ν θ [] ( )       ( 41
)
where ψ (0) (x) is the digamma function [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF].

2) Measurement World: The reduced variable is

X = 0 ν ξ θ , (42)
and

E(X) = ν. The PDF is p e (| )= 2( 2) 0 2 0 2 2 22 1 0 θξ νξ νθ ν ν νξ θ νν / / // / - + () /( ) . Γ (43)
The mean (ν > 2) is

E()= 2 0 θ ν ν ξ - . ( 44
)
The logarithmic mean is

e E[ln ] exp. () (0) 0 = 22 θ ν ψ ν ξ - ( )       (45)
Table II and Fig. 4 show the biases and the 95% confidence interval versus EDF.

III. D

A. A log-Unbiased Estimator

We know that σ τ y 2

() is an unbiased estimator of σ τ y 2

() . However, E [ln(ξ )] ≠ ln(θ 0 ) and e E [ln(ξ )] ≠ θ 0 . Similarly, E [ln(θ )] ≠ ln(ξ 0 ) and e E [ln(θ )] ≠ ξ 0 .

L θ ν ψ ν ξ = 22 (0) - ( )       { } (46)
and, with ν = 1,

ˆln ln L θ ξξ =( 3.5625 )= () 1.2705. ⋅+ (47)
The mathematical expectation of Lθ is then

〈 〉 ˆln . L θ θ =( ) 0 (48)
Lθ is then an unbiased estimator of ln(θ 0 ).

On the other hand, from (45), for measurement data ξ 0 :

E[ln ]l ne xp ˆ, () = 22 = (0) 0 θ ν ψ ν ξ θ - ( )       { } L (49)
and from (24), with ν = 1,

E[ln ]l n () =( 3.5625 )= . 0 θξ θ ⋅ L (50)
The mathematical expectation of ln(θ ) is equal to ˆ. L θ Therefore, e Lθ is a log-unbiased estimator of θ. For fitting a log-log plot, this log-unbiased estimator is then optimal.

B. Fitting an Allan Variance Curve

Fig. 5 displays the bias estimates (crosses) as well as the log-unbiased estimates (circles) of a simulated oscillator. Fitting the AVAR estimates rather than the logunbiased estimates leads to a significant underestimation of the long-term asymptotes (τ and τ2 asymptotes) and therefore of the very-long-term stability of the oscillator.

C. Simulation Results

1) Pure Random-Walk FM: We first simulated pure random-walk FM, because it is often the dominating noise for the very-long-term. a) Principle of the simulation: We performed one million simulations of a sequence of 1024 frequency deviation samples of pure random-walk FM with a noise level (model parameter) within the interval: 10 -3 < h -2 < 10 +9 with a uniform distribution of the exponent 2 between -3 and +9. Thus, this variation range can cover all of the model parameters which may give a measurement data ĥ-2 equal to 1. It is then possible to select either the realizations for which the model parameter h -2 is close to 1 at 5% (i.e., within the interval [0.952, 1.05]), or the realizations for which the measurement data ĥ-2 is close to 1 at 5%. b) Study of AVAR for the higher τ value: It is well known that an AVAR measurement data for the higher τ value is χ 1 2 distributed. Because each sequence contains 1024 samples, the higher τ argument of AVAR is τ max = 1024/2 = 512 arbitrary units (a.u.). From Table I, the theoretical AVAR result (model parameter) for h -2 = 1 a.u. and τ = 512 a.u. is surement data. To more easily compare these values, we normalize them by the coefficient K n given by (51):

σ τ πτ π yn h K
θ σ π ξ σ = (512) = 1024 3 = (512) = 2 2 2 2 2 21 2 y nn y nn K h K K YY K - -        ˆ() ,        (52)
where Y 1 and Y 2 are, respectively, the average of the first 512 frequency deviation samples and the average of the last 512 frequency deviation samples of a considered simulated sequence. Thus, we can select either the sequences for which θ = θ 0 = 1 ± 5% (model world) or the sequences for which ξ = ξ 0 = 1 ± 5% (measurement world). We can then easily compare the histogram of the ξ measurement data, knowing that θ 0 = 1 ± 5%, to the theoretical conditional probability given by (18) as well as the histogram of the θ parameter, knowing that ξ 0 = 1 ± 5%, to the theoretical conditional probability given by ( 22). Fig. 6 shows the very good agreement between histograms and theoretical conditional probability for both the model world and measurement world.

For the model world, we obtained 3568 sequences for which θ 0 = 1 ± 5% a.u., 〈ξ 〉 = 0.977 a.u., exp(〈ln(ξ )〉) = 0.276 a.u., 〈 〉 exp( ˆ)

L θ = 3.48 a.u., exp( ˆ) 〈〉 L θ = 0.983 a.u. For the measurement world, we obtained 3523 sequences for which ξ 0 = 1 ± 5% a.u., 〈θ 〉 = 19 800 a.u., exp(〈ln(θ )〉) = 3.84 a.u., and 3479 sequences for which exp( ˆ)

L θ = 1 ± 5% a.u., 〈θ 〉 = 5630 a.u., exp(〈ln(θ )〉) = 1.04 a.u.

Again, the results of this simulation are in perfect agreement with the theoretical considerations we have developed previously.

c) Results in the model world:

The model parameter is now the random-walk FM level h -2 and the measurement data are its estimation ˆ.

h -2 Among the million simulated sequences, 3568 have a model parameter h -2 within the interval [0.952, 1.05]: model parameter, 〈h -2 〉 = 1.0007 a.u.; fit over classical estimates, 〈 〉 -ĥ 2 = 0.911 a.u. (-8.9%); and fit over log-unbiased estimates, 〈 〉 -ĥ 2 = 0.911 a.u. (-0.9%).

The estimation of the random-walk noise level is clearly better with the fit over the log-unbiased estimates, but the estimation of h -2 over the classical estimates is only slightly underestimated because it is based on all estimates, including those obtained for low τ values that are little biased (see Fig. 7).

2) White FM + Random Walk FM: We simulated 10 000 sequences of a combination of white FM and random-walk FM in such a way that the random-walk FM is the dominating noise over the last decade. For the simulations, we used model parameter: h 0 = 1 a.u., h -2 = 1.5 • 10 -5 a.u.; fit over classical estimates: 〈 〉 ĥ0 = 0.969 a.u. (-3.1%), 〈 〉 -ĥ 2 = 1.01 • 10 -5 a.u. (-33%); and fit over log-unbiased esti- Allan deviation (ADEV) curve for a pure random-walk FM with a unit noise level (obtained using SigmaTheta-1.2 [START_REF] Vernotte | Sig-maTheta software package[END_REF]). The crosses represent the ADEV estimates, whereas the circles represent the square root of the log-unbiased estimates. Fig. 8. Example of one realization among 10 000 simulations of an Allan deviation (ADEV) curve for a combination of white FM and randomwalk FM (obtained using SigmaTheta-1.2 [START_REF] Vernotte | Sig-maTheta software package[END_REF]). As in Fig. 1, the graph shows ADEV but the fit was performed over the Allan variance (AVAR) estimates. The crosses represent the ADEV estimates, whereas the circles represent the square root of the log-unbiased estimates. For comparison, Theo1 estimates [START_REF] Howe | dTheo1: Characterization of very long-term frequency stability[END_REF] are represented as asterisks.

y 2 ( 2 () is a measurement data of σ τ y 2 (

 222 = 10 s) = h 0 /20 (where h 0 is the white FM level) as a parameter; σ τ y = 10 s).

Fig. 2 .

 2 Fig. 2. Dispersion of Allan variance (AVAR) estimates for a pure white frequency noise and their equivalent degrees of freedom, ν (65 536 samples at 1 s).

Fig. 3 .

 3 Fig. 3. (a) Model parameter and measurement data for a χ 1 2 distribution (10 6 realizations of θ • χ 1 2 for θ = 10 k and -3 < k < +9). (b) Measurement dispersion for a parameter θ 0 ≈ 1. (c) Model parameter values for a measurement data ξ 0 ≈ 1.

  back into the model world for k = 1 and using (8):

  and E(X) = ν. The PDF is

Let us now consider the theoretical σ y 2 ( 512 )

 2512 as the model parameter and the σy 2 (512) computed value as the mea-

Fig. 4 .

 4 Fig. 4. Biases and 95% confidence intervals over the model parameter θ versus equivalent degrees of freedom, ν, in the measurement world where the measurement data are ξ 0 = 1. The circles represent the exponential of the mathematical expectation of the log of the model parameter and the crosses represent the mathematical expectation of the model parameter (which does not converge for ν ≤ 2).

Fig. 5 .

 5 Fig.5. Example of log-log fit (obtained with SigmaTheta-1.2[START_REF] Vernotte | Sig-maTheta software package[END_REF]). As in Fig.1, the graph shows the Allan deviation (ADEV) but the fit was performed over the Allan variance (AVAR) estimates. The crosses represent the ADEV estimates, whereas the circles represent the square root of the log-unbiased estimates.

Fig. 6 .

 6 Fig. 6. Comparison of the histogram of the simulations and the theoretical conditional probability distribution function (PDF) (a) in the model world, for which the model parameter is θ 0 = 1, and (b) in the measurement world, for which the measurement data are ξ 0 = 1.

Fig. 7 .

 7 Fig.7. Example of one realization among one million simulations of an Allan deviation (ADEV) curve for a pure random-walk FM with a unit noise level (obtained using SigmaTheta-1.2[START_REF] Vernotte | Sig-maTheta software package[END_REF]). The crosses represent the ADEV estimates, whereas the circles represent the square root of the log-unbiased estimates.

TABLE I .

 I T R   A V   D N T   P L M    L F D.

TABLE II .

 II B  95% C I B V EDF ν.

	ν	E(θ )/ξ 0	e E[ln(θ )] /ξ 0	B 2.5%	B 97.5%
	1	∞	3.56213	0.19905	1018.26086
	2	∞	1.78108	0.27108	39.49783
	3	3.00000	1.44625	0.32091	13.90209
	4	2.00000	1.31043	0.35896	8.25735
	5	1.66667	1.23754	0.38964	6.01530
	6	1.50000	1.19222	0.41524	4.84911
	7	1.40000	1.16137	0.43715	4.14232
	8	1.33333	1.13902	0.45624	3.67017
	9	1.28571	1.12210	0.47312	3.33285
	10	1.25000	1.10885	0.48821	3.07978

Let us define a new estimator Lθ as ˆln exp,

It is generally accepted that the frequency deviation samples follow statistics that can be approximated by a Gaussian model[START_REF] Lesage | Characterization of frequency stability: Uncertainty due to the finite number of measurements[END_REF], although the non-stationary effects of low-frequency noises could limit the relevance of this model.

This distribution represents, at best, the total ignorance of a value, i.e., a uniform law for the order of magnitude of this value. It was formalized as the Newcomb-Benford law[START_REF] Benford | The law of anomalous numbers[END_REF], or first-digit law, and in particular explains the reason for the uneven wear of logarithm tables[START_REF] Newcomb | Note on the frequency of use of the different digits in natural numbers[END_REF]. It is the basis of the concept of prior ignorance in Bayesian statistics[START_REF] Bernardo | Wiley Series in Probability and Statistics)[END_REF].

mates: 〈 〉 ĥ0 = 0.983 a.u. (-1.7%), 〈 〉 -ĥ 2 = 1.42 • 10 -5 a.u. (-5.5%).

In this case, the estimation of the random-walk FM noise level is drastically improved by the use of the logunbiased estimates because the random walk asymptote is fitted over the very last τ values, which are significantly affected by the log-log bias (see Fig. 8). For comparison, Theo1 estimates [START_REF] Howe | dTheo1: Characterization of very long-term frequency stability[END_REF] were added on Fig. 8.

IV. C

AVAR estimates are not biased and are easily fitted by weighted relative least squares. However, ln(AVAR) estimates are biased, and this log-log bias affects the fits performed on log-log plots.

A log-unbiased estimator has been defined in such a way that its log is an unbiased estimator of the log of the model parameter. On a log-log AVAR plot, it is then recommended to fit the log-unbiased estimates rather than the classical AVAR estimates.

For this purpose, we developed the SigmaTheta software package [START_REF] Vernotte | Sig-maTheta software package[END_REF] which performs fits according to the criterion (4) including log-unbiased estimator. It is free software developed under a GNU license for UNIX and Windows (Microsoft Corp., Redmond, WA) systems, available at http://theta.obs-besancon.fr/spip.php?article103.

R