Thibaut Morel-Journel
email: thibaut.morel@uclouvain.be

Claire Rais Assa

Ludovic Mailleret

Elodie Vercken

It's all about connections: hubs and invasion in habitat networks

Keywords: Allee effect, connectivity, establishment, hub, individual-based model, invasion, microcosm, network, simulation, spread

During the early stages of invasion, the interaction between the features of the invaded landscape, notably its spatial structure, and the internal dynamics of an introduced population, has a crucial impact on establishment and spread. By approximating introduction areas as networks of patches linked by dispersal, we characterized their spatial structure with specific metrics and tested their impact on two essential steps of the invasion process: establishment and spread. By combining simulations with experimental introductions of Trichogramma chilonis (Hymenoptera: Trichogrammatidae) in artificial laboratory microcosms, we demonstrated that spread was hindered by clusters and accelerated by hubs but was also affected by smallpopulation mechanisms prevalent for invasions, such as Allee effects. Establishment was also affected by demographic mechanisms, in interaction with network metrics. These results highlight the importance of considering the demography of invaders as well as the structure of the invaded area to predict the outcome of invasions.

Introduction

Managing invasions becomes increasingly costly and decreasingly efficient with the expansion of the invaded area [START_REF] Simberloff | Impacts of biological invasions: what's what and the way forward[END_REF]. In invasion biology, focus is placed on the processes occurring during the first generations after the introduction of an exotic species, before its proliferation and spread besides the introduction site. The dispersal abilities of individuals during invasions and the way they evolve over the course of the range expansion have been documented theoretically (Travis et al. 2009;Burton et al. 2010) and experimentally, notably among protists [START_REF] Fronhofer | Eco-evolutionary feedbacks during experimental range expansions[END_REF] and arthropods [START_REF] Ochocki | Rapid evolution of dispersal ability makes biological invasions faster and more variable[END_REF][START_REF] Weiss-Lehman | Rapid trait evolution drives increased speed and variance in experimental range expansions[END_REF]. Although the role of individual characteristics on the spread patterns is well known, the way individuals disperse across a landscape also depends on the interaction between individuals and the features of the environment (Calabrese & Fagan 2004), which can themselves affect the evolution of dispersal [START_REF] Baguette | Landscape connectivity and animal behavior: functional grain as a key determinant for dispersal[END_REF]. Therefore, understanding the structure of the invaded landscape is also essential to understand the patterns of colonization observed during the first stages of an invasion.

Networks have gained popularity in the last decades as a method to represent the spatial structure of such landscapes (Urban & [START_REF] Keitt | Allee effects, invasion pinning, and species' borders[END_REF]Minor & Urban 2007). They are used to represent habitat patches suitable for the species considered (the vertices of the network) and the way dispersal connects them (the edges of the network). Networks are a powerful tool associated with various metrics and used across a variety of scientific fields, but they are still seldom used in invasion biology. In this field, they are classically used to describe food webs, to study how they are impacted by non-native species or to assess the invasibility of a community based on its trophic structure [START_REF] Romanuk | Predicting invasion success in complex ecological networks[END_REF][START_REF] Lurgi | Network complexity and species traits mediate the effects of biological invasions on dynamic food webs[END_REF]. A few studies also represent the spatial structure of riverine systems as networks, to study the impact of the network structure on the composition of a community invading a new environment [START_REF] Seymour | Active colonization dynamics and diversity patterns are influenced by dendritic network connectivity and species interactions[END_REF][START_REF] Alther | Fluvial network topology shapes communities of native and non-native amphipods[END_REF]. This study aims at using networks in another way: to investigate the establishment and spread of the invader, with a special focus on small-population demographic processes.

Although few usages of networks are recorded in invasion biology, other fields provide insight about the impact of the spatial structure of the introduction area on the first stages of an invasion. Notably, several epidemiology studies correlate specific metrics with patterns of spread. A prime example is the combination of networks with models derived from classical Susceptible-Infected (SI) models [START_REF] Kermack | A Contribution to the Mathematical Theory of Epidemics[END_REF], e.g. to study the spread of nosocomial infections in hospitals [START_REF] Ueno | Controlling nosocomial infection based on structure of hospital social networks[END_REF] or the transmission of pathogens through grooming among macaques [START_REF] Romano | Modeling infection transmission in primate networks to predict centrality-based risk[END_REF]. These studies usually consider networks made up of hosts (the vertices) connected by social interactions (the edges). Other studies consider larger-scale networks, in which vertices are populations of hosts and edges are host movements between populations [START_REF] Arino | Disease spread in metapopulations[END_REF]. Such a framework has for example been used to explain the extent of plague during the fourteenth century [START_REF] Gómez | Network theory may explain the vulnerability of medieval human settlements to the Black Death pandemic[END_REF] or to assess epidemic risks in the Japanese airline network [START_REF] Tanaka | Random and targeted interventions for epidemic control in metapopulation models[END_REF].

The conceptual similarities between the spread of disease outbreaks and the spread of invading organisms have led studies in invasion biology to use models originating from epidemiological modelling. They notably have been used to study the impact of landscape structure on the risks of invasive spread through forests [START_REF] Ferrari | Initial conditions and their effect on invasion velocity across heterogeneous landscapes[END_REF] or across marinas via ballast waters [START_REF] Floerl | The importance of transport hubs in steppingstone invasions[END_REF]. However, these studies fail to consider the internal dynamics of each patch. Moreover, this impact of network structure on invasion remains to be experimentally tested.

Networks are also often used in conservation biology, notably to map the structure of existing metapopulations, in terrestrial (Bunn et al. 2000;Urban & Keitt 2001), marine (Treml et al. 2008;Almany et al. 2009) and riverine contexts [START_REF] Schick | Directed connectivity among fish populations in a riverine network[END_REF]. Network metrics are notably used to describe the connectivity between sub-populations, to study the mean lifetime of metapopulations (e.g. [START_REF] Bode | Using complex network metrics to predict the persistence of metapopulations with asymmetric connectivity patterns[END_REF][START_REF] Drechsler | Predicting metapopulation lifetime from macroscopic network properties[END_REF][START_REF] Kininmonth | Metapopulation mean life time within complex networks[END_REF][START_REF] Shtilerman | The effects of connectivity on metapopulation persistence: network symmetry and degree correlations[END_REF] or the extinction risks of sub-populations (e.g. [START_REF] Gilarranz | Spatial network structure and metapopulation persistence[END_REF][START_REF] Peck | Networks of habitat patches in tsetse fly control: Implications of metapopulation structure on assessing local extinction probabilities[END_REF][START_REF] Webb | How does network structure and complexity in river systems affect population abundance and persistence?[END_REF]. They are also used to identify essential populations to maintain connectivity in a landscape [START_REF] Bodin | Ranking individual habitat patches as connectivity providers: integrating network analysis and patch removal experiments[END_REF][START_REF] Baranyi | Contribution of habitat patches to network connectivity: redundancy and uniqueness of topological indices[END_REF][START_REF] Watson | Identifying critical regions in small-world marine metapopulations[END_REF]). Yet, those studies often focus on the dispersal capabilities of individuals [START_REF] Bodin | Ranking individual habitat patches as connectivity providers: integrating network analysis and patch removal experiments[END_REF][START_REF] Baranyi | Contribution of habitat patches to network connectivity: redundancy and uniqueness of topological indices[END_REF] or consider local population dynamics only through global colonization and extinction probabilities, without considering explicitly the internal dynamics of the sub-populations [START_REF] Bode | Using complex network metrics to predict the persistence of metapopulations with asymmetric connectivity patterns[END_REF][START_REF] Gilarranz | Spatial network structure and metapopulation persistence[END_REF][START_REF] Shtilerman | The effects of connectivity on metapopulation persistence: network symmetry and degree correlations[END_REF]. The studies considering the internal dynamics of populations more precisely and its interaction with network features usually concern large, already established populations, and their susceptibility to disturbances, biotic (Mari et al. 2014) or abiotic [START_REF] Gilarranz | Effects of network modularity on the spread of perturbation impact in experimental metapopulations[END_REF]. Conversely, the present study focuses on small populations during their growth phase, and on specific associated mechanisms that can interact with the network structure of the landscape.

Invaders usually experience multiple demographic bottlenecks, firstly at the initial introduction and subsequently at each colonization event, which can strongly impact the success of an invasion. Mechanisms associated with small population sizes are known to decrease establishment rates if the number of individuals introduced is too low [START_REF] Simberloff | The role of propagule pressure in biological invasions[END_REF]) and to block the spread of invaders to unoccupied patches [START_REF] Keitt | Allee effects, invasion pinning, and species' borders[END_REF][START_REF] Johnson | Allee effects and pulsed invasion by the gypsy moth[END_REF].

Because of their easily tractable effect on population dynamics, Allee effects are classically considered to account for small population dynamics [START_REF] Courchamp | Allee effects in ecology and conservation[END_REF]. However, other mechanisms affect small introduced populations. For instance, they have high risks of going extinct at random, a phenomenon known as demographic stochasticity [START_REF] Lande | Stochastic population models in ecology and conservation[END_REF].

Their probability of not producing any individual able to successfully colonize other patches is also greater, a phenomenon referred to as dispersal stochasticity (Morel-Journel et al. 2016a). These demographic mechanisms have been shown to interact with dispersal and impact establishment. Strong emigration notably tends to decrease the establishment rate of introduced populations by decreasing individual density in the introduction site [START_REF] Kean | Effects of dispersal on local population increase[END_REF][START_REF] Robinet | Dispersion in time and space affect mating success and Allee effects in invading gypsy moth populations[END_REF]Morel-Journel et al. 2016b). This study adresses the interplay between the local dynamics of introduced populations and the network structure, during the first stages of an invasion. Network structure will be characterized using two metrics related to the distribution of edges, i.e. the way patches are connected to one another: clustering and centralization.

