Javier González-Delgado

Amin Sagar

Christophe Zanon

Kresten Lindorff-Larsen

Pau Bernadó

Pierre Neuvial

Juan Cortés

WASCO: A Wasserstein-based statistical tool to compare conformational ensembles of intrinsically disordered proteins *

The structural investigation of intrinsically disordered proteins (IDPs) requires ensemble models describing the diversity of the conformational states of the molecule. Due to their probabilistic nature, there is a need for new paradigms that understand and treat IDPs from a purely statistical point of view, considering their conformational ensembles as well-defined probability distributions. In this work, we define a conformational ensemble as an ordered set of probability distributions and provide a suitable metric to detect differences between two given ensembles at the residue level, both locally and globally. The underlying geometry of the conformational space is properly integrated, one ensemble being characterized by a set of probability distributions supported on the three-dimensional Euclidean space (for global-scale comparisons) and on the two-dimensional flat torus (for local-scale comparisons). The inherent uncertainty of the data is also taken into account to provide finer estimations of the differences between ensembles. Additionally, an overall distance between ensembles is defined from the differences at the residue level. We illustrate the interest of the approach with several examples of applications for the comparison of conformational ensembles: (i) produced from molecular dynamics (MD) simulations using different force fields, and (ii) before and after refinement with experimental data. We also show the usefulness of the method to assess the convergence of MD simulations, and discuss other potential applications such as in machine-learning-based approaches.

Introduction

The comparison of protein structures is a crucial problem in structural biology. In the early works [START_REF] Rao | Comparison of super-secondary structures in proteins[END_REF][START_REF] Maiorov | Significance of root-mean-square deviation in comparing three-dimensional structures of globular proteins[END_REF], the use of root-mean-square deviation (RMSD) was introduced and discussed as a metric between conformations of folded proteins, and later extended to its ensemble version [START_REF] Brüschweiler | Efficient RMSD measures for the comparison of two molecular ensembles[END_REF]. More recently, Lindorff-Larsen and Ferkinghoff-Borg [START_REF] Lindorff-Larsen | Similarity measures for protein ensembles[END_REF] defined three metrics that allow overall comparison between ensembles of ordered/structured systems, with stronger mathematical guarantees, but using RMSD as a distance between individual conformations, which complicates its extension to disordered structures. Cazals et al. [START_REF] Cazals | Conformational ensembles and sampled energy landscapes: Analysis and comparison[END_REF] used a graph-based representation of the conformational space based on a set of low-energy conformations (i.e. local minima of the potential energy landscape) and compared them with the more suitable Wasserstein distance. To do so, they used the least-RMSD as ground metric between conformations. The methods presented in [START_REF] Lindorff-Larsen | Similarity measures for protein ensembles[END_REF] and [START_REF] Cazals | Conformational ensembles and sampled energy landscapes: Analysis and comparison[END_REF] are well suited to examine conformational ensembles of molecules that present a well-characterized energy landscape. However, their application to molecules with energy landscapes where low-energy conformations are difficult to identify, as it is the case of IDPs, is inappropriate.

A few recent works have dealt with the comparison of conformational ensembles of IDPs. Huihui and Ghosh [START_REF] Huihui | Intrachain interaction topology can identify functionally similar intrinsically disordered proteins[END_REF] focused on averaged conformational properties over ensembles as informative descriptors of their function. They proposed a sequence-decoration metric that classifies IDPs using only their primary structure together with their charge configuration. The same idea of comparing average descriptors was applied by Lazar et al. [START_REF] Lazar | Distancebased metrics for comparing conformational ensembles of intrinsically disordered proteins[END_REF], who proposed an ensemble comparison tool based on differences between average pairwise distances. Due to the huge conformational variability of IDPs, it is, however, important to take into account both the average properties as well as the distribution around those averages. Describing IDP conformations as being drawn from probability distributions determining their structure may yield to an important loss of information (or even misleading results) if the whole distribution is reduced to its mean. Even when comparing two (possibly multivariate) Gaussian distributions, the difference between the two depends both on the means and variances [START_REF] Kullback | An application of information theory to multivariate analysis[END_REF][START_REF] Zhou | From sample similarity to ensemble similarity: Probabilistic distance measures in reproducing kernel Hilbert space[END_REF]; thus, methods for comparing ensembles should ideally include also higher order moments of the probability distributions. This is why a statistical approach that integrates the entire probability law defining an ensemble is crucial to correctly capture the existing differences between disordered ensembles.

The probability distributions describing the ensembles need to be compared using a suitable metric, well-adapted to the geometric features of the underlying spaces. The Wasserstein distance [START_REF] Villani | Optimal Transport: Old and New[END_REF], sometimes called "earth mover's distance", integrates the geometry of the space where the distributions are supported and provides strong mathematical guarantees. Moreover, it has a physical interpretation, as it is defined as the minimum transportation cost needed to reconfigure the mass of one probability distribution to recover the other. All this makes Wasserstein distance substantially preferable to other metrics currently used in the literature (e.g. Kullback-Leibler divergence, Helliger distance), as discussed in Section 2.

In this work, we define a set of probability distributions that characterize at local and global level the highly variable conformations in an ensemble of disordered proteins, and to which we can have access in practice. These probability laws can then be compared using the Wasserstein distance, allowing the identification of residue-specific and overall discrepancies. We also propose an approach to integrate the intrinsic uncertainty of the data within the metric, which enables a more clear identification of the relevant differences between the ensembles. The method has been implemented inside a purely non-parametric framework, avoiding model assumptions, dimensionality reduction or further simplifications that may yield significant loss of information.

In the following sections, we provide an overall description of the proposed methodology, which is further detailed in the Supplementary Information (SI), together with several cases of applications that illustrate how our method identifies residue-specific and overall discrepancies between conformational ensembles of IDPs or flexible peptides generated for example by molecular dynamics simulations or stochastic sampling techniques. Finally, we discuss current limitations and possible extensions of WASCO, as well as the great potential interest of this type of metric for its integration in machine-learning-based (ML-based) methods applied to generate or to refine conformational ensembles of IDPs.

Methods

Due to the intrinsic probabilistic nature of IDPs, descriptors of their conformational ensembles should be conceived from a purely statistical point of view. To do so, we seek to locally and globally describe conformational ensembles using well-defined probability distributions and to develop statistical tools allowing their comparison. The main questions to answer are therefore: (1) which is the best way to define those probability distributions? and (2) how these distributions have to be compared to provide quantitative information about similarities and differences between ensembles?

Defining conformational ensembles as a set of probability distributions

IDP ensembles can be described at both local and global scales, providing complementary information. We aim at defining an ordered set of probability distributions that account for the highly variable structure of the ensemble and, above all, that can be estimated in practice from a set of sampled conformations.

The most important aspects of the local structure can be described by the dihedral angles (φ, ψ) for each amino acid residue along the sequence. Therefore, for each residue, the ensemble is locally characterized by a two-dimensional random variable (φ, ψ) or, in other words, by a probability distribution supported on the two-dimensional flat torus T 2 [START_REF] Mardia | Protein bioinformatics and mixtures of bivariate von Mises distributions for angular data[END_REF][START_REF] Boomsma | A generative, probabilistic model of local protein structure[END_REF]. If we denote such distribution as P l i , for the residue at the i-th position, we define the local structural descriptor of an ensemble as the L-tuple

(P l 1 , . . . , P l L ), P l i ∈ P(T 2 ) for all i = 1, . . . , L, ( 1
)
where L is the sequence length and P(T 2 ) denotes the space of probability distributions supported on

T 2 .
Describing the global structure is a less trivial task. The use the absolute positions of the atoms and an absolute reference frame for the entire ensemble is not an appropriate description as it is sensitive to rigid-body motions. Therefore, our approach uses the relative positions of all pairs of residues along the sequence, which are invariant under rigid-body motion. More precisely, we define the position of a given residue as the position of its C β atom when it exists and of its C α atom otherwise. If i, j ∈ {1, . . . , L}, i ̸ = j, denote two different sequence positions, let --→ R i,j be the three-dimensional random variable determining the relative position of j-th residue with respect to the i-th one. If we denote P g i,j

the probability distribution associated to --→ R i,j , we define the global structural descriptor of an ensemble as the (L(L -1)/2)-tuple (P g 1,2 , P g 1,3 , . . . , P g L-1,L ), P g i,j ∈ P(R 3 ) for all i = 1, . . . , L -1,

j = i + 1, . . . , L, ( 2
)
where L is the sequence length and P(R 3 ) denotes the space of probability distributions supported on the three-dimensional Euclidean space.

Accessing empirical probability distributions from sampled conformations

Estimating the local structural descriptor (1) is immediate as we have direct access to dihedral angles (φ, ψ) from the sample of conformations. Therefore, the local structural descriptor will be estimated by its empirical counterpart (P l 1;n , . . . , P l L;n ),

where each P l i,n , i = 1, . . . , L, is the empirical probability distribution of P l i , and n is the number of conformations constituting the sample. Such empirical probability distributions are commonly represented through Ramachandran maps [START_REF] Ramachandran | Stereochemistry of polypeptide chain configurations[END_REF].

Obtaining a sample of --→ R i,j from the set of conformations is less direct. To compute a set of comparable --→ R i,j vectors from all conformations, their coordinates must be expressed on the same reference system. To do so, we first define a reference frame at the i-th residue, using only the positions of the i-th C ′ , C α i, j along the sequence. In (a), the reference frame is built for every conformation at residue i. In (b), all the frames are superimposed using this reference frame. Then, for any j ̸ = i, the vectors -→ r i,j constitute a sample of --→ R i,j .

