Géraldine Patey
email: geraldine.patey@gmail.com

Catherine M Couillard
email: catherine.couillard@dfo-mpo.gc.ca

Hilaire Drouineau

Guy Verreault
email: guy.verreault@mffp.gouv.qc.ca

Fabien Pierron
email: f.pierron@epoc.u-bordeaux1.fr

Patrick Lambert
email: patrick.lambert@irstea.fr

Magalie Baudrimont
email: magalie.baudrimont@u-bordeaux.fr

Patrice Couture
email: patrice.couture@ete.inrs.ca

Saint François Lake

Early back-calculated size-at-age of Atlantic yellow eels sampled

Keywords: Anguilla rostrata, Anguilla anguilla, back-calculated length, growth, habitat, St. Lawrence River system, Gironde River system

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

A decline of North European and American eel populations (here called "Atlantic" eels) has been observed over the last 30 years [START_REF] Castonguay | Why is recruitment of the American eel, Anguilla rostrata, declining in the St. Lawrence River and Gulf?[END_REF][START_REF] Cosewic | Assessment and status report on the American eel Anguilla rostrata in Canada[END_REF][START_REF] Dekker | Status of the European eel stock and fisheries[END_REF]FAO/ICES 2009). Multiple factors have been evoked, including barriers to migration, overfishing, parasites, pollutants and oceanographic and climate changes [START_REF] Castonguay | Why is recruitment of the American eel, Anguilla rostrata, declining in the St. Lawrence River and Gulf?[END_REF][START_REF] Couillard | Assessment of Fat Reserves Adequacy in the First Migrant Silver American Eels of a Large-Scale Stocking Experiment[END_REF][START_REF] Geeraerts | The effects of contaminants in European eel: a review[END_REF], but it remains difficult to disentangle their relative contributions and they have probably acted in synergy [START_REF] Jacoby | Synergistic patterns of threat and the challenges facing global anguillid eel conservation[END_REF][START_REF] Miller | Did a "perfect storm" of oceanic changes and continental anthropogenic impacts cause northern hemisphere anguillid recruitment reductions?[END_REF].

Atlantic eels are two panmictic and catadromous species with very similar life cycles [START_REF] Avise | Catadromous eels of the North Atlantic: a review of molecular genetic findings relevant to natural history, population structure, speciation, and phylogeny[END_REF]. After a long trans-oceanic larval drift from their spawning grounds in the Sargasso Sea to the continental shelf, glass eels penetrate into coastal areas, then they grow into elvers and later into yellow eels [START_REF] Tesch | Migration from continental waters to the spawning grounds[END_REF]. Yellow eels complete their somatic growth phase in a large range of habitats (estuaries, rivers and lakes) [START_REF] Daverat | Tactics and demographic attributes in the European eel Anguilla anguilla in the Gironde watershed, SW France[END_REF][START_REF] Jessop | Geographic effects on American eel (Anguilla rostrata) life history characteristics and strategies[END_REF].

At the end of the yellow phase, they undergo a last metamorphosis, turn into silver eels and begin their long migration to the Sargasso Sea where they reproduce and die [START_REF] Tesch | Migration from continental waters to the spawning grounds[END_REF].

The growth phase associated with lipid accumulation is a key determinant of trans-oceanic migration, fecundity, and consequently of reproductive success [START_REF] Belpaire | Decreasing eel stocks: survival of the fattest?[END_REF][START_REF] Svedäng | Low fat contents in female silver eels: indications of insufficient energetic stores for migration and gonadal development[END_REF]. Since the metamorphosis into silver eels and the success of their reproductive migration are heavily dependent on body length and fat content at migration, growth rate of the yellow phase directly affects reproductive success [START_REF] Belpaire | Decreasing eel stocks: survival of the fattest?[END_REF][START_REF] Couillard | Assessment of Fat Reserves Adequacy in the First Migrant Silver American Eels of a Large-Scale Stocking Experiment[END_REF]Van Ginneken and Van den Thillart 2000).

Female American and European eels exhibit a latitudinal cline in growth rate, partly explained by temperature, with lower growth rate (at the silver stage) and higher size-atmetamorphosis at higher latitudes and/or longer distances from the spawning site [START_REF] Jessop | Geographic effects on American eel (Anguilla rostrata) life history characteristics and strategies[END_REF]Vélez-Espino and Koops 2009). Our current understanding is that within the St. Lawrence system (Canada), the growth rate of female eels decreases and size-at-metamorphosis increases from the Gulf of St. Lawrence (downstream section) to Lake Ontario (upstream section) [START_REF] Jessop | Geographic effects on American eel (Anguilla rostrata) life history characteristics and strategies[END_REF].

While eel's growth rate is known to be largely plastic [START_REF] Côté | In absence of local adaptation, plasticity and spatially varying selection rule: a view from genomic reaction norms in a panmictic species (Anguilla rostrata)[END_REF][START_REF] Drouineau | Assessing the impact of discretisation assumptions in a length-structured population growth model[END_REF], genetic/epigenetic differences linked to differential migratory patterns or to selection are increasingly supported as contributing factors for the early differential growth rates observed in elvers at different latitudes [START_REF] Côté | Differential effects of origin and salinity rearing conditions on growth of glass eels of the American eel Anguilla rostrata: implications for stocking programmes[END_REF][START_REF] Côté | In absence of local adaptation, plasticity and spatially varying selection rule: a view from genomic reaction norms in a panmictic species (Anguilla rostrata)[END_REF].

At the river basin scale, eels can use a wide range of habitats as demonstrated by studies using otolith microchemistry [START_REF] Daverat | Tactics and demographic attributes in the European eel Anguilla anguilla in the Gironde watershed, SW France[END_REF][START_REF] Thibault | Facultative catadromy in American eels: testing the conditional strategy hypothesis[END_REF]). This also results in variability in growth rates at this basin scale. Yellow eels with different migratory tactics have different growth rates, as shown in the Gironde System (France) for A. anguilla [START_REF] Daverat | One century of eel growth: changes and implications[END_REF][START_REF] Daverat | Tactics and demographic attributes in the European eel Anguilla anguilla in the Gironde watershed, SW France[END_REF][START_REF] Edeline | Dispersal in metamorphosing juvenile eel Anguilla anguilla[END_REF]Edeline et al. 2005a) and in various sites in the St. Lawrence system for A. rostrata [START_REF] Boivin | How does salinity influence habitat selection and growth in juvenile American eels Anguilla rostrata?[END_REF][START_REF] Cairns | Salinity-linked growth in anguillid eels and the paradox of temperate-zone catadromy[END_REF][START_REF] Jessop | Effects of inter-habitat migration on the evaluation of growth rate and habitat residence of American eels Anguilla rostrata[END_REF][START_REF] Thibault | Facultative catadromy in American eels: testing the conditional strategy hypothesis[END_REF]. Resident freshwater eels display lower growth rates than resident brackish Page 4 of 37 Can. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18

For personal use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version of record.

Author-produced version of the article published in Can. J. Fish. Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025

water eels, and the growth rates of nomadic migrants vary according to habitat salinity [START_REF] Daverat | Tactics and demographic attributes in the European eel Anguilla anguilla in the Gironde watershed, SW France[END_REF][START_REF] Thibault | Facultative catadromy in American eels: testing the conditional strategy hypothesis[END_REF]. As observed for the latitudinal variations of growth rate, differential growth patterns between salinity ecotypes are detected at an early age (glass eels and elvers) and appear to be related to genetic/epigenetic differences [START_REF] Côté | Differential effects of origin and salinity rearing conditions on growth of glass eels of the American eel Anguilla rostrata: implications for stocking programmes[END_REF][START_REF] Edeline | Dispersal in metamorphosing juvenile eel Anguilla anguilla[END_REF][START_REF] Pavey | RAD Sequencing Highlights Polygenic Discrimination of Habitat Ecotypes in the Panmictic American Eel[END_REF].

The present study is part of an international project examining the influence of anthropogenic and environmental factors on the decline of Atlantic eels. This large-scale sampling program provided a unique opportunity to investigate early size-at-age of yellow

Atlantic eels captured at the same time with similar protocols in different habitats of the Gironde and St. Lawrence systems. Previous studies demonstrated that both species exhibited, at an early age, extreme inter-individual phenotypic variance in terms of size-at-age related to genetic/epigenetic differences and to the use of highly heterogeneous environments [START_REF] Côté | Differential effects of origin and salinity rearing conditions on growth of glass eels of the American eel Anguilla rostrata: implications for stocking programmes[END_REF][START_REF] Côté | In absence of local adaptation, plasticity and spatially varying selection rule: a view from genomic reaction norms in a panmictic species (Anguilla rostrata)[END_REF][START_REF] Daverat | Tactics and demographic attributes in the European eel Anguilla anguilla in the Gironde watershed, SW France[END_REF][START_REF] Gagnaire | The genetic consequences of spatially varying selection in the panmictic American eel (Anguilla rostrata)[END_REF]. The objective of our study was to test if the current spatial patterns of early size-at-age of yellow Anguilla sp. in these two systems are consistent with the current paradigm of higher early growth rates of eels growing in habitats more saline and/or less distant from the Sargasso Sea spawning area. Early growth rates of yellow eels were compared within species among habitats more or less saline and/or distant from the Sargasso Sea, using back-calculated length from 1 to 5 years old. Then, we carried out linear mixed effect models for each site including random effects to account for possible inter-cohort's variation in lengths and individual random effects. Finally, the early sizeat-age of American and European eels were compared with historical data to examine the temporal evolution of eel early growth in recent years.