Clustering indicates the presence of clusters, i.e. groups of vertices well connected to one another [START_REF] Watts | Collective dynamics of "small-world[END_REF][START_REF] Jordán | Characterizing the importance of habitat patches and corridors in maintaining the landscape connectivity of a Pholidoptera transsylvanica (Orthoptera) metapopulation[END_REF]. Clusters have notably been shown to hinder the spread of pathogens [START_REF] Badham | The impact of network clustering and assortativity on epidemic behaviour[END_REF], reduce the extinction risks of already colonized patches [START_REF] Kininmonth | Metapopulation mean life time within complex networks[END_REF], and buffer the propagation of disturbances across metapopulations [START_REF] Gilarranz | Effects of network modularity on the spread of perturbation impact in experimental metapopulations[END_REF]). In the context of invasion biology, they are expected to slow down the spread by "trapping" invaders. Invasions in highly clustered landscape are therefore expected to create fewer, larger populations, which could therefore also be less vulnerable to extinction because of small population mechanisms. Centralization indicates the presence of hubs, i.e. patches concentrating many dispersal flows across the landscape. Hubs are usually central elements in metapopulations, increasing connectivity and the persistence of remote populations otherwise isolated in fragmented landscapes [START_REF] Watson | Identifying critical regions in small-world marine metapopulations[END_REF]. They have been shown to increase the speed of infection spread [START_REF] Ueno | Controlling nosocomial infection based on structure of hospital social networks[END_REF][START_REF] Romano | Modeling infection transmission in primate networks to predict centrality-based risk[END_REF] and to be more likely to generate epidemic outbreaks [START_REF] Da Silva | Predicting epidemic outbreak from individual features of the spreaders[END_REF]. In invasion biology, hubs are expected to facilitate the rapid spread of invaders across the landscape. Therefore, they may also increase extinction risks of small introduced populations by decreasing the population density locally, especially if the hub is the introduction site itself.

We developed an individual-based model to simulate the first generations following introduction, in landscapes whose structure was based on networks with known centralization and clustering levels. In addition, we performed introductions of the parasitoid Trichogramma chilonis, in artificial landscapes, and followed the invasion dynamics during ten generations.

We chose the spatial structures of these artificial landscapes among networks used in the simulations to specifically test the impacts of centralization. We performed simulations in the presence of Allee effects or with demographic stochasticity alone, to investigate the interaction between landscape structure and small population mechanisms. T. chilonis does not suffer from Allee effects (Morel-Journel et al. 2016a) but exhibits a strong stochasticity in reproduction and in dispersal. Therefore, it provided a "null biological model" of an introduced population for the experiment. The combination of simulations and experiment is a powerful tool to test hypotheses, and this approach has been at the root of landmark results that have become classics in ecology, e.g. the competitive exclusion principle [START_REF] Gause | The Struggle for Existence[END_REF], chaotic population dynamics [START_REF] Cushing | Chaos in ecology: experimental nonlinear dynamics[END_REF], stochasticity in spread rates [START_REF] Melbourne | Highly variable spread rates in replicated biological invasions: fundamental limits to predictability[END_REF] or population extinction risks [START_REF] Drake | Experimental demonstration of a two-phase population extinction hazard[END_REF]. Nevertheless, invasion biologyas well as epidemiology and conservation biology, two other fields using spatial networkis heavily dominated by theoretical development and lack experimental validations of classical theoretical predictions, which are essential for the emergence of new hypotheses and the building of a comprehensive ecological theory.

Our results provide the first experimental pieces of evidence of the impacts of landscape structure on the spread of invasive species. The consistency between these results and others in different fields demonstrates that the impacts of centralization and clustering on spreading speed are robust to the ecological context. In addition, our simulations highlight the impact of Allee effects, not only on the persistence of invasive populations, but also on the colonization speed itself. Another interesting result is that high-density mechanisms, such as overcompetition, can operate even at an early stage and underpin invasion failures depending on the network structure of the landscape. Besides highlighting the importance of landscape structure to predict the outcome of invasions, our study demonstrates its interplay with internal population dynamics.

Methods

The network structure of the landscapes

Using the igraph package [START_REF] Csardi | The igraph software package for complex network research[END_REF] of the R software (R CoreTeam 2018), we generated 112618 networks with a number of vertices Nv = 10 and a number of edges Ne = 15 (see Supporting information 1 about the generation and selection of these networks). This set of networks represented all the possible different structure of connected, undirected and unlabeled networks. We computed two metrics to characterize each network: their clustering coefficient T indicating the presence of clusters, and their centralization level C indicating the presence of hubs (see Supporting information 1 about the co mputation of these indices). The value of C is based on two measures of centrality, which encompasses several, conceptually distinct measures, such as the degree, the eigenvector centrality, the betweenness and the closeness. Although all these measures were highly correlated in the set of networks considered (see Supporting information 1), we chose to consider here the closeness and the betweenness.

Both were chosen because they carry information about the structure of the paths linking vertices, the closeness accounting for the length of these paths, and the betweenness for their tendency to lie on these paths [START_REF] Freeman | Centrality in social networks conceptual clarification[END_REF]. Because of the correlation between these metrics, we based our measure of centralization on the product between betweenness and closeness, rather than considering them separately. We also identified the most central vertex of a network (called the hub thereafter) based on the product between its betweenness and its closeness. If several vertices shared the highest centrality value, the hub was selected at random among them.

Simulations

We developed an individual-based model to simulate invasions in landscape whose spatial structure was described in the previous section (see Supporting information 2 about the description of the model). This model was in discrete time, with two successive phases dispersal between patches and local population growth. In this model, individuals had no intrinsic dispersal propensity, and dispersal probability only depends on di (the degree of the vertex the individual is in) and pe (the dispersal rate of individuals in a vertex with di = 1). Local population growth was itself divided in two phases. The first concerns mating, affected by m, a density-independent mating probability, and a, a parameter describing potential mating Allee effects. The second is a reproduction phase, affected by α, the intra-specific competition, β the fecundity rate of individuals and s, the survival rate of juveniles.

Using this model, we simulated 20 time steps after a single initial introduction of 15 individuals in one patch of the landscape, for pe = 0.03, m = 0.8, α = 0.008, β = 15 and s = 0.1 (see Supporting information 4 for sensitivity analyses of these parameters). The initial population size was chosen so that extinction because of low numbers was possible but not systematic. We simulated 2x2 invasions for each network in our set, with either a = 0 (without Allee effects) or a = 2.5 (with Allee effects), and with an introduction site which was either the hub (the patch corresponding to the most central vertex) or another patch selected at random.

Experiment

To complement the simulations, we monitored artificial invasions of the egg parasitoid Trichogramma chilonis (Hymenoptera: Trichogrammatidae) in laboratory microcosm landscapes, for ten non-overlapping generations following an initial introduction of 15 individuals (see Supporting information 3 for details about the experimental setup). This model species is especially suited for our experiment, because of its small size (approx. 1mm), its short and regular developmental cycle (consistently 9 days to reach the adult stage), and its ability to parasitize its laboratory host Ephestia kuehniella (Lepidoptera: Pyralidae). The latter allowed us to focus on T.chilonis by suppressing the host's dynamic, and estimate population sizes throughout the experiment by counting the parasitized at each generation (see Supporting information 3). Moreover, we could control for impacts of the genetic makeup of the introduced populations by introducing individuals from the same inbred line, founded by a single female and maintained in the laboratory for over 100 generations.

All the microcosm landscapes were made up of 10 patches connected by 15 corridors, and their spatial structure were selected at random in two subsets of the set of networks previously established (see Supporting information 1): the "high-C" and "low-C" landscapes, whose centralization values were respectively in the 10% highest and 10% lowest of all the centralization values in the network set. We performed 63 artificial introductions equally distributed across three treatments: (i) in a patch selected at random in a low-C landscape, (ii) in a hub in a high-C landscape, (iii) in another patch than the hub in a high-C landscape. The 63 artificial landscapes were distributed equally across 7 experimental blocks, i.e. 3 replicates of each treatment in each block.

Analysis of the results

We considered three response variables: the global extinction rate, the introduction site's extinction rate and the colonization speed. The global extinction rate was computed as the proportion of invasions during which all the patches went extinct. The introduction site's extinction rate was computed as the proportion of invasions during which the introduction site went extinct at least once. Both extinction rates followed binomial distributions and were analyzed with logistic regressions. The colonization speed was computed as the ratio between the maximal number of patches colonized and the number of generations to reach this maximal extent. Considering the maximum extent allowed us to differentiate failed establishments from the start from "boom-and-bust" invasions. The colonization speed had a continuous and positive distribution clustered towards low values and was therefore analyzed with a gamma regression.