and N H atoms. This frame, whose construction is detailed in the Supplementary Information (SI), is a meaningful representation of the spatial pose of each residue.

The reference frame associated to each residue i ∈ {1, . . . , L} allows to express the relative positions of all residues j ̸ = i with respect to i. Moreover, the definition of a reference system allows the superposition of all the conformations in the ensemble. This is illustrated in Figure 1, for three conformations. Consequently, for every j ̸ = i, we will have access to n realizations of the random variable --→ R i,j or, in other words, to a point cloud in the three-dimensional Euclidean space, representing a sample drawn from the distribution of P g i,j . Therefore, the global structural descriptor of the ensemble (2) will be estimated by its empirical counterpart (P g 1,2;n , P g 1,3;n , . . . ,

P g L-1,L;n ), (4)
where P g i,j,n is the empirical probability distribution of P g i,j , for all i = 1, . . . , L -1, j = i + 1, . . . , L. An example of a pair of samples of --→ R i,j is presented in Figure S8.

Distances between local and global structural descriptors

After defining the local and global structural descriptors of an ensemble as an ordered set of probability distributions, the choice of a suitable metric allowing inter-ensemble comparisons becomes the subsequent question to deal with. The basic properties that such a metric should have are:

1. Satisfying the mathematical properties that define a distance (i.e. being 0 if an only if the two compared distributions are identical, symmetry and triangle inequality), 2. Integrating the geometry of the underlying space.

The use of metrics between probability distributions is not new in structural biology. For instance, Ting et al. [START_REF] Ting | Neighbordependent ramachandran probability distributions of amino acids developed from a hierarchical dirichlet process model[END_REF] used Hellinger distance to detect differences between (φ, ψ) distributions. However, this metric does not take into account the geometry of the underlying space (in particular here, its periodicity). A symmetrized Kullback-Leibler (KL) or the Jensen-Shannon (JS) divergence was used in [START_REF] Lindorff-Larsen | Similarity measures for protein ensembles[END_REF][START_REF] Tiberti | ENCORE: software for quantitative ensemble comparison[END_REF] to compare ensembles of ordered systems. This metric has a firm interpretation, based on information theory (in particular the JS divergence is the square of a metric). However, it still misses the geometrical reliability and does not satisfy triangle inequality, which makes comparisons between multiple ensembles difficult to interpret.

Besides satisfying conditions 1 and 2, the Wasserstein distance, derived from the theory of Optimal Transport (OT), provides both strong theoretical guarantees [START_REF] Villani | Optimal Transport: Old and New[END_REF] and attractive empirical performance [START_REF] Peyré | Computational optimal transport: With applications to data science[END_REF]. Informally, it represents the minimum transportation cost needed to reconfigure the mass of one probability distribution to recover the other. A more technical definition is provided in SI, and we refer to [START_REF] Peyré | Computational optimal transport: With applications to data science[END_REF] for an in-depth introduction to OT. Most of the applications of OT are related to the very active field of machine learning (ML), notably in the framework of generative networks [START_REF] Arjovsky | Wasserstein generative adversarial networks[END_REF], robustness [START_REF] Serrurier | Achieving robustness in classification using optimal transport with hinge regularization[END_REF] or fairness [START_REF] Del Barrio | A central limit theorem for lp transportation cost on the real line with application to fairness assessment in machine learning[END_REF], among others. With some notable exemptions [START_REF] Cazals | Conformational ensembles and sampled energy landscapes: Analysis and comparison[END_REF][START_REF] Berg | Towards a molecular basis of ubiquitin signaling: A dual-scale simulation study of ubiquitin dimers[END_REF][START_REF] Rosenbaum | Inferring a continuous distribution of atom coordinates from cryo-em images using vaes[END_REF][START_REF] Damjanovic | Catboss: Cluster analysis of trajectories based on segment splitting[END_REF][START_REF] González-Delgado | Two-sample goodnessof-fit tests on the flat torus based on wasserstein distance and their relevance to structural biology[END_REF], Wasserstein distance has not been widely used in structural biology. More related to our work, in [START_REF] Cazals | Conformational ensembles and sampled energy landscapes: Analysis and comparison[END_REF], Cazals et al. used Wasserstein distance to compare energy landscapes sampled from conformational ensembles. Recently, it was used in [START_REF] González-Delgado | Two-sample goodnessof-fit tests on the flat torus based on wasserstein distance and their relevance to structural biology[END_REF] to define statistical tests assessing differences between (ϕ, ψ) distributions. The incorporation of the underlying geometry to its definition makes it a well-adapted metric to measure distances between local and global structural descriptors of the ensembles. Details and important considerations regarding its practical computation are given in SI.

The comparison tool

Consider two ensembles A, B, associated to two protein sequences of equal length L, and let n A , n B be their number of conformations, respectively. We define the differences between local structural descriptors of A and B as the L-tuple of Wasserstein distances

(W l,A,B 1 , . . . , W l,A,B L ) = W(P l,A 1;n A , P l,B 1;n B ), . . . , W(P l,A L;n A , P l,B L;n B ) , (5)
where P l,A i;n A (resp. P l,B i;n B ) denotes the i-th distribution of the empirical local structural descriptor (3) of ensemble A (resp. B). Statistical tests to assess whether any W l,A,B i is significantly different from zero have been recently defined in [START_REF] González-Delgado | Two-sample goodnessof-fit tests on the flat torus based on wasserstein distance and their relevance to structural biology[END_REF]. The second of the introduced techniques is better adapted to our problem, as it only detects the more important discrepancies and accepts slight differences that may arise from experimental or computational procedures. This is discussed in detail in [START_REF] González-Delgado | Two-sample goodnessof-fit tests on the flat torus based on wasserstein distance and their relevance to structural biology[END_REF]. Consequently, together with the L-tuple (5) of distances comparing local structural descriptors, we are able to supply a L-tuple of p-values (corrected for multiplicity [START_REF] Holm | A simple sequentially rejective multiple test procedure[END_REF]) accounting for the statistical significance of the corresponding distances:

(p A,B 1 , . . . , p A,B L ). (6)
Recall that a small p-value p A,B i indicates strong evidence that the true distance that W l,A,B i estimates is different from zero. In other words, small p-values show significant differences between the corresponding local structural descriptors. Therefore, the vector (6) enables the identification of those residues where the differences are more important, and those residues for which differences can be assigned as non-significant.

Analogously, the difference between global structural descriptors of A and B is defined as the (L(L -

1)/2)- tuple (W g,A,B 1,2 , . . . , W g,A,B L-1,L ) = W(P g,A 1,2;n A , P g,B 1,2;n B ), . . . , W(P g,A L-1,L;n A , P g,B L-1,L;n B ) , ( 7
)
where P g,A i,j;n A (resp. P g,B i,j;n B ) denotes the i, j distribution of the empirical global structural descriptor (4) of ensemble A (resp. B). In this case, we are not able to provide a vector of p-values assessing the significance of the global differences. This is due to the intrinsic limitations of the underlying mathematical theory when the ground space has dimension d ≥ 3. Note that (7) can be more naturally represented as a triangular (L -1) × (L -1) matrix W g,A,B , whose elements are given by (W g,A,B ) ij = W g,A,B i,j . Graphically, the matrix W g,A,B is represented using a color scale to fill the coefficients according to distance values. As the diagonal will remain empty, it will be filled with the local distances [START_REF] Cazals | Conformational ensembles and sampled energy landscapes: Analysis and comparison[END_REF]. This will also allow to assess whether changes on local structural descriptors are related with changes in global structural descriptors and to compare both scales within the same representation.

Accounting for uncertainty

The variability in experimental and simulated structures causes uncertainties and statistical noise that may substantially bias the distance estimation. For example, when running a MD simulation, independent replicas of the same simulation setup may results in non-negligible differences that distort the analysis of the comparison matrices. The same may occur when comparing two uniformly chosen subsets of conformations from an ensemble generated by stochastic sampling techniques [START_REF] Ozenne | Flexible-meccano: a tool for the generation of explicit ensemble descriptions of intrinsically disordered proteins and their associated experimental observables[END_REF][START_REF] Estaña | Realistic ensemble models of intrinsically disordered proteins using a structure-encoding coil database[END_REF]. In order to soften the effect of uncertainty and to obtain net estimates of the differences between a pair of ensembles, we will use (if available) independent replicas of the same ensemble. These replicas may also be produced by uniform subsampling of the set of conformations. However, special care must be taken when subsampling MD trajectories as the convergence of the simulation must be ensured for the subsamples to be representative of the entire ensemble.

Let A 1 , . . . , A n I (resp. B 1 , . . . , B n I ) be n I independent replicas of ensemble A (resp. B). The corrected difference between local structural descriptors of A and B is defined as the L-tuple

( W l,A,B 1 , . . . , W l,A,B L ), (8)
where each corrected distance is defined as

