Keywords: Reroutage, process number, vertex separation, largeur de chemins

Reconfiguration de routage à l'aide du jeu des gendarmes et du voleur Résumé : Le process number est le nombre de requêtes simultanément perturbées lors du phase de reconfiguration de routage dans un réseau orienté connexions. Du point de vue de la théorie des graphes, ce problème est similaire à la pathwidth, mais pas toujours égal. Déterminer ce paramètre est en général NP-complet. Dans ce rapport, nous proposons un algorithme heuristique estimant en temps polynomial le process number d'un graphe orienté. Cet algorithme a de meilleures performances que les algorithmes existant.

Introduction

When designing a wavelength division multiplexing (WDM) network, designers have to take into account expected traffic patterns, class of services proposed for failure resilience, routing algorithms, various functionality of the components, and many other parameters. Such designs are extremely hard problems to solve, in particular due to unpredictable evolution of the traffic pattern during the lifetime of the network. Also, links are generally oversized to simplify the design of the network. However, with the rapid increase of bandwidth requirements it is more and more difficult to adapt the network usage to traffic variation and so to maintain a near-optimal usage of resources.

In WDM networks without wavelength conversion, the acceptation of a new connection request is subject to the availability of a lightpath (and so resources) from end to end in the network. Whatever the routing algorithm chosen, successive addition and removal of connections may lead to a poor usage of resources. Thus new connection requests might be rejected although the network has enough resources to handle the traffic pattern, up to the rerouting of some requests. For example, in Fig. 2(a), the network is a 9 nodes grid with one wavelength and 5 requests a, b, c, d, e. A new connection request r will be rejected in Fig. 2(a) although the routing of Fig. 2(b) is possible. Therefore, it is necessary to change from time to time the routing of established connections to improve the usage of resources, and so accept more traffic.

In this paper we concentrate on the reconfiguration phase, that is the problem of switching the set of connections from current routing, R 1 , to a new pre-computed routing, R 2 .

To solve this problem, several strategies are possible. For example, one may interrupt all connections requests that have to be rerouted and after restart them with their new routes in R 2 . This strategy is very simple, but requires to store a huge amount of traffic during the reconfiguration phase and so requires large and costly buffers.

Another strategy is to reroute connection requests one after the other, as soon as destination resources (lightpath in the new routing) are available, minimizing the number of simultaneous interrupted requests during the reconfiguration. It corresponds to compute the process number (Sec. 2) in the dependency digraph [START_REF] Jose | Connection rerouting/network reconfiguration. Design of Reliable Communication Networks (DRCN)[END_REF]. This requires to compute the scheduling of the rerouting, taking into account that destination resources might be currently used by other connections. More precisely, resources assigned to request r in R 2 might be used by some request r in R 1 , thus request r has to be rerouted before r. We represent these constraints by a digraph D = (V, A), the dependency digraph [START_REF] Jose | Connection rerouting/network reconfiguration. Design of Reliable Communication Networks (DRCN)[END_REF], in which each node corresponds to a request, and there is an arc from vertex u to vertex v if v must be rerouted before u.

For example to switch from routing R 1 of Fig. 2(a) to routing R 2 of Fig. 2(b), we construct the dependency digraph D of Fig. 2(c). It has one node per connection that has to be rerouted, and so, connection e is not represented in D since it is not affected. It has an arc from a to b since connection b has to be switched before connection a. Similarly, it has arcs from b to c and from b to d since c and d must be switched before b, and so on.

When the digraph D is a dag (direct acyclic graph), the scheduling is straightforward, but in general, it may contain cycles. To break them, some requests have to be temporarily interrupted, thus removing some incident arcs in D. So, the optimization problem is to find a scheduling minimizing the number of requests simultaneously interrupted. This problem has first been considered in [START_REF] Jose | Connection rerouting/network reconfiguration. Design of Reliable Communication Networks (DRCN)[END_REF]. They have proposed a heuristic algorithm based on the minimum feedback vertex set (mfvs) of D (i.e. the smallest subset S of nodes of D such that D -S is a dag) that performs well on small instances. This problem has also been considered in [START_REF] Coudert | Rerouting requests in wdm networks[END_REF][START_REF] Coudert | Characterization of graphs and digraphs with small process number[END_REF] with a modeling in terms of cops-and-robber game [START_REF] Kirousis | Searching and pebbling[END_REF]: the process number, denoted pn. It allows to use more advance theoretical concepts and so to design more efficient heuristic algorithms. For example, the digraph of Fig. 1 is such that |mfvs(D)| = n/2 and pn(D) = 2. So it is possible to reroute corresponding requests with only two requests simultaneously interrupted.

This paper starts in Sec. 2 with the modeling of our reconfiguration problem as a copsand-robber game in the dependencies digraph, and so the definition of the process number. Then, in Sec. 3 we present some previous results linking pathwidth, vertex separation and process number. We also recall the heuristic algorithm proposed in [START_REF] Jose | Connection rerouting/network reconfiguration. Design of Reliable Communication Networks (DRCN)[END_REF]. In Sec. 4, we propose a heuristic algorithm to compute the process number of digraphs (i.e. maximum number of simultaneously interrupted requests) and the corresponding network reconfiguration strategy. Finally, in Sec. 5 we present simulations results and analyze the efficiency of our heuristic algorithm compared to others.

Modeling

It has been proved in [START_REF] Coudert | Rerouting requests in wdm networks[END_REF][START_REF] Coudert | Characterization of graphs and digraphs with small process number[END_REF] that the routing reconfiguration problem can be expressed as a cops-and-robber game [START_REF] Kirousis | Searching and pebbling[END_REF], as for the pathwidth [START_REF] Robertson | Graph minors. I. Excluding a forest[END_REF]. An interruption is represented by placing an agent on the corresponding node in D. A node is said processed when the corresponding request has been rerouted. If the node was occupied by an agent, then it can be reused. We call a process strategy a series of the three following actions (rules) allowing to reroute all requests with respect to the constraints represented by the digraph. R 1 . Put an agent on a node (interrupt a connection). R 2 . Remove an agent from a node if all its out-neighbors are either processed or occupied by an agent (restart a connection on its final route when destination resources are available). The node is now processed (connection has been rerouted ).

R 3 . Process a node if all its out-neighbors are either processed or occupied by an agent1 (connection has been rerouted because destination resources are available).

A p-process strategy is a strategy which process the digraph using p agents and the process number, pn(D), is the smallest p such that a p-process strategy exists.

Clearly, pn(D) is upper bounded by the minimum feedback vertex set (mfvs) of D, that is the smallest subset S of nodes of D such that D -S is a dag. However, this bound is very large. For example, the digraph represented in Fig. 1 is such that |mfvs(D)| = n/2, but it can be 2-processed. For that, we put an agent on node u (rule R 1 ). Then we apply rule R 3 to process node u . After, we put a second agent on node v (rule R 1 ), process node v (rule R 3 ) and then process node v (rule R 2 ). We repeat on the predecessor of v until processing of the out-neighbor w of u. Finally, we apply rule R 2 on u. Notice also that when D is a dag the scheduling is straightforward using rule R 3 from the leaves (nodes with out-degree 0), and we have pn(D) = 0.

Previous work

First of all, it has been proved in [START_REF] Coudert | Characterization of graphs and digraphs with small process number[END_REF] that vs(D) ≤ pn(D) ≤ vs(D) + 1, where vs(D) is the vertex separation of D. This implies that determining pn(D) is NP-complete in general, since it is also the case for vs(D). Furthermore, when D is a symmetric digraph, the same result holds for the underlying graph G. We also have vs(G) = pw(G), where pw(G) is the pathwidth of G [START_REF] Kinnersley | The vertex separation number of a graph equals its pathwidth[END_REF][START_REF] Robertson | Graph minors. I. Excluding a forest[END_REF], and determining the pathwidth is NP-complete in general, approximable within O(log 2 |V (G)|), and there is no polynomial-time algorithm with an absolute error guarantee of |V (G)| 1-ε for any ε > 0 [START_REF] Kirousis | Searching and pebbling[END_REF][START_REF] Bodlaender | Approximating treewidth, pathwidth, frontsize, and shortest elimination tree[END_REF]. The same holds for the process number and so for the routing reconfiguration problem considered in [START_REF] Jose | Connection rerouting/network reconfiguration. Design of Reliable Communication Networks (DRCN)[END_REF][START_REF] Coudert | Characterization of graphs and digraphs with small process number[END_REF] and in this paper, thus motivating the development of an efficient heuristic algorithm.

Notice that for specific topologies it is possible to determine the process number in polynomial time. In particular, a characterization of digraphs with process number 1 and 2 is given in [START_REF] Coudert | Characterization of graphs and digraphs with small process number[END_REF] as well as O(n+m) and O(n 2 (n+m)) respectively time complexity recognition algorithms.