we considered a fourth response variable for the experiment: the size reached by a local population at the generation preceding its extinction. These values were distributed as an overdispersed Poisson distribution, that were analyzed using a negative binomial regression.

We used AIC comparisons to assess the respective importance of C and T as explanatory variables in the simulations. We considered a model including both variables (CT model), one with only C (C model) and one with only T (T model). We computed ΔC as the difference between the AIC of the CT model and the T model, and ΔT as the difference between the AIC of the CT model and the C model. We considered that values of ΔT or ΔC greater than 10 indicated that the support for the CT model was unequivocally higher than for the other model [START_REF] Burnham | Model selection and multimodel inference: a practical information-theoretic approach[END_REF]. To analyze the experimental results, we used generalized linear mixed models, to account for potential differences created by the experimental blocks. The experimental block was treated as a crossed random effect, as the treatments were distributed equally across all blocks.

Results

Simulations

The global extinction rate was always negligible in the simulations when a = 0 (no Allee effect), regardless of the introduction site or the landscape characteristics (Figure 1A,B). It was on average higher when a = 2.5 but was independent from clustering whether the introduction was performed in the hub (ΔT = -1.02) or elsewhere in the landscape (ΔT = -1.87). The centralization level had an impact, but only when a = 2.5 and for an introduction in the hub (ΔC = 602.37). In this case, the global extinction rate increased with C (Figure 1A).

The introduction site's extinction rate was also overall higher when a = 2.5 than when a = 0 but was impacted by centralization regardless of Allee effects when individuals were introduced in the hub (Table 1). In this case, centralization increased the extinction rate, with a large majority of introduction sites going extinct at high C levels (Figure 1C). The impact of clustering appeared only for introductions in the hub and a = 2.5 (ΔT = 9.70). In this case, it decreased slightly the extinction rate of the introduction site (Figure 1D). Centralization had a substantial positive impact on colonization speeds when individuals were introduced in the hub itself (Figure 1E), whether a = 0 (ΔC = 1578.84) or a = 2.5 (ΔC = 3250.45). Its positive impact was weaker for introductions outside of the hub, although it was still significant when a = 0 (ΔC = 130.50, Figure 1E). Clustering markedly decreased the colonization speeds when a = 0, when introductions were performed in the hub (ΔT = 865.15) or elsewhere (ΔT = 648.80). Its impact was however limited when a = 2.5 for introductions in the hub (ΔT = 29.20) and negligible for introductions outside of the hub (ΔT = -0.58). Overall, the colonization speeds were higher when a = 0 compared to a = 2.5.

The sensitivity analyses performed on the parameters of the model (see Supporting information 4) showed that, although the parameters could have an impact on the values of the response variables (the global extinction rate, the introduction site's extinction rate or the colonization speed), they did not affect qualitatively the relationships between centralization or clustering, and these response variables.

Experiment

We were not able to evidence differences between the three experimental treatments concerning the extinction rates (Figure 2), neither at the level of the whole landscape (Likelihood Ratio Test, df = 2; p = 0.344) nor at the level of the introduction site (Likelihood Ratio Test, df = 2; p = 0.747). However, the sizes of the local populations going extinct during the invasions of high-C landscapes after an introduction outside of the hub were significantly higher than those in low-C landscapes (Wald test; z = 2.325, p = 0.021) and those in high-C landscapes with introductions in the hub (Wald test; z = 2.951, p = 0.004).

The colonization speeds recorded during the experiment were overall low (way less than one patch per generation on average) but they were significantly higher in the high-C landscapes compared to the low-C landscapes (Wald test; z = 2.980, p = 0.0014), but only when the hub was the introduction site. Otherwise, there were no discernable differences in the colonization speeds between high-C and low-C landscapes (Wald test z = 0.594; p = 0.552).

Discussion

Impacts of landscape structure on spread

Both the simulations and the experiments evidenced an impact of landscape structure on spread.

The colonization of new patches by invaders was slowed down by the presence of clusters and accelerated by the presence of hubs. The effect of centralization was dominant when individuals were initially introduced in the hub itself, although it was still significant when the introduction site was chosen at random in the simulations. Our results are consistent with other studies in other fields, and thus confirm the effect of the network structure itself on spread, regardless of context.

Centrality was repeatedly shown to facilitate the transmission of infections (e.g. [START_REF] Christley | Infection in Social Networks: Using Network Analysis to Identify High-Risk Individuals[END_REF][START_REF] Ueno | Controlling nosocomial infection based on structure of hospital social networks[END_REF][START_REF] Da Silva | Predicting epidemic outbreak from individual features of the spreaders[END_REF][START_REF] Romano | Modeling infection transmission in primate networks to predict centrality-based risk[END_REF].

Epidemics starting in hubs were shown to reach greater sizes, thus underlining the influence of the hub at the start of spread [START_REF] Da Silva | Predicting epidemic outbreak from individual features of the spreaders[END_REF]), which we confirm in the context of invasions simulations and experiments. Other studies showed that the most central agents also have higher chances of being infected [START_REF] Christley | Infection in Social Networks: Using Network Analysis to Identify High-Risk Individuals[END_REF][START_REF] Romano | Modeling infection transmission in primate networks to predict centrality-based risk[END_REF]. In studies in epidemiology considering networks of populations, the most central vertices are susceptible to re-infections, even after the initial epidemic outbreak. This phenomenon, likely behind some of the greatest epidemics (e.g. [START_REF] Gómez | Network theory may explain the vulnerability of medieval human settlements to the Black Death pandemic[END_REF], also underpins major challenges for the control of invasive populations: the most central patches do not only increase the speed of invasive spread but are also more likely to be invaded.

Our simulation results concerning clustering are also consistent with epidemiological studies showing its role in limiting the speed and extent of infections [START_REF] Keeling | The implications of network structure for epidemic dynamics[END_REF][START_REF] Miller | Spread of infectious disease through clustered populations[END_REF][START_REF] Badham | The impact of network clustering and assortativity on epidemic behaviour[END_REF]. Moreover, [START_REF] Gilarranz | Effects of network modularity on the spread of perturbation impact in experimental metapopulations[END_REF] provide experimental evidence that modularity also prevents spread in a similar context: the propagation of disturbances. Modularity, which is particularly used in trophic [START_REF] Dormann | A method for detecting modules in quantitative bipartite networks[END_REF][START_REF] Beckett | Improved community detection in weighted bipartite networks[END_REF]) and genetic networks [START_REF] Fletcher | Network modularity reveals critical scales for connectivity in ecology and evolution[END_REF][START_REF] Peterman | Assessing modularity in genetic networks to manage spatially structured metapopulations[END_REF], describes the belonging of vertices to modules, i.e. subparts of the network highly connected [START_REF] Newman | Finding and evaluating community structure in networks[END_REF]. In a metapopulation context, [START_REF] Gilarranz | Effects of network modularity on the spread of perturbation impact in experimental metapopulations[END_REF] show that disturbances occurring in one module mostly affect other populations in this module, similarly to the way clusters "trap" introduced individuals in our simulation results.

Although they did not dramatically change the impacts of centralization or clustering, Allee effects reduced the colonization speed in the simulations. This result is consistent with the theoretical predictions of [START_REF] Keitt | Allee effects, invasion pinning, and species' borders[END_REF], which suggest that Allee effects can act as a supplementary hurdle to colonization, preventing the small populations at the margin of the invaded area to produce enough dispersing individuals to successfully colonizing new patches.

This impact of population dynamics on colonization speed echoes theoretical results on the variation of speed rates with density-dependent growth [START_REF] Sullivan | Density dependence in demography and dispersal generates fluctuating invasion speeds[END_REF]) and on pushed invasion [START_REF] Lewis | Allee dynamics and the spread of invading organisms[END_REF][START_REF] Roques | Allee effect promotes diversity in traveling waves of colonization[END_REF]. The low colonization speeds observed during our microcosm invasions indicated that colonization failures also occurred during the experiment, creating a "pushed-like" invasion pattern likely caused by dispersal stochasticity, as evidenced in Morel-Journel et al. (2016a) in the same experimental system.

Impacts of landscape structure on establishment

In addition to affecting the spread of invaders, the network structure of the introduction area also impacted the dynamics of the introduced populations. While the results concerning spread suggest that the same mechanisms were at play in the experiment and the simulations, extinction appears to be underpinned by different causes. In the simulations, the centralization of the landscape increased the introduction site's extinction rate when individuals were introduced in the hub. This result highlights a tradeoff between colonization and local persistence at the beginning of invasions, which is consistent with previous works in invasion biology. Theoretical [START_REF] Lewis | Allee dynamics and the spread of invading organisms[END_REF][START_REF] Kanarek | Allee effects, aggregation, and invasion success[END_REF] and empirical [START_REF] Robinet | Dispersion in time and space affect mating success and Allee effects in invading gypsy moth populations[END_REF][START_REF] Vercken | Critical patch size generated by Allee effect in gypsy moth, Lymantria dispar (L.)[END_REF]) studies underline the negative impact of dispersal soon after the introduction on the persistence of introduced populations. However, the mechanism invoked by these studies to explain these extinctions is the Allee effect, while our simulation results