Materials and methods

The St. Lawrence Estuary is one of the largest and most productive ecosystems in Canada [START_REF] Dufour | Estuary and Gulf of St. Lawrence marine ecosystem overview and assessment report[END_REF]. It is 5 times longer and 9.5 times wider than the Gironde Estuary (GE) (Fig. 1 and Table S1). The GE is one of the largest estuaries in Western Europe and it originates from the confluence of two freshwater rivers, the Garonne and the Dordogne rivers, at the Bec d'Ambès. Both systems are composed of marine, brackish and freshwater sections (Fig. 1).

For each hydrographical system during early summer in 2011 and 2012, a total of 120 yellow eels were collected at 4 sites along brackish and freshwater sections (Table 1). Our scientific collector permits restricted the number of fish to 30 eels per site because of their threatened status, a major constraint of this study. As further constraint, fish were captured in two different years (15 eels per site per year) due to the low abundance of eels in several sites and the requirement of our sampling design for selecting specific size ranges (see below). In Quebec, eels were captured in two fluvial lakes, St. François Lake (SFL) and St. Pierre Lake (SPL) located in the freshwater section and in the mouth of two freshwater tributaries of the Estuary and the Gulf of St. Lawrence respectively, the Sud-Ouest River (SOR) and St. Jean River (SJR), (Fig. 1a and Table 1). SFL and SPL are shallow freshwater fluvial lakes near Montréal in an urbanized and industrialized region [START_REF] Béguer-Pon | Large-scale, seasonal habitat use and movements of yellow American eels in the St. Lawrence River revealed by acoustic telemetry[END_REF][START_REF] Hodson | Mercury concentrations in amphipods and fish of the Saint Lawrence River (Canada) are unrelated to concentrations of legacy mercury in sediments[END_REF]. Previous studies have shown that the estimated average transit times of the upstream migration of early life stages of eels between the Gulf and these lakes are 3 years and 5-6 years for SPL and SFL respectively (G.

Verreault, unpublished data; [START_REF] Benchetrit | Using otolith microchemistry to reconstruct habitat use of American eels Anguilla rostrata in the St. Lawrence River-Lake Ontario system[END_REF]. Yellow eels from SOR and SJR typically arrive at about 1 to 2 Years old (Y) and grow in running freshwater habitats, which are tributaries opening into a brackish estuary (G. Verreault, unpublished data). In France, eels were captured in two freshwater rivers, the Dordogne (DOR) and Garonne (GAR) Rivers, one brackish site, the Gironde estuary (GIR) with a salinity gradient moving up and downstream as a function of both tide and freshwater discharge [START_REF] Larrose | Assessing the current state of the Gironde Estuary by mapping priority contaminant distribution and risk potential in surface sediment[END_REF]) and one salt marsh site, Certes (CER). This site is located in the Certes marsh connected to the Arcachon basin, a nearby coastal environment not influenced by water from the Gironde watershed (Fig. 1b and Table 1). The GIR site is known to be contaminated by anthropic and agricultural activities [START_REF] Delpech | Development of a fish-based index to assess the ecological quality of transitional waters: the case of French estuaries[END_REF][START_REF] Tableau | A fish-based index of estuarine ecological quality incorporating information from both scientific fish survey and experts knowledge[END_REF]. The DOR site is considered as the least contaminated site in the watershed and the GAR site is located near the city of Bordeaux in a contaminated region [START_REF] Tapie | PBDE and PCB contamination of eels from the Gironde estuary: from glass eels to silver eels[END_REF]. The DOR and GAR sites are located in the tidally influenced part of the rivers. In contrast to A. rostrata from the St. Lawrence system, the estimated transit time of A.

anguilla elvers from the coast to their growth habitat is very short, an average of 2-3 weeks to reach the GAR and DOR sites [START_REF] Beaulaton | The efficiency of selective tidal stream transport in glass eel entering the Gironde (France)[END_REF].

European and American yellow eels were captured using trap nets (mesh size of 6.4 mm), trawl (mesh size of 31 mm) and electro-fishing by professional fishers or technical staff (Table 1). In 2011 and 2012, all eels were from a single capture in one day, except for SOR and SJR sites in Quebec where eels were captured in two days due to a small sample the first day. Based on existing information according to [START_REF] Nilo | Synthèse des connaissances et établissement d'une programmation de recherche sur l'anguille d'Amérique (Anguilla rostrata)[END_REF], yellow eels from a pre-defined size were targeted, in order to have the same range of ages (5-10 years old) for all sites. In Quebec, the selected fish ranged between 495 to 655 mm in length and in France, the selected fish ranged between 433 to 578 mm. After each capture, eels were kept in tanks with aerated water from the site for few hours until processing.

Morphometric measurements and tissue sampling

American and European eels were processed using the same standardized protocols in France and in Quebec. Total body length (± 1 mm) and total body mass (± 0.1 g) were measured and then eels were killed by decapitation. Eel manipulations were made in agreement with approved protocols from animal ethics committees in both countries. Sex (undifferentiated: U; Male: M; Female: F) was determined by macroscopic observation of the gonads using previously described criteria [START_REF] Jones | On the sex and distribution of the freshwater eel (Anguilla anguilla)[END_REF]. Based on calculated gonadal and pectoral fin indices [START_REF] Durif | Analyse des paramètres morphologiques et physiologiques lors de la préparation à la migration de dévalaison chez l'anguille européenne (Anguilla anguilla) du lac de Grand-Lieu (Loire-Atlantique)[END_REF][START_REF] Durif | Seasonal evolution and individual differences in silvering eels from different locations[END_REF]), all eels retained were immature yellow eels.

Otolith processing

Sagittal otoliths were extracted, cleaned with distilled water, dried and stored in plastic Eppendorf tubes at ambient room temperature until processing. Otoliths were processed for age determination according to Verreault et al. ( 2009) using a five-stage procedure: 1) cleaning with sodium hypochlorite (2.5%) and rinsing with demineralized water and drying with ethanol (70%); 2) embedding in epoxy resin (24h); 3) sanding on 800-grit paper and polishing with aluminium oxide powder (1µm); 4) etching with EDTA solution, and finally; 5) staining with 0.01% toluidine blue. Otoliths were observed with a binocular microscope (60X) under transmitted light for aging.

Back-calculation of length

Back calculation was undertaken using a relationship between body size and otolith radius [START_REF] Francis | Back-calculation of fish length: a critical review[END_REF]). The biological intercept procedure proposed by [START_REF] Campana | How reliable are growth back-calculations based on otoliths?[END_REF] uses a biologically-determined intercept in the back-calculation equation using the mean size of the body and the otolith at the elver stage. This method is robust to any variations in the fish-otolith relationship and provided from [START_REF] Tremblay | Reproductive strategy of female American Eels among five subpopulations in the St. Lawrence River watershed[END_REF]. Back-calculated body length (L t ) at age t are calculated as follows:

Page 8 of 37

Can. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18

For personal use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version of record.

Author-produced version of the article published in Can. J. Fish. Aquat. Sci., 2018, 75 (8), 1270-1279

The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025

L t = L c + (L c -L bi ) x (O t -O c ) / (O c -O bi )
where O t is the measured otolith radius at age t, L c and O c are the body length and otolith radius at capture, respectively, and L bi and O bi are the body length and otolith radius at the biological intercept, respectively. Measurements of the radius (mm) of the t th annuli (R t ), which is the distance from the elver mark (the first growth check outside the nucleus from where continental age determination commences) of the otolith to the t th annuli (meaning: O t -O bi ) and of the radius (mm) of the otolith at capture (R), which is the distance from the elver mark of the otolith to the edge (O c -O bi ). These measurements were always made along the longest axis on the ventral side of the sagittae of the otolith. In this study, the mean total body length of elvers when they recruit to coasts (L bi ) was set to 67 mm (SD=4.5), for A. anguilla and 62 mm (SD=3.2) for A. rostrata, which corresponds to the age 0 year old or continental age (ICES 2009;[START_REF] Lambert | Variabilités morphologique et du taux d'hormone de croissance des civelles d'anguilles européennes dans l'estuaire de la Gironde au cours de la saison 1997-1998[END_REF].

Statistical analyses

Statistical analyses were conducted using Prism v. 5.0 (GraphPad Software Inc., San Diego, CA, USA) and R (R. 2011). Mixed models were fitted using the package lme4 [START_REF] Horton | A multistate capture-recapture modeling strategy to separate true survival from permanent emigration for a passive integrated transponder tagged population of stream fish[END_REF]) and Tukey tests were carried out using package multcomp [START_REF] Bates | Fitting linear mixed-effects models using lme4[END_REF]. For all statistical tests, a probability of p≤0.05 was considered significant. It was not possible to collect same-sized eels at each sampling site. To increase the comparability of eels among sampling sites, only eels between 2 and 20 years old were considered in this analysis. The total number of eels retained for statistical analyses were 102 in 2011 and 112 in 2012 (Tables 2 and3).