W l,A,B i = 1 n I n I s=1 W l,As,Bs i - 1 2(n I -1) n I s=2 W l,A1,As i + W l,B1,Bs i +
, for all i = 1, . . . , L, (

where, for any real number x, (x) + = x if x > 0 and (x) + = 0 otherwise. The first term in ( 9) is an average of n I Wasserstein distances between n I paired independent replicas of A and B. As it was shown in [START_REF] Sommerfeld | Optimal transport: Fast probabilistic approximation with exact solvers[END_REF], an average of Wasserstein distances between sub-samples of the same population is a pertinent estimate of the Wasserstein distance between the two entire populations that, in addition, conserves the properties that mathematically define a distance. Therefore, this first term estimates the Wasserstein distance between the entire populations of A and B (conceived as the union of all independent replicas), softening the variability. To this brutto inter-ensemble difference, we subtract an average of the Wasserstein distances between independent replicas of the same population (intra-ensemble). Note that, for the sake of computational simplicity, we just compared the first replica of each ensemble with the subsequent ones. This alignment is arbitrary and can be set otherwise. Of course, distances between all pairs of replicas can be added to this term. The more combinations are added to [START_REF] Zhou | From sample similarity to ensemble similarity: Probabilistic distance measures in reproducing kernel Hilbert space[END_REF], the finer will be the estimate of the (unknown) true Wasserstein distance between the ensembles but, as replicas are independent, different alignments for a given number of combinations should not yield substantial discrepancies on the quality of this estimate. The same applies if n I is different for A and B; both terms in (9) can be accordingly adapted. As it is illustrated in Section 3, the use of corrected distances (9) contribute to reduce the noise coming from structural uncertainty and help to emphasize residue-specific differences in the matrix representation. For the distances between global structural descriptors, the correction is performed analogously.

Setting an interpretable scale

When defining an absolute distance or score between conformational ensembles, providing the clues to ease its interpretation is crucial. The problem of interpreting unbounded metrics with no intrinsic reference values has been widely discussed since the introduction of RMSD for the comparison of pairs of conformations [START_REF] Rao | Comparison of super-secondary structures in proteins[END_REF][START_REF] Maiorov | Significance of root-mean-square deviation in comparing three-dimensional structures of globular proteins[END_REF]. Here, we do not seek to define any cutoff to binarize the resulting matrices, but to provide a more informative continuous scale. To do so, we aim at quantifying the magnitude of the inter-ensemble distances compared to the intra-ensemble ones, using the uncertainty estimate as a reference. If we denote as W l,A,B inter (resp. W l,A,B intra ) the first (resp. second) term in ( 9), the score

W l,A,B i W l,A,B intra = W l,A,B inter -W l,A,B intra + W l,A,B intra , ( 10
)
corresponds to the relative difference between the inter-ensemble and intra-ensemble differences. Once again, this score is analogously defined for differences between global structural descriptors.

An overall distance between ensembles

In some situations, it may be of interest to perform overall comparisons between multiple ensembles.

To do so, moving from a residue-specific analysis to a comparison at the whole structure level might be preferable. The definition of a score for the overall ensemble has been addressed for ordered systems [START_REF] Lindorff-Larsen | Similarity measures for protein ensembles[END_REF].

Here, we propose to define such a score by aggregating all the residue-specific distances computed using the above-described methods. We recall that if d 1 , . . . , d L are L distances defined on L metric spaces X 1 , . . . , X L , the function

d 2 1 + • • • + d 2 L is a distance on the product space X 1 × • • • × X L . Consequently, OW l,A,B = L i=1 W l,A,B i 2 1/2 (11)
is a distance on the product space of all dihedral angles along the sequence and, therefore, serves to quantify the overall local discrepancy between a pair of ensembles. Analogously,

OW g,A,B =   L-1 i=1 L j=i+1 w ij W g,A,B i,j 2   1/2
, with w ij > 0 for all i, j ∈ {1, . . . , L}, (

is a distance on the product space of all pairwise relative positions of the residues in both ensembles, and serves to quantify the overall global discrepancy. Note that we have assigned a positive weight w ij to each global distance in [START_REF] Boomsma | A generative, probabilistic model of local protein structure[END_REF]. This allows to consider distances between specific residue pairs as more relevant than the others when computing the overall discrepancy [START_REF] Holm | Protein structure comparison by alignment of distance matrices[END_REF]. For instance, we can highlight differences between global structural descriptors that appear for residue pairs that are far from each other in the sequence, i.e. large |i -j|, and neglect distances between neighboring residue pairs, i.e small |i -j|. This can be done by choosing w ij as an appropriate increasing function of |i -j|, as

w ij = w(i, j) = 1 tanh 1 tanh |i -j| L -1 1 2 , ( 13
)
which satisfies w i,i = 0 for all i and w 1,L = w L,1 = 1.

The drawback of this definition of the overall distance is that it does not take into account the uncertainty discussed in Section 2.4.1. To solve this problem, the same strategy to define a global score can be performed by replacing each W l,A,B i (resp. W g,A,B i,j

) by its corresponding corrected distance W l,A,B i (resp. W g,A,B i,j

) in (11) (resp. ( 12)). However, this strategy makes the triangle inequality for the overall metric no longer satisfied. Both scores can be implemented by the practitioner and used depending on the specific comparison context.

The Jupyter Notebook

The WASCO comparison tool has been implemented through an easy-to-use Jupyter Notebook. It is available at https://gitlab.laas.fr/moma/WASCO, together with its installation guidelines and detailed implementation instructions. The notebook takes a pair of ensembles as input and returns the comparison results through the matrix defined in Section 2.4, containing global and local differences. Users can choose to correct the computed distances by uncertainty, as proposed in Section 2.4.1. When independent replicas are not provided as input, subsampling is used to emulate them. If this correction is performed, results are displayed in the interpretable scale defined in Section 2.4.2. The overall score defined in Section 2.4.3, aggregating the corrected distances, is also returned by the tool.

Ensembles can be provided as input in several of the most common data formats. WASCO accepts one .xtc file per replica, together with a .pdb file including the topology information of the molecule, one multiframe .pdb file per replica or a folder per replica containing one .pdb file per conformation. The user can also choose to compare ensembles for sequence segments (of equal length) instead of the entire sequence. Details are provided in the notebook documentation.

Due to the large number of Wasserstein distances to be computed (L(L -1)/2 + L per pair of replicas), the computation time might be considerably high. The number of conformations constituting the ensemble also has a significant impact, due to computational limitations of the existing OT algorithms when sample sizes and dimension increase. In order to return results within a reasonable amount of time, WASCO computes Wasserstein distances in parallel. The required CPU time depends on the number of conformations, replicas and sequence length of the ensembles. For small proteins of L ∼ 30 and ensembles of reasonable size n A , n B ∼ 10 4 , the CPU time using 20 threads is less than 15 minutes using a standard computing server. However, for larger proteins of L ∼ 150 and large ensembles with n A , n B ∼ 10 5 , the CPU time using 20 threads goes up to some hours. Additionally, comparing large ensembles of substantially longer sequences (L ≫ 150) might cause memory problems, as all pairwise relative positions for every conformation need to be stocked. Therefore, the suitability of the sizes of the ensembles must be considered before launching WASCO. Adapting WASCO to longer sequences with large conformational ensembles remains an objective for future work.

The output of WASCO is given through a matrix, whose entries are the values of the score [START_REF] Villani | Optimal Transport: Old and New[END_REF] computed for local and global distances, when independent replicas are provided. Otherwise, the matrix depicts the values of the non-corrected inter-ensemble distances ( 5), [START_REF] Lazar | Distancebased metrics for comparing conformational ensembles of intrinsically disordered proteins[END_REF]. The values for the discrepancies between the global structural descriptors (4) are provided in the lower triangle. The differences between the local structural descriptors (3) are displayed along the diagonal. Details on the interpretation of the matrix are given in Section S1.3 and illustrated in Figure S4. These guidelines are also presented in the software documentation.

Results

In this section, we present several applications to illustrate the different possibilities enabled by WASCO. In all cases, the distances between local and global structural descriptors were corrected for uncertainty using (9), as independent replicas were available. The results are depicted through the score [START_REF] Villani | Optimal Transport: Old and New[END_REF], representing the relative difference between the inter-ensemble distances and the uncertainty. Both overall local and global discrepancies between pairs of ensembles were computed plugging the corrected distances in [START_REF] Mardia | Protein bioinformatics and mixtures of bivariate von Mises distributions for angular data[END_REF] and [START_REF] Boomsma | A generative, probabilistic model of local protein structure[END_REF], as discussed in Section 2.4.3. The weight function [START_REF] Ramachandran | Stereochemistry of polypeptide chain configurations[END_REF] was used to highlight differences between residue pairs far from each other in the sequence and reduce differences between neighboring amino acids. Note that this weighting is considered only to compute the overall distance [START_REF] Boomsma | A generative, probabilistic model of local protein structure[END_REF], and not to depict distance values in the matrix representation, which correspond to the interpretable scale defined in Section 2.4.2. An additional analysis illustrating the application of WASCO to assess the convergence of MD simulations is included in SI.

Comparison of ensembles produced by MD simulations using different force-fields

We applied WASCO to compare the results of MD simulations using different force-fields presented in [START_REF] Jephthah | Force field effects in simulations of flexible peptides with varying polyproline II propensity[END_REF] for two flexible peptides showing a significant propensity to form poly-l-proline type II (PPII) structures. Four different force-fields, having demonstrated relatively good performances to simulate IDPs were applied: AMBER ff99SB-disp, AMBER ff99SB-ILDN, CHARMM36IDPFF, and CHARMM36m (details and references to these force-fields can be found in [START_REF] Jephthah | Force field effects in simulations of flexible peptides with varying polyproline II propensity[END_REF]). For simplicity, we will refer to these force-fields as disp, ildn, c36idp and c36m, respectively. As independent replicas for each simulation were available, we could perform the correction for uncertainty [START_REF] Zhou | From sample similarity to ensemble similarity: Probabilistic distance measures in reproducing kernel Hilbert space[END_REF]. Figure 2 presents the output of WASCO for several pairwise comparisons of conformational ensembles of Histatin-5 (Hst5) obtained with the different force-fields. The matrices and the overall dissimilarities suggest that the generated structures are closer (in Wasserstein distance) when they are simulated using c36idp and c36m (which we can define as group-I), or disp and ildn (group-II). This is not surprising as group-I are versions of CHARMM and group-II are versions of AMBER. Indeed, matrices (a) and (b), comparing force fields inside group-I and inside group-II respectively, present overall global differences which are small compared to those of panels (c) and (d), which compare force-fields of different groups. The same conclusion can be reached by comparing the magnitude of the scales of both pairs of matrices. The two remaining comparisons (ildn vs. ildn and c36m vs. disp) are not included in Figure 2 as the corresponding matrices are qualitatively equivalent to (c) and (d). Similar observations have been made when comparing ensembles of folded proteins generated using related force-fields [START_REF] Tiberti | ENCORE: software for quantitative ensemble comparison[END_REF][START_REF] Martín-García | Comparing molecular dynamics force fields in the essential subspace[END_REF].

Matrices returned by WASCO also allow a residue-specific analysis of the distances. In Figure 2, panels (c) and (d) show that the most relevant global differences appear in regions close to the diagonal (i.e. between residue pairs close in the sequence), where the inter-ensemble corrected distances rise up to 6-7 times the intra-ensemble ones. This is not the case when comparing force-fields inside the same group, as the largest differences appear in more internal matrix regions (i.e. between residue pairs more distant in the sequence). However, these corrected differences represent less than the half of the intraensemble distances. The information displayed on the diagonal allows the detection of the residues where the local conformation change more abruptly between force-fields. These local changes are restricted to smaller regions, contrary to the observed behaviour of global differences, which appeared for more extent regions inside the lower triangle and not for isolated pairs of amino acids. In some cases, substantial local distances appear in residues where global structure also changes (see, for example, residues next to the N-terminus in (a,c)). However, this correspondence is not observed in all matrices. We repeated the same analysis for MD simulations of PEP3 with the same force-fields. Results are presented in SI.

Structural impact of SAXS ensemble refinement

Using Hst5 as an example, we applied WASCO to evaluate the structural impact of SAXS refinement with the Ensemble Optimization Method (EOM) [START_REF] Bernadó | Structural characterization of flexible proteins using small-angle X-ray scattering[END_REF] on the resulting ensemble. We first compared the Hst5 ensemble simulated with Flexible-Meccano [START_REF] Ozenne | Flexible-meccano: a tool for the generation of explicit ensemble descriptions of intrinsically disordered proteins and their associated experimental observables[END_REF][START_REF] Bernadó | A structural model for unfolded proteins from residual dipolar couplings and small-angle x-ray scattering[END_REF] with the refined one using previously reported SAXS data [START_REF] Sagar | Comment on the optimal parameters to derive intrinsically disordered protein conformational ensembles from smallangle X-ray scattering data using the ensemble optimization method[END_REF]. The results are presented in Figure 3. Note that a previous EOM analysis of these data suggested that Hst5 in solution is slightly more extended that the random coil model generated with Flexible-Meccano [START_REF] Sagar | Comment on the optimal parameters to derive intrinsically disordered protein conformational ensembles from smallangle X-ray scattering data using the ensemble optimization method[END_REF]. Small but significant differences were observed at the central part of the peptide (from residues 6 to 13). Most probably, the SAXS-based refinement selected conformations with an extended central region to account for the overall expansion of peptide in the solution [START_REF] Bernadó | Structural characterization of flexible proteins using small-angle X-ray scattering[END_REF]. Moreover, we observed highly significant local distances that propagate towards the interior of the matrix. In other words, these residues with large local distances conformationally influence their closest neighbours. Intriguingly, this propagation seems to only occur towards the C-terminus.

We next assessed whether the direction in which conformations are built have a structural effect and change the refined ensemble. To do so, we generated two Hst5 ensembles using a stochastic sampling method similar to Flexible-Meccano but using a different strategy [START_REF] Estaña | Realistic ensemble models of intrinsically disordered proteins using a structure-encoding coil database[END_REF], where the chains were built either from N-to-C or from C-to-N. When using these two ensembles to fit the experimental curve, the resulting distance matrices displayed very similar features for local and global distances (Figures S9a andS9b), suggesting that the chain-building direction does not have a relevant effect. In both cases, a systematic increase in the distances is observed for the central residues, as observed in the previous analysis (Figure 3).

In a recent study, ENCORE was used to show that refined ensembles were closer to each other than different input ensembles [START_REF] Tiberti | ENCORE: software for quantitative ensemble comparison[END_REF]. This can also be illustrated using WASCO, by comparing the Hst5 ensembles generated in both directions before and after the filtering with SAXS data (Figures S9c andS9d). These comparisons clearly showed that both global and local differences were smaller for the refined ensembles than for the input ones, as observed when comparing the maximum values of the corresponding color scales. As we were comparing very similar ensembles, we expected the distances to be small. Nevertheless, we observe one significant local difference on the diagonal in Figure S9c that disappeared after refinement.

Discussion

We have presented a novel method to compare conformational ensemble models of highly flexible proteins. WASCO is based on a non-parametric framework: local and global structural descriptors of the conformational space are defined as distributions and do not rely on probabilistic or statistic models. This allows capturing the entire variability of the ensemble without information loss. The distributions are compared using the Wasserstein distance, which has strong mathematical guarantees and respects the geometry of the underlying space. To this metric, we incorporated the structural uncertainty presented in experimental and simulated ensembles. Using this strategy, WASCO highlights the relevant differences between ensembles. We have illustrated several possible applications of WASCO as an additional tool for the investigation of IDPs and flexible peptides. It provides complementary information with respect to other tools to analyze and compare conformational ensembles based on global descriptors, such as the radius of gyration [START_REF] Tria | Advanced ensemble modelling of flexible macromolecules using X-ray solution scattering[END_REF] or secondary structure propensities [START_REF] Jephthah | Force field effects in simulations of flexible peptides with varying polyproline II propensity[END_REF]. Besides, the presented approach is advantageous with respect to simpler comparison techniques based on average descriptors, such as the difference of median distance matrices introduced in [START_REF] Lazar | Distancebased metrics for comparing conformational ensembles of intrinsically disordered proteins[END_REF]. This is illustrated with an example in SI (Section S2.3). Thanks to its accuracy to identify differences between ensembles, WASCO has great potential interest for integration into ML-based methods for generating or refining conformational ensembles of IDPs [START_REF] Lindorff-Larsen | On the potential of machine learning to examine the relationship between sequence, structure, dynamics and function of intrinsically disordered proteins[END_REF][START_REF] Janson | Direct generation of protein conformational ensembles via machine learning[END_REF][START_REF] Zhang | Learning to evolve structural ensembles of unfolded and disordered proteins using experimental solution data[END_REF]. More precisely, metrics based on WASCO can be used to evaluate the performance of these methods, or as a loss function when training neural network models.

WASCO has been implemented in an open-source Jupyter Notebook, which enables an easy use of the methods as well as their adaptation or extension to particular needs. The main drawback of the current implementation is its limitation to deal with considerably large ensembles of long IDPs. Adapting WASCO to larger chains remains for future work. Other interesting directions for future work will be the extension of WASCO to compare ensembles of multi-domain proteins, and to operate with coarsegrained models. The extension of WASCO to compare ensembles for chains of different length is also an interesting but challenging work. Note however that the Jupyter Notebook enables the user to select sequence fragments of equal length for the comparison.

Figure 1 :

 1 Figure 1: Illustration of how samples of global structural descriptors are obtained, for a pair of positionsi, j along the sequence. In (a), the reference frame is built for every conformation at residue i. In (b), all the frames are superimposed using this reference frame. Then, for any j ̸ = i, the vectors -→ r i,j constitute a sample of --→ R i,j .

Figure 2 :

 2 Figure 2: Comparison of Molecular Dynamics simulations of Hst5 ensemble using different force-fields. The color scale W /W intra corresponds to the score (10), representing the relative difference between the inter-ensemble distances and the uncertainty. The coefficients in the lower-triangle (in red) correspond to the global differences. The coefficients along the diagonal (in blue) correspond to the local differences. Blue stars indicate that the corresponding local corrected distance is significantly different from zero (the associated p-value (6) is smaller than α = 0.05). Note the different scales used in the different plots.

Figure 3 :

 3 Figure 3: Comparison of Hst5 ensemble before and after filtering with experimental SAXS data. The color scale W /W intra corresponds to the score (10), representing the relative difference between the inter-ensemble distances and the uncertainty. The coefficients in the lower triangle (in red) correspond to the global differences. The coefficients along the diagonal (in blue) correspond to the local differences. Blue stars indicate that the corresponding local corrected distance is significantly different from zero (the associated p-value (6) is smaller than α = 0.05).

Acknowledgements

We are grateful to Francesco Pesce, Sthitadhi Maiti and Matthias Heyden for providing useful data. We thank Gabriella Gerlach, Frederik Emil Thomasen and Philipp Schake for their helpful discussions and valuable feedback on WASCO implementation. This work has been partially supported by the French National Research Agency (ANR) through grant ANR-19-P3IA-0004, the LabEx CIMI (ANR-11-LABX-0040) and EpiGenMed (ANR-10-LABX-12-01) within the French State Programme "Investissements d'Avenir", by the European Research Council under the European Union's H2020 Framework Programme (2014-2020)/ ERC Grant agreement n°[648030] awarded to PB and by the Lundbeck Foundation BRAINSTRUC initiative (R155-2015-2666). The CBS is a member of France-BioImaging (FBI) and the French Infrastructure for Integrated Structural Biology (FRISBI), 2 national infrastructures supported by the French National Research Agency (ANR-10-INBS-04-01 and ANR-10-INBS-05, respectively).


Robin P M Gauff
email: gauff.robin@yahoo.de

Christophe Lejeusne

Stéphane Greff

Stéphane Loisel

Olivier Bohner

Dominique Davoult

Getting used to the heat -Success of non-indigenous species facing climate change might depend on prior climate context

Keywords: Introduced species, Climate Change, in situ warming, marina, Hot Plates, Metabolomics

Climate change constitutes a major challenge for marine urban ecosystems and ocean warming will likely strongly affect local communities. Non-Indigenous Species (NIS) have been shown to often have higher heat resistance than natives, but studies investigating how forthcoming global warming might affect them in marine urban environments remain scarce, especially in situ studies.

Here we used an in situ warming experiment in a NW Mediterranean (warm temperate) and a NE Atlantic (cold temperate) marina to see how global warming might affect recruited communities in the near future. In both marinas, warming resulted in significantly different community structure, lower biomass, and more empty space compared to control. However, while in the warm temperate marina, NIS showed an increased surface cover, it was reduced in the cold temperate one. Metabolomic analyses in the Atlantic marina revealed potential heat stress experienced by the introduced bryozoan Bugula neritina and identified a potential link between heat stress and the expression of a halogenated alkaloid, Caelestine A. The present results might indicate that global warming might not cause identical effects everywhere on the prevalence of NIS, and that present local climate may play a role in their response.

Introduction

Climate change constitutes a major challenge for humanity and the situation has incited to strong political engagements to counter it (UN 2015;IPCC 2021). It is however unlikely that efforts suffice to limit global warming to +1.5°C, aimed for by the Paris agreement, and a scenario with a global temperature increase of +3.2°C seems most likely [START_REF] Raftery | Less than 2 °C warming by 2100 unlikely[END_REF]). Associated to this, sea surface temperature is also projected to increase with climate change, with mean estimates ranging from +0.7°C to +3°C above preindustrial levels depending on the considered scenario (IPCC 2019). Climate change already poses a significant threat on natural ecosystems [START_REF] Hanson | Drought disturbance from climate change: Response of United States forests[END_REF][START_REF] Ayres | Assessing the consequences of global change for forest disturbance from herbivores and pathogens[END_REF][START_REF] Hughes | Climate Change , Human Impacts , and the Resilience of Coral Reefs[END_REF][START_REF] Lejeusne | Climate change effects on a miniature ocean: the highly diverse, highly impacted Mediterranean Sea[END_REF][START_REF] Waldock | Overview Articles Temperature-Driven Biodiversity Change: Disentangling Space and Time[END_REF][START_REF] Lonhart | Shifts in the distribution and abundance of coastal marine species along the eastern Pacific Ocean during marine heatwaves from 2013 to 2018[END_REF]), but as urban areas constitute islands of increased heat due to human activities and building materials, climate change effects may be even exacerbated [START_REF] Chapman | The impact of urbanization and climate change on urban temperatures: a systematic review[END_REF]. These combined effects of climate and urbanization can be observed in terrestrial, aquatic and marine urban environments and are affecting species characteristics [START_REF] Luo | Phenological responses of plants to climate change in an urban environment[END_REF][START_REF] Chapman | The impact of urbanization and climate change on urban temperatures: a systematic review[END_REF][START_REF] Brans | The heat is on: genetic adaptation to urbanization mediated by thermal tolerance and body size[END_REF][START_REF] Todd | Towards an urban marine ecology: characterizing the drivers, patterns and processes of marine ecosystems in coastal cities[END_REF]. Understanding how climate change influences urban ecosystems, which are already subjected to cumulative anthropic impacts, constitutes a grand challenge of high priority [START_REF] Borja | Past and Future Grand Challenges in Marine Ecosystem Ecology[END_REF].