When the digraph is composed of several strongly connected components (sccs), it is possible to process each scc independently, and then obtain the processing of D. For that, let {scc i } be the set of sccs of D. Let also D scc be the dag of sccs of D, thus containing one node per component of {scc i } and one arc from node i to node j iff D contains an arc from some node u ∈ scc i to some node v ∈ scc j . Then the processing of D consists in processing the sccs sequentially following the order given by D scc . Since sccs are processed independently from each other the Lemma 1 follows. Notice that the computation of {scc i } and D scc takes time O(n + m) using standard algorithms. We now recall the principle of the heuristic algorithm proposed in [START_REF] Jose | Connection rerouting/network reconfiguration. Design of Reliable Communication Networks (DRCN)[END_REF], HeurJS, to obtain the rerouting ordering of a digraph D = (V, A). It uses Lemma 1, and so we may assume that D is strongly connected. The main idea of HeurJS is to break elementary cycles. For that, it first computes the c elementary cycles of D using Johnson's algorithm [START_REF] Johnson | Finding all the elementary circuits of a directed graph[END_REF] and assigns to each node u the number α(u) of cycles it belongs to. Then it places an agent on the node u 1 of maximum weight, thus breaking α(u 1 ) cycles, updates the values of each remaining nodes and repeats if α(u i ) > 0, where u i is the new node of maximum weight. Let L be the set of nodes covered by an agent. When all cycles have been broken, the digraph D -L is a dag and so can be 0-processed. Thus it processes D -L in appropriate order and then the nodes of L. The number of agents simultaneously used is |L|, and the time complexity of HeurJS is dominated by the computation of elementary cycles which takes time O((n + m)(c + 1)) [START_REF] Jose | Connection rerouting/network reconfiguration. Design of Reliable Communication Networks (DRCN)[END_REF]. Notice that this algorithm is exponential with the number c of elementary cycles, and so with the number n of nodes. Thus, HeurJS can only permit to solve the problem for small digraphs, as it has been designed for.

To the best of our understanding, HeurJS is in fact a heuristic algorithm for mfvs applied on each scc of D, following Lemma 1. Although shown to be quite efficient in practice through a large number of experimentations [START_REF] Jose | Connection rerouting/network reconfiguration. Design of Reliable Communication Networks (DRCN)[END_REF], its result could be far from optimality (see the example of Fig. 1).

Process strategy based heuristic algorithm

We now propose a heuristic algorithm HeurCM for solving the reconfiguration problem, taking advantages of the modeling with the process number proposed in [START_REF] Coudert | Characterization of graphs and digraphs with small process number[END_REF]. This heuristic algorithm has polynomial time complexity, and with high probability improves upon HeurJS on the number of agents needed simultaneously (corresponding to the number of requests disturbed at the same time) in resulting strategy.

Flow circulation algorithm

Our algorithm is based on a flow circulation algorithm which has the objective of choosing the best candidate node to receive an agent and so to break a large set of cycles.

The principle is the following. Given a dependency digraph D = (V, A), each node u ∈ V is initially assigned a weight q 0 (u) = 1/ |V |. Then, at each round t, each node u sends (1 -µ)q t (u)/d + (u) to each of its out-neighbors, where d + (u) = |Γ + (u)| with Γ + (u) the set of out-neighbors of u in D, and µ ∈ (0, 1). Furthermore u keeps µq t (u) in order to satisfy discrete Markov chain conditions. Indeed we transform the digraph

D = (V, A) into D µ = (V µ , A µ ), such that V µ = V and A µ = A ∪ {(u, u), ∀u ∈ V }.
In other words, we add loops to each node of D. Clearly, the value of µ will influence the convergence time (i.e. number of rounds) of the flow circulation algorithm without significant impact on the final weights. Thus, we choose it close to 0, and so we will only slightly slow down the convergence time. We have for all u ∈ D µ :

q t+1 (u) = µq t (u) + v∈Γ -(u) (1 -µ)q t (v)/d + (v) (1)
where Γ -(u) is the set of predecessors of u.

After k rounds (the value of k will be discussed in Lemma 3), the node with maximum weight is the best candidate. With flow circulation algorithm, we expect that the node with maximum weight belongs to many cycles. So its removal from the digraph may break many cycles.

We will now evaluate the number k of rounds needed in the flow circulation algorithm. For that, we will first prove in Lemma 2 the convergence of the algorithm using a discrete Markov chain for which we define the set of states {u 1 , u 2 , ..., u |Vµ| } corresponding to the nodes in V µ and the transition matrix M associated to D µ . We have for all x, y ∈

V µ M [x, y] =    (1 -µ)/d + (x) when y ∈ Γ + (x) µ when x = y 0 otherwise (2)
where Γ + (x) is the set of out-neighbors of x and d + (x) = |Γ + (x)|. We also define the vector of probability q t (V µ ) = (q t (u 1 ), q t (u 2 ), ..., q t (u |Vµ| )).

q t (u i ) (1 ≤ i ≤ |V µ |)
is the probability that the state of the system is u i at step t ≥ 0.

Lemma 2. Let D µ be a strongly connected digraph with loops. Whatever q 0 (V µ ), we have lim t→+∞ q t (V µ ) = q(V µ ), where q(V µ ) is the unique stationary weights vector.

Proof. The proof relies on the discrete time Markov chain described previously, where the vector q t (V µ ) represents the probability of each state at step t ≥ 0 and M , associated to D µ , represents the transition matrix (see Equation 2).

From the Markov chains theory, we know that if the chain is irreducible and aperiodic (ergodic), then there exists a unique stationary distribution q(V µ ) whatever the initial state q 0 (V µ ) [START_REF] Meyn | Markov Chains and Stochastic Stability[END_REF].

• The irreducibility of the Markov chain is obvious because D µ is a strongly connected digraph. Indeed from each state x, there exists a positive probability to move to state y since there exists a path from x to y in D µ .

• The Markov chain is aperiodic since D µ is a strongly connected digraph with loops [START_REF] Meyn | Markov Chains and Stochastic Stability[END_REF].

Fig. 2(c) consists on a digraph D = (V, A) of 4 nodes a, b, c, d. We construct D µ = (V µ , A µ ) (Fig. 2(d)), adding loops in order to satisfy the Markov conditions described before, even if it is not necessary in this example. The probabilities on the edges in Fig. 2(d) represent the Markov transition matrix (see Equation 2). We have q t (V µ ) = (q t (a), q t (b), q t (c), q t (d)) such that q 0 (V µ ) = (0.25, 0.25, 0.25, 0.25). After 4 steps, we have INRIA q 4 (V µ ) ≈ (0.1257, 0.2473, 0.2503, 0.3767), that is to say almost the stable vector q(V µ ) = (0.125, 0.25, 0.25, 0.375). Thus node d has maximum weight.

The Markov chain admits always a stable state due to the artificial loop transition on each node. We now have to choose k sufficiently large to get the stable vector q(V µ ) with an allowed error ε, that is

ε k ≤ ε, where ε k = ||q k (V µ ) -q(V µ )|| ∞ = O(|λ 2 | k ),
denoting λ 2 the second largest eigenvalue of the transition matrix M . Thus we will obtain the node with the maximum weight.

Lemma 3. If |λ 2 | < ε 1/k ,
then flow circulation algorithm computes q(V µ ) with an allowed error ε.

Proof. From Perron-Frobenius Theorem [START_REF] Axelsson | Iterative solution methods[END_REF], we know that the convergence time from q t (V µ ) to the stationary vector q(V µ ) depends only of the second largest eigenvalue λ 2 of M (the first largest is always 1). Thus when

|λ 2 | k < ε, the stable state is obtained with ||q k (V µ ) - q(V µ )|| ∞ ≤ ε.
Since computing the second largest eigenvalue value λ 2 of M is time consuming, in our experiments we choose to fix arbitrarily k = n. For the example of Fig. 2(d), we thus have ||q 4 (V µ ) -q(V µ )|| ∞ ≤ 0.003. Another example is with n = 500 requests and an allowed error of ε = 0.01. Then, the flow circulation algorithm computes the stationary vector of the Markov chain if |λ 2 | < 0.99 (remember that the second largest eigenvalue is strictly lower than 1, and that 0.99 500 < 0.01). Thus with a very high probability, we get the right stable vector and so a good candidate node for breaking cycles. See Sec. 5.1 for more details.

Heuristic algorithm

We assume that we are given a strongly connected digraph D. Otherwise, we apply our algorithm on each sccs, according to Lemma 1. We also assume that D is different than a single node digraph (otherwise we can process it easily) and it has been transformed into D µ , but we say D for simplicity.

The heuristic algorithm, HeurCM, is described in Algorithm 1. It consists in choosing a node to place an agent using the flow circulation algorithm, and then it decomposes the digraph into sccs and repeats on each scc i . Algorithm HeurCM maintains the set of nodes covered by an agent, S, sorted by insertion dates. Also S.last is the latest inserted node that has not yet been processed. After any step of HeurCM, S may change: insertion of a new node covered by an agent (rule R 1 ), and removal of nodes that can be processed (rule R 2 ). The number of agents used by HeurCM is thus the maximum size reached by S during the algorithm, S. max.