Most analyses were conducted for each country/species separately. Firstly, biological characteristics of eels at capture were statistically compared among sites within a species and year. For each year, data normally distributed (using the Shapiro-Wilk test) were compared with one-way ANOVA followed by a Tukey's multiple comparison test. Data not normally distributed was compared among sites using non-parametric Kruskal-Wallis test (KW) followed by the Dunn's multiple comparison test. The relative proportions of female and undifferentiated eels were compared among sites or between sampling years in each country using a Chi-square test.

To compare size among sites within systems for eels of the same age, back-calculated lengths from 1 to 5 Years old (Y) were compared among sites within a species using the KW test as described above. Then, we carried out linear mixed effect models for each site including random effects to account for possible inter-cohort variations in length and individual random effects. First, we fitted a linear mixed effect model for each site to assess whether length at 1Y significantly differed among sites:

(1) L 1 ~Site + (1|cohort)

with L 1 , the length-at-age at 1Y, Site, the site of capture and (1|cohort), a random effect to account for possible inter-cohort variations in length. Tukey multiple comparison tests were then carried out to detect significant differences among sites.

A second analysis was performed to check whether growth rates were different among sites. We first computed yearly growth increments for each fish, i.e. the difference in backcalculated lengths between two successive ages. Growth increments are a linear function of initial length when fish follow a Von Bertalanffy growth curve (i.e.

I t =L t -L t-1 =(L ∞ -L t-1 )•(1-exp(-K))
with I t , the growth increment between year t-1 and year t, Lt the size-at-age at age t, K the Brody coefficient and L ∞ the asymptotic length). In view of this, we fitted a generalized linear mixed model on growth increments for each species separately:

(2) I t ~Site+L t-1 +(1|eel) +(1|cohort) with I t , the vector of yearly growth increments (we used the 4 increments for each from 1Y to 5Y) and L t-1 , the vector containing the corresponding initial length. We added a cohort random effect to account for inter-cohort differences and an individual random effect since each eel was considered 4 times in the analysis. We used a gamma family distribution to account for dissymmetry in growth increment distributions [START_REF] Delong | Estimating growth and mortality of juvenile winter flounder, Pseudopleuronectes americanus, with a length-based model[END_REF][START_REF] Drouineau | Assessing the impact of discretisation assumptions in a length-structured population growth model[END_REF], with an identity link to preserve the linear relationship. Similarly, to the first model, we carried out Tukey multiple tests to compare yearly growth increments among sites.

Results

Eel biological characteristics in the St. Lawrence and Gironde systems

For A. anguilla sampling, male eels were removed due to small sample size (n=3 in the Gironde system) and for A. rostrata sampling, no males have been detected (n=0 in the St.

Lawrence system). In 2011 and 2012, body length and age of A. anguilla from the Gironde system ranged from 278 to 825 mm and from 3 to 18 years old respectively, and differed among sites (H=38.5, p<0.0001 in 2011 and H=12.8, p=0.005 in 2012 for body length; H=33.8, p<0.0001 in 2011 and H =25.3, p<0.0001 in 2012 for age, Tables 2 and3). Eels from DOR were smaller and lighter than eels captured in GAR, GIR and CER with median body length and body mass 1.3 and 2.8 times lower, respectively (H=37.7, p<0.0001 in 2011 and H=21.8, p<0.0001 in 2012 for body mass, Tables 2 and3). In both sampling years, the proportions of undifferentiated eels increased from the salt marsh habitat (CER) to freshwater habitats and the highest proportions of undifferentiated eels were observed in DOR and GAR (Tables 2 and3). No significant difference was observed in biological characteristics between female and undifferentiated European eels in GAR and DOR (Student t-test, Tables S2).

In the St. 2 and3). Eels from SOR were consistently smaller and lighter than eels from SPL and SFL (H=44.4, p<0.0001 in 2011 and H=40.7, p<0.0001 in 2012, for body mass, Tables 2 and3). In 2011, body length and age were higher in eels from SPL and SFL compared to SJR (Table 2). In 2011 and 2012, body mass was higher in SPL and SFL than in SJR (Tables 2 and3).

In 2011, yellow eels captured in the Quebec sites were older than those sampled in France (U=357.5, p<0.0001, Mann Whitney test).

Comparison of back-calculated length at age 1Y to 5Y of eels from different habitats in each species

In the Gironde system, at 1Y, back-calculated lengths of eels sampled in the salt marsh CER site (median of 219 mm) were higher than those of eels sampled in the estuarine GIR site and in the fluvial GAR and DOR sites (medians were respectively 152, 136 and 134 mm, H=51.4, p<0.0001, Fig. 2a and Table S3). At 2Y, back-calculated lengths of eels from the estuarine GIR site became higher than those of eels from fluvial GAR and DOR sites and the difference was maintained until 5Y (F=51.2, p<0.0001 for 2Y; F=48.7, p<0.0001 for 3Y; F=30.6, p<0.0001 for 4Y and F=16.7, p<0.0001 for 5Y, Fig. 2a). From 2Y to 5Y, eels from the marine CER site had the higher median of back-calculated length compared to the other sites and eels from GAR and DOR sites remained with the lowest median of back-calculated length (Fig. 2a).

In the St. Lawrence system, at 1Y, back-calculated lengths of eels captured in fluvial SPL and SFL sites (with a median length of 184 and 171 mm respectively) were higher than those of eels from the estuarine SOR and SJR sites (with a median length of 134 and 142 mm respectively, (F=31, p<0.0001, Fig. 2b and Table S3). The significant difference between higher back-calculated length of eels from SPL and SFL and lower back-calculated length of eels from SOR and SJR was maintained until 5Y (H=46.4, p<0.0001 for 2Y; H=50.2, p<0.0001 for 3Y;

F=31.2, p<0.0001 for 4Y and F=40.6, p<0.0001 for 5Y, Fig. 2b and Table S3).

Linear mixed effect models

In both years, for A. anguilla, this analysis confirmed that the length at age 1Y was higher in CER than in any other site. Length at age 1Y in GIR was slightly higher than in DOR, and no difference was observed between GAR and DOR. This confirmed the decreasing gradient in length at age 1Y from downstream to upstream habitats (Table 4). For A. rostrata, length at age 1Y was higher in the two upstream lakes than in downstream rivers but no differences were detected between the two lakes or between the two rivers (Table 4).

Results were rather similar regarding subsequent growth-increments (Table 5): increments were higher in downstream sites than in upstream sites for A. anguilla, although CER and GIR did not differ. Yearly growth increments from 1Y to 5Y were consistently higher in the two upstream lakes than in downstream rivers for A. rostrata.

Discussion

This study confirms that the current spatial pattern of early growth rates of yellow A.

anguilla in the Gironde system is consistent with the paradigm of higher early growth rates of eels growing in habitats more saline and/or less distant from the Sargasso Sea spawning area.

Moreover, comparisons with historical data in the Gironde system indicate that this spatial pattern remains stable despite temporal variation in the overall growth rates. In contrast, our study reveals that early growth rates of yellow A. rostrata in the St. Lawrence system are highest in eels captured in the less saline habitats, more distant from the Sargasso Sea, in contradiction with the current paradigm. This new information is useful to guide further studies on the mechanisms leading to the different ecotypes of Atlantic eels in different rearing habitats.

Spatial variations in early back-calculated length in Anguilla anguilla

At 1Y, eels from the salt marsh site CER had higher back-calculated length than eels from the estuarine site GIR and from the fluvial sites GAR and DOR. From 2Y to 5Y, while eels from CER maintained their higher back-calculated length, eels from the estuarine site GIR displayed an increasingly faster growth rate with higher back-calculated length than those of eels from freshwater sites. These results are consistent with previous studies, which demonstrated that eels sampled in brackish habitats grew faster than those from riverine freshwater habitats.

Specifically, [START_REF] Daverat | Tactics and demographic attributes in the European eel Anguilla anguilla in the Gironde watershed, SW France[END_REF] have reached similar conclusions on the spatial variations of early growth rates of eels sampled in 2001-2004 in the Gironde system. Since the transit time of A. anguilla glass eels from the coast to their growth habitat was estimated to last 2-3 weeks according to [START_REF] Beaulaton | The efficiency of selective tidal stream transport in glass eel entering the Gironde (France)[END_REF], it is likely that their growth was already influenced by the characteristics (biological, chemical and/or physical) of the selected habitat at the age of 1Y. However, laboratory [START_REF] Edeline | Dispersal in metamorphosing juvenile eel Anguilla anguilla[END_REF]) and genetic studies [START_REF] Pujolar | Growth rate correlates to individual heterozygosity in the European eel, Anguilla anguilla L[END_REF] in A. anguilla suggest that the ecotypic differences among eels occupying different habitats are not only the effect of phenotypic plasticity. They may also be caused by functional genetic differences resulting from intra-generational spatially varying selection and/or genotypedependent habitat choice (or both) of ecologically divergent habitats [START_REF] Pavey | RAD Sequencing Highlights Polygenic Discrimination of Habitat Ecotypes in the Panmictic American Eel[END_REF].