show that the extinction of the introduction site can occur because of demographic stochasticity alone.

We also observed local extinctions during the experiment, although they occurred regardless of the treatment considered. If some extinctions occurring at small population sizes can be attributed to demographic stochasticity, other populations went extinct while they were close to the carrying capacity, which would rather suggest that they suffered from overcompetition. This would be consistent with the biology of T. chilonis, which can be subject to strong over-competition through superparasitism [START_REF] Suzuki | Sex allocation and effects of superparasitism on secondary sex ratios in the gregarious parasitoid, Trichogramma chilonis (Hymenoptera: Trichogrammatidae)[END_REF]. This phenomenon occurred most often for introductions outside of the hub in high-C landscapes, i.e. in poorly connected patches. Because of their low dispersal rate, individuals mostly remained in the isolated introduction site, where they rapidly suffered from superparasitism and eventually went extinct.

Complete establishment failures did not only require the extinction of the introduction site, but also that individuals fail to colonize other patches in the landscape. While we recorded extinctions of the introduction site regardless of Allee effects in the simulations, most of them did no lead to global extinctions for a = 0. Interestingly, the most important impact of Allee effects on establishment was to prevent the colonization of new patches in the landscape.

During the experiment, colonization failures resulted from an extinction of the introduction site, either because of demographic stochasticity or over-competition, combined with colonization failures because of low dispersal rate and dispersal stochasticity.

Impact of the size of the landscape

The number of vertices of the landscapes used in this study (Nv = 10) was not only smaller than the values used for epidemiology studies, but also smaller than the sizes used for network describing entire landscapes. As we aimed at describing invasions dynamics just after the introduction, we created landscapes corresponding to the direct surroundings of the introduction site only. To test the validity of the results presented in larger landscapes, we performed additional simulations, with networks with a larger number of edges (Ne varying between 20 and 40) and with a larger number of vertices (Nv varying between 20 and 100) (see Supporting information 5). The results indicate that the effects of clustering were robust to variations in the number of vertices, while the effects of centralization became weaker, although they remained qualitatively the same. The effects of these metrics were also robust to variations in the number of edges, up to a point. For high values of Ne, such network-level metrics become irrelevant, as almost every patch can be considered a hub and belonging to a cluster.

Impact of dispersal behavior of the individuals

In our model, the dispersal of individuals was treated as a diffusion process. Hence, individuals had no intrinsic dispersal propensity, and dispersed only depending on the way the introduction site was connected to the rest of the landscape. Previous results by Morel-Journel et al. (2016b) demonstrated that this approximation is valid to describe the dispersal of T. chilonis across our experimental system. Yet, we completed our results by considering emigration rates dependent on the choice of individuals, rather than the structure of the landscape (see Supporting information 6). Dispersal choice based on environmental quality would not have been relevant to our study, as the sameness of habitat quality across all patches of the landscape is a central aspect of the study. We rather focused on the impact of the presence of conspecifics, i.e. densitydependent dispersal. The results we obtained were qualitatively similar to those previously obtained on spread and on the extinction rate of the introduction site, with little impact of the strength of the density-dependence of dispersal. These additional simulations therefore suggest that the patterns observed in our simulations are robust to variations in the dispersal propensity of individuals. However, we did not test the potential role of the patch quality on this propensity, as all the patches in all our landscapes were considered having the same quality.

Besides these external factors, the dispersal behavior of individuals can also depend on intrinsic differences between individuals. During spread, dispersal is not only impacted by the genetic background of the introduced individuals, but also affects the spatial distribution of genotypes across space. This feedback loop was not considered in our model because all the individuals were considered identical, but such evolution during range expansion has been documented experimentally [START_REF] Fronhofer | Eco-evolutionary feedbacks during experimental range expansions[END_REF][START_REF] Ochocki | Rapid evolution of dispersal ability makes biological invasions faster and more variable[END_REF][START_REF] Weiss-Lehman | Rapid trait evolution drives increased speed and variance in experimental range expansions[END_REF]). Yet, we expected the impact of the genetic background of individuals or potential variations between populations to be minimal during our experiment, because of the very low genetic variability in the populations initially introduced (15 individuals from the same inbred line).

Conclusion

Predicting the fate of introduced species remains a central objective of invasion biology. This study is a first demonstration of the use of network theory in this context to characterize the structure of landscapes and predict their invasibility. We built upon previous results in epidemiology or conservation to investigate the robustness of the influence of network structure on spread patterns. In addition, our study highlights interactions between the two network metrics studiedcentralization and clusteringand small population dynamics characteristic of early stages of invasion. Among the involved small population mechanisms, Allee effects often provide a simple and elegant way to describe the positive density dependence occurring in small introduced populations, but they lack generality and empirical support [START_REF] Kramer | The evidence for Allee effects[END_REF][START_REF] Gregory | Limited evidence for the demographic Allee effect from numerous species across taxa[END_REF]. Our results demonstrate that demographic and dispersal stochasticities can create similar patterns, by bringing small, well connected populations to extinction, or by preventing colonization. Moreover, negative density-dependence and competitive interactions can also interact with landscape structure and affect the outcome of potential invasions. Table 1: Differences in AIC between the M-model (ΔC) or the C-model (ΔM) and the CM-model for the simulation results. All the values greater than 10 (indicating unequivocal higher support for the CM-model compared to the other one) are in bold.

Response variable

FiguresFigure 1 :

 1 Figures

Figure 2 :

 2 Figure 2: Mean experimental values for global extinction rates (A), the introduction site's

Acknowledgements

We gratefully acknowledge the support of the Department Santé des Plantes et Environnement from the INRA. This research was supported in part by the European Commission through the 7th Framework Program (PURE project, contract number 265865). This work is a part of the thesis of TMJ, financed by the University Nice Sophia Antipolis. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.


Niklas Felix König

Dr Abdelaziz Al Ouahabi

Dr Salomé Poyer

Laurence Charles
email: laurence.charles@univ-amu.fr

Jean-François Lutz
email: jflutz@unistra.fr

Dr Msc N F König

J.-F Lutz

Prof L Charles

A simple post-polymerization modification method for controlling sidechain information in digital polymers

Keywords: Sequence-controlled polymers, information-containing macromolecules, digitally-encoded polymers, polymer modification, solid-phase synthesis

A three-step post-polymerization modification procedure was studied for the design of digitally-encoded poly(phosphodiester)s with controllable side-groups. Sequence-defined precursors were synthesized, either manually on polystyrene resins or automatically on controlled pore glass supports, using two phosphoramidite monomers containing either terminal alkynes or triisopropylsilyl (TIPS)protected alkyne side-groups. Afterwards, these polymers were modified by stepwise copper-catalyzed azide-alkyne cycloaddition (CuAAC). The terminal alkynes were first reacted with a model azide compound and, after deprotection of the TIPS groups, the remaining alkynes were reacted with another organic azide. This simple method allows quantitative side-chain modification, thus opening up interesting avenues for the preparation of a wide variety of digital polymers. Yet, optimized protocols reported herein shall be followed. For instance, the use of iodine in iterative phosphoramidite chemistry may result in the formation of unwished iodo-acetylenes. Practical solutions to bypass that problem are reported herein.

Digital polymers are natural or non-natural information-containing macromolecules that store binary information in their monomer sequences. [1] For instance, DNA, which is the natural storage medium of genetic information, can be engineered to store digital data. [2] Alternatively, digital information can be stored in a wide variety of synthetic polymers, including poly(triazole amide)s, [3] poly(phosphodiester)s, [4] poly(alkoxyamine phosphodiester)s, [5] poly(alkoxyamine amide)s, [6] and polyurethanes. [7] Digital polymers are currently studied for different applications including data storage, [8] long-term storage [9] and anti-counterfeiting tags. [10] In all cases, information is usually written in the polymer chains using a step-by-step iterative process and is read using a sequencing technology. [1c] Numerous sequencing methods are available for deciphering the sequences of information-containing macromolecules. [11] For instance, binary information stored on non-natural polymer backbones may be potentially deciphered using tandem mass spectrometry (MS/MS), protein nanopore sequencing, solidstate nanopore sequencing and controlled depolymerization. It was pointed out that the molecular structure of non-natural polymers can be intentionally tuned to facilitate sequencing. [1c] For instance, it was observed that the presence of repeated alkoxyamine and carbamate linkages facilitates MS/MS sequencing of information-containing macromolecules. [5][6][7] However, polymers that are tailor-made for one particular sequencing technique are not necessarily optimal for another one. For instance, nanopore sequencing, [12] which is a non-destructive analytical technique, may require specific readable side-chain motifs rather than breakable main-chain motifs. In this approach, sequence-defined polyelectrolytes are subjected to an ionic current and thread through nanopores of defined diameter. The interactions of the monomers with the pore lead to voltage variations that can be monitored and related to a sequence. Yet, this technique is not trivial and DNA pore sequencing required about twenty years of optimization. [12] This is, in part, due to the fact that DNA nucleobases have closely-related structures (i.e. purine and pyrimidine heterocycles) that are difficult to distinguish from another. In comparison, the molecular bits of non-natural sequence-coded polymers could be intentionally engineered to interact to a higher or smaller extent with a pore, thus rendering the analytes more "readable" than DNA. More generally speaking, the molecular structure of monomer bits may greatly influence the storage density and readability of digital polymers. [1c] In this context, it seems crucial to develop synthetic methods that permit to tune easily side-chain information in digital polymers. Here, we describe a universal approach, which allows side-chain modification in digitally-encoded poly(phosphodiester)s.