The dependency digraph D of Fig. 2(c) is strongly connected and we apply HeurCM on it. We first apply flow circulation algorithm to choose the node to be covered by an agent, i.e. the node with maximum weight. So we choose d and add it to S (Fig. 2(c)). Since D -{d} is a dag, we can now process all remaining nodes (line 1) and finally process d. Thus HeurCM uses 1 agent in this example.

RR n°6694

Algorithm 1 HeurCM Require: A strongly connected digraph D, the set of nodes covered by an agent S and the set of active nodes A Ensure: The number of agents needed 1: Apply flow circulation on D.

2: Let u be one of the nodes of maximum weight. In case of equality we choose the closest node to S.last. 3: Place an agent on u: add it to S and remove it from A. 4: Process all the nodes whose can be processed updating S and A.

HeurCM(scc i , S), i = 1.. |{scc i }| 8: end if
With the dependency digraph of Fig. 1, andw.l.o.g., node u is chosen among the set of nodes of maximum weights (all belonging to the main cycle) and added to S. Now, u is processed and then the digraph is decomposed into sccs, each being a 2-cycle, on which we apply recursively HeurCM (line 1). Thus, HeurCM will use only 2 agents.

When choosing the node of maximum weight at line 1 of Algorithm 1, it may happen that among the set of candidates, one of them allows to process a node of S. This is typically the case when the digraph is a bidirectional path u 0 , u 1 , u 2 , u 3 , . . . , u r , u r+1 . W.l.o.g., we may assume that u 1 is chosen at line 1 and so added to S. We process node u 0 , and we repeat the algorithm on u 2 , u 3 , . . . , u r+1 . Now, candidates will be u 3 and u r (the symmetry of the graph give same weights for u 3 and u r ), but u 3 is clearly a better choice since is allows to process u 2 and after u 1 thus releasing one agent. It is why in case of equality, we choose the closest node to S.last. Finally, note that the rerouting strategy follows directly our algorithm. With HeurCM, it is sufficient to know entering and leaving dates in S and in A to know when a connection is suspended (enter in S) or switched to its new route (leave S or A).

Simulations

In this section, we first analyze further the convergence time of the flow circulation algorithm. Then, we analyze the performance of our heuristic algorithm HeurCM.

INRIA

Convergence time of flow circulation

We present here some simulations showing that the second eigenvalue λ 2 of M respects almost surely the conditions given in Lemma 3. For that, we have designed different kinds of random Markov matrices M n×n .

• A fully random matrix, that is to say for each line of M (i.e. for each node), the outdegree and the out-neighbors are chosen randomly and uniformly. Fig. 3(d) and Fig. 3(b) show the distribution of λ 2 for n = 100 and 1000 tests. For these two simulations, λ 2 is always lower than 0.7. Furthermore, Fig. 3(c) describes that with high probability, λ 2 is small.

• We have designed a random matrix with constant average out-degree d, that is to say for each line of M , d is a binomial distribution of parameters d/n and n. Then the out-neighbors are chosen randomly and uniformly. Fig. 3(d) shows the value λ 2 for each possible average out-degree for a digraph of 100 nodes. Note that for each average outdegree, 100 tests has been done. λ 2 is very small (lower than 0.5) except for very small average out-degree.

Simulations of HeurJS and HeurCM

We have implemented and analysed the efficiency of HeurJS and our algorithm HeurCM (in term of number of agents simultaneously required and in term of computation time).

• In [START_REF] Coudert | Characterization of graphs and digraphs with small process number[END_REF], a characterization of the graphs with process number 2 is done and Fig. 4(a) shows the approximate process number computed by HeurCM for these graphs. The number of agents increases very slowly with the number of nodes (requests). For a digraph of 1000 nodes, the number of agents required is almost 5.

• For a symmetric grid G n×n , Fig. 4(c) shows that the approximate process number computed by HeurCM is closed to the exact value (n + 1 if n > 2) whereas the number of agents required by HeurJS increases exponentially. Furthermore, the computation time is very smaller for HeurCM compared to HeurJS (Fig. 4(d)), because a grid is a strongly connected symmetric digraph and there is an exponential number of cycles.

• For a 20 nodes random graph, starting with an average out-degree d closed to 0, increasing it until d = 6, Fig. 4(b) shows that HeurCM requires less agents that HeurJS.

Conclusion

In this paper, we have proposed a new heuristic algorithm for the reconfiguration problem of switching the connections from one routing to another with objective of minimizing the number of connections simultaneously interrupted. This new heuristic algorithm is based on a modeling with cops-and-robbers games and improves upon previous proposal in both efficiency and computation time, as shown by our simulation results.

The next step is to design an exact and efficient exponential algorithm viable in practice for digraphs with a small number of nodes.

Figure 1 :

 1 Figure 1: Digraph D such that pn(D) = 2 and mfvs(D) = n/2.

  Initial routing R 1 of requests a, b, c, d, e.The new request r can not be satisfied.

  The new routing R 2 with requests a, b, c, d, e and r.

  Dependency digraph D for switching the routing from R 1 to R 2 . Dµ from D to apply flow circulation algorithm.

Figure 2 :

 2 Figure 2: Dependency digraph D (c) and corresponding D µ (d) to switch from routing R 1 (a) to routing R 2 (b).

  λ 2 for a number of nodes n = 100, 1000 tests, and µ = 0.25. λ 2 for a number of nodes n = 100, 1000 tests, and µ = 0.25. Repartition of λ 2 for n = 100 and 100 tests. λ 2 according to the average out-degree for graphs with n = 100 nodes and µ = 0.25. 100 tests for each average out-degree.

Figure 3 :

 3 Figure 3: Simulations: second largest eigenvalue.

Figure 4 :

 4 Figure 4: Simulations: efficiency of HeurCM.

Lemma 4 .

 4 The worst case time complexity of HeurCM is O(kn(n + m)). Proof. Each flow circulation round takes time O(k(n + m)) and in the worst case D ∩ {A} = D -{u} is a single scc.

5 :

 5 Decompose D ∩ {A} into sccs: {scc i } 6: if {scc i } is not empty then

	7:

RR n°6694

In terms of cops and robber games in undirected graphs, rule R 3 expresses that the fugitive is forced to move at each step.INRIA

http://www.inria.fr ISSN 0249-6399

Acknowledgments

This work has been partially supported by ANR JC OSERA, région PACA, European projects IST FET AEOLUS and COST 293 Graal.


Veronica Roman-Reyna

Enora Dupas

Sophie Cesbron

Guido Marchi

Sara Campigli

Mary Ann Hansen

Elizabeth Bush

Melanie Prarat

Katherine Shiplett

Melanie L Lewis

Melanie L Lewis Ivey

Joy Pierzynski

Sally A Miller

Francesca Peduto Hand

Marie-Agnes Jacques

Jonathan M Jacobs

Metagenomic sequencing for rapid identification of Xylella fastidiosa from leaf samples

Keywords: Xylella fastidiosa, metagenomics, diagnostics, short-read sequencing

à 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

Xylella fastidiosa (Xf), is a globally distributed insect-transmitted plant pathogenic bacterium, causing diseases on a large hosts range. To date, 595 plant species grouped belonging to 85 botanical families have been reported as Xf hosts [START_REF] Efsa | Update of the Xylella spp. host plant databasesystematic literature search up to 30 June 2019[END_REF], some of which are of major socio-economic interest, such as grapevine, olive, citrus, coffee and almond [START_REF] Eppo | PM 7/24 (4) Xylella fastidiosa[END_REF]. Xf colonizes the xylem vessels of plants where it forms biofilms (3) that, together with tyloses and gums produced by the plant in response to the infection (4), limit water translocation. Infected hosts display symptoms of leaf scorches and plant dieback finally followed by plant death [START_REF] Chatterjee | A cell-cell signaling sensor is required for virulence and insect transmission of Xylella fastidiosa[END_REF].

Xf was first described in and limited to the Americas but recently emerged in Europe, highlighting the pathogen's capacity to expand its host range and geographical . CC-BY-NC-ND 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint (which this version posted May 15, 2021. ; https://doi.org/10.1101/2021.05.12.443947 doi: bioRxiv preprint 4 distribution [START_REF] Eppo | PM 7/24 (4) Xylella fastidiosa[END_REF][START_REF] Pierce | The California vine disease : a preliminary report of investigations[END_REF]. The pathogen was reported in Italy in 2013, where is currently devastating Apulian olive production, then detected in France in 2015, Spain in 2016 and Portugal in 2018, both on cultivated as well as spontaneous Mediterranean plant species [START_REF] Eppo | PM 7/24 (4) Xylella fastidiosa[END_REF]. The primary control strategy for Xf diseases includes eradication of hosts; therefore, fast and accurate detection is necessary to prevent major losses to growers and further pathogen dispersal.