At 1Y, no significant difference in length was detected between eels from the estuarine site GIR and eels from the freshwater sites GAR and DOR. This observation could be related to a high within-site inter-individual variability in back-calculated length. One factor of variability is the older age of some of the eels sampled in GIR and GAR. However, results of analyses using a linear mixed model integrating individual variability and cohort effect were consistent with the results of the KW test indicating these sources of variability had no significant impact on the results. Overall these data and the similarity with the pattern described a decade earlier by [START_REF] Daverat | Tactics and demographic attributes in the European eel Anguilla anguilla in the Gironde watershed, SW France[END_REF] indicate a very consistent spatial pattern of higher growth at an early age (from 1Y to 5Y) in yellow eels sampled in the more saline site (CER) compared to eels sampled in the upstream fluvial section of the Gironde system.

Spatial variations in early back-calculated length in Anguilla rostrata

Anguilla rostrata sampled in the fluvial freshwater sites (SPL and SFL) in the Upper St.

Lawrence River had a higher back-calculated length from 1Y to 5Y compared to same-aged A.

rostrata sampled downstream in tributaries opening into a brackish estuary (SJR and SOR). In our study, the difference in back-calculated length in A. rostrata sampled in the upstream St.

Lawrence system and those sampled in downstream tributaries was present as early as 1Y. At this early age, A. rostrata of the fast-growing group (LSP, LSF)

have not yet reached their growing habitat. It is likely that they were still migrating in the St. Lawrence system and therefore their growth rate could not yet be influenced by the local conditions in the fluvial lakes [START_REF] Benchetrit | Using otolith microchemistry to reconstruct habitat use of American eels Anguilla rostrata in the St. Lawrence River-Lake Ontario system[END_REF]).

The higher back-calculated length at 1Y of eels from the upstream portion of the St.

Lawrence system compared to more downstream populations that we report here contrasts with recent literature. Based on previous studies on A. rostrata captured in the Gulf of St. Lawrence and in tributaries of the marine portion of the St. Lawrence system, a slower growth at 1Y was expected with increasing distance of the capture sites from the Sargasso Sea [START_REF] Boivin | How does salinity influence habitat selection and growth in juvenile American eels Anguilla rostrata?[END_REF].

A common garden experiment by [START_REF] Côté | Differential effects of origin and salinity rearing conditions on growth of glass eels of the American eel Anguilla rostrata: implications for stocking programmes[END_REF] suggested a genetic basis for growth differences among A. rostrata early life stages from different geographical locations. They found that A. rostrata elvers from the Mira River (Nova Scotia; shortest distance to the Sargasso Sea) grew faster than those from Grande-Rivière-Blanche (Québec, longest distance to the Sargasso Sea), a result confirmed by the study of [START_REF] Boivin | How does salinity influence habitat selection and growth in juvenile American eels Anguilla rostrata?[END_REF]. Molecular tools were used to support the hypothesis that life-history traits of A. rostrata were influenced by spatially variable genetic selection along latitudinal gradients of temperature and/or increasing distance from the Sargasso Sea and by genotype-dependent habitat choice [START_REF] Boivin | How does salinity influence habitat selection and growth in juvenile American eels Anguilla rostrata?[END_REF][START_REF] Côté | Differential effects of origin and salinity rearing conditions on growth of glass eels of the American eel Anguilla rostrata: implications for stocking programmes[END_REF][START_REF] Gagnaire | The genetic consequences of spatially varying selection in the panmictic American eel (Anguilla rostrata)[END_REF][START_REF] Pavey | RAD Sequencing Highlights Polygenic Discrimination of Habitat Ecotypes in the Panmictic American Eel[END_REF]. Moreover, glass eels of A. rostrata from New Brunswick and Nova Scotia sites (which are closer to the Sargasso Sea) translocated to Lake Ontario exhibited faster annual growth and matured at smaller sizes and ages at reproductive migration than wild migrating silver eels from Lake Ontario [START_REF] Couillard | Assessment of Fat Reserves Adequacy in the First Migrant Silver American Eels of a Large-Scale Stocking Experiment[END_REF][START_REF] Stacey | A caution for conservation stocking as an approach for recovering Atlantic eels[END_REF]Verreault et al. 2010). Despite some limitations in fish samplings, our study clearly demonstrates an inverse pattern of early growth rate from 1Y to 5Y for A. rostrata with fast growing eels reaching the upstream St. Lawrence system and slow growing eels remaining in the downstream section. These observations question our current understanding of eels adaptation to environmental heterogeneity [START_REF] Mateo | Cause or consequence? Exploring the role of phenotypic plasticity and genetic polymorphism in the emergence of phenotypic spatial patterns of the European eel[END_REF].

Similarities between A. anguilla and A. rostrata and comparisons with previous studies

This study is the first to use a common methodology to compare the back-calculated length of the yellow stage of the two species of Atlantic eels. Given the wide range of habitats covered in this study (a larger hydro-geographical scale associated with longer upstream migration of early life stages, longer distance to the Sargasso Sea for Gironde system eels and different temperature regimes for the two river systems), we found that from 1Y to 5Y, the backcalculated length of A. rostrata collected in the St. Lawrence system in the northern part of their distribution range presented some similarities with values from A. anguilla collected in the Author-produced version of the article published in Can. J. Fish. Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025

Gironde system. The back-calculated length of A. rostrata from the upper St. Lawrence system (SPL and SFL) was similar to that of eels from the downstream sites GIR (at 2Y and 3Y) and CER (at 4Y and 5Y). From 1Y to 5Y, the back-calculated length of eels from the downstream St.

Lawrence system sites, SJR and SOR, were comparable to those of slow-growing eels from the upstream sites GAR and DOR. Our study highlighted an inverse pattern of early growth rate according to the upstream-downstream gradient of the location of habitats between A. rostrata with faster growers from the upstream section of the St. Lawrence system and A. anguilla with faster growers from the downstream section of the Gironde system.

In our study, median of back-calculated length at 1Y and 5Y of A. anguilla collected in the Gironde in 2011-2012 (Table S3) were approximately twice those observed in eels captured 10 years earlier in the same sites in the Gironde system by [START_REF] Daverat | Tactics and demographic attributes in the European eel Anguilla anguilla in the Gironde watershed, SW France[END_REF]. This suggests an improvement in growth performance possibly related to abiotic factors (salinity, temperature, productivity and habitat quality) and/or biotic factors (low density of eels, weak coastal recruitment of glass eels) [START_REF] Bonhommeau | Impact of climate on eel populations of the Northern Hemisphere[END_REF][START_REF] Daverat | One century of eel growth: changes and implications[END_REF][START_REF] Knights | A review of the possible impacts of long-term oceanic and climate changes and fishing mortality on recruitment of anguillid eels of the Northern Hemisphere[END_REF].

Despite the observed changes in growth performance in A. anguilla, the spatial pattern of early growth rate did not change with time. However, the back-calculated length at 5Y of slow growing A. rostrata collected in SJR in 2011-2012 are similar to those calculated for yellow eels captured at the same site in 2004 by [START_REF] Thibault | Facultative catadromy in American eels: testing the conditional strategy hypothesis[END_REF]. These data suggest that, at least for the SJR site, no significant change in growth rate up to 5Y has occurred in the last 7-8 years. Unfortunately, to our knowledge, no historic data for early GR of yellow eels collected in the upper part of the St.

Lawrence system is available for comparison.

Limitations and future directions

The unique biology of eels, and especially their panmictic nature, the variability in their migratory patterns during their growth phase and their large distribution ranges from tropical to subarctic climates, makes their study a challenge in any circumstance. Population declines and severe restrictions on commercial and scientific collection impose further difficulties. Sampling sites in this international study were selected according to salinity, metal and organic contamination gradients found in sediments and biota along the Saint Lawrence and the Gironde systems [START_REF] Laporte | RAD sequencing reveals within generation polygenic selection in response to anthropogenic organic and metal contamination in North Atlantic Eels[END_REF]). Nevertheless, these sites cannot be considered to represent fairly the wide range of habitats that the species occupy. Some limitations related to the sampling of eels (30 eels captured in each site and in two years and the use of different fishing gear adapted to each habitat) and the wide range of estimated age (2 to 18 years) may constitute biases.