Scheme 1. General concept herein for the selective functionalization of sequence-defined poly(phosphodiester)s: (1) CuAAC modification of the 0-bit units, (2) deprotection of the 1-bit units; (3) CuAAC modification of the 1-bit units.

Scheme 1 shows the strategy studied in this work. Several interesting routes have been described during the last few years for the synthesis of sequence-defined polymers. [13] Among them, iterative phosphoramidite chemistry, which was initially introduced for the synthesis of oligonucleotides, has been shown to be also a robust method for the preparation of non-natural sequence-defined poly(phosphodiester)s. [4a, 14] This method can be performed manually but can also be automated, thus allowing preparation of long coded chains. [4b] Phosphoramidite chemistry was therefore selected herein for the synthesis of digitally-encoded polymers of different lengths. These precursors contained two reactive bits 0 and 1 as shown in Scheme 1; one of them bears terminal alkynes and the other TIPS-protected alkynes. Although a variety of reactions can be performed on acetylenes, CuAAC was chosen for modifying these precursors. Indeed, this reaction was shown to be versatile for stepwise modification of DNA and synthetic polymers. [15] The terminal alkynes in 0-bits were first reacted with an organic azide R1-N3 to afford substituted triazoles. Afterwards, the protecting TIPS groups in 1-bits were removed and the remaining alkynes were reacted by CuAAC with another azide R2-N3. Two phosphoramidite monomers were used to build the information-containing polymers: 5a that bears two pendant terminal alkynes and 5b which contains two TIPS-protected alkyne functions (Figure 1a). The synthesis of 5a begins with the reaction of diethylmalonate and propargyl bromide (i) and is followed by a reduction to the corresponding diol 3a (ii). For 5b, the corresponding TIPS-functionalized propargyl bromide was employed in step (i). For both monomers, one alcohol function of the diols is then mono-protected by a 4,4'-dimethoxytrityl (DMT) group (iii) and the remaining alcohol is converted into a reactive phosphoramidite (iv). Both 5a and 5b were obtained in excellent yields. The mono-protected intermediates 4a and 4b were also used to prepare polystyrene supports 7a and 7b for manual iterative synthesis (Figure 1b). They were first transformed into base-labile succinidyl ester linkers 6a and 6b (v) and then coupled to an amino-methylated polystyrene resin (vi). Besides modified polystyrene resins, a commercial controlled pore glass support 8 was also used for automated synthesis (Figure 1c).

The sequence-coded precursors were synthesized by iterative phosphoramidite chemistry. Classical cycles involving DMT-deprotection (vii), alcohol-phosphoramidite coupling (viii) and oxidation of the formed phosphites into phosphates (ix) were used. In order to illustrate the versatility of our approach, some polymers were prepared manually on polystyrene supports 7a and 7b (Figure 1d), whereas some others were prepared on 8 using an automated DNA synthesizer. In all cases, readable binary sequences were implemented through the use of 5a and 5b. However, robot-made sequences contain a thymine residue that comes from the use of thymine-loaded support 8. After reaching the desired sequence and length, the 2-cyanoethyl phosphate protective groups (x) and the ester linker were cleaved (xii) (Figure 1d). All the digital polymers were characterized by negative mode high resolution electrospray mass spectrometry (ESI-HR-MS) (Table 1 and Figures S1-S8). In all cases, the targeted structures were detected as deprotonated molecules under different charge states. Yet, it should be noted that the polymers prepared manually on 7a and 7b (P1-P4 in Table 1) exhibited impurities and slight polydispersity due to missed coupling steps, whereas those prepared on 8 using the automated synthesizer (P5-P8 in Table 1) were monodisperse. This is probably due to the fact that the automated procedure employs two capping states per cycle that were not performed during the manual synthesis of short oligomers. Moreover, it should be remarked that conventional inverse-phase chromatography purification cannot be applied to these sequence-coded precursors because of the apolar character of the TIPS-protecting groups.

Table 1. ESI-HR-MS characterization of the sequence-coded precursors. [a] Supp.

Sequence Yield m/zth m/zexp P1 7a 01100110 [b-c] 53% 966.1367 [d] 966.1375 P2 7a 0110 [b-c] 73% 708.8620 [e] 708.8626 P3 7b 1001 [b-c] 68% 708.8620 [e] 708.8638 P4 7b 1001 [b] 69% 708.8620 [e] 708.8628 P5 8 T 01000001 54% 1288.9641 [e] 1288.9659 P6 8 T 01001110 01001011 70% 880.4022 [f] 880.4030 P7 8 T 01000001 01000001 01000011 49% 945.2324 [g] 945.2326

P8 8
T 01001000 01000001 01010000 01010000 01011001 01000000 01010000 01001000 01000100 1079.0947 [h] 1079.0873 [i] [a] The numbers 0 and 1 denote monomers 5a and 5b, respectively. [b] Contains missed-steps. [c] Contains iodo-acetylene units. [d] [M-3H] 3-. [e] [M-2H] 2-. [f] [M-7H] 7-. [g] [M-8H] 8-. [h] [M-20H] 20-. [i] Lack of standards at the appropriate charge state accounts for the larger error (-7 ppm) associated with measurement of this highly charged oligomer. All ions were measured at the isotopic maximum. For some samples prepared on polystyrene resins, ESI-HR-MS also evidenced the presence of several species with masses corresponding to increments of 126 Da higher than the expected mass of the poly(phosphodiester)s (Figure 2 and Figures S1-S3). Such signals are most probably due to the fact that some terminal alkynes get randomly iodinated during the iterative synthesis. Since this side-reaction is not observed for polymers prepared on 8, a logical explanation would be that traces of iodine used in (ix)

stay entrapped in the resins despite extensive washing. These traces may accumulate and react with terminal alkynes when ammonia is used in step (xii). [16] Even if acetylene iodination does not prevent post-modification by alkyne-azide cycloaddition, [17] the concomitant use of both terminal alkynes and iodo-alkynes may lead to heterogeneously-modified, thus unreadable, polymers. In order to avoid this side-reaction, a sodium thiosulfate wash (xi) was conducted to reduce residual iodine before performing (xii) (Figure 1d). Sample P4, which was prepared following that procedure, did not exhibit any iodinated species in ESI-HR-MS (Figure S4), whereas comparable sample P3, which was cleaved in the presence of iodine, contained iodo-acetylene units (Figure S3). Altogether, the results of Table 1 suggest that any binary sequence can be synthesized using 5a and 5b.

Sequential CuAAC modification was then performed on some non-iodinated sequence-defined precursors. Although modification steps (xiii), (xiv) and (xv) could be potentially directly performed on solid-supports, step (xiv) may interfere with the base-labile linkers used in the present work. Therefore, sequential modification of the cleaved polymers was investigated in solution. Figure 3 shows for example the characterization of the post-polymerization modification of resin-made polymer P4. As displayed in Figure 1e, the terminal alkynes of the resin-made precursors were first reacted by CuAAC with a model azide, namely 1-azido-2-(2-(2-methoxyethoxy)ethoxy)ethane 10 (xiii). This reaction was performed in a mixture of DMSO and tert-butanol using Cu(I)Br as copper source, tris(benzyltriazolyl-methyl)amine (TBTA) as ligand following conditions previously optimized by Carrell and coworkers for DNA modification. [15b] The resulting mono-modified polymer P4' was purified by dialysis in methanol and characterized by ESI-HR-MS (Figure S9) and 1 H NMR (Figure 3c). Both methods evidenced successful CuAAC modification. In particular, quantitative modification was supported by the complete disappearance of terminal alkyne protons signals in the 1 H NMR spectrum at 2.3 ppm and the apparition of new peaks due to triazole proton (8.1 ppm), methylene protons adjacent to the triazole rings (4.5 ppm and 2.8 ppm) and methylenoxy protons (3.75-3.4 ppm). Afterwards, the TIPS protecting groups of the 1units in P4' were cleaved with TBAF in THF solution (xiv) to obtain P4" (Figure S10). Quantitative deprotection was evidenced by the complete disappearance of the triisopropylsilyl protons at 1 ppm in 1 H NMR and the re-appearance of terminal alkyne protons at 2.3 ppm (Figure 3c). The second CuAAC modification (xv) was then achieved in the same conditions as in (xiii) using N- (2-(2-(2-(2azidoethoxy)ethoxy)ethoxy)ethyl)-2,2,2-trifluoroacetamide 11 as a model azide. After methanol dialysis, the fully modified digitally-encoded polymer P4''' was characterized by ESI-HR-MS (Figure 3a), 19 F NMR (Figure 3b) and 1 H NMR (Figure 3c). The latter technique evidenced the quantitative disappearance of terminal alkyne protons and a significant increase of the intensity of the methylenoxy protons broad peak. Successful modification was also confirmed by the appearance of a signal due to trifluoroacetamide moieties at -77 ppm in 19 F NMR. Perhaps more importantly, ESI-HR-MS evidenced the formation of polymer P4''' in which all 0 and 1 units were quantitatively modified in steps (xiii) and (xv), respectively.