The diagnostic of diseases caused by fastidious pathogens like Xf is difficult. This difficulty is increased as infected plants may remain asymptomatic for very long periods of time, which is associated with low bacterial concentrations, and by an irregular distribution of the pathogens in the plants (6). It is of major interest to develop reliable and highly sensitive tools for detection and detailed identification that can be used directly on plant extracts. Current standards for Xf diagnostics primarily rely on quantitative real-time PCR (qPCR) assays to detect and sometimes identify the bacterium [START_REF] Francis | Genome-based PCR Primers for Specific and Sensitive Detection and Quantification of Xylella fastidiosa[END_REF][START_REF] Harper | Development of LAMP and real-time PCR methods for the rapid detection of Xylella fastidiosa for quarantine and field applications[END_REF][START_REF] Li | Development and systematic validation of qPCR assays for rapid and reliable differentiation of Xylella fastidiosa strains causing citrus variegated chlorosis[END_REF][START_REF] Ouyang | Enhanced Reliability and Accuracy for Field Deployable Bioforensic Detection and Discrimination of Xylella fastidiosa subsp. pauca, Causal Agent of Citrus Variegated Chlorosis Using Razor Ex Technology and TaqMan Quantitative PCR[END_REF][START_REF] Ito | Universal detection of phytoplasmas and Xylella spp. by TaqMan singleplex and multiplex real-time PCR with dual priming oligonucleotides[END_REF][START_REF] Bonants | Development and Evaluation of a Triplex TaqMan Assay and Next-Generation Sequence Analysis for Improved Detection of Xylella in Plant Material[END_REF], followed by the amplification and sequencing of two, for subspecies identification, to seven housekeeping genes (cysG, gltT, holC, leuA, malF, nuoL and petC) for is Sequence Type (ST) determination and phylogeny reconstruction (2) (Fig. 1A). Five subspecies are proposed in X. fastidiosa, ie. fastidiosa, multiplex, pauca, morus, and sandyi [START_REF] Schaad | Xylella fastidiosa subspecies: X. fastidiosa subsp piercei, subsp. nov., X. fastidiosa subsp. multiplex subsp. nov., and X. fastidiosa subsp. pauca subsp. nov[END_REF][START_REF] Nunney | Large-Scale Intersubspecific Recombination in the Plant-Pathogenic Bacterium Xylella fastidiosa Is Associated with the Host Shift to Mulberry[END_REF][START_REF] Denancé | Identification of genetic relationships and subspecies signatures in Xylella fastidiosa[END_REF]. However, whole genome analyses revealed similarities of the subspecies fastidiosa, morus and sandyi, which cluster into one clade. Moreover, genome analysis indicated high frequency of horizontal gene transfer and recombination among Xf subspecies [START_REF] Nunney | Large-Scale Intersubspecific Recombination in the Plant-Pathogenic Bacterium Xylella fastidiosa Is Associated with the Host Shift to Mulberry[END_REF][START_REF] Denancé | Identification of genetic relationships and subspecies signatures in Xylella fastidiosa[END_REF][START_REF] Marcelletti | Genome-wide comparison and taxonomic relatedness of multiple Xylella fastidiosa strains reveal the occurrence of three subspecies and a new Xylella species[END_REF].

Plant samples infected by more than one Xf strain belonging to several subspecies are not uncommon and are not easy to detect, [START_REF] Dupas | Novel Tetraplex Quantitative PCR Assays for Simultaneous Detection and Identification of Xylella fastidiosa Subspecies in Plant Tissues[END_REF][START_REF] Denancé | Several subspecies and sequence types are associated with the emergence of Xylella fastidiosa in natural settings in France[END_REF]. Nevertheless, current methods do not consider new or emerging variants resulting from pathogen genetic recombination [START_REF] Nunney | Large-Scale Intersubspecific Recombination in the Plant-Pathogenic Bacterium Xylella fastidiosa Is Associated with the Host Shift to Mulberry[END_REF]. For example, qPCR with high Cq values (>35) are considered inconclusive [START_REF] Eppo | PM 7/24 (4) Xylella fastidiosa[END_REF], making decisions about disease control difficult. A complementary tool for diagnostics is the use of Next-Generation Sequencing [START_REF] Jones | Viral Diagnostics in Plants Using Next Generation Sequencing: Computational Analysis in Practice[END_REF] (Fig. 1). Because this approach can be directly used on plant extracts, it is not biased towards known sequences and provides more information about the pathogen genome, such as virulence traits. Metagenomics, the study of genetic material from environmental samples, beyond whole genome sequencing allows for the detection of strains from several subspecies and ST at the same time from the host [START_REF] Llontop | Strain-Level Identification of Bacterial Tomato Pathogens Directly from Metagenomic Sequences[END_REF]. Recently, the use of long-read sequencing as diagnostic tool identified Xf subspecies and ST from infected samples [START_REF] Bonants | Development and Evaluation of a Triplex TaqMan Assay and Next-Generation Sequence Analysis for Improved Detection of Xylella in Plant Material[END_REF][START_REF] Faino | Nanopore sequencing for the detection and the identification of Xylella fastidiosa subspecies and sequence types from naturally infected plant material[END_REF].

In this study, we developed and tested a metagenomics pipeline using in-house shortread sequencing as a complementary approach for affordable and accurate Xf detection. We were able to use metagenomics to identify Xf to strain level in single and mixed infected plant samples, at concentrations as low as one picogram of bacterial DNA per gram of tissue. In addition, we tested naturally infected field samples from Europe and the United States. We identified Xf subspecies in samples with Cq values equal to and greater than 37, which is beyond the threshold of detection for the standard and certified qPCR methods [START_REF] Eppo | PM 7/24 (4) Xylella fastidiosa[END_REF]. Overall, we developed a robust diagnostics pipeline that could be easily extended to other pathogens and implemented for surveillance of emerging agricultural threats.

6

RESULTS

Metagenomics for diagnostics pipeline

We developed and tested a metagenomics pipeline for Xf detection and subspecies identification (Fig. 1). We tested this pipeline based on three types of DNA samples: from bacterial colonies in culture, spiked plant samples, and naturally infected plant samples (Fig. 1B). To recover and identify Xf subspecies and compare it to the already sequenced genomes, we developed a pipeline that uses six different tools and custommade databases [START_REF] Wood | Kraken: ultrafast metagenomic sequence classification using exact alignments[END_REF] (Fig. 1C). The pipeline recovers Xf reads with the software Kraken2 and a custom-made database [START_REF] Wood | Kraken: ultrafast metagenomic sequence classification using exact alignments[END_REF]. The database has user-specified genomes for Xf reads identification. The user-specified genomes belonged to Xylella (n=81), Xanthomonas sp. (n=10), E. coli (n=1) and several plant sequences from NCBI (Table S1). The database had plant sequences because some Xf genomes from NCBI contained plant genomic DNA sequences. We could not clean all 18S sequences, plant plastids, or chloroplast reads from the NCBI Xf genomes. Therefore, the plant reads in the database serve as a filter to ensure plant reads were not misidentifying as Xf reads.

After Kraken2, the pipeline de-novo assembled the recovered Xf reads into contigs with the program SPAdes [START_REF] Bankevich | SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing[END_REF]. The pipeline used the Xf contigs for four different analyses: 1) subspecies identification, 2) phylogeny reconstruction, 3) identification of the already sequenced genetically closest strains, and 4) alleles for MLST profile and virulencerelated genes determination. The pipeline used Xf contigs, and the tool SendSketch from the BBMap software to identify subspecies. Then it used Pyani and LINbase software to reconstruct phylogeny by ANI [START_REF] Tian | LINbase: a web server for genome-based identification of prokaryotes as members of crowdsourced taxa[END_REF][START_REF] Pritchard | Genomics and taxonomy in diagnostics for food security: soft-rotting enterobacterial plant pathogens[END_REF]. Next, it assigned Xf strains to each Xf contig to identify the closest strains with the tool BBSplit from BBMap software. Finally, to identify specific genes or alleles, the pipeline used local BLAST with two type of subjects: a) one subject was the reported MLST allele genes and b) the other subject was the protein sequences from genes associated with virulence.

Xf identification from in-silico prepared samples

To test the pipeline sensitivity, we used in-silico samples with target (e.g. Xf) and nontarget DNAs (e.g. non-host plant). The samples included variable amounts of non-host barley (Hordeum vulgare) sequenced reads in silico spiked with Xff CFBP 7970 reads.

We obtained a strong linear correlation between the Kraken2 results and the proportion of spiked Xf sequence reads (y =103.21x -0. 0127; R² = 1) (Fig. S1). We recovered Xf reads and assembled them as contigs using SPAdes. With the Xf contigs, we performed BLAST analysis to identify MLST alleles and virulence genes. We were able to identify one to four MLST-related gene for samples spiked with 0.5 to 2.4% Xf reads (Table S3). This result indicated that we cannot capture the full MLST gene set for ST identification with less than 2.4% Xf reads (Table S3). We calculated, for all samples, the percentage of gene similarity to the virulence-related genes (Table S4).

The percentage of gene similarity increased with the higher number of spiked Xf reads.

