J J Lembrechts
email: jonas.lembrechts@uantwerpen.be.

De Gruyter

Ramirez-Rojas

Lenoir

Jonas Lembrechts

A framework to bridge scales in distribution modelling of soil microbiota

Keywords: species distribution models, soil microbes, biotic interactions, joint distribution models, niche modelling, microclimate

Creating accurate habitat suitability and distribution models (HSDMs) for soil microbiota is far more challenging than for aboveground organism groups. In this perspective paper, we propose a conceptual framework that addresses several of the critical issues holding back further applications. Most importantly, we tackle the mismatch between the broad-scale, long-term averages of environmental variables traditionally used, and the environment as experienced by soil microbiota themselves. We suggest using nested sampling designs across environmental gradients and objectively integrating spatially hierarchic heterogeneity as covariates in HSDMs.

Secondly, to incorporate the crucial role of taxa co-occurrence as driver of soil microbial distributions, we promote the use of joint species distribution models, a class of models that jointly analyze multiple species' distributions, quantifying both species-specific environmental responses (i.e. the environmental niche) and covariance among species (i.e. biotic interactions). Our approach allows incorporating the environmental niche and its associated distribution across multiple spatial scales. The proposed framework facilitates the inclusion of the true relationships between soil organisms and their abiotic and biotic environment in distribution models, which is crucial to improve predictions of soil microbial redistributions as a result of global change.

Introduction

Habitat suitability and distribution models (HSDMs) have turned into irreplaceable assets to study the spatial distribution and redistribution of species under given environmental circumstances, and changes thereof [START_REF] Elith | Species Distribution Models: ecological explanation and prediction across space and time[END_REF]. This broad group of statistical models relates known species occurrences (or presence-absence) with information about the environmental conditions at these locations [START_REF] Guisan | Predicting species distribution: offering more than simple habitat models[END_REF][START_REF] Elith | Species Distribution Models: ecological explanation and prediction across space and time[END_REF][START_REF] Jiménez-Valverde | Use of niche models in invasive species risk assessments[END_REF]. Based on statistically or theoretically derived response curves [START_REF] Guisan | Habitat suitability and distribution models: with applications in R[END_REF], they aim to define the environmental niche and related geographical range in which organisms can operate. Applications range from studies on the effects of anthropogenic climate change to predictions of biological invasions, and they have proven valuable across a wide range of organism groups both above and below the soil surface [START_REF] Martiny | Microbial biogeography: putting microorganisms on the map[END_REF][START_REF] Elith | Species Distribution Models: ecological explanation and prediction across space and time[END_REF][START_REF] Sato | Modelling the global distribution of fungal species: new insights into microbial cosmopolitanism[END_REF][START_REF] Bahram | Structure and function of the global topsoil microbiome[END_REF]. For belowground organisms in particular, however, the successful application of HSDMs needs to overcome important limitations concerning the scale, resolution and accuracy that can reliably and routinely be obtained. Indeed, many of the available global, high-resolution gridded environmental variables cover the above-ground world only, with many of the key drivers of soil microbial distributions (like soil pH and soil organic matter) insufficiently recorded, or at coarse resolutions only [START_REF] Fierer | The diversity and biogeography of soil bacterial communities[END_REF][START_REF] Mod | What we use is not what we know: environmental predictors in plant distribution models[END_REF]. This spatial mismatch between the environment as experienced by organisms on the one hand, and the data we have available to model their distribution on the other, is particularly pronounced for soil microbes. Their distribution, more than for most aboveground species groups, indeed depends largely on processes happening at a small sub-meter scale in the hard-tomeasure belowground environment.

To understand and predict soil functioning in a rapidly changing world, we need the ability to accurately model the spatial distribution patterns of soil microorganisms. To achieve this, we need to incorporate local-scale drivers and processes relevant to these species groups in our HSDMs, in order to improve predictions of large-scale distributional patterns, especially in response to global change. For example, accurately attributing local variance to abiotic drivers versus biotic interactions can aid our understanding of the importance of species interactions involving microbes as drivers of spatial distributions [START_REF] De Mesquita | Incorporating biotic factors in species distribution modeling: are interactions with soil microbes important?[END_REF], Ovaskainen et al., 2017).

Such an approach can help improve our predictions of the faith of ecologically important, but so far elusive groups including for example plant pathogens, mycorrhizal fungi and nitrogen fixing bacteria and their effects on ecosystem processes [START_REF] Martiny | Microbial biogeography: putting microorganisms on the map[END_REF]. While many of these questions have been covered in microbial distribution studies before, a clear framework is still lacking to incorporate the relevant processes at the local scale of sub-centimeters to metersin which microbes operatewith the large-scale spatial patterns across several kilometers that are the key focus in spatial ecology, without the need for an ever-increasing higher resolution of the used spatial data.

Belowground microbial distributions: the importance of the local scale

It has long been assumed that many soil microbiota have cosmopolitan distributions. This theory stems from the observation that they have large population sizes and short generation times resulting in high dispersal, as well as the capacity to disperse over long distances if they get into the air or water streams [START_REF] Frey | The spatial distribution of soil biota[END_REF]. Additionally, the actual drivers of their occurrence might vary more on a sub-centimeter scale than on the kilometer-scale traditionally used, as is the case for e.g. pore size, or biotic interactions with other soil microbiota. If soil microbial distributions would indeed relate only little to global-level environmental gradients, the validity of HSDMapproaches would be hampered, as it would be impossible to identify an environmental niche at the regional or global scale. However, recent research has demonstrated that spatial patterns of microbial diversity qualitatively similar to those observed for plants and animals exist, suggesting that soil microbes indeed "have a biogeography" [START_REF] Martiny | Microbial biogeography: putting microorganisms on the map[END_REF][START_REF] Hanson | Beyond biogeographic patterns: processes shaping the microbial landscape[END_REF][START_REF] Sato | Modelling the global distribution of fungal species: new insights into microbial cosmopolitanism[END_REF][START_REF] Peay | Dimensions of biodiversity in the Earth mycobiome[END_REF]. While the existence of such spatial structuring should not come as a surprise for host-dependent soil microbes like plant-associated mycorrhizae, the presence of a spatial distributional pattern is more telling for microbes like saprotrophs that operate independent of a plant host. Nevertheless, their high dispersal flexibility implies thateven more so than for most aboveground organisms -the range of environmental conditions and the nature of biotic interactions at the local scale is likely to be of particular importance for their distribution [START_REF] Sato | Modelling the global distribution of fungal species: new insights into microbial cosmopolitanism[END_REF][START_REF] Vos | Micro-scale determinants of bacterial diversity in soil[END_REF], and an efficient way to include this local environmental variation in HSDMs is thus crucial.

While existing HSDMs often rely heavily on air temperature averages obtained from global climate models (e.g. [START_REF] Aguilar | Ecological niche models reveal the importance of climate variability for the biogeography of protosteloid amoebae[END_REF][START_REF] Tedersoo | Global diversity and geography of soil fungi[END_REF], the use of disconnected environmental variables overlooking the small-scale nature of the community processes involved is unlikely to provide meaningful descriptions of the regional and global distribution of soil biota, and thus not fully exploit the potential of HSDMs for soil organisms [START_REF] Fierer | The diversity and biogeography of soil bacterial communities[END_REF]. Recent HSDMs studies indeed increasingly acknowledge the need for more environmental data at the scale and especially the location relevant to the studied organisms, as opposed to the traditional use of coarse spatiotemporal averages of environmental conditions [START_REF] Lembrechts | Incorporating microclimate into species distribution models[END_REF][START_REF] Lembrechts | Comparing temperature data sources for use in species distribution models: From in-situ logging to remote sensing[END_REF][START_REF] Zellweger | Advances in microclimate ecology arising from remote sensing[END_REF]. This issue becomes particularly important when one aims to use HSDMs predictively to investigate the effects of global (or local) change on species (re)distributions. Incorporating local conditions is for example critical in regard to microrefugia, where isolated populations can persist in a favorable microclimate amid deteriorating climatic conditions until the latter become favorable again [START_REF] Hannah | Fine-grain modeling of species' response to climate change: holdouts, stepping-stones, and microrefugia[END_REF], Lenoir et al., 2017). For many soil microbiota, these microrefugia might exist within centimeters to meters from where the species currently occurs due to the large local heterogeneity in conditions [START_REF] Veresoglou | Extinction risk of soil biota[END_REF], and HSDMs averaging out these local gradients are thus likely to largely overestimate local extinction incidences, as has been shown for other organism groups [START_REF] Lenoir | Local temperatures inferred from plant communities suggest strong spatial buffering of climate warming across Northern Europe[END_REF]. Similarly, this local heterogeneity can create stepping stones of favorable conditions that facilitate range expansion beyond that predicted at a coarse scale [START_REF] Lembrechts | Microclimate variability in alpine ecosystems as stepping stones for non-native plant establishment above their current elevational limit[END_REF]. The relevant data to counter these biased estimates is nevertheless rarely available.

The way forward to improve HSDMs for soil microbiota

In this paper, we introduce a conceptual framework that provides the missing link between the local-scale (abiotic and biotic) belowground environment and large-scale distribution patterns, building on the recent advancements in tools, measurement techniques and models (Fig. 1, see also [START_REF] Lembrechts | A framework to bridge scales in distribution modelling of soil microbiotaa modelled example[END_REF] for a simple overview of the followed steps in R). Keeping the strengths and limitations of HSDMs in mind, this framework builds on 1) smart in-situ measurements of local environmental heterogeneity and community distributions, 2) recent developments in spatial modelling of environmental conditions, and 3) the hierarchical integration of the local with the regional scale. It thus aims for a relatively simple way to incorporate the complexity of the belowground world into statistical HSDMs without dramatically increasing the sampling effort, and is especially suited for distributions driven by a complex interaction between local environmental conditions and biotic interactions which mechanistic distribution models cannot (yet) handle.