In addition to spatially variable selection and/or genotype-dependant habitat choice, several other processes could contribute to the observed growth variability in 1Y-old eels, including epigenetic changes, early environment influences with carry-over effects on growth and interactions between genotype and environment through phenotypic plasticity. A reduction in the abundance of glass eels arriving in continental habitats associated with the low density of eels in their growth habitats could have contributed to increase their length at early age [START_REF] Castonguay | Why is recruitment of the American eel, Anguilla rostrata, declining in the St. Lawrence River and Gulf?[END_REF][START_REF] Moriarty | Management of the European eel[END_REF]. Temporal variations in hatching times [START_REF] Kuroki | Hatching time and larval growth of Atlantic eels in the Sargasso Sea[END_REF] and in environmental conditions during development [START_REF] Jonsson | Early environment influences later performance in fishes[END_REF]) could also affect early growth rates and behaviour. The season in which glass eels arrived on the continental coasts may also affect their size [START_REF] Desaunay | Seasonal and long-term changes in biometrics of eel larvae: a possible relationship between recruitment variation and North Atlantic ecosystem productivity[END_REF]. Local environmental factors such as tidal height, water temperature, water discharges and the difference in temperatures between estuarine and river habitats are known to have effects on glass eel's abundance and upstream migratory behavior of migrants [START_REF] Jessop | Annual variability in the effects of water temperature, discharge, and tidal stage on the migration of American eel elvers from estuary to river[END_REF][START_REF] Lambert | Effects of body condition and water temperature on Anguilla anguilla glass eels migratory behavior[END_REF]). In addition, physiological traits of eels such as body condition, endocrine system (thyroid hormones) and health status may affect the locomotor activity, the migratory behavior and growth performance (Edeline et al.

Page 18 of 37 Can. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18

For personal use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version of record.

Author-produced version of the article published in Can. J. Fish. Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025

2005b; [START_REF] Lefebvre | Is the continental life of the European eel Anguilla anguilla affected by the parasitic invader Anguillicoloides crassus?[END_REF]. Further studies are needed on the interactions between genetic variation and phenotypic plasticity in response to environmental variation during early development of Atlantic eels and their impact on migration patterns and later-life history performance. For that purpose, longer term studies focussing on variations of growth and behaviour of 0-1 Y eels in relation to genotypes and/or environmental factors are recommended, with larger sample sizes and greater number of sites.

In spite of these limitations, our study provides a new indication that for A. rostrata sampled along the St. Lawrence system, faster growing individuals appear more likely to colonize the head of the system than slower growing individuals. We have also highlighted that the median of back-calculated length of A. anguilla at 1Y and 5Y is twice higher than that reported 10 years earlier. Further studies and monitoring are clearly needed to test the hypothesis of a possible increase of early growth rate in relation to the modifications of abiotic and biotic factors and its consequences for age, length and lipid reserves at migration in A. anguilla and A.

rostrata.

Finally, whereas several studies have investigated the relationships between genetic selection and growth rate characteristics of glass eels in A. rostrata, these studies are limited to eels from tributaries of the lower St. Lawrence Estuary and Gulf and there is a need to confirm these relationships in A. anguilla. Comparative studies of these two related Atlantic eel species would help to better understand the relationship between differential growth rate, size-atmetamorphosis, lipid accumulation, sexual differentiation and their genetic and environmental bases in these two species. This information would facilitate prediction of the effects of global environmental changes on the recruitment of these two declining species. 2 Distance between the sampling sites and the Sargasso Sea spawning ground were estimated via the Google Earth web tool. For personal use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version of record.

Author-produced version of the article published in Can. J. Fish. Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025

  Freshwater; BW = Brackish water. Salinity values of each sampling site are detailed inBaillon et al. (2015).

a

  Reference: Tapie et al. (2011); b Reference: Jessop (2010).

  . Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18