Samples with a lower number of Xf reads had a low genome coverage to recover and analyze complete gene sequences (Table S4, S7).

We then identified Xf subspecies using the Xf contigs. Since the in-silico samples only identified Xff reads, we expected that SendSketch assigned all contigs to Xff. However, we found that 9 to 15% of Xff contigs were instead assigned to Xfm. Based on these results, we did two additional analyses to determine the best approach to analyze the Xf subspecies composition. For the first analysis, we hypothesized that complete assembled genomes would reduce the percentage of reads assigned to other subspecies. To test this, we created a smaller Kraken2 database with 30 genomes instead of 81. These 30 Xf genomes had complete assemblies. We recovered 1 to 2% fewer Xf reads with the new database, and the subspecies distribution remained the same (data not shown). The results indicated that Xf subspecies classification is not related to the level of genome assembly. Therefore, the original Kraken2 database with 81 NCBI Xf genomes was retained for all further analyses.

For the second analysis, we manually assessed the SendSketch sensitivity to mixed infections with new in-silico samples. The new samples included a set amount of barley reads in silico spiked with variable amounts of Xff CFBP 7970 and Xfm CFBP8418 reads (Table S2). For these new in-silico samples, we recovered Xf reads with Kraken2, assembled the reads as contigs and run SendSketch to identify subspecies. When using BLAST, a certain number of Xf contigs mapped equally (100% identity) to Xff and Xfm (Xf core contigs). We observed that Xf core contigs are directly proportional to the total of Xf recovered reads and samples with a higher Xff to Xfm spiked reads ratio (Table S2, Fig. S2). Moreover, the tool SendSketch randomly assigned the subspecies to Xf contigs with 100% identity to Xff and Xfm. Consequently, we developed a manual correction to separate single from mixed infections. We only used samples with either Xfm or Xff; Xfp was not part of the analysis. The correction consists of calculating the logarithm of the Xfm : Xff contigs ratio. Corrected log-ratios from 0.081 to 0.4 are considered a mixed infection. Log-ratios below 0.08 will be a single Xff infection, and higher than 0.43 will be Xfm single infection (Fig. S2).

To evaluate the pipeline with samples free of Xf, we used extracted DNAs of two healthy plant samples and a non-Xylella controls (i.e., barley leaves infiltrated with Xanthomonas). For the artificially inoculate barley samples, Kraken2 software recovered 20 to 30% of the total reads as Xanthomonas. For all four Xf-free samples, Kraken2 recovered six to 19 Xf reads (Table S2). All these Xf reads corresponded to plant reads based on the BLAST webtool from NCBI. Based on these results, the pipeline considers a sample Xf-free when it cannot recover more than 19 Xf reads (Fig. S3).

Xf identification from isolated bacteria

To test the capacity of iSeq 100 sequencing, we used six gDNAs from isolated bacteria and two known Xf gDNAs as control. The six Xf gDNAs were isolated from Italian field samples (See Material and Methods). The two control Xf gDNA samples were Xff CFBP 7970 (CFBP 7970 iSeq100) and Xfm CFBP 8418 (CFBP 8418 iSeq100). All eight gDNA samples were sequenced with the iSeq 100 System. For all eight samples, Kraken2 recovered 99% of the total reads as Xf reads (Table S2). We assembled the Xf reads as contigs and classified them into subspecies. For CFBP 7970 and CFBP 8418 for which a genome was already available, 54-78% of the contigs corresponded to Xf core contigs, 44-65% to their respective subspecies, and 2% to the closest subspecies. On average, within the six Italian samples, 20% of the contigs corresponded to Xf core contigs, 25% to Xfm, and 2% to Xff.

The ANI values were consistent with the Xf contig abundance (Fig. 2, Table S2). All six Italian samples and CFBP 8418 iSeq100 had 99 to 100% identity to Xfm and less than 97% identity to Xff and Xfp. The control sample, CFBP 7970 iseq100, had 100% identity to Xff and less than 98% identity to Xfm and Xfp (Fig. 2, Table S2).

For strain identification, we used the program BBSplit and the Harvest suite. For each sample, we selected the top three closest strains based on the program BBSplit output.

Then, we used these closest strains and the sample to compare the number of single nucleotide polymorphisms (SNPs) with the Harvest suite. For the sample CFBP 8418 iSeq100, the closest strain with 30 SNPs was Xfm CFBP 8418 (Table S5, S6). For the sample CFBP 7970 iSeq100, the closest strain with fewer SNPs was Xff CFBP 7970.

The six Italian samples had the same three closest Xfm strains, TOS5, TOS4, and TOS14. All six samples had fewer SNPs when compared against the strain Xfm TOS4.

The three TOS strains and the Italian samples were isolated from the outbreak area of Monte Argentario, Tuscany, Italy [START_REF] Giampetruzzi | Draft Genome Sequence Resources of Three Strains (TOS4, TOS5, and TOS14) of Xylella fastidiosa Infecting Different Host Plants in the Newly Discovered Outbreak in Tuscany, Italy[END_REF].

We performed BLAST analysis to identify MLST alleles and virulence genes for all eight isolated bacteria with the assembled Xf contigs. For virulence genes, the sample CFBP 7970 iSeq100 had 100% similarity to all Xff virulence genes except for rpfE (96.4%) (Table S4). The sample CFBP 8418 iSeq100 had 100% similarity to all Xfm virulence genes except for pilB (99.8%). The six Italian samples had the same similarity percentages for all Xfm virulence genes except for hemagglutinin (95.7 to 100%). We were able to identify all virulence genes and to complete the allelic profiles for ST identification (Table S3). As expected, the ST identified for the sample CFBP 7970 iSeq100 was ST2, and the one for CFBP 8418 iSeq100 was ST6. The six Italian samples had the same ST87 number.

Xf identification from spiked plant samples

We tested the pipeline with DNA extracted from grapevine petioles and midribs artificially inoculated with known bacterial concentrations of the strain Xff CFBP 7970, Xfm CFBP 8418, or an equal mix of both strains. Kraken2 output recovered 0.01% to 74.8% of the total sequences as Xf reads (Table S2). The percentage of recovered Xf reads had a positive correlation with Log10 CFU values (R 2 = 0.9876) (Fig. 3A) and Cq values (R 2 = 0.9141) (data not shown). The pipeline detected Xf with lowest bacterial concentration tested in this study (1× 10 4 CFU/ml), equivalent to a Cq 28.85 and 1.62 pg.µL -1 .

After Xf contig assembly, we were able to identify the subspecies for all samples, and the log-ratio separated single from mixed infections (Fig. 3B,3C). The log-ratio for Xffsingle infections varied between -1.56 and -0.23. The log-ratio for Xfm-single infection was 1.15, while for mixed-strains infection was 0.08. The ANI values confirmed the Xff and Xfm subspecies for single infected samples (Fig. S4, Table S2). The mixed sample (Xff CFBP 7970 + Xfm CFBP 8418) had a higher ANI value for Xfm. This result was consistent with a higher number of Xfm contigs for the mixed-strain sample (Fig. 3C, Table S2).

Based on BBSplit results, the genetically closest strains sequenced for most of the Xffsingle infections, was CFBP 7970, followed by ATCC 35879, GV230, and TPD4 (Table S5). For Xfm-single infection, the genetically closest Xfm strains were Dixon and CFBP 8418. For the mixed infection, 90% of Xf contigs were assigned to Xfm Dixon and 6% to Xff strains.

With the assembled Xf contigs, we performed BLAST analysis to identify MLST alleles and virulence genes. We identified the ST number for three of the seven artificially inoculated samples. The mixed infected sample of grapevine inoculated with the strain CFBP 8418 were identified as ST6 (Tables S4). The grapevine sample inoculated with CFBP 7970 (10 8 CFU/ml) was ST2. We could not assign an ST the sample for grapevine sample inoculated with the strain CFBP 7970 (10 7 CFU/ml) because we only identified six of the seven MLST alleles. In contrast to MLST analysis, we detected at least two virulence genes per sample (Table S4). The sample with the lowest CFU values, CFBP 7970 (10 4 CFU/ml), had 42% similarity to Xff hemagglutinin and 41% similarity to Xfm pilQ. For the remaining Xff-single infected samples, the percentage of similarity to a single subspecies increased with the higher CFU number, which is also associated with higher genome coverage (Table S7). For the Xfm-single infection, all the virulence genes had 100% similarity to Xfm. For the mixed infection, all the virulence genes had 100% similarity to Xfm and, on average, 98% to Xff.

Xf identification from field-collected samples

Finally, we tested the iSeq 100 sequencing capacity with European and American field samples (Table S2). We used 24 samples with Cq values ranging from 21 to 40 based on Harper's qPCR assay. We used three samples that were negative based on the same qPCR assay. The DNAs from the 27 samples were extracted from six different hosts: Olea europaea, Polygala myrtifolia (France and Italy), Quercus ilex, Spartium junceum, Rhamnus alaternus, and Vitis vinifera.

