Riccardo Testolin
email: riccardo.testolin@gmail.com

Carlos Pérez Carmona

Fabio Attore

Peter Borchardt

Helge Bruelheide

Jiri Dolezal

Manfred Finckh

Sylvia Haider

Andreas Hemp

Ute Jandt

Andrei Yu Korolyuk

Jonathan Lenoir

Natalia Makunina

George P Malanson

Ladislav Mucina

Jalil Noroozi

Arkadiusz Nowak

Robert K Peet

Gwendolyn Peyre

Francesco Maria Sabatini

Jozef Šibík

Petr Sklenář

Kiril Vassilev

Risto Virtanen

Susan K Wiser

Evgeny G Zibzeev

Borja Jiménez-Alfaro

Global functional variation in alpine vegetation

Keywords: Alpine biomes, alpine vegetation, evolutionary history, functional convergence, macroclimate, phylogenetic dissimilarity, trait pools, trait probability density

Questions: What are the functional trade-offs of vascular plant species in global alpine ecosystems? How is functional variation related to vegetation zones, climatic groups and biogeographic realms? What is the relative contribution of macroclimate and evolutionary history in shaping the functional variation of alpine plant communities? Location: Global.

Methods:

We compiled a dataset of alpine vegetation with 5,532 geo-referenced plots, 1,933 species and six plant functional traits. We used principal component analysis to quantify functional trade-offs among species and trait probability density to assess the functional dissimilarity of alpine vegetation in different vegetation zones, climatic groups and biogeographic realms. We used multiple regression on distance matrices to model community functional dissimilarity against environmental and phylogenetic dissimilarity, controlling for geographic distance.

Results:

The first two PCA axes explained 66% of the species' functional variation and were related to the leaf and stem economic spectra, respectively. Trait probability density was largely independent of vegetation zone and macroclimate but differed across biogeographic realms. The same pattern emerged for both species pool and community levels. The effects of environmental and phylogenetic dissimilarities on community functional dissimilarity had similar magnitude, while the effect of geographic distance was negligible.

Conclusions: Plant species in alpine areas reflect the global variation of plant function, but with a predominant role of resource-use strategies. Current macroclimate exerts a limited effect on alpine vegetation, mostly acting at the community level in combination with evolutionary 6 history. Global alpine vegetation is functionally unrelated to the vegetation zones in which it is embedded, exhibiting strong functional convergence across regions.

A. was funded by the Marie Curie Clarín-COFUND program of the Principality of Asturias-EU (ACB17-26) and the Spanish Research Agency (AEI/ 10.13039/501100011033)". sPlot was funded by the German Research Foundation as one of iDiv's (DFG FZT 118, 202548816) research platforms.

Introduction

Alpine environments (i.e. high-elevation habitats above the climatic treeline) cover about 3% of land outside Antarctica [START_REF] Körner | A definition of mountains and their bioclimatic belts for global comparisons of biodiversity data[END_REF][START_REF] Testolin | Global distribution and bioclimatic characterization of alpine biomes[END_REF] and can be found in all continents and at all latitudes [START_REF] Körner | Alpine plant life. Functional plant ecology of high mountain ecosystems[END_REF]. These habitats include global biodiversity hotspots [START_REF] Myers | Biodiversity hotspots for conservation priorities[END_REF] and support about 10,000 plant species worldwide, many of which are endemics [START_REF] Körner | Alpine plant life. Functional plant ecology of high mountain ecosystems[END_REF]. Globally, the vegetation of alpine environments is dominated by few growth forms (e.g. dwarf shrubs, graminoids, herbaceous rosettes and cushions), reflecting functional adaptations to the characteristics of high-mountain ecosystems, such as low temperatures, short growing season and limited nutrient availability [START_REF] Körner | Alpine plant diversity: a global survey and functional interpretations[END_REF][START_REF] Körner | Alpine plant life. Functional plant ecology of high mountain ecosystems[END_REF][START_REF] Körner | Plant adaptations to alpine environments[END_REF]Dolezal et al., 2016;[START_REF] Stanisci | Functional composition and diversity of leaf traits in subalpine versus alpine vegetation in the Apennines[END_REF]. Some growth forms, however, are more abundant in certain regions (e.g. evergreen dwarf shrubs in boreal ranges, succulents in semi-arid zones, sclerophyllous species in mediterranean-type climates) or are unique to specific areas, such as giant rosettes in tropical mountains (e.g. Espletia and Dendrosenecio) (Nagy and Grabherr, 2009). Nevertheless, growth forms are poor descriptors of the functional adaptations of alpine vegetation, with several features of alpine plants found to vary widely within a single growth form [START_REF] Körner | Alpine plant diversity: a global survey and functional interpretations[END_REF] or show no variation among different growth forms [START_REF] Körner | Carbon and nitrogen stable isotope signals for an entire alpine flora, based on herbarium samples[END_REF].

In alpine environments, plants have adapted to low temperatures and low nutrient supply (Nagy and Grabherr, 2009). In comparison to lowland species, alpine plants are normally shorter, with smaller leaves and lighter seeds [START_REF] Körner | Alpine plant life. Functional plant ecology of high mountain ecosystems[END_REF][START_REF] Pellissier | Plant traits co-vary with altitude in grasslands and forests in the European Alps[END_REF]. These traits increase frost tolerance, photosynthetic efficiency and dispersal success, so they are globally ubiquitous in alpine vegetation [START_REF] Körner | Alpine plant life. Functional plant ecology of high mountain ecosystems[END_REF]. Yet, considerable variation remains among species from different alpine regions, e.g., differences in leaf traits [START_REF] Halloy | Comparative leaf morphology spectra of plant communities in New Zealand, the Andes and the European Alps[END_REF][START_REF] Pyankov | Leaf structure and specific leaf mass: The alpine desert plants of the Eastern Pamirs, Tadjikistan[END_REF], suggesting a response to environmental and evolutionary drivers. Most research comparing plants across alpine regions, however, has only focused on individual traits.

Assessing how multiple traits vary simultaneously may allow identification of the trade-offs of plant form and function, i.e., the different strategies used by alpine plants for resource acquisition, growth and reproduction [START_REF] Grime | Vegetation classification by reference to strategies[END_REF][START_REF] Díaz | The global spectrum of plant form and function[END_REF][START_REF] Bruelheide | Global trait-environment relationships of plant communities[END_REF].

Global alpine areas can be grouped according to their macroclimate [START_REF] Testolin | Global distribution and bioclimatic characterization of alpine biomes[END_REF], and they are linked to different vegetation zones [START_REF] Walter | Global classification of natural terrestrial ecosystems[END_REF] and biomes characterised by distinct evolutionary history and own species pools [START_REF] Mucina | Biome: evolution of a crucial ecological and biogeographical concept[END_REF]. Indeed, present-day alpine floras are the result of upward shifts of species undergoing regional radiations and long-distance migrations associated with the displacement of cold-climate biomes, such as during the Pleistocene glacial cycles [START_REF] Billings | Adaptations and origins of alpine plants[END_REF][START_REF] Horandl | The evolution and biogeography of alpine species in Ranunculus (Ranunculaceae): a global comparison[END_REF][START_REF] Jiménez-Alfaro | Postglacial determinants of regional species pools in alpine grasslands[END_REF]. The historical legacy of ancestral species, which may have belonged to different vegetation zones and biogeographic realms, together with the environmental filtering of the current macroclimate, have determined the diversity of alpine trait pools, i.e. the total set of plant trait values found in an alpine region today.

Factors selecting for favourable combinations of traits are generally scale-dependent (Garnier et al., 2016). At continental scales, trait pools are defined by the interplay of macroclimate and evolutionary history [START_REF] Moncrieff | Revising the biome concept for understanding and predicting global change impacts[END_REF][START_REF] Mucina | Biome: evolution of a crucial ecological and biogeographical concept[END_REF], with the latter constrained by the long-term isolation of major landforms [START_REF] Chaboureau | Tectonic-driven climate change and the diversification of angiosperms[END_REF] and by the phylogenetic origin of species occurring in a biogeographic realm [START_REF] Holt | An update of Wallace's zoogeographic regions of the world[END_REF][START_REF] Daru | Understanding the processes underpinning patterns of phylogenetic regionalization[END_REF][START_REF] Daru | Unravelling the evolutionary origins of biogeographic assemblages[END_REF]. At the scale of local plant communities, trait pools are further constrained by biotic and abiotic filters that select species assemblages with favourable trait syndromes [START_REF] Lavorel | Plant functional classifications: from general groups to specific groups based on response to disturbance[END_REF][START_REF] Zobel | The species pool concept as a framework for studying patterns of plant diversity[END_REF][START_REF] Mucina | Biome: evolution of a crucial ecological and biogeographical concept[END_REF]. As a consequence, the trait values in communities might deviate from those of trait pools [START_REF] Grime | Trait convergence and trait divergence in herbaceous plant communities: Mechanisms and consequences[END_REF][START_REF] Marks | Alternative designs and the evolution of functional diversity[END_REF] and may depend on local conditions (e.g. soil properties and topoclimate) rather than macroclimate [START_REF] Bruelheide | Global trait-environment relationships of plant communities[END_REF]. However, considering the varied origin of plants across global alpine environments [START_REF] Billings | Adaptations and origins of alpine plants[END_REF], an evolutionary mark on functionality might still be detectable at the level of communities [START_REF] Srivastava | Phylogenetic diversity and the functioning of ecosystems[END_REF]. Linking local filtering to evolutionary and biogeographic history remains a major challenge in macroecology and new approaches that incorporate different facets of diversity are required to understand patterns and processes across scales [START_REF] Pärtel | Macroecology of biodiversity: disentangling local and regional effects[END_REF][START_REF] Ladouceur | The functional trait spectrum of European temperate grasslands[END_REF]. Disentangling the effect of macroclimate and evolutionary history might therefore open new prospects for understanding, and possibly predicting, biodiversity patterns in alpine regions.

Here, we provide the first overview of the functional variation of alpine vegetation and an attempt to infer possible drivers across spatial scales. Specifically, we aim to 1) describe the functional trade-offs of vascular plant species in global alpine ecosystems; 2) assess the functional variation of trait pools and local communities among vegetation zones, climatic groups and biogeographic realms; and 3) quantify the relative contribution of macroclimate and evolutionary history in shaping the functional variation of alpine plant communities.

Methods

Study system and data selection

We used data featuring alpine vegetation defined as any vascular plant community above the climatic treeline [START_REF] Körner | Alpine plant life. Functional plant ecology of high mountain ecosystems[END_REF]. In addition to strictly zonal habitats dominated by graminoids, forbs and dwarf shrubs, we also included snow-patch plant communities and vegetation on rocks and screes, as they are also found ubiquitously across the alpine belt. The plot data collected by the authors, compiled from the literature, or stored in the sPlot database (v2.1) [START_REF] Bruelheide | sPlot -A new tool for global vegetation analyses[END_REF], was first filtered using habitat classifications of the data sources (Table S1 Appendix S1), and then further reduced by excluding plots with tree species or incomplete taxonomic identification. We standardised datasets from different sources by identifying a minimum common set of plot attributes including plot size, elevation and geographic coordinates. Species names were harmonised using the Taxonomic Name Resolution Service [START_REF] Boyle | The taxonomic name resolution service: an online tool for automated standardization of plant names[END_REF] (http://tnrs.iplantcollaborative.org) with default settings. Species cover values coded with discrete scales were transformed to the mean value of the corresponding percentage interval.

Subspecies and varieties were merged at the species resolution by summing the respective percentage cover values. At this point, the dataset consisted of 8,419 plots of alpine vascular vegetation with 4,651 plant species recorded.

Each plot was assigned to the vegetation zone dominating the same ecoregion, i.e. montane grasslands and shrublands, temperate broadleaf and mixed forests, temperate coniferous forests, tropical and subtropical moist broadleaf forests, and tundra [START_REF] Olson | Terrestrial ecoregions of the world: a new map of life on earth[END_REF]. These physiognomic types encompassing large areas are presumed to contribute ancestral clades with potential impact on current alpine trait pools. We also assigned the plots to one of three groups summarizing the climatic variability of global alpine areas, representing regional alpine biomes in the classification scheme of [START_REF] Testolin | Global distribution and bioclimatic characterization of alpine biomes[END_REF]: 1) oceanic, characterised by greater precipitation and relative temperature stability; 2) continental, defined by low precipitation and large annual temperature amplitudes; 3) subtropical, encompassing both tropical and subtropical alpine areas and characterised by low annual precipitation and contrasting diurnal temperature cycles. Single plots falling slightly outside the boundaries of the commonest vegetation zone or climatic group for a given region were manually assigned to those. Finally, each plot was assigned to a biogeographic realm: Afrotropic, Australasia, Nearctic, Neotropic or Palearctic.

Each realm represents a broadly defined geographic region characterised by typical flora and fauna and supposed to have a distinct evolutionary history.

For each species, we extracted the gap-filled trait information from the TRY database (v5.0) [START_REF] Shan | Gap filling in the plant kingdom -Trait prediction using hierarchical probabilistic matrix factorization[END_REF]Fazayeli et al., 2014;[START_REF] Schrodt | BHPMF -a hierarchical Bayesian approach to gap-filling and trait prediction for macroecology and functional biogeography[END_REF][START_REF] Kattge | TRY plant trait databaseenhanced coverage and open access[END_REF], provided by the sPlot database as species average values [START_REF] Bruelheide | sPlot -A new tool for global vegetation analyses[END_REF]. We selected six plant functional traits: leaf area (one-sided surface of the fresh leaf), specific leaf area (leaf area per leaf dry mass; SLA), leaf dry matter content (leaf dry mass per leaf fresh mass; LDMC), leaf nitrogen (N) (N per leaf dry mass), plant height (maximum total height of the plant) and seed mass (dry mass of the seed). We chose these traits because they are commonly used to characterise tundra and alpine vegetation [START_REF] Bjorkman | Tundra Trait Team: A database of plant traits spanning the tundra biome[END_REF][START_REF] Thomas | Traditional plant functional groups explain variation in economic but not size-related traits across the tundra biome[END_REF][START_REF] Liancourt | Plant's-eye view of temperature governs elevational distributions[END_REF], and they are fully representative of plant ecological strategies [START_REF] Díaz | The global spectrum of plant form and function[END_REF].

The gap-filling process employed hierarchical Bayesian modelling to estimate missing trait values based on other traits available in TRY for individuals of the same species [START_REF] Schrodt | BHPMF -a hierarchical Bayesian approach to gap-filling and trait prediction for macroecology and functional biogeography[END_REF]. Only traits of those species having at least one measured trait observation were imputed.

Of all species that were selected based on trait data availability, 99% had at least one measurement in TRY for leaf area, 96% for SLA, 98% for LDMC, 97% for leaf N, 91% for plant height and 85% for seed mass. The values of plant height, leaf area and seed mass were log10transformed to reduce skewness. Species for which trait information was not available (n = 2,517) were removed. At the community level, we only considered plots with at least 50% cumulative cover of species with trait data. We chose 50% cover as a trade-off between the inclusion of plots for which trait data were scarce versus the representativeness of dominant vegetation in each community. An alternative set of results obtained choosing more conservative thresholds of 75% and 90% cumulative cover of species with trait data (Appendix S2) showed minor differences with the results presented here. The final dataset consisted of 5,532 vegetation plots between 0.25 and 400 m 2 in size sampled between 1923 and 2019, with 1,933 species belonging to five vegetation zones, three climatic groups and five biogeographic realms (Figure 1). All the following analyses have been carried out using R 3.6.3 (R Core Team, 2020).

Functional trade-offs and variation of trait pools

To analyse the relationships among traits of the plant species in our dataset, we performed a principal component analysis (PCA) of the standardised values of the six traits. The loadings of the individual traits were then used to identify the main axes of variation and possible trade-offs of plant strategies [START_REF] Díaz | The global spectrum of plant form and function[END_REF].

To compare the trait pools across the vegetation zones, climatic groups and biogeographic realms, we employed trait probability densities, a scale-independent framework that implements the concept of the niche hypervolume while accounting for the probabilistic nature of traits [START_REF] Carmona | Traits without borders: integrating functional diversity across scales[END_REF]. This method requires both the mean and the standard deviation of each trait for all the species. As reliable information on the standard deviations was not available, we assumed it to be constant across species and estimated it as 50% of the standard deviation of all species' mean values for each trait [START_REF] Lamanna | Functional trait space and the latitudinal diversity gradient[END_REF][START_REF] Carmona | TPD: Methods for Measuring Functional Diversity Based on Trait Probability Density[END_REF]. Then, we calculated the individual trait pools as the probability densities for each vegetation zone, climatic group and biogeographic realm using the package "TPD" [START_REF] Carmona | TPD: Methods for Measuring Functional Diversity Based on Trait Probability Density[END_REF], accounting for species frequencies (i.e. the number of plots in each group where a certain species was recorded). We assessed the functional variation of trait pools using kernel density plots and calculated pair-wise functional dissimilarities among trait pools using the "dissim" function of the "TPD" package. The significance of the pair-wise dissimilarities was evaluated in a nullmodelling framework (Geange et al., 2011;[START_REF] Traba | Realised niche changes in a native herbivore assemblage associated with the presence of livestock[END_REF] by pooling the observations from each pair, randomizing the species' labels 999 times while keeping the number of species constant for each group and ranking the pair-wise dissimilarity values among the simulated trait probability densities. This allowed us to calculate the Bonferroni-corrected p-values for each comparison as:

𝑝 = (1 - 𝑟 𝑜𝑏𝑠 𝑖 + 1 ) × 𝑛
where 𝑟 𝑜𝑏𝑠 is the rank of the observed dissimilarity value among the simulated ones, 𝑖 is the number of simulations, and 𝑛 is the number of pair-wise comparisons (Legendre and Legendre, 2012;[START_REF] Traba | Realised niche changes in a native herbivore assemblage associated with the presence of livestock[END_REF]. To assess the overall functional variation among vegetation zones, climatic groups and biogeographic realms, while excluding potentially redundant information, we calculated multi-trait probability densities by using the first two axes of the PCA of the six traits and repeated the same analyses described above for individual traits.

Functional variation of communities

To analyse the variation in trait values across plant communities, we calculated the multi-trait functional dissimilarities between all vegetation plots as described above, accounting for the cover of the species within each plot. The pair-wise dissimilarities were displayed in the same PCA space as the individual species by calculating the community weighted means of the first two PCA axes for each plot. This allowed the visualisation of the functional variation of plots belonging to different groups. Significant differences among plots belonging to different vegetation zones, climatic groups and biogeographic realms were tested using PERMANOVA [START_REF] Anderson | A new method for non-parametric multivariate analysis of variance[END_REF], implemented by the "adonis" function of the R package "vegan" [START_REF] Oksanen | Community Ecology Package. R package version 2[END_REF], with 999 permutations.

To quantify the relative contribution of climate and evolutionary history in determining functional variation among communities, we modelled community functional dissimilarity as a function of environmental and phylogenetic dissimilarity while controlling for geographic distance. Phylogenetic data were provided by the sPlot database based on the phylogeny of [START_REF] Qian | An updated megaphylogeny of plants, a tool for generating plant phylogenies and an analysis of phylogenetic community structure[END_REF]. Species present in our dataset but missing from this phylogeny were added next to a randomly selected congener, if available [START_REF] Bruelheide | sPlot -A new tool for global vegetation analyses[END_REF]. First, we selected the set of species for which both trait and phylogenetic data were available (n = 1,674) and further subset the vegetation plots by keeping those with at least 50% cumulative cover of these species.

Thus, we obtained a subset of 5,047 plots and calculated the multi-trait functional dissimilarities between all possible pairs of plots as described above. We also performed alternative selections of plots with 75% and 90% cumulative cover of species with trait and phylogenetic data to assess the effect of a more conservative cumulative cover threshold on the model results (Appendix S2). Then, we built a set of climatic variables known to affect alpine vegetation [START_REF] Körner | Alpine plant life. Functional plant ecology of high mountain ecosystems[END_REF][START_REF] Moser | Environmental determinants of vascular plant species richness in the Austrian Alps[END_REF]Nagy and Grabherr, 2009) using data from CHELSA bioclimatic database at ~1 km spatial resolution [START_REF] Karger | Climatologies at high resolution for the earth's land surface areas[END_REF]. The included variables were mean temperature, precipitation, growing degree days and mean potential evapotranspiration. Each variable was calculated within the time frame of the growing season, defined as days with mean temperature >0.9 °C [START_REF] Paulsen | A climate-based model to predict potential treeline position around the globe[END_REF]. Growing degree days (i.e. the sum of monthly temperatures > 0.9 °C multiplied by the total number of days) were calculated using the "growingDegDays" function of the R package "envirem" [START_REF] Title | ENVIREM: an expanded set of bioclimatic and topographic variables increases flexibility and improves performance of ecological niche modeling[END_REF]. Mean potential evapotranspiration of the growing season was estimated with the "hargreaves" function of the R package "SPEI" (Beguería and Vicente-Serrano, 2017), using maximum and minimum monthly values of temperature and monthly precipitation. The monthly values of potential evapotranspiration obtained were then averaged across months with mean temperature above 0.9 °C. We standardised the four climatic variables and calculated the Euclidean distance among each pair of plots as a measure of environmental dissimilarity.

To account for the evolutionary history of plant species in different communities, we also calculated the pair-wise phylogenetic dissimilarity between plots (Ives and Helmus, 2010) with the "pcd" function of the R package "picante" [START_REF] Kembel | Picante: R tools for integrating phylogenies and ecology[END_REF]. To account for the spatial aggregation of plots and unmeasured regional effects on the estimated functional dissimilarity, we calculated the pair-wise geographical distances between plots. Finally, we modelled functional community dissimilarity against these three distanced-based predictors using multiple regression on distance matrices (MRM) with the "lm" function. Despite our measure of functional dissimilarity is constrained between 0 and 1, our dataset mainly encompassed intermediate levels of functional turnover (Figure S1 Appendix S1), allowing us to treat it as approximately linear (Ferrier et al., 2007). Further, a linear modelling approach allowed us to calculate the adjusted R 2 of all the sub-models necessary to perform variance partitioning [START_REF] Borcard | Partialling out the spatial component of ecological variation[END_REF][START_REF] Swenson | Functional and Phylogenetic Ecology in R. Use R! New York[END_REF].

Results

The first two PCA axes accounted for 66% of the total trait variation among species. The other axes explained less variation than expected by chance and were not considered further. The first axis (PC1; 35% of variation) was mainly related to variations in LDMC, leaf N and SLA, while the second axis (PC2; 31% of variation) was linked to leaf area, plant height and seed mass (Figure 2, Table S2 Appendix S1).

When focusing at the level of vegetation zones, we observed negligible differences in trait probability density and low functional dissimilarity among trait pools. The only exception was the alpine vegetation related to tropical and subtropical moist broadleaf forests, which exhibited slightly greater plant height values compared to other vegetation zones (Figure 3; Table S3 Appendix S1). Among climatic groups, subtropical alpine areas also exhibited greater plant height values when compared to oceanic and continental ones, with minor variation in the distribution of other traits (Figure 3; Table S3 Appendix S1). However, we observed considerable variability in trait probability density among biogeographic realms. The alpine vegetation of the Australasia and Neotropic realms had lower SLA compared to that of the others and similar values of plant height to the Afrotropic, which were greater than those of the Palearctic and Nearctic. As for leaf area and seed mass, the Neotropics generally showed higher values compared to the Palearctic and Nearctic, which in turn presented higher leaf N and lower LDMC than Australasia (Figure 3; Table S3 Appendix S1). Multi-trait patterns seemingly reflected those observed at the single trait level. Among vegetation zones and climatic groups, multi-trait functional dissimilarities were not significant or very modest (Table 1). Conversely, Finally, the MRM model fit on a subset of plots with available phylogenetic information explained 16.6% of the communities' functional dissimilarity. Environmental and phylogenetic dissimilarities both explained 6.2% individually, while 4% was shared between the two of them.

Geographic distance exhibited a marginal effect, explaining only 0.3% (Figure 5). Again, adopting more conservative thresholds of cumulative cover of species with trait and phylogenetic data did not significantly affect the results (Appendix 2).

Discussion

Functional trade-offs of alpine plant species

We selected six traits linked to resource use, growth and reproduction of plants and used PCA to describe the functional trade-offs of 1,933 vascular plant species in global alpine ecosystems.

PC1 differentiated strategies in terms of investments of nutrients and dry mass in leaves and hence the leaf economics spectrum [START_REF] Wright | The worldwide leaf economics spectrum[END_REF]. This spectrum discriminates between species, those with high leaf construction costs (high LDMC, low SLA) and low leaf nutrient concentrations (low leaf N) related to slower vegetative development rates versus fast-growing species with high leaf nutrient concentration and cheaper construction costs that promote a quick return of the investments in nutrients and carbon [START_REF] Wright | The worldwide leaf economics spectrum[END_REF]. PC2 reflected differences in plant size, conforming to the stem economics spectrum [START_REF] Baraloto | Decoupled leaf and stem economics in rain forest trees[END_REF] that separates taller plants able to carry larger leaves and seeds (large plant height, leaf area and seed mass) from smaller plants. These results agree with previous analyses of alpine and tundra vascular plants (Dolezal et al., 2016;[START_REF] Thomas | Traditional plant functional groups explain variation in economic but not size-related traits across the tundra biome[END_REF] and are consistent with directions of variation in the global spectrum of plant form and function [START_REF] Díaz | The global spectrum of plant form and function[END_REF]. Nevertheless, the predominance of variation in resource-use strategies rather than size reflects the absence of trees and tall shrubs in alpine vegetation and the general abundance of prostrate species which are mainly differentiated by local conditions. Small size allows alpine plants to respond to and modify the microclimate near the ground [START_REF] Geiger | The climate near the ground[END_REF] by accumulating heat under the leaf canopy regardless of fluctuations of the macroenvironment [START_REF] Körner | Functional morphology of mountain plants[END_REF][START_REF] Körner | Alpine plant life. Functional plant ecology of high mountain ecosystems[END_REF]. Additionally, the main variation observed along PC1 could be explained by the greater variability of leaf construction costs of alpine plants, which depend on local temperature, frost stress and prolonged exposure to light [START_REF] Körner | Functional morphology of mountain plants[END_REF]. Any of these may vary widely even within a single mountain range, hence the greater variation of the related traits [START_REF] Stanisci | Functional composition and diversity of leaf traits in subalpine versus alpine vegetation in the Apennines[END_REF].

Variation of alpine trait pools

Trait pools of alpine plants were largely independent of the vegetation zone, suggesting that alpine vegetation is functionally different from the surrounding flora in which it is embedded.

Thus, the convergence of growth forms that characterises alpine vegetation [START_REF] Körner | Alpine plant life. Functional plant ecology of high mountain ecosystems[END_REF][START_REF] Körner | Plant adaptations to alpine environments[END_REF][START_REF] Aubert | A revised worldwide catalogue of cushion plants 100 years after Hauri and Schröter[END_REF] follows adaptation to similar ecological conditions [START_REF] Givnish | Giant lobelias exemplify convergent evolution[END_REF][START_REF] Givnish | Convergent evolution, adaptive radiation, and species diversification in plants[END_REF][START_REF] Horandl | The evolution and biogeography of alpine species in Ranunculus (Ranunculaceae): a global comparison[END_REF][START_REF] Hughes | The ubiquity of alpine plant radiations: From the Andes to the Hengduan Mountains[END_REF]. This finding contrasts with the view of alpine areas as elevational "orobiomes" closely related to the "zonobiomes" they originate from [START_REF] Walter | Global classification of natural terrestrial ecosystems[END_REF] but agrees with the distinction of alpine ecosystems from other terrestrial biomes [START_REF] Testolin | Global distribution and bioclimatic characterization of alpine biomes[END_REF]. Trait pools were also convergent among climatic groups, indicating that macroclimatic differences above the treeline have little influence on the functional features of alpine vegetation, which is consistent with the similar patterns of primary productivity found across global alpine biomes [START_REF] Testolin | Global distribution and bioclimatic characterization of alpine biomes[END_REF].

However, we observed some divergence of trait pools across biogeographic realms. Such functional differences might emerge even among structurally similar plant groups when these are compared across areas with distinct evolutionary histories [START_REF] Alvarado-Cárdenas | To converge or not to converge in environmental space: Testing for similar environments between analogous succulent plants of North America and Africa[END_REF].

Specifically, we observed a distinction between the trait pool of the Holarctic realm and those of the Neotropic and Australasia, while the trait pool of the Afrotropic realm occupied an intermediate position. This pattern likely reflects different evolutionary histories and adaptations of alpine vegetation in the Northern and Southern Hemispheres [START_REF] Billings | Adaptations and origins of alpine plants[END_REF]. Indeed, much of the ancestral alpine vascular flora originated during the Miocene (23-5 Ma) from Arcto-Tertiary and Antarcto-Tertiary floras through upward migration and evolution of lowland taxa [START_REF] Billings | Adaptations and origins of alpine plants[END_REF]. Consequently, Holarctic alpine vegetation shares many species with the Arctic [START_REF] Billings | Adaptations and origins of alpine plants[END_REF] and has major links with Afrotropic alpine species [START_REF] Linder | The evolution of African plant diversity[END_REF][START_REF] Carbutt | Reconciling ecological and phytogeographical spatial boundaries to clarify the limits of the montane and alpine regions of sub-Sahelian Africa[END_REF]. In contrast, a large part of Neotropic alpine plants originated locally through migration and adaptation of Neotropical lowland species [START_REF] Sklenář | Tropical and temperate: Evolutionary history of páramo flora[END_REF], some of which also contributed to Afrotropic lineages [START_REF] Linder | The evolution of African plant diversity[END_REF]. Finally, the functional similarity of Neotropics and Australasia probably derived from both migration [START_REF] Raven | Plate tectonics and Australasian paleobiogeography[END_REF] and convergent evolution during the Pliocene (5 Ma) and the Pleistocene (2.5 Ma), when further mountain uplift and repeated glaciations led to the diversification of the respective alpine floras [START_REF] Mcglone | Endemism, species selection and the origin and distribution of the vascular plant flora of New Zealand[END_REF][START_REF] Winkworth | Evolution of the New Zealand mountain flora: Origins, diversification and dispersal[END_REF][START_REF] Sklenář | Tropical and temperate: Evolutionary history of páramo flora[END_REF][START_REF] Madriñán | Páramo is the world's fastest evolving and coolest biodiversity hotspot[END_REF].

Functional variation of alpine communities

Communities were not functionally distinct among vegetation zones or climatic groups, while biogeographic realms exhibited greater discriminatory power, as they did for the trait pools.

Australasian communities form an isolated group characterised by leaves with high construction costs, which agrees with the greater abundance of sclerophyllous dwarf shrubs in the Australasian alpine flora relative to other global alpine regions [START_REF] Ballantyne | Shrub facilitation is an important driver of alpine plant community diversity and functional composition[END_REF].

The functional distinctness of Australasian alpine communities from Holarctic and Tropical ones may also reflect differences in trait pools between the two hemispheres and could be related to the long-time isolation (45 -49 Ma) of Australasian landforms from other biogeographic realms of Gondwanan origin [START_REF] Raven | Plate tectonics and Australasian paleobiogeography[END_REF]. Holarctic and Tropical communities, however,

were not as functionally distinct as their trait pools, indicating that other processes apart from regional evolutionary history are involved at the local scale. Indeed, although our model highlighted the presence of a phylogenetic signal in functional dissimilarity, environmental dissimilarity explained an equal amount of variance. This is consistent with the process of niche conservatism in highly heterogeneous areas, where the retention of the ancestral niche characteristics could lead to both conservatism and divergence of the realised niche (i.e. the functional characteristics; [START_REF] Pyron | Phylogenetic niche conservatism and the evolutionary basis of ecological speciation[END_REF]. Nevertheless, the negligible effect of geographic distance and the large amount of unexplained variance point to fine-scale environmental factors (e.g. soil properties and topoclimate), disturbance and biotic interactions as the main drivers of community trait composition in alpine ecosystems [START_REF] Grime | Trait convergence and trait divergence in herbaceous plant communities: Mechanisms and consequences[END_REF][START_REF] Dolezal | Functionally distinct assembly of vascular plants colonizing alpine cushions suggests their vulnerability to climate change[END_REF].

Assumptions and caveats

Even though we used the largest dataset of alpine vegetation ever collected, our study does not come without uncertainties. First, many mountain regions, including outstanding centres of alpine plant diversity such as the Himalayas and Hengduan Mountains (Favre et al., 2015;[START_REF] Xing | Uplift-driven diversification in the Hengduan Mountains, a temperate biodiversity hotspot[END_REF][START_REF] Muellner-Riehl | Origins of global mountain plant biodiversity: testing the 'mountain-geobiodiversity hypothesis[END_REF][START_REF] Ding | Ancient orogenic and monsoondriven assembly of the world's richest temperate alpine flora[END_REF], were not represented in our dataset, preventing us to provide a complete global picture of alpine plant functional variation.

Still, our plots encompass alpine vegetation in six continents from boreal to tropical latitudes, allowing meaningful global comparisons that could be further refined by the future inclusion of additional alpine regions, especially in the tropical and subtropical belts. Second, when comparing functional dissimilarities across geographical units and spatial scales, we presumed that the species for which trait data were available were also representative of the dominant vegetation in our study areas. For several tropical species, however, such data were not available, and we had to exclude many plots in Africa and South America. Although we recognise that this could have led to the exclusion of unusual combinations of traits and that even rare species can drive trait divergence among communities at the regional scale [START_REF] Richardson | Rare species drive local trait diversity in two geographically disjunct examples of a naturally rare alpine ecosystem in New Zealand[END_REF], this is probably less relevant at the global level. Third, we note that our dataset encompasses vegetation plots of very different sizes (0.25 -400 m 2 ). As species richness generally increases with area [START_REF] Lomolino | Ecology's most general, yet protean pattern: the species-area relationship[END_REF], larger plots might be functionally richer than smaller ones [START_REF] Smith | Characterizing scale-dependent community assembly using the functional-diversity-area relationship[END_REF][START_REF] Wang | Phylogenetic and functional diversity area relationships in two temperate forests[END_REF], biasing the comparison among plots. However, when accounting for species abundancesor, in our case, coverthe relationship between functional diversity indices and plot size tends to weaken or disappear because of species dominance and functional redundancy [START_REF] Karadimou | Functional diversity exhibits a diverse relationship with area, even a decreasing one[END_REF]. Therefore, as the trait probability density framework accounts for the distribution of trait values in plant communities, plot size likely had a minor effect in the estimation of functional dissimilarity among alpine communities. Although we collected most of the plot data in alpine vegetation currently available, much effort is still needed to collect data with a consistent sampling protocol, including functional traits and a proper representation of species and vegetation types from disparate global regions.

Conclusions

This study provides the first overview of the global functional trait variation in alpine vegetation.

While alpine species exhibit the same trade-offs observed in vascular plants globally, the absence of trees in alpine ecosystems leads to a greater variety of traits related to resource-use strategies rather than size. We found that alpine vegetation is scarcely related to the vegetation zones in which it is embedded and is largely independent of macroclimatic patterns, at least for the traits analysed in this study. However, evolutionary history seemingly affected current trait pools, and phylogenetic constraints and macroclimate equally determine the functional dissimilarity of communities. Overall, our results indicate a strong functional convergence of adult plant traits in global alpine vegetation, with implications at both regional and community level. This finding agrees with the functional convergence observed for regeneration traits in alpine plants across continents (Fernández-Pascual et al., 2020), further supporting a distinct delineation of alpine ecosystems in the context of the global biomes. Yet, other factors not accounted for in this study (e.g. soil properties, topoclimatic gradients) are likely influencing functional traits of alpine vegetation locally. In this respect, future work should be oriented toward the inclusion of additional fine-scale environmental characteristics, as well as trait data from tropical and subtropical species currently underrepresented in global datasets.

Tables Table 1. Multi-trait pair-wise dissimilarities (Diss) of alpine vegetation between vegetation zones, climatic groups and biogeographic realms. Significant dissimilarities (P < 0.05) are in bold.

Diss

  among biogeographic realms, Palearctic and Nearctic were similar to one another and differentiated from Neotropic and Australasia, with the Afrotropic pool taking an intermediate position between the two groups.Multi-trait dissimilarities of alpine plant communities revealed distinct patterns among biogeographic realms (Figure4), with Australasian, Afrotropic and Neotropic plots characterised by larger values of LDMC and smaller SLA and leaf N. PERMANOVA showed that biogeographic realms explained 19% of the functional variation (R 2 = 0.19, p = 0.001), vegetation zones explained 11% (R 2 = 0.11, p = 0.001) and climatic groups explained 5% (R 2 = 0.05, p = 0.001). The same patterns emerged when considering more conservative thresholds of cumulative cover of species with trait data (Appendix S2).

FiguresFigure 1 .

 1 Figures

  Figure 1. Spatial distribution of 5,532 alpine vegetation plots across vegetation zones, climatic groups and biogeographic realms.

Figure 2 .

 2 Figure 2. Functional variation of 1,933 vascular plant species in global alpine areas along the first two principal components of six traits representing main functional trade-offs. SLA = Specific leaf area; LDMC = Leaf dry matter content.

Figure 3 .Figure 4 .

 34 Figure 3. Kernel density plots of alpine trait pools estimated using trait probability density for six individual plant traits among vegetation zones, climatic groups and biogeographic realms. SLA = Specific leaf area; LDMC = Leaf dry matter content.

  Vegetation zones Montane grasslands and shrublands -Temperate broadleaf and mixed forests 0.16 ns Montane grasslands and shrublands -Temperate coniferous forests 0.27 ** Montane grasslands and shrublands -Tropical and subtropical moist broadleaf forests 0.29 ns Montane grasslands and shrublands -Tundra 0.27 ns Temperate broadleaf and mixed forests -Temperate coniferous forests 0.24 ** Temperate broadleaf and mixed forests -Tropical and subtropical moist broadleaf forests 0.34 * Temperate broadleaf and mixed forests -Tundra 0.23 ns Temperate coniferous forests -Tropical and subtropical moist broadleaf forests 0.28 ns Temperate coniferous forests -Tundra 0.25 ns Tropical and subtropical moist broadleaf forests -Tundra 0.38 ns Significance codes: ***: p<0.001; **: p< 0.01; *: p<0.05; ns: p ≥ 0.05.

	Climatic groups
	Continental -Oceanic	0.19 **
	Continental -Subtropical	0.22 *
	Oceanic -Subtropical	0.31 **
	Biogeographic realms
	Afrotropic -Australasia	0.40 ns
	Afrotropic -Nearctic	0.26 ns
	Afrotropic -Neotropic	0.29 ns
	Afrotropic -Palearctic	0.23 ns
	Australasia -Nearctic	0.47 **
	Australasia -Neotropic	0.36 ns
	Australasia -Palearctic	0.45 **
	Nearctic -Neotropic	0.43 **
	Nearctic -Palearctic	0.13 ns
	Neotropic -Palearctic	0.41 **

Acknowledgments: The authors would like to thank Keith McDougall (Department of Planning, Industry and Environment, Queanbeyan NSW, Australia) for providing the vegetation data for the Australian Alps, Wolfgang Willner (Department of Botany and Biodiversity Research, University of Vienna, Vienna, Austria) for providing data for the Austrian Alps, and Miska Luoto (Department of Geosciences and Geography, University of Helsinki, Helsinki, Finland) for sampling part of the vegetation data in the Scandinavian Mountains, along with all the other data collectors. We are also grateful to sPlot, the global vegetation-plot database, which was funded by the German Research Foundation as one of iDiv (DFG FZT 118, 202548816) research platforms. This study has been supported by the TRY initiative on plant traits (http://www.trydb.org). The TRY initiative and database is hosted, developed and maintained by J. Kattge and G. Bönisch (Max Planck Institute for Biogeochemistry). TRY is currently supported by DIVERSITAS/Future Earth and the German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig. This work is part of the ALPVEG network (www.alpveg.com/).

Vegetation dynamics at the upper elevational limit of vascular plants in Himalaya. Scientific Reports, 6, 1-13. https://doi.org/10.1038/srep24881 Favre, A., Päckert, M., Pauls, S. U., Jähnig, S. C., Uhl, D., Michalak, I., and Muellner-Riehl, A. N. (2015) The role of the uplift of the Qinghai-Tibetan Plateau for the evolution of Tibetan biotas. Biological Reviews of the Cambridge Philosophical Society, 90, 236-253. https://doi.org/10.1111/brv.12107 Fazayeli, F., Banerjee, A., Kattge, J., Schrodt, F., and Reich, P. B. (2014) Uncertainty quantified matrix completion using bayesian hierarchical matrix factorization. In: 2014 13th International Conference on Machine Learning and Applications IEEE, pp. 312-317. https://doi.org/10.1109/ICMLA.2014.56 Fernández-Pascual, E., Carta, A., Mondoni, A., Cavieres, L., Rosbakh, S., Venn, S., et al. (2020) The seed germination spectrum of alpine plants: a global meta-analysis. New Phytologist (early view). https://doi.org/10.1111/nph.17086 Ferrier, S., Manion, G., Elith, J., and Richardson, K. (2007) Using generalized dissimilarity modelling to analyse and predict patterns of beta diversity in regional biodiversity assessment. Diversity and Distributions, 13, 252-264. https://doi.org/10.1111/j.1472-4642.2007.00341.x Garnier, E., Navas, M.-L., and Grigulis, K. (Eds) (2016) Plant functional diversity. Plant Functional Diversity. Oxford: Oxford University Press. https://doi.org/10.1093/acprof:oso/9780198757368.001.0001 Geange, S. W., Pledger, S., Burns, K. C., and Shima, J. S. (2011) A unified analysis of niche

Author contributions: B. J.-A., L. M., R. T. and C. P. C. conceived the study and developed the methodology. P. B., J. D., M. F., A. H., A. Y. K., J. L., N. M., G. P. M., J. N., A. N., R. K. P., G. P., J. Š., P. S., K. V., R. V., S. K. W. and E. G. Z. provided the vegetation plot data. F. M. S. and H. B. facilitated access to the sPlot database. R.T. analysed the data and produced the outputs. R.T and B.J.-A wrote the first manuscript draft. B.J.-A supervised the study. All the authors discussed the methodology and commented on various versions of the manuscript.

Data availability:

The vegetation plot data is stored and managed by sPlot and the corresponding author and is available upon request. The R code used to carry out the analyses is available at: https://dx.doi.org/10.6084/m9.figshare.14040152


Benoit Tallon

Philippe Roux
email: philippe.roux@univ-grenoble-alpes.fr

Guillaume Matte

Jean Guillard

John H Page

Sergey E Skipetrov

Ultra-Slow Acoustic Energy Transport in Dense Fish Aggregates Ultra-Slow Acoustic Energy Transport in Dense Fish Aggregates

A dramatic slowing down of acoustic wave transport in dense fish shoals is observed in open-sea fish cages. By employing a multi-beam ultrasonic antenna, we observe the coherent backscattering (CBS) phenomenon. We extract key parameters of wave transport such as the transport mean free path and the energy transport velocity of diffusive waves from diffusion theory fits to the experimental data. The energy transport velocity is found to be about 10 times smaller than the speed of sound in water, a value that is exceptionally low compared with most observations in acoustics. By studying different models of the fish body, we explain the basic mechanism responsible for the observed very slow transport of ultrasonic waves in dense fish shoals. Our results show that, while the fish swim bladder plays an important role in wave scattering, other organs have to be considered to explain ultra-low energy transport velocities.

INTRODUCTION

Because of their swim bladder (analogous to an immersed air bubble), bony fish (Osteichthyes) are very strong scatterers for underwater acoustic waves. Thus, (multi-beam) sonar techniques are very efficient to locate and characterize shoals, aggregates or even isolated fish. Most of the fisheries acoustics methods are developed under the single scattering approximation [START_REF] Simmonds | Fisheries Acoustics: Theory and Practices[END_REF] , i.e. for low fish concentrations when the wave is scattered at most once during propagation. For such low densities (∼ 1 fish/m 3 ), fish counting is straightforward and efficient using traditional methods such as echo-counting or echo-integration 2 . However, fish aggregates can be very dense (∼ 50 fish/m 3 ) for aquaculture purposes, for example, and in naturally occurring fish schools. In these cases, the backscattered signal received by the antenna is scattered by several fish during the wave's propagation, which makes fish counting much more challenging. This multiple scattering regime is problematic for the aquaculture industry, for which nonintrusive biomass estimation is one of the most important issues in regular practice 3 . In a previous work 4 , we suggested the use of mesoscopic physics and, in particular, of multiple scattering theory based on the diffusion approximation to deal with wave propagation in dense fish shoals.

In the diffusion approximation, one assumes that after a propagation distance corresponding to the transport mean free path * , the average intensity of multiply scattered waves follows a diffusion process (such as in heat diffusion) with a characteristic diffusivity D = v e * /3. The diffusivity involves the energy transport velocity of diffusive waves v e = 3D/ * . The energy velocity is thus a key parameter describing wave transport, and is proportional to the ratio of energy flux to energy density. For diluted or non-resonant systems, v e v 0 (where v 0 is the sound speed in water). However, as demonstrated with previous model systems, the energy velocity can be highly impacted by resonant phenomena that are typically encountered when the acoustic wavelength becomes similar to the typical size of the scatterers [5][6][7] .

In this paper, we first report an observation of an ultra-low value of the energy transport velocity of diffusive acoustic waves in a dense fish shoal through coherent backscattering (CBS) measurements [8][9][10][START_REF] Akkermans | Mesoscopic Physics of Electrons and Photons[END_REF] . CBS is a wave interference phenomenon that manifests itself by an enhancement (by a factor of 2) of the average backscattered intensity measured in the direction opposite to that of the incident wave. The angular profile of backscattered intensity has a cone shape in the stationary (continuous wave) limit, with a width that depends on * . In the dynamic case, accessed using a pulsed source, the temporal evolution of the backscattering peak's width depends on D and hence on v e . In order to scan the angular dependence of the backscattered intensity, we employ a multi-beam sonar probe (based on the Seapix technological brick [START_REF] Mosca | [END_REF] ) in a large fish cage anchored in open sea. In the second part of the paper, we present a comparative study of energy transport velocity calculations based on Mie theory 13,14 . The comparison of four scattering models for individual fish reveals that, while fish are usually approximated as air bubbles in water for scattering of acoustic waves, their complex structure can play an important role for wave transport in dense shoals. This comparison enables us to identify the essential features that need to be accounted for, and to propose a simple model to replicate the scattering properties of fish in dense shoals.

RESULTS

Experiments

Coherent backscattering (or weak localization) is a mesoscopic phenomenon that has been observed for light 8,9 , ultrasound 10 , matter waves 15 and seismic waves 16 . This effect is due to the constructive interference of waves following time-reversed pairs of paths. This interference produces a peak, centred on the exact backscattering direction, in the angular profile of backscattered intensity. In order to observe CBS, we use a Mills Cross multi-beam antenna made of two perpendicular ultrasonic arrays (2×64 transducers with a central wavelength λ = 1 cm). Experiments involve sending a short acoustic pulse (central frequency 150 kHz) into the fish cage by firing all transducers at the same time and performing an angular scan of the backscattered intensity using the beamforming method 17 . The cage contains a shoal of gilthead sea breams (Sparus aurata) with an average mass of 150 g and concentration η ∼ 50 fish/m 3 .

Time integration of the backscattered intensity yields the stationary CBS profile, which has an angular width of ∆θ ∼ λ / * (where λ is the wavelength in water). We fit the entire profile with the predictions of diffusion theory 11 ; see Fig. 1a. The fitting parameters are the transport mean free path * = (6.0 ± 0.2) cm and the absorption length a = (100 ± 3) cm, which characterizes the exponential decay of acoustic energy due to losses. The value k * = 36 1 (with k = 2π/λ ) indicates that even if sound is strongly scattered in this shoal, no complex phenomena such as strong localization (occurring for k * ∼ 1) impact the CBS peak shape.

The time-resolved dynamic CBS peak narrows with time (Fig. 1b). From diffusion theory [START_REF] Akkermans | Mesoscopic Physics of Electrons and Photons[END_REF] , the CBS width depends on diffusivity as ∆θ ∝ 1/ √ Dt. Using the measured value of * and the diffusion theory, the fitting of the dynamic CBS peak gives us D = (1.3 ± 0.1) m 2 /s. In this way, the simultaneous measurement of * and D leads to the energy velocity v e = 3D/ * = (65 ± 5) m/s. This value of v e is surprising low (an order of magnitude lower than the sound speed in water v 0 = 1500 m/s) , which is a very rare observation for acoustic waves 4,6,7,18,19 . Similarly slow diffusion has been observed previously with CBS measurements in dense shoals of sea breams, sea basses (Dicentrarchus labrax) or croakers (Argyrosomus regius) 4 , but no such extreme behaviour has been found for acoustic waves in other multiply scattering media 6,7,18,19 . However, unlike this past work on dense fish shoals, the cage considered in this paper has a sufficiently low fish concentration to allow us to neglect mesoscopic interferences that might impact the energy velocity, thereby facilitating a quantitative interpretation of the current results. In the following section, we employ a model 14 that takes into account the scattering delay induced by the fish and explains this ultra low value of energy velocity.

Energy velocity calculation

Here, we employ a microscopic model of the energy velocity 14 in order to identify the mechanisms responsible for the ultra low v e value. This description is based on a model that accounts for the delay ∆t ave induced by an immersed scatterer, and predicts a result for v e that may be approximated as 19 1

v e ≈ 1 v gr + ησ T ∆t ave . (1)
Here v gr is the group velocity of the average wave field (v gr ∼ v 0 in the absence of dispersion effects) and η is the scatterer concentration. ∆t ave is the scattering delay time for a single scattering event, calculated from the intensity-weighted angle-averaged phase derivative with frequency of the scattering amplitude [equation ( 4)], and σ T is the total scattering cross section. Thus, to obtain slow diffusive waves, the system has to be dense (high η) with strong scatterers (high σ T ) and large scattering delay (large ∆t ave ).

In order to identify the parts of the fish that are responsible for the ultra-low energy velocity, we use equation ( 1) to compute v e for four idealized spherically symmetric scattering models representing different simplifications of the complex fish body (Fig. 2). Since this theory enables analytic expressions to be obtained for relatively simple scattering geometries, the goal here is to capitalise on this capability to search for the key features that need to be considered in order to understand the origin of the remarkably slow energy velocity. Thus, rather than attempting a complex simulation that might obscure the basic scientific mechanism(s) at play, we focus on very simple models to reveal the basic scattering mechanisms involved:

• Model a: an air bubble representing the fish swim bladder with radius R 1 = 10 mm (Fig. 2a), longitudinal wave speed v l1 = 340 m/s and density ρ 1 = 0.001 g/cm 3 .

• Model b: a homogeneous soft sphere representing the fish flesh with radius R 1 = 76 mm, longitudinal wave speed v l2 = 1600 m/s, shear wave speed v t2 = 10 m/s and density ρ 2 = 1.1 g/cm 3 (Fig. 2b).

• Model c: a combination of models a and b representing the swim bladder surrounded by a flesh layer (Fig. 2c).

• Model d: similar to model c with an additional hard thin layer representing fish scales and bones with radius R 3 = 78 mm, longitudinal wave speed v l3 = 1600 m/s, shear wave speed v t3 = 900 m/s and density ρ 3 = 1.4 g/cm 3 (Fig. 2d).

Some scattering theories allow calculations for spheroids [START_REF] Flammer | Spheroidal wave functions[END_REF] that might be closer to the actual fish shape, but none of those calculations are developed for multilayer scatterers. However, the spherical approximation is suitable in the present case because of the randomized fish orientation in the azimuthal plane. Thus on average, the effective scatterer shape seen by the incident plane wave can be approximated as a sphere. Model a represents the usual assumption in fisheries acoustics (in the single scattering regime): since, at least near resonance, the swimbladder is the most reflective organ for acoustic waves, fish shoals are often seen as clouds of air bubbles in water. However, as shown in Fig. 3a, Eq. ( 1) applied to model a fails to explain the ultra-low value of energy velocity found in our experiments. The same conclusion can be drawn for models b and c, for which we also obtain v e ∼ v 0 . On the other hand, model d predicts a very slow energy velocity, v e ≈ 100 m/s, at the frequency f = 150 kHz, and ranges from about 50 to 150 m/s over the bandwidth of the transducers. In comparison with the other models, these values are very close to our measurement (v e = 65 ± 5 m/s). The phase and group velocities (see Supplementary Information, Fig. S2) can also be calculated for these models. Both of these velocities are found to be very close to v 0 for all four models. These results indicate that the complex fish structure mostly impacts only the diffusive waves; the ballistic wave velocities (the "line-of-sight" propagation through the shoal) are not significantly affected by the fish scattering and concentration. Thus, the hard thin layer surrounding the soft solid representing fish flesh seems to play an important role for the slowing down of diffusive acoustic waves. Figure 3(b) shows the calculation of energy velocity for a range of fish concentrations η and as a function of size variation ∆R/R (here the size variation ∆R/R is a factor that is equally applied to the three radii R 1 , R 2 and R 3 ). It is important to note the weak size dependence of energy velocity on ∆R/R that proves that the drop of v e /v 0 in model d with respect to model c is not due to size differences (outer radius R 3 versus R 2 ). In contrast, the fish concentration η seems to have a significant impact on v e . This effect could be interesting for enabling a new method of biomass assessment with acoustic waves compared with more traditional acoustic methods (see, for example, Ref. [START_REF] Foote | [END_REF] ).

DISCUSSION

In this section, we interpret the role of each part of the fish model on the scattering delay. The double core-shell structure of model d seems to explain the observed ultra-low energy velocity. As has been observed in the past [START_REF] Liu | [END_REF] , core-shell scatterers can indeed exhibit very strong scattering. However, no slow diffusive waves have been measured for such scatterers in the past. Resonant mechanisms have been identified [5][6][7] as being responsible for the decrease of v e in both optics and acoustics, but only for homogeneous scatterers. These effects result in large energy velocity variations, with v e ∼ v 0 far from resonant frequencies and v e v 0 around resonances. In the present case, such frequency fluctuations (expected over several tens of kHz) cannot be observed in our experiments, since the sonar bandwidth is too narrow; hence, frequency-resolved measurements were not feasible, and our experimental determination of v e corresponds to a narrow bandwidth-limited average around the central frequency (150 kHz) of the transducers. Thus, from experimental observations, we can measure conclusively the low v e value but we cannot conclude anything about its potential variations over a larger frequency range.

A way to interpret the scattering delay impacting the energy velocity is the calculation of the acoustic energy density inside and outside the scatterer 19,22 . High values of energy density suggest that waves are "stored" in the scatterer. This energy is then released in the surrounding medium with a certain delay resulting in a slowing down of diffusive wave transport. Figure 4 shows the energy density calculations for longitudinal waves at the frequency f = 150 kHz for the four different model scatterers.

The energy density calculation for an air bubble (model a) exhibits predominantly forward scattering as expected at high frequency for small scatterers with large acoustic contrast [START_REF] Van De Hulst | Light scattering by small particles[END_REF] . Since the frequency being considered in this study is far from the resonance of the bubble ( f res ∼ 0.5 kHz), the scattering is very weak for this model, and the energy density inside the scatterer is small. Thus, one might expect only a modest scattering delay, and the scattering strength (or total cross section σ T ) induced by the isolated swim bladder is not high enough for the product of scattering delay and cross section to cause a significant decrease of v e . For models b and c, both the energy density inside the scatterers and the forward scattering are somewhat larger than for model a, suggesting that there could be a bigger difference between group and energy velocities, but additional information would be needed to assess if they could be interesting candidates for predicting a slow energy velocity. However, for the double core-shell system (model d), the result of the energy density calculation is much more striking, as we observe a very strong increase of scattered energy density (the large increase in wave energy stored in the fish body is due to the hard scales and bones layer). The subwavelength-scale outer layer helps the generation of slow shear waves via mode conversion and the trapping of both longitudinal and shear waves in the fish flesh. This stored acoustic energy is then re-radiated into the surrounding water with a large delay 14 . The link between the stored energy and the large scattering delay becomes clear when comparing Fig. 4 and Fig. 5, which shows the angular dependence of the scattering delay ∆t(θ ) for all models [Fig. 5(a)-(d)] (see the Supplementary Information for calculation details). In particular, this figure shows that the angle-resolved scattering delays are much larger, typically by a couple of orders of magnitude, for model d than for the other models, for which the energy densities inside the scatterers are much less. Specifically, for model a, the delay is relatively small negative at all scattering angles, whereas for models b and c, the delays are negative for most angles and have, typically, somewhat larger magnitudes. In all three cases, these results suggest that the total angle averaged delay will be fairly small and certainly negative (indicating a slight enhancement of the energy velocity relative to the group velocity). For model d, however, very large delays, both positive and negative, are seen, with the positive delays dominating. These observations are confirmed by doing the angular integration of ∆t(θ ), yielding average scattering delays per wave period ∆t ave /T 0 for models a through d of -0.4, -1.4, -1.4 and 80.8, respectively. Thus, we find that adding the hard coating increases the magnitude of the angle-averaged delay by a factor of approximately 100 or more in comparison with the other models.

To further illustrate this large scattering delay, we calculate the dynamic energy density of the scattered waves W s outside the scatterers (Fig. 6). W s is obtained by integrating W ( f ) (c.f., Fig. 4 at f = 150 kHz) over the bandwidth of the incident pulse. The resulting temporal evolution of W is spatially averaged on the region R 3 < r < b (where r = x 2 + y 2 and b = 85 mm is the average inter-scatterer distance for the concentration η considered here. In this way, we obtain the dynamics of the acoustic energy density in the surroundings of each scatterer. The higher overall energy obtained with model d confirms the strong influence on the scattering induced by the presence of scales and bones. Furthermore, while models a, b and c all predict that W d decreases quickly (which explains v e ∼ v 0 for these systems), a much slower decay of W d with time is obtained with model d. This slow decay for model d demonstrates clearly that the large stored acoustic energy in the fish body (Fig. 4(d)) is slowly radiated into the surrounding medium (Fig. 6). To summarize, the core-shell model system with swim bladder + flesh + scalesand-bones layer leads to a large slowly decaying energy density associated with a large scattering delay. These observations and the quite large fish concentration η successfully explain the slow diffusion of acoustic waves observed in the sea bream shoal.

1 0 -4 1 0 -2 1 0 0 1 0 2 1 0 4 1 0 4 1 0 2 1 0 0 1 0 -2 1 0 -4 1 0 -2 1 0 0 1 0 2 1 0 4 0 0 . 2 5 0 . 5 0 . 7 5 1 1 0 4 1 0 2 1 0 0 1 0 -2 0 0 . 2 5 0 . 5 0 . 7 5 1 | ∆t ( q ) | ( µs ) f o r ∆t ( q ) > 0 ∆t ( q ) > 0 ( a ) | ∆t ( q ) | ( µs ) f o r ∆t ( q ) < 0 ∆t ( q ) < 0 ∆t ( q ) > 0 ( b ) ∆t ( q ) < 0 | ∆t ( q ) | ( µs ) f o r ∆t ( q ) > 0 ∆t ( q ) > 0 ( c ) | ∆t ( q ) | ( µs ) f o r ∆t ( q ) < 0 q ( r a d / π) ∆t ( q ) < 0 ∆t ( q ) > 0 ( d ) q ( r a d / π) ∆t ( q ) < 0

CONCLUSION

In conclusion, while the assumption that fish shoals are equivalent to air bubble clouds is adequate for most acoustic studies 1 , some cases with dense fish shoals have to be considered more carefully. Indeed, high fish concentration (η ∼ 50-100 fish/m 3 ) can lead to an accumulation of scattered energy, a resulting increase in scattering delay and a drastic slowing down of sound diffusion. Thus, it is essential to properly estimate v e in fish shoals in order to extend the range of application of acoustic fisheries techniques to very dense shoals. Several improvements (such as FIG. 6: Temporal evolution of scattered energy density W s (plots are normalized by the total incident energy density E i ). The solid blue line represents the energy density of incident pulse (without a scatterer). The time axis is normalized by the width of the incident pulse T 0 = 0.1 ms. scatterer shape and structure) can, in principle, be implemented in the model in order to calculate more accurate v e values. Nonetheless, as a first study, the simple models used here are sufficient to reveal the essential wave physics behind our observations, and enable us to demonstrate the strong impact of coating layers surrounding the swim bladder on wave transport in dense shoals.

As shown in Fig. 3b, the strong energy velocity dependence on the fish concentration could be useful for fish counting purposes. The method is particularly interesting since v e depends strongly on concentration and only weakly on fish size variation (over a reasonable size range for aquaculture conditions: ∆R/R ± 15%). In particular, such a tool could help the aquaculture industry for which large fish concentrations make acoustic biomass estimation impossible using traditional acoustic approaches based on the assumption of single scattering. Noninvasive monitoring of fish farms using diffusive waves 4,[START_REF] Tallon | [END_REF] is currently under study with long-term experiments (several weeks) in order to investigate the impact of fish size and biomass variations on the diffusive transport of sound waves.

METHODS

The sea cages are located in the Mediterranean Sea (Cannes, France) where the water temperature is about 20 • and salinity 3.6%. The cage in which experiments were conducted is cubic with a volume of 125 m 3 . The distance from the bottom of the cage to the sea bottom is z = 6.5 m. To maintain the organic label of the farm and to avoid the need for drug treatments, the fish densities in these cages are lower than those in intensive farming facilities (where mass densities can reach 100 kg/m 3 ). The feeding procedures are controlled to obtain a calibrated fish size.

The theory used to fit the experimental CBS data is derived from a diffusion equation for the average intensity I(r,t) [START_REF] Akkermans | Mesoscopic Physics of Electrons and Photons[END_REF] that is solved for a disordered medium occupying the half-space z > 0, with a delta-function source at z = z = * 25 :

I(t) = I 0 2π +∞ -∞ z 0 exp(-γ 0 z ) D(1 + γ 0 z 0 ) exp(-iΩt)dΩ, (2)

INTRODUCTION

This document provides additional information on the experiments that were performed to investigate dense fish shoals in open-sea fish cages, and on the theory used to interpret the results of these measurements.

EXPERIMENTAL SET-UP

The CBS experiments were performed using a Mills cross-shaped antenna that consists of two perpendicular linear arrays, each with 64 ultrasonic piezoelectric transducers. Each transducer is narrow band (bandwidth, 10 kHz) with a central frequency of 150 kHz, which corresponds to a wavelength λ of 1 cm in water. The transducers are square in shape and have sides of 0.5 cm λ /2. The antenna (expressly designed by the iXblue company for aquaculture monitoring) is placed just below the water surface, facing the sea bottom.

At full power, the radiation pressure of the ultrasonic antenna is 133 dB•µPa (measured at a distance d = 1 m from the source and for the frequency f = 150 kHz). The hearing threshold for sciaenidae fish is about 60 dB•µPa in the frequency range f ∈ [0.1 -1] kHz S1 . No avoidance behavior of the fish has been observed during our experiments. For comparison, some vessel noise signatures are reported in Ref. S2 and turn out to be around 130 dB•µPa in the hearing frequency range of sciaenidae fish.

DIFFUSION THEORY

As indicated in the Methods section of the main paper, the theory used to interpret the experimental CBS data is obtained from a diffusion model for the average intensity S3 I(r,t) in a disordered medium occupying the half-space z > 0 when the source is a delta-function at z = z = * S4 . The expression for I(r,t) is repeated here for convenience:

I(t) = I 0 2π +∞ -∞ z 0 exp(-γ 0 z ) D(1 + γ 0 z 0 ) exp(-iΩt)dΩ, (S1)
where γ 2 0 (Ω) = -iΩ D + 1 Dτ a , τ a is the characteristic absorption time and z 0 = 2 3 1+R 1-R * is the extrapolation length, with R = 0.99 as the reflection coefficient for the water/air interface. The theoretical expression used to fit the dynamic CBS profile I(θ ,t) follows from the same diffusion theory, and is given by

I(θ ,t) = I 0 2π +∞ -∞ z 0 D e -γ 0 z 1 + γ 0 z 0 + e -γz 1 + γz 0 exp(-iΩt)dΩ, ( S2
)
where γ 2 (θ , Ω) = -iΩ D + k 2 0 sin 2 (θ ) + 1 Dτ a and γ 0 = γ(θ = 0, Ω). The expression for the stationary CBS profile is obtained by integrating equation (S2) over time S5 .

PHASE AND GROUP VELOCITIES

We outline the calculation of the wave number k for the average field. This wave number, useful for the energy velocity calculation, is expressed as :

k 2 = k 2 0 + 4πηF(θ = 0), ( S3
)
where η is the scatterer concentration and k 0 = ω/v 0 is the wave number of the incident wave at frequency f = ω/2π in water. In the far field approximation, the scattering function F(θ ) = |F(θ )|exp[iϕ(θ )], with magnitude |F(θ )| and phase ϕ(θ ), represents the scattering amplitude in the direction given by the angle θ with respect to the incident wave vector. F(θ ) is given by the following expression

F(θ ) = 1 ik 0 ∑ n (2n + 1)A n P n (cos θ ), (S4)
where P n represents the Legendre polynomials and A n the scattering amplitude coefficients of the scattered field. The A n coefficients are obtained by solving the Mie problem S6 : by imposing stress and displacement continuity conditions at the interfaces of a spherical scatterer (model a), one obtains a system of 3 linear equations with 3 coefficients S7 . The solution of this system of equations requires inverting a 3 × 3 matrix, and leads to all of the 3 partial wave amplitudes (the scattered longitudinal wave, the refracted longitudinal wave and the refracted shear wave). The calculation of scattering amplitude of a coated sphere (model c) is obtained by solving a system of 7 linear equations with 7 coefficients, as detailed in Ref. S8. Following the same procedure, the last case of a double coated sphere (model d) requires the solution of a system of 10 linear equations with 10 coefficients given by stress and displacement continuity conditions at the 3 interfaces.

FIG. S1: Scattering cross section calculated for the four model scatterers.

The scattering function allows the calculation of the scattering cross-section σ to be performed:

σ = 2π dθ sin θ |F(θ )| 2 . ( S5
)
The scattering cross-section of our four spherical scatterers (models a-d) is represented in figure S1. As expected from the energy density calculations (Fig. 4), the scattering from the simple air bubble is weak in comparison with the more elaborate spherical models b, c and d. The significantly larger scattering strengths for the three more complicated spherical models indicate that they are more promising candidates for achieving slower wave transport, especially model d. However it is evident that scattering strength alone is not sufficient to explain the energy velocity predictions obtained from these models (see Fig. 3(a) of the main paper). The knowledge of the wave number k leads to the calculation of the phase velocity v p = ω/Re[k] and the group velocity v gr = ∂ ω/∂ Re[k]. Both these velocities are plotted in figure S2 for the four model scatterers. In all cases, both of these velocities are very close to the water velocity v 0 . These calculations reveal the weak phase shift of the forward scattered waves (since v p and v gr only depend on F(θ = 0) and ∂ Re[F(θ = 0)]/∂ ω, respectively). In terms of energy transport, these results indicate that the transport velocity of the average wave field (v gr ) is barely affected by the scattering. Thus, the ultra-slow transport of diffusive waves must result only from the very large scattering delay of multiply scattered waves that is described in the next section.

ENERGY VELOCITY

The energy velocity calculation is based on the renormalization of the sound speed in pure water due to the scattering delay S10,S11 Under the approximation v 2 0 /v p v gr ≈ 1 (for weak dispersion of the average wave field), the energy velocity can be rewritten as S12 1 v e = 1 v gr + ησ ∆t ave , (S8)

where the average scattering delay ∆t ave is equal to ∂ ω /σ . While the average-field velocities only involve the scattering function in the forward direction, the scattering delay of diffusive waves depends on the frequency derivative of phase shift in all scattering directions ∂ ϕ(θ )/∂ ω. The group velocity being similar to the sound speed in pure water (v gr ∼ v 0 ), it is this scattering delay that is responsible for the remarkably slow values of the energy velocity obtained for the multi-layer sphere (model d). Thus, the adding of the hard scales layer significantly slows down the diffusive wave transport by facilitating the storing of acoustic energy within the scatterer, and the releasing of this energy with a large delay.

FIG. 1 :

 1 FIG. 1: (a) Stationary CBS profile (blue open circles) fitted with the diffusion theory (red dashed line). (b) Dynamic CBS profile fitted with the diffusion theory (black solid lines). In (b), the intensity is normalized by its peak value at each time to show more clearly the temporal evolution of the width.

FIG. 2 :

 2 FIG. 2: On the left: Schematic representation of a sea bream body. On the right: Models a, b, c and d used for the energy velocity calculation.

FIG. 3 :

 3 FIG. 3: (a) Normalized velocity v e /v 0 calculations versus frequency for the four different models. (b) Normalized velocity v e /v 0 calculations versus fish size variation ∆R/R and concentration η at a given frequency f = 150 kHz for model d. The inset represents v e /v 0 versus fish concentration for ∆R/R = 0 and f = 150 kHz.

FIG. 4 :

 4 FIG. 4: Calculations of the longitudinal acoustic energy density W for the four different model scatterers at f = 150 kHz. W 0 is the incident energy density. The white dotted line represents the limits of the swim bladder, the dashed pink line the flesh layer and the solid black line the thin scales and bones layer. White arrows indicates the direction of the incident plane wave.

FIG. 5 :

 5 FIG. 5: (a) to (d) Angular dependence of the normalized scattering delays ∆t θ calculated for the four models (a, b, c and d), at a frequency of 150 kHz.

FIG. S2 :

 S2 FIG.S2: Phase (upper panel) and group (lower panel) velocities calculated for the four different models.

  the delay parameter δ is equal toδ = 2πηv gr v p ω ∂ ReF(0) ∂ ω + π 0 sin(θ )|F(θ )| 2 ∂ ϕ(θ ) ∂ ω dθ . (S7)

  ∆t ave = dθ sin θ |F(θ )| 2 ∂ ϕ(θ ) ∂ ω dθ sin θ |F(θ )| 2 = 2π dθ sin θ ∆t(θ ). (S9)The right side of this equation follows from the definition of the total scattering cross section [equation (S5)], and is expressed in terms of the angle-dependent scattering delay ∆t(θ) = |F(θ )| 2 ∂ ϕ(θ )

ACKNOWLEDGMENTS

The authors wish to thank Olivier Lerda from iXblue for Seapix signal processing routines. JHP is grateful for support from the Natural Sciences and Engineering Research Council of Canada's Discovery Grant Program (RGPIN-2016-06042).

where γ 2 0 (Ω) = -iΩ D + 1 Dτ a , τ a is the characteristic absorption time and z 0 = 2 3 1+R 1-R * is the extrapolation length, with R = 0.99 as the reflection coefficient for the water/air interface. Theoretical predictions for both the dynamic and stationary CBS peaks can then be obtained using this theory 4 .

The energy velocity theory 5 (Eq. 1) is based on the calculation of the scattering function F(θ ) = |F(θ )|e iϕ(θ ) with magnitude |F(θ )| and phase ϕ(θ ). F(θ ) represents the scattering amplitude in the direction given by the angle θ with respect to the incident wavevector, and is given by the following expression:

where k 0 represents the wavenumber of incident wave, P n the Legendre polynomials and A n the scattering amplitude coefficients of the scattered field. The A n coefficients are obtained by solving the Mie problem 13 , invoking the calculation of stress and displacement continuity conditions (for both longitudinal and shear waves) at the boundaries R 1 , R 2 and R 3 [START_REF] Faran | [END_REF] . F(θ ) is then used to obtain the scattering cross-section σ and delay ∆t ave that are needed for the calculation of equation ( 1):

Additional details on the experimental method and theoretical calculations can be found in the Supplementary Information. To whom correspondence should be addressed; E-mail: philippe.roux@univ-grenoble-alpes.fr