Figure 2 .

 2 Figure 2. Back-calculated length in relation to age from 1Y to 5Y for a) Anguilla anguilla

  Can. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18For personal use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version of record. -produced version of the article published in Can. J. Fish.Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025

	Author
	Page 10 of 37

  Lawrence system, the body length and age of A. rostrata captured in 2011 and 2012 ranged from 265 to 963 mm and from 5 to 19 years old respectively and both parameters Can. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18For personal use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version of record. -produced version of the article published in Can. J. Fish.Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025 differed among sites (F=174.0, p<0.0001 in 2011 and H=41.5 p<0.0001 in 2012, for body length and F=13.1, p<0.0001 in 2011 and H=6.1, p<0.0001 in 2012 for age, Tables

	Author
	Page 11 of

  Page 12 of 37 Can. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18For personal use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version of record. -produced version of the article published in Can. J. Fish.Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025

	Author

  Can. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18For personal use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version of record. -produced version of the article published in Can. J. Fish.Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025

	Author
	Page 13 of

  Can. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18For personal use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version of record.

	Author-produced version of the article published in Can. J. Fish. Aquat. Sci., 2018, 75 (8), 1270-1279
	The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025	Page 14 of 37

  Can. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18For personal use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version of record.

	Page 16 of 37

Table 1 .

 1 Characteristics of sampling sites in France and in Quebec.

			FRANCE (A. anguilla)			QUEBEC (A. rostrata)
	Site	DOR	GAR	GIR	CER	SFL	SPL	SOR	SJR
	Latitude (°N)	44°54'30" 44°43'51"	45°12'23" 44°41'18"	48°51'40"	48° 21' 17"	46°19'50" 45°09'18"
	Longitude (°W)	0°15'01"	0°28'05"	0°43'13"	1°1'39"	64°28'47" 68°45'28"	72°32'06" 74°23'04"
	Sampling
	period (month
	and years)
	Fishing gear	Trap nets Trap nets	Trawl	Trap nets	Trap nets	Electro-fishing	Trap nets Trap nets
	Habitat 1 (Salinity)	FW	FW	BW	BW	FW	FW	BW	FW
	Mean length at silvering (mm)	752

Table 2 .

 2 Biological characteristics for Anguilla anguilla and Anguilla rostrata for each sampling site in Can. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18 al use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version oAuthor-produced version of the article published in Can. J. Fish.Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025

													Page 30 of 37
	2011. Values represent the median (Q1, Q3) of each biological variable (BL: Body Length; BM: Body
	Mass. Sex ratios are presented in % of the total number of eels within a site (Chi-square test). Different
	superscript lower case letters indicate statistically significant differences among sites within species
	(ANOVA or Kruskall Wallis).
										2011
				A. anguilla					A. rostrata
	Site	DOR	GAR	GIR	CER	SFL	SPL	SOR	SJR
	N total
	(Females and	13	14	15	11	15	7	15	15
	undifferentiated)
	Sex ratio (%):
	Female	62	86	100	100	100	100	100	100
	Undifferentiated	38	a	14	ab	0	b	0	b	0	0	0	0
	BL (mm)	325 (302, 354) c	434 (393, 503) b	565 (530, 614) a	440 (395, 478) b	832 (801, 856) a	894 (876, 929) a	323 (290, 382) c	b (465, 606) 535
	BM (g)	51 (39, 69) c	127 (103, 202) b	300 (249, 424) a	128 (102, 167) b	1254 (1147, 1298) a	1565 (1435, 1642) a	50 (32, 92) c	b (172, 440) 275
	Age (year)	5 (4, 5) b	6 (5, 7) b	8 (8, 12) a	5 (4, 5) b	12 (12, 13) a	13 (11, 16) a	8 (7, 10) b	a (8, 12) 10

Table 3 .

 3 Biological characteristics for Anguilla anguilla and Anguilla rostrata for each sampling site in 2012. Values represent the median (Q1, Q3) of each biological variable (BL: Body Length; BM: Body Mass). Sex ratios are presented in % of the total number of eels within a site (Chi-square test). Different superscript lower case letters indicate statistically significant differences among sites within species (ANOVA or Kruskall Wallis). Can. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18 al use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version oAuthor-produced version of the article published in Can. J. Fish.Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025

	2012

Table 4 :

 4 Tukey multiple comparisons tests on the site effects in the linear

	mixed models on length-at-age 1. Statistical significance codes: ≤ 0.00001
	'***'; ≤ 0.001 '**'; ≤ 0.01 '*'; ≤ 0.05 '.'; ≤ 0.1 'ns'
	Species	Site 1	comparison	Site 2	significance level
	A. anguilla
		CER	>	GIR	***
		CER	>	DOR	***
		CER	>	GAR	***
		GIR	>	GAR	.
		GIR	>	DOR	*
		GAR	>	DOR	ns
	A. rostrata
		SFL	>	SPL	ns
		SFL	>	SJR	***
		SFL	>	SOR	***
		SPL	>	SJR	***
		SPL	>	SOR	***
		SJR	>	SOR

Table 5 :

 5 Tukey multiple comparisons tests on the site effects inCan. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18For personal use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version of record. -produced version of the article published in Can. J. Fish.Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025 Figure 1. Study areas and sampling sites for (a) yellow eels captured in Quebec (North Eastern Canada) and (b) yellow eels captured in the Gironde region (South West France)showing the separation into 3 salinity sectors (freshwater, brackish and marine sectors). Black squares refer to sampling sites and black circles refer to the closest cities. The symbols (★)indicate the limits of tidal influence. The map was made using QGIS project, a free and open source geographic information system (www.qgis.org).Can. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18For personal use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version of record.Author-produced version of the article published in Can. J. Fish.Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025

	the linear mixed models on growth increments. Statistical
	significance codes: ≤ 0.00001 '***'; ≤ 0.001 '**'; ≤ 0.01 '*'; ≤
	0.05 '.'; ≤ 0.1 'ns'
	Species	Site 1	comparison Site 2	significance level
	A. anguilla
		CER	>	GIR	ns
		CER	>	DOR	***
		CER	>	GAR	.
		GIR	>	GAR	*
		GIR	>	DOR	**
		GAR	>	DOR
	A. rostrata
		SFL	<	SPL	ns
		SFL	>	SJR	***
		SFL	>	SOR	***
		SPL	>	SJR	***
		SPL	>	SOR	***
		SJR	>	SOR	.

Author-produced version of the article published in Can. J. Fish.Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025

Acknowledgments

Page 19 of Can. J. Fish. Aquat. Sci. Downloaded from www.nrcresearchpress.com by CEMAGREF on 01/17/18 For personal use only. This Just-IN manuscript is the accepted manuscript prior to copy editing and page composition. It may differ from the final official version of record. The authors wish to acknowledge all the personnel of the INRS-ETE, UMR EPOC (university of Bordeaux) and IRSTEA (Bordeaux-Cestas) who participated in the field sampling. We would like to thank the technical staff from the Ministère des Forêts, de la Faune et des Parcs for otolith expertise and for providing us with back-calculated length data. This work was supported by the joint ANR Blanc (ANR-10-INTB-1704)-NSERC SGP grant IMMORTEEL and by the Aquitaine region (GAGILAU and MICANG projects). The laboratory work of Géraldine Patey was in part supported by Fisheries and Oceans Canada.

June 2011 May 2012 June 2011 May 2012 June 2011 May 2012 May 2011 May 2012 June 2011 May 2012 June 2011 June 2012 July 2011 June 2012 June 2011 June 2012

Author-produced version of the article published in Can. J. Fish. Aquat. Sci., 2018, 75 (8), 1270-1279 The original publication is available at http://www.nrcresearchpress.com doi : 10.1139/cjfas-2017-0025

Author-produced version of the article published in Can. J. Fish. Aquat. Sci., 2018, 75 (8)


Arturo Laurenzi

Enrico Mingo Hoffman

Matteo Parigi Polverini

Nikos G Tsagarakis

Balancing Control through Post-Optimization of Contact Forces

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Balancing represents a crucial requirement for legged robots expected to cope with a variety of unstructured terrains and environments. Despite a deep knowledge on the dynamics governing the balancing of articulated bodies, balancing is still considered a challenging control problem, especially when dealing with torque controlled humanoids and legged robots in general. In this respect, the existing control approaches are generally classified in two categories.

The first category achieves balancing through a two-stage methodology. An optimal contact force distribution problem is first solved with respect to the robot centroidal dynamics. This phase will be hereafter referred to as pre-optimization of contact forces. The computed contact forces are then mapped to joint torques under quasi-static assumptions, see Hyon et al. [START_REF] Hyon | Full-body compliant humanhumanoid interaction: Balancing in the presence of unknown external forces[END_REF], [START_REF] Hyon | Compliant terrain adaptation for biped humanoids without measuring ground surface and contact forces[END_REF] and Ott et al. [START_REF] Ott | Posture and balance control for biped robots based on contact force optimization[END_REF], [START_REF] Henze | An approach to combine balancing with hierarchical whole-body control for legged humanoid robots[END_REF], [START_REF] Henze | Multi-contact balancing of humanoid robots in confined spaces: Utilizing knee contacts[END_REF], or through inverse dynamics, as in [START_REF] Lee | Ground reaction force control at each foot: A momentum-based humanoid balance controller for non-level and non-stationary ground[END_REF], [START_REF]A momentum-based balance controller for humanoid robots on non-level and non-stationary ground[END_REF] and [START_REF] Stephens | Dynamic balance force control for compliant humanoid robots[END_REF].

In opposition to the first category, the humanoid balancing problem can be alternatively addressed in a single-stage fashion by entirely exploiting the full-body inverse dynamics. Several inverse dynamics controllers, e.g. [START_REF] Khatib | Whole-body dynamic behavior and control of human-like robots[END_REF], [START_REF] Sentis | Synthesis of whole-body behaviors through hierarchical control of behavioral primitives[END_REF], [START_REF] Sentis | Compliant control of multicontact and center-of-mass behaviors in humanoid robots[END_REF], compute joint torques by modeling contacts as rigid constraints and projecting the dynamics into a constraint free space. In this way an explicit solution of the contact force distribution problem is not required, although optimality of the problem is not guaranteed. Nevertheless, Righetti et al. [START_REF] Righetti | Inverse dynamics control of floating-base robots with external constraints: A unified view[END_REF], [START_REF] Righetti | Optimal distribution of contact forces with inverse-dynamics control[END_REF] showed that is possible to design inverse dynamics controllers and operational space controllers that are optimal with respect to any combination of linear and quadratic 1 Advanced Robotics Department (ADVR), Istituto Italiano di Tecnologia, Genova, Italy {arturo.laurenzi, enrico.mingo, matteo.parigi, nikolaos.tsagarakis}@iit.it cost in the contact forces and in the torque commands. On the other hand, methods based on a hierarchical Quadratic Programming (QP) formulation of the full-body inverse dynamics, see [START_REF] Mansard | A dedicated solver for fast operational-space inverse dynamics[END_REF], [START_REF] Saab | Dynamic whole-body motion generation under rigid contacts and other unilateral constraints[END_REF], [START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF], and [START_REF] Herzog | Balancing experiments on a torque-controlled humanoid with hierarchical inverse dynamics[END_REF], [START_REF] Herzog | Momentum control with hierarchical inverse dynamics on a torque-controlled humanoid[END_REF], explicitly consider contact forces as variables for the resulting optimization problem.

It is worth pointing out that, despite a clear advantage in terms of required computation time of the first category of methods over the second, the pre-optimization of contact forces raises a major concern about the role of the momentum of rotation in balance control. It is known in fact that the kinetic momentum of rotation is not directly related to the actual orientation of an articulated system [START_REF] Wieber | Holonomy and nonholonomy in the dynamics of articulated motion[END_REF], [START_REF] Wieber | Modeling and control of legged robots[END_REF]. As a consequence, controlling the momentum of rotation for a balancing task may end up in a body rotation which is incompatible with the task itself.

In a recent work by the authors [START_REF] Hoffman | Multi-priority cartesian impedance control based on quadratic programming optimization[END_REF] a prioritized Cartesian impedance controller for redundant fixed-base robots has been proposed, based on a hierarchical QP formulation. Extending this approach to floating-base legged robots to subsequently address the humanoid balancing problem is the main focus of the present paper. The proposed balancing controller belongs to the first category of methods, i.e. it adopts a two-stage technique. In this respect, aiming to overcome the main limitation inherent in the pre-optimization of contact forces, i.e. the control of the momentum of rotation, we propose to treat a legged robot as a fully actuated fixed-base system in order to compute the desired joint torques according to [START_REF] Hoffman | Multi-priority cartesian impedance control based on quadratic programming optimization[END_REF]. Only at this point, the underactuated component of the obtained torques can be mapped into contact forces through an optimal distribution problem, hereafter referred to as post-optimization of contact forces. This way the tricky specification of a desired momentum of rotation is circumvented. In addition to the aforementioned advantage over pre-optimization methods and the reduced number of decision variables compared to single-stage methods, the add-on nature of the proposed control approach entail no modification of the original algorithm in [START_REF] Hoffman | Multi-priority cartesian impedance control based on quadratic programming optimization[END_REF].

The paper is organized as follows. The proposed torque control approach to humanoid balancing is introduced and discussed in Sec. II. Simulation results performed on both CENTAURO and COMAN+ humanoid platforms are presented in Sec. III-A, together with preliminary experimental results. Finally, concluding remarks and future work directions can be found in Sec. IV.

II. MATHEMATICAL FORMULATION

This section introduces the mathematical formulation of our prioritized force controller. In Section II-A we introduce the notation, and we give some brief background on our previous work [START_REF] Hoffman | Multi-priority cartesian impedance control based on quadratic programming optimization[END_REF]. Then, in Section II-B we present the mathematical modeling for the floating-base case; we highlight the major differences that forbid direct use of our previous work. Finally, in Section II-C we present the main contribution of this paper, which is the extension to the floating-base case.

A. The fixed-base case

Let us consider a generic task x ∈ R m expressed in terms of a function x = f (q); the corresponding task velocity is obtained by differentiation:

ẋ = J (q) q, (1)
where J (q) ∈ R m×n is the task Jacobian matrix. The dependence on q will be omitted in the following for brevity. We are interested in finding joint torques such that our task shows a desired behaviour with respect to reference signals, and external perturbations as well. To this aim, let us introduce the robot dynamics in contact with the environment:

B(q)q + h(q, q) = τ + J T F x ; (2)
in ( 2), B ∈ R n×n represents the joint-space inertia matrix, h ∈ R n is the vector of bias torques which are needed to produce zero joint acceleration, and τ ∈ R n is the vector of joint torques. External wrenches F x ∈ R 6 applied to our task are also considered. By further differentiation of (1) and substituting (2) we obtain the task dynamics as follows [START_REF] Siciliano | Robotics: modelling, planning and control[END_REF]:

Λ(q)ẍ + µ(q, q) = F τ |x + F x (3)
Where:

• Λ ∈ R m×m is the task-space inertia matrix; it is given by

Λ = J B -1 J T -1 ; (4)
• µ ∈ R m is the bias force vector (i.e. the force required to produce zero task acceleration); • F τ ∈ R m is the task space force that is produced by the joint torques; its expression is the following:

F τ |x = J T τ , (5)
where J ∈ R n×m is known as the dynamicallyconsistent pseudo-inverse of J , which has the following expression:

J = B -1 J T Λ. (6)
From (3) it can be seen that an arbitrary task behaviour can be achieved through the term F τ , which is controlled through proper choice of the joint torques according to [START_REF] Henze | An approach to combine balancing with hierarchical whole-body control for legged humanoid robots[END_REF]. Moreover, for a redundant robot the matrix J T has a non-empty null space, enabling a prioritized multiple task formulation. Finally, the prioritized formulation is solved as cascade of QPs of the following form:

min τi A i τ i -b i 2 + τ i 2 s.t. b l ≤ D τ i ≤ b u u l ≤ τ i ≤ u u A i-1 τ * i-1 = A i-1 τ i . . . A 1 τ * 1 = A 1 τ i , (7)
for values of i ranging from i = 1 to the number of priority levels. Hierarchical force control is obtained by choosing

A i = J i B -1 and b i = Λ -1 i F i .
By using this formulation, it was possible to implement mixed a stiff/compliant behaviour with priority enforcement on the bi-manual upper body of our CENTAURO and CO-MAN+ robots.

B. The floating-base case

The main focus of the present work is the extension of the algorithm that was presented in the previous section to the floating-base case. From a modeling point of view, a legged robot shows two structural differences with respect to a fixed-base one:

• under-actuation: floating-base robots can be described in terms of n degrees of freedom (one for each joint) plus six additional coordinates describing the pose of some robot link w.r.t. to an inertial world frame. Such a link is usually called floating-base. These additional DoFs are usually modelled by introducing a virtual sixdof chain of passive (unactuated) joints, which is known as virtual chain. • Contact forces: legged robots must always interact with the environment in order to be controlled in their full (6 + n)-dimensional coordinate space. Indeed, from the so-called centroidal dynamics equation, we know that the global motion of the robot is entirely given by the contact forces. The virtual chain formulation allows to easily extend the fixed-base dynamics (2) to the floating-base case. We just augment the generalized coordinate vector with six virtual joints as in the following equation1 :

q = q u q a , (8)
where q ∈ R 6+n is obtained by stacking the configuration vector of virtual joints q u ∈ R 6 with the one corresponding to the n actuated joints q a ∈ R n . Then, the floating-base dynamics equation is given by

B(q)q + h(q, q) = Sτ + J T C F C + J T F x ; (9)
compared to (2), the joint torques vector is pre-multiplied by a matrix S ∈ R (n+6)×n that maps actuation torques into torques for the full floating-base robot:

S = 0 6×n I n×n . (10)
Finally, contact forces are taken into account by introducing the Jacobian of all support links J C ∈ R k×(n+6) and the corresponding overall contact wrench F C ∈ R k , with k equal to the contact constraint dimension (e.g. k = 12 for a humanoid in double support). By repeating the same steps as for the fixed-base case, the task dynamics is obtained as follows:

Λ(q)ẍ + µ(q, q) = F τ |x + F C|x + F x . (11)
It can be seen that the contribution from actuated joint torques is changed slightly w.r.t. the fixed-base case, and it has the following expression:

F τ |x = J T S τ ; (12)
however, this does not affect the mathematical formulation of our prioritized controller. The most important change is given by the coupling between the task dynamics and the contact wrenches, through the term

F C|x = J T J C F C = Λ J B -1 J T C F C ; (13)
on this regard, it is worth noticing that coupling matrix between square brackets in ( 13) is always non-zero regardless of the specific task x: indeed, all tasks are coupled through the virtual chain connecting the floating-base to the world frame. Such a coupling hinders the direct application of our algorithm [START_REF] Hoffman | Multi-priority cartesian impedance control based on quadratic programming optimization[END_REF]; in order to tackle the problem, the authors see three possible methods, as explained

• pre-optimization of contact forces. This approach has been successfully used in [START_REF] Ott | Posture and balance control for biped robots based on contact force optimization[END_REF] for balancing control applied to the lower-body of the TORO robot. From the robot centroidal dynamics, the authors compute contact forces that achieve the desired center-of-mass and angular momentum behavior. More specifically, the angular momentum is used to perform orientation control of the robot base link.

Once that such forces have been obtained, they can be made to disappear from the dynamics formulation by, for instance, redefining the bias torque vector as

ĥ = h -J T C F C . (14)
However, it is the authors' belief that such an idea contains a pitfall, namely that the centroidal dynamics of the robot is constrained to the value obtained during the pre-optimization phase (see remark in Section III-A.2). Consequently, the center-of-mass behaviour is always a first-priority task for the resulting controller, and the same applies to the angular momentum. This is undesirable for a prioritized controller; moreover, it is not clear which reference should be assigned to the angular momentum, given its non-holonomy as stated in [START_REF] Wieber | Holonomy and nonholonomy in the dynamics of articulated motion[END_REF]. • Joint optimization of joint torque and contact forces.

This approach is conceptually similar to the one of [START_REF] Herzog | Momentum control with hierarchical inverse dynamics on a torque-controlled humanoid[END_REF], where the optimization is carried out over contact forces and joint accelerations. Full control over the task hierarchy and system momentum is retained at the cost of an increased number of optimization variables. • Post-optimization of contact forces. A third option, which is the main contribution of the present work, consists in treating the floating-base robot as a fixedbase one, by mathematically replacing the sum of under-actuated joint torques and contact torques with an equivalent completely-actuated torque vector. Then, a post-optimization phase is set up in order to map back the obtained virtual joint torques and forces to equivalent contact wrenches. At the best of the authors' knowledge, this approach has not been explored before.

A full description of such method is the subject of the following section.

C. Post-optimization of contact forces

Starting from ( 9), let us define an equivalent fully-actuated torque vector

τ τ = Sτ + J T C F C ; (15)
with such a definition, the floating-base dynamics formally resembles the fixed-base one:

B(q)q + h(q, q) = τ + J T F x ; (16)
consequently, our algorithm of [START_REF] Hoffman | Multi-priority cartesian impedance control based on quadratic programming optimization[END_REF] can be applied without any modification, yielding some optimized value τ * for the fully-actuated torque vector (15) that permits to achieve the desired hierarchical motion and interaction. Our problem is then to recover the contact force information, by taking into account the under-actuated nature of the system. Essentially, this amounts to solving the system of equations ( 15) for τ and F C . On this regard, notice that the number of equations is N eq = n + 6, while the number of variables is N var = n + k; this means that, in general, the problem admits ∞ Nr many solutions, with N r = k -6.

We can reduce the number of unknowns by exploiting the structure of the actuation matrix [START_REF] Khatib | Whole-body dynamic behavior and control of human-like robots[END_REF]. By focusing on the first six rows of (15) we obtain

τu = J T C,u F C , (17)
where the subscript "u" indicates the sub-matrix corresponding to the unactuated virtual joints. Equation ( 17) is torque-independent, and contains only the contact forces as variables. Once that these have been determined, joint torques can be recovered by looking at the bottom n rows of (15) and solving for τ :

τ = τa -J T C,a F C . ( 18
)
The transpose of the unactuated part of the contact Jacobian acts as a grasp matrix G ∈ R 6×k :

G = J T C,u ; (19)
hence, we can draw inspiration from [START_REF] Ott | Posture and balance control for biped robots based on contact force optimization[END_REF] and obtain a leastsquares solution for F C :

F * C = G † τu , (20)
where the dagger symbol † denotes the pseudo-inverse. However, we can also exploit a QP formulation in order to enforce inequality constraints, as for instance friction cones, as follows:

min FC GF C -τu 2 s.t. b l ≤ D F C ≤ b u u l ≤ F C ≤ u u . (21)
Torque constraints can be introduced as well by considering the dependency on the contact forces as given by [START_REF] Herzog | Balancing experiments on a torque-controlled humanoid with hierarchical inverse dynamics[END_REF].

D. Discussion

The main advantage of the proposed post-optimization formulation with respect to the pre-optimization is clearly given by the proper handling of priorities between tasks, since neither the Center of Mass (CoM) motion nor the robot angular momentum need to be set a-priori. Moreover, the user is relieved from specifying a target angular momentum to the robot, which seems to be problematic as discussed in Section II-B. With our approach, target values for the centroidal dynamics are obtained from the fixed-base solution τ * , and then the corresponding contact wrenches are optimized.

The main drawback is that the first optimization stage may give back a solution that is not feasible under force constraints (e.g. friction cones), which means that the objective value of ( 21) will be greater than zero. This could lead to loss of performance and, in the worst case, instability of the closed loop system. However, this situation can eventually be used to detect the need to perform a recovery action, e.g. a step.

On the other hand, a joint torque-force optimization strategy could perform better in such a case, since constraints on forces are taken into account for motion/force control as well. The price to pay is an increased number of decision variables.

As a final consideration, our post-optimization approach gives us the possibility of adapting our fixed-base formulation to the floating-base case without any need for modification, since the under-actuation is dealt with at a separate stage that is completely decoupled from force/motion control.

III. SIMULATIONS AND PRELIMINARY EXPERIMENTS

In order to validate the proposed extension of our controller [START_REF] Hoffman | Multi-priority cartesian impedance control based on quadratic programming optimization[END_REF] to the floating-base case, we set up two simulation scenarios in Gazebo involving two different legged platforms: CENTAURO [START_REF] Baccelliere | Development of a human size and strength compliant bi-manual platform for realistic heavy manipulation tasks[END_REF], a 39 degrees-of-freedom (DoF) hybrid wheeled-legged quadruped equipped with a humanoid upper-body and COMAN+, which is a 28 DoF, 1.70 meters tall humanoid robot. Both robots are fully torquecontrolled by feeding back the measured link-side joint torques. Moreover, our control architecture XBotCore [START_REF] Muratore | Xbotcore: A real-time cross-robot software platform[END_REF] allows for transparent switch between simulation and real hardware, and hard real-time control on the actual robot. We conclude the section with a brief overview of some preliminary experimental results on the COMAN+ robot.

A. Gazebo simulations

The controller's first stage is formalized in the following way:

         i World TFoot i / World TWaist/ World TLHand + World TRHand / T Posture          << CJoint Torque Limits , (22)
where the symbol A T B denotes a Cartesian impedance task of the frame B relative to the frame A; such a task is implemented by commanding simple virtual wrenches as in the following expression:

F d = K d (x d -x) + D d (x d -ẋ) , (23)
where x, x d , ẋ and ẋd represent the actual and desired Cartesian poses and twists, respectively, whereas K d and D d denote the desired Cartesian impedance. These are combined by means of the operators "+" and "/", which are used to set aggregation and null-space relations, respectively. Finally, the symbol "<<" denotes insertion of constraints into the problem.

Recall from Section II-C that the first stage of our method computes a fully actuated torque vector, that is then mapped to an under-actuated torque vector via post-optimization of contact forces. We implement the post-optimization stage as in [START_REF] Wieber | Modeling and control of legged robots[END_REF], considering linearized friction cones as inequality constraints:

|F t | ≤ √ 2µ 2 F n , F n ≥ 0; (24)
where F t and F n are the tangential and normal components of contact forces, respectively, and µ = 0.3 is the considered friction coefficient. As a final observation, for all our simulation experiments the state of the floating-base link is directly taken from the simulator. 1) CENTAURO balancing under external disturbances: The first simulation scenario consists on a balancing task for the CENTAURO robot under the disturbance of external forces applied on the robot waist. Screen shots from the performed simulations are reported in Fig. 1. The desired Cartesian impedance for the waist position task has been set equal to: K d = 500 500 500 N m and D d = 200 200 200 N s m . A constant force of 200 N is first applied downward along the z-direction for 2 s. The obtained contact forces along the z-direction are shown in Fig. 2(a), while the corresponding waist position error is shown in Fig. 2(b). As expected, the maximum value of the actual waist position error (blue solid line), approximately 0.4 m, is consistent with the expected value (black solid line) obtained through [START_REF] Baccelliere | Development of a human size and strength compliant bi-manual platform for realistic heavy manipulation tasks[END_REF].

To further highlight the impact of friction cones, a second simulation involves the application of a constant force of 90 N sideways along the y-direction, see Fig. 3. Note that, according to [START_REF] Muratore | Xbotcore: A real-time cross-robot software platform[END_REF], the x and y components of the contact 2) COMAN+ picking a box: In the second simulation scenario we consider the COMAN+ robot picking a two kilogram box from the knee level. However, the box weight is not known to the controller. The task has to be performed using whole-body motions to reach the box and pick it up.