Kraken2 recovered 0.004% to 1.43% of the total reads as Xf (Table S2). We assembled the Xf reads into one to 2896 contigs. We found all samples had at least one contig with at least 400 bp. The limit of Xf detection with Harper and tetraplex qPCR correspond to 30-37 Cq values [START_REF] Dupas | Novel Tetraplex Quantitative PCR Assays for Simultaneous Detection and Identification of Xylella fastidiosa Subspecies in Plant Tissues[END_REF]. Therefore, we evaluated 16 samples that either had less than 30 Xf contigs, were classified as Xf-negative, or had Cq higher than 30 [START_REF] Eppo | PM 7/24 (4) Xylella fastidiosa[END_REF]. We used each Xf contig from the 16 samples as query for a Nucleotide BLAST search using webtool from NCBI. Eleven of the 16 samples gave 100% identity to Xf genomes. Hence, these 11 samples were considered Xf-positive. All contigs from the other five samples FR1-Pm, FR1-Oe, IT6-Sj, IT11-Sj, and US1-Vv, had 100% identity to chloroplast and 18S plant sequences but none to Xf. Therefore, these five samples were considered Xfnegative (Fig. 4). With our pipeline, we were able to detect Xf in samples considered inconclusive by qPCR according to Harper's (EPPO 2019).

We then used the contigs from the 22 Xf-positive samples for subspecies classification.

Overall, the samples had 50 to 100% of contigs classified as Xf core contigs (Fig. 4).

Three French samples (FR2-Qi, FR2-Pm, FR4-Oe) and five Italian samples (IT5-Sj, IT7-Sj, IT8-Sj, IT9-Sj, IT6-Ra) had 1 to 6 contigs assigned as Xfm. The French sample FR3-Oe, seven Italian samples (IT2-Sj, IT3-Sj, IT4-Sj, IT12-Sj, IT4-Ra, IT2-Pm, IT3-Pm) and the USA sample US2-Vv, were Xfm-single infected based on the manual log correction (log-ratio > 0.43) (Table S2). The sample FR2-Oe had 17% contigs assigned to Xfp and the sample US3-Vv was Xff-single infected (log-ratio < 0.08).

Then, we determined the ANI values and Xf strain composition for the 22 Xf-positive samples. Both results were consistent with the subspecies identification. The Italian samples had 99-100% ANI to Xfm. Five French samples (FR1-Qi, FR2-Pm, FR3-Oe, FR4-Oe, FR2-Qi) had 99-100% ANI to Xfm and FR2-Oe 98% ANI to Xfp (Fig. 5A, Table S2). The sample US2-Vv had 99% ANI to Xfm, and US3-Vv had 100% to Xff (Fig. 5A).

For strain distribution, all the French, USA, and three Italian samples (IT2-Sj, IT4-Sj and IT5-Sj) had more than 40% Xf contigs had 100% identity to one strain (Table S5, Fig. 5B). Except for IT2-Sj, Italian samples had most of the contigs assigned to the three strains TOS4, TOS5, and TOS14. The sample IT2-Sj had more reads assigned to Xfm RAAR14. Overall, the tools SendSketch, Pyani and Bbsplit validated the qPCR subspecies results for field samples.

We performed BLAST analysis to identify MLST alleles and virulence genes for the 22 infected samples. For 16 of 22 Xf-positive samples, we found the percentage of gene similarity to be 26 to 100% for at least one virulence gene (Table S4). Distinct from single gene analysis, we only identified some MLST-related alleles for four samples (US3-Vv, FR2-Oe, IT2-Pm, IT3-Pm); consequently, we could not identify the ST number (Table S3).

To compare some of our results with a high-performance, deep-sequencing Illumina platform as a control, we selected nine samples for re-sequencing with the MiSeq platform using the same iSeq 100 libraries from this study. The nine samples were IT3-Pm, IT5-Sj, FR2-Oe, FR3-Oe, US1-Vv, US2-Vv, FR1-Pm, FR2-Pm, IT9-Sj (Table S8, Fig. S5). We analyzed the MiSeq sequences with our pipeline and recovered 0.005%-0.792% of total reads as Xf reads with Kraken2. We assembled the Xf reads into contigs and manually assessed all samples with less than 30 Xf contigs. The NCBI Blastn analysis indicated that the samples US1-Vv and FR1-Pm, had contigs with 100% identity to plant reads; therefore, we confirmed they were Xf-negative samples. The other seven samples were considered Xf-positive. We followed the pipeline to identify and determine subspecies, phylogeny, genetically closest sequenced genome strains, MLST profile, and virulence-related genes. The results for subspecies and phylogeny identification were the same between MiSeq and iSeq100 sequencers (Fig. S5), but there were some differences for the other three analysis results (Table S2, Table S8).

For the genetically closest sequenced genome strains analysis, all samples gave the same strain distribution as iSeq100 results, except FR2-Oe which showed Xfp OLS0478 instead of Xfp COF0407 as the most abundant strain. These two Xfp strains are phylogenetically close. For MLST analysis, we identified four more alleles for the sample IT3-Pm with the MiSeq platform than with iSeq100 (Table S3), while we only detected two alleles in the sample US2-Vv sequenced with Miseq. We were only able to detect MLST alleles for the sample FR2-Oe with the iSeq100 platform. We were able to calculate the percentage of gene similarity for more virulence genes with the MiSeq platform than with the iSeq100 platform. The variation between Illumina platforms was not consistent.

DISCUSSION

In this study, we developed a user-friendly metagenomic pipeline to identify and determine Xylella fastidiosa subspecies from field-collected samples without the need for pathogen isolation. We demonstrated the flexibility of the pipeline by using seven different plant hosts and three DNA extraction methods. We recovered and assembled Xf reads into contigs from total DNA samples. We used percentage of similarity to a single subspecies to identify Xf subspecies and validated the results through phylogeny and strain proximity. Finally, we examined potential virulence-related genes among all sequenced samples.

To recover Xf reads from field samples with the tool Kraken2, we used Xf genomes available on NCBI. We found plant plastid reads in all genomes obtained from pure Xf cultures, except for Xff CFBP7970. We decided to add plant reads in the database to filter out potential plant contamination. We still recovered Xf reads for some plant samples reported as Xf-negative. Therefore, we manually examined the contigs of samples with less than 30 Xf contigs to determine if they have a low Xf concentration or are negative samples. More than the Cq values, the contig evaluation will be necessary to determine if a sample is truly Xf -negative.

We observed that the Kraken2 tool not only recovered Xf reads but also identified Xf subspecies. We decided not to use Kraken2 to identify subspecies because we found that the recovery is affected by incomplete NCBI Xf genomes subspecies information. Kraken2 uses by default the NCBI taxonomy to classify reads, if the genomes used to build the database does not have a subspecies information it will keep most of the reads at the species level. To improve the subspecies resolution, we decided to use contigs instead of reads and the tool SendSketch (BBMap tool). Contigs or assembled reads increases the coverage and reduces false-positive reads [START_REF] Ayling | New approaches for metagenome assembly with short reads[END_REF]. We used Sendskecth because it uses MinHash algorithm to be fast and it takes into account whole genomes and do not uses taxonomy.

After we identified samples as Xf-positive, we defined Xf subspecies. We observed that some samples had mapped contigs to both Xff and Xfm. These are contigs most likely associated with core sequences as only 3% of Xff and Xfm genomes are different. However, it is also possible to have a percentage of annotation error due to sequencing contamination [START_REF] Nunney | The Importance of Multilocus Sequence Typing: Cautionary Tales from the Bacterium Xylella fastidiosa[END_REF][START_REF] Kircher | Double indexing overcomes inaccuracies in multiplex sequencing on the Illumina platform[END_REF]. To determine if the samples were single-or mixed-infected, we corrected the results by calculating the Log Xfm: Xff contigs ratio (see Materials and methods). Once we defined the Log values for Xff, Xfm, and Xfm+Xff (mixed) infections, we validated the presence of single Xfm infections in all the tested European samples.

Identifying Xf subspecies in ornamental and crop plants is essential for the correct application of eradication measures or for plant movements within the EU territory according to regulation (EU) 2020/1201 [START_REF] Nunney | Recent Evolutionary Radiation and Host Plant Specialization in the Xylella fastidiosa Subspecies Native to the United States[END_REF].

The number of reads generated by the iSeq100 sequencer highlighted some limitations of the pipeline. For example, the ST were not determined for any field sample because we could not recover the complete sequences of the seven genes. This is probably associated with low genome coverage. The low number of reads also hampers deep SNP diversity and intersubspecific homologous recombination analyses. Other sequencing systems, with higher reads output than the iseq100, can also be used with this pipeline as the input is fastq files. With a high number of reads, the pipeline will provide better resolution to recover the MSLT genes. Some of the diagnostic tools for Xf diagnostics are qPCR and Sanger sequencing.