In addition, incompletely-modified sequences could not be detected, thus suggesting that the stepwise CuAAC protocol is valid for modification of information-containing macromolecules. A slightly modified protocol was used for the modification of polymers prepared on 8 that were typically obtained in smaller quantities that those prepared on 7a and 7b. The first CuAAC step was not performed in solution but directly on the controlled-pore glass cartridge. Subsequently, the monomodified polymer was cleaved from the support and TIPS deprotection and the second CuAAC steps In summary, stepwise CuAAC modification is an efficient strategy for the side-chain modification of digitally-encoded polymers. Importantly, all modification steps studied herein appeared to be quantitative, thus allowing synthesis of defect-free readable sequences. Consequently, this work opens up interesting avenues for the design of "universally" readable non-natural information-containing macromolecules.

Figure 1 .

 1 Figure 1. (a) Monomer synthesis: i) NaH, (protected) propargylbromide; ii) LiAlH4; iii) DMTr chloride; iv) 2-cyanoethyl-N,N-diisopropylchlorophosphoramidite. (b) Functionalization of polystyrene resin: v) succinic anhydride, DMAP; vi) aminomethylated polystyrene, DMAP, DCC. (c) Commercial dT-CPG support modified with deoxythymidine. (d) Example of iterative synthesis on polystyrene support: vii) DMT deprotection with TCA; viii) coupling: 5a or 5b, 1H-tetrazole; ix) oxidation: I2; x) phosphate deprotection with piperidine; xi) iodine reduction with sodium thiosulfate; xii) cleavage from resin with aqueous methylamine/ammonia. (e) Example of post-polymerization modification: xiii) R1-N3 10, CuBr, TBTA, sodium ascorbate; xiv) TBAF; xv) R2-N3 11, CuBr, TBTA, sodium ascorbate.

Figure 2 .

 2 Figure 2. ESI-HR-MS spectrum (zoom of the 950-1200 m/z range) recorded for the sequence-coded precursor P1 with the sequence 01100110, detected as a triply deprotonated molecule (m/z 966.2) together with iodinated species.

Figure 3 .

 3 Figure 3. Characterization of the binary-coded poly(phosphodiester)s P4 at different stages of its modification: after solid-phase synthesis (light grey), after CuAAC transformation of the 0-units P4' (blue), after TIPS-removal on the 1-units P4" (dark grey) and after CuAAC modification of the 1-units P4''' (red). (a) ESI-HR-MS spectrum of the fully modified polymer P4'''.(b) 19 F NMR spectrum of the fully modified polymer P4'''.(c) 1 H NMR spectra recorded in methanol-d4 for P4, P4', P4", and P4'''. The displayed regions are8.5-7.5 ppm and 4.7-0.5 ppm. The cut region 7.5-4.7 ppm contains the water peak of the deuterated solvent.

  solution. Following this methodology, sample P7 with sequence T 01010000 01010000 01000100 was modified with azides 10 and 11. Figures S11 and S12 show the ESI-HR-MS and 1 H NMR characterization of the intermediates P7' and P7". Ultimately, FigureS13shows the 1 H NMR of the final polymer P7'''. The complete disappearance of terminal alkyne protons in that spectrum confirmed the successful stepwise CuAAC modification of the digital polymer

Acknowledgements J.F.L. thanks the H2020 program of the European Union (project Euro-Sequences, H2020-MSCA-ITN-2014, grant agreement n°642083), the European Research Council (ERC Proofs-of-Concept project Sequence Barcodes, grant agreement n°680097) and the CNRS for financial support. The PhD position of N.F.K. is supported by the ITN Euro-Sequences. The manager position of A.A.O. is supported by the ERC.

Entry for the Table of Contents

COMMUNICATION

In digital polymers, binary information is written using two comonomers, with slightly different molecular structures, that are defined as 0 and 1 bits. Here, we describe a simple stepwise modification method that permits to tune the molecular structure of the bits after polymerization.

Niklas


Jean-Arthur Amalian

Benoit Salomé Poyer

Eric Petit

Sofia Telitel

Valérie Monnier

Denise Karamessini

Didier Gigmes

Jean-François Lutz
email: jflutz@unistra.fr.

Laurence Charles
email: laurence.charles@univ-amu.fr

Eric Benoit

Petit

Negative mode MS/MS to read digital information encoded in sequence-defined oligo

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

Introduction

Synthetic digital polymers have recently emerged as an interesting new class of macromolecules that can be used to store information at the molecular level [START_REF] Lutz | Sequence-controlled polymers[END_REF][START_REF] Lutz | From precision polymers to complex materials and systems[END_REF]. By using a set of two co-monomers arbitrarily defined as the 0-and 1-bits of the ASCII alphabet [START_REF] Lutz | Writing on polymer chains[END_REF][START_REF] Colquhoun | Information-containing macromolecules[END_REF], our group has reported various synthesis methodologies involving different chemistries to produce so-called digitally encoded polymers such as poly(alkoxyamine amide)s [START_REF] Roy | Design and synthesis of digitally encoded polymers that can be decoded and erased[END_REF][START_REF] Roy | Convergent synthesis of digitally-encoded poly(alkoxyamine amide)s[END_REF], poly(phosphodiester)s [START_REF] Al Ouahabi | Synthesis of non-natural sequence-encoded polymers using phosphoramidite chemistry[END_REF][START_REF] Al Ouahabi | Synthesis of monodisperse sequencecoded polymers with chain lengths above DP100[END_REF], poly(alkoxyamine phospho-diester)s [START_REF] Cavallo | Orthogonal synthesis of easyto-read information-containing polymers using phosphoramidite and radical coupling steps[END_REF], poly(triazole amide)s [START_REF] Trinh | Synthesis of monodisperse sequence-defined polymers using protecting-group-free iterative strategies[END_REF][START_REF] Trinh | Preparation of informationcontaining macromolecules by ligation of dyad-encoded oligomers[END_REF], and polyurethanes [START_REF] Gunay | Chemoselective synthesis of uniform sequence-coded polyurethanes and their use as molecular tags[END_REF][START_REF] Karamessini | Identification-tagging of methacrylate-based intraocular implants using sequence defined polyurethane barcodes[END_REF]. In all these species, the chemical code was based on mass, with a simple H/CH3molecular variation between the two co-monomers. As a result, MS/MS sequencing appeared as a key technique to "read" information that was "written" in the chains. Indeed, similarly to biopolymers such as peptides [START_REF] Biemann | Laying the ground work for proteomicsmass spectrometry from1958 to 1988[END_REF], oligonucleotides [START_REF] Mcluckey | Tandem mass-spectrometry of small:multiply charged oligonucleotides[END_REF], or carbohydrates [START_REF] Domon | A systematic nomenclature for carbohydrate fragmentations in FAB-MS/MS spectra of glycoconjugates[END_REF], collision-induced dissociation (CID) of sequence-defined synthetic polymers proceeds via back-bone cleavages independently of comonomer arrangement along the chain, yielding MS/MS data that allow their sequence to be reconstructed. However, complexity of CID spectra (which affects MS/MS readability of encoded information) highly depends on the chemistry of the polymer backbone, as commonly observed for synthetic polymers [START_REF] Wesdemiotis | Fragmentation pathways of polymer ions[END_REF]. Most often, and similar to the case of biopolymers, multiple fragment series are generated upon activation, which leads to redundant information when used for sequencing purposes. Moreover, spectral complexity further increases with the size of dissociating chains: because longer chains promptly adopt high charge states upon electrospray ionization(ESI), their dissociation products also exhibit a charge state distribution. Typically, a chain-length of at least 100 residues is wished for such digital polymers, in order to prepare libraries of moderate size macromolecules to compose long messages rather than using very high molecular weight polymers, as learnt from studies involving DNA chains as a medium to store information [START_REF] Church | Next-generation digital information storage inDNA[END_REF][START_REF] Goldman | Towards practical high-capacity, low-maintenance information storage in synthesized DNA[END_REF].In this context, simplicity of MS/MS pattern obtained for small oligomers (n < 10) prepared in the earlier development stages of new sequence-defined polymers is a key requirement before dedicating further efforts in optimizing synthesis protocols to produce longer congeners with perfect monodispersity. So far, candidates that best fulfilled this requirement were polyurethanes (PUs) [START_REF] Gunay | Chemoselective synthesis of uniform sequence-coded polyurethanes and their use as molecular tags[END_REF].As depicted in Scheme 1, information-encoded PUs are built from two co-monomeric units composed of a carbamate function linked to an alkyl segment, the latter being used to implement a binary 0/1code with R = H for the 0-bit and R = CH3for the 1-bit. Cooligomeric chains are capped by a carboxypentyl-end-group and a hydroxyl-termination. Although their structure allowed coded PU oligomers to readily ionize in both polarity modes, dissociation of deprotonated molecules was found to offer an extremely simple MS/MS pattern, with a single series of fragments spaced by the mass of one or the other co-monomer (vide infra). This peculiar pattern allows any binary messages to be literally read from CID spectra. In order to figure out whether such a simple MS/MS sequencing could also be anticipated for longer chains, mechanisms underlying dissociation reactions observed for PU co-oligomers were studied in details.