Notice that, in order to be able to reach the box, the z-translation task is removed from the waist control. As a side-note, the orientation is managed using a quaternion formulation as in [START_REF] Nakanishi | Operational space control: A theoretical and empirical comparison[END_REF]. The outcome of the simulation experiment can be seen in the accompanying video2 . Remark: Notice that the obtained whole-body motion would be impossible under a pre-optimization framework, in which the CoM trajectory is decided beforehand. In this case, the only solution is to explicitly decrease the CoM height in order for the hands to reach the box.

B. Preliminary Experiments on the COMAN+ platform

In this preliminary experiment we test the presented controller, with the real humanoid robot COMAN+. As the only difference compared to our simulations, in this case we need to estimate the state of our base-link in terms of pose and twist w.r.t. an inertial world frame. In order to do so, we simply fuse kinematic information from the lower-body with IMU measurements, assuming fixed contacts.

After some controller tuning, we managed to have the robot balance on two feet, even when subject to small external perturbations (see Figure 4b). Despite the encouraging results, we experienced instabilities and strong vibrations, in a very similar way as described in [START_REF] Englsberger | Torque-based dynamic walking -A long way from simulation to experiment[END_REF]; indeed, implementation of full-torque controllers on legged robots appears to be challenging, and it needs further investigation from the authors in order to improve the robustness of the controller and successfully deploy it to real systems.