These tools require amplification of known Xf genome regions, but as is the case of MLST, do not consider new or emerging variants. Moreover, these tools introduce bias due to primer design, have unresolved results with high Cq values, and may take longer since they follow a multistep process. Our pipeline complements these conventional tools by obtaining metagenomic data directly from symptomatic or asymptomatic samples and increasing the detection power. We found some discrepancies between the number of recovered Xf reads and Cq values. These differences could be caused by PCR inhibitors or the genomic target region for qPCR that underestimates the bacterial concentration [START_REF] Modesti | Experience of the Latium region (Central Italy) as a pest-free area for monitoring of Xylella fastidiosa: distinctive features of molecular diagnostic methods[END_REF][START_REF] Schrader | PCR inhibitorsoccurrence, properties and removal[END_REF]. Metagenomics sequencing is becoming a more affordable and faster approach for diagnostics. For example, the whole detection/identification with qPCR and MLST scheme could have an estimated cost of 52-54USD per sample and takes three to four days to detect one to seven genes. With the iSeq100, it could cost 50-70USD (when having 12 samples in the same run) and take two days but while also allowing to get a complete genomic analysis of the plant and its pathogenic and commensal microbiota.

In conclusion, our pipeline provides Xf taxonomy and functional information for diagnostics without extensive knowledge of the host or pathogen. The pipeline databases used for the analysis can be public repositories or privately collected gDNA and could be adapted by the user and tailored to different plant pathogens. The sequencing can be adapted to be an in-house system as the library preparation and sequencing are user-friendly and not limited by the DNA quality or quantity. The analysis can be adjusted to detect several pathogens simultaneously. We hope this pipeline can be used for early detection of Xf or other crop pathogens and incorporated as part of management strategies.

MATERIAL AND METHODS

Xylella fastidiosa strains

The strain Xf subsp. fastidiosa (Xff) CFBP 7970 isolated in the United States (Florida) in 2013 from Vitis vinifera and Xf subsp. multiplex (Xfm) CFBP 8418 isolated in France in 2015 from Spartium junceum were provided by the French Collection of Plant-Associated Bacteria (CIRM-CFBP (CIRM-CFBP https://www6. inra. fr/cirm_eng/CFBP -Plant-Associated-Bacteria) and used as controls for whole-genome sequencing. Both strains were cultivated on modified PWG medium (Gelrite 12 g. L -1 ; soytone 4 g. L -1 ; bacto tryptone 1 g. L -1 ; MgSO4. 7H2O 0, 4 g. L -1 ; K2HPO4 1. 2 g. L -1 ; KH2PO4 1 g. L -1 ; hemin chloride (0. 1% in NaOH 0. 05 M) 10 ml. L -1 ; BSA (7. 5%) 24 ml. L -1 ; L-glutamine 4 g. L -1 ) at 28°C for one week.

Artificially inoculated plant samples

For artificial inoculations, 10 ml of a calibrated CFBP 7970 or CFBP 8418 strain suspension was spiked in 2 g of detached V. vinifera leaves. Sterile water was used for negative controls. The DNA extraction was performed using a CTAB-based extraction protocol [START_REF] Eppo | PM 7/24 (4) Xylella fastidiosa[END_REF] with slight modifications in order to concentrate bacterial DNA. After a 20min centrifugation at 20,000 g of the sample, the pellet was resuspended in 1ml of CTAB buffer. At the end of the extraction protocol, the pellet was resuspended in 50 µl of sterile demineralized water. Xf presence in the infected samples was checked using Harper's qPCR assay [START_REF] Harper | Development of LAMP and real-time PCR methods for the rapid detection of Xylella fastidiosa for quarantine and field applications[END_REF].

Plant material and bacterial gDNA

Healthy plant material of Vitis vinifera (2 g of leaf petioles) was spiked with 10 ml of a calibrated CFBP 7970 or CFBP 8418 strains suspension. Sterile water was used as negative control. The DNA extraction was performed using a CTAB-based extraction protocol [START_REF] Eppo | PM 7/24 (4) Xylella fastidiosa[END_REF] with slight modifications in order to concentrate bacterial DNA. After a 20-min centrifugation at 20,000 g of the plant macerate sample was resuspended in 1ml of CTAB buffer. At the end of the extraction protocol, the pellet was resuspended in 50 µl of sterile demineralized water. Xf presence in the infected samples was checked using Harper's qPCR assay [START_REF] Harper | Development of LAMP and real-time PCR methods for the rapid detection of Xylella fastidiosa for quarantine and field applications[END_REF].

Naturally infected samples were collected in Europe and the USA. Symptomatic samples of Olea europaea, Polygala myrtifolia and Quercus ilex were collected in October 2018 in Corsica (France) and in September 2019 in the French Riviera. Xf detection and DNA extraction of whole infected-plant tissue were performed as mentioned above. Xf subspecies were identified using the tetraplex qPCR [START_REF] Dupas | Novel Tetraplex Quantitative PCR Assays for Simultaneous Detection and Identification of Xylella fastidiosa Subspecies in Plant Tissues[END_REF].

Twig tissues, leaf petioles or green shoots of Rhamnus alaternus, Spartium junceum and P. myrtifolia growing in the Xf outbreak zone of Monte Argentario (Grosseto, Tuscany, Italy) were collected during 2019 and 2020 [START_REF] Marchi | First detection of Xylella fastidiosa subsp. multiplex DNA in Tuscany (Italy) | Phytopathologia Mediterranea[END_REF]. Xf was detected using Harper's qPCR assay [START_REF] Harper | Development of LAMP and real-time PCR methods for the rapid detection of Xylella fastidiosa for quarantine and field applications[END_REF]. For samples with a Cq value lower than 30, Xf isolation was attempted on Buffered Charcoal Yeast Extract (BCYE) agar according to PM 7/24-4 (2).

Bacterial isolates that became visible to the unaided eye within three days of incubation at 28°C were discarded; those that became visible thereafter were streaked twice for purity on BCYE agar and identified as Xf based on qPCR results [START_REF] Harper | Development of LAMP and real-time PCR methods for the rapid detection of Xylella fastidiosa for quarantine and field applications[END_REF]. Reactions were carried out after boiling bacterial suspension for 10 min. The DNA of one isolate among those that tested positive by qPCR from each plant, was extracted using the CTAB based protocol and further characterized to subspecies and ST level following the Multi Locus Sequence Typing (MLST) approach [START_REF] Yuan | Multilocus sequence typing of Xylella fastidiosa causing Pierce's disease and oleander leaf scorch in the United States[END_REF]. The GoTaq probe qPCR Master Mix (Promega, A6102) and GoTaq G2 (Promega, M784B) polymerase were used for qPCR and conventional PCR experiments, respectively. The bacterial DNA of three isolates from R. alaternus (IT1-Ra to IT3-Ra), one from S. junceum (IT1-Sj), one from P. myrtifolia (IT1-Pm) and one from Prunus dulcis (IT1-Pd), were sequenced in this study.

The assembled genomes of the six isolates were deposited on NCBI under the Bioproject PRJNA728043. Approximately 50-100 mg of grape leaf or petiole tissue was excised from each sample using a razor blade. The excised tissue was transferred to lysing Matrix A tubes (MP, 6910-500) and ground using a FastPrep 24 (MP, 116004500). DNA was extracted using ISOLATE II Plant DNA extraction kit (Bioline, 52070) following manufacturer's recommendations and CTAB lysis buffer. For Xf detection, qPCR Harper's was performed on the StepOnePlusTM system (Life Technologies, 4376600) with Sensi-FAST Probe Hi-ROX qPCR kit (Bioline, 82005) [START_REF] Harper | Development of LAMP and real-time PCR methods for the rapid detection of Xylella fastidiosa for quarantine and field applications[END_REF].

Pipeline controls

To test the pipeline, seven samples were used as negative controls. The controls were two DNA samples from healthy barley and wild grass leaves that were grown in a greenhouse, two DNA samples from barley leaves were infiltrated with a bacterial suspension (10 8 CFU/ml) of Xanthomonas translucens pv. translucens UPB886. Three healthy samples were collected in France in 2020 and in the USA in 2019. Petioles and midribs were collected from healthy Olea europaea plants in a non-Xf infected area (Angers, France) and from Polygala myrtifolia plants that were purchased form a local nursery. V. vinifera leaves were collected from the vineyard in Greene County (Va, USA).The DNA from these samples was extracted using CTAB method as mentioned above. The absence of Xf in the healthy plants was confirmed using Harper's qPCR assay [START_REF] Harper | Development of LAMP and real-time PCR methods for the rapid detection of Xylella fastidiosa for quarantine and field applications[END_REF].

For in-silico pipeline controls, two types of positive controls were used. To validate the detection of different concentrations of Xf, barley fasta sequence files were in-silico mixed with reads of the Xff CFBP 7970 sequenced in this study. The final proportion of Xff reads in the sample ranged from 0.2 to 2.4% of total reads. To validate the detection limits for single and mixed infections, the barley fasta sequence files were in silico mixed with different proportion of Xff CFBP 7970 and Xfm CFBP 8418 reads, to get from % CFBP7970 : 99% CFBP8418 to 99% CFBP7970 : 1% CFBP8418. iSeq 100 sequencing iSeq 100 sequencing libraries were prepared according to the Illumina reference guide for Nextera DNA Flex Library Prep and Nextera DNA CD Indexes. In brief, 200 to 500 ng of DNA were quantified by spectrophotometry and used for library preparation. Then, the libraries were diluted to have the same starting concentration prior to sample pooling. Eight to 12 libraries were mixed together, and 1 nM of pooled library was used per run. The sequencing settings were paired-ended (PE) read type, 151 read cycles and eight index cycles. In the iSeq 100 System, the illumina GenerateFASTQ Analysis Module for base calling and demultiplexing was selected. After sequencing for 17 h, fastq paired end read files were extracted from the machine for subsequent analysis.

Pipeline for Xylella sp. detection, classification and quantification via metagenomic analysis

Fastq files and the program Kraken2 were used to recover Xylella reads [START_REF] Wood | Kraken: ultrafast metagenomic sequence classification using exact alignments[END_REF]. The Kraken2 command options were: --paired, --minimum-hit-groups 5, --report and -db.

The database was created with 92 NCBI genomes: 79 from Xf, two from Xylella taiwanensis, ten from Xanthomonas sp. and one from Escherichia coli (Table S1). The 81 Xylella genomes were used to recover Xylella reads. The last 11 genomes were added to remove reads common to Proteobacteria that might give a false positive match. The tool SendSketch with nt server was used to make sure the 81 NCBI Xf genomes were not contaminated with plant reads (BBMap -Bushnell B.sourceforge.net/projects/bbmap/; October 30, 2019). Forty-nine NCBI sequences from plant 18S and chloroplast were added to the Kraken2 customized database to avoid extracting reads annotated as plant reads (Table S1).

The reads classified as Xylella were extracted with the script extract_kraken_reads from the KrakenTools suite (GitHub jenniferlu717/KrakenTools). The extracted Xf reads were used for downstream analysis. First Xf reads were de-novo assembled with the software SPAdes (23) using defaults settings and the option --only-assembler. Second, the Xf contigs were the query sequences in the Basic Local Alignment Search Tool website (NCBI) to confirm if they were Xf reads, or misclassified plant reads. The Blastn parameters were Nucleotide collection (nr/nt) as database and megablast program selection.

The Xf-positive contigs were used in four different analyses: 1) to identify subspecies, 2) to reconstruct phylogeny, 3) to identify the genetically closest strains already sequenced and 4) to identify alleles from specific genes and the MLST profile.