Climate change may substantially modify biological invasion patterns through niche modifications inducing new areas to become available for invasive species [START_REF] Zhang | Impacts of climate change on geographical distributions of invasive ascidians[END_REF][START_REF] Vilizzi | A global-scale screening of non-native aquatic organisms to identify potentially invasive species under current and future climate conditions[END_REF]. NIS sometimes completely restructure the ecosystem leading to the loss of biodiversity and ecosystem services [START_REF] Pejchar | Invasive species, ecosystem services and human wellbeing[END_REF][START_REF] Johnston | Measuring Bahamian lionfish impacts to marine ecological services using habitat equivalency analysis[END_REF][START_REF] Walsh | Invasive species triggers a massive loss of ecosystem services through a trophic cascade[END_REF]. Besides the long term warming, climate change is projected to further increase the frequency of extreme weather events such as heat waves which can cause mass mortalities and restructure communities in natural [START_REF] Garrabou | Mass mortality in Northwestern Mediterranean rocky benthic communities: Effects of the 2003 heat wave[END_REF][START_REF] Lejeusne | Climate change effects on a miniature ocean: the highly diverse, highly impacted Mediterranean Sea[END_REF][START_REF] Collier | Temperature extremes reduce seagrass growth and induce mortality[END_REF] and in urban ecosystems [START_REF] Smale | Disentangling the impacts of heat wave magnitude, duration and timing on the structure and diversity of sessile marine assemblages[END_REF][START_REF] Castro | Winners and losers : prevalence of non-indigenous species under simulated marine heatwaves and high propagule pressure[END_REF]. This effect may directly benefit NIS as many have been shown to have a higher temperature tolerance [START_REF] Zerebecki | Temperature tolerance and stress proteins as mechanisms of invasive species success[END_REF][START_REF] Lejeusne | Do invaders always perform better? Comparing the response of native and invasive shrimps to temperature and salinity gradients in south-west Spain[END_REF]Kenworthy et al. 2018a;[START_REF] Gewing | Temperature resilience facilitates invasion success of the solitary ascidian Herdmania momus[END_REF]. Increased warming can also exacerbate the effect of pollution on local fauna by accelerating metabolism [START_REF] Sokolova | Interactive effects of metal pollution and temperature on metabolism in aquatic ectotherms: Implications of global climate change[END_REF][START_REF] Cabral | Synergistic effects of climate change and marine pollution: An overlooked interaction in coastal and estuarine areas[END_REF][START_REF] Todd | Towards an urban marine ecology: characterizing the drivers, patterns and processes of marine ecosystems in coastal cities[END_REF]. This may even further advantage NIS as they appear more tolerant to pollution as well [START_REF] Piola | Pollution reduces native diversity and increases invader dominance in marine hard-substrate communities[END_REF][START_REF] Osborne | Copper pollution enhances the competitive advantage of invasive ascidians[END_REF]. In the context of invasion, climate change might thus be a particular management challenge as it may increase the impact of NIS. For now, many NIS are established in their introduced range but are not considered invasive as they do not proliferate and do not cause adverse effects on the environment [START_REF] Occhipinti-Ambrogi | A uniform terminology on bioinvasions: A chimera or an operative tool?[END_REF][START_REF] Blackburn | A proposed unified framework for biological invasions[END_REF][START_REF] Ojaveer | Classification of Non-Indigenous Species Based on Their Impacts: Considerations for Application in Marine Management[END_REF]). Yet some of these NIS may be sleeper species, which are currently impeded in their proliferation by environmental conditions and may become invasive if these conditions change due to climate change [START_REF] Bradley | Regional Invasive Species & Climate Change Management Challenge: Preparing for sleeper species[END_REF][START_REF] Spear | The Invasion Ecology of Sleeper Populations: Prevalence, Persistence, and Abrupt Shifts[END_REF]. Furthermore, climate exerts strong selective pressures on native and introduced populations and some species, including NIS, might rapidly adapt to climate change through genetically mediated differences in thermal tolerances [START_REF] Brans | The heat is on: genetic adaptation to urbanization mediated by thermal tolerance and body size[END_REF][START_REF] Sork | Genomic Studies of Local Adaptation in Natural Plant Populations[END_REF] or through plastic responses [START_REF] Zerebecki | Temperature tolerance and stress proteins as mechanisms of invasive species success[END_REF][START_REF] Sork | Genomic Studies of Local Adaptation in Natural Plant Populations[END_REF]. Divergences in thermal tolerances of species from different geographical locations may be good indicators of their capacity to adapt to climate change [START_REF] Sorte | Geographic variation in temperature tolerance as an indicator of potential population responses to climate change[END_REF]. It is thus crucial to foster the understanding of how climate change might impact marine urban communities, especially considering the important place the NIS occupy.