IV. CONCLUSION AND FUTURE WORKS

In this work, a novel method to to formulate task-space inverse dynamics of floating-base robots has been developed as a follow up of the work in [START_REF] Hoffman | Multi-priority cartesian impedance control based on quadratic programming optimization[END_REF]. The method consists in a two step optimization: during a first optimization stage, a fully-actuated torque vector that realizes the desired tasks is found using [START_REF] Hoffman | Multi-priority cartesian impedance control based on quadratic programming optimization[END_REF]; such a vector contains non-zero torques at the virtual joints, and is not directly applicable to the robot. Hence, during a second optimization stage this wrench is mapped back to the available contacts. The proposed method, in comparison to [START_REF] Ott | Posture and balance control for biped robots based on contact force optimization[END_REF], has the main advantage of avoiding to constrain the robot centroidal dynamics.

The method has been tested first in simulation using two floating base robots with substantially different kinematic structure: the CENTAURO and the COMAN+ robot. The first is a wheeled, quadruped robot with humanoid torso, the second a bipedal humanoid robot. Furthermore, preliminary results are reported in the real COMAN+ robot which is subject to small perturbations. Future works will address first the implementation on the real hardware which, as stated previously, at the moment presents some stability problems, very similar to the one in [START_REF] Englsberger | Torque-based dynamic walking -A long way from simulation to experiment[END_REF]. Second we would like to explore a single-stage torque-force optimization, and compare it to [START_REF] Herzog | Momentum control with hierarchical inverse dynamics on a torque-controlled humanoid[END_REF]. Finally, it is the authors' belief that this method can be extended to other problem that can be treated by introducing virtual kinematic chains, for example when contact forces to pick a box has to be computed.

Fig. 1 .

 1 Fig.1. Screen shots from CENTAURO simulations in Gazebo. In the upper plots a constant force of 200 N is applied downward on the robot waist, while in the lower plots a constant force of 90 N is applied sideways.

  (a) z-component of contact forces. (b) Actual (blue solid line) vs. expected (black dashed line) waist position error along the pushing direction, i.e. the z-direction.

Fig. 2 .

 2 Fig. 2. Time histories from CENTAURO simulation: an external constant force of 200 N is applied downward (z-direction) on the robot waist for 2 s (shaded area).

Fig. 3 .

 3 Fig. 3.Contact forces' time history from CENTAURO simulation: an external constant force of 90 N is applied sideways (y-direction) on the robot waist for 2 s (shaded area).

Fig. 4 .

 4 Fig. 4. On the left, COMAN+ in simulation. During the squat motion needed to pick up the box, the knees of the robot (red dot) enter in contact with the table generating unwanted forces which perturb the CoM. Green dots represent intentional contacts. On the right preliminary experiments with the real COMAN+ platform

DIBRIS, Universitá di Genova, Italy * These authors have contributed equally to this work.

We use the subscripts "u" for unactuated, and "a" for actuated, respectively.

The attached video is also available at https://youtu.be/ p8fwwV_zZa8

ACKNOWLEDGMENT

The research leading to these results has received funding from the European Unions Horizon 2020 research and innovation program under grant agreements No. 644839 (CENTAURO) and No. 644727 (CogIMon).