To identify subspecies, the tool SendSketch was run with the parameters, mode=sequence, records=2, and format=3, minani=100, minhit=1, address=ref, level=0 (BBMap -Bushnell B.sourceforge.net/projects/bbmap/). Only contigs with 100% average nucleotide identity (ANI) were used to identify Xf subspecies. The Xf contigs with no hits were considered as core sequences. For visualization, results were plotted using stacked bars.

To reconstruct phylogeny, ANI was calculated using the software Pyani (v. 0. 2. 10) and LINbase [START_REF] Tian | LINbase: a web server for genome-based identification of prokaryotes as members of crowdsourced taxa[END_REF][START_REF] Pritchard | Genomics and taxonomy in diagnostics for food security: soft-rotting enterobacterial plant pathogens[END_REF]. For Pyani, the option -ANIm was set and for LINbase the "Identify using a gene sequence" was set as identification method. The R package ComplexHeatmap was used to visualize the ANI cluster analysis with the parameters clustering_distance_rows = robust_dist and clustering_method_rows = "average".

Robust_dist was a function suggested by the ComplexHeatmap tutorial.

To identify the genetically closest already sequenced Xf strains, the tool BBSplit was run with Xf-contigs (BBMap -Bushnell B.sourceforge.net/projects/bbmap/). The tool used 81 Xf genomes from NCBI and the options minratio=1 ambig=best. For comparisons, each sample was normalized to their total Xf contigs and plotted using the R package ComplexHeatmap. For the isolated bacterial genomes, the most abundant strains were used as reference to identify genomics variants using Harvest suite tools [START_REF] Treangen | The Harvest suite for rapid core-genome alignment and visualization of thousands of intraspecific microbial genomes[END_REF].

To determine specific gene alleles, percentage of identity was calculated using Xf contigs as query and the Blastn algorithm (Nucleotide-Nucleotide BLAST 2. 8. 1+). The database contained complete nucleotide sequences for all alleles for the seven genes used for ST identification (cysG, gltT, leuA, malF, nuoL, holC and petC) [START_REF] Yuan | Multilocus sequence typing of Xylella fastidiosa causing Pierce's disease and oleander leaf scorch in the United States[END_REF]. All 147 alleles were downloaded from the website PubMLST (36) (Last updated: 2019-03-06).

To determine the presence of reported Xf virulence-related or common to several plant pathogenic bacterial genes [START_REF] Da Silva | Fastidian gum: the Xylella fastidiosa exopolysaccharide possibly involved in bacterial pathogenicity[END_REF][START_REF] Chatterjee | Living in two Worlds: The Plant and Insect Lifestyles of Xylella fastidiosa[END_REF], percentage of similarity was calculated using Xf contigs as query using Blastx algorithm (Nucleotide-Nucleotide BLAST 2. 8. 1+. The databases contained complete amino acid sequences for gumBCDE, pilBMQTVW, rpfCEFG, tolC, 6-phosphogluconolactonase (pgl), and hemagglutinin from the Xfm M12 and Xff M23 NCBI genomes.

MiSeq sequencing

To validate the iSeq 100 results, the same iSeq 100 libraries were used for MiSeq deep sequencing. Nine samples were selected, at least one from each iSeq 100 run, making sure not to use the same i5 and i7 tags. These nine libraries were sent to the Animal Disease Diagnostic Laboratory (Ohio Department of Agriculture, Reynoldsburg, Ohio) for sequencing. Library preparation was performed using an Illumina DNA Flex kit, and 2x250 sequencing was performed on the MiSeq platform using V3 chemistry. The pipeline described above was used to analyze the MiSeq fastq files. fastq reads from all samples were then used for metagenomic analysis. We created a database to recover Xf reads using Kraken2. The database contained Xylella, Xanthomonas and Escherichia coli genomes. We also added plant plastid genomes to remove false positive results. Xf reads were recovered from the fastq files. The Xf recovered reads were de-novo assembled to obtain Xf contigs, using SPADes. The Xf contigs were used in four different analyses, subspecies identification, phylogeny reconstruction, identification of the genetically closest strain with a sequenced genome and alleles from specific genes. To determine subspecies, we used the tool SendSketch. To reconstruct phylogeny, we calculated ANI using Pyani and the website tool LINbase (https://linbase. org/). To determine the genetically closest known Xf strain, we detected the number of hits to each Xf strain using BBSplit. To identify specific genes alleles, we calculated the percentage of identity to the seven MLST genes (cysG, . CC-BY-NC-ND 4.0 International license available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

The copyright holder for this preprint (which this version posted May 15, 2021. ; https://doi.org/10.1101/2021.05.12.443947 doi: bioRxiv preprint 36 gltT, holC, leuA, malF, nuoL , petC), and the percentage of similarity to 17 virulencerelated proteins using local BLAST+. Graphics were created with BioRender.

  Vitis vinifera DNA samples were received from the Virginia Tech Plant Disease Clinic (VA, USA). All samples were collected from Virginia vineyards in 2019. US2-Vv sample was collected from the vineyard in Greene County and US3-Vv sample was from a vineyard in Isle of Wight County. The DNA extraction and Xf detection protocols are based on work instructions from the Virginia Tech Plant Disease Clinic (VTPDC).

Fig. 1 .

 1 Fig. 1. Metagenomics for diagnostic pipeline. A) Sequenced-based detection. Two

Fig. 2 :

 2 Fig. 2: Phylogenetic reconstruction of isolated bacteria used in this study. The cluster

Fig. 3 :

 3 Fig. 3: Spiked samples with different dilutions and mixed samples. A) scatter plot

Fig. 4 .

 4 Fig. 4. Xf from field-collected samples and read mapping to subspecies from database.

Fig. 5 :

 5 Fig. 5: Metagenomics analyses of field samples identifies bacterial subspecies. A) The

available under a was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint (which this version posted May 15, 2021. ; https://doi.org/10.1101/2021.05.12.443947 doi: bioRxiv preprint

ACKNOWLEDGEMENTS

The authors are grateful for funding support from the Ohio Department of Agriculture Specialty Crops Block Grant (AGR-SCG-19-03) and USDA NIFA FACT (2021-67021-34343) to JMJ. Sara Campigli was financed by a grant from the Phytosanitary Service of the Tuscany Region (Italy). We would like to thank The Ohio Supercomputer Center for providing High Performance Computing resources.

We would also like to thank Dr. Stephen Cohen, Dr. Jeff Chang, and Dr. Alexandra Weinsberg for their valuable input for the experiment development and editing.