Yves Cherel
email: cherel@cebc.cnrs.fr

Christophe Barbraud

Maxime Lahournat

Audrey Jaeger

Sébastien Jaquemet

Ross M Wanless

Richard T Phillips

David R Thompson

Paco Bustamante

Richard A Phillips

Accumulate or eliminate? Seasonal mercury dynamics in albatrosses, the most contaminated family of birds

Keywords: Antarctica, Feathers, Foraging, Methylmercury, Moult, Seabirds, Southern Ocean, Stable isotopes

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

An adult light-mantled sooty albatross at rest in a cliff on Possession Island, Crozet Archipelago, southern Indian Ocean

INTRODUCTION

Mercury (Hg) is a persistent and non-essential trace element that ranks third on the priority list of hazardous substances, based on toxicity and prevalence at contaminated sites (ATSDR; https://www.atsdr.cdc.gov/spl/). Hg is mobilized from geological deposits through both natural and anthropogenic processes and travels long distances within the atmosphere to reach even the most remote regions on Earth located far from emission sources. After atmospheric deposition, microorganisms methylate inorganic Hg (iHg) into methylmercury (MeHg), a powerful neurotoxin that bioaccumulates in organisms and biomagnifies in food webs to levels that pose major health risks to humans and wildlife [START_REF] Driscoll | Mercury as a global pollutant: sources, pathways, and effects[END_REF]. Among consumers, birds exhibit varying levels of trophic, spatial and temporal integration of contaminants and so are effective sentinels of MeHg bioavailability [START_REF] Furness | Birds as monitors of pollutants[END_REF][START_REF] Evers | Patterns and interpretation of mercury exposure in freshwater avian communities in northeastern North America[END_REF]. Since MeHg contamination increases from terrestrial to aquatic ecosystems, consumption of freshwater and marine foods constitutes the major sources of MeHg exposure in humans and animals [START_REF] Evers | Patterns and interpretation of mercury exposure in freshwater avian communities in northeastern North America[END_REF][START_REF] Driscoll | Mercury as a global pollutant: sources, pathways, and effects[END_REF]. Hence, seabirds have been frequently used as indicators of MeHg contamination in the marine environment [START_REF] Monteiro | Seabirds as monitors of mercury in the marine environment[END_REF][START_REF] Burger | Marine birds as sentinels of environmental pollution[END_REF].

Albatrosses are iconic pelagic seabirds with life-history traits that make them potentially at high risk of MeHg exposure. They are extremely long-lived animals with some individuals living >50 years [START_REF] Wasser | Avian longevities and their interpretation under evolutionary theories of senescence[END_REF], and they are at the top of the food web, typically preying on, or scavenging large carnivorous nekton [START_REF] Cherel | A review of the food of albatrosses[END_REF][START_REF] Cherel | Feeding ecology, isotopic niche, and ingestion of fishery-related items of the wandering albatross Diomedea exulans at Kerguelen and Crozet Islands[END_REF]. The liver is the principal tissue for long-term Hg storage [START_REF] Thompson | Metal levels in marine vertebrates[END_REF], and pioneering studies have shown that some albatross livers contained amongst the highest total Hg (THg) levels documented for apparently healthy free-living vertebrates, albeit lower than those recorded in some marine mammals [START_REF] Muirhead | Heavy metal concentrations in the tissues of seabirds from Gough Island, South Atlantic Ocean[END_REF][START_REF] Honda | High residue levels and detoxification processes of mercury in pelagic seabirds from the North Pacific Ocean[END_REF][START_REF] Stewart | Influence of phylogeny, diet, moult schedule and sex on heavy metal concentrations in New Zealand Procellariiformes[END_REF]. Subsequent measurements generally confirmed substantial Hg contamination of albatrosses, but with large variation depending on tissue type, individual, population and species [START_REF] Thompson | Temporal and spatial variation in mercury concentrations in some albatrosses and petrels from the sub-Antarctic[END_REF][START_REF] Kim | Mercury levels and its chemical form in tissues and organs of seabirds[END_REF][START_REF] Hindell | Mercury and cadmium concentrations in the tissues of three species of southern albatrosses[END_REF][START_REF] Stewart | Influence of phylogeny, diet, moult schedule and sex on heavy metal concentrations in New Zealand Procellariiformes[END_REF]. Comparison between contemporary and historical specimens showed slight or no increase in albatross THg levels over time, thus suggesting that contamination results primarily from natural -not anthropogenic -processes [START_REF] Thompson | Temporal and spatial variation in mercury concentrations in some albatrosses and petrels from the sub-Antarctic[END_REF], but see [START_REF] Vo | Temporal increase in organic mercury in an endangered pelagic seabird assessed by century-old museum specimens[END_REF].

The inference is that albatrosses are amongst the organisms most contaminated by Hg [START_REF] Thompson | Temporal and spatial variation in mercury concentrations in some albatrosses and petrels from the sub-Antarctic[END_REF], but the available data set remained largely incomplete and suffered from several limitations: (i) many taxa and populations have not been sampled, (ii) sample sizes were often low (n < 5), and (iii) relating published data to a species under the current taxonomy can be difficult because of splitting of species in the last decade [START_REF] Phillips | The conservation status and priorities for albatrosses and large petrels[END_REF]. This, together with the remoteness of most colonies, has resulted in a critical lack of comprehensive knowledge on Hg contamination in albatrosses.

Here, we report THg concentrations in 516 individual albatrosses from35 populations, including all 20 species and subspecies that breed in the Southern Ocean, thus providing new information on THg levels in 11 taxa and 22 populations (Table 1). THg concentration was measured in body feathers for practical, ethical and scientific reasons: (i) body feathers can be collected easily and noninvasively from live birds and they allow comparisons of metal exposure over time using museum specimens [START_REF] Thompson | Temporal and spatial variation in mercury concentrations in some albatrosses and petrels from the sub-Antarctic[END_REF][START_REF] Vo | Temporal increase in organic mercury in an endangered pelagic seabird assessed by century-old museum specimens[END_REF][START_REF] Carravieri | Penguins as bioindicators of mercury contamination in the southern Indian Ocean: geographical and temporal trends[END_REF], (ii) body feathers present less THg variation than do wing and tail feathers, and their collection does not impair flying ability [START_REF] Furness | Using bird feathers to measure mercury in the environment: relationships between mercury content and moult[END_REF][START_REF] Thompson | Temporal and spatial variation in mercury concentrations in some albatrosses and petrels from the sub-Antarctic[END_REF], (iii) most THg in feathers is organic, thus providing a mean of measuring MeHg exposure of birds (Thompson and Furness,1989b;[START_REF] Renedo | Combined impacts of longline fisheries and climate on the persistence of the Amsterdam albatross Diomedea amsterdamensis[END_REF], and (iv) plumage is a major pathway for MeHg elimination in avian species [START_REF] Burger | Metals in avian feathers: bioindicators of environmental pollution[END_REF][START_REF] Monteiro | Seabirds as monitors of mercury in the marine environment[END_REF]. We tested for the importance of several potential factors that might drive Hg exposure, including breeding frequency, geographical location of the breeding colonies and isotopic proxies of the foraging habitats and trophic levels during moult (feather d13C and d15N, respectively). In addition, we examined a potential latitudinal effect on Hg exposure in the Southern Ocean, because recent investigations indicate that birds foraging in colder, southerly waters present lower tissue Hg concentrations than those feeding further north in temperate waters (Carravieri et al., 2014a[START_REF] Carravieri | Penguins as bioindicators of mercury contamination in the southern Indian Ocean: geographical and temporal trends[END_REF][START_REF] Carravieri | From Antarctica to the subtropics: contrasted geographical concentrations of selenium, mercury, and persistent organic pollutants in skua chicks (Catharacta spp.)[END_REF].

MATERIALS AND METHODS

Study sites, birds and sampling

Body feathers were collected at each breeding site from 7 to 33 randomly chosen adult albatrosses over the period 2004-2013. Fieldwork was carried out on 14 islands and archipelagos that are scattered within the Southern Ocean (south of the Subtropical Front, STF) or in warmer fringing waters (Table S1). Two, five and seven islands were located in the Atlantic, Indian and Pacific oceans, respectively. The Subtropical (STZ), Subantarctic (SAZ) and Antarctic (AZ) Zones are here defined as the zones north of the STF, between the STF and the Polar Front (PF), and south of the PF, respectively (Fig. 1). The 14 sampling sites are located within the STZ (Amsterdam, Tasmanian and Chatham islands), the SAZ (Gough, Prince Edward, Crozet, Kerguelen, Auckland, Snares, Campbell, Antipodes and Bounty islands), and the AZ (South Georgia and Heard Island) (Table S1). Based on feather δ 13 C isoscapes (Jaeger et al., 2010), values of > -18.3‰, -21.2 to -18.3‰, and < -21.2‰ were considered to correspond to STZ, SAZ and AZ, respectively.

Moult in albatrosses and feather sampling

To test for the potential effects of foraging habitat and trophic position on feather THg concentrations, THg and isotopic values were measured on the same body feather taken from each individual albatross. In the Southern Ocean, d13C values of seabirds can be used to infer their foraging habitats [START_REF] Cherel | Geographical variation in carbon stable isotope signatures of marine predators: a tool to investigate their foraging areas in the Southern Ocean[END_REF][START_REF] Anderson | Influence of trophic position and foraging range on mercury levels within a seabird community[END_REF][START_REF] Quillfeldt | Moving polewards in winter: a recent change in the migratory strategy of a pelagic seabird?[END_REF] and δ 15 N values increase with trophic level [START_REF] Cherel | Isotopic niches and trophic levels of myctophid fishes and their predators in the Southern Ocean[END_REF]. Because keratin is a metabolically inactive molecule that is inert following synthesis, isotope values reflect the diet at the time of feather growth [START_REF] Hobson | Assessing avian diets using stable isotopes I: turnover of 13 C in tissues[END_REF][START_REF] Bearhop | Factors that influence assimilation rates and fractionation of nitrogen and carbon stable isotopes in avian blood and feathers[END_REF]. Hence, stable isotope analysis of feathers documents the feeding ecology of albatrosses during moult [START_REF] Cherel | Using stable-isotope analysis of feathers to distinguish moulting and breeding origins of seabirds[END_REF][START_REF] Anderson | Influence of trophic position and foraging range on mercury levels within a seabird community[END_REF]. Three potential limitations of the method are notable. (i) The exact timing of growth of body feathers in albatrosses is unknown as they are replaced gradually over the long inter-breeding period, with only ~7% of body feathers being moulted and regrown at any one time [START_REF] Battam | Chemical composition and tissue energy density of the cuttlefish (Sepia apama) and its assimilation efficiency by Diomedea albatrosses[END_REF]. The temporal window represented in the composition of body feathers depends on albatross breeding frequency and duration of the breeding cycle.

Consequently, the interbreeding (moulting) period spans a full year plus a winter (~16 months) in small biennial breeders (the grey-headed, sooty and light-mantled albatrosses), a full year (~12 months) in large biennial breeders of the genus Diomedea (seven taxa), and one winter (~4 months) in annual breeders (10 taxa) (Table S1). (ii) Moult of body feathers in albatrosses rarely occurs during the breeding period [START_REF] Catry | Differences in the timing and extent of annual moult of blackbrowed albatrosses Thalassarche melanophris living in contrasting environments[END_REF]. In the present study, most feather δ 13 C and δ 15 N values were different from the values that characterize feeding near the breeding sites [START_REF] Cherel | A comprehensive isotopic investigation of habitat preferences in nonbreeding albatrosses from the Southern Ocean[END_REF], thus verifying that feather moult took place during the inter-breeding period. (iii) More importantly, there is a temporal mismatch between isotopic values and THg concentrations in feathers of adult birds. The latter is considered to represent Hg accumulated between two consecutive moulting cycles, i.e. a much longer integration period than that corresponding to feather growth [START_REF] Bond | Relationships between stable isotopes and metal contaminants in feathers are spurious and biologically uninformative[END_REF]. However, a preliminary analysis of our data indicated a more complex picture; feather THg concentrations were correlated with δ 13 C in feathers, which represents the carbon source in moulting habitat (see Results). The very low THg concentrations in feathers that were grown in Antarctic waters (but not further north) prompted further investigations using three albatross populations that partly moult in Antarctic waters; grey-headed, sooty and light-mantled albatrosses from the Prince Edward, Gough, and Kerguelen Islands, respectively [START_REF] Jaeger | Stable isotopes document inter-and intra-specific variation in feeding ecology of nine large southern Procellariiformes[END_REF][START_REF] Cherel | A comprehensive isotopic investigation of habitat preferences in nonbreeding albatrosses from the Southern Ocean[END_REF]. In those birds, THg concentrations and δ 13 C and δ 15 N values were measured on three additional body feathers from each individual (for a total of four feathers per bird) to better investigate potential relationships between Hg exposure and the foraging ecology of albatrosses during the inter-breeding period.

Mercury and stable isotope analyses

THg concentrations, δ 13 C and δ 15 N values were measured on the same individual body feathers.

Feathers were first cleaned to remove surface contaminants using a 2:1 chloroform:methanol solution followed by two methanol rinses. They were then oven dried for 48 h at 50 °C. Feather samples (1-5mg dry weight, dw) were analyzed for THg with an Advanced Mercury Analyzer spectrophotometer (Altec AMA 254) [START_REF] Bustamante | Total and organic Hg concentrations in cephalopods from the North East Atlantic waters: influence of geographical origin and feeding ecology[END_REF]. The analyses were repeated in duplicate or triplicate until the relative standard deviation was <10% for multiple aliquots. Tort-2 Lobster Hepatopancreas, NRC, Canada (THg concentration: 0.27 ± 0.06 µg/g dw) was used as Certified Reference Material (CRM) to check the accuracy of the method. After adjusting the CRM mass to represent the same amount of THg in reference samples as in those from feathers, our measured values for CRM averaged 0.273 ± 0.008 µg/g dw (n = 53), thus corresponding to a recovery rate of 101 ± 3%. Blanks were analyzed at the beginning of each set of samples, and the detection limit of the method was 0.005 µg/g dw. Feather THg concentrations are presented relative to the dry weight.

For isotopic analyses, ~0.4 mg subsamples of the feather homogenates were weighed in tin cups. A continuous flow mass spectrometer (Thermo Scientific Delta V Advantage) was coupled to an elemental analyzer (Thermo Scientific Flash EA 1112) to measure feather δ 13 C and δ 15 N values, and feather carbon and nitrogen contents, respectively. Stable isotope ratios are expressed using standard d notation relative to carbonate Vienna PeeDee Belemnite and atmospheric nitrogen. The internal laboratory standards were acetanilide. Observed analytical errors were <0.10‰ for both δ 13 C and δ 15 N values. The data on feather THg concentrations of albatrosses from the Kerguelen Islands, and most of the data presented here on feather δ 13 C and δ 15 N values are already published, but in distinct scientific contexts [START_REF] Cherel | A comprehensive isotopic investigation of habitat preferences in nonbreeding albatrosses from the Southern Ocean[END_REF]Carravieri et al., 2014b). Values are means ± SD.

Statistical analyses

The effects of species, phylogeny (Diomedea, large Thalassarche, small Thalassarche, Phoebetria), breeding frequency, breeding location, oceanographic zone (Table S1), ocean basin (Indian, Pacific, Atlantic), and foraging habitat and trophic level during moult (δ 13 C and δ 15 N, respectively) on feather THg concentrations (516 feathers from 516 individuals from 35 populations) were investigated using non-parametric smoothing regression techniques. While body condition is a main driving factor of avian THg concentrations (Fort et al., 2015), its effect was not tested because moult occurs at sea, thus precluding handling albatrosses at that time.

We used Generalised Additive Mixed Models (GAMM) in order to model potential non-linear relationships between THg concentrations and covariates, and to account for pseudoreplication. Indeed, our sampling unit was the population (n = 35, Table S1), represented by several individuals (mean of 15 ± 6, range 7-33). Population was thus modelled as a random effect in the GAMM. Data exploration was first carried out following [START_REF] Zuur | Mixed effects models and extensions in ecology with R[END_REF]. We then used variance inflation factors to assess which explanatory variables were collinear and should be dropped, based on a cut-off value of 4 [START_REF] Zuur | Mixed effects models and extensions in ecology with R[END_REF]. The following was then considered as a starting model: THg ~ s(δ 13 C) + s(δ 15 N) + species + phylogeny + breeding frequency + breeding location + oceanographic zone + ocean basin + two by two interactions between all categorical variables + random(population).

We then searched for the optimal model in terms of the explanatory variables and their interactions using a top-down strategy removing terms with non-significant p-values (significance level of  < 0.01). The GAMM models ('gamm4' library, version 0.2-4) were specified with a Gaussian family and a penalized thin plate regression spline, the optimal span for smoothing was determined by AIC, and the smoothing parameters estimated using restricted maximum likelihood [START_REF] Wood | Thin plate regression splines[END_REF][START_REF] Zuur | Mixed effects models and extensions in ecology with R[END_REF]. The adjusted R-squared for the model was defined as the proportion of variance explained, where original variance and residual variance were both estimated using unbiased estimators. Model validation was performed by inspecting the relationship between residuals and fitted values [START_REF] Zuur | Mixed effects models and extensions in ecology with R[END_REF]. GAMM analyses were performed using R 3.3.2 (R Development Core Team, 2016) and other statistical analyses using SYSTAT 13 for WINDOWS (Systat Software, Chicago).

RESULTS

THg concentrations varied 49-fold within the whole data set (618 body feathers), with the lowest value in a light-mantled albatross from the Kerguelen Islands and the highest in an Amsterdam albatross from Amsterdam Island (1.4 and 70.6 µg/g dw, respectively). All measurements of feather THg were >1 µg/g, of which 74.1% and 1.1% were >5 µg/g and >40 µg/g, respectively, which are the two most widely used toxicity thresholds for feather THg [START_REF] Burger | Risk, mercury levels, and birds: relating adverse laboratory effects to field biomonitoring[END_REF][START_REF] Evers | Adverse effects from environmental mercury loads on breeding common loons[END_REF]. Notably, all of the seven most THgcontaminated feathers (>40 mg/g) were from Diomedea albatrosses, including four Amsterdam albatrosses. At the population level (n = 35), THg concentrations varied nine-fold, with mean values ranging from 3.8 µg/g in Atlantic yellow-nosed albatross to 34.6 µg/g in Amsterdam albatross (Fig. 2, Table S1). Accordingly, feather THg concentrations differed significantly between populations (Kruskal-Wallis, H = 353.3, p < 0.0001).

Results from the optimal GAMM model are presented in Table 2. Of the fixed factors, species, breeding frequency, breeding location, oceanographic zone and ocean basin were not retained, but phylogeny and feather δ 13 C and δ 15 N values had significant effects on feather THg concentrations. THg concentrations differed between albatross taxa, from 6.4 ± 4.4 µg/g in small Thalassarche species to 21.9 ± 11.0 µg/g in Diomedea spp. (Fig. 3). Feather THg concentrations were positively related to feather δ 13 C values (Fig. 4) and the estimated degree of freedom of the smooth term (2.74) indicated a near-linear effect. In contrast, the relationship between feather δ 15 N and THg concentrations was clearly non-linear (Fig. 4). The overall deviance explained by the optimal model was 49.5%. The population random effect explained 48.9% of the total variance of the residuals. There was no pattern between the residuals and the fitted values.

As expected, feather δ 13 C values ranged widely in the three populations for which isotopic values and THg concentrations were determined on four body feathers per bird. Grey-headed, sooty and light-mantled albatrosses moulted over a large latitudinal gradient from Antarctica to the subtropics, with considerable variation both within and between individuals (Table 3, Fig. 5). Two features are notable: (i) all feathers that were grown in the AZ contained less THg than those moulted and regrown further north, and (ii) THg population concentrations of feathers grown in the AZ remained low over a wide range of δ 15 N values. For example, feathers of lightmantled albatrosses contained 8.8 times more THg if they were grown in the STZ than in the AZ, with no overlap in feather THg concentrations between the two groups. There was no change in THg concentration in feathers grown in the AZ, while their δ 15 N values ranged from 7.3 to 11.7‰ (a 4.4‰ difference). When feathers grew in warmer waters, THg concentrations were high, ranged widely and were associated with higher δ 15 N values (Fig. 5). Conversely, feathers of grey-headed albatrosses that were grown south and north of the PF showed no significant difference in mean δ 15 N values (n = 8 and 36, 11.0 ± 1.1 and 11.8 ± 1.5‰, respectively, Mann-Whitney test, U = 86.0, p = 0.078), but those grown in the AZ contained 2.4 times less THg than those grown further north (4.5 ± 1.2 and 10.9 ± 6.4 µg/g; U = 22.0, p < 0.0001).

DISCUSSION

Our key finding is that albatrosses are the most contaminated by Hg of all families of birds. The present work provides extensive new data on MeHg exposure in albatrosses from the Southern Ocean, where most species breed, thus allowing a comprehensive synthesis of the available information on albatrosses worldwide (Table 1). Two features are notable: (i) THg measurements are urgently needed from two Pacific Phoebastria species, waved albatross and short-tailed albatross (but see [START_REF] Shinsuke | Trace element accumulation in shorttailed albatrosses (Diomedea albatrus) and black-footed albatrosses (Diomedea nigripes) from Torishima Island, Japan[END_REF], and (ii) seven (32%) of the albatross species and subspecies present mean body feather THg concentrations > 20 mg/g, including all five taxa of the wandering albatross species complex (wandering, Antipodean, Gibson's, Amsterdam and Tristan albatrosses), sooty albatross, and black-footed albatross. A review of THg body feathers emphasizes the exceptionally high contamination levels of albatrosses when compared to adult seabirds from the more investigated Northern Atlantic and Northern Pacific (>35°N; from 0.6 to 22.6 µg/g, and from 0.9 to 10.3 µg/g, respectively) [START_REF] Kim | Mercury levels and its chemical form in tissues and organs of seabirds[END_REF][START_REF] Monteiro | Contemporary patterns of mercury contamination in the Portuguese Atlantic inferred from mercury concentrations in seabird tissues[END_REF]Bond and Diamond, 2008;Eagles-Smith et al., 2009). Only four populations of Procellariidae (closely-related to Diomedeidae) showed THg levels > 20 µg/g, northern giant petrel Macronectes halli and great-winged petrel Pterodroma macroptera from Prince Edward Islands, Atlantic petrel P. incerta from Gough, and Bulwer's petrel Bulweria bulwerii from the Azores, Madeira and Salvages [START_REF] Monteiro | Contemporary patterns of mercury contamination in the Portuguese Atlantic inferred from mercury concentrations in seabird tissues[END_REF][START_REF] Becker | Feather mercury concentrations in Southern Ocean seabirds: variation by species, site and time[END_REF]. In contrast, average THg body feather levels >20 µg/g are very rare in other avian species, occurring only in a few diurnal and nocturnal raptors from the more contaminated Northern Hemisphere [START_REF] Burger | Metals in avian feathers: bioindicators of environmental pollution[END_REF][START_REF] Bowerman | Using feathers to assess risk of mercury and selenium to bald eagle reproduction in the Great Lakes Region[END_REF], and in tern chicks in two highly polluted areas [START_REF] Furness | Spatial and temporal variation in mercury contamination of seabirds in the North Sea[END_REF][START_REF] Herring | Mercury exposure may suppress baseline corticosterone levels in juvenile birds[END_REF].

Phylogeny explained much of the variation in feather THg concentrations in albatrosses;

Diomedea species ranked first, and their body feathers contained on average 3.4 more THg than smaller Thalassarche species (Fig. 3). The wandering albatross species complex were the most Hg contaminated taxa, whereas northern and southern royal albatrosses showed moderate THg levels. These two groups of Diomedea albatrosses have essentially the same life history traits, differing ecologically only in their main feeding grounds. Royal albatrosses forage primarily in neritic waters, whereas the wandering albatross and related taxa are essentially oceanic (Nicholls et al., 2002;Waugh and Weimerskirch, 2003), which fits well with a pattern of increasing feather THg values from coastal to oceanic seabirds [START_REF] Ochoa-Acuña | Mercury in feathers from Chilean birds: influence of location, feeding strategy, and taxonomic affiliation[END_REF].

Phoebetria albatrosses ranked second, with the differences between the two species in feather THg concentrations likely explained by contrasting distributions; light-mantled albatrosses partly moult in the AZ, whereas sooty albatrosses do so further north (see below). Finally, the lowest feather THg concentrations were found in Thalassarche albatrosses; these are small to medium-sized albatrosses that are mostly annual breeders and, because they tend to eat more crustaceans, do not feed at as high trophic levels as the Diomedea albatrosses (Prince et al., 1993;[START_REF] Cherel | A review of the food of albatrosses[END_REF][START_REF] Jaeger | Stable isotopes document inter-and intra-specific variation in feeding ecology of nine large southern Procellariiformes[END_REF]. They thus differ from Diomedea and Phoebetria albatrosses in several life-history traits, which likely explain the lower THg concentrations in their feathers (this study), and also in internal tissues [START_REF] Stewart | Influence of phylogeny, diet, moult schedule and sex on heavy metal concentrations in New Zealand Procellariiformes[END_REF].

Feather δ 15 N was also related to THg concentrations, which was to some extent expected given δ 15 N values is a proxy of trophic level and hence linked to biomagnification processes within ecosystems (the higher the trophic position, the higher the MeHg exposure). However, the relationship was nonlinear (Fig. 4), indicating a more complex relationship involving other confounding factors, including varying δ 15 N baselines at moulting grounds. The large range of feather δ 13 C values showed that plumage replacement took place in habitats that varied from Antarctic oceanic waters to northern coastal upwelling zones, depending on species, population and individual [START_REF] Cherel | A comprehensive isotopic investigation of habitat preferences in nonbreeding albatrosses from the Southern Ocean[END_REF]. The different moulting habitats are likely marked by contrasting δ 15 N baselines that propagate up the food webs and thus raw δ 15 N values do not directly indicate relative trophic level. For example, the δ 15 N baseline varies little within the Southern Ocean but increases sharply at the STF, remaining elevated further north [START_REF] Altabet | Sedimentary nitrogen isotopic ratio as a recorder for surface ocean nitrate utilization[END_REF], and the δ 15 N baseline varies significantly within the Humboldt Current where many albatrosses moult [START_REF] Mollier-Vogel | Nitrogen isotope gradients off Peru and Ecuador related to upwelling, productivity, nutrient uptake and oxygen deficiency[END_REF][START_REF] Cherel | A comprehensive isotopic investigation of habitat preferences in nonbreeding albatrosses from the Southern Ocean[END_REF].

Another key finding of the present study was the major influence of feather δ 13 C values, which reflect moulting habitat, on feather THg concentration in albatrosses. Conversely, breeding locality, oceanographic zone and ocean basin during breeding did not show significant relationships with feather THg concentration. This is presumably related largely to the timing of body feather moult, which takes place during the inter-breeding period. Importantly, our results indicate that in albatrosses, body feather THg concentration integrates predominantly very short-term MeHg intake at the temporal scale of a single growing body feather (Fig. 5).

Consequently, there is no temporal mismatch between feather stable isotopic values and THg concentrations, all three parameters reflecting the feeding habits of albatrosses during feather growth. Hence, unlike feathers from most birds [START_REF] Furness | Using bird feathers to measure mercury in the environment: relationships between mercury content and moult[END_REF][START_REF] Braune | Mercury levels in Bonaparte's gulls (Larus philadelphia) during autumn molt in the Quoddy Region, New Brunswick, Canada[END_REF][START_REF] Bond | Relationships between stable isotopes and metal contaminants in feathers are spurious and biologically uninformative[END_REF], body feathers of albatrosses carry information about MeHg availability on the moulting grounds, with differences in the carbon source reflected in δ 13 C values. Furthermore, it follows that because albatross body feather THg concentrations reflect MeHg intake over the time of feather growth, feather THg concentrations do not fully reflect MeHg accumulation over the long-term. Thus, plumage renewal in albatrosses does not appear to be a mechanism by which much of MeHg accumulated during the inter-moult period can be eliminated, unlike many other bird species that excrete most ingested and stored MeHg into newly grown feathers [START_REF] Burger | Metals in avian feathers: bioindicators of environmental pollution[END_REF][START_REF] Monteiro | Seabirds as monitors of mercury in the marine environment[END_REF]. Consequently, a critical issue for albatrosses is the fate of ingested MeHg, specifically in order to cope with potentially toxic effects. Since albatrosses are known to efficiently demethylate MeHg, we propose that they use two complementary detoxification mechanisms: (i) during moult, dietary MeHg is primarily eliminated directly into the keratin of growing feathers, and (ii) outside the moulting period, MeHg is continuously converted to iHg in the liver where it accumulates. These two mechanisms explain: (i) why feather THg does not increase with age in albatrosses [START_REF] Thompson | Temporal and spatial variation in mercury concentrations in some albatrosses and petrels from the sub-Antarctic[END_REF]Tavares et al., 2013;[START_REF] Bustamante | High feather mercury concentrations in the wandering albatross are related to sex, breeding status and trophic ecology with no demographic consequences[END_REF], (ii) the records of exceptionally high liver THg concentrations, with a large predominance of iHg over MeHg, in some (presumably old) individual albatrosses [START_REF] Honda | High residue levels and detoxification processes of mercury in pelagic seabirds from the North Pacific Ocean[END_REF]Thompson and Furness, 1989a;[START_REF] Kim | Mercury levels and its chemical form in tissues and organs of seabirds[END_REF][START_REF] Stewart | Influence of phylogeny, diet, moult schedule and sex on heavy metal concentrations in New Zealand Procellariiformes[END_REF], and (iii) the much lower THg burden in feathers (<10%) than in the whole body, indicating that THg excretion during the moult is negligible in albatrosses [START_REF] Honda | High residue levels and detoxification processes of mercury in pelagic seabirds from the North Pacific Ocean[END_REF][START_REF] Kim | Mercury levels and its chemical form in tissues and organs of seabirds[END_REF].

We conclude that two strategies to deal with ingested MeHg operate in birds and we propose to call them the income and capital strategies in reference to the previous definition of those terms in the management of energy for breeding [START_REF] Jonsson | Capital and income breeding as alternative tactics of resource use in reproduction[END_REF]. Birds adopting the income strategy (e.g. albatrosses) convert ingested MeHg to iHg in the liver outside the moulting period and only eliminate MeHg ingested during feather growth into newly growing feathers. In species using the capital strategy (e.g. gulls; [START_REF] Braune | Mercury levels in Bonaparte's gulls (Larus philadelphia) during autumn molt in the Quoddy Region, New Brunswick, Canada[END_REF], most feather MeHg comes from the accumulation of dietary MeHg into the body reservoir during the period between successive moults. A consequence of the capital strategy is that, unlike in albatrosses, feathers grown by these other species contain a major part (56-94%) of the THg burden of the whole body, indicating that they are the endpoint in a primary pathway for THg elimination [START_REF] Burger | Metals in avian feathers: bioindicators of environmental pollution[END_REF][START_REF] Monteiro | Seabirds as monitors of mercury in the marine environment[END_REF]. The concept of income and capital strategies requires further investigation, especially because some birds probably fall somewhere along an income-capital continuum. Albatrosses stands at the income end of the spectrum, and the income strategy is likely to be the prominent strategy in seabirds for which the moulting habitat affects feather THg concentrations [START_REF] Fort | Spatial ecotoxicology: migratory Arctic seabirds are exposed to mercury contamination while overwintering in the northwest Atlantic[END_REF][START_REF] Watanuki | Feather mercury concentration in streaked shearwaters wintering in separate areas of southeast Asia[END_REF].

Overall, feather d13C values were near-linearly and positively related to feather THg concentrations (Fig. 4). Hence, Hg exposure is higher for albatrosses moulting in warmer, more northerly, than colder, more southerly waters, because both baseline and consumer d13C values vary with latitude; the higher the latitude, the lower the δ 13 C values [START_REF] Cherel | Geographical variation in carbon stable isotope signatures of marine predators: a tool to investigate their foraging areas in the Southern Ocean[END_REF][START_REF] Anderson | Influence of trophic position and foraging range on mercury levels within a seabird community[END_REF][START_REF] Jaeger | Seabird satellite tracking validates the use of latitudinal isoscapes to depict predators' foraging areas in the Southern Ocean[END_REF][START_REF] Quillfeldt | Moving polewards in winter: a recent change in the migratory strategy of a pelagic seabird?[END_REF]. This relationship held both between and within individual birds. While all body feathers of grey headed, sooty and lightmantled albatrosses grown in the STZ and the SAZ had elevated THg levels, feathers grown in the AZ showed consistently low THg concentrations (Fig. 5). Surprisingly, feather THg concentrations remained low whatever their corresponding δ 15 N values, indicating that habitat is much more important than trophic level when moulting in the AZ. The large feather δ 15 N range suggests that albatrosses fed on prey from low trophic levels (probably including Antarctic krill Euphausia superba) in southern Antarctic waters to high trophic levels (fish and squid) south of the PF, but with no concomitant feather THg increase. The primary importance of feeding latitude on Hg levels is highlighted by intra and inter-individual variation in moulting grounds. Some individuals (e.g. light-mantled albatrosses shown in green and pink on Fig. 5) replaced their plumage in various oceanic zones, resulting in contrasting THg concentrations.

Other individuals moulted over a more restricted latitudinal range, and consequently, levels of Hg contamination of birds moulting in the AZ were approximately ten times lower than in those moulting further north (e.g. the cyan vs. red individuals, 2.1 vs. 19.9 µg/g; Fig. 5).

Our data confirm that exposure to MeHg of top predators decreased with increasing latitude in the Southern Ocean (Carravieri et al., 2014a[START_REF] Carravieri | Penguins as bioindicators of mercury contamination in the southern Indian Ocean: geographical and temporal trends[END_REF][START_REF] Carravieri | From Antarctica to the subtropics: contrasted geographical concentrations of selenium, mercury, and persistent organic pollutants in skua chicks (Catharacta spp.)[END_REF]. Moreover, they show a stepwise drop in MeHg availability at the PF, thus suggesting that the PF acts as an oceanographic barrier [START_REF] Guynn | Mercury concentrations in the Patagonian toothfish, Dissostichus eleginoides Smitt 1898, among three distinct stocks[END_REF]. Lower Hg contamination in predators feeding in the AZ in our study contrasts with earlier work on Hg in the Southern Ocean, which showed higher MeHg concentrations available at shallower depths in waters to the south, compared to the north, of the PF [START_REF] Cossa | Mercury in the Southern Ocean[END_REF]. Two non-exclusive explanations for this apparent mismatch between the base and the top of the trophic web merit further investigations: (i) a more complex oceanic food web involving more trophic levels and thus higher MeHg biomagnification processes in warmer (northern) than colder (southern) waters, and (ii) a higher MeHg enrichment within the trophic web due to higher production of MeHg by pelagic organisms [START_REF] Pucko | Transformation of mercury at the bottom of the Arctic food web: an overlooked puzzle in the mercury exposure narrative[END_REF] to the north of the PF.

Conclusions

Overall, albatrosses clearly illustrate how elevated THg levels in avian species are governed by a combination of several factors, specifically (i) foraging grounds marked by high MeHg bioavailability (temperate waters) and (ii) feeding at the top of trophic webs (maximal biomagnification). The Diomedeidae is among the world's most threatened families of birds, will all the 22 species classified from Near Threatened to Critically Endangered in the IUCN Red List (Table S1). They face a number of direct and indirect threats that are causes of population declines, including climate change, fishery-related mortality, infectious diseases, introduced predators, and pollutants [START_REF] Phillips | The conservation status and priorities for albatrosses and large petrels[END_REF]. In this context, it is noteworthy that two of the three Critically Endangered species were also the most contaminated by Hg (Amsterdam and Tristan albatrosses), and that information is lacking on the third species (waved albatross). It is especially relevant that most individuals with THg levels above the upper feather toxicity threshold (40 µg/g; [START_REF] Burger | Risk, mercury levels, and birds: relating adverse laboratory effects to field biomonitoring[END_REF][START_REF] Evers | Adverse effects from environmental mercury loads on breeding common loons[END_REF] were breeding Amsterdam albatrosses, the rarest albatross on Earth with a single population of approximately 165 adults (C. Barbraud et al., unpublished data). The whole population of Amsterdam albatrosses is at risk, because chicks are also heavily contaminated, with body feather concentrations averaging 7.8 µg/g (Y. Cherel et al., unpublished data).

One of the main problems in environmental toxicology is to interpret the impact of observed contaminant concentrations on wildlife. Most adult albatrosses have THg levels above the lower feather toxicity threshold (5 µg/g; Fig. 2), but seabirds are expected to have higher toxicity resistance than other avian taxa since they cope with relatively high MeHg ingestion through efficient hepatic demethylation processes (Thompson and Furness, 1989a;[START_REF] Ikemoto | Detoxification mechanism of heavy metals in marine mammals and seabirds: interaction of selenium with mercury, silver, copper, zinc, and cadmium in liver[END_REF]. However, pioneering investigations on the biological consequences of Hg contamination pointed out (i) negative effects of blood THg on the immune function of black-footed albatrosses [START_REF] Finkelstein | Contaminantassociated alteration of immune function in black-footed albatross (Phoebastria nigripes), a North Pacific predator[END_REF], and (ii) negative effects of blood THg on demographic parameters of wandering albatrosses [START_REF] Goutte | Demographic consequences of heavy metals and persistent organic pollutants in a vulnerable long-lived bird, the wandering albatross[END_REF]. These data emphasize the urgent need for robust assessment of the impact of Hg exposure on the physiology and demography of albatrosses at a species, population and individual level, because Hg loading to the World Ocean is likely to continue rising [START_REF] Lamborg | Mercury in the anthropocene ocean[END_REF], thus driving an increase accumulation of Hg in the biota, including albatrosses [START_REF] Vo | Temporal increase in organic mercury in an endangered pelagic seabird assessed by century-old museum specimens[END_REF][START_REF] Becker | Feather mercury concentrations in Southern Ocean seabirds: variation by species, site and time[END_REF]. [START_REF] Burger | Risk, mercury levels, and birds: relating adverse laboratory effects to field biomonitoring[END_REF][START_REF] Evers | Adverse effects from environmental mercury loads on breeding common loons[END_REF]. Colour and symbols: red circles (Diomedea albatrosses), grey triangles up (large Thalassarche spp.), blue triangles down (small Thalassarche spp.), green diamonds (Phoebetria spp.). Abbreviations are detailed in Table S1.

Values are means ± SD. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.). Table 1. Review of total mercury (THg) concentrations in body feathers of adult albatrosses. All the existing literature was reviewed up to date (December 2017). A few data were deleted due to various problems, including (i) too small sample sizes (n <5; [START_REF] Ochoa-Acuña | Mercury in feathers from Chilean birds: influence of location, feeding strategy, and taxonomic affiliation[END_REF][START_REF] Shinsuke | Trace element accumulation in shorttailed albatrosses (Diomedea albatrus) and black-footed albatrosses (Diomedea nigripes) from Torishima Island, Japan[END_REF], (ii) taxonomic uncertainties (e.g. the Diomedea exulans species complex; Thompson & Furness (1989b), [START_REF] Thompson | Temporal and spatial variation in mercury concentrations in some albatrosses and petrels from the sub-Antarctic[END_REF], and (iii) geographical uncertainties together with other feather types than body feathers (e.g. North Pacific and flight feathers; [START_REF] Honda | High residue levels and detoxification processes of mercury in pelagic seabirds from the North Pacific Ocean[END_REF][START_REF] Kim | Mercury levels and its chemical form in tissues and organs of seabirds[END_REF]

Figure 1 .

 1 Figure 1. Location of albatross breeding localities and of the main oceanic fronts and zones of the Southern Ocean and fringing waters. Abbreviations: AZ, Antarctic Zone; PF, Polar Front; SAZ, Subantarctic Zone; STF, Subtropical Front; STZ, Subtropical Zone.

Figure 2 .

 2 Figure 2. Concentrations of total mercury (THg) in body feathers of albatrosses from the Southern Ocean and fringing subtropical waters. Albatross taxa were deliberately placed in a sequence of increasing feather THg concentrations. Limited data for THg concentrations in feathers indicate a range of 5 to 40 µg/g (dotted lines) as being associated with adverse effects in birds[START_REF] Burger | Risk, mercury levels, and birds: relating adverse laboratory effects to field biomonitoring[END_REF][START_REF] Evers | Adverse effects from environmental mercury loads on breeding common loons[END_REF]. Colour and symbols: red circles

Figure 3 .

 3 Figure 3. Histograms of feather THg concentrations in relation to albatross phylogeny. Kruskal-Wallis, H = 219.89, p <0.001; superscript letters indicate results from Conover-Inman tests for pairwise comparisons, all p <0.0001; see also text and Table2for GAMM statistics.Values are means ± SD.

Figure 4 .

 4 Figure 4. Fitted GAMM results showing the relationship between feather THg concentration and feather δ 13 C and δ 15 N values for the 516 sampled albatrosses. The y-axis gives the smooth transformed value, and the little vertical lines along the x-axis indicate the feather THg values of the observations. Dashed lines indicate the 95% confidence intervals.

Fig. 5 .

 5 Fig. 5. Total mercury (THg) versus δ 13 C (a proxy of moulting habitat) and δ 15 N (a proxy of trophic position) in feathers of light-mantled sooty albatrosses from Kerguelen Islands. THg concentrations and δ 13 C and δ 15 N values were measured in four different body feathers per individual (one different colour per bird). Following Jaeger et al. (2010), dashed lines correspond to the isotopic estimation of the Polar Front (PF) and of the Subtropical Front (STF), which delimit the Antarctic (AZ), Subantarctic (SAZ) and Subtropical (STZ) Zones; LMSA, light-mantled albatross. Diamonds and circles indicate feathers that grew in the Antarctic Zone and further north, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Table 2 .

 2 . Values are means ± SD with ranges in parentheses. GAMM results for feather THg concentrations as a function of covariates for the optimal model (516 sampled individuals from 35 populations). Variance components plus their SD are shown for random effects; edf indicates the estimated degrees of freedom for the smooth terms.

	Taxa	Location	Individuals (n)	THg (µg/g)	References
	Diomedea spp.
	Wandering albatross	South Georgia	66	19.59 ± 10.12 (5.44 -73.42)	Thompson et al. (1993)
			14	27.43 ± 8.14 (15.40 -45.36)	Anderson et al. (2009)
			34	20.14 ± 7.64 (8.09 -39.82)	Tavares et al. (2013)
			10	20.31 ± 5.88 (14.69 -29.97)	This study
		Prince Edward	29	24.83 ± 12.35 (11.84 -58.58)	Thompson et al. (1993)
			12	24.80 ± 8.61 (11.87-41.95)	This study
		Crozet	201	23.87 ± 11.83. (5.86 -94.72)	Bustamante et al. (2016)
			12	23.75 ± 6.84 (13.79-38.30)	This study
		Kerguelen	12	16.59 ± 3.78 (9.84 -24.13) Carraviéri et al. (2014b), this study
	Antipodean albatross	Antipodes	15	21.87 ± 9.27 (9.22 -41.93)	This study
	Gibson's albatross	Auckland	15	29.50 ± 8.17 (16.85 -50.09)	This study
	Amsterdam albatross	Amsterdam	18	34.60 ± 12.50 (18.57-70.63)	This study
	Tristan albatross	Gough	59	30.2 ± 15.6	Thompson et al. (1990)
			27	28.00 ± 14.31 (10.80 -65.50)	Thompson et al. (1993)
			12	24.96 ± 6.85 (15.74-35.74)	This study
	Southern royal albatross	Campbell	22	11.52 ± 13.86 (1.95 -45.46)	Thompson et al. (1993)
		Auckland	15	11.78 ± 8.14 (2.33 -27.57)	This study
	Northern royal albatross	Chatham	20	11.40 ± 6.46 (3.50 -26.27)	This study
	Phoebastria spp.
	Black-footed albatross	Midway	17	19.60 ± 1.75	Burger & Gochfeld (2000a,b)
	Laysan albatross	Midway	13	3.46 ± 0.39	Burger & Gochfeld (2000a,b)
			10	1.86 ± 0.35	Burger & Gochfeld (2000c)
	Large Thalassarche spp.
	Shy albatross	Tasmania	8	7.42 ± 2.82 (5.26-13.12)	This study
	White-capped albatross	Auckland	31	11.99 ± 4.67 (3.87 -20.80)	This study
	Salvin's albatross	Snares	15	5.59 ± 3.36 (2.04-15.60)	This study
		Bounty	16	7.39 ± 4.96 (2.68-22.16)	This study
	Chatham albatross	Chatham	11	11.46 ± 6.65 (2.32 -23.10)	This study
	Small Thalassarche spp.
	Black-browed albatross	Falklands	30	2.68 ± 1.22 (1.00 -6.64)	Thompson et al. (1993)
		South Georgia	20	4.57 ± 1.98 (1.41 -9.14)	Thompson et al. (1993)
			16	5.39 ± 2.05	Becker et al. (2002)
			16	8.35 ± 2.63 (4.24 -12.97)	Anderson et al. (2009)
			10	6.86 ± 2.87 (4.31 -12.76)	This study
		Kerguelen	33	4.07 ± 1.78 (1.71 -7.75)	Carraviéri et al. (2014b), this study
		Heard	8	9.93 ± 6.61 (2.60-19.63)	This study
	Campbell albatross	Campbell	35	10.06 ± 4.40 (3.52 -19.54)	Thompson et al. (1993)
			20	12.03 ± 5.33 (3.59 -22.31)	This study

Table 3 .

 3 Feather δ 13 C and δ 15 N values and total mercury (THg) concentrations in relation to albatross moulting zones (according to the estimated isotopic positions of oceanic fronts in Jaeger et al. 2010). Stables isotopes and THg were measured on four different body feathers per individual bird. Kruskal-Wallis H tests were performed to compare moulting zones and species. Significant statistics (p <0.05) are highlighted in bold. Values are means ± SD.

		Grey-headed albatross	Sooty albatross	Light-mantled sooty	Kruskal-Wallis
		Prince Edward	Gough	albatross, Kerguelen	H tests
		(n = 11)	(n = 13)	(n = 10)
	Antarctic Zone
	Feathers (n)	8	9	13
	Feather δ 13 C (‰)	-22.2 ± 0.6	-22.5 ± 0.4	-25.3 ± 1.5	H = 17.3, p <0.0001
	Feather δ 15 N (‰)	11.0 ± 1.1	8.7 ± 0.7	9.1 ± 1.3	H = 12.2, p = 0.002
	Feather THg (µg/g dry mass)	4.5 ± 1.2	4.6 ± 2.8	2.2 ± 0.6	H = 14.8, p = 0.001
	Subantarctic Zone
	Feathers (n)	25	34	20
	Feather δ 13 C (‰)	-19.7 ± 0.7	-19.0 ± 0.4	-19.6 ± 0.6	H = 21.6, p <0.0001
	Feather δ 15 N (‰)	11.6 ± 1.6	13.4 ± 0.7	13.3 ± 0.4	H = 21.6, p <0.0001
	Feather THg (µg/g dry mass)	11.9 ± 7.3	14.7 ± 3.9	14.4 ± 5.6	H = 7.0, p = 0.031
	Subtropical Zone
	Feathers (n)	11	9	7
	Feather δ 13 C (‰)	-17.2 ± 0.5	-17.8 ± 0.3	-17.5 ± 0.3	H = 7.4, p = 0.025
	Feather δ 15 N (‰)	12.3 ± 1.1	14.2 ± 0.5	14.4 ± 0.3	H = 12.9, p = 0.002
	Feather THg (µg/g dry mass)	8.8 ± 2.7	14.0 ± 4.2	19.4 ± 7.5	H = 11.7, p = 0.003
	Kruskal-Wallis H tests
	Feather δ 13 C (‰)	H =

34.2, p <0.0001 H = 36.2, p <0.0001 H = 32.6, p <0.0001

	Feather δ 15 N (‰)	H = 4.3, p = 0.116	H =

29.9, p <0.0001 H = 32.4, p <0.0001 Feather

  THg (µg/g dry mass) H

= 14.01, p = 0.001 H = 19.5, p <0.0001 H = 26.6, p <0.0001

	All three oceanic zones
	Feathers (n)	44	52	40
	Feather δ 13 C (‰)	-19.5 ± 1.8	-19.4 ± 1.6	-21.1 ± 3.2	H = 8.2, p = 0.016
	Feather δ 15 N (‰)	11.7 ± 1.4	12.7 ± 2.0	12.1 ± 2.3	H = 14.2, p = 0.001
	Feather THg (µg/g dry mass)	9.8 ± 6.3	12.8 ± 5.3	11.3 ± 8.3	H = 8.9, p = 0.012

Acknowledgements. The authors thank all those colleagues who helped with collecting albatross feathers at the different breeding sites, G. Guillou from the Plateforme Analyses Isotopiques of LIENSs for running stable isotope samples and C. Churlaud and M. Brault-Favrou from the Plateforme Analyses Elémentaires of LIENSs for their technical advices during Hg analyses. The present work was partly supported financially and logistically by the Institut Polaire Français Paul Emile Victor (IPEV, programme no. 109, H. Weimerskirch), the Terres Australes et Antarctiques Françaises (TAAF), and the British Antarctic Survey. This paper represents a contribution to the Ecosystems programme of the British Antarctic Survey. S. Jaquemet benefited from the National Research

Supplementary material

Table S1. Conservation status, breeding frequency, sampling sites, stable isotope values (δ 13 C and δ 15 N) and total mercury (THg) concentrations in body feathers of breeding albatrosses from the Southern Ocean and fringing subtropical waters. Values are means ± SD.


Sebastien Valeyre
email: sebastien.valeyre@jl-investments.com

Denis S Grebenkov
email: denis.grebenkov@polytechnique.edu

Sofiane Aboura
email: sofiane.aboura@univ-paris13.fr

Emergence of correlations between securities at short time scales

The correlation matrix is the key element in optimal portfolio allocation and risk management. In particular, the eigenvectors of the correlation matrix corresponding to large eigenvalues can be used to identify the market mode, sectors and style factors. We investigate how these eigenvalues depend on the time scale of securities returns in the U.S. market. For this purpose, oneminute returns of the largest 533 U.S. stocks are aggregated at different time scales and used to estimate the correlation matrix and its spectral properties. We propose a simple lead-lag factor model to capture and reproduce the observed time-scale dependence of eigenvalues. We reveal the emergence of several dominant eigenvalues as the time scale increases. This important finding evidences that the underlying economic and financial mechanisms determining the correlation structure of securities depend as well on time scales.

Introduction

How do the eigenvalues of securities correlation matrices emerge at different time scales? This fundamental question is important because crosscorrelations change over different investment horizons while a reliable empirical determination of the correlation matrix remains difficult due to its time and frequency dependence. This was first evidenced by Epps, who demonstrated the decay of correlations among U.S. stocks when shifting from daily to intra-daily time scales (or frequencies) [START_REF] Epps | Comovements in Stock Prices in the Very Short Run[END_REF]. In other words, the price correlation decreases with the duration of the time interval over which price changes are measured. The economic argument behind the Epps effect is that the information is not instantaneously transmitted at shorter time intervals, where the average adjustment lag in response of prices lies approximately between 10 and 60 minutes. This appears to reduce the scope of the Efficient Market Hypothesis [START_REF] Fama | Efficient capital markets: A review of theory and empirical work[END_REF] at short time scales given that tick data prices seem to adjust to new information only after a lag time, thus do not reflect all available information. Since its inception, the Epps effect has been confirmed by several studies, although its impact has been progressively declined in the NYSE, indicating that the market becomes increasingly more efficient [START_REF] Tóth | Increasing market efficiency: Evolution of crosscorrelations of stock returns[END_REF].

The dependence of securities cross-correlations on time scales can be captured via the eigenvalues of the correlation matrix. In particular, the largest eigenvalue reflects changes in the average correlation between stocks, whereas the corresponding eigenvector is associated to the "market mode". Kwapien et al. showed a significant elevation of the largest eigenvalue with increasing time scale using data from 1 minute to 2 days from NYSE, NASDAQ and Deutsche Börse (1997)(1998)(1999) [START_REF] Kwapien | Time scales involved in emergent market coherence[END_REF]. Using high-frequency stock returns from NYSE, AMEX and NASDAQ (1994)(1995)(1996)(1997), Plerou et al. supported the idea that the largest eigenvalue and its eigenvector reflect the collective response of the entire market to stimuli such as certain news breaks (e.g., central bank interest rates hikes) [START_REF] Plerou | Random matrix approach to cross correlations in financial data[END_REF]. This is particularly true during periods of high volatility when the collective behavior is enhanced. Coronnello et al. confirmed that the largest eigenvalue, computed from 5-minute data, describes the common behavior of the stocks composing the LSE stock index (2002) [START_REF] Coronnello | Sector Identification in a Set of Stock Return Time Series Traded at the London Stock Exchange[END_REF].

As firms having similar business activities are correlated, some other eigenvectors can economically be interpreted as business sectors [START_REF] Gopikrishnan | Quantifying and interpreting collective behavior in financial markets[END_REF]. So, Gopikrishnan et al. computed the eigenvectors of cross-correlation matrices of 1000 U.S. stocks at a 30-minute scale (1994)(1995) and a 1-day scale [START_REF] Gopikrishnan | Quantifying and interpreting collective behavior in financial markets[END_REF]. They found that the correlations in a business sector, captured via an eigenvector, were stable in time and could be used for the construction of optimal portfolios with a stable Sharpe ratio. In the same vein, as similar trading strategies induce cross-correlations in stocks, some eigenvectors can be financially interpreted as style factors. The corresponding eigenvalues are thus expected to exhibit non-trivial dependence on time scales. However, an accurate statistical analysis of multiple eigenvalues at different time scales is challenging due to measurement noises. In fact, as the correlation matrix is estimated from time series of stocks' returns, its elements are unavoidably random and thus prone to fluctuations. These fluctuations become larger as the length of time series is reduced, i.e., when the time scale is increased. While the largest eigenvalue typically exceeds the level of fluctuations by two orders of magnitude, the other eigenvalues rapidly reach this level and become non-informative. Several researchers employed the random matrix theory to distinguish economically significant eigenvalues from noise [START_REF] Laloux | Noise Dressing of Financial Correlation Matrices[END_REF][START_REF] Pafka | Exponential Weighting and Random-Matrix-Theory-Based Filtering of Financial Covariance Matrices for Portfolio Optimization[END_REF][START_REF] Burda | Signal and noise in correlation matrix[END_REF][START_REF] Potters | Financial applications of random matrix theory: Old laces and new pieces[END_REF][START_REF] Conlon | Random matrix theory and fund of funds portfolio optimisation[END_REF]. In particular, Laloux et al. showed that only 6% of the eigenvalues carried some information of the S&P 500 (1991)(1992)(1993)(1994)(1995)(1996), while the remaining 94% eigenvalues were hidden by noise [START_REF] Laloux | Noise Dressing of Financial Correlation Matrices[END_REF]. Guhr and Kalber proposed an alternative statistical approach to reduce noise that they called "power mapping" [START_REF] Guhr | A new method to estimate the noise in financial correlation matrices[END_REF]. Andersson et al. extended this work by comparing the power mapping approach to a standard filtering method discarding noisy eigenvalues for Markowitz portfolio optimization using daily Swedish stock market returns (1999)(2000)(2001)(2002)(2003) [START_REF] Andersson | Power Mapping and Noise Reduction for Financial Correlations[END_REF].

In this paper, we consider the correlation matrix of financial securities and investigate the emergence of its eigenvalues at small time scales. As the financial literature on this critical issue remains sparse, this research fills the gap by investigating the eigenvalues at intraday time scales using 1-min returns. We propose a simple model, coined the "lead-lag factor model", as an adaptation of the well-known "one-factor marker model" [START_REF] Ross | The Arbitrage Theory of Capital Asset Pricing[END_REF] to smaller time scales and to multiple sectors and style factors. In this model, stock returns are correlated to the returns of selected factors at earlier time steps. A detailed description of the eigenvalues as functions of the time scale is then derived. An empirical validation is performed on long time series of 1-min returns of a large universe of U.S. stocks. To get several significant eigenvalues at time scales from 1 minute to 2 hours, the correlation matrix was estimated over the whole available period (2013-2017) so that variations of cross-correlations over time were ignored (note that the dynamics of the eigenvalues and eigenvectors over time has been investigated elsewhere [START_REF] Conlon | Cross-correlation dynamics in financial time series[END_REF][START_REF] Allez | Eigenvector dynamics: General theory and some applications[END_REF][START_REF] Buccheri | Evolution of correlation structure of industrial indices of U.S. equity markets[END_REF]). In spite of its simple character, the lead-lag factor model is shown to be able to reproduce the dependence of large eigenvalues on the time scale.

The paper is organized as follows. In Sec. 2, we estimate the correlation matrix of U.S. stocks' returns at different time scales and present the empirical dependence of large eigenvalues on the time scale. To rationalize the observed behavior, we develop in Sec. 3 the lead-lag factor model and compare it to empirical results. Section 4 summarizes and concludes. Some derivations and more technical analysis of the lead-lag factor model are presented in Appendices.

Empirical results

Data description

We study the correlation structure of a universe that includes 533 U.S. stocks whose capitalization exceeded 1 billion dollars in 2013. For the considered period from 1st of January 2013 to 28th of June 2017, our database contains 338 176 1-min returns for each stock. We have also verified that the arithmetic aggregation of returns, r i (1) + . . . + r i (τ ), is almost identical to considering the product (1 + r i (1)) . . . (1 + r i (τ )) -1, given that the 1-min returns r i (t) are very small.

From the time series of 1-min returns, we estimate the correlation matrix over the whole available period, and then compute its eigenvalues. Then we aggregate the returns into 2-min, 4-min, ..., 128-min returns, producing time series with 169 088, 84 544, ..., 2 642 points, respectively. At each time scale τ , we repeat the computation to investigate the dependence of the eigenvalues on τ .

Empirical results

Figure 1a shows the four largest eigenvalues of the covariance matrix of 533 U.S. stocks' returns, computed by aggregating 1-min returns with the time scale τ , ranging from 1 minutes to 128 minutes (2 hours). The first two eigenvalues exhibit almost linear growth with τ , the others show minor deviations from linearity at small τ but scale linearly with τ at large τ . This behavior reflects the diffusion-like growth of the variance of aggregated returns; in particular, if the returns were independent, the eigenvalues of the corresponding covariance matrix, C ij = τ σ 2 i δ ij , would be just λ i = τ σ 2 i , and thus proportional to τ . Although correlations affect this linear growth, their effect is subdominant, at least for large eigenvalues, as witnessed by Fig. 1a. To highlight the effect of correlations, we focus on the eigenvalues of the correlation matrix. This choice is also justified from the financial point of view to level off the variability of stocks volatilities.

Figure 1b shows the four largest eigenvalues of the correlation matrix of the same 533 U.S. stocks' returns. If the returns were independent, the correlation matrix would be the identity, and thus all its eigenvalues would be equal to 1. The growth of these eigenvalues with the time scale τ indicates strong cross-correlations between stocks. The largest eigenvalue can be naturally attributed to the market mode, whereas the next eigenvalues correspond to different sectors and style factors.

After a sharp growth at short time scales (few minutes), the eigenvalues slowly approach to their long-time limits. The existence of these upper bounds is expected because the sum of eigenvalues of a correlation matrix is equal to its size (i.e., to the number of stocks, N). This saturation effect contrasts with the unlimited growth of eigenvalues of the covariance matrix (Fig. 1a). Finding the functional form of this approach and identifying its characteristic time scales present the main aim of our work. Recently, Benzaquen et al. proposed a multivariate linear propagator model for dissecting crossimpact on stock markets and revealing their dynamics [START_REF] Benzaquen | Dissecting cross-impact on stock markets: an empirical analysis[END_REF]. Due to its very general form accounting for both cross-correlations and auto-correlations of stocks, the proposed model contains too many parameters, while the resulting formulas are not explicit. Our ambition is rather the opposite and consists in suggesting an explicit model, as simple as possible, that would capture the empirical results shown in Fig. 1b and thus provide a minimalistic framework for their financial interpretation.

The lead-lag factor model

Basic lead-lag one-factor model

We consider a trading universe with N assets. In a conventional onefactor model, the return of the i-th asset at time t, r i (t), is modeled as a combination of a specific, asset-dependent random fluctuation, ε i (t), and an overall market contribution, R(t),

r i (t) = ε i (t) + βR(t), (1)
with a market sensitivity β (that we generalize below to other factors). The asset-specific random fluctuations ε i (t) are typically modeled as independent centered Gaussian variables with volatilities σ i .

We propose a modification of this conventional model by incorporating the lead-lag effect, in which the i-th asset return at time t is influenced by a common factor R(t -k) at earlier times t -k, with progressively decaying weights:

r i (t) = ε i (t) + β ∞ k=0 α k R(t -k), (2)
where 0 ≤ α < 1 characterizes the relaxation time of the memory decay. Note that the upper limit of the sum in Eq. ( 2) is formally extended to infinity, bearing in mind that contributions for very large k are exponentially small. We will analyze the model in the stationary regime as t → ∞ in order to eliminate transient effects.

The common term R(t) can be interpreted as an idealized factor without auto-correlations in an efficient market that most stocks follow with a lead lag delay. We model therefore R(t) by independent centered Gaussian variables with volatility Σ. The term R(t) can represent the market mode but also sectors or style factors, or any popular trading portfolio. Moreover, R(t) can also be interpreted as being linked to the market order transactions for a particular strategy (market, sector or styles). In this light, our model can be seen as an extension of the Kyle model [START_REF] Kyle | Continuous Auctions and Insider Trading[END_REF] that explains the impact of transactions on price for a single stock and without delay. Here, we consider multiple stocks and include an exponential decay of the impact. While more sophisticated models with a power law decay of the impact were proposed [START_REF] Benzaquen | Dissecting cross-impact on stock markets: an empirical analysis[END_REF][START_REF] Bouchaud | How markets slowly digest changes in supply and demand[END_REF], we will show that our minimalistic model is enough to reproduce a slow growth of the eigenvalues of the correlation matrix. For the sake of clarify, we first analyze this basic lead-lag one-factor model and then discuss its several straightforward extensions.

The one-factor relation ( 2) is the basic model for returns at the smallest time scale. We then consider the returns aggregated on the time scale τ :

r τ i (t) = τ -1 ℓ=0 r i (t -ℓ), (3)
with t being a multiple of τ . Under the former Gaussian assumptions, the covariance function of the aggregated returns reads (see Appendix A):

C τ ij = r τ i (t)r τ j (t) = τ σ 2 δ ij + β 2 Σ 2 τ (1 -α 2 ) -2α(1 -α τ ) (1 -α 2 )(1 -α) 2 , ( 4
)
where • • • denotes the expectation, and δ ij = 1 for i = j, and 0 otherwise. Note that we set here σ i = σ for all assets for simplicity (this simplification will be relaxed below). As we consider the stationary regime, the covariance function does not depend on time t. Denoting

κ α (τ ) = τ (1 -α 2 ) -2α(1 -α τ ) 1 -α 2 , (5)
one gets the correlation matrix

C τ ij = C τ ij C τ ii C τ jj = 1 i = j, ρ 2 (τ ) i = j, (6)
with

ρ(τ ) = 1 + η(τ )/γ -1/2 , (7)
where

γ = Σ 2 β 2 σ 2 (8)
and

η(τ ) = τ κ α (τ )/(1 -α) 2 = (1 -α) 2 1 -2α 1-α 2 (1 -α τ )/τ . ( 9
)
The function η(τ ), that will play the central role in our analysis, monotonously decreases from

η(1) = 1 -α 2 to η(∞) = (1 -α) 2 .
Since the matrix C τ -(1 -ρ 2 (τ ))I has rank 1 (I being the identity N × N matrix), there are N -1 eigenvalues λ i = 1 -ρ 2 (τ ). In turn, the single largest eigenvalue of the correlation matrix C τ can be obtained as follows: N = Tr(C τ ) = λ 1 + (N -1)λ i , from which λ 1 = 1 + (N -1)ρ 2 (τ ). We get thus the complete description of the eigenvalues as functions of the time scale τ :

λ 1 = 1 + (N -1)ρ 2 (τ ), ( 10
) λ i = 1 -ρ 2 (τ ) (i = 2, 3, . . . , N). (11)
In the limit of very large τ , one finds

ρ 2 (∞) = 1 + (1 -α) 2 /γ -1 . ( 12
)
This simplest lead-lag one-factor model predicts a monotonous growth of the largest eigenvalue (corresponding to the market mode) with the time scale τ , up to a saturation plateau. In turn, the other eigenvalues exhibit a monotonous decrease to a plateau. In spite of the exponential decay of the lead-lag memory effect in Eq. ( 2), the approach to the plateau is governed by a slow, 1/τ power law, in a qualitative agreement with the empirical observation (see Sec. 3.5 for quantitative comparison). In particular, this approach has no well-defined time scale. While the basic model can potentially capture the behavior of the largest eigenvalue, it clearly fails to distinguish other eigenvalues. One needs therefore to relax some simplifying assumptions to render the model more realistic.

General lead-lag one-factor model

We start by introducing arbitrary volatilities σ i and sensitivities β i of the i-th asset to the common factor R(t):

r i (t) = ε i (t) + β i ∞ k=0 α k R(t -k). (13)
In this case, the computation is precisely the same, the only difference is that

C τ ij = τ σ 2 i δ ij + Σ 2 β i β j κ α (τ ). ( 14
)
As a consequence, the structure of the correlation matrix is fully determined by β i , whereas the dependence on the time scale τ is still represented by κ α (τ ). The correlation matrix reads

C τ ij = 1 (i = j), ρ i (τ )ρ j (τ ) (i = j), (15)
with

ρ i (τ ) = 1 + η(τ )/γ i -1/2 , γ i = Σ 2 β 2 i σ 2 i . ( 16
)
The eigenvalues of this correlation matrix can be computed as follows.

If all γ i are distinct 1 , the components of an eigenvector are

v i = ρ i Q λ -1 + ρ 2 i (i = 1, . . . , N), with Q = N i=1 ρ i v i , (17)
from which one gets the equation on the eigenvalues λ

N i=1 ρ 2 i λ -1 + ρ 2 i = 1. (18)
This equation has N distinct solutions that can be characterized in terms of ρ 2 i (see Appendix B). When N is large, the largest eigenvalue is expected to be large, and the asymptotic expansion of Eq. ( 18) yields

λ 1 ≈ N i=1 ρ 2 i = N i=1 1 + η(τ )/γ i -1 . (19)
In turn, the other eigenvalues are below 1 (see Appendix B). As a consequence, such a lead-lag one-factor model cannot reproduce several eigenvalues larger than 1. For this purpose, one needs to consider multiple factors.

General lead-lag multi-factor model

Now we consider a general lead-lag multi-factor model

r i (t) = ε i (t) + ∞ k=0 α k F f =1 β i,f R f (t -k), (20)
1 When some γ i are identical, the analysis of eigenvalues becomes more involved (see Appendix B), but the largest eigenvalue still satisfies Eq. ( 18) and can thus be approximated by Eq. [START_REF] Benzaquen | Dissecting cross-impact on stock markets: an empirical analysis[END_REF]. In particular, if all γ i = γ, one gets λ 1 ≈ N ρ 2 (τ ), which is close to the exact solution [START_REF] Burda | Signal and noise in correlation matrix[END_REF].

where ε i (t) are independent centered Gaussian variables (representing random fluctuations specific to the stock i) with variance σ 2 i , F is the number of factors, R f (t) are independent centered Gaussian returns of the factor f with variance Σ 2 f , β i,f is the sensitivity of the stock i to the factor f , and α sets the relaxation time. Repeating the computation from Appendix A, one gets

C τ ij = δ ij + (1 -δ ij ) F f =1 ρ i,f ρ j,f , (21)
where

ρ i,f (τ ) = Σ f β i,f β i ρ i (τ ) (22)
and

ρ i (τ ) = 1 + η(τ )/γ i -1/2 , γ i = β 2 i σ 2 i , β 2 i = F f =1 Σ 2 f β 2 i,f . (23)
Considering ρ i,f as the elements of an N × F matrix ρ, one can rewrite Eq. ( 21) in a matrix form

C τ = (I -P ) + ρρ † , ( 24
)
where P is the diagonal matrix formed by ρ 2 i , and † denotes the matrix transpose.

The matrix ρ of size N ×F plays the central role in the following analysis. As the elements of the matrix ρ are real, ρρ † , as well as ρ † ρ, are positive semidefinite matrices which have nonnegative eigenvalues. The rank of the matrix ρ is equal to that of matrices ρρ † and ρ † ρ and thus cannot exceed min{F, N}. Given that F ≪ N, the correlation matrix C τ appears as the perturbation of a diagonal matrix by a low-rank matrix.

The eigenvalues of the correlation matrix are the zeros of the determinant

0 = det λI -C τ = det λI -I + P -ρρ † . ( 25
)
Since ρρ † is a low-rank perturbation, one can expect, as in the one-factor case of Appendix B, that most eigenvalues coincide with that of the unperturbed diagonal matrix I -P , i.e., they are given by 1 -ρ 2 i for some indices i. These eigenvalues are essentially hidden by noise and non-exploitable in practice. We are interested in large eigenvalues that (significantly) exceed 1.

Practical approximation

As we will discuss in detail in Sec. 3.5, empirical data exhibit the shortrange memory effect (α is small) and the relatively small impact of the factors onto the variance of individual stocks as compared to the stock-specific fluctuations (γ i are small). In this situation, which is particular to the time series of securities returns at the considered time scales, one has η(τ )/γ i ≫ 1 so that ρ i (τ ) in Eq. ( 23) can be approximated as

ρ 2 i (τ ) ≃ γ i η(τ ) . ( 31
)
This approximation greatly simplifies the elements of the matrix ρ † ρ:

(ρ † ρ) f,g = N i=1 Σ f β i,f ρ i (τ ) β i =ρ i,f Σ g β i,g ρ i (τ ) β i =ρ i,g ≈ N η(τ ) Γ f,g , (32)
where the matrix elements Γ f,g do not depend on the time scale:

Γ f,g = Σ f Σ g N N i=1 β i,f β i,g σ 2 i . ( 33
)
As a consequence, all the elements of the matrix ρ † ρ and thus its eigenvalues exhibit the same dependence on the time scale τ , expressed via the explicit function η(τ ) given by Eq. ( 9). Denoting the eigenvalues of the matrix Γ as γ f (f = 1, . . . , F ), one gets the following approximation for large eigenvalues of the correlation matrix:

λ f ≈ Nγ f η(τ ) (f = 1, . . . , F ) . ( 34
)
From the explicit form ( 16) of η(τ ), one deduces a slow, 1/τ , power law approach of the eigenvalue to the saturation level as the time scale τ increases. Within this approximate computation, all large eigenvalues exhibit the same dependence on the time scale. In practice, one aims at constructing the factors R f to capture independent features of cross-correlations in the market. The sensitivies β i,f and β i,g of the stock i to factors R f and R g are thus expected to be "orthogonal", and this property can be formally expressed by requiring that the nondiagonal elements of the matrix Γ are negligible. In this case, the eigenvalues γ f are given by the diagonal elements

γ f = Γ f,f = 1 N N i=1 Σ 2 f β 2 i,f σ 2 i . ( 35
)
This is a kind of empirical mean of the squared sensitivities β 2 i,f , normalized by the squared volatilities σ 2 i .

Application to empirical data

We aim at applying the lead-lag factor model to fit the eigenvalues of the empirical correlation matrix of U.S. stocks' returns. The fitting formula (34) has two adjustable parameters: the relaxation time α in the function η(τ ) and the amplitude Nγ f . Using the least square fitting algorithm implemented as the routine lsqcurvefit in Matlab, we apply the formula (34) separately to each empirical eigenvalue.

Figure 2 shows the fitting of the four largest eigenvalues. The good quality of the fit by the lead-lag factor model indicates that, in spite of numerous simplifying assumptions on which the model was built, it captures the overall behavior qualitatively well. In particular, the eigenvalues converge to limiting values, at least for the considered short-time scales (up to 2 hours). Moreover, this saturation level is approached slowly, with the characteristic 1/τ power law dependence. The adjustable parameters are summarized in Table 1. Rewriting the attenuation factor α k in the lead-lag factor model (2) as exp(-t/t α ) with t = kτ 0 and t α = τ 0 / ln(1/α), where τ 0 = 1 min is the finest time scale of the time series used, one gets the relaxation time t α in minutes. One can see that the relaxation times α (or t α ) for four eigenvalues are close to each other. In other words, all the dominant eigenmodes evolve 34) of the four largest eigenvalues of the correlation matrix of N = 533 U.S. stocks' returns, computed by aggregating 1-min returns with the time scale τ . The adjustable parameters α and γ are summarized in Table 1.

at comparable time scales. This is an important conclusion which refutes a common belief that the market mode (corresponding to the largest eigenvalue) evolves at a time scale that is significantly different from other modes (sectors and style factors). The values of t α are of the order of one minute, in agreement with predictions by Benzaquen et al. [START_REF] Benzaquen | Dissecting cross-impact on stock markets: an empirical analysis[END_REF]. Remarkably, while the lead-lag memory effects vanish so rapidly, they impact the behavior of the eigenvalues at much longer time scales. In particular, if the lead-lag was ignored (by setting α = 0), the largest eigenvalue would be ≃ Nγ 1 and independent of the time scale τ . For instance, using the estimated value γ 1 = 0.17 and setting α = 0, one would get the largest eigenvalue to be 90, which is significantly smaller than the expected limit 128 for α = 0.16 or the observed value 130 at τ = 128 min.

Conclusion

We investigated the dependence of the eigenvalues of the correlation matrix on the time scale τ . Aggregating 1-min returns of the largest 533 U.S. stocks (2013-2017) to estimate the correlation matrix at different time scales, we showed that its large eigenvalues grow with τ and apparently saturate to limiting values. This growth reflects the important phenomenon that interstock correlations accumulate over time scales.

To rationalize this phenomenon and to interpret empirical observations, we developed the lead-lag factor model. In the one-factor case, each stock is considered to be partly correlated to a given lead-lag factor. Under several simplifying assumptions, we derived a simple formula for large relevant eigenvalues. This formula containing just two easily interpretable adjustable parameters, was then validated on empirical data.

The relaxation time of the stock market was estimated to be around 1 minute. A possible interpretion of this observation can be that a transaction can generate a cascade of transactions that decays in 1 minute so that the impact of transaction on price decays in 1 minute. As correlations emerge from the cross-impact of transactions on prices, we model this effect by extending the Kyle model to the impact of transaction on preferential portfolios with a lead lag effect.

The small value of the observed relaxation time suggests that correlation measurements based on 5 minutes returns should provide a good proxy of correlation of daily returns for risk management, in line with the conclusion by Liu et al. on volatility estimation [START_REF] Liu | Does anything beat 5-minute RV? A comparison of realized measures across multiple asset classes[END_REF]. However, other phenomena are likely to occur at much larger time scales (from day to month), e.g., autocorrelations of returns of financial factors (book, size, momentum) due to herding effect, or lack of liquidity. An accurate estimation of correlations at larger time scales remains a challenging problem because of a limited number of the available returns and thus higher impact of noise in the estimated correlation matrix. To overcome this limitation, one can either consider time horizons over several decades (in which case neglecting variations of corrections over time becomes debatable), or reduce the number of considered securities and thus the dimension of the correlation matrix (in which case financial meaning of estimated correlations may be debatable). A possible solution consists in constructing relevant financial factors and investigating how their correlations change with the time scale, as suggested by our factor-based model.

Appendix B. Analysis of the lead-lag one-factor model

We study in more detail the model [START_REF] Ross | The Arbitrage Theory of Capital Asset Pricing[END_REF] of the correlation matrix C, with ρ i (τ ) given by Eq. ( 16). This matrix is a perturbation of the identity matrix by a rank one matrix, for which many spectral properties are known (see, e.g., [START_REF] Golub | Some modified matrix eigenvalue problems[END_REF]). This matrix combines both effects: the correlation coefficient ρ and the impact of the exponential moving average (with the coefficient α). We search for an eigenvector of this matrix as

v = (v 1 , v 2 , . . . , v n ) † . Writing explicitly Cv = λv, we get v i (1 -ρ 2 i ) + ρ i Q = λv i (i = 1, . . . , N), (B.1)
where

Q = N i=1 v i ρ i . (B.2)
First, we note that if ρ i = 0 for some i, then the above equation is reduced to v i = λv i that has two solutions: either λ = 1 and v i can be arbitrary; or v i = 0. One can check that if ρ i 1 = . . . = ρ i k = 0 for k stocks, then the correlation matrix has the eigenvalue λ = 1 with the multiplicity k. The corresponding eigenvectors can be chosen as an orthogonal basis in the subspace R k . In turn, the remaining n -k eigenvalues are nontrivial, and can be determined as discussed below. In what follows, we focus on these nontrivial eigenvalues, i.e., we assume that all ρ i = 0.

The equation (B.1) has two solutions:

(i) either λ = 1 -ρ 2 i and Q = 0; or (ii) λ = 1 -ρ 2 i and v i = ρ i Q λ -1 + ρ 2 i . (B.3)
In the latter case, one can substitute this expression into Eq. (B.2) to get an equation on the eigenvalue λ:

N i=1 ρ 2 i λ -1 + ρ 2 i = 1. (B.4)
This equation can be seen as a polynomial of degree N which has N (a priori complex-valued) zeros. Finally, Q can be fixed by setting the normalization condition on v:

1 = N i=1 v 2 i = Q 2 N i=1 ρ 2 i (λ -1 + ρ 2 i ) 2 . (B.5)
This is a generic situation.

Let us return to the first option, namely, we suppose that λ = 1 -ρ 2 k for some index k that implies that Q = 0. If all ρ i are distinct, i.e., ρ 1 = ρ 2 = . . . = ρ N , so that v i = 0 for all i = k, but, due to Q = 0, it would also imply that v k = 0. As a consequence, v = 0 but this is not an eigenvector. We conclude that, if all ρ i are distinct, then λ cannot be given by 1 -ρ 2 i , and this option is excluded. Now, we consider the case when two or more values ρ i are identical. For instance, let us assume that ρ 1 = ρ 2 = ρ 3 = . . . = ρ N . In this case, λ = 1 -ρ 2 1 is indeed an eigenvalue. In fact, one gets Q = 0 and thus v i = 0 for i > 2. However, one has

Q = ρ 1 v 1 + ρ 2 v 2 = ρ 1 (v 1 + v 2 ) = 0, implying that v 1 = -v 2 .
The normalization condition implies thus v 1 = -v 2 = 1/ √ 2. We conclude that λ = 1 -ρ 2 1 is then a single eigenvalue. More generally, if ρ 1 = ρ 2 = . . . = ρ k = ρ k+1 = . . . = ρ N , then the eigenvalue λ = 1 -ρ 2 1 has the multiplicity k -1.

In general, it is convenient to denote z i = 1 -ρ 2 i and to order them in an increasing order:

z 1 ≤ z 2 ≤ z 3 ≤ . . . ≤ z N (B.6)
or, equivalently, by grouping the eventual identical values: In other words, there are i 1 identical values z 1 = . . . = z i 1 ; i 2 identical values z i 1 +1 = . . . = z i 1 +i 2 , etc. (note that when all z i are distinct, one has i 1 = i 2 = . . . = 1). In this configuration, the correlation matrix has: the eigenvalue z 1 with the multiplicity i 1 -1 (if i 1 > 1); the eigenvalue z 2 with the multiplicity i 2 -1 (if i 2 > 1); etc. If for some k, z i k = 1, then this eigenvalues has the multiplicity i k . Finally, the remaining eigenvalues are determined as solutions of Eq. (B.4) that can be written as f (z) = 1, with

f (z) = N i=1 ρ 2 i z -1 + ρ 2 i = N i=1 1 -z i z -z i . (B.8)
The terms with z i = 1 (resulting in the eigenvalue λ = 1) are excluded from this sum. Moreover, if some z i are identical, the corresponding terms are just grouped together. As a consequence, the equation f (z) = 1 is reduced to a polynomial of degree at most N (the degree N corresponding to the case when all z i are distinct).

It is worth noting that the function f (z) is decreasing everywhere:

f ′ (z) = - N i=1 ρ 2 i (z -z i ) 2 < 0. (B.9)
As a consequence, one gets immediately that each interval (z i , z i+1 ) (with z i < z i+1 and z N +1 = ∞) has exactly one solution of the equation f (z) = 0, i.e., one eigenvalue. In particular, one gets the following bounds for the smallest eigenvalue

z 1 ≤ min 1≤i≤N {λ i } ≤ z 2 . (B.10)
We conclude that all eigenvalues are positive if and only if z 1 ≥ 0, i.e., ρ 2 i ≤ 1 for all i. In other words, the inequalities ρ 2 i ≤ 1 for all i present the necessary and sufficient condition for the positive definiteness of the matrix. These conditions are evidently satisfied in our setting.

Since f (1) ≥ 1, one gets the following bound for the largest eigenvalue (note that the eigenvalues are ordered in descending order, λ 1 ≥ λ 2 ≥ . . ., in contrast to z k ). However, this bound is rather weak. In turn, since λ 2 ≤ z N = 1 -ρ 2 N < 1, all other eigenvalues are below 1: λ i < 1 (i = 2, 3, . . . , N). (B.12)

Figure 1 :

 1 Figure 1: Four largest eigenvalues of the covariance matrix (a) and of the correlation matrix (b) for returns of 533 U.S. stocks, computed by aggregating 1-min returns with the time scale τ , varying from 1 minute to 128 minutes (2 hours).

Figure 2 :

 2 Figure 2: Fitting by Eq. (34) of the four largest eigenvalues of the correlation matrix of N = 533 U.S. stocks' returns, computed by aggregating 1-min returns with the time scale τ . The adjustable parameters α and γ are summarized in Table1.

z 1 =

 1 z 2 = . . . = z i 1 < z i 1 +1 = z i 1 +2 = . . . = z i 1 +i 2 < . . . < z i 1 +...+im = z i 1 +...+im+1 = . . . = z N .(B.7)

Table 1 :

 1 Two adjustable parameters of the fitting formula (34) applied to four largest eigenvalues of the correlation matrix of N = 533 U.S. stocks' returns. The corresponding relaxation time t α in minutes is obtained as 1 min/ ln(1/α).

	f	1	2	3	4
	γ	0.17 0.03 0.02 0.01
	α	0.16 0.25 0.18 0.26
	t α (min) 0.55 0.72 0.58 0.74

If λ exceeds 1, it cannot be equal to 1 -ρ 2 i for all i, the matrix λI -I + P is nonsingular, its inverse exists, so that one can rewrite Eq. (25) as 0 = det λI -I + P det I -ρ † (λI -I + P ) -1 ρ ,

from which one gets a new equation on eigenvalues:

(here we used a general property: if A ∈ C m×m is nonsingular matrix and U, V ∈ C m×r , then det(A + UV * ) = det(A)det(I + V * A -1 U), see [START_REF] Xiong | On the eigenvalues of a specially rank-r updated complex matrix[END_REF]). Denoting the F × F matrix in the determinant as φ(λ), one can write explicitly its elements as

The solutions of Eq. (27) determine some eigenvalues λ of the correlation matrix in Eq. [START_REF] Bouchaud | How markets slowly digest changes in supply and demand[END_REF]. As one typically deals with the situation N ≫ F , the reduction of the original determinant equation ( 25) for a matrix of size N ×N to Eq. ( 27) for a matrix of size F × F is a significant numerical simplification of the problem. Most importantly, this formal solution allows one to get analytical insights onto the eigenvalues, as we did in the one-factor case in Appendix B. Note that in the one-factor case (F = 1), the determinant equation ( 27) is simply reduced to

i.e., we retrieve Eq. ( 18). If one searches for large eigenvalues, λ ≫ 1, one can neglect the matrix P -I in comparison to λI in Eq. ( 27), that yields

In other words, the large eigenvalues of the correlation matrix can be approximated by the eigenvalues of the matrix ρ † ρ of size F × F . This symmetric positive semi-definite matrix has F nonnegative eigenvalues that correspond to F factors.

Appendix A. Computation of the covariance matrix

The covariance matrix of aggregated centered Gaussian returns r τ i (t) defined by Eq. ( 3) is

The first term in this expression comes from the uncorrelated stock-dependent fluctuations. The independence of returns R(k) implies

To calculate these four sums, it is convenient to consider separately various terms depending on ℓ 1 and ℓ 2 :

• there are τ terms with

(A.3)

• there are τ -1 terms with

Moreover, the same contribution comes from ℓ 1 = ℓ 2 -1 and k 1 = k 2 + 1.

• similarly, there are τ -j terms with

and this contribution is doubled by the symmetry argument.

• finally, there is one term with ℓ 1 = ℓ 2 +(τ -1) and thus k 1 = k 2 -(τ -1) whose contribution is α τ -1 /(1 -α 2 ).

Combining all these terms, one gets after simplifications Eq. ( 4).