Scheme 1. Structure of information-containing oligourethane with a defined sequence of 0 (R=H, m0=115.1 Da) and 1 (R=CH3, m1=129.1 Da) co-monomeric units.

Experimental 2.1. Chemicals

Digitally-encoded polyurethanes studied in this work were synthesized using an orthogonal iterative solid-phase approach described in a recent publication [START_REF] Gunay | Chemoselective synthesis of uniform sequence-coded polyurethanes and their use as molecular tags[END_REF]. The oligourethane with ancarboxylic acid moiety was prepared using ethyl 6-aminohexanoate in the last step of the iterative synthesis. The terminal ester function was then acidified during cleavage and the formed ethanol was evaporated along with cleavage reagents, thus affording the desired compound. The sample containing a brominated end-group was prepared using a reported post-oligomerization modification procedure [START_REF] Telitel | Sequence-coded ATRP macroinitiators[END_REF].All oligomer samples (a few mg) were first solubilized in methanol (SDS, Peypin, France), then diluted (1/10 3 -1/10 4 , v/v)in a methanolic solution of ammonium acetate (3 mM) purchased from Sigma Aldrich (Saint Louis, MO), prior being subjected to electrospray ionization (ESI). Poly(ethylene glycol) (PEG) and poly(methylmethacrylate) (PMMA) used as internal standards for accurate mass measurements were from Sigma Aldrich. Deuterated methanol (CD3OD, 99.80%D) used for H/D exchange experiments was from Eurisotop (Saint-Aubin, France). In these experiments, a 50 µL sample solution was evaporated to dryness and dissolved in100 µL CD3OD (3 times) prior injection in the ESI source.

Mass spectrometry

High resolution MS and MS/MS experiments were performed using a QStar Elite mass spectrometer (Applied Biosystems SCIEX, Concord, ON, Canada). Ions were accurately mass measured after internal calibration of the orthogonal acceleration time-of-flight (oa-TOF) mass analyzer, using two anionic adducts of PEG or PMMA to bracket the targeted analyte m/z value in the MS mode [START_REF] Charles | Influence of internal standard charge state on the accuracy of mass measurements in orthogonal acceleration time-of-flight mass spectrometers[END_REF]or the precursor ion signal in the MS/MS mode. In CID experiments, precursor ions were selected in a quadrupole mass analyzer, injected into the collision cell, and product ions were measured in the oa-TOF. Low resolution MS3experiments were performed usinga3200 Q-TRAP mass spectrometer (Applied Biosystems SCIEX).This instrument is equipped with two quadrupole mass analyzers separated by a collision cell, hence allowing classical MS/MS experiments. However, applying appropriate voltages at both the entrance and the exit of the second mass analyzer allows this quadrupole to be operated as a linear ion trap. As a result, MS3experiments can also be performed, with the first dissociation stage occurring in the collision cell while the second one takes place in the linear ion trap. Both instruments were equipped with an electrospray ionization source operated in the negative ion mode(capillary voltage: -4200 V; cone voltage: -75 V), using air as the nebulizing gas (10 psi) and nitrogen as the curtain gas (20 psi) and as the collision gas. Sample solutions were introduced in the ESI source with a syringe pump at a 5 µL min -1 flow rate. Instrument control, data acquisition and data processing of all experiments were achieved using Analyst software (QS 2.0 for QStar, 1.6.3 forQ-TRAP) provided by Applied Biosystems SCIEX.A Waters Synapt G2 HDMS mass spectrometer was used to perform pseudo-MS3experiments, also known as post-IMS parallelMS3-type acquisition. Typically, the quadrupole mass analyzer was used for selection of precursor ions that were further activated in an ion trap (collision energy: 35 V; argon pressure: 0.5 bar). So-obtained product ions were then all ejected from the trap, passed through a cooling cell (helium flow: 180 mL min -1 ) and reached the Travelling Wave Ion Mobility Spectrometry (TWIMS) cell operated at 2.94 mbar of N2, with a wave velocity (m s -1 )/wave height (V) couple of 1100/40. After IMS separation, fragments were activated in a collision cell (collision energy: 25 V; argon pressure: 0.5 bar) and ions generated during this second activation stage were mass measured using the oa-TOF mass analyzer in the 50-1500 m/z range. The ESI source of this instrument was operated in the negative ion mode (capillary voltage: -2.2 kV; sampling cone voltage: -50 V; desolvation gas (N2) flow: 300 L h -1 at 150•C; source temperature: 50•C). Instrument control, data acquisition and data processing of all experiments were achieved with the MassLynx4.1 programs, while 2D-plots were generated with Driftscope 2.1, both provided by Waters.

Results and discussion

Sequence-defined PUs were ionized in the negative ion mode ESI after deprotonation of the carboxylic acid group in the α moiety (Scheme 1). Single stage MS experiments confirmed that the requested uniform structure has been achieved for these sequence-defined polymers, with main detection of targeted species validating that all chains exhibit the same length (Sup-porting Fig. S1). As mentioned in the Introduction, MS/MS spectra obtained for deprotonated PU oligomers were extremely simple and hence readily allowed full coverage of their sequence, as described in details hereafter.

CID of deprotonated PUs

Collisional activation of deprotonated PU oligomers generated a single fragment series, with the distance between two consecutive congeners being equal to the mass of one or the other comonomeric unit [START_REF] Gunay | Chemoselective synthesis of uniform sequence-coded polyurethanes and their use as molecular tags[END_REF]. Such CID spectra easily allowed reading of any message in digitally encoded PUs, as exemplified with the two isomeric chains 101010 (Fig. 1a) and 000111 (Fig. 1b) of the same co-monomeric composition but distinct sequences. This very simple dissociation pattern resulted from the occurrence of a single bond cleavage in each repeating unit, that is, the O (CO) linkage of the carbamate moiety, combined with the fact that the negative charge remained located at the α chain end, hence preventing any complementary products to be detected. This greatly contrasts with the complex dissociation behavior of protonated PUs formed in positive mode ESI, observed to yield four ion series together with abundant secondary fragments (Sup-porting Fig. S2), as previously reported for related species [START_REF] Warrass | Oligocarbamates as MHC class I ligands[END_REF][START_REF] Gies | Collision induced dissociation study of ester-based polyurethane fragmentation reactions[END_REF]. According to the nomenclature established for synthetic polymer fragments [START_REF] Wesdemiotis | Fragmentation pathways of polymer ions[END_REF], product ions generated from deprotonated Pus were named ai -as they were formed upon cleavage of the first bond in the monomeric unit, with the subscripted i designating the total number of partial or entire repeating units they contain. Owing to their m/z values and as supported by accurate mass measurements (Supporting Table S1), ai -fragments were formed after a proton has been transferred from the righthand side of the dissociating center to the ether oxygen atom of the carbamate moiety. In order to find out which proton was transferred and whether ai -fragments were formed according to competitive or consecutive processes, this reaction was mechanistically studied using the 101010 oligomer as a model.

Competitive vs consecutive mechanisms.

Both consecutive and competitive processes can be envisaged to account for the O (CO) bond cleavage in the carbamate moiety. As previously reported during CID of ester-based polyurethanes in the positive ion mode [START_REF] Gies | Collision induced dissociation study of ester-based polyurethane fragmentation reactions[END_REF], proton transfer can proceed via a charge-remote reaction from the nearby amide NH group. Accordingly, ai -fragments could be formed either from the precursor ion (in a competitive manner) or upon dissociation of higher ai -congeners (in a consecutive manner). The alternative proton transfer proceeding from the terminal OH can also be envisaged to occur consecutively and competitively. However, the efficiency of the latter process is expected to decrease when considering formation of small ai -ions that would imply a long distance proton transfer. Monitoring the abundance of ai -ions as a function of the activation energy imparted to their precursor (Fig. 2a) confirmed that both competitive and consecutive processes were at work. On the one hand, a similar energy threshold of about 0.60 eV was measured for all product ions (except a1-), hence in favor of competitive processes. On the other hand, activation energy allowing maximum fragment abundance (highlighted by a square in each data series of Fig. 2a) was noticed to increase as the fragment size decreased, as typically observed for consecutive processes. The smallest ai -congener at m/z 131.1 was a noticeable exception, as its detection required significantly higher activation energy (at least 0.85 eV, center-of-mass frame) for the m/z 863.5 precursor (red plot in Fig. 2a). This particular behavior suggested that a1-was formed according to a different reaction as compared to larger fragments. Data recorded in MS3experiments (Fig. 2b-f) confirmed these assumptions. On the one hand, all ai -fragments (with i = 2-5) of the deprotonated 101010 6-mer were formed upon activation of any an -ions (with n >i), hence in a consecutive manner. On the other hand, the smallest a1-ion was never observed in any of theseMS3spectra, indicating that cleavage of the first carbamate bond only proceeded from the precursor ion. Alternatively, this last result could be explained by the fact that energy imparted to ions in the linear ion trap during the second activation stage might not be high enough to induce this particular reaction.

Mechanism for O(CO) bond cleavage in carbamate functions