The effect of global warming and/or heatwaves on a variety of coastal ecosystems has most of the time been conducted through in aquaria and mesocosm studies (Sorte et al. 2010a;b;[START_REF] Smale | Disentangling the impacts of heat wave magnitude, duration and timing on the structure and diversity of sessile marine assemblages[END_REF][START_REF] Wahl | A mesocosm concept for the simulation of near-natural shallow underwater climates: The Kiel Outdoor Benthocosms (KOB)[END_REF][START_REF] Pansch | Heat waves and their significance for a temperate benthic community: A near-natural experimental approach[END_REF][START_REF] Castro | Winners and losers : prevalence of non-indigenous species under simulated marine heatwaves and high propagule pressure[END_REF]) but completely in situ experiments remain scarce (examples are [START_REF] Smale | Turning on the Heat : Ecological Response to Simulated Warming in the Sea[END_REF][START_REF] Smale | Community responses to seawater warming are conserved across diverse biological groupings and taxonomic resolutions[END_REF][START_REF] Ashton | Warming by 1°C Drives Species and Assemblage Level Responses in Antarctica's Marine Shallows[END_REF], due to technical and experimental difficulties of in situ water warming. Marinas are the major model for studying marine urban communities, since they combine the different disturbances exerted by urbanization and have a high abundance of NIS, but they also offer logistic advantages making it possible to conduct in situ water warming experiments in an anthropogenic context.

In the present study we thus aim to assess the in situ effect of +3°C sea surface warming (RCP8.5; IPCC 2019) in two marinas, one in the NW Mediterranean and one in the NE Atlantic. We deployed indirectly heated hot plates in open systems in two marinas and assessed differences in sessile community characteristics from not heated control panels. We also focused on the metabolism of a common marina species (reported as NIS in both surveyed marinas) to see if it was affected under increased heat. The observed results could help to make predictions about the effect of climate change on these two marinas, especially in regard to the prevalence of NIS, which are supposed to be advantaged by increased heat.

Materials and Methods

Study sites

The two sites selected for our study were the Marina du Chateau in Brest (henceforth "Chateau"; Finistere, France; 48°22'45.4"N, 4°29'21.2"W) and the Port de Carry-le-Rouet (henceforth "Carry";

Bouches-du-Rhône, France; 43°19'47.0"N, 5°09'09.5"E). Both marinas have similar boating capacities (650 and 590 berths respectively), are semi-enclosed, and subject to anthropic impacts due to their position in their urban area and water runoffs at the inner parts of the marinas. We deployed experimental setups at the entrance of both marinas to avoid a synergetic and confounding effect of strong pollution in the inner parts of marinas [START_REF] Je | Changes in benthic communities along a presumed pollution gradient in Vancouver Harbour[END_REF]Kenworthy et al. 2018b).

Experimental setup

For each marina we deployed a structure with 10 settlement panels (Correx Polyethylene, 20 cm x 20 cm; Fig. 1). These recruitment panels were randomly assigned to be either control panels (5 per marina) or being constantly heated at +3°C above the recorded ambient water temperature (5 per marina). The experiment lasted 50 days in both marinas and started in mid-August 2020 in Carry and early September for Chateau. To reach a stable water temperature difference with ambient water, the settlement panels were placed at the top of transparent bottom opened acryl cylinders. Water column warming of the cylinders was achieved through 200W aquaria resistors (approx. 20 cm away from the panels) regulated with a built-in computer (Raspberry Pi3+) coupled to thermometers (DS18B20) measuring the outside water temperature (Fig. 1, Sup. Fig. 1). The water temperature around the heating elements and directly around the panels was regularly monitored to avoid overheating.

Passthroughs at the top of the cylinders aimed to increase connectivity with the surrounding water and to avoid excessive stratification of the water within the cylinder.

Fauna identification and community analysis

At the end of the experiment all panels were transported to the laboratory where they were maintained in open circuit aquaria until identification by a single observer in less than 36 h. A surface of 15 cm x 15 cm was analyzed by superimposing a 100-point grid allowing for a resolution of more than 0.4 points per cm² necessary for a reliable estimation of the surface cover [START_REF] Taormina | Optimizing image-based protocol to monitor macroepibenthic communities colonizing artificial structures[END_REF].

At each intersection of the grid all living individuals of any species were counted. Species present on a panel that did not fall within an intersection were counted with 1% cover. Identification went to the lowest possible taxonomical level [START_REF] Zibrowius | Étude Morphologique, Systématique et Écologique, des Serpulidae (Annelida Polychaeta) de la Région de Marseille[END_REF][START_REF] Hayward | British ascophoran bryozoans[END_REF][START_REF] Hayward | Cyclostome Bryozoans (Synopses of the British Fauna 34)[END_REF], 1995[START_REF] Hayward | Cheilostomatous bryozoa (Part 1)[END_REF][START_REF] Bianchi | Guide per il riconoscimento delle specie animali delle acque lagunari e costiere italiane -Policheti Serpuloidei[END_REF][START_REF] Riedl | Fauna und Flora des Mittelmeeres: Ein systematischer Meeresführer für Biologen und Naturfreunde[END_REF][START_REF] Zabala | Illustrated keys for the classification of Mediterranean Bryozoa[END_REF][START_REF] Brunetti | Ascidiacea of the European waters[END_REF]. After analysis, the panels were scraped, and biomass was dried at 60°C for one week, then weighted. To determine ash free dry mass (AFDM) the dry communities were burned at 520°C for 6 h.

The community matrix has been analyzed using the 'vegan' R package (version 2.4-6;[START_REF] Oksanen | Vegan: community ecology package[END_REF]. The Shannon (H') and Simpson (λ) diversity indices as well as Pielou (J') evenness and species richness were assessed for each panel. They were compared between the control and heated treatment with an ANOVA. A Poisson GLM was applied for species richness since it constitutes count data. The community matrix was transformed into a Bray-Curtis dissimilarity matrix and multivariate homogeneity of groups was verified with the 'betadisper' function. A Principal Coordinate Analysis (PCoA) was used to visualize the community of both treatments for the two marinas. A subsequent PERMANOVA (10⁵ perm.) was used to assess for significant differences between treatments of each marina. To identify the species contributing to the differences between treatments a SIMPER analysis followed by an Multipattern analysis using the 'indicSpecies' package was conducted (10 4 perm.; ver. 1.7.1;[START_REF] Cáceres | A framework for estimating niche metrics using the resemblance between qualitative resources[END_REF][START_REF] De Cáceres | R package: indicspecies[END_REF]. An ANOVA was used on the AFDM to determine if significant differences were present between treatments for each marina. The cover occupied by empty space without any species and the cover exerted by NIS, was analyzed via a Poisson glm since they originated from count data (One count = 1% cover).

Metabolome identification and analysis

The invasive bryozoan Bugula neritina was used for metabolomic analyses, as it usually occurs in both marinas due to its worldwide introduction and having been previously studied in the Chateau marina including for metabolomics (Gauff et al. in prep a). Unfortunately, biomass in the experiment from the Carry marina was too low for metabolomic sampling, so analyses were only conducted in the Chateau marina. For each panel, one large colony of B. neritina was sampled. Samples were immediately frozen in dry ice, then stored at -80°C, freeze dried and ground to dust. The total metabolome was extracted out of 5 mg of dry powder with 3 times 1 mL MeOH/DCM (CH2Cl2; 1:1) and ultrasonication at room temperature during 2 min. Supernatants were filtered on PTFE syringe filters (13 mm, 0.22 µm, Restek®). Two samples (one of each treatment) were excluded due to strong contamination. Before injection into Liquid Chromatography Mass Spectrometry (LC-MS) the samples were resuspended in 500 µL MeOH and a pool sample was created by mixing 30 µL of each. Two blank samples, which were submitted to the entirety of the extraction process were additionally produced. LC-MS were performed with UHPLC (Dionex Ultimate 3000 with RS Pump, autosampler, thermostat column and UV diode array, Thermo Scientific®) coupled to a quadrupole Time of Flight spectrometer (qToF) with an ESI source (Impact II, Bruker Daltonics®). We conducted UHPLC separation on an Acclaim TM RSLC 120C18 column (150 mm x 2.1 mm, 2.2 µm, ThermoScientific®) at 40°C at 0.5 mL min -1 elution rate. Water (Polar, LC/MS grade, Carlo Erba®) with 0.1% formic acid (A) and acetonitrile (Apolar, LC/MS grade) with 0.1% formic acid (B) were used as chromatographic solvents. Chromatography cycles were set as follows: 1) mix of 5:95 of A:B during 2 min, 2) linear increase up to 100% B during 8 min followed by 100% B during 4 min, 3) return to initial conditions (5:95, A:B) for 3 min (17 min total). We injected 10 µL in MS1+ and MS2+ with N2 at 3.5 bars as nebulizer gas, N2 at 12 L min -1 as dry gas, capillary temp. at 200 °C and voltage at 3500 V for a scan from 50 to 1200 mass units at 2 Hz. A spectrometer calibration with a formate:acetate solution was performed before each sample and the pool sample was injected every 4 samples to correct for source drift. Blank samples and pure solvent were injected to filter for ions linked to the extraction process.

LC-MS profiles were automatically recalibrated with the internal calibration in 'Bruker Compass DataAnalysis' (version 4.3) resulting in a m/z precision of 2 ppm. The netCDF data was processed with the 'XCMS' R package (version 3.12; [START_REF] Smith | XCMS: processing mass spectrometry data for metabolite profiling using nonlinear peak alignment, matching, and identification[END_REF]) with peak picking for the detection of features, obiwarp retention time correction, feature grouping (time frame 5 s) and peak filling. Ions were manually normalized in Microsoft Excel [START_REF] Van Der Kloet | Analytical error reduction using single point calibration for accurate and precise metabolomic phenotyping[END_REF]) and were filtered to exclude ions originating from the blanks, based on signal/noise ratio (S/N = 10), ions with high variability in the pool sample (> 0.2) and ions that were auto-correlated (> 0.8). To increase relevancy, we excluded the first 100 seconds of chromatography of our analysis.

