Irene Calderón-Sanou

Marc Ohlmann

Tamara Münkemüller

Lucie Zinger

Mickael Hedde

Clément Lionnet

Camille Martinez-Almoyna

Amélie Saillard

Julien Renaud

Nicolas Le Guillarme

Ludovic Gielly

Orchamp Consortium

Wilfried Thuiller

Mountain soil multitrophic networks shaped by the interplay between habitat and pedoclimatic conditions

Keywords: network dissimilarity, soil biodiversity, trophic interactions, environmental DNA metabarcoding, ecological filters, metaweb

WT, TM, MO and IC-S conceptualised the study. AS and WT led the field work, with participation of MO, CM-A, JR, CL, LG, IC-S and several members of the Orchamp consortium. LG executed eDNA extractions and all laboratory samples processing

Introduction

Soil biodiversity encompasses a complex network of interactions among functionally and trophically diverse organisms, playing a vital role in supporting ecosystem functions and services such as carbon sequestration, organic matter decomposition, and enhancing plant performance and resistance to pests and stress [START_REF] Bardgett | Belowground biodiversity and ecosystem functioning[END_REF][START_REF] Delgado-Baquerizo | Multiple elements of soil biodiversity drive ecosystem functions across biomes[END_REF]. Over the past decades, our understanding of the spatial distribution of soil taxa has improved, particularly for specific groups like bacteria [START_REF] Fierer | The diversity and biogeography of soil bacterial communities[END_REF][START_REF] Delgado-Baquerizo | A global atlas of the dominant bacteria found in soil[END_REF] and fungi [START_REF] Tedersoo | Global diversity and geography of soil fungi[END_REF][START_REF] Davison | Global assessment of arbuscular mycorrhizal fungus diversity reveals very low endemism[END_REF]. However, our current knowledge primarily relies on a few specific taxa, while the broader ecological processes shaping the whole soil communities and their organisation into interaction networks remain poorly understood. This knowledge gap stem from challenges associated with capturing the intricate spatial heterogeneity inherent to the soil environment [START_REF] Ettema | Spatial soil ecology[END_REF], and the methodological complexities of sampling and analysing the vast diversity of the soil biota [START_REF] Geisen | A methodological framework to embrace soil biodiversity[END_REF]. These limitations have hindered the application of existing theoretical frameworks on soil biota [START_REF] Thakur | Towards an integrative understanding of soil biodiversity[END_REF]. Advancing research in this direction promises to unravel the ecological processes that structure soil biodiversity and to predict the impacts of global change on terrestrial ecosystems [START_REF] Soliveres | Biodiversity at multiple trophic levels is needed for ecosystem multifunctionality[END_REF][START_REF] Eisenhauer | The multidimensionality of soil macroecology[END_REF].

Gaining insights into the organization of soil communities necessitates studying spatial, environmental, and biotic drivers, and this poses challenges [START_REF] Decaëns | Macroecological patterns in soil communities[END_REF]Münkemüller et al., 2020;[START_REF] Eisenhauer | Frontiers in soil ecology-Insights from the World Biodiversity Forum 2022[END_REF]. First, capturing the full complexity of soil systems requires accounting for their inherent spatial heterogeneity from landscape to local scales [START_REF] Ettema | Spatial soil ecology[END_REF][START_REF] Thakur | Towards an integrative understanding of soil biodiversity[END_REF]. At the landscape scale, dominant habitat types, such as grasslands or forests, are crucial for soil taxonomic and functional composition [START_REF] Fiore-Donno | Contrasting Responses of Protistan Plant Parasites and Phagotrophs to Ecosystems, Land Management and Soil Properties[END_REF][START_REF] Arribas | The limited spatial scale of dispersal in soil arthropods revealed with whole-community haplotype-level metabarcoding[END_REF][START_REF] Sepp | Woody encroachment in grassland elicits complex changes in the functional structure of above-and belowground biota[END_REF]. Consequently, the heterogeneity between habitats within a landscape can give rise to mosaic patterns of soil communities. Moreover, the spatial connectivity of habitats and the varying dispersal capacities of soil organisms can influence the local composition and structure of soil communities [START_REF] Arribas | The limited spatial scale of dispersal in soil arthropods revealed with whole-community haplotype-level metabarcoding[END_REF]. But even at the local scale of a few meters, variable abiotic conditions, including microclimate and soil physico-chemical properties, along with stochastic processes, can lead to pronounced differences in soil communities [START_REF] Ramirez | Biogeographic patterns in below-ground diversity in New York City's Central Park are similar to those observed globally[END_REF][START_REF] O'brien | Spatial scale drives patterns in soil bacterial diversity[END_REF][START_REF] Zinger | Body size determines soil community assembly in a tropical forest[END_REF].

The structure of the vegetation within an habitat or its taxonomy or functional composition can also affect the abundance and diversity of different soil taxa or trophic groups [START_REF] Fitzpatrick | Relative importance of tree species richness, tree functional type, and microenvironment for soil macrofauna communities in European forests[END_REF][START_REF] Noguerales | Community metabarcoding reveals the relative role of environmental filtering and spatial processes in metacommunity dynamics of soil microarthropods across a mosaic of montane forests[END_REF][START_REF] Calderón-Sanou | Energy and physiological tolerance explain multi-trophic soil diversity in temperate mountains[END_REF]. Effectively disentangling the effects of spatial, environmental and biotic processes necessitates a sampling design that encompasses multiple spatial scales, ranging from the heterogeneity between habitats to smallscale soil variations.

Second, comprehensive sampling and cross-taxa analysis are required to capture the high diversity in soils. Recent advancements in monitoring techniques, such as environmental DNA metabarcoding (eDNA) [START_REF] Taberlet | Environmental DNA[END_REF], offer cost-effective means to obtain extensive data on soil biota at large spatial extent [START_REF] Taberlet | Environmental DNA: for biodiversity research and monitoring[END_REF][START_REF] Geisen | A methodological framework to embrace soil biodiversity[END_REF]. To more concisely present the different taxa they can be classified into trophic groups that share similar resources or prey (Eltonian niche; Elton 1927), e.g. fungivorous nematodes or photoautotrophic bacteria [START_REF] Louca | Decoupling function and taxonomy in the global ocean microbiome[END_REF][START_REF] Potapov | Feeding habits and multifunctional classification of soil-associated consumers from protists to vertebrates[END_REF]. To better comprehend the complex interactions, we can construct multitrophic networks that consider multiple trophic groups (nodes) and their linkages across trophic levels (links) [START_REF] Thompson | Food webs: reconciling the structure and function of biodiversity[END_REF][START_REF] Gravel | Bringing Elton and Grinnell Together: a quantitative framework to represent the biogeography of ecological interaction networks[END_REF]. However, constructing multitrophic networks from large eDNA datasets is not without challenges. It entails handling and classifying thousands of sequences into trophic groups, and subsequently linking these groups based on known interactions. Obstacles include the resolution limitations of DNA markers for species-level assignment, the incompleteness of sequence reference databases, inconsistent terminology in trophic ecology across soil taxa [START_REF] Hedde | A common framework for developing robust soil fauna classifications[END_REF], and the sparse trophic information available for soil organisms in literature and repositories (but see [START_REF] Potapov | Feeding habits and multifunctional classification of soil-associated consumers from protists to vertebrates[END_REF]. Nevertheless, the application of artificial intelligence tools now offers the potential for automating classification tasks [START_REF] Compson | Chapter Two -Linking DNA Metabarcoding and Text Mining to Create Network-Based Biomonitoring Tools: A Case Study on Boreal Wetland Macroinvertebrate Communities[END_REF]Le Guillarme and Thuiller, 2023).

Third, if we are to understand soil biodiversity, not as a group of independent taxa, but rather as multitrophic networks, we need the appropriate methods and frameworks to investigate how ecological processes act on soil trophic groups and their interactions. The mono-trophic community assembly framework [START_REF] Keddy | Assembly and Response Rules: Two Goals for Predictive Community Ecology[END_REF][START_REF] Thuiller | A road map for integrating eco-evolutionary processes into biodiversity models[END_REF] can be extended to multitrophic networks using the metaweb concept [START_REF] Dunne | The network structure of food webs[END_REF]. At the regional scale, the metaweb represents the potential ecological network encompassing all species (or trophic groups) from the regional pool and their potential interactions (Fig. 1c). At the local scale, realised local networks emerge as subsets of the metaweb due to ecological filtering of species [START_REF] Bauer | Biotic filtering by species' interactions constrains food-web variability across spatial and abiotic gradients[END_REF]. This filtering process can be influenced by spatial filters like landscape barriers, limiting dispersal [START_REF] Peay | Evidence of dispersal limitation in soil microorganisms: Isolation reduces species richness on mycorrhizal tree islands[END_REF][START_REF] Gobbi | Hay meadows' overriding effect shapes ground beetle functional diversity in mountainous landscapes[END_REF][START_REF] Blasbichler | Abrupt boundaries between mountain meadows and forests separate ground-dwelling invertebrate communities: a case study from South Tyrol, Italy[END_REF], abiotic factors filtering out species lacking physiological adaptations [START_REF] Maaß | Environmental filtering vs. resourcebased niche partitioning in diverse soil animal assemblages[END_REF][START_REF] Glassman | Environmental filtering by pH and soil nutrients drives community assembly in fungi at fine spatial scales[END_REF][START_REF] Steinwandter | Ground-dwelling invertebrates of the high alpine: Changes in diversity and community composition along elevation (1500-3000 m)[END_REF], and biotic interactions such as mutualisms, herbivory and predation [START_REF] Vályi | Community assembly and coexistence in communities of arbuscular mycorrhizal fungi[END_REF]. By using the metaweb concept, we can evaluate the effects of ecological filters not only on the taxonomic composition of the local community but also on its network structure (i.e. composition of interactions). Important biotic interactions in the soil are via plant-soil feedbacks, where the taxonomic or functional composition of plant communities influences soil multitrophic networks and vice-versa (Kardol and De Long, 2018;Kardol et al., 2018). Moreover, trophic interactions within the network can elucidate its local structure and composition [START_REF] Thakur | Trophic Regulations of the Soil Microbiome[END_REF]. We anticipate that interactions between taxa or groups within the network should result in covariation across environmental and spatial gradients. For instance, environments characterised by decomposer communities dominated by bacteria (rather than fungi) should promote the dominance of bacteria consumers and related higher-level consumers within the multitrophic network [START_REF] Moore | Temporal and spatial heterogeneity of trophic interactions within below-ground food webs[END_REF], see e.g. [START_REF] Martinez-Almoyna | Differential effects of soil trophic networks on microbial decomposition activity in mountain ecosystems[END_REF]. Null models can be used to test the significance of such co-variation in the assembly of ecological communities [START_REF] Caruso | Fluctuating ecological networks: A synthesis of maximum-entropy approaches for pattern detection and process inference[END_REF].

Here, we analysed the spatial variation in soil multitrophic network structure, encompassing trophic group and interaction composition. We used data from all over the French Alps including lowland forests and high-elevation grassland ecosystems. To represent their spatial heterogeneity, we employed a stratified and nested sampling design, consisting of 24 elevation gradients with multiple plots at varying elevations (Fig. 1a). We measured soil biodiversity in 418 soil samples, using environmental DNA (eDNA). Using an ontology-based data integration pipeline, combining multiple databases with existing knowledge on the trophic habits of soil organisms (Le Guillarme and Thuiller, 2023;Le Guillarme et al., 2023), we constructed a metaweb comprising 55 soil trophic groups with 383 potential trophic interactions (Fig. 1b).

Finally, we applied the metaweb framework to characterise the importance of ecological filters on soil multitrophic networks (Fig. 1c), using a combination of network diversity indices based on Hill numbers [START_REF] Ohlmann | Diversity indices for ecological networks: a unifying framework using Hill numbers[END_REF] and associated null models. We quantified local variations in soil multitrophic network structure both between and within habitats, addressing three specific questions: (Q1) Do different habitats differ in soil multitrophic networks? (Q2) Do trophic interactions contribute to variations in group abundances between habitats?, and, (Q3) How do abiotic conditions, spatial factors, and plant communities shape soil multitrophic networks within habitats?

Material and Methods

Study site

The data come from the Orchamp long-term biodiversity observatory (www.orchamp.osug.fr, Appendix 1). It encompasses 24 elevation gradients (at the time of this study) across the French Alps sampled between 2016 and 2020 (Fig. 1a). These gradients represent diverse climatic, vegetation and soil conditions. Each gradient consists of four to nine 30 x 30 m permanent plots, spaced approximately 200 m apart in elevation. The elevation range covered by each gradient varies between 585 and 1380 m (Fig. S1). We worked with 127 plots, ranging in total from 280 m and 2780 m above sea level.

Habitat type (forest, grassland, shrubland) was determined on-site. Most gradients covered all three types of habitats, with some exceptions where only one or two types were present (e.g. grassland and shrubland in high-elevation mountains). Plant species abundances were estimated at the peak of the vegetation along a linear transect traversing each plot, using the pin-point method [START_REF] Jonasson | Evaluation of the Point Intercept Method for the Estimation of Plant Biomass[END_REF]. Soil sampling was performed at the end of the summer season in three subplots (2 x 2 m) located along a separate 4 m wide transect within each plot.

In each subplot, we collected ten soil cores, 5 cm in diameter and ~15 cm in depth (excluding litter), pooled together to a biological sample. Some elevation gradients (37 plots) were sampled twice (i.e. in different years).

Soil biota sampling

Soil biodiversity was measured using environmental DNA (eDNA) metabarcoding from the soil samples. In the field, we extracted eDNA from a 15 g aliquot following [START_REF] Taberlet | Environmental DNA[END_REF][START_REF] Taberlet | Environmental DNA: for biodiversity research and monitoring[END_REF]. The remaining soil was sieved at 2 mm and used to measure soil physicochemical properties (Table 1).