H/D exchange (HDX) experiments were performed in order to identify which proton transfer induced the O (CO) bond cleavage in carbamate moieties. When prepared in CD3OD, the studied101010 6-mer was detected as a deprotonated molecule with mostly up to 6 D atoms at m/z 869.5, while full H/D exchange con-ducting to a 7 D containing species was detected to a small extent at m/z 870.5 (Supporting Fig. S3a). Similarly, HDX experiments con-ducted for different PU molecules systematically showed their m/z values to mostly shift by a number equal to their polymerization degree (data not shown). Although surprising, this result suggested that H/D exchanges mainly occurred in all carbamate groups but hardly in the hydroxyl end-group. This assumption was further supported by HDX data obtained for PU oligomers with an ester moiety instead of OH as the ω termination: as shown in Supporting Fig. S3b, m/z value measured for the ester derivative 4-mer 0001clearly indicated deuteration of the four carbamate groups. CID data recorded for the 101010-d6 oligomer at m/z 869.5 are presented in Fig. 3. It should be noted that, to perform HDX experiments, no ammonium acetate was added in PU samples prior to ESI in order to avoid H/D back exchanges in solution. As a result, depro-tonation yield of PUs was poor, and the number of accumulated spectra had to be increased by a factor of 10 in MS/MS to compensate for the quite low absolute abundance of product ions. In this MS/MS spectrum, two peaks were observed for all ai -ions, hence formed after transfer of either D (m/z in pink) or H (m/z in purple) according to mechanisms proposed in Scheme 2a-b for cleavage of the O(CO) bond. Variation of relative abundance of these two peaks as a function of the fragment size showed that formation of product ions via H-transfer became less favored as their size decreased. This result suggests that H-transfer mostly occurred from the terminal group (Scheme 2b): this would be consistent with the lower efficiency expected for such long range proton transfers and would further support that the 101010-d6 precursor ion mostly contained an OH (instead of OD) terminal group. To account for the more energy demanding reaction leading to a1 -, a particular conformation of the α termination was proposed (Scheme 2c), in which folding of the pentyl segment would allow the carboxylate moiety to strongly interact with the NH proton of the first unit. As a result, this proton could no longer be transferred to induce the cleavage of the nearby O (CO) bond according to the pathway shown in Scheme 2a. Alternatively, activation of the pre-cursor ion would lead to protonation of the carboxylate group, and subsequent electron transfers would conduct to the formation of a1 -(Scheme 2c). This mechanism is consistent with CID of 101010-d6 showing a preferential form of a1 -at m/z 132.1 as a result of deuterium transfer from the first carbamate moiety (Fig. 3). More-over, existence of this hydrogen bonded conformation might also prevent release of CO2, which was never observed from either the precursor or the product ions. However, our assumption of too low activated ai -fragments to account for the absence of a1 -in all MS 3 spectra of Fig. 3 (although these reactions could mechanistically be envisaged according to Scheme 2a) was still to be validated. For this purpose, pseudo-MS 3 experiments were performed in a Synapt G2 instrument to achieve more energetic collisions during the second activation stage. In this post-IMS parallel MS 3 -type acquisition, the m/z 863.5 precursor ion was selected in a quadrupole mass analyzer and then dissociated in a quadrupole ion trap, while a collision cell was used to activate product ions that were separated in the mobility cell. Using constant instrumental settings for the second activation stage (i.e.25 eV, laboratory frame) allowed the energy imparted to ai-fragments to increase as their size decreased, that is, 1.27 eV for m/z 748.4, 1.51 eV for m/z 619.4, 1.84 eV for m/z 504.3, 2.41 eV for m/z 375.2, and 3.33 eV for m/z 260.1. In this collision regime, the a1-fragment could indeed be generated from all higher ai -congeners, as clearly shown by signals measured for m/z 131.1 at all drift times of ai -in the 2D-plot of Fig. 4. This result hence confirmed the limited activation energy available in the linear ion trap of the instrument used to acquire MS3data of Fig. 3.

CID of longer PU chains

The PU synthesis protocol was then optimized to produce longer sequence-controlled chains. As illustrated for the 2-byte (i.e., 16 residues) PU coding for 1101001011011000 (Supporting Fig. S4a), monodispersity was achieved for this sample but the oligomer was detected as a trifluoroacetate adduct rather than a deprotonated molecule, indicating a strong decrease of deprotonation yield of the α group as the chain length increased. Nevertheless, activating this[M+TFA] -gave rise to the elimination of one trifluoroacetic acid molecule to form [M-H] -, which further dissociation yielded the same simple MS/MS pattern as smaller oligomers, with the whole ai -ion series allowing a full coverage of the digitally encoded PU sequence (Supporting Fig. S4b). Although not influencing the dissociation behavior of PUs, the chain length was observed to have a strong impact on sensitivity, due to a decreased ionization yield combined with a lower transmission efficiency of high m/z ions by the quadrupole mass analyzer used for precursor selection in MS/MS. As a result, a large number of spectra had to be accumulated (50 in the MS mode, 500 in the MS/MS mode) to obtain data shown in Supporting Fig. S4.To address this issue, the ω termination was modified to include a second carboxylic acid function. Although operated on the small 00111 coded PU, this structural modification allowed doubly deprotonated molecules to be generated in the negative ion mode ESI. As illustrated in Fig. 5, this charge state increase did not modify the CID behavior of PU, with main fragmentation reactions still observed to be O-(CO) bond cleavage in each repeating unit. Due to both end-groups being deprotonated in these di-acidic species, a pair of complementary ai -/xi -fragments was detected upon cleavage at each carbamate moiety, allowing the 00111 sequence to be reconstructed from one or the other PU termination while not introducing too much complexity in the MS/MS data at least for the studied 5-mer. This derivatization protocol still needs to further be optimized to enable production of monodisperse samples of long sequence-controlled polymers.

Conclusion.

As shown in this mechanistic study, most simple MS/MS pat-terns do not necessarily imply that dissociating ions experienced a single reaction. Cleavage of O (CO) bond in carbamate groups was indeed evidenced to proceed via three mechanisms of different efficiencies. The bottom-up approach implemented here to anticipate the CID behavior of long chains from that observed for very short oligomers showed the same dissociation reactions applied for Pus with a chain length up to 16 residues. However, this size limitation was not related to fragmentation issues but to a decreasing limit of detection of PUs as their size increased, due to a decreasing ionization yield combined to a lower transmission efficiency of high m/z ions in the precursor ion selection stage. Introduction of a second acidic moiety to increase ion charge state (and hence decrease their m/z) was shown to be a promising strategy to address this issue, as this structural modification did not modify the CID behavior of digitally encoded PUs, as demonstrated for short oligomers. It should be acknowledged, however, that this approach will allow at most doubly charged species to be obtained: owing to the upper m/z limit for precursor ion selection in the quadrupole analyzer used here, the upper limit to readable data size would be about 4 bytes of information. A more efficient strategy to store larger amounts of information in PU chains would consist of making use of language higher than binary.

Fig. 1 .

 1 Fig. 1.Negative mode ESI-MS/MS recorded at 1.10 eV (center-of-mass frame) for a)101010 and b) 000111 PU oligomers, both as deprotonated molecules (m/z 863.5), with corresponding fragmentation scheme on top. * designates secondary fragments.

Fig. 2 .

 2 Fig. 2. a) Breakdown curves of ai -fragments from [101010 -H] -at m/z 863.5. MS 3 spectra obtained after activation of b) a6 -at m/z748.4, c) a5 -at m/z619.3, d) a4 -at m/z504.3, e) a3 -at m/z375.2,and f) a2 -at m/z260.2. * designates secondary fragments and # electronic noise.

Fig. 3 .

 3 Fig. 3. Negative mode ESI-MS/MS of the 101010-d6 oligomer (shown in inset) at m/z 869.5. Colored m/z values indicate which ai -form was generated after transfer of D in carbamate moieties (pink) or H from hydroxy group (purple). * designates secondary fragments.

Fig. 4 .

 4 Fig. 4.Pseudo-MS3 data recorded on the Synapt G2 HDMS and represented as a 2D-plot, where fragments generated in the second activation stage are detected at the same drift time (hence vertically aligned) as compared to their respective precursor.

Fig. 5 .

 5 Fig. 5. Negative mode ESI-MS/MS (0.96 eV, center-of-mass frame) of the di-acid terminated 00111 PU detected as a deprotonated molecule at m/z 452.3. As indicated in the corresponding fragmentation scheme (top), fragments containing the original α moiety were named ai -(in grey) while ions containing the  termination were named as xi -(in green). * designates secondary fragments.

Acknowledgments. L.C. acknowledges support from Spectropole, the Analytical Facility of Aix-Marseille University, by allowing a special access to the instruments purchased with European Funding (FEDEROBJ2142-3341). J.F.L. thanks the H2020 program of the European Union (project Euro-Sequences, H2020-MSCA-ITN-2014, grant agreement n•642083), and the Cluster of Excellence Chemistry of Complex Systems (LabEx CSC) for financial support. The PhD positions of D. K. and B. E. P. are supported by the ITN Euro-Sequences and LabEX CSC, respectively.

Appendix A. Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ijms.2017.07.006.