Metabolome analyses were mainly performed in 'Metaboanalyst' [START_REF] Xia | MetaboAnalyst: A web server for metabolomic data analysis and interpretation[END_REF]). Our matrix had less than 5000 features and was log transformed and auto scaled. A Principal Correspondence Analysis (PCA) was used to visualize our sample structure. A PPLS-DA driven permutation Model Validation Analysis (MVA) (10 5 perm., 6 components) using the 'RVAideMemoire' package in R (version 0.9-79; Hervé 2021) tested for statistical differences between the treatments. A Volcano plot from 'Metaboanalyst' identified features with a significant foldchange (t-Test, p < 0.05; FC > 1.5). These molecules were analyzed with fragmentation mass spectrometry (MS2) spectra. Molecular networking via Global Natural Products Social Molecular Networking (GNPS; [START_REF] Wang | Sharing and community curation of mass spectrometry data with Global Natural Products Social Molecular Networking[END_REF]) and visualized in 'cytoscape' (version 3.7.2) helped ion identification.

Molecules of particular interest, requiring further confirmation upon identification were compared against a commercially produced standards (Ambinter®, Greenpharma) by injecting the sample with the highest intensity and the commercial extract in a short injection series following the same parameters as described earlier.

Results

Community structure

None of the diversity indices had significant differences between control and heated panels in any marina (Carry: H' J' λ, ANOVA, |R²| < 0.27, p > 0.07; S, Poisson glm, z = -0.52, p = 0.62; Chateau: H' J' λ, ANOVA, |R²| < 0.12; p > 0.33; S, Poisson glm, z = 0.25, p = 0.81). The species list for each marina can be found in Sup. Tab. 1. On the PCoA of the community structure (Fig. 2), panels were distributed according to their treatment (Control vs Heated) on the PCo1 axis, explaining 40.6% of the variability between samples for Carry (Fig. 2 A) and 70.7% for Chateau (Fig. 2 B). The PCo2 axis (< 26.3%) corresponded to intra treatment variability. For both marinas, significant effects of the warming were observable on the communities (PERMANOVA; Carry: R² = 0.34, p = 0.008; Chateau: R ²= 0.54, p = 0.008). The SIMPER analysis revealed that in both cases, 5 species contributed to more than 70% of the contrast between the treatments (over 85% for Chateau; Fig. 3; Tab. 1). At the Carry marina, only the bryozoan Aetea sica has been identified as indicator species for the control (Tab. 1). In the Château marina however, most of the contributors identified by SIMPER analysis were also indicator species of the control treatment (Tab. 1). For this marina two additional species were identified as indicator species for the heated treatment: the bryozoan Cryptosula pallasiana and the annelid Spirobranchus triqueter (Tab. 1). Tab. 1: Species most contributing to the contrast between the control and heated (+3°C) communities and indicator species for the treatments. Species in bold with an asterisk are non-indigenous species. SIMPER analysis score indicates the percentage a species contributes to the contrast between treatments (see also Fig. 3). The direction of the warming effect on the percent cover occupied by these species is indicated. Indicator species for each treatment identified through Multipattern analysis are indicated with their association statistic (Stat.) and the p-value. SIMPER scores in parentheses are scores that are not part of the most important contributors. Ash free dry mass (AFDM) varied significantly between treatments for both marinas (ANOVA, R² > 0.56, p < 0.005; Fig. 4) with lower biomass in the heated treatment. Accordingly, empty space was significantly higher in the heated treatment for both marinas (Poisson glm, z > 2.28, p < 0.022; Fig. 4 A B). The absolute cover exerted by NIS as well as the ratio of NIS cover on total cover (proportional cover) was higher in the heated treatment for Carry-le-Rouet (Poisson glm, z = 2.149, p = 0.0317; Fig.

Carry-le-Rouet

Species

Metabolome analysis

The metabolome of Bugula neritina showed a marginally significant difference between the control and the heated treatment (Fig. 5; MVA: CER = 0.181, p-value = 0.052). A total of 12 molecules were found to vary significantly between treatments. Seven had a significant 1.5-fold under expression while 5 had a significant 1.5-fold overexpression (t-Test, p < 0.05; Fig. 6). From these molecules, we were able to identify 6 (Tab. 2). A molecule identified as a Bromo-Oxo-Quinoline-Carboxylic Acid (BQA) had a fragmentation pattern (Sup. Tab. 2: Identified molecules with a significant fold change (FC) between heated vs control treatments in Bugula neritina in the Chateau marina. Identification confidence levels are indicated [START_REF] Schymanski | Identifying small molecules via high resolution mass spectrometry: Communicating confidence[END_REF]

Discussion

The present results indicate that long-term warming due to climate change may affect marina communities and the metabolism of certain species. We chose to simulate in situ a seawater warming of +3°C compared to the control in accordance with the mean estimate of global surface temperature increase projected in RCP8.5 (IPCC 2019). Simulated warming in mesocosms (Sorte et al. 2010b) and in situ [START_REF] Smale | Turning on the Heat : Ecological Response to Simulated Warming in the Sea[END_REF][START_REF] Smale | Community responses to seawater warming are conserved across diverse biological groupings and taxonomic resolutions[END_REF][START_REF] Ashton | Warming by 1°C Drives Species and Assemblage Level Responses in Antarctica's Marine Shallows[END_REF]) has been shown to strongly impact communities favoring opportunistic species and NIS. In accordance with these previous observations, our present case revealed significant differences in the recruited communities between heated and control treatments in both marinas, which was associated to lower biomass (AFDM) and higher empty space in the heated treatment. At Carry-le-Rouet (Carry, Mediterranean Sea), this coincided with higher absolute and proportional abundance of non-indigenous species in the heated treatment. In the marina du Chateau de Brest (Chateau, Atlantic) however, NIS were negatively impacted by warming, just as were natives as shown by the constant proportional abundance. Yet, Cryptosula pallasiana, native from the Mediterranean but introduced in the Atlantic, has been identified as an indicator species of the heated treatment in the Chateau marina. No data on the thermal tolerance of this species is available in the literature. However, it seems reasonable to assume that this tolerance should be above the temperatures observed in the Atlantic marina due to its origin (Mediterranean Sea). In its introduced habitat, it has been shown to almost double its growth rate at 18°C compared to 4°C [START_REF] Amui-Vedel | Zooid size and growth rate of the bryozoan Cryptosula pallasiana Moll in relation to temperature, in culture and in its natural environment[END_REF], which may explain our observations. While anecdotal, this might indicate that this species could have profited from our simulated warming.

The higher temperature did however not only influence the settlement dynamics on the panels, but also had an impact on the metabolism of our model species. In Carry we were unfortunately unable to sample B. neritina since too few individuals were present (mainly in the control). In the Chateau marina however, a slight difference between the metabolome of the heated and the control treatment is visible in the PCA and MVA test.

A total of 12 molecules had significant fold changes (FC > 1.5) between the treatments. Most notably two phospholipids (Glycerophosphoethanolamine, Lyso-PE (22:6) and Glycerophosphocholine, Lyso-PC (22:6)) and a Bromo-Oxo-Quinolinecarboxylic Acid (BQA), were identified as significantly affected molecules. We were able to confirm the BQA structure and to match it with a commercially available standard for Caelestine A, a molecule described in the ascidian Aplidium caelestis [START_REF] Yin | Caelestines A-D, Brominated Quinolinecarboxylic Acids from the Australian Ascidian Aplidium caelestis[END_REF]. QA derivates (especially brominated ones) have rarely been identified in nature and their exact function is thus difficult to establish [START_REF] Yin | Caelestines A-D, Brominated Quinolinecarboxylic Acids from the Australian Ascidian Aplidium caelestis[END_REF]. Caelestines have been shown to have minor cyto-and bacterio-toxicity at high concentrations [START_REF] Yin | Caelestines A-D, Brominated Quinolinecarboxylic Acids from the Australian Ascidian Aplidium caelestis[END_REF]. Quinolinecarboxylic acids (QA) seem to inhibit the activity of one of the most common protein kinases, Protein Kinase CK2 [START_REF] Syniugin | Design, synthesis and evaluation of 3-quinoline carboxylic acids as new inhibitors of protein kinase CK2 Design, synthesis and evaluation of 3quinoline carboxylic acids as new inhibitors of protein kinase CK2[END_REF]), which has numerous functions including inhibiting apoptosis and whose activity seems to increase with cellular stress [START_REF] Manni | Human Cancer Biology Protein Kinase CK2 Protects Multiple Myeloma Cells from ER Stress-Induced Apoptosis and from the Cytotoxic Effect of HSP90 Inhibition through Regulation of the Unfolded Protein Response[END_REF][START_REF] Piazza | Protein kinase CK2 in hematologic malignancies: reliance on a pivotal cell survival regulator by oncogenic signaling pathways[END_REF]. Considering the very high overexpression of Two phospholipids varied between the control and the heated treatment (Lyso-PC (22:6) and Lyso-PE (22:6)). Phospholipids play a major role in cell membrane fluidity, which changes as a response to ambient temperature [START_REF] Pruitt | Adaptations to temperature in the cellular membranes of crustacea: membrane structure and metabolism[END_REF][START_REF] Cuculescu | Temperature acclimation of marine crabs: Changes in plasma membrane fluidity and lipid composition[END_REF][START_REF] Sung | Acquired tolerance to temperature extremes[END_REF][START_REF] Chen | Profiling lipidome changes of Pyropia haitanensis in short-term response to high-temperature stress[END_REF]. The type of lipids used as well as the number of unsaturation is a key component in this response [START_REF] Pruitt | Adaptations to temperature in the cellular membranes of crustacea: membrane structure and metabolism[END_REF][START_REF] Cuculescu | Temperature acclimation of marine crabs: Changes in plasma membrane fluidity and lipid composition[END_REF]. It is interesting to note that in our case, both phospholipids have the same adducts. Lyso-PC (22:6) might results directly from the conversion of Lyso-PE (22:6) trough successive methylation by the phosphatidylethanolamine N-methyltransferase [START_REF] Gibellini | The Kennedy pathway-de novo synthesis of phosphatidylethanolamine and phosphatidylcholine[END_REF],

indicating that this conversion is associated to the increased heat. In the marine seaweed Pyropia sp.

long term warming and short term heat shocks modified Lyso-PC and Lyso-PE concentrations generally increasing the first and decreasing the latter [START_REF] Chen | Profiling lipidome changes of Pyropia haitanensis in short-term response to high-temperature stress[END_REF][START_REF] Song | Heat shockinduced metabolic conversion of membrane lipids, fatty acids and volatile organic compounds of Pyropia haitanensis under different heat shock time[END_REF], which corresponds to our observed results. Furthermore, PC and PE of different length including PC (22:6) and PE (22:6) have been shown to vary in the crabs Carcinus maenas and Cancer pagurus in respect to temperature [START_REF] Cuculescu | Temperature acclimation of marine crabs: Changes in plasma membrane fluidity and lipid composition[END_REF]. The same Lyso-PE has also been previously described in B. neritina, varying according to environmental conditions in the Chateau marina (Gauff et al. in prep a). It may constitute a marker of good environmental conditions as Lyso-PE may play a role in gametogenesis in certain marine organisms [START_REF] Ivanisevic | Lysophospholipids in the Mediterranean Sponge Oscarella tuberculata: Seasonal Variability and Putative Biological Role[END_REF].

It is interesting to note that despite apparent similarities of the effect of +3°C warming in both marinas (Community differences, lower biomass, more empty space), the effect of warming on NIS prevalence was opposite, including B. neritina: in the Chateau marina the higher temperature led to lower abundance of NIS indicating that the generated heat stress may cause a disadvantage. This is further supported by the metabolomic differences between heated and control treatment in the Chateau marina, potentially indicating heat stress in B. neritina. Yet, in Carry-le-Rouet NIS, including B. neritina, had an increased abundance in the heated treatment. This last result is more in adequation with anterior results from the literature, with opportunistic species and NIS having higher biomass in warm treatments [START_REF] Smale | Turning on the Heat : Ecological Response to Simulated Warming in the Sea[END_REF][START_REF] Smale | Community responses to seawater warming are conserved across diverse biological groupings and taxonomic resolutions[END_REF][START_REF] Ashton | Warming by 1°C Drives Species and Assemblage Level Responses in Antarctica's Marine Shallows[END_REF]. Moreover, an increased growth of B. neritina under higher temperatures has also been reported [START_REF] Lord | Impact of seawater temperature on growth and recruitment of invasive fouling species at the global scale[END_REF]) and numerous NIS, including B.

neritina, have been shown to have high thermal tolerance [START_REF] Hitoshi | Growth of the Bryozoan Bugula neritina in the Sea at Various Water Temperatures[END_REF][START_REF] Zerebecki | Temperature tolerance and stress proteins as mechanisms of invasive species success[END_REF][START_REF] Lejeusne | Do invaders always perform better? Comparing the response of native and invasive shrimps to temperature and salinity gradients in south-west Spain[END_REF][START_REF] Lord | Impact of seawater temperature on growth and recruitment of invasive fouling species at the global scale[END_REF]Kenworthy et al. 2018a). However, strong variation (up to 2°C) in temperature tolerance among geographically separated populations of several fouling species, including B. neritina, has been reported [START_REF] Sorte | Geographic variation in temperature tolerance as an indicator of potential population responses to climate change[END_REF][START_REF] Zerebecki | Temperature tolerance and stress proteins as mechanisms of invasive species success[END_REF]. This may be due to adaptation and/or acclimation of populations to their local climate and its variability [START_REF] Sorte | Geographic variation in temperature tolerance as an indicator of potential population responses to climate change[END_REF]. Tolerance to temperature-related stress may be acquired through numerous cellular pathways [START_REF] Sung | Acquired tolerance to temperature extremes[END_REF]. Such tolerance can be increased in an invasive species by the exposure to stress [START_REF] Clutton | Early developmental stages of native populations of Ciona intestinalis under increased temperature are affected by local habitat history[END_REF]) by modifying epigenetics [START_REF] Ardura | Stress related epigenetic changes may explain opportunistic success in biological invasions in Antipode mussels[END_REF] or be modulated by maternal effects [START_REF] Burgess | Temperature-induced maternal effects and environmental predictability[END_REF]. High temperature heatwaves in a marina close to Carry-le-Rouet have been shown to lead to summer mass mortality for numerous species including B. neritina (Gauff et al. in prep b). It may be possible that the regular exposure of Mediterranean sessile communities to high temperatures and heat waves [START_REF] Garrabou | Mass mortality in Northwestern Mediterranean rocky benthic communities: Effects of the 2003 heat wave[END_REF][START_REF] Lejeusne | Climate change effects on a miniature ocean: the highly diverse, highly impacted Mediterranean Sea[END_REF] increases the thermal tolerance of present NIS, thus allowing their spread in warmer water. Conversely, the water temperature in the vicinity of the Chateau marina is more constant due to rarer heatwaves and higher water mixing due to higher tides and wave action (no stratification). This may explain why the simulated warming strongly impacted the present NIS which would not be accustomed to such disturbance.