Before expanding further, it is noteworthy that the definition of a 'species', as traditionally used in HSDMs, is far less straightforward for soil microbes than it is for aboveground organisms, and that approaches to identify micro-organism identities differ for different microbial groups. The following framework focuses specifically on soil micro-organisms that are routinely identified as operational taxonomic units (OTUs) or amplicon sequence variants (ASVs). These OTUs are commonly used as proxies for species (e.g. [START_REF] Bahram | Stochastic distribution of small soil eukaryotes resulting from high dispersal and drift in a local environment[END_REF], yet recently, delineation of taxa has often been abandoned in favor of ASVs [START_REF] Knight | Best practices for analysing microbiomes[END_REF] as the unit of analysis, in particular when sequence heterogeneity is not too high. The latter avoids the arbitrary binning of sequences and allows comparisons across different studies, yet it would increase complexity of the models due to a higher number of identities many of which would be intraspecific, as would be the case for e.g. fungal ITS-sequences [START_REF] Heeger | Long-read DNA metabarcoding of ribosomal RNA in the analysis of fungi from aquatic environments[END_REF]. If needed, meaningful HSDMs can also still be produced by focusing on keystone taxa or higher taxonomic resolutions only.

Nested sampling of abiotic conditions

Global databases and gridded datasets are currently being developed for many soil variables (e.g. SoilGrids, [START_REF] Hengl | SoilGrids250m: Global gridded soil information based on machine learning[END_REF]. Even though in most cases database resolutions are still relatively coarse for application to the world of microbiota, the accuracy of the available data on many local environmental drivers is rapidly rising. For soil moisture, for instance, the first attempts to obtain global remotely sensed data at a coarse scale (3 and 9 km) are now appearing [START_REF] He | Will remote sensing shapte the next generation of species distribution models[END_REF]. For soil microbial distributions, however, large-scale gridded datasets are still several orders of magnitude too coarse in spatial resolution. For many other relevant variables, e.g. soil climate, users often still have to rely on the above-ground equivalent for higher resolutions (e.g. free-air temperatures from CHELSA at 1 km² resolution [START_REF] Karger | Climatologies at high resolution for the earth's land surface areas[END_REF] vs soil temperatures from ERA5 at 9 km² resolution (Copernicus Climate Change Service, 2019)).

The success of any application of HSDMs to soil microbial distributions thus fundamentally needs to rely on a sampling design that allows linking the local heterogeneity, inherently present in soils, with the microbial distribution at regional scales (Fig. 2). We therefore recommend a nested sampling approach, measuring the relevant environmental variables right where it matters for each specific study organism (e.g. measuring soil temperatures in-situ instead of relying on freeair large-scale interpolations [START_REF] Lembrechts | Incorporating microclimate into species distribution models[END_REF]. The idea of a nested sampling design is to capture the spatial autocorrelation signal that may exist both at a very fine spatial resolution as well as across larger spatial extents, due to the modifying impact of both climatic and biotic interactions on soil conditions [START_REF] King | High levels of microbial biomass and activity in unvegetated tropical and temperate alpine soils[END_REF][START_REF] King | Biogeography and habitat modelling of high-alpine bacteria[END_REF], Ovaskainen et al., 2017). Such nested sampling can not only provide a better mechanistic understanding of drivers of the local and regional distribution of a focal study organism, it is also a crucial requirement to use HSDMs for predictions under global change [START_REF] Mateo | Looking for an optimal hierarchical approach for ecologically meaningful niche modelling[END_REF].

These measurements should cover both the local (e.g. meters or even sub-centimeters apart, or at different depths in the soil) and regional scales (e.g. up to kilometers apart) [START_REF] King | Biogeography and habitat modelling of high-alpine bacteria[END_REF].

Importantly, one should not aim for the finest possible resolution of environmental variables across the whole landscape. Such a run for an increasing refinement in the spatial resolution would greatly reduce the efficiency of HSDMs over larger spatial extents, while not necessarily improving their accuracy [START_REF] Bennie | Seeing the woods for the trees -when is microclimate important in species distribution models?[END_REF][START_REF] Lembrechts | Incorporating microclimate into species distribution models[END_REF]. Instead, we suggest to explicitly keep local-scale heterogeneity as a variable in regional-scale distribution models when aggregating to a coarser resolution. This can for example be achieved by using environmental data with a fine spatiotemporal resolution from a selection of sites to estimate the variation around the mean within a certain measurement location both in space and time (Fig. 2a). For example, one can take into account the differences in the variability in soil moisture levels occurring within different environments along large-scale environmental gradients (e.g in mountains vs. flatlands).

Such estimates of within-pixel heterogeneity can be obtained with relatively simple aggregating techniques, for example by using the local standard deviation in addition to the mean, or by fitting correlative models that include within-pixel values as a covariate to incorporate uncertainty. If desired, more complex approaches can be used, for example using mechanistic models that describe the spatiotemporal variability in abiotic conditions in the system, or by calibrating HSDMs at the local scale using Hierarchical Niche Models (HNMs)distribution models that allow formally including information from different spatial scales [START_REF] Mateo | Looking for an optimal hierarchical approach for ecologically meaningful niche modelling[END_REF].

Converting nested sampling into gridded datasets

Importantly, measured variables at the sampling locations are only one part of the story. To allow spatial interpolations of a species' distribution between the sampling locations, gridded environmental datasetsboth of the means and the variation around these meansare necessary. The absence of such gridded datasets for belowground environmental variables with sufficient resolution and reliability is one of the fundamental issues currently hampering further improvements in HSDMs, especially so for belowground organisms. In what follows, we describe how the measured environmental conditions can be converted into gridded products with a spatial resolution far finer than what is currently used, including the even finer-scaled local heterogeneity at the sub-centimeter to meter scale.

Recent work has shown that remotely sensed drivers and proxies can accurately be linked as covariates to measured environmental variables in the soil by using hybrid models combining statistical correlations with mechanistic knowledge of the drivers of environmental variation (Fig. 2, [START_REF] Lembrechts | Incorporating microclimate into species distribution models[END_REF][START_REF] Zellweger | Advances in microclimate ecology arising from remote sensing[END_REF]. For example, soil moisture correlates strongly with local topography, for which highly accurate data at meter to centimeter-resolution can be obtained through remote sensing (satellite-or LiDAR-based DEMs) [START_REF] Sörensen | On the calculation of the topographic wetness index: evaluation of different methods based on field observations[END_REF]. By combining such remotely-sensed data with in-situ measurements, soil moisture can be modelled with increasingly high accuracy, e.g. up to 1 m² resolution across a 3 km² regional extent [START_REF] Kemppinen | Modelling soil moisture in a highlatitude landscape using LiDAR and soil data[END_REF]. Often, a resolution of 20 x 20 m up to 1 x 1 km is achievable using freely available satellite-based remotely-sensed datasets, like digital elevation models. If one applies these interpolation techniques to both the local averages and the within-pixel heterogeneity, both can be converted into gridded datasets with the desired spatial resolution (Fig. 2b). To accomplish this interpolation of the means and variation in these in-situ measurements, one can either use simple statistical modelling methods -making sure that model residuals are not spatially autocorrelated -such as generalized linear mixed-effects models, generalized additive mixedeffects models, boosted regression trees, or random forests (e.g. [START_REF] Kemppinen | Modelling soil moisture in a highlatitude landscape using LiDAR and soil data[END_REF], as well as spatially explicit geostatistical approaches such as spatial kriging or geographically weighted regressions (GWR) [START_REF] Fotheringham | Geographically weighted regression: the analysis of spatially varying relationships[END_REF][START_REF] Lembrechts | Comparing temperature data sources for use in species distribution models: From in-situ logging to remote sensing[END_REF]. The latter technique estimate model parameters at each geographical location by implementing a kernel and weighing explanatory variables by distance. Although such geostatistical approaches have a more explicit way of incorporating spatial structure than traditional generalized linear models or generalized linear mixed-effects models for spatial interpolation of environmental variables, they cannot be used to extrapolate microclimate outside the spatiotemporal extent covered by the data [START_REF] Lembrechts | Incorporating microclimate into species distribution models[END_REF]. While doing this at the scale of single regional studies has thus proven productive, globally coordinated efforts measuring and compiling such data on belowground conditionspreferably using a nested framework -are needed more urgently than ever [START_REF] Robock | The Global Soil Moisture Data Bank[END_REF][START_REF] Slessarev | Water balance creates a threshold in soil pH at the global scale[END_REF][START_REF] Hengl | SoilGrids250m: Global gridded soil information based on machine learning[END_REF]. Such efforts are indeed required to make a wider range of soil environmental variables available for interpolation at large spatial scales, and to calibrate emerging mechanistic models that allow predicting variables at every relevant scale [START_REF] Kearney | microclim: Global estimates of hourly microclimate based on long-term monthly climate averages[END_REF].

Using these techniques shown to be effective in studies on aboveground organisms, in combination with a nested sampling approach can help solve the persisting mismatch between the available gridded environmental data products and the environment as it is perceived by soil organisms themselves, and will greatly enhance the reliability of soil microbial HSDMs. It would indeed allow to significantly improve predictive distribution models, as it permits for microrefugia to be identified within a local pixel through the variation in observed values [START_REF] Hattab | Towards a better understanding of potential impacts of climate change on marine species distribution: a multiscale modelling approach[END_REF].

The number of samples needed to achieve this heavily depends on the local heterogeneity present in the studied variable (e.g. a higher sampling density might be needed for soil moisture in mountainous environments than in flat terrain).

Even if we obtain sufficiently detailed environmental variables to allow the calibration of accurate and fine-grained HSDMs, the question remains which of these environmental variables are driving the regional distribution of soil microbes. Importantly, which variables are relevant for a particular group of microbes should be decided on a case-by-case basis, and requires a thorough understanding of the microbial group under study. Indeed, these driversand the scale on which they operatewill depend among others on the size of the studied organism groups. For example, [START_REF] Birkhofer | General relationships between abiotic soil properties and soil biota across spatial scales and different land-use types[END_REF] showed that abiotic soil properties explained significant amounts of variation in fungal diversity, but not in yeasts or bacteria. The use of HSDMs should thus always be seen in parallel with other experimental approaches in which the mechanisms behind the spatial drivers are further disentangled [START_REF] Ettema | Spatial soil ecology[END_REF].

Nevertheless, we argue here that many available datasets of soil microbial diversity would be suitable for our approach. The idea of nested sampling per se is not new in HSDMs [START_REF] Diez | Hierarchical analysis of species distributions and abundance across environmental gradients[END_REF][START_REF] Elith | Species Distribution Models: ecological explanation and prediction across space and time[END_REF], and applications of similar sampling designs have been used in studies of soil microbial diversity (i.e. repeated local sampling along environmental gradients, or hierarchically structured designs to capture variation at different scales) (e.g. [START_REF] King | High levels of microbial biomass and activity in unvegetated tropical and temperate alpine soils[END_REF][START_REF] Bahram | Stochastic distribution of small soil eukaryotes resulting from high dispersal and drift in a local environment[END_REF][START_REF] Zhou | Temperature mediates continental-scale diversity of microbes in forest soils[END_REF]. The current scale of most sampling campaigns (collecting soil diversity and environmental information with an accuracy of meters) would indeed already provide sufficient detail to answer a myriad of questions regarding the true relationship between soil microbial spatial distributions and the local environment, especially if one would refrain from the common practice to pool all samples out of a plot together to create an 'average' sample, and thus allow objective integration of local-scale heterogeneity into models [START_REF] Manter | Negative effects of sample pooling on PCR-based estimates of soil microbial richness and community structure[END_REF]. Global database efforts compiling data on soil microbial distributions and/or their abiotic drivers can for example explicitly keep this hierarchy in their structure and, as much as possible, link up the diversity data with in-situ measured local-scale environmental conditions.

Emerging efforts to synchronize such sampling efforts across the globe (e.g. [START_REF] Kao | NEON terrestrial field observations: designing continental-scale, standardized sampling[END_REF][START_REF] Thompson | A communal catalogue reveals Earth's multiscale microbial diversity[END_REF] are thus pivotal to move this forward, as are efforts to standardize sampling protocols between studies [START_REF] Halbritter | The handbook for standardized field and laboratory measurements in terrestrial climate change experiments and observational studies (ClimEx)[END_REF] to reduce uncertainties when merging studies [START_REF] Ramirez | Detecting macroecological patterns in bacterial communities across independent studies of global soils[END_REF].

Pooling soil microbial communities

For the sampling of soil microbial communities, practical restrictions (time and money) often preclude the analysis of large numbers of samples. However, as summarized in Fig. 2c, repeated sampling of local communities within a plot is not a prerequisite for the successful application of our framework. Indeed, soil microbial studies traditionally measure the pooled community structure by analyzing a sample of mixed soil from a random set of locations within a plot and expressing OTUs by their (normalized) read number after sequencing [START_REF] Staddon | Microbial diversity and community structure of postdisturbance forest soils as determined by sole-carbon-source utilization patterns[END_REF][START_REF] Cleary | Long-term antibiotic exposure in soil is associated with changes in microbial community structure and prevalence of class 1 integrons[END_REF]. This approach on its own already provides an integration of the variability in local communities driven by the local environmental heterogeneity. Linking the pooled community datawhich is thus weighed based on the local distribution of environmental conditions -to the average and heterogeneity in local environmental conditions effectively takes into account the necessary fine-scale variation. It is however recommended to validate the local relationship between microbial communities as displayed in Fig. 2b and local environmental conditions in a few sites.

Importantly, note that OTU abundances are semi-quantitative relative abundances only. They are indeed only meaningful in relation to the community in which they are measured, but they do give information on shifts in dominance of a certain OTU between different samples. Additionally, the proposed approach still suffers from false absences in the community data, as do all spatial sampling schemes balancing output with resource input [START_REF] Jiménez-Valverde | The ghost of unbalanced species distribution data in geographical model predictions[END_REF]. The latter issue is not solved here, and would not necessarily benefit proportionally from increasing sample sizes.

Including biotic interactions in HSDMs of soil microbiota

A recurring pattern in analyses of the spatial distribution of soil microbes is that the local spatial distribution of many soil microbiota is not only determined by the abiotic environment, but also by the presence or absence of other taxa [START_REF] Ettema | Spatial soil ecology[END_REF][START_REF] Wardle | Ecological linkages between aboveground and belowground biota[END_REF][START_REF] De Vries | Abiotic drivers and plant traits explain landscape-scale patterns in soil microbial communities[END_REF][START_REF] Bahram | Structure and function of the global topsoil microbiome[END_REF]. Indeed, for many soil microbes spatial distribution largely depends on where competitors, facilitators, mutualists, predators and/or other interactors reside, especially at the local scale. This is for example the case for mycorrhizal fungi, which rely on the presence of relevant host plants (and vice versa), but also strongly interact with each other, other fungi, bacteria, and protists [START_REF] Sato | Modelling the global distribution of fungal species: new insights into microbial cosmopolitanism[END_REF]. To get reliable predictions of their spatial distributions, it will thus be important to include these local biotic interactions, within and between trophic levels, in our HSDMs.

Again, there are a variety of ways in which one can integrate these community effects. One could for example include plot-level metrics of the cohesion of interaction networks [START_REF] Herren | Cohesion: a method for quantifying the connectivity of microbial communities[END_REF] or metrics for food web complexity [START_REF] Pellissier | Combining food web and species distribution models for improved community projections[END_REF] as covariates in HSDMs.

Alternatively, a recently developed species distribution modelling tool now also allows to integrate the need for both an environmental and community-based approach in our HSDM-framework: joint species distribution modelling (jSDMs, see e.g. [START_REF] Ovaskainen | How to make more out of community data? A conceptual framework and its implementation as models and software[END_REF], Tikhonov et al., 2017, Fig. 1). The so-called jSDMs are a class of models that jointly analyze multiple species' distributions, quantifying both species-specific environmental responses (i.e. the environmental niche of the focal species) and covariance among species (i.e. positive or negative cooccurrences that can potentially capture biotic interactions). The use of jSDMs has proven a major step forward regarding the modelling of the distribution of a wide range of aboveground species groups, especially when incorporating a hierarchical study design [START_REF] Ovaskainen | How to make more out of community data? A conceptual framework and its implementation as models and software[END_REF]. It is important to remember, however, that jSDMs are not suitable for highly diverse communities of more than 100 species (with example studies commonly using up to around 50 co-occurring species, [START_REF] Ovaskainen | How to make more out of community data? A conceptual framework and its implementation as models and software[END_REF], Tikhonov et al., 2017). One should thus focus on specific groups of soil microbes, like nitrogen fixers, mycorrhizae or pathogens, and/or restrict analyses to the most common OTUs only. The latter is realistic, as the structure of belowground communities is often very uneven, with the top 2% of microbial taxa worldwide shown to make up 41% of microbial abundance [START_REF] Delgado-Baquerizo | A global atlas of the dominant bacteria found in soil[END_REF]. If one does want to model whole communities of thousands of taxa, other network based approaches, e.g. using plot-level interaction network metrics like the cohesion factor of interaction networks as mentioned above can be recommended.

For optimization of HSDMs for soil microbiota, we strongly recommend this integration of both abiotic and biotic drivers. Interestingly, this can work both for the interactions between soil organisms among each other, and between soil organisms and the aboveground world (e.g. between plant species and their mycorrhiza, if such data is available, see also Fig. 1). Vice versa, jSDMs can also be used to include the role of soil microbes as drivers of the distribution of aboveground organisms like plants [START_REF] De Mesquita | Incorporating biotic factors in species distribution modeling: are interactions with soil microbes important?[END_REF]. The result is a hybrid model that allows giving equal importance to specific species and local environmental conditions as drivers of species distributions.

Such integration of biotic interactions could also help bypass the issue of functional redundancy in soil microbial ecology, i.e. the idea that many soil microbial taxa would share the same function and thus are often only absent (or rare) because their redundant counterpart is common [START_REF] Allison | Resistance, resilience, and redundancy in microbial communities[END_REF][START_REF] Mori | Low multifunctional redundancy of soil fungal diversity at multiple scales[END_REF]. This would result in an apparent unfilling of their ecological niche, i.e. the absence of the species in an environment with suitable environmental conditions. However, if the presence of a certain redundant species affects the distribution of the focal species, this will be picked up by the species correlation matrix used in the jSDM and result in a high relative importance of biotic co-occurrences.

Conclusion

In this perspective paper, we argued that we should aim to model the current and projected distribution of soil microbiota at a relevant scale, resolution and measurement location, and with the relevant variables that drive their distributions. We propose using a nested sampling and hierarchical modelling approach for the spatial niche and associated distribution, with the possibility of including different drivers of the distribution of the studied microbiota at different scales, and accounting for (above-and belowground) biotic interactions, and refer to the recent advancements in high-resolution spatial environmental modelling. Applying this framework to existing datasets of soil microbial distributions -and keeping it in mind when setting up sampling networks and globally coordinated efforts -promises to resolve many remaining questions about the local, regional and global distribution of soil microbes. Eventually this approach will allow answering fundamental questions concerning microbial distribution and community assembly, but will also be relevant in applied work such as predicting abundance and resilience of soil-dwelling pests or beneficial soil microbes in agriculture. Importantly, this framework paves the way towards the prediction of changes in soil microbial distributions as a result of different global change drivers, for which first the true relationships between organismal distributions and the environment need to be unraveled. The strength of the approach lies in the nested sampling of above-and belowground communities and environmental data, and the consequent inclusion of those at relevant scales in HSDMs. Microbial OTUs and environmental conditions are sampled both at a high resolution in a selection of plots, and at a coarse resolution across a region (see Fig. 2). The OTU-distribution data is then used to obtain information on the distribution of a specific focal species of interest, as well as on its co-occurrences with other microbes in the microbial community. The environmental variables can provide information on interactions with aboveground organisms (e.g. between plants and their associated mycorrhizas), and on the species' abiotic niche. The latter can be converted into gridded data, both for coarser-scale averages and their fine-scale heterogeneity, by interpolating them using remotely sensed gridded environmental proxies. Then, the regional distribution of specific focal species of interest can be modelled using habitat suitability and distribution models (HSDMs), as a function of both these biotic interactions and the abiotic drivers, for example using joint species distribution models (see main text). Note that displayed details regarding sample sizes, resolution and spatial distribution of sample points are examples only, with actual requirements depending on the study system, study question and practical limitations. abiotic conditions and the specific abiotic niche of each OTU (visualized here using the same color gradient as for the environmental gradient). A pooled community composition sample will then have OTU abundances dependent on the relative frequency of each local environmental condition within a plot. As this pooled community data thus inherently contains the underlying local-scale variation in the microbial community, one can calibrate distribution models using pooled community data measured in the red and black dots depicted in (b). Finally, the distribution of each OTU can be modelled as a function of the predicted regional averages and heterogeneity in abiotic conditions using HSDMs.

Figure 1 .

 1 Figure 1. General overview of the proposed roadmap for HSDMs of soil microbiota at relevant spatial scales.

  18

Figure 2 .

 2 Figure 2. Detailed view of how the proposed nested sampling approach can be used to improve HSDMs, here exemplified for the link between microbial communities and local abiotic conditions (e.g. soil moisture).

  (a) Fine-scale measurements of the abiotic conditions of interest in a selection of small plots, distributed evenly along large-scale environmental gradients (e.g. in similar vegetation types, across topographic gradients in landscapes of 100 x 100 km), can be used to get an estimate of the local heterogeneity in environmental conditions. (b) These local relationships in a selection of plots (red dots) can be extrapolated across a whole region using coarser-grained averages of the same environmental variables (measured in the sites represented by black dots) and the relationships between these local conditions (averages and heterogeneity) and the large-scale environmental gradients. For example, one can extrapolate the 19 relationships in (a) using high-resolution digital elevation models and, if repeated in different vegetation types, gridded vegetation maps. (c) The locally present microbial community depends on the local range of

Acknowledgements

We gratefully acknowledge the support from the Research Foundation -Flanders (FWO) to JJL and the comments from the 3 reviewers.

Funding

This work was supported by a postdoctoral fellowship from the Research Foundation -Flanders (FWO) to JJL.

Figure legends


domain: sdu sdv

text:
Jean Goedert
email: jean.goedert@protonmail.fr

Romain Amiot
email: romain.amiot@univ-lyon1.fr

Didier Berthet
email: didier.berthet@museedesconfluences.fr

François Fourel
email: francois.fourel@univ-lyon1.fr

Laurent Simon
email: laurent.simon@unvi-lyon1.fr

Christophe Lécuyer
email: christophe.lecuyer@univ-lyon1.fr

Combined oxygen and sulphur isotope analysis-a new tool to unravel vertebrate (paleo)-ecology

Keywords: Geochemistry, Stable isotope, Biogenic apatite, Ecology, Fossil (Gingerich et al

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

U N C O R R E C T E D P R O O F Introduction

Background information

Vertebrate evolution has been many times punctuated by ecological transitions between terrestrial and aquatic (freshwater vs seawater) environments resulting in major radiation events: during the Late Devonian-Early Carboniferous, early tetrapods left and Gordon 1965;Gat 1984;[START_REF] Gat | The stable isotope composition of waters of the eastern Mediterranean Sea[END_REF]. Hypersaline lagoons or sabkhas (but also inland lakes, e.g. in East Africa) can also reach δ 18 O w values higher than + 2‰ (e.g., Gat and Levy 1978).

The δ 18 O w values of freshwaters mainly derive from those of meteoric waters (groundwater contributions being possible) whose ultimate source is seawater. Evaporation of seawater at low latitudes, distillation and cooling of the humid air mass during its transport towards high latitudes are responsible for the negative δ 18 O values of meteoric waters (Dansgaard 1964). At the global scale, the higher the latitude and altitude, the lower the δ 18 O values of rainfall and snow. These values are comprised between -6 and -2‰ at low latitudes and decrease down to about -15‰ at high latitudes, polar caps excluded. Oxygen isotope compositions of vertebrate biogenic apatite phosphate (δ 18 O p ) are linearly correlated with the oxygen isotope composition of their environmental waters [START_REF] Longinelli | Oxygen isotopes in mammal bone phosphate: a new tool for paleohydrological and paleoclimatological research?[END_REF][START_REF] Luz | Oxygen isotope variations in phosphate of biogenic apatites, III. Conodonts[END_REF]. Consequently, vertebrates living or ingesting different environmental waters will record in their bones distinct oxygen isotope compositions. Nonetheless, it is worth to note that physiological factors such as evaporative transcutaneous water loss and thermo-metabolism, which are species-specific, also impact the oxygen isotope compositions recorded in bioapatites (e.g. [START_REF] Kohn | Predicting animal δ18O: accounting for diet and physiological adaptation[END_REF][START_REF] Levin | A stable isotope aridity index for terrestrial environments[END_REF].

Sulphur isotope composition of vertebrate apatite

Sulphur isotope composition of sulphates (δ 34 S) is highly variable in modern aquatic environments. Marine environments have high and relatively uniform sulphate δ 34 S values close to + 21.0‰ (Böttcher et al. 2007). Most freshwater environments (e.g. rivers, lakes, ponds, precipitations) have comparatively lower sulphate δ 34 S values, ranging from -20.0 to + 20.0‰ [START_REF] Krouse | Sulphur isotopes in our environment[END_REF][START_REF] Kaplan | Stable isotopes of sulphur, nitrogen and deuterium in recent marine environments[END_REF]Nehlich 2015). It has been shown that the sulphur isotope composition of food is recorded in vertebrate organic tissues (e.g. muscles, hairs) or molecules (e.g. bone collagen) with low isotopic fractionation (+ 0.5‰ ± 2.4‰, Nehlich 2015), especially when compared to the oxygen isotopic system.

A recent study also measured very low sulphur isotope fractionation values between the collagen of sub-fossil red fox and that of its preys (ranging from -0.54 to + 0.03‰, with a mean analytical error of ± 0.4; [START_REF] Krajcarz | Prey-to-fox isotopic enrichment of 34S in bone collagen: implications for palaeoecological studies[END_REF]. Notably, this study further allows such low sulphur isotope fractionation to apply for carnivores.

Sulphur isotope analysis of vertebrate organic tissues is, therefore, particularly relevant to differentiate between freshwater and seawater environments. In particular, this method has been successfully used to determine the living environment exploited by

U N C O R R E C T E D P R O O F 140 fish

Oxygen isotope analysis

Bone apatite samples were treated following the wet chemistry protocol described by (Crowson et al. 1991) and slightly modified by [START_REF] Lécuyer | Thermal excursions in the ocean at the Cretaceous-Tertiary boundary (northern Morocco): δ18O record of phosphatic fish debris[END_REF]). This protocol consists in the isolation of phosphate (PO 4 3-) from apatite using acid dissolution and anion-exchange resin. For each sample, 30 mg of enamel powder was dissolved in 2 mL of 2 M HF overnight. The CaF 2 residue was separated by centrifugation, and the solution was neutralised by adding 2.2 mL of 2 M KOH. 2.5 mL of Amberlite™ anionexchange resin was added to the solution to separate the PO 4 3-ions. After 24 h, the solution was removed and the resin was eluted with 27.5 mL of 0.5 M NH 4 NO 3 . After (V-CDT), (Fourel et al. 2015;[START_REF] Goedert | High-precision 34S/32S measurements in vertebrate bioapatites using purge-and-trap EA-IRMS technology[END_REF]) to ensure that analytical conditions were optimal to perform sulphur isotope analyses of samples with low-S content.

The sample average standard deviation for δ 34 S measurements is 0.34‰ ± 0.34‰. Data are reported as δ 34 S in ‰ vs V-CDT.

271

Results

272

Oxygen isotope

273

The different vertebrates analysed had oxygen isotope com-274 positions ranging from + 12.1 to + 24.2‰ V-SMOW (Online 275 Information 1; Fig. 1a), which mainly reflect the variability of 276 oxygen isotope compositions of environmental waters. On the Fig. 2 Reconstructed oxygen and sulphur isotope composition of the environmental waters (δ 18 O w , δ 34 S w ) of the modern vertebrates. For oxygen, the isotopic composition of water was calculated using published isotopic fractionation equations for different groups of vertebrates (Online Information 3). For sulphur, the isotopic composition of water is very close to that recorded in bone apatite (i.e., almost no isotopic fractionation) and was calculated using published values of sulphur isotope composition of bone apatite and associated environmental water measured in present-day vertebrates [START_REF] Goedert | Euryhaline ecology of early tetrapods revealed by stable isotopes[END_REF]; Online Information 4). Each data point represents a biologically independent animal (n = 64) and corresponds to the average value of five and three repeated measurements for oxygen and sulphur isotope analysis, respectively (see "Material and Methods"). Each dot is numbered according to the species it represents (cf. Table 1). Error bars of each individual data point are given in Table S2 andS3 . 2003[START_REF] Levin | A stable isotope aridity index for terrestrial environments[END_REF]Tütken et al. 2006;Amiot et al. 2015Amiot et al. , 2009Amiot et al. , 315 2010;;Pouech et al. 2014;[START_REF] Guy | Tooth oxygen isotopes reveal Late Bronze Age origin of Mediterranean fish aquaculture and trade[END_REF]. As illustrated by 316 our results, this analysis can be particularly useful to distin- 2008). Although the different vertebrates sampled come from different region of the world, it should be noted for instance that the Eurasian otter (Lutra lutra) and the Eurasian beaver (Castor fibre), both having a semi-aquatic lifestyle, have recorded lower oxygen isotope ratios in their bones than fully terrestrial mammals (Online Information 1 and Fig. 1a). This is also the case for the semi-aquatic mallard duck (Anas platyrhynchos), which recorded in its bones lower oxygen isotope ratios than the common buzzard (Buteo buteo) (Online Information 1 and Fig. 1). In the latter case, it is worthy to note that both specimens come from the same geographic area and therefore rely on environmental waters of comparable oxygen isotope compositions.

On the contrary, it can be used to detect desert lifestyle [START_REF] Lécuyer | Oxygen isotope compositions of reptile bones and teeth: a potential record of terrestrial and marine paleoenvironments[END_REF]. For instance, the horned desert viper (Cerastes cerastes) and the Kleinmann's tortoise (Testudo kleinmanni), had both recorded high oxygen isotope ratios in their bones.

Nonetheless, for low-latitude environments, oxygen isotope compositions of freshwater and marine environments can display significant overlap. Consequently, water oxygen isotope compositions recorded in vertebrate apatites may not always be a diagnostic tracer of their living environment (e.g.

Pouech et al. 2014).

Sulphur isotope composition

Compared to oxygen, sulphur isotopes have been less applied to question the ecology of extinct vertebrates, principally due to technical difficulties. Due to the large amplitude of natural isotopic variations, particularly observed between terrestrial and marine environments, it remains a particularly relevant environmental tracer (cf. Background information).

However, as discussed in the "Introduction" section, the 'sea spray' effect may complicate interpretation concerning the living environment of vertebrates for terrestrial environment located in the influenced of marine ones. Moreover, some freshwater settings may have sulphur isotope compositions close to that of marine environments. For instance, rivers draining basins in which marine evaporites are exposed may have elevated dissolved sulphate content (more than 200 mg/L for the Colorado River system (Shope and Gerner 2014)) and δ 34 S values (up to seawater-like 19.5‰ for the Mackenzie River system [START_REF] Hitchon | Hydrogeochemistry of the surface waters of the Mackenzie River drainage basin, Canada-III. Stable isotopes of oxygen, carbon and sulphur[END_REF]). Therefore, vertebrates living in such environments are expected to have high sulphur isotope compositions that could be misinterpreted as reflecting an aqueous environment at least submitted to some marine influences. Finally, vertebrate species living in aquatic environments submitted to the influences of both fresh and

422

The general picture we have of major ecological transitions 423 that took place during vertebrate evolution are incomplete and 424 potentially biased as it corresponds to the final stages of these 425 transitions. For instance, the colonisation of terrestrial environments by early tetrapods at the beginning of the Carboniferous gave rise to a wide evolutionary radiation of terrestrial tetrapods that are still present on lands today. Similarly, the multiple iterations of secondary adaptation to the aquatic environment are well illustrated by the numerous species of vertebrates belonging to different groups (crocodiles, snakes, turtles, lizards, birds and mammals), which live in present-day aquatic environments. All these vertebrates testify that different groups adapted to new environments from a common ancestor.

However, the way these major ecological transitions proceeded, especially during their early stages, is difficult to infer and often remained elusive. Indeed, morpho-functional adaptations to a specific environment can be diachronous with its effective use (exaptation); the diagnose of living environment of vertebrates from morpho-functional analysis is thereby limited. Therefore, the combined use of 18 O/ 16 O and 34 S/ 32 S ratios of skeletal apatite should be particularly promising and powerful to document major ecological transitions in the fossil record for any phylogenetic group of vertebrates. For instance, this method has already been successfully applied to determine the aquatic environment of some Devonian early tetrapods and their associated vertebrate fauna [START_REF] Goedert | Euryhaline ecology of early tetrapods revealed by stable isotopes[END_REF]. Furthermore, it could also help to precise the ecology of some present-day aquatic vertebrates and shed light on the modalities of transition between terrestrial and aquatic environments during the course of vertebrate evolution over the Phanerozoic. It is also worthy to note that this method has the potential to shed light on the ecology of numerous present-day vertebrates living in transitional environments, and for which the ecology remains unclear.

30F

  aquatic environments and colonised terrestrial ones (Ahlberg and Milner 1994); during the Jurassic-Cretaceous, various 32 crocodylomorphs belonging to the thalattosuchians, the 33 pholidosaurids, the dyrosaurids and the eusuchians, radiated in 34 the marine environments (Martin et al. 2014); One hundred mil-35 lion years later, during the Cenozoic (Eocene), early cetaceans 36 also experienced a secondary adaptation to aquatic environments controlled by evaporation and condensation processes during 86 which isotopic fractionation takes place (Craig and Gordon 87 1965; Dansgaard 1964). Marine environments have relatively 88 uniform δ 18 O w values of 0 ± 1‰ except at high latitudes, where δ 18 O w values are lower, ranging from -3 to -1‰ due to mixing with ice melt, and at tropical latitudes where high evaporation rates result in positive δ 18 O w values ranging from + 1 to + 2‰, especially in closed tropical and subtropical seas like the Red Sea, the Dead Sea, Mediterranean Sea or Caribbean Sea (Craig

177

  Figs. 1 and 2 were plotted using KaleidaGraph 4.5.3 software. Figures were drawn using Inkskape 0.92.3.

  4 h, 0.5 mL of NH 4 OH and 15 mL of an ammoniacal solution of AgNO 3 were added, and the samples were placed in a thermostated bath at 70 °C during 7 h, allowing the precipitation of silver phosphate (Ag 3 PO 4 ) crystals.When only a few mg of apatite powders could be collected, the wet chemistry procedure was adapted following(Bernard et al. 2009) for small sample weights (about 3 mg).Oxygen isotope compositions were measured using a high-temperature pyrolysis (Py) technique involving a VarioPYROcube™ elemental analyser (EA) interfaced in continuous flow (CF) mode to an Isoprime™ isotopic ratio mass spectrometer (IRMS) (EA-Py-CF-IRMS technique(Fourel et al. 2011;[START_REF] Lécuyer | High-precision determination of 18O/16O ratios of silver phosphate by EA-pyrolysis-IRMS continuous flow technique[END_REF]) at the Laboratoire de Géologie de Lyon (UMR 5276, Université Claude Bernard Lyon 1). For each sample, 5 aliquots of 300 μg of Ag 3 PO 4 were mixed with 300 μg of pure graphite powder and loaded in silver foil capsules. Pyrolysis was performed at 1450 °C. Measurements were calibrated against the NBS120c (natural Miocene phosphorite from Florida: δ 18 O = 21.7‰ (V-SMOW),[START_REF] Lécuyer | Thermal excursions in the ocean at the Cretaceous-Tertiary boundary (northern Morocco): δ18O record of phosphatic fish debris[END_REF] and the NBS127 (barium sulphate, BaSO 4 : δ 18 O = 9.3‰ (V-SMOW),[START_REF] Hut | Stable isotope reference samples for geochemical and hydrological investigations[END_REF]. Silver phosphate samples precipitated from standard NBS120c were repeatedly analysed (δ 18 O p = 21.6‰; 1σ = 0.4; n = 16) along with the silver phosphate samples derived from fossil bioapatites to ensure that no isotopic fractionation took place during the 29 ± 0.14‰. Data are reported as δ 18 O p in ‰ 243 values vs V-SMOW. 244 Sulphur isotope analysis 245 Sulphur isotope compositions were measured using a 246 VarioPYROcube™ elemental analyser in NCS combustion 247 mode interfaced in continuous-flow mode with an Isoprime 248 100™ isotope ratio mass spectrometer hosted by the plat-249 form 'Ecologie Isotopique' of the 'Laboratoire d'Ecologie des Hydrosystèmes Naturels et Anthropisés' (LEHNA, UMR 5023, Villeurbanne, France). For each bone apatite sample, 3 aliquots of 7 mg of bioapatite powder were mixed with 20 mg of pure tungsten oxide (WO 3 ) powder and loaded in tin foil capsules. Tungsten oxide is a powerful oxidant ensuring the full thermal decomposition of apatite sulphate into sulphur dioxide (SO 2 ) gas (Goedert et al. 2016). Measurements have been calibrated against the NBS127 (barium sulphate, BaSO 4 δ 34 S = + 20.3‰ (V-CDT), (Halas and Szaran 2001) and S1 (silver sulphide, Ag 2 S δ 34 S = -0.3‰ (V-CDT), (Robinson 1995) international standards. For each analytical run of bone samples, we have also analysed BCR32 samples as a compositional

Fig. 1

 1 Fig. 1 δ 18 O p and δ 34 S apatite values

  Fig.2Reconstructed oxygen and sulphur isotope composition of the environmental waters (δ 18 O w , δ 34 S w ) of the modern vertebrates. For oxygen, the isotopic composition of water was calculated using published isotopic fractionation equations for different groups of vertebrates (Online Information 3). For sulphur, the isotopic composition of water is very close to that recorded in bone apatite (i.e., almost no isotopic fractionation) and was calculated using published values of sulphur isotope composition of bone apatite and associated environmental water measured in present-day vertebrates[START_REF] Goedert | Euryhaline ecology of early tetrapods revealed by stable isotopes[END_REF]; Online Information 4). Each data point represents a biologically independent animal (n = 64) and corresponds to the average value of five and three repeated measurements for oxygen and sulphur isotope analysis, respectively (see "Material and Methods"). Each dot is numbered according to the species it represents (cf. Table1). Error bars of each individual data point are given in TableS2 and S3for oxygen and sulphur respectively. Results are given as variations in parts per mille from the ratio of 18 O/ 16 O in Vienna Mean Oean Water (‰ VSMOW) for oxygen and 34 S/ 32 S in Vienna Canyon Diablo Troilite (‰ VCDT) for sulphur. Species living in freshwater are represented by light blue dots; those living in seawater are represented by dark blue dots, and green dots are used for terrestrial species.(1): Crocodylus niloticus; (2): Crocodylus siamensis; (3): Gavialis gangeticus; (4): Crocodylus porosus; (5): Crocodylus suchus; (6): Crocodylus suchus; (7): Acrochordus

  to note that vertebrates which 283 live in both freshwater to seawater environment had interme-284 diate median δ 18 O p values (δ 18 O p = + 17.7‰, 1σ = 0.9, n = 6), 285 although the difference was only significant compared to con-286 tinental environments and not seawater ones (P = 0.01255 and 287 P = 0.1611, respectively). One exception concerns the horned 288 desert viper (Cerastes cerastes) and the Kleinmann's tortoise 289 (Testudo kleinmanni), which had both recorded high oxygen 290 isotope ratios in their bones due to their desert lifestyle. 291 Sulphur isotope 292 The different vertebrates analysed had sulphur isotope com-293 positions apatite (δ 34 S apatite ) ranging from + 1.1 to + 22.9‰ V-294 CDT (Online Information 1; Fig. 1b). On the whole, verte-295 brates living or foraging in marine environments had signifi-296 cantly higher δ 34 S values than those living or foraging in con-297 tinental (freshwater or terrestrial) environments (median 298 δ 34 S apatite = + 16.9‰, 1σ = 4.4, n = 18 vs δ 34 S apatite = + 299 10.4‰, 1σ = 4.4, n = 40; P = 0.0001357). This isotopic pat-300 tern reflects an almost systematic 34 S-enrichment of marine 301 environments compared to continental ones. It is again worth 302 to note that vertebrates living in freshwater to seawater envi-303 ronment had intermediate median δ 34 S apatite values (δ 34 S = + 304 13.8‰, 1σ = 6.0, n = 6), although the difference was not sig-305 nificant with that of continental or marine environments (P = 306 0.1063 and P = 0.5264). Sulphur isotope analysis of fossilised 307 apatite can, therefore, help to detect the presence or proximity 308 of seawater in the living environments of extinct vertebrates. of vertebrate biogenic apatite has 312 been widely applied to fossilised apatite of extinct vertebrates 313 to get information on their living environment (e.g. Clementz 314 et al

317

  guish vertebrates living or foraging in marine environments 318 from those living or foraging in continental (freshwater or 319 terrestrial) ones (e.g. sharks: Gates Q6 2019; mosasaurs: Makádi 320 et al. 2012; coelacanths: Simon et al. 2003). 321 It can also be used to further differentiate aquatic or semi-322 aquatic lifestyle from a terrestrial one in the case of sympatric 323 vertebrates (e.g. Amiot et al. 2010). Indeed, terrestrial animals lose more water than semi-aquatic animals through transcutaneous evaporation or sweat. Water lost during this process as vapour is preferentially 16 O-enriched, resulting in a relative 18 O-enrichment of the residual body water (Cerling et al.

  like in estuaries, may record a sulphur isotope 376 composition in their bioapatite difficult to correctly interpret in 377 terms of living environment. 378 Combined oxygen and sulphur isotope composition 379 On the whole, the combined use of oxygen and sulphur isotope 380 compositions of bone apatite allows, in most cases, environmen-381 tal identification for the present-day vertebrates after the conver-382 sion of the measured δ 18 O and δ 34 S values of apatite into envi-383 ronmental water δ 18 O value and dissolved environmental sul-384 phate δ 34 S values using known isotopic fractionation equations 385 (Fig. 2; Online Information 3 and 4).386 The complementarity of these two isotopic systems lies in the 387 different abundance ratios of oxygen and sulphur, respectively, in 388 seawater and freshwater bodies. Indeed, oxygen is equally pres-389 ent (as H 2 O) in both marine and freshwater reservoirs whereas 390 sulphur content (as SO 4 2-) of seawater is generally 100 to 1000 391 higher than in freshwater (Fry and Chumchal 2011). 392 Consequently, sulphur isotopes will be particularly relevant to 393 detect the presence of seawater in the environment, even if only 394 a small quantity of seawater intrudes freshwater environment, 395 and oxygen isotopes will be relevant to quantify the amount of 396 freshwater in the environment, in particular in aquatic environ-397 ments where freshwater and seawater are mixing, like in deltas or 398 estuaries (Goedert et al. 2018). 399 Vertebrates living or foraging in marine environments tend to 400 have higher oxygen and sulphur isotope compositions recorded 401 in their bone apatite than those from freshwater and terrestrial 402 habitats. This rule is especially valid when we compare verte-403 brates of close phylogenetic affinity. For instance, the wild gha-404 rial (Gavialis gangeticus), living in freshwater streams, and the 405 two captive specimens of desert crocodiles (Crocodylus suchus), 406 kept in freshwater at the Zoo of Lyon, have recorded in their bone 407 apatite δ 18 O p and δ 34 S apatite values (+ 12.1‰ and + 15.9‰, + 408 14.5‰ and + 8.2‰, and + 13.9‰ and + 8.5‰, respectively) 409 lower than those measured in bones of the wild Nile crocodile 410 (Crocodilus niloticus; + 17.5‰ and + 20.1‰) and saltwater 411 crocodile (Crocodilus porosus; + 16.8‰ and + 16.3‰), both 412 known to undertake incursions in brackish waters to seawaters 413 (cf. Supplementary Information). Similarly, the sea otter 414 (Enhydra lutris), fully adapted to life in seawater, has higher 415 δ 18 O p and δ 34 S apatite values (+ 17.4‰ and + 18.1‰) than those 416 of the Eurasian otter (Lutra lutra) (δ 18 O p = + 14.0‰ and 417 δ 34 S = + 12.8‰), inhabiting freshwater environments. In a simi-418 lar way, the marine narwhal (Monodon monoceros) has higher δ 18 O p and δ 34 S apatite values (+ 17.0‰ and + 16.0‰) than those of 420 the South Asian river dolphin (Platanista gangetica; + 14.3‰ and + 10.2‰).

and methods 172 Sixty-four vertebrate bone apatite samples have been col- lected and analysed in this study (Online Information 1). 174 Samples were selected to encompass a broad ecological

  at the species and population levels (Fry 2002; Fry and 141 Chumchal 2011;[START_REF] Hesslein | Stable isotopes of sulphur, carbon, and nitrogen as indicators of trophic level and fish migration in the lower Mackenzie River basin, Canada[END_REF] Nehlich et al. 2013;

		142	Trembaczowski 2011) or in archaeological studies to know if
		143	ancient human populations relied on freshwater or marine food
		144	resources (e.g. Bocherens et al. 2016). More generally, terrestrial
		145	environments (including freshwater ones) and animals living
		146	there have generally relatively low δ 34 S values compared to ma-
		147	rine environments. Nonetheless, it is worth to note that coastal or
		148	island environments may be substantially influenced by sulphate
		149	from marine environments, which can be redeposited as rain or
		150	aerosols (the so-called 'sea spray' effect) with sulphate δ 34 S
	Q2	151	values close to those of marine environment (+ 20.3‰; Nielsen
		152	1974; Norman et al. 2006). Consequently, the δ 34 S values of
		153	vertebrates living in those terrestrial environments submitted to
		154	sea spray effect can be relatively high and may complicate inter-
		155	pretation concerning the living environment.
		156	Due to technical difficulties, sulphur isotope analyses have
		157	been only applied to organic tissues that easily degrade after
		158	animal death and are rarely preserved in the fossil record. A
		159	new method has been recently developed to measure the sul-
		160	phur isotope ratios ( 34 S/ 32 S) of sulphate compound in calcium
		161	phosphate minerals (analytical precision equals 0.5‰ (1σ))
		162	with a low-S concentration (0.14% to 1.19%) such as verte-
		163	brate bioapatites (Fourel et al. 2015; Goedert et al. 2016).
		164	Previous results indicated that sulphur isotope compositions
		165	of environmental waters are recorded in vertebrate inorganic
		166	tissues (bone apatite) with low isotopic fractionation (0.8‰ ±
		167	0.8‰, n = 5; Goedert et al. 2018). Therefore, sulphur isotope
		168	analysis of bioapatite from extinct vertebrates can provide
		169	estimates of the salinity of their aqueous environments
		170	(Goedert et al. 2018).
		175	and taxonomic spectrum of vertebrates (crocodiles, snakes,
		176	turtles, mammals, birds, lizards, fish and amphibians). For

171

Material

JrnlID 114_ArtID 1664_Proof# 1 -26/12/2019

Bocherens H, Drucker DG, Haidle MN et al (2016) Isotopic evidence (C,

Acknowledgements We thank the Musée des Confluences de Lyon, M. Creuzé des Châtelliers (Centre de Ressources pour les Sciences de l'Evolution (CERESE, FED 4271, Université de Lyon, Université Claude Bernard Lyon 1)), G. Douay (Zoo de Lyon), O. de Lataillade (Ferme du Ciron), P. François (Pierrelatte), E. Liatout (Maison Liatout) and the fishery (Maison Pupier) for providing bone material. Sulphur and oxygen isotope compositions were measured at the platforms 'Ecologie Isotopique' (LEHNA) 'Isotopes Stables' (LGLTPE), respectively.

Funding information This study was supported by the CNRS INSU program InterrVie, and the Institut Universitaire de France (C.L.). Data availability All data generated or analysed during this study are included in this published article [and its online information files].

Compliance with ethical standards

Author contribution statement All authors contributed to the study conception and design. Material preparation and data collection were performed by J. Goedert, D. Berthet and R. Amiot. Material analysis were performed by J. Goedert, F. Fourel and L. Simon. The first draft of the manuscript was written by J. Goedert, R. Amiot and C. Lécuyer, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Conflict of interest


Natkamon Tovanich
email: natkamon.tovanich@irt-systemx.fr

Nicolas Soulié
email: nicolas.soulie@imt-bs.eu

Petra Isenberg
email: petra.isenberg@inria.fr

Visual Analytics of Bitcoin Mining Pool Evolution: On the Road Toward Stability?

Keywords: Bitcoin, Mining pools, Visual analytics

We present our work on visual analytics tools to support the analysis of Bitcoin mining pool evolution. Mining blocks are a critical component of the Bitcoin ecosystem, helping to keep the system secure, valid, and stable. At the same time, mining is a resource-intensive activity that continues to get more and more difficult. Mining pools have emerged to address this issue and to ensure a more stable and predictable income by sharing computing power. Yet, increased centralization of the mining power is also not without dangers (e. g., the 51% attack), and, thus, it is important to better understand and analyze mining pool activities in Bitcoin. Here, we report three contributions: our extensive data collection on Bitcoin mining pools, our development of two custom visualizations, and our first exploratory data analysis leading to hypotheses and documented activities about pools' main features such as market share, reward rules, or location.

I. INTRODUCTION

The sharp increase of Bitcoin's value has attracted great attention on cryptocurrencies and their underlying technology, the Blockchain. In the Bitcoin ecosystem, miners are major actors from cryptographic as well as economic viewpoints. Miners ensure the validity and security of the Blockchain, which is fundamental for people's trust in the cryptocurrency. Mining is a computing and energy-intensive activity that is rewarded by newly created bitcoins in case of successful mining [START_REF] Böhme | Bitcoin: Economics, technology, and governance[END_REF]. While this incentive scheme is fundamental for a proper functioning and growth of Bitcoin, the required resources and unpredictability of mining success has led miners to join their resources in pools. Mining pools gather the mining power of all members to maximize the probability of successful mining, and so, expected rewards [START_REF] Schrijvers | Incentive compatibility of bitcoin mining pool reward functions[END_REF]. However, concentration of mining power in the hands of very few pools raises security issues, such as an increased likelihood of a 51% attack [START_REF] Eyal | Majority is not enough: Bitcoin mining is vulnerable[END_REF]. Despite the potential danger, mining pools have grown steadily and become critical supporters of the Bitcoin ecosystem.

Understanding the sustainability of the organization of mining activities in pools is critical to assess the future of Bitcoin, and more generally, of blockchain-based technology. A major component of the sustainability of mining systems is the economic efficiency of pools. Yet, it is always difficult to assess the efficiency of an organization. In a competitive context, the stability of organizations that operate on a market tends to illustrate their relative efficiency. In our case, the stability of mining pools and their distributions can illustrate the efficiency of the system. If new pools easily enter the market, they generate visible turbulence in the system and indicate that the mining ecosystem is not yet efficiently working.

In this article, we make three contributions: (a) we contribute an extensive data collection on Bitcoin mining, (b) we developed custom visualizations dedicated to analyzing the emergence and evolution of mining pools from our data, and (c) we report on findings from our first exploratory analysis using our data and visualizations. Our exploratory analysis shows that mining pools tend to converge toward a limited set of configurations (reward rule, location, hash rate, etc.). In particular, we show that mining market share distributions have become more stable across time and are now operated by a few numbers of pools. Moreover, pools have converged towards a main specific reward payout rule called Pay-Per-Share (PPS). Concerning the geographic location of mining pools, historically, the dominant pools have been hosted in China but tend now to become more globally operated. All this evidence suggests that mining activities have undergone a rationalization process and might currently provide a solid foundation for the stability of the Bitcoin ecosystem.

II. RELATED WORKS

As mining is central to how Bitcoin and other cryptocurrencies work, researchers have started to analyze it more closely. Much of the research on Bitcoin mining and mining pools is theoretical in nature and focuses on analyzing reward rules [START_REF] Schrijvers | Incentive compatibility of bitcoin mining pool reward functions[END_REF], [START_REF] Rosenfeld | Analysis of bitcoin pooled mining reward systems[END_REF], [START_REF] Fisch | Socially optimal mining pools[END_REF], mining strategies [START_REF] Houy | The bitcoin mining game[END_REF], [START_REF] Lewenberg | Bitcoin mining pools: A cooperative game theoretic analysis[END_REF], and attacks [START_REF] Eyal | Majority is not enough: Bitcoin mining is vulnerable[END_REF], [START_REF] Kroll | The economics of bitcoin mining, or bitcoin in the presence of adversaries[END_REF]- [START_REF] Laszka | When bitcoin mining pools run dry[END_REF].

Yet, empirical analyses on mining pools have also recently emerged. Wang and Liu [START_REF] Wang | Exploring miner evolution in bitcoin network[END_REF] provided evidence that the top mining pools became larger as the hash rate grew exponentially between Mar. 2013 and Mar. 2014. The authors analyzed the profit of mining regarding hardware cost and electricity price and concluded that the profit became negative when the hash rate increased faster than the Bitcoin price. Romiti et al. [START_REF] Romiti | A deep dive into bitcoin mining pools: An empirical analysis of mining shares[END_REF] analyzed the distribution of mining pools from Dec. 2013 to Dec. 2018 and found that 3-4 mining pools controlled >50% of the hash rate. The authors further analyzed the reward payout among the top-3 mining pools and found that a small number of members received a total of >50% of the reward from the pool. In addition, the authors detected cross-pool miners who received rewards from multiple mining pool and provided evidence that miners tend to transfer their rewards to exchange services and wallet providers. Wang et al. [START_REF] Wang | Measurement and analysis of the bitcoin networks: A view from mining pools[END_REF] analyzed the daily hash rate of the top mining pools from Feb. 2016 to Jan. 2019. They found that mining pools increased their hash rate exponentially in order to maintain their market share. Mining pools were caught in the Prisoner's Dilemma-if they increase their hash rate to compete with other pools their mining profit diminishes. Finally, Wang et al.'s data confirmed that pools tended to collect transactions that gave large transaction fees to maximize their profit.

Instead of focusing on just a few mining pools and a limited timeframe, our work allows to study mining pool distribution more broadly using a number of different parameters. As such, our work has some overlap with previous works (e. g., Romiti et al.'s study of mining pool distribution) but also extends the work with custom visualizations that show not only the distributions but also help to detect changes in the mining pool ranks over time. We also provide a first exploratory analysis of the evolution of reward rules and locations which has not been studied in previous empirical works. Our results provide new information toward more realistic theoretical models of the Bitcoin mining organization.

III. DATA DESCRIPTION

To analyze the evolution of Bitcoin mining over time, we collected coinbase transactions (the first transaction in each block) from the Bitcoin genesis block (the first ever mined block) until the last block mined in 2019. From the coinbase transactions, we extracted a number of different metrics as described next. We also collected monthly Bitcoin statistics from Blockchain.info and historic data about mining pools from the Bitcoin Wiki to give context to the mining activity.

A. Mining pool reward distribution

We tagged each coinbase transaction we collected with one mining pool (or "unknown") based on the data from Romiti et al. [START_REF] Romiti | A deep dive into bitcoin mining pools: An empirical analysis of mining shares[END_REF] who collected and manually resolved pool names from various sources (i. e. Blockchain.info, BTC.com, and Blocktrail.com). For each transaction we calculated a mining reward as the summation of the pre-defined block reward and the variably assigned transaction fees for the block. We converted the mining reward from Bitcoin currency to US Dollar using the daily market price data from Blockchain.info. Finally, we aggregated the number of blocks and mining rewards for each mining pool per month and calculated each pool's market share-the percentage of blocks mined by the pool compared to the total blocks in the month.

B. External mining pool data

Some important information about mining pools is not available for extraction from the Bitcoin blockchain. For example, data on reward payout schemes (how are reward payments distributed to members), mining pool locations, or fees are specific to how each pool operates. To obtain additional information, we collected tabular data available in the Bitcoin Wiki 1 . The wiki maintains an overview table of mining pools that includes information regarding pool name, location, size, reward type, transaction fees, among others. The table is maintained by the community and regularly changed. We crawled the editing history of the construct panel table about mining pools. We simplified the panel data to monthly granularity using the most recent edit for each month. For the mining pools that did not exist in the table, we manually searched for information from the pool's website and mining pool comparison articles from the internet 2 . The reward rules are categorized into main groups (see Fig. 2). The location of mining pools are grouped by continent-level (see Fig. 3). We defined "global" for pools that operated in more than two continents and kept China separately.

Finally, we collected Bitcoin statistics on hash rate and mining difficulty, transactions (number per block, mempool size, etc.), and market price from Blockchain.info and aggregated the data to obtain averages for each month.

IV. VISUAL ANALYTICS OF BITCOIN MINING POOLS

We built two custom visualizations to analyze mining pool evolution based on a user-centered design process. In this process two computer scientists collaborated closely with an economist (all authors of this paper). We identified design requirements and goals and iteratively revised prototype visualization designs. The two main analysis tasks we identified in the process are centered on: (1) the emergence and evolution of mining pool market shares and (2) external factors influencing the evolution of mining pool market shares and their distribution.

A. Visualizing the emergence and evolution of mining pool market shares

In order to learn when mining pools started to gain market share and which pools dominated the market over time, we focused our first visualization on temporal mining activity. Fig. 1 (a) provides information on the historic evolution of the top 25 mining pools based on the average market share over the entire period. This visualization allows the viewer to see an overview of active mining pools row-by-row. For each mining pool, we see when it started to mine blocks, how it accumulated or lost mining power, and whether it stopped its operations. For each pool and each month in the timeline, we draw a circle representing its market share and use a heated color scale to represent mining power domination (0%->50%) to highlight the pools that were likely to dominate the Bitcoin network at a time. We chose this representation over a simple heat map alternative because it allows us to represent two quantitative variables.

B. Visualizing external factors influencing the evolution of mining pool market shares and their distribution

Understanding when and why mining pool market shares change is important to understanding the broader mining ecology in Bitcoin. To address this analysis challenge, our second visualization presents the temporal distribution of mining pools in connection to the panel data we collected (reward rules, location, etc.). Fig. 2 and3 show the temporal evolution of reward rules and pool locations weighted by pools' market shares. For the visualization, we used a ribbon chart designa stacked bar chart showing a quantitative measurement over time with ribbons connecting the same data category. Each mining pool is represented as a bar whose height is relative to the pool's market share in a month and the color indicates categorical data (i. e. mining pool, reward rule, or location). For each month, the bars are sorted by the total market share of the category from the highest to the lowest value and within each category pools are sorted by their market shares. This sorting helps us to identify large mining pools or categories that dominated mining power. Ribbons help to track a mining pool's ranking connection between months. In contrast to our first visualization, here, we can see overall distributions and detect the top mining pools while the first visualization is better for observing the development of individual pools.

V. FINDINGS: MAIN TRENDS IN MINING POOL EVOLUTION

Next, we highlight exploratory findings we made using the visualizations we created. In particular, we saw that mining has undergone important changes from its beginning to the way people engage in it currently. Three main characteristics allow to highlight different stages of mining pool evolution: market share concentration, reward rules, and location.

A. Mining pool market share

In economy, considering the evolution of market share distributions gives a good overview of an activity's stability. The analysis is typically done through longitudinal studies of the market concentration index. Beyond economic interests, the question of mining activity concentration in Bitcoin is important to understand potential security issues (consensus manipulation). The Herfindahl index 3 provides a basic measure of the concentration of an activity or market. This aggregated measure does not, however, provide information about the cause of index variation. Fig. 1 mitigates this drawback by combining the Herfindahl index with the market share of each mining pool, allowing for more relevant analysis of the mining market evolution.

Until 2011, most mining activities were performed by small independent miners following Bitcoin's original spirit 4 . The 3 Herfindahl index (H) in Fig. 1 (b) is computed on a monthly basis. For each month, let's n be the number of active mining pools and x i pool i market share (i. e. the number of blocks mined by pool i divided by total number of blocks found in the month), then:

H = n i=1 x 2
i . H is equal to 1 for perfectly concentrated market and converges toward 0 for very fragmented ones. 4 It is known that some pools existed between 2009 and 2011, however they are not identifiable in our data and were very small. mining market was then very fragmented as illustrated by the Herfindahl index close to 0 (Fig. 1 (b)). According to our data, four cycles of concentration of market share can be observed. All these cycles are associated with a rise of bitcoin's value, as investments in mining hardware are indeed strongly correlated with bitcoin's value [START_REF] Prat | An equilibrium model of the market for bitcoin mining[END_REF]. By increasing the expected revenue of mining, a rise of bitcoin value gives an incentive to miners to enhance their investments and might also attract new miners into pools [START_REF] Wang | Measurement and analysis of the bitcoin networks: A view from mining pools[END_REF]. A primary cycle of concentration started mechanically in the early 2011 with the emergence of first large and identifiable pools, namely DeepBit and Eligius. This rise got amplified by the increase of DeepBit's hash power during the first peak of bitcoin's value that occurred in Jul. 2011. At this moment bitcoin's value reached more than $30, compared to roughly $1 few months ago. The growing expected revenue gained from mining lead probably miners to create new large pools such as BTCGuild and SlushPool. The emergence of these two pools participated to the rise of market share concentration until the beginning of 2012. The second cycle began in 2013. Two peaks occurred in Apr. and Nov., during which bitcoin's value respectively reached more than $260 and $1,200. During this period, existing pools increased their hash power (e. g., BTCGuild, SlushPool or 50BTC), and new pools experienced significant growth (F2Pool and GHash.IO), especially in 2014. After the decrease of bitcoin's value in 2015, the market price raised again and reached peaks at roughly $450 in Dec. 2015 and even $750 in Jun. 2016. During this period, F2Pool and AntPool increased largely their hash power. Newly created pools (Bitfury, BTCC Pool, and BW.COM) also enhanced their hash rate. The fourth cycle was driven by the sharp rise of Bitcoin market price at the end of 2017. Most existing large pools, and in particular BTC.com and AntPool, increased their hash rates. Both pools are owned by BitMain, the Chinese Bitcoin mining hardware manufacturer.

It is noticeable that the four cycles have a decreasing amplitude along time, indicating that the market might progressively reach an equilibrium situation. This is due to the exponential increase of hash power which involves very large investment for an incumbent or a new pool to get a significant market share. Since mid-2018, the concentration of mining activities has been rather stable, with 4 main pools accounting each for 10% or more of market share: F2Pool, AntPool, BTC.com, and Poolin.

B. Reward rules and location

Another striking evidence of the rationalization of mining lies in the evolution of reward rules applied by pools. A wide variety of reward rules exist ranging from simple methods (Proportional (Prop.), Pay-Per-Share (PPS), etc.) to more complex ones (Double Geometric Method (DGM), etc.) [START_REF] Rosenfeld | Analysis of bitcoin pooled mining reward systems[END_REF]. This rule variety leaves rooms for strategic decisions by pool managers in order to adopt a reward scheme that attracts miners and makes them behaving loyally, e. g., to avoid block withholding, multi-pooling or pool hopping [START_REF] Schrijvers | Incentive compatibility of bitcoin mining pool reward functions[END_REF], [START_REF] Belotti | Bitcoin pool-hopping detection[END_REF]. Ribbon charts are useful tools to observe the evolution of an organizational feature such as reward rules. In Fig. 2, we can observe a large heterogeneity of applied reward rules until 2016 5 . Progressively, two rewards schemes became dominant patterns with PPS and PPS+ (i. e. PPS for block rewards and PPLNS for transaction fees in most cases). Compared to proportional schemes (Prop., Score, etc.), PPS and PPS+ are less favorable to block withholding [START_REF] Schrijvers | Incentive compatibility of bitcoin mining pool reward functions[END_REF] and offer resistance to pool-hopping [START_REF] Belotti | Bitcoin pool-hopping detection[END_REF]. These reward rules seems to provide an attractive and stable incentive scheme for miners. Pools using PPS or PPS+ can differ however in the fee they apply on reward and transaction fees.

Fig. 3 shows the large evolution in mining pool locations. Until 2015, most pools were located in Europe and the US. The rapid growth of two important Chinese pools (F2Pool and AntPool) profoundly modified this landscape and made China the largest pool hosting nation from 2015 to 2018. According to Hileman and Rauchs [START_REF] Hileman | Global cryptocurrency benchmarking study[END_REF], cheap electricity and land costs in remote Chinese areas (e. g., Sichuan) are major 5 For detailed descriptions of the reward rules, see: https://en.bitcoin.it/wiki/ Comparison of mining pools. drivers of this location pattern. Between Mar. 2015 and Feb. 2017, the combination of Chinese pools exceeded 51% mining power threshold and therefore posed the risk of majority attack. This situation persisted until 2019 when global pools became dominant. This growth resulted from both the growth of originally global pools (BTC.com and Poolin), and also from Chinese or European pools that turned into global ones (F2Pool and SlushPool).

VI. DISCUSSION AND CONCLUSION

Mining pools are the heart of the Bitcoin ecosystem's security and growth. Their evolution toward rational and stable organizations is critical for Bitcoin's future. Using visual analytics approaches, we collected and joined a large number of datasets about mining pools and developed two custom visualizations to investigate their evolution according to a number of different parameters. In particular, we combined the aggregated measure of market share with pool-level information across time. Our representations allow us to document the rise of mining activity concentration and give rise to insights and further hypotheses as to the cause of this evolution. Important events regarding mining activity profitability (e. g., halving days and bitcoin value) seem to be major drivers of changes to the Bitcoin ecosystem.

We found the use of ribbon charts insightful to show two other important trends regarding pool organization. First, we can observe the rise of PPS and PPS+ as the main reward rules used by pools. These rules seem to provide attractive and stable incentives for miners. Secondly, ribbon charts are very useful to highlight the emergence of China as the main hosting country for pools and show that global pools have become dominant in recent times. These findings raise questions worth further examination, in particular regarding factors that stimulate the stability of Bitcoin mining.

The visualizations we developed are relatively simple but already provide first useful tools for people who want to analyze the evolution of mining pools by highlighting relevant elements that affect these organizations. We are working on a more encompassing visualization tool to help economists develop a more realistic model of Bitcoin mining.

Fig. 1 .

 1 Fig. 1. (a) The emergence and evolution of the top 25 mining pools. The figure shows the market share of each mining pool as the size of the circle over time. The color scale encodes the mining power domination in the Bitcoin network from no domination (0% mining power) in yellow to completely domination (>50%) in red. Two grey vertical lines in both charts indicates two halving days on Nov. 28, 2012 and Jul. 9, 2016. (b) The Herfindahl index measures the concentration of mining pools in the Bitcoin network. (c-e) Bitcoin mining statistics data from Blockchain.info. (c) Total hash rate in tera hashes per second. (d) Average market price in US Dollar. (e) Total mining reward in Bitcoin currency divided into block reward (blue line) and transaction fees (orange line).

Fig. 2 .

 2 Fig.2. Ribbon chart of market share distribution according to reward rules. Each bar represents the market share of a mining pool while its color stands for a specific reward. The stacked bars are sorted monthly by the total market share of each reward rule and by the market share of the mining pools within each group. The legend on the right is sorted by the market share over the entire period. The red horizontal line indicates a 50% mining power threshold.

Fig. 3 .

 3 Fig.3. Ribbon chart of market share distribution according to mining pool locations. Each bar represents the market share of a mining pool while its color stands for a location. The stacked bars are sorted monthly by the total market share of each location and by the market share of the mining pools within each group. The legend is sorted by the market share of each location over the entire period. The red horizontal line indicates a 50% mining power threshold.

https://en.bitcoin.it/wiki/Comparison of mining pools

e. g., https://www.cryptocompare.com/mining/#/pools and https://www. bitcoinmarketjournal.com/best-bitcoin-mining-pool/