We employed six DNA markers, including two universal markers (Euka02 for eukaryote, Bact01 for bacteria) and four clade-specific markers (Fung02 for fungi, Inse01 for insect, Olig01 for oligochaeta, and Coll02 for collembola) [START_REF] Taberlet | Environmental DNA: for biodiversity research and monitoring[END_REF]. Details regarding the markers and molecular analyses are detailed in Appendix 2. A standardised bioinformatic pipeline was applied to curate the retrieved sequences from contaminants and errors [START_REF] Calderón-Sanou | From environmental DNA sequences to ecological conclusions: How strong is the influence of methodological choices?[END_REF], using the OBITools software [START_REF] Boyer | OBITOOLS : a UNIXinspired software package for DNA metabarcoding[END_REF] and the 'metabaR' R package (R Core Team, 2020;[START_REF] Zinger | metabaR: An r package for the evaluation and improvement of DNA metabarcoding data quality[END_REF]. Sequences were clustered into Molecular Operational Taxonomic Unit (MOTU) with SUMACLUST [START_REF] Mercier | SUMATRA and SUMACLUST: fast and exact comparison and clustering of sequences[END_REF] using a clustering threshold of 97% for Euka02, Fung01 and Bacte01, 85% for Coll01, 88% for Olig01 and 95% for Inse01, following the recommendations in [START_REF] Bonin | Optimal sequence similarity thresholds for clustering of molecular operational taxonomic units in DNA metabarcoding studies[END_REF]. Taxonomic annotations were performed on MOTUs with the SILVAngs pipeline [START_REF] Quast | The SILVA ribosomal RNA gene database project: improved data processing and webbased tools[END_REF] (SILVA version 138.1) for ribosomal universal markers. For clade specific markers, we used the ecotag program from the OBITools, and marker-specific databases built with the ecoPCR program [START_REF] Ficetola | An In silico approach for the evaluation of DNA barcodes[END_REF], from the EMBL database version 136. Metazoan taxa not registered in the European region were removed from the data (<1% of total reads) using the GBIFfilter tool (https://github.com/nleguillarme/gbif-filter-python).

Multitrophic network construction

Trophic groups and classes: We assigned the retrieved taxonomically annotated MOTUs to 10 trophic classes including autotrophs, decomposers, detritivores, phytophages or phytoparasites, plant mutualists, bacterivores, fungivores, omnivores, predators and zooparasites. We then refined those trophic classes by subdividing them into 55 trophic groups in total. These were defined by separating phylogenetic distant groups that could have a different set of prey/predators (e.g. bacterivorous mites vs. bacterivorous nematodes) and groups differing in their resources acquisition strategy (e.g. different types of mycorrhiza and saprotrophs). The taxonomic rank used to delimitate phylogenetic distant groups comprised Bacteria, Fungi, Protista, the different phyla within Metazoa, and the different classes or orders within Arthropoda and Annelida. (Appendix 2).

To facilitate the assignment of the taxonomically annotated MOTUs to trophic classes, we first built a knowledge graph integrating information about the trophic interactions and feeding habits of soil-associated consumers from multiple data sources using the ontology-based data integration pipeline described in Le Guillarme & Thuiller (2023). This trophic knowledge graph uses the NCBITaxon ontology and the Soil Food Web Ontology (SFWO) (Le Guillarme et al., 2023) as global schemas for reconciling taxonomic and semantic heterogeneities between the data sources. It provides a unified access to multisource trophic information across taxonomic groups and trophic levels. The list of data sources and details on the assignment criteria can be found in Appendix 2.

Building the metaweb: A metaweb is a potential network containing all trophic groups and their potential interactions [START_REF] Dunne | The network structure of food webs[END_REF]. We added three basal resources to construct our metaweb: light (or chemical energy), plants, and organic matter. Trophic links were then added between groups based on their main feeding preferences. Therefore, plant mutualists and phytoparasites were linked with the plant resource, detritivores and decomposers were linked with the organic matter resource and autotrophs were associated with the light resource.

Bacterivores were linked with all trophic groups containing bacteria, and fungivores with all trophic groups containing fungi. The trophic interactions of omnivores, predators and zooparasites, were determined through a literature review of dietary preferences of taxa representing over 90% of the read abundance of the group (Appendix 2).

Local networks: From the metaweb, the composition and structure of the local soil multitrophic networks were deduced based on locally present trophic groups, i.e. local networks are subgraphs of the metaweb. This assumes that co-occurring groups interact locally as defined in the metaweb. We approximated local relative abundance of each trophic group using a double-transformation, where first, total read counts were transformed into proportions within the sample, and second, the resulting proportions were standardised by the largest proportion observed across samples for each trophic group ('eDNA index' in [START_REF] Kelly | Understanding PCR Processes to Draw Meaningful Conclusions from Environmental DNA Studies[END_REF]. The relative abundance of each trophic group varied thus between 0 (absent) to 1 (largest observed proportion), allowing to have a comparable measure across trophic groups. Abundance values of the basal resources (light, plant, organic matter) were set to non-null (0.001) to not affect diversity measurements. For trophic class abundance estimates, we summed group relative abundances.

Network analyses

With the aim of analysing spatial variation in soil multitrophic network structure, we used beta diversity metrics that account for compositional changes in the network across samples. We used a set of network diversity indices based on the generalisation of Hill numbers proposed by [START_REF] Ohlmann | Diversity indices for ecological networks: a unifying framework using Hill numbers[END_REF]. These indices allow to quantify diversity in trophic groups and trophic interactions, varying the weight of the relative abundance of groups and links. Here, we focussed on Shannon diversity (q=1, see Calderón-Sanou et al. 2020) but also replicated the analyses using species richness (q=0) and Simpson diversity (q=2) to assess the extent to which the observed changes are compositional versus structural. Since interactions in the metaweb are binary (either presence or absent), we approximated the abundance of interactions as the product of the relative abundances of the interacting groups. These diversity indices can be further decomposed into alpha, gamma, and ß-diversity (ßT), and be used to calculate the pairwise dissimilarities for both groups (ßpP) and links (ßpL). ß-diversity (ßT) quantifies the effective number of different communities, based on shared groups (ßTP) or interactions (ßTL), ranging from 1 (indicating identical group or interaction abundance distribution) to the total number of communities, when networks do not share any common group or interaction. The diversity metrics were computed using the R package 'econetwork' [START_REF] Miele | econetwork: Analyzing Ecological Networks[END_REF].

To quantify changes in the spatial structure of soil multitrophic networks across habitats, we calculated pairwise dissimilarity measures for both groups (ßpP) and interactions (ßpL) among all pairs of local networks. Subsequently, we used UMAP, a non-linear dimension reduction algorithm, to visualize the local networks in a two-dimensional space based on their group and interaction dissimilarities [START_REF] Mcinnes | UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction[END_REF]. Given the strong structural differences observed in local soil multitrophic networks between forests and grasslands (Fig. 2,S3), and the high variability in network composition in shrublands due their intermediate state along the elevation gradient, our analysis focused exclusively on grasslands and forests for further comparison. We compared the relative abundances of trophic groups and classes between forests and grasslands with a Wilcoxon test and used the 'group-TL-tsne' network layout from the R package 'metanetwork' [START_REF] Ohlmann | metanetwork': a R package dedicated to handling and representing trophic metanetworks[END_REF] together with edge bundling provided by 'edgebundle' R package [START_REF] Schoch | edgebundle: Algorithms for Bundling Edges in Networks and Visualizing Flow and Metro Maps[END_REF] to represent the difference network between average networks in forests vs. grasslands.

We tested the imprint of trophic interactions on community assembly with null models. Firstly, we constructed an average soil multitrophic network for each habitat by using the average relative abundances of trophic groups per habitat. We then calculated the ß-diversity (ßT) for both group and link interactions between the two habitats. Next, we built the null model. Since group and interaction diversities are strongly correlated, we aimed at keeping group diversity constant in the null model while modifying interaction diversity. This was achieved by permuting the node labels within the metaweb. Consequently, group ß-diversity remains unchanged as this diversity index is invariant to label permutation. However, interaction ßdiversity was affected since it relies on the product of relative abundances of interacting groups.

If dominant interacting groups differ between the two habitats, the observed interaction ßdiversity would be higher than expected under a random distribution of group abundances. To evaluate the significance of our findings, we conducted 3000 permutations and compared the observed diversities with the corresponding null distribution.

Finally, to quantify the relative importance of environmental and spatial distances in explaining network dissimilarity within habitats (i.e. grasslands and forests), we used Generalised Dissimilarity Models (GDM) [START_REF] Ferrier | Using generalized dissimilarity modelling to analyse and predict patterns of beta diversity in regional biodiversity assessment[END_REF]. We built two GDM per habitat using trophic group turnover (ßpP) and interaction turnover (ßpL) as response variables, with the R package 'gdm' [START_REF] Fitzpatrick | Relative importance of tree species richness, tree functional type, and microenvironment for soil macrofauna communities in European forests[END_REF]. For the environmental predictors we selected a set of weakly correlated variables (Pearson's r < 0.5, Fig. S2) representing climatic, soil and vegetation conditions that could influence soil organisms (Table 1). The spatial coordinates of the samples and all the selected environmental variables were used as predictors. The environmental distance between samples was directly calculated by the 'gdm' package, but for plant composition, we used our own dissimilarity indices to represent changes in vegetation composition from both a taxonomic and functional point of view (Table 1). Plant taxonomic dissimilarity between plots was calculated using the beta.pair function from the R package 'betapart' [START_REF] Baselga | betapart: Partitioning Beta Diversity into Turnover and Nestedness Components[END_REF]. To estimate plant functional dissimilarity we retrieved trait values for species of which cumulative abundances represented at least 90% of the plot coverage. Missing values (<15% of the total data) were estimated using the imputation method offered by the R package 'mice' (van Buuren and Groothuis-Oudshoorn, 2011). Plant functional dissimilarity was estimated using the beta.fd.multidim function from the R package 'mFD' [START_REF] Magneville | mFD: Compute and Illustrate the Multiple Facets of Functional Diversity[END_REF]. In the analysis, soil samples from the same subplot collected in different years (i.e. 95 of the 418 samples) were treated as separate samples as they differed in environmental conditions, except for climate and NDVI that were averaged across a period of 10 years. Their spatial dependency was considered through the spatial coordinates of the plot.

Results

The metaweb constructed for the French Alps region using the taxa identified in 418 soil samples consisted of 55 trophic groups, three basal resources and a total of 383 potential interactions (Fig. S3). On average, local networks were composed of 41 ± 4 trophic groups (and resources), with 206 ± 37 interactions. All 10 trophic classes were present in all local networks, with a single exception where no group from the fungivore class was detected. The mean dissimilarity between local soil multitrophic networks was primarily influenced by changes in the dominance of trophic groups, rather than by changes in the presence or absence of groups (as indicated by increasing dissimilarity when increasing parameter q, which represents the weight of trophic group abundance and interaction, Fig. S4).

Overall, we found that multitrophic network composition varied across the environmental range covered in this study. The habitat type, particularly when comparing grassland and forest habitats, strongly structured the network composition, both in terms of trophic group and interaction composition (Q1) (Fig. 2, Fig. S5). The influence of habitat type on interaction composition was more pronounced compared to trophic group composition (Table S1). When examining shrubland areas, which represent ongoing shifts from grassland to forest, we found highly variable network structures (Fig. 2). These structures exhibited similarities to both forest and grassland networks, potentially due to the spatial proximity of shrubland sampling sites to forests or grasslands, or the specific stage of succession in the shrubland areas. By contrasting the average networks of forests and grasslands, we identified compositional differences in the abundance of trophic groups and classes between the two habitats (Fig. 3). Forest soil multitrophic networks were enriched in ectomycorrhizal fungi, litter and wood saprotrophic fungi, macro-detritivores like earthworms and diplopods at basal level, fungivores (particularly eudaphic collembola and protura), soil top predators including predatory mites, centipedes, predatory coleoptera and pseudoscorpions, as well as zooparasites. In contrast, grassland soil multitrophic networks had a higher proportion of decomposers and detritivores, particularly soil saprotrophic fungi, enchytraeids and coprophagous coleoptera (e.g. dung beetles), arbuscular mycorrhizal fungi and root endophytes, and most groups of plant phytoparasites and autotrophs. Grasslands also exhibited a higher proportion of bacterivores and omnivores compared to forests, notably bacterivorous nematodes, rotifers and epigeic collembola, at higher trophic levels. We further examined the ß-diversity of interactions between habitats and compared it to null expectations derived from randomising trophic group abundances while preserving trophic group diversity and network structure (see methods). The observed ßdiversity was significantly higher than expected by chance (Fig. 4), indicating that trophic groups enriched in each habitat were not randomly distributed in the network, but instead formed pairs of trophically linked groups. This suggests that differences in trophic group abundance between habitats (Fig. 3) are influenced by the structure of trophic links within the network (Q2).

The composition of forest and grassland soil multitrophic networks was influenced by distinct ecological filters (Q3). In both habitats, predictors associated with these filters accounted for approximately 10% of the variance in soil multitrophic network dissimilarity (Fig. 5, Fig. S6).

In forests, spatial and environmental predictors played a more important role, particularly soil C/N ratio (related to organic matter degradability), pH, and frozen degree days (FDD) as indicators of abiotic stress (Fig. 5a). Dissimilarities in trophic groups and interactions were more pronounced between sites separated by large spatial distances and remained consistent along the gradients of FDD and pH (Fig. 5b). Changes in soil C/N ratio primarily influenced interaction dissimilarity, particularly at higher values of the C/N ratio gradient (>20).

In contrast, grassland soil multitrophic networks were predominantly driven by plant taxonomic and functional composition, as well as soil C/N ratio (Fig. 5a). Soil C/N ratio and plant taxonomic/functional composition better explained interaction dissimilarity (Fig. 5c), while plant taxonomic composition was more effective in predicting trophic group dissimilarity. Dissimilarities in trophic groups and interactions remained relatively constant along the plant dissimilarity gradient but were more pronounced in the lower part of the FDD gradient (> 40 cumulative degrees) and the upper part of the C/N ratio gradient (>15).

Discussion

Our study employed high-resolution data spanning contrasted habitats and extensive environmental gradients to quantify the variability of soil multitrophic network structure. This approach revealed new insights into how the environment, space and biotic interactions shape soil multitrophic network assembly in the French Alps. We observed significant differences in dominant trophic groups between habitats, that was associated with the structure of the trophic links in the network (i.e. the abundance of trophically linked groups co-vary across habitats).

Furthermore, we uncovered notable discrepancies in the relative importance of the ecological filters between forests and grasslands.

Habitat filtering in soil multitrophic networks

This study encompassed diverse climatic regions in the French Alps (e.g. alpine, continental, mediterranean) and an elevation range of 2500m, resulting in significant environmental variability within and across habitat. As anticipated, habitat filtering strongly influenced soil multitrophic network structure, especially when comparing forests and grasslands [START_REF] Crotty | Divergence of feeding channels within the soil food web determined by ecosystem type[END_REF][START_REF] Arribas | The limited spatial scale of dispersal in soil arthropods revealed with whole-community haplotype-level metabarcoding[END_REF][START_REF] Sepp | Woody encroachment in grassland elicits complex changes in the functional structure of above-and belowground biota[END_REF]. Differences in abiotic conditions and resource availability between habitats could explain the observed differences in the proportions of trophic groups. Grassland soils exhibit lower acidity, higher P and N availability and a lower C/N ratio [START_REF] Joimel | Physico-chemical characteristics of topsoil for contrasted forest, agricultural, urban and industrial land uses in France[END_REF], Fig. S1), and have a higher below-ground plant biomass, which is more easily decomposed, leading to significant inputs of readily decomposable organic matter into the soil [START_REF] Mason | GRASSLAND SOILS[END_REF][START_REF] Heděnec | Global distribution of soil fauna functional groups and their estimated litter consumption across biomes[END_REF]. These factors likely enhanced resource availability for detritus-based trophic groups, resulting in the higher dominance of copiotroph decomposers (e.g. saprophytic fungi and bacteria), micro-or mesodetritivores (i.e. Enchytraeidae) and microbivores in grasslands. In contrast, plant detritus in forests, which is more challenging to decompose, serves as a direct resource for litter and wood saprotrophic fungi and provide a habitat for macro-detritivores, which feed on plant litter and associated fungi [START_REF] David | The ecology of saprophagous macroarthropods (millipedes, woodlice) in the context of global change[END_REF][START_REF] Zuo | Diversity of macro-detritivores in dead wood is influenced by tree species, decay stage and environment[END_REF]. Light availability in grasslands may favour the presence of photoautotrophic bacteria that use light as their primary resource, and facilitate interactions, such as arbuscular mycorrhiza [START_REF] Konvalinková | Lights Off for Arbuscular Mycorrhiza: On Its Symbiotic Functioning under Light Deprivation[END_REF]. Moreover, grasslands had a higher proportion of phytoparasites, likely due to the higher biomass of fine roots in this habitat [START_REF] Jackson | A global budget for fine root biomass, surface area, and nutrient contents[END_REF][START_REF] Heděnec | Global distribution of soil fauna functional groups and their estimated litter consumption across biomes[END_REF], which are more easy to colonise by symbionts, while forests had a higher proportion of zooparasites, benefiting from a greater abundance of potential hosts, predominantly arthropods, molluscs or earthworms, for the taxa belonging to our zooparasite groups. These findings align with our results showing that trophic interactions within the soil multitrophic network contributed to explain network dissimilarities between habitats (Fig. 4) and are consistent with previous studies examining changes in the trophic composition of protists [START_REF] De Araujo | Protist species richness and soil microbiome complexity increase towards climax vegetation in the Brazilian Cerrado[END_REF][START_REF] Fiore-Donno | Contrasting Responses of Protistan Plant Parasites and Phagotrophs to Ecosystems, Land Management and Soil Properties[END_REF] and nematodes [START_REF] Zhao | Soil energy pathways of different ecosystems using nematode trophic group analysis: a meta analysis[END_REF] across habitats. The distinct identities of soil antagonists in the soil multitrophic networks from grassland to forests warrant further investigation as they could provide insights for land management strategies in mountain ecosystems [START_REF] Wall | Soil biodiversity and human health[END_REF].

Impact of ecological filters within habitats

Our study design aimed to investigate the impact of ecological processes on soil multitrophic network structure across elevation gradients and diverse environmental conditions. It is important to note that the studied mountain system included various forest types (e.g. lowelevation deciduous, high-elevation coniferous) and grassland types (e.g. subalpine, alpine).

To capture this ecological heterogeneity within habitats, our research incorporated variations in abiotic conditions, plant composition, and geographic distance, offering a comprehensive analysis of the complex environmental variation characterising our mountainous systems.

Although the GDM explained a relatively low percentage of variance when quantifying changes in the structure of local networks (Fig. 5a), it is important to note that typical deviance explained for models of compositional dissimilarity ranges from 20% to 50%, and this percentage can decrease with an increasing number of sites, which was relatively high in our study [START_REF] Ferrier | Using generalized dissimilarity modelling to analyse and predict patterns of beta diversity in regional biodiversity assessment[END_REF][START_REF] Mokany | A working guide to harnessing generalized dissimilarity modelling for biodiversity analysis and conservation assessment[END_REF]. Furthermore, the resolution of our soil multitrophic networks using trophic groups instead of species led to low average dissimilarities (0.27 and 0.48 in forests, and 0.30 and 0.53 in grasslands for groups and interactions respectively), reducing model predictability [START_REF] Mokany | A working guide to harnessing generalized dissimilarity modelling for biodiversity analysis and conservation assessment[END_REF]. Importantly, we found that trophic link structure significantly influenced changes in trophic groups abundance composition (Fig. 4). However, current models cannot account for this factor, and incorporating network structure in diversity models may enhance predictability but requires further implementation [START_REF] Poggiato | Integrating food webs in species distribution models improves ecological niche estimation and predictions[END_REF]. Spatial processes, environmental filtering, and changes in plant community composition explained variation in soil multitrophic networks within habitats. In forests, spatial distance played a significant role, explaining a substantial portion of trophic group (31%) and interaction dissimilarity (58%), even when accounting for confounding environmental changes.

Limited dispersal opportunities in closed forest habitats, particularly for organisms dispersed through wind (anemochory), might influence soil multitrophic network structure in mountain forests. Conversely, open grassland habitats subject to strong winds [START_REF] Tackenberg | Wind dispersal of alpine plant species: A comparison with lowland species[END_REF] and high levels of herbivory, both domestic and wild, can facilitate the dispersal of soil organisms. Furthermore, mountain massifs and large distances between forested areas (up to 250 km) can impact meso-and macrofauna, susceptible to ecological drift due to their limited dispersal capacity or smaller population sizes [START_REF] Arribas | The limited spatial scale of dispersal in soil arthropods revealed with whole-community haplotype-level metabarcoding[END_REF][START_REF] Kang | Stochastic processes regulate belowground community assembly in alpine grasslands on the Tibetan Plateau[END_REF]. This is in line with our results showing that network dissimilarity increased with greater spatial distances (Fig. 5b). Yet, the spatial distance used in the model may not fully capture the complexity of the spatial structure, and further investigation using experimental setups to clarify spatial connectivity and dispersal limitation effects is required.

Environmental filtering played a significant role in both habitats, primarily influenced by soil edaphic properties. In the French Alps, the soil C/N ratio emerged as the most influential abiotic factor in shaping soil multitrophic network composition. It reflects organic matter decomposability and resource availability, as discussed above, impacting decomposer groups, their consumers, and the entire feeding channel. Conversely, variables related to total energy or resources, such as the NDVI, GDD and SOM, had a smaller impact on soil multitrophic network dissimilarity. While this contrasts with previous studies emphasizing the role of these variables in diversity patterns of soil groups in mountains [START_REF] Zinger | Contrasting Diversity Patterns of Crenarchaeal, Bacterial and Fungal Soil Communities in an Alpine Landscape[END_REF][START_REF] Calderón-Sanou | Energy and physiological tolerance explain multi-trophic soil diversity in temperate mountains[END_REF][START_REF] Steinwandter | Ground-dwelling invertebrates of the high alpine: Changes in diversity and community composition along elevation (1500-3000 m)[END_REF], it highlights the emergence of different ecological drivers at different levels of biodiversity organization. In our study, resource quality, particularly litter quality, seemed to dominate network composition rather than quantity.

Additionally, environmental factors associated with harsh conditions, such as pH and frozen degree days, known drivers of soil taxonomic diversity and composition [START_REF] Fierer | The diversity and biogeography of soil bacterial communities[END_REF][START_REF] Decaëns | Macroecological patterns in soil communities[END_REF][START_REF] Gobbi | Hay meadows' overriding effect shapes ground beetle functional diversity in mountainous landscapes[END_REF], also played a major role in shaping soil multitrophic networks.

Plant-soil interactions likely play a crucial role in shaping grasslands soil multitrophic networks in the French Alps. Although we did not directly measure these interactions, the taxonomic and/or functional composition of plant communities had a major impact on soil multitrophic network structure. Mountain grasslands exhibit high variability in plant species richness and the presence of functionally important taxa such as Fabaceae, Brassicaceae or Poaceae across different elevations and massifs. Changes in plant diversity can influence the composition of soil multitrophic networks [START_REF] Eisenhauer | Plant diversity effects on soil food webs are stronger than those of elevated CO2 and N deposition in a long-term grassland experiment[END_REF] , and altering plant identity can affect different soil taxa [START_REF] Zinger | Contrasting Diversity Patterns of Crenarchaeal, Bacterial and Fungal Soil Communities in an Alpine Landscape[END_REF]. Interestingly, taxonomic dissimilarity had a larger effect than plant functional dissimilarity, possibly because above-ground traits used in our study may have a lesser role in plant-soil interactions compared to root traits (Kardol and De Long, 2018;[START_REF] Wilschut | Root traits and belowground herbivores relate to plant-soil feedback variation among congeners[END_REF]. Surprisingly, in forests, plant composition had a minor impact on soil multitrophic network structure. Soil variables like pH and C/N ratio may better capture variations in forest type (e.g. broadleaf vs coniferous forests), directly influencing soil multitrophic networks.

Conclusions

This study enhances our understanding of how ecological processes operate in below-ground communities and networks, thereby extending our knowledge of ecological communities as a whole. Furthermore, our research offers valuable insights into the impacts of global change on soil biodiversity and terrestrial ecosystem functioning, particularly in vulnerable mountainous systems experiencing climate change and land-use modifications. Land management decisions, such as forestry, agriculture, and tourism, can alter local land-use and habitats, necessitating consideration of their effects on soil multitrophic networks and ecosystem functioning.

However, further investigations are needed to explore the underlying mechanisms driving spatial changes in network structure. Additionally, studying the recovery and resilience of soil multitrophic networks to land-use changes is essential, considering the dynamic nature of ecosystems in the face of a constantly changing climate. These future studies will enhance our understanding of soil multitrophic networks and their responses to environmental changes, aiding informed land management decisions and sustainable practices.

Figure 1 .

 1 Figure 1. Study design to investigate soil multitrophic network structure in the French Alps (a) Study sites across the French Alps were selected using a stratified sampling design across elevations and nested design with subplots within gradients. The collected 418 soil samples were distributed across nested spatial distances ranging from 8 m to 250 km. (b) eDNA and bioinformatic analysis were used to assign detected organisms to taxonomic and functional information. (c) The metaweb framework was used to quantify spatial changes in the structure of soil multitrophic networks across the study area and examine their associations with environmental and spatial factors.

Figure 2 .

 2 Figure 2. 2-D plane representing the dissimilarity between soil multitrophic networks based on the composition of trophic groups and trophic interactions with a colour scale representing the habitat type. Each dot is a local network from a single sub-plot. Dots that are close to each other in the 2D plane have similar trophic group compositions or interactions.

Figure 3 .

 3 Figure 3. Differences in trophic group and class abundances between soil multitrophic networks in forests and grasslands. Colours indicate the significance of Wilcoxon tests calculated for the difference in abundances of each group between habitats (p<0.05) while node sizes are proportional to significant abundance differences. A Wilcoxon test was also applied to compare differences in trophic class abundances (dashed circles) between habitats. Abundance of trophic classes correspond to the sum of trophic groups relative abundances.

Figure 4 .

 4 Figure 4. Null distribution of the ß-diversity of interactions (ßTL) at fixed trophic group diversity.The values of the histogram show the ßTL between the average soil multitrophic networks of forest and grassland computed using 3000 permutations of trophic group abundances at fixed network structure, leading to various ßTL values while keeping trophic group ß-diversity constant. The vertical line shows the observed ßTL across habitats. The value of ßTL varies from one when the average network of forest and grassland have identical interaction abundances distribution, to two, when the networks do not share any common group or interaction.

Figure 5 .

 5 Figure 5. Main predictors of soil multitrophic network dissimilarity within forest and grassland habitats. The relative importance of each predictor was estimated using Generalised Dissimilarity Models (a). It is measured as the sum of the coefficients of the three I-splines of the focus predictor and represents the total amount of change along the predictor gradient. Partial curves for three of the most important predictors of network dissimilarity in forest (b) and grassland (c). Each panel shows the network dissimilarity as a function of an environmental predictor when holding all other variables to their mean. The slope at any point on the curve indicates the rate of network dissimilarity at that position along the environmental gradient (x-axis), while the total height reached by the function indicates the total amount of network dissimilarity due to that environmental predictor. SOM: Soil Organic Matter, FDD: Frost Degree Days, GDD: Growing Degree Days, AP: Annual Precipitation, NDVI: Normalised Difference Vegetation Index.

Table 1 .

 1 Environmental variables used in this study.

	524
		Environmental variable	Abbr. Description	Temporal range	Resolution Source
		Growing degree days	GDD	Annual sum of degree-days over a 0°C soil temperature
	Climate	Annual precipitation	AP	The sum of all the monthly precipitation estimates Annual sum of average daily	Average across	Plot	SAFRAN-SURFEX/ISBA-Crocus-MEPRA reanalysis
				degrees below zero in soil	2009-2019		(Vernay et al., 2021)
		Frost degree days	FDD	temperature. Summarizes the
				duration and intensity of ground
				freezing events.
		Soil pH	pH	Top-soil wahter pH	Same year	Sample
		Soil Organic Matter	SOM	Soil organic matter content	than the	Sample
	Soil	Soil C/N	C/N	Soil carbon to nitrogen ration in organic matter	sample was collected (2016-	Sample	Measured in the lab
					2020)
	Vegetation	Plant taxonomic dissimilarity Plant functional dissimilarity		Jaccard pairwise dissimilarity index between subplots. Jaccard-like functional index. Traits: specific leaf area (SLA), leaf C/N, root depth (ordered variable), vegetative plant height, and woodyness index (categorical variable ranging from herbaceous to woody).	Same year than the sample was collected (2016-2020)	Plot Plot	Botanical surveys Botanical surveys and our own trait measurement values for species (median values across individuals). Root depth was , extracted from (Landolt et al., 2010).
				Estimated from the surface spectral			MODIS (Moderate
		Normalised Difference Vegetation Index	NDVI	reflectance at a resolution of 250 m. Raw NDVI times series were pre-processed following (Choler, 2015), and we kept the mean yearly sum of	Average across 2010-2020	Plot	Resolution Imaging Spectroradiometer), available online: https://lpdaac.usgs.gov/produ
				NDVI greater than 0.2			cts/mod09q1v006/.

This project was possible thanks to the many people that participated in the field work and in the lab procedures, including students, interns, postdoctoral fellows, field assistants, park managers and rangers. We thank Cindy Arnoldi for conducting the soil physicochemical analyses.

This research was funded by the French Agence Nationale de la Recherche (ANR) through the GlobNet & EcoNet projects (ANR-16-CE02-0009, ANR-18-CE02-0010-01), from two 'Investissement d'Avenir' grant (Idex UGA: ANR-15-IDEX-02; Montane: OSUG@2020: ANR-10-LAB-56, MIAI@Grenoble Alpes: ANR-19-P3IA-0003), the METRO Grenoble Alpes, the Isere departement and the Office Français de la Biodiversité (OFB). Support was also received from the Lautaret Garden-UMS 3370 (Univ. Grenoble Alpes, CNRS, SAJF, 38000 Grenoble, France), member of AnaEE-France (ANR-11-INBS-0001AnaEE-Services, Investissements d'Avenir frame) and of the LTSER Zone Atelier Alpes, a member of the eLTER-Europe network.


Alice Berthelot
email: alice.berthelot@univ-lyon1.fr

Elodie Romeo

Xavier Dagany

Thierry Deschamps

Eliot Herry

Robin Kineider

Yannis De Leon

Emmeline Brassac

Guillaume Kermouche

Etienne Barthel

Christine Martinet
email: christine.martinet@univ-lyon1.fr

Strength of Brillouin spectroscopy to identify spatial densification model in indented silica

Keywords: mechanical of glass, vibrationnal spectroscopies, plasticity, indentation

Introduction

States of Art

Due to its amorphous nature, mechanical transformations of glass, especially plastic deformation, are very difficult to model. The need to predict the response of glass answer to external strains implies the development of closely related experiments and simulations. One model experiment is the study of glasses under micro-indentation, in particular the analysis of the imprint and the cracks formed during the process, which is strongly dependent on the nature of the glass as described by Rouxel et al [START_REF] Rouxel | Indentation of glasses[END_REF]. A complementary experiment is the application of pure hydrostatic pressure into the glass sample, where plastic deformation of silicate glasses occurs above the elastic limit, between 6 and 9 GPa, depending on the silicate glass composition [START_REF] Hemley | Raman spectroscopy of siO 2 glass at high pressure[END_REF][START_REF] Huang | Amorphous-amorphous transitions in silica glass. II. Irreversible transitions and densification limit[END_REF][START_REF] Champagnon | High pressure elastic and plastic deformations of silica : In situ diamond anvil cell raman experiments[END_REF][START_REF] Benmore | Structural and topological changes in silica glass at pressure[END_REF][START_REF] Deschamps | Soda-lime silicate glass under hydrostatic pressure and indentation : A micro-raman study[END_REF][START_REF] Tran | Micro-Brillouin spectroscopy mapping of the residual density field induced by Vickers indentation in a soda-lime silicate glass[END_REF]. For hydrostatic and non-hydrostatic stresses, densification rate depends on the initial structure but also on the thermo-mechanical history and the chemical composition of the glass [START_REF] Barlet | Hardness and toughness of sodium borosilicate glasses via vickers's indentations[END_REF][START_REF] Limbach | Elasticity, deformation and fracture of mixed fluoride phosphate glasses[END_REF]. In particular, non-isotropic stress can induce another plastic deformation process in competition with densification : shear flow [START_REF] Rouxel | Poisson's ratio and the densification of glass under high pressure[END_REF][START_REF] Kermouche | Mechanical modelling of indentationinduced densification in amorphous silica[END_REF][START_REF] Kermouche | Perfectly plastic flow in silica glass[END_REF].

Models attempt to determine the effects of each process on the evolution of the glass structure and the emergence of cracks.

In silicate glasses, plastic deformation usually results in a permanent increase in density [START_REF] Ernsberger | Role of densification in deformation of glasses under point loading[END_REF] which has been well correlated with Poisson's ratio by AFM measurements [START_REF] Yoshida | Indentation induced densification of sodium borate glasses[END_REF]. As an anomalous glass with low Poisson's ratio, fused silica, is more sensitive to densification process than to shear flow. This material has a high capacity for permanent densification up to a density increase of 21 % [START_REF] Sonneville | Progressive transformations of silica glass upon densification[END_REF][START_REF] Deschamps | Elastic anomalous behavior of silica glass under high-pressure : In-situ raman study[END_REF]. Soda-lime silicates, with a higher Poisson's ratio and a lower ability to densify compared to silica, are considered as normal glasses and show large plastic shear flows [START_REF] Rouxel | Indentation deformation mechanism in glass : Densification versus shear flow[END_REF]. The plastic behavior of silicate glasses can be described by a hybrid law involving pressure and shear stress [START_REF] Kermouche | Mechanical modelling of indentationinduced densification in amorphous silica[END_REF][START_REF] Lambropoulos | Surface microroughness of optical glasses under deterministic microgrinding[END_REF].

Classical Raman Spectrometry for indents study

In order to validate these indentation models, micro-spectroscopic studies were carried out on indent. To this end, the first measurements have allowed to highlight the normal/anomalous characteristics of different silicate glasses [START_REF] Arora | Indentation deformation/fracture of normal and anomalous glasses[END_REF]. Micro-Raman spectroscopy studies were then carried out to probe the densification inside the imprint. For the float glass with higher Poisson's ratio, the Raman mapping of the barycentre of the Raman bands 500-700 cm -1 , gave a density cartography [START_REF] Deschamps | Soda-lime silicate glass under hydrostatic pressure and indentation : A micro-raman study[END_REF]. Deschamps et al observed a progressive evolution of densification as a function of the distance from the centre of the imprint. The fracture lines that separate the different iso-density regions are well reproduced by finite element simulations [START_REF] Kassir-Bodon | Raman mapping of the indentation-induced densification of a soda-lime-silicate glass[END_REF].

Then, for silica with a low Poisson's ratio, such as pure silica, cartographies of the Raman shift of the D2 band, which corresponds to the breathing mode of 3-fold rings, gave an indentation-induced densification map, which allowed the mechanical behaviour of glasse under indentation to be implemented [START_REF] Perriot | Raman microspectroscopic characterization of amorphous silica plastic behavior[END_REF]. For this last experiment, a constitutive law derived from finite element simulations reproduced quite well the experimental results of the indentation-induced densification rate map [START_REF] Kermouche | Mechanical modelling of indentationinduced densification in amorphous silica[END_REF] by also taking into account the hydrostatic hardening behaviour [START_REF] Vandembroucq | Density hardening plasticity and mechanical ageing of silica glass under pressure : A raman spectroscopic study[END_REF]. However, we have recently shown that the D2 band is not only sensitive to densification but also to shear stress and is not a reliable marker of the densification process [START_REF] Martinet | Highlighting the impact of shear strain on the SiO2 glass structure : From experiments to atomistic simulations[END_REF].

Chemical dissolution experiments on silicate imprints after 10mN and 50mN indentations for silica and soda-lime glasses, have allowed the observation of a quasi-homogeneously densified region just below the surface in the indented region [START_REF] Niu | Locally enhanced dissolution rate as a probe for nanocontact-induced densification in oxide glasses[END_REF]. Very recent work of Bruns et al presents results from a classical Raman study of indented pristine silica with loads varying from 0.3 N to 10 N [START_REF] Bruns | Indentation densification of fused silica assessed by raman spectroscopy and constitutive finite element analysis[END_REF]. They study the densification process using Raman spectroscopy of the main band (250 cm -1 -550 cm -1 ) which provides an average analysis of the Si-O-Si angle distribution in glasses and has so far been a good probe of the densification process for hydrostatic compression [START_REF] Deschamps | Permanent densification of compressed silica glass : A raman-density calibration curve[END_REF]. By measuring the evolution of the main band centroid with the depth of laser focus, Bruns et al have compared their experimental results with a constitutive finite element analysis. They found that the saturation of densification under indentation is closed to 18 % for different size of indent.

Nevertheless, this classical Raman analysis spectrum is problematic if we want to study the spatial evolution of densification. The shape of the laser beam used for Raman spectrometry does not allow to remove the answer of the pristine silica that is located below the indented area. In fact, since the main bands of the Raman spectra of the two glasses (densified and non-densified) overlap strongly, it is difficult to separate the imprint contribution from the far field non-densified response.

Brillouin Spectrometry for indents study

As an alternative technique to classical Raman spectroscopy, Brillouin scattering spectroscopy can measure the energy of acoustic phonons of transparent materials and provides elastic constants if density, longitudinal and transverse acoustic modes have been measured. The evolution of these coupled physical quantities can then be studied by measuring the Brillouin shift. In the back-scattering geometry, only the longitudinal acoustic mode is measured and Brillouin shift ν B can be described by the following formula :

ν B = 2nV AL λ
where n is the optical refactive index, λ the wavelength of the laser,V AL the longitudinal acoustic wave velocity. V AL depends on Young's modulus, density, Poisson's ratio and the optical refractive index which is difficult to measure. A calibration curve relating the Brillouin shift to the densification factor was then deduced experimentally from macroscopic densified silica samples where the density has been measured by the Archimedian method [START_REF] Deschamps | Elastic moduli of permanently densified silica glasses[END_REF]. Using Brillouin spectrometry, Tran et al studied the densification of the imprint of float glass under 20 N load indentation [START_REF] Tran | Micro-Brillouin spectroscopy mapping of the residual density field induced by Vickers indentation in a soda-lime silicate glass[END_REF]. Although the Brillouin shift between the indent center and pristine float glass is small, this work has shown that the evolution of the densification rate is progressive within the imprint. Despite the fact that these calibration curves do not take into account the effects of shear stresses, which are far from negligible in an indentation, this technique is very promising to complement other spectroscopic techniques, as Raman and luminescence spectroscopies [START_REF] Perriot | Raman microspectroscopic characterization of amorphous silica plastic behavior[END_REF] and will allow to roughly estimate the density evolution of glass inside the imprints for silicate glasses. Then, our present work proposes to compare the well-known behaviour of a normal silicate glass, float glass, with an anomalous glass, silica, under indentation and to analyse the indented volume by micro-Brillouin spectrometry.

Experiments

Samples description

Samples used in this study correspond to pure SiO 2 glass (commercial "Suprasil 300"), with a low OH content, less than 1 ppm, and a commercial float glass (mainly chemical components : 72% SiO 2 , 15% N a 2 O, 8% CaO, 4% M gO and about 1% (Al 2 O 3 + F eO 3 + SO 3 ) in mole %). On flat and polished face of these two glasses, Vickers micro-indentations of different loads have been performed from 1 N to 10 N. 1 N load imprints does not exhibit any microscopic fracture initiation whereas 3 N and 10 N load imprints show microscopic fractures below the indent and laterally close to the surface of the glass.

Brillouin spectrometric measurements

We performed micro-Brillouin experiments at room temperature, using a Sandercock tandem Fabry Perot [START_REF] Lindsay | Construction and alignment of a high performance multipass vernier tandem Fabry-Perot interferometer[END_REF], coupled with a microscope in back-scattering geometry with a x100 long working distance objective (diameter of the laser spot and its Rayleigh length are about 0.520 µm and 0.410 µm respectively), using a 532 nm N d 3+ :YAG laser excitation. The micro-positioning is ensured by a Marzhauser X-Y-Z motorized micro-stage with a spatial minimal step of 1 µm. Each spectrum have been recorded after a long accumulation time (about 90 minutes by spectrum), in order to obtain a good signal to noise ratio. For each data spectrum, band, this band widens to reveal a broader band at higher frequency corresponding to the contribution of the plastically deformed region under the indent. These broadening and shift of the Brillouin line are in agreement with the previous work of Niu et al [START_REF] Tran | Micro-Brillouin spectroscopy mapping of the residual density field induced by Vickers indentation in a soda-lime silicate glass[END_REF]. On the contrary, for silica, two distinct bands appear at the center of the indent, one corresponding to undensified silica glass and one corresponding to the contribution of highly densified glass. This last band can be attributed to the glass just below the imprint. These two distinct bands have never been seen before and needs to be understood.

Results

Vertical dependance of densification for Silica and window glass imprints

The discontinuity of the Brillouin spectra, observed for silica glass and not for float glass, has to examined in the light of the nature of the optical probe used. This geometry allows us to have a high transversal resolution (in X and Y directions), but a poor vertical resolution because the laser beam is focused over length of few microns. When the laser focus is in the indented region, we observe the Brillouin spectra of the densified glasses, but also the Brillouin signal scattered by the lower layers. This second contribution is reduced due to the confocal geometry, but is still been observed [START_REF] Bruns | Indentation densification of fused silica assessed by raman spectroscopy and constitutive finite element analysis[END_REF]. In other words, the observed spectrum is mainly due to scattering within the Rayleigh region, but also from the region below and above the focused one. In order to separate the different possible contributions at different depths inside the sample, we have performed X-Y and Z scans of Brillouin measurements inside the indent. Z-Scan at the centre of the indent were obtained for both glasses by progressively varying the depth of the focus point. In order to have access to the response of the surface we performed measurements with a laser focalisation above 5 µm from the surface to 28 µm for the depth. The evolution of the normalized spectrum with depth is plotted in Figure 1.c for float glass and in Figure .1.d for silica glass. Even tough the obtained mapsare raw results, these graphs allow us to access to the depth of the indented area and its location, and to study the evolution of densification as a function of depth of focus. The results are then compared with existing models.

Evolution of Brillouin shift with radial distance for silica glass 10 N imprint.

Figure 2 shows the evolution of the Brillouin spectra with lateral positions from the centre of the indent to a distance of 15 µm along a line passing from the indent centre to a pyramid face centre at a fixed depth (depth of the imprint surface centre). We observe a large change in the Brillouin spectrum as we move away from the centre, the densified peak shifts continuously from 38 GHz to the undensified peak. The relative area of this peak is almost constant for distances less than 7.5 µm. For measurement points at distance greater than 10 µm, we observe in the Brillouin spectrum an additional component in the low energy part for the undensified peak, wich can be attributed to the scattering by glass cracks.

In order to compare the depth and lateral evolutions of the compression processes, we have plotted in respectively. We observe a very weak depth evolution of the densified peak (peak D), which never reaches the the undensified peak (peak U). the densification process is not progressive along the depth, unlike the evolution along the radial direction. Then, in a first approach, we can question any hemispherical symmetry for the densification process.

Evolution of Brillouin shift at the imprint center with loads.

Figure 4 shows the evolution the Brillouin spectra after for different loads (from 1 N to 10 N). spectra were measured with a laser focus at the surface of the indent centre. For comparison, the Brillouin spectrum of undensified silica has been added to the figure. The two peaks described above are observed for all loads. It appears that the position and intensity of the peak D depends on the load of the indenter. Figure 5 presents the evolution of the relative weight of peak D, as a function of the depth of focus, for the different indenter loads. The relative weight of peak D is calculated by dividing the area of peak D by the sum of the areas of each peak. These plots give access to the location of the modified area, under the imprint and its evolution with load. As expected, we observed an increase in the densified area with load indenter. Although this study is hampered by the Gaussian structure of the laser beam which is focused oover a length equal to twice the Rayleigh zone, in our case just over 1 µm, it is possible to study the evolution of the densification location at the µm scale. The evolution of the Brillouin shift as a function of the depth of focus is different from 1 N to 10 N mass load. In fact, for the 1 N curve, the densified peak (peak D) is significant only at the surface of the sample (for a focus between -5 to 5 µm) and disappears just below it. Peak D is large and centred around 36.2 GHz, which means that we are already observing densified silica at this loading. From the results of Deschamps results [START_REF] Deschamps | Elastic moduli of permanently densified silica glasses[END_REF], with the maximum of peak D, we can roughly estimate a densification ratio of 14% for 1 N sample for a depth equal to 0 µm. For the 3 N indented sample, the maximum of the Brillouin shift is about 36.8 GHz, corresponding to a densification of about 15 %, and is localised at the surface of the sample. For this sample, the relative area of the peak D is less than 10% for a depth of focus of 10 µm. For the 10 N sample, the densification appears to be higher, 17%, than for the 3 N sample and on a larger volume. Since the strain distribution is expected to be invariant due to the equivalent cone shape of the indenter, this variation of the maximum density with load may illustrate the difficulties of the spectroscopic mapping at small sizes.

Discussion

Brillouin spectrometry capacity to differentiate the undensified and densified glass answers

For float glass, as a normal glass, the densification process appears to be progressive in all with a constant decreasing in densification as one moves away from the area [START_REF] Tran | Micro-Brillouin spectroscopy mapping of the residual density field induced by Vickers indentation in a soda-lime silicate glass[END_REF]. This progressive density change has already been confirmed by Raman experiments and FE simulations [START_REF] Kassir-Bodon | Raman mapping of the indentation-induced densification of a soda-lime-silicate glass[END_REF]. Brillouin spectrometry shows a broadening a slight shift of the line towards higher Nevertheless, the Brillouin shift of the longitudinal mode have a small amplitude variation from 33.5 GHz to 34 GHz between pristine and fully densified float glass, the limit of the resolution of the spectrometer, to 0. 5 GHz. Our spectrometer does not us to resolve a possible substructure of the Brillouin spectra. This variation is too small to conclude on the presence or the absence of a sublayer with uniform densification. Nevertheless, the well-documented theoretical and experimental results show that this progressive density change the imprint centre perfectly matches the cracks pattern of float glass.

In the literature, the densification ratio of silica can be determined from Raman measurements using the position of the barycentre, called the centroid, of the main band [START_REF] Deschamps | Permanent densification of compressed silica glass : A raman-density calibration curve[END_REF]. Studying of the z-dependence of densification after indentation and using main band Raman spectrometry become an hazardous challenge due to the large spectral overlap of densified and undensified components [START_REF] Bruns | Indentation densification of fused silica assessed by raman spectroscopy and constitutive finite element analysis[END_REF]. This measure can, at best, propose a mean spectrum between densified region where the laser is focused and the far-field response of the undensified glass as suggested by Bruns. Our results show that Brillouin spectrometry for silica indentation has the advantage of spectrally separating the two components and allowing us to follow the evolution of densification with z. For silica, as an anomalous glass, the amplitude of densification is large enough to observe modifications of the Brillouin peak position and structure, which varies from 31.5 GHz to 38 GHz.

Our Brillouin spectrometer is able to distinguish the substructure of the Brillouin spectrum produced by different densified regions (Fig. 1).

As we have seen few intermediate components between the densified and undensified peaks, we can conclude that the densification does not decrease continuously from the surface of the indent. However, in Fig. 3a we see a densified peak with a FWHM around 2 GHz, which is broader than the undensified peak.

So we have an inhomogeneous broadening of the Brillouin line due to a partial gradient of densification within the indented area. We do not observe a fully densified silica, but an intermediate densified silica with a densification ratio varying between 14 and 16.6 %. This result suggests that the permanent densification under the indentation does not decrease smoothly with depth from the indentation point surface as it has been observed on float glass and previously in the literature on silica [START_REF] Perriot | Raman microspectroscopic characterization of amorphous silica plastic behavior[END_REF][START_REF] Bruns | Indentation densification of fused silica assessed by raman spectroscopy and constitutive finite element analysis[END_REF]. This vertical discontinuity in the densification process could be understood by considering the difference in compressibility and availability of each material to shear flow : silica, as anomalous glass, is a highly compressible glass with a maximum density change of about 21% and float glass, as a normal glass, has a lower maximum density change of about 6%. As Januchta et al explains that, depending on densification ability and Poisson's ratio, two mechanical behaviours occur under indentation [START_REF] Januchta | Indentation deformation in oxide glasses : Quantification, structural changes, and relation to cracking[END_REF]. For sufficiently high Poisson's ratios, for example float glass with a Poisson's ratio of 0.23, which ensures a large shear flow, we observe a progressive decrease of densification in both the lateral and longitudinal directions of the imprint [START_REF] Tran | Micro-Brillouin spectroscopy mapping of the residual density field induced by Vickers indentation in a soda-lime silicate glass[END_REF][START_REF] Kassir-Bodon | Raman mapping of the indentation-induced densification of a soda-lime-silicate glass[END_REF]. For low Poisson's ratio materials, e.g. silica glass with a Poisson's ratio equal of 0.15, the density pattern under the imprint is not clearly defined. Our results do not agree with the experimental Raman results obtained by [START_REF] Perriot | Raman microspectroscopic characterization of amorphous silica plastic behavior[END_REF] and shown in Fig. 6a, where the densification is progressive along the vertical and lateral directions. In this case, these results were obtained from Raman scattering, using the evolution of the D2 shift as a function of density. However, we have recently demonstrated that the D2 line (position and intensity) is highly dependent on the thermomechanical history the densified silica [START_REF] Martinet | Permanently densified SiO2 glasses : A structural approach[END_REF]. Then, the reliability of a clear relations between the D2 line and density is not so obvious. Nevertheless, our present results are consistent with the FE simulation of a 2D-axisymmetric cone, where a quasi-homogenous densified region was observed, followed by a rapid decrease in densification below [START_REF] Perriot | Raman microspectroscopic characterization of amorphous silica plastic behavior[END_REF]. Rouxel et al has hyptothesised the creation of a densified region below the indented surface with a sharp transition to the far-field undensified region [START_REF] Rouxel | Indentation deformation mechanism in glass : Densification versus shear flow[END_REF]. To our knowledge, only the experimental work of Niu et al, shown in Fig. 6b), has demonstrated such behaviour for silica glass but also for soda-lime glass [START_REF] Niu | Locally enhanced dissolution rate as a probe for nanocontact-induced densification in oxide glasses[END_REF]. Our present results, resumed in scheme fig. 6c), confirm and refine this "discontinuous" model along the vertical direction for silica glass : we clearly demonstrate a small thickness, dz for the transition region from densified, extending from the surface to z imprint , to undensified glass. Niu's description and our work diverge in the lateral direction : unlike their work, we observe a progressive decrease in densification from the centre of the imprint. This last result seems to agree at least qualitatively with the predictions of Kermouche et al [START_REF] Kermouche | Mechanical modelling of indentationinduced densification in amorphous silica[END_REF]. As a rapid decrease in density between the densified area and undensified regions induces high local stresses, it could explain the origin of the large lateral crack, under the plastic core of the indent, currently observed for silica [START_REF] Arora | Indentation deformation/fracture of normal and anomalous glasses[END_REF].

Test of density evolution

In order to compare our experimental results with indentation simulations, we have calculated Brillouin spectra scattered by two model imprints : one in which the densification ratio evolves linearly from its saturation value (21%), the other in which the densification is constant up to a certain depth and then decreases rapidly and linearly with depth and then falls to a value close to 0. In a first approximation, we have also assumed invariance of the the indent characteristics in the planar directions (X and Y). Since we have measured the evolution of the experimental Brillouin spectrum focus depth, we have calculated and compared Brillouin spectra emitted by the sample for different focus depths. For the first densification model, the slope is the free parameter, for the second model, the densification ratio of the homogeneous region, the beginning and the end of the slope are free parameters.

The calculated Brillouin spectra result from the laser scattered at all depths of the sample. We assume that the number of Brillouin photons is proportional to the incident laser intensity which is calculated for each glass slice, at a distance z from the laser focus point. We have assumed that each glass slice, with a densification rate given by the densification evolution curve to be tested, produces a Brillouin spectrum with a Lorentzian profile L(ν -ν B (z)) with a central position ν B (z) derived from the paper of Deschamps et al densification factor-Brillouin shift calibration curve [START_REF] Deschamps | Elastic moduli of permanently densified silica glasses[END_REF]. The FWHM of the Lorentzian profile was determined experimentally using a homogeneous glass plate (0.5 GHz). Before integrating the response of each layer over the entire thickness of the sample, it is necessary to consider and estimate the confocal filtering of the light for each slide of the sample.

The estimation of the confocal filtering is based on the work of Grauw et al [START_REF] Grauw | Axial resolution of confocal raman microscopes : Gaussian beam theory and practice[END_REF] and is structured around 2 assumptions : the Gaussian structure of the scattered beam and the confocal detection. For each sample slice, located at a distance z from the focus point, we have calculated the light scattered by the slice and transmitted through the pinhole with a diameter R ph . We then obtained a confocal filtering factor, f R ph (z)

for each z. This factor represents the ratio between the intensity scattered by the slice collected with the pinhole and that collected without.

We have defined r, the radius of the incident laser beam at a distance z from the laser focus point. This radius is given by

r(z) = r 0 1 + (z/Z r ) 2
where r 0 is the radius of the laser beam point at the focus point and Z r the Rayleigh length of the Gaussian incident laser beam.

We have defined r as the radius of the image with the detection system of the irradiated area in the sample located at z , z as the image shift offset from the pinhole position. The expressions of r and z were calculated by conjugation relation for microscope from r and z. We have assumed that the structure of the scattered light collected by the miscroscope has a gaussian beam structure with Z rimage , the Rayleigh length associated with the image gaussian beam focused at z . To estimate the confocal filtering factor f R ph (z), we had to calculate r at ph the radius of the image Gaussian beam calculated at the pinhole position :

r at ph (z) = r 1 + (z)/Z rimage ) 2
Finally, rhz expression for f R ph (z) for a slice at z from the laser focus point is :

f R ph (z) = 1 -exp( -2R 2 ph r 2 at ph (z)
)

Once we have obtained the evolution of f R ph (z) with depth shift z, the total scattered Brillouin spectrum generated by sample, is obtained by integration :

S Brillouin (ν) = +∞ -∞ L(ν -ν B (z)) • f R ph (z) • dz
From these expressions, we can calculate the evolution of the Brillouin spectra as a function of the laser focus depth for the two models : the "continuous" one and the "step" one and compare them with the experimental spectra.

By optimising the overlap between the experimental spectra and the modelled spectra, we have attempted to determine the best densification profile that verifies the experimental values for continuous densification ratio and "step" densification ratio evolutions. The results of these optimisations are plotted in is much better at reproducing the evolution of the Brillouin Spectra than the continuous one. We can also see that our model does not reproduce well the linewidth of the densified peak D, which can be a signature of non-uniform densification in the planar direction and a lack of lateral resolution even with our x100 microscope objective, as seen in Fig. 3d and by Bruns et al [START_REF] Bruns | Indentation densification of fused silica assessed by raman spectroscopy and constitutive finite element analysis[END_REF]. The best values for this step model are a maximum densification ratio equal to 16.2%, a thickness z imprint for the densified zone equal to 7.3 µm and a thickness dz for the transition part equal to 12.3 µm. The uncertainty in determining the thickness of the transition is about 1.5 µm. In a reasonable way, we can conclude that the progressive model in z profile is not valid, because it cannot in any way reproduce the two peaks profile and its evolution with depth that we observe experimentally. Our new experimental results therefore provide a pertinent test to challenge theoretical and simulated models of densification within the indented region.

Conclusion

Density mapping from micro-Brillouin scattering micro-indentation imprints had revealed large differences between the Brillouin spectra of silica and float glass. Concerning float glass, our present results confirm a progressive density evolution but do not allow us to conclude on the presence or not of a highly densified core inside the imprint. On the contrary, the depth profile of the imprint on silica shows a quasihomogeneous densification region whose size depends on the load, followed by a smaller transition region to undensified glass. Our experimental results have allowed us to observe and to confirm the localisation of this boundary region below the highly densified glass.A more progressive evolution of density is observed in the direction. To date, and to our knowledge, no model predicts this profile of densification. Our For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

Figure 1 :

 1 Figure 1: Comparison of Brillouin spectra of pristine glass and in the centre of indent of 10 N load for float glass (a) and pure silica (b) for focalisation depth equal to 0 µm (top views). Evolution of area intensity, normalized with all spectrum, of the Brillouin longitudinal acoustic mode, with depth focalisation for float glass (c) and pure silica (d). Size of spectrometer diaphragm is 300 µm

Figures 1 .

 1 Figures 1.a and 1.b show Brillouin spectra of silica and float glasses respectively for pristine glasses and at the centre of the indent surface for 10 N load. The spectrum of pristine float glass shows a well defined

Figure 3 Figure 2 :Figure 3 :

 323 Figure 3 the evolution of the position of the two peaks at the intensity maximum and their FWHM extension as a function of radial distance (3.a) and depth (3.b), for 10 N indent. The intensity of the colour, for each measurement is proportional to the relative weight of each fitted peak. As it had been demonstrated by Deschamps et al[START_REF] Deschamps | Elastic moduli of permanently densified silica glasses[END_REF], from ex-situ samples recovered after irrversible high pressure experiences and using Raman and Brillouin spectroscopies, it is possible to obtain the densification ratio as a function of the maximum pressure for silica glass. Using the literal expressions of Brillouin shift and densification ratio from this last paper, we can deduce the evolution of the densification ratio in the imprint as a function of

Figure 4 :

 4 Figure4: Evolution of Brillouin spectra of pristine silica in the centre of indent and at the surface of indent for various loads : 1N (blue 3N (orange curve) and 10 N (green curve). Brillouin spectrum of non-indented silica is also plotted (red line). Size of spectrometer diaphragm is 300 µm.

Figure 5 :

 5 Figure 5: Evolution of densified peak D area divided by sum of peak U and peak D area with focalisation depth for various load : 1 N (blue curve) 3N g (orange curve) and 10 N (green curve). Size of spectrometer diaphragm is 300 µm.

Figure 6 :

 6 Figure 6: Densification profiles of indented silica glass from various experimental works.a) cross-section densification map from Raman experiments panel, Perriot et al [21]) and corresponding map from finite element calculation (right panel,Kermouche et al [11]) b) shematic description of density profile beneath the indenter from chemical dissolution experiments probed by AFM (Niu c) Schematic iso-density map reflection the Brillouin spectroscopy results from this present study

Figure. 7 for 10 N

 710 imprint. In the lower right figure, the optimised evolutions of densification fonctions for the two models are plotted. In Figure7, three experimental Brillouin spectra for different focus depths (-2 µm (a),4 µm (b), and 10 µm(c)) and the modelled spectra are shown. It appears that the step model for densification

Figure 7 :

 7 Figure 7: Evolution of modeled Brillouin spectra (continuous densification in blue and step densification in orange) for focus laser depth : -2 µm (a),4 µm(b) and 10 µm (c) compared to experimental data. For experimental spectra, size of spectrometer diaphragm is 150 µm.

  experimental results therefore highlight the need to propose a new model for densification under indentation that explains the results observed and detailed in this paper. This work has also demonstrated the strenght of Brillouin spectrometry to refine our understanding of plastic deformation in glasses. As improving the mechanical strength of silicate glasses remains a key challenge, a step change in our understanding of the mechanisms involved is essential.