The present results indicate that in the short term, climate change related ocean surface warming might not have the same effect on the prevalence of non-indigenous species everywhere. Certain locations, where regular heat events are already experienced, might have a higher risk at favoring them compared to cooler and more stable climates. This might be attributed to heat resistant individuals being favored by the intense selection exerted by heat waves. The impact of climate change might thus also depend on how global warming modifies cooler temperate climates and if this increases extreme weather events like heat waves. It is crucial that further studies investigate in situ how heat waves and global warming influence marine urban communities and impact of NIS within. These studies might also integrate the conjoint effect of ocean acidification due to increased CO2 concentrations and temperature as they both are likely to highly influence, which species profit from climate change.

Fig. 1 :

 1 Fig. 1: in situ warming structure deployed in the marinas. Photo credit: Wilfried Thomas (CNRS Station Biologique de Roscoff)

Fig. 2 :

 2 Fig. 2: Principal Coordinate Analysis (PCoA) of the control and heated (+3°C) communities in the A) Carry-le-Rouet Marina (Mediterranean) and in B) Marina du Château (Atlantic). Significant differences in community structure: PERMANOVA: Carry (R² = 0.34, p = 0.008); Chateau (R ²= 0.54, p = 0.008)

  4 C; ANOVA, R² = 0.29, p value = 0.063). In the Chateau marina the absolute cover was significantly lower in the heated treatment while the proportional cover remained unaffected (Poisson glm, z = -5.36, p < 0.001; Fig. 4 D; ANOVA, NS).

Fig. 4 :

 4 Fig. 4: Boxplots of the percent cover of ash free dry mass (AFDM; A, B) unoccupied empty space (C, D) and the percent cover occupied by Non Indigenous Species (NIS; E, F) for the control and heated (+3°C) treatment in Carry-le-Rouet (Mediterranean; A, C, E) an d the Marina du Château (Atlantic; B, D, F).A: ANOVA, R² = 0.56, p = 0.008; B: ANOVA, R² = 0.60, p = 0.005; C:Poisson glm, z = 2.28, p = 0.022; D: Poisson glm, z = 7.57, p < 0.0001; E: Poisson glm, z = 2.149, p = 0.0317; F: Poisson glm, p < 0.001).

Fig. 6 :

 6 Fig. 5: Principal Component Analysis (PCA) of the

  Caelestine A in the heated treatment and the implication of QAs with the inhibition of CK2, this molecule might be linked to the stress response of Bugula neritina to heat stress. While yet speculative the implication of Caelestine A in heat stress response should be further investigated in future experiments with B. neritina and other species like A. caelestis.

warmly thank the harbor offices of Brest and Carry-le-Rouet (especially the harbor master Henri Salati), and the Aix-Marseille-Provence Métropole for letting us conduct experiments in the marinas. Special thanks also to Centre de Ressources Biologiques Marines (CRBM) and the Roscoff Aquarium Services (RAS), specifically Gaëtan Schires and Sébastien Henry Metabolomics analyses were conducted at Chemical Ecology and Metabolomics facilities at METABOMER Roscoff (Cedric Leroux) and IMBE (funded by ANR, French National Research Agency, INEE the Institute of Ecology and Environment of the CNRS, the Total Foundation and the French Sud PACA regional council). We also want to acknowledge Jean-Charles Martin (C2VN BIOMET platform, Marseille) for creating in house metabolome filtering scripts for R. Thank you Elsa Davies for English proofreading. Funding for this project has been provided through a PhD grant from the Sorbonne Université -Museum National d'Histoire Naturelle (Ecole Doctorale 227). This work was also supported by the French National Program EC2CO (Ecosphère Continentale et Côtière) and the program 'InPor' of the INEE-CNRS' PEPS EcoMob (PIs: C. Lejeusne and D. Davoult).

Funding for this project has been provided through a PhD grant from the Sorbonne Université -Museum National d'Histoire Naturelle (Ecole Doctorale 227). This work was also supported by the French National ProgramEC2CO (Ecosphère Continentale et Côtière) and the program 'InPor' of the INEE-CNRS' PEPS EcoMob (PIs: C. Lejeusne and D. Davoult).

Declarations

Funding

Competing Interests

The authors have no relevant financial or non-financial interests to disclose.

Author Contributions