Mohamed Wahbi
email: mohamed.wahbi@emn.fr

Redouane Ezzahir
email: red.ezzahir@gmail.com

Christian Bessiere
email: bessiere@lirmm.fr

Asynchronous Forward Bounding Revisited

The Distributed Constraint Optimization Problem (DCOP) is a powerful framework for modeling and solving applications in multi-agent coordination. Asynchronous Forward Bounding (AFB BJ) is one of the best algorithms to solve DCOPs. We propose AFB BJ + , a revisited version of AFB BJ in which we refine the lower bound computations. We also propose to compute lower bounds for the whole domain of the last assigned agent instead of only doing this for its current assignment. This reduces both the number of messages needed and the time future agents remain idle. In addition, these lower bounds can be used as a value ordering heuristic in AFB BJ + . The experimental evaluation on standard benchmark problems shows the efficiency of AFB BJ + compared to other algorithms for DCOPs.

Introduction

Distributed Constraint Optimization Problem (DCOP) is a powerful framework to model a wide range of applications in multi-agent coordination such as distributed scheduling [START_REF] Maheswaran | Taking DCOP to the real world: Efficient complete solutions for distributed multi-event scheduling[END_REF], distributed planning [START_REF] Bonnet-Torrés | Multiply-constrained dcop for distributed planning and scheduling[END_REF], distributed resource allocation [START_REF] Petcu | A Value Ordering Heuristic for Distributed Resource Allocation[END_REF], target tracking in sensor networks [START_REF] Modi | ADOPT: Asynchronous Distributed Constraint Optimization with Quality Guarantees[END_REF] distributed vehicle routing [START_REF] Léauté | Coordinating Logistics Operations with Privacy Guarantees[END_REF], etc. A DCOP consists of a group of autonomous agents, where each agent has an independent computing power. Each agent owns a local constraint network. Variables owned by different agents are connected by constraints. These constraints specify a non-negative constraint cost for combinations of values assigned to the variables they connect. In general, constraints or value assignments may be strategic information or private choice that should not be revealed or delegated to other agents. Thus, each agent only has control on its variables and only knows constraints that involve them. DCOP addresses problems in which agents must, in a distributed manner, assign values to their variables such that the sum of the constraint costs of all constraints is minimized.

Several complete algorithms for solving DCOPs have been proposed in the last decade. The pioneer complete asynchronous algorithm is Adopt [START_REF] Modi | ADOPT: Asynchronous Distributed Constraint Optimization with Quality Guarantees[END_REF]. Later on, the closely related BnB-Adopt [START_REF] Yeoh | BnB-ADOPT: An Asynchronous Branch-and-Bound DCOP Algorithm[END_REF] was presented. BnB-Adopt changes the nature of the search from Adopt best-first search to a depth-first branch-and-bound strategy, obtaining better performance. Gutierrez and Meseguer show that some of the messages exchanged by Adopt and BnB-Adopt turned out to be redundant [START_REF] Gutierrez | Removing redundant messages in n-ary bnb-adopt[END_REF]. By removing these redundant messages they obtain more efficient algorithms: Adopt + and BnB-Adopt + . The algorithms mentioned so far perform assignments concurrently and asynchronously. Thereby, the perception of agents on the variable assignments of other agents is in general inconsistent.

Another category of algorithms for solving DCOPs is that of algorithms performing assignments sequentially and synchronously. The synchronous branch and bound (SyncBB) [START_REF] Hirayama | Distributed partial constraint satisfaction problem[END_REF] is the basic systematic search algorithm in this category. In SyncBB, only the agent holding the token is allowed to perform an assignment while the other agents remain idle. Once it assigns its variables, it passes on the token and then remains idle. Thus, SyncBB does not make any use of concurrent computation. No-Commitment Branch and Bound (NCBB) is another synchronous polynomial-space search algorithm for solving DCOPs [START_REF] Chechetka | No-Commitment Branch and Bound Search for Distributed Constraint Optimization[END_REF]. To capture independent sub-problems, NCBB arranges agents in constraint tree ordering. NCBB incorporates, in a synchronous search, a concurrent computation of lower bounds in non-intersecting areas of the search space based on the constraint tree structure.

Another attempt to incorporate a concurrent computation in a synchronous search was applied in Asynchronous Forward Bounding (AFB) [START_REF] Gershman | Asynchronous Forward-Bounding for Distributed Constraints Optimization[END_REF]. AFB can be seen as an improvement of SyncBB where agents extend a partial assignment as long as the lower bound on its cost does not exceed the global upper bound (i.e., the cost of the best solution found so far). In AFB, the lower bounds are computed concurrently by unassigned agents. Thus, each synchronous extension of the current partial assignment is followed by an asynchronous forward bounding phase. Forward bounding propagates the bounds on the cost of the partial assignment by sending to all unassigned agents copies of the extended partial assignment. When the lower bound of all assignments of an agent exceeds the upper bound, it performs a simple backtrack to the previous assigned agent. Later, the AFB has been enhanced by the addition of a backjumping mechanism, resulting in the AFB BJ algorithm [START_REF] Gershman | Asynchronous Forward Bounding for Distributed COPs[END_REF]. The authors report that AFB BJ, especially combined with the minimal local cost value ordering heuristic performs significantly better than other DCOP algorithms.

In this paper, we propose AFB BJ + , a revisited version of AFB BJ in which we refine the lower bound computations. We also propose to compute lower bounds for the whole domain of the last assigned agent instead of only doing this for its current assignment. Thus, an unassigned agent computes the lower bound for each value in the domain of the agent requesting it. This reduces both the number of messages needed and the time future agents remain idle. Hence, we take all possible advantage from the asynchronicity of the system. In addition, these lower bounds are used as a value ordering heuristic in AFB BJ + . Thus, an agent assigns first values with minimal lower bound.

This paper is structured as follows. Section 2 gives the necessary background on DCOP and a short description of the AFB algorithms. We present the AFB BJ + algorithm in Section 3. Correctness proofs are given in Section 4. We report experimental results in Section 5. Finally, we conclude in Section 6.

Background

Basic definitions and notations

The Distributed Constraint Optimization Problem (DCOP) has been formalized in [START_REF] Modi | ADOPT: Asynchronous Distributed Constraint Optimization with Quality Guarantees[END_REF] as a quadruple (A, X , D, C), where A is a set of p agents {A 1 , . . . , A p }, X is a set of n variables {x 1 , . . . , x n }, where each variable x j is controlled by one agent in A. D = {D 1 , . . . , D n } is a set of n domains, where D j is the set of possible values to which variable x j may be assigned. Only the agent controlling a variable can assign a value to it and has knowledge of its domain. C = {c ij : D i × D j → R + } is a set of binary utility constraints (i.e., soft constraints). Each utility constraint c ij ∈ C is defined over the pair of variables {x i , x j } ⊆ X . We say that x i and x j are neighbors.

For simplicity purposes, we consider a restricted version of DCOP where each agent holds exactly one variable (p = n). Thus, we use the terms agent (A j ) and variable (x j ) interchangeably and we identify the agent ID with its variable index (j). Furthermore, all agents store a unique total order ≺ on agents. Agents appearing before an agent A j ∈ A in the total order are the higher priority agents and those appearing after A j are the lower priority agents. The order ≺ divides the set Γ(x j ) of neighbors of A j into higher priority neighbors Γ -(x j ), and lower priority neighbors Γ + (x j ). For sake of clarity, we assume that the total order is the lexicographic ordering

[A 1 , A 2 , . . . , A n ].
In the rest of the paper, we consider a generic agent A j ∈ A. Thus, j is the level of agent A j .

An assignment for agent A j is a tuple (x j , v j ), where v j is a value from D j . A j maintains a counter t j and increments it whenever it changes its value. The value of the counter tags each generated assignment. When comparing two assignments for the same agent, the most up to date is the one with the greatest tag. A current partial assignment (CPA) is an ordered set of assignments, e.g., Y = [(x 1 , v 1 ), . . . ,(x j , v j )] s.t. x 1 ≺ . . . ≺ x j . The set of all variables assigned in Y is denoted by var(Y ) = {x 1 , . . . , x j }. A timestamp associated to a CPA Y is an ordered list of counters [t 1 , t 2 , . . . , t j ] where t i is the tag of the variable x i s.t. x i ∈ var(Y ) [START_REF] Nguyen | Dynamic Distributed BackJumping[END_REF][START_REF] Wahbi | Nogood-Based Asynchronous Forward-Checking Algorithms[END_REF]. When comparing two CPAs, the strongest one is that associated with the lexicographically greater time-stamp. That is, the CPA with greatest value on the first counter on which they differ, if any, otherwise the longest one. Let Y = [(x 1 , v 1 ), . . . ,(x i , v i ), . . . ,(x j , v j )] be a CPA, the subset of Y including all variables down to x i is denoted by

Y i = [(x 1 , v 1 ), . . . ,(x i , v i )].
The guaranteed cost of a CPA Y , denoted by gc(Y ), is the sum of all utility constraints c ij s.t. x i and x j are assigned in Y (Eq. 1).

gc(Y ) = cij ∈C c ij (v i , v j ) | (x i , v i ),(x j , v j ) ∈ Y . (1)
A full assignment Y is a CPA that involves all variables of the problem, i.e., var(Y ) = X . The goal of a DCOP solver is to distributively find a full assignment Y * with minimal cost, that is,

Y * = arg min Y {gc(Y ) | var(Y ) = X }.
In the following, we will present standard AFB algorithms. We will use a nomenclature of messages and data structures different from those used in the original paper [START_REF] Gershman | Asynchronous Forward Bounding for Distributed COPs[END_REF] in order to be closer to those used in our approach.

lb: Y = [ (x 1 , v 1 ), . . . ,(x j , v j )], lb k (Y ) back: Y = [ (x 1 , v 1 ), . . . ,(x j-1 , v j-1 )]
Fig. 1: The messages exchanged by the AFB algorithm.

Asynchronous Forward Bounding (AFB) algorithm

In Asynchronous Forward Bounding (AFB) [START_REF] Gershman | Asynchronous Forward-Bounding for Distributed Constraints Optimization[END_REF], agents assign their variables sequentially and unassigned agents asynchronously try to compute lower bounds on the CPA, say Y . Agents perform assignments of their variables only when they hold the current partial assignments Y (i.e., Y is the token). Each extension of the CPA Y , is followed by a Forward Bounding (FB) phase. The FB phase is performed by sending forward copies of Y to all unassigned agents. In the FB phase, it is required from unassigned agents to compute a lower bound on the cost increment caused by an assignment of their variables on Y . Once computed, the lower bounds are sent back to the agent that sent the request, i.e., the last assigned agent in Y . Due to the asynchronous nature of the FB phase, multiple CPAs may be present at a given moment in time. However, the time-stamp mechanism is used by agents to discard obsolete ones. The lower bounds collected from unassigned agents are used to compute a lower bound on the CPA. When the computed lower bound becomes larger than the current upper bound (i.e., the cost of the best full CPA found so far), the CPA is pruned. Concretely, whenever agent A j receives a valid lower bound from an unassigned agent, it adds it to that received from other agents and checks if the cumulative lower bound exceeds the upper bound. In such a case, A j tries to assign an alternative value to its variable. If such value is not available, it needs to backtrack. When agent A j takes the decision to backtrack, it sends the CPA (Y j-1 ) backwards to the last agent assigned on it (i.e., A j-1 ). However, if A j is the first agent in the ordering, it ends the search process after claiming this to other agents. The AFB algorithm then reports that the optimal solution is the best full CPA found so far.

Fig. 1 shows the messages exchanged by the AFB algorithm. 4 AFB agents exchange the following types of messages: ok? : a message which contains the CPA Y with its cost gc(Y ). When A j assigns its variable, it sends this message to the next agent in the ordering (A j+1 ). back : a message which contains an inconsistent CPA. It is sent back to agent A j-1 requiring it to change its assignment. fb? : a message which contains a copy of an ok? message. It is sent by A j to unassigned agents to compute a lower bound on the CPA it carries. lb : a message which contains a lower bound on the current partial assignment. It is sent as response to a fb? message.

The computation of lower bounds on AFB is performed as follows. In a preprocessing step, A j computes the minimal future cost estimation for each possible value in its domain incurred by every lower priority agent A k (Eq. 2). f c j (v j ) is a lower bound on the cost of constraints involving the assignment (x j , v j ) and all its lower priority neighbors. Agents compute these estimations only once and store them.

f c j (v j ) = x k ∈Γ + (xj ) min v k ∈D k {c jk (v j , v k )}. (2)
Given a current partial assignment Y i , and an unassigned agent A j , the local cost of assigning a value v j to A j is the sum of the constraints costs of this value with all assignments in Y i s.t. i < j (Eq. 3).

lc j (Y i , v j ) = (x h ,v h )∈ Y i s.t. h≤i<j c hj (v h , v j ) (3)
Summing the local cost of an assignment (x j , v j ) and the minimal future cost incurred by lower priority neighbors provides a future cost of an eventual extension of Y i with(x j , v j ). The lower bound of Y i on a unassigned agent (A j ) is the minimal future cost over all its values (Eq. 4). Thus, whenever a higher priority agent A i requires from A j to compute a lower bound on a CPA Y i , it responds by sending back the minimal lower bound of Y i over all values in its domain, i.e., lb j (Y i ).

lb j (Y i ) = min vj ∈Dj {lc j (Y i , v j ) + f c j (v j )}. (4)
By collecting lower bounds from lower priority agents, agent A j can compute a lower bound on its CPA Y j . The lower bounds on a CPA Y j reported by lower priority agents are accumulated and summed up with the guaranteed cost of Y j to provide a lower bound on the cost of a complete extension of Y j (Eq. 5).

lb(Y

j ) = gc(Y j ) + A k Aj lb k (Y j ) (5)
If the computed lower bound lb(Y j ) exceeds the current known upper bound (U B j ), A j needs to change its current value v j on Y j by a new value v j generating a stronger CPA. Then, search continues with the generated CPA. If A j has already tested all possible values for its variable, it backtracks, asking the previous agent in the ordering (A j-1 ) to assign a new value to its variable through a back message.

Asynchronous Forward Bounding with CBJ (AFB BJ)

The Asynchronous Forward Bounding with backjumping (AFB BJ) was obtained by adding a backjumping mechanism to standard AFB [START_REF] Gershman | Asynchronous Forward Bounding for Distributed COPs[END_REF]. When the lower bounds of all values exceed the upper bound, instead of backtracking to the most recently assigned variable, AFB BJ tries to jump to the last assigned agent such that its re-assignment could possibly lead to a solution. To this end, agents in AFB BJ use some maintained data structures that we introduce in the following. Another feature of AFB BJ is that the agent that performs an assignment, uses the minimal local cost (Eq. 3) as value ordering heuristic.

When an agent A j assigns its variable, it sends an ok? message to the next agent in the ordering. In AFB BJ, the ok? message contains the current partial assignment Y j and an array of guaranteed costs, one for each level. gc(Y

j )[j] is the cost of Y j = [(x 1 , v 1 ), . . . ,(x j , v j )] where gc(Y j )[0] = 0. For each i ∈ 1..j -1, gc(Y j )[i] equals that received from A j-1 , i.e., gc(Y j-1 )[i]. gc(Y j )[j] = gc(Y j ) = gc(Y j-1 ) + lc j (Y j-1 , v j ) (6)
In order to perform backjumping, AFB BJ agents compute a lower bound for each level on the CPA. Concretely, instead of computing the lower bound for the whole received CPA Y i , A j computes it for each Y h where h ≤ i < j. To this end, A j first computes the local cost for each level h (Eq. 7).

lc j (Y i , v j )[h] = lc j (Y h , v j ) (7)
The lower bound at level h (Eq. 8) is then the minimal lower bound of Y h over all values in D j . When a higher agent A i requests from A j to compute its lower bound on a CPA Y i , it answers by sending an array of lower bounds, one for each level h where

h ≤ i < j. lb j (Y i )[h] = min vj ∈Dj {lc j (Y i , v j )[h] + f c j (v j )}. (8)
When A j successfully assigns a value v j to its variable, it sends forward copies of the extended CPA, Y j , to each unassigned agent A k and awaits for receiving from them the array of lower bounds. The lower bounds denoted by lb k (Y j ), is an array in which the i th element (1 ≤ i ≤ j) contains a lower bound on the cost of assigning a value to A k with respect to the assignment on Y i (Eq. 8). Once A j receives lower bounds arrays, it computes a lower bound on the cost of any full assignment (Eqs. 9).

lb(Y j )[i] = gc(Y j )[i] + A k Aj lb k (Y j )[i]. (9)
These lower bounds are used by the AFB BJ to determinate the backjumping level. For more details about the way in which the level of backjumping is calculated we refer the reader to [START_REF] Gershman | Asynchronous Forward Bounding for Distributed COPs[END_REF].

3 Asynchronous Forward Bounding revisited AFB BJ + is a revisited version of AFB BJ in which we propose a refinement of the lower bound computations. We also propose to compute lower bounds for the whole domain of the last assigned agent instead of only doing this for its current assignment. Thus, an unassigned agent computes the lower bound for each value in the domain of the agent requesting it. In addition, these lower bounds are used as a value ordering heuristic.

Lower bound refinement

When an agent A i successfully assigns a value v i to its variable, it sends forward copies of the extended CPA, Y i , to each unassigned agent and awaits for receiving from them the array of lower bounds. When agent A j receives this CPA Y i (through a fb? message), it computes the lower bound for each level h where h ≤ i < j (Eq. 8). When computing the lower bound of level h only assignments on Y h are considered. We also add the cost of assigning v i to x i (i.e.,(x i , v i )) to this lower bound. Moreover, we also add the minimal cost of constraints with variables (x m ) between x h and x i . Thus, instead of being a lower bound on a possible extension of Y h by a possible assignment of A j , it will be a lower bound on a possible extension of Y h by both A i and A j and agents between x h and x i . Hence, we revise Eq. 8 to get Eq. 10 where the first and the last terms remain as in the original equation.

lb j (Y i )[h] = min vj ∈Dj lc j (Y i , v j )[h] + i-1 m=h+1 min vm∈Dm {c mj (v m , v j )} + c ij (v i , v j ) + f c j (v j ) (10)
At first glance, it seems that this will require more computational effort from unassigned agents, however it is not the case. One can simply compute the array of lower bounds, as is already done in AFB BJ, and at the end it adds to each level the cost with variable x i . We obtain the addition of the third term, i.e., c ij (v i , v j ). To get the quantity to be added by the second term (i.e., i-1 m=h+1 min vm∈Dm {c mj (v m , v j )}), we use the same principle used in Eq. 2. Agents compute for each value the estimations of each level only once and store them.

The refinement of the lower bounds computation allows agents to get more accurate lower bounds on their assignments. Thus, the accumulated lower bound at each level is increased. This mechanism allows earlier detection of CPAs with lower bound larger than the upper bound. In addition, by doing this, the back message will be sent as high as possible in the agent ordering, thus saving unnecessary search effort.

Lower bounds for the whole domain

In AFB BJ, the forward bounding phase is very expensive in term of communication load. FB requires for each value in D j , 2 × (nj) messages (a fb? and a lb message for each lower agent). Thus, FB needs, for each CPA Y j-1 , 2×|D j |×(n-j) messages.

In AFB BJ + , we propose to compute lower bounds for the whole domain of the last assigned agent instead of only computing this for its current assignment. When an agent receives a fb? message it answers by sending back a two-dimensional array, an array for each value in the domain of the receiver agent. Hence, the forward bounding phase will need, for each CPA Y j-1 , only 2(nj) messages, 2 messages for each lower agent. When agent A j receives a fb? from agent A i , instead of computing lb j (Y i )[h] only for the current value of x i , A j computes it for each value

v i in D i , lb j (Y i-1 )[h][v i ] using Eq. 11. 5 ∀h ∈ 1..i-1, ∀v i ∈ D i , lb j (Y i-1 )[h][v i ] = lb j (Y i-1 ∪(x i , v i ))[h] (11)

Avoiding redundancy

Another feature of the AFB BJ + algorithm is that agents retain and maintain the received lower bounds to avoid redundant messages. When an agent A j receives a fb? message from agent A i , instead of clearing all information it stores (namely the collected lower bounds and the computed ones with their local costs), it clears only irrelevant information w.r.t the received CPA Y i . Concretely, A j compares the time-stamp of the received CPA with its CPA. If its CPA is stronger than the received one, the message is discarded. Otherwise, A j gets the index h ≤ i of the first counter on which they differ. All local costs and lower bounds on the current partial assignment Y h-1 remain valid. Thus, agent A j will not re-compute lower bounds for this part. The same thing is done for ok? messages. Whenever agent A j receives a CPA Y j-1 , it updates all stored information by only removing parts that are not compatible with Y j-1 . When A j succeeds in assigning its variable, it sends forward copy of the extended CPA Y j in fb? messages to its lower agents. However, some of these messages are redundant. To avoid this, each agent A j stores, for each lower priority agent A k , the agent which is the closest to A j in the neighbors of A k higher than A j . As long as the assignment of such agent or agents higher than him were not updated, there is no need to send fb? message to A k . Thus, redundant messages and computations are saved. Moreover, the agent assigning its variable has more accurate lower bounds for all values in its domain. As long as the new complete array of lower bounds has not yet been received, the remaining valid part can be used as a lower bound estimation for each value in the current domain using Eq. 12. h is the lowest valid level for lower bound received from A k .

lb(Y j ) = gc(Y j ) + k j lb k (Y j-1 )[h][v j ] s.t. (x j , v j ) ∈ Y j (12)

Promising value ordering heuristic

Unlike AFB BJ that uses minimal local cost as value ordering heuristic, an AFB BJ + agent uses a different strategy for reordering values in its current domain. All computations performed so far by unassigned agents to calculate lower bounds are used to reorder values in the current domain. Thus, when receiving an ok? message, A j computes the lower bounds for all values in its domain using Eq. 12. A j chooses to assign first values with minimal lower bound. Then, instead of considering only costs with past variable, both costs with past variables and estimations of costs on future variables are considered. We mimic an informed memory-bounded version of A * , instead of simulating an uninformed memory-bounded version of A * .

AFB BJ + description

Fig. 2 presents the pseudo-code of AFB BJ + executed by every agent A j . Agent A j maintains a variable U B j that stores the current upper-bound (the cost of the best solution found so far) initialized to +∞, v * j that stores the value of A j on the solution, Y that stores the strongest received CPA, GC an array of size j -1 that stores the guaranteed costs where GC from a lower agent A k . Since lb k (Y )[0][v j ] depends only on the assignments of x j and x k , it is initialized to min

[i] = gc(Y i ),
v k ∈D k {c jk (v j , v k )}.
Thus, it is a valid lower bound for all CPAs that contains(x j , v j ). Eq. 2 is obtained by summing lb k (Y )[0][v j ] for each lower agent A k .

AFB BJ + starts by initializing the local data structures of A j (lines 1-4). A j then enters in the waiting and processing message loop (line 6). Each received message holds a CPA msg.Y and its corresponding guaranteed costs msg.GC. Due to the asynchronous nature of the algorithm, some messages may be obsolete. A j uses the time-stamping mechanism to discard those messages (line 9). If the received CPA (msg.Y ) is stronger than Y , A j updates Y and GC by the received ones (line 10). Then, A j clears all irrelevant lower bounds computed or received so far (line 11). Agent A j attaches to each message it sends its U B j . The upper bound U B j and v * j are updated when a received message carries a new upper bound smaller than the stored one (line 8).

Upon receiving an ok? message, A j marks that it must send fb? messages by setting mustSendF B to true. Next, it attempts to extend the received CPA by calling procedure ExtendCPA() (line 13).

When calling ExtendCPA(), A j tries to find a value with the minimum lower bound (Eq. 12) without exceeding U B j (lines [START_REF] Wahbi | DisChoco 2: A Platform for Distributed Constraint Reasoning[END_REF][START_REF] Wahbi | Nogood-Based Asynchronous Forward-Checking Algorithms[END_REF]. If such value does not exist, A j backtracks (Backtrack() call, line 19). Otherwise, A j extends the CPA by adding its new assignment and increments its counter t j . If the resulting CPA includes assignments of all agents (line 22), a solution is found and then the upper bound is updated. Instead of broadcasting the new solution and its associated upper-bound, A j calls ExtendCPA() to continue the search (line 26). Since U B j is always attached to the exchanged messages, other agents will be informed of this new upper bound when continuing the search. At the end of the search, the best assignment of A j is stored in v * j . If A j is not the last agent on the ordering, it sends the extended CPA to the next agent (line 28). Afterwards, A j sends fb? messages to all lower priority agents (lines 30-31). 6When A j receives a fb? message, it computes for each value from the domain of the sender a lower bound on the cost increment caused by adding an assignment to its variable using Eq. 10. These lower bounds are sent back to the agent who sent the fb? message through a lb message.

When A j receives a valid lb message, it saves the attached lower bounds (line 37). It checks if this new information causes the current partial assignment to exceed the upper-bound. In such a case, A j calls ExtendCPA() in order to change its assignment (line 38).

Agent A j calls procedure Backtrack() whenever the lower bounds of all its values exceed the upper-bound. When this occurs, A j computes to which agent the CPA Y should be sent to (the backtracking target). A j goes over all candidates, from j -1 down to 1, looking for the first agent it finds that its reassignment could lead to a full assignment with a cost lower than U B j . This agent is the latest assigned agent A i where lb(Y )[i-1] < U B j (line 33). If such an agent exists, A j sends him a back message (line 34). Otherwise, A j reports this to other agents through stp messages (line 35) and terminates its execution. Lemma 1. AFB BJ + is guaranteed to terminate.

Proof. (Sketch) The proof is close to the one given in [START_REF] Wahbi | Nogood-Based Asynchronous Forward-Checking Algorithms[END_REF]. It can easily be obtained by induction on the agent ordering that there will be a finite number of new generated CPAs (at most d n , where n is the number of variables and d is the maximum domain size), and that agents can never fall into an infinite loop for a given CPA.

To prove that AFB BJ + is correct, we need to prove that the correctness inherent to AFB BJ is not violated by the lower bound refinements and the non-broadcasting of solution messages.

(Sketch) Assuming the correctness of AFB BJ, the lower bounds without refinement terms are consistent. It is thus enough to prove that the costs included in the refinement terms are not redundant. All constraints considered in the calculation of the second and third terms of Eq. 10 have not been included in the first and fourth terms. Moreover, these constraints are not included in the lower bounds computed by other lower priority agents. Therefore, costs added by refinement terms are not redundant and then Eq. 12 is a lower bound on Y j . Lemma 2. By the end of AFB BJ + , each agent stores in its U B j the cost of the optimal solution Y * and in v * j its value on Y * . Proof. (Sketch) In agent A n , lb(Y n ) equals gc(Y n ) because it does not have lower priority agents (Eq. 12). A n updates its U B n and v * n only when it generates a full CPA (Y n ) with gc(Y n ) smaller than its current upper bound (lines 23-24). Thus, U B n only decreases. Let σ be the smallest generated U B n , i.e., σ is the cost of the latest generated full CPA Y n . In AFB BJ + , (i) each agent A j attaches to each message it sends its U B j . Agent A j only updates its U B j and its v * j when the upper bound carried in a received message is smaller than the stored one (line 8, Fig. 2). (ii) All agents will receive at least one message after the generation of σ (at least they will receive stp messages before they stop their execution). (iii) Messages are only sent after receiving and processing other messages. A n attaches σ to each message it sends. Hence, all messages that follow the generation of σ will contain it. Because σ is the smallest generated upper bound and following (i), (ii) and (iii), when the search is ended, U B j of each agent A j equals σ and its v * j equals that assigned to x j in Y n . Now, one needs to prove that σ is the cost of the optimal solution Y * (i.e., σ equals gc(Y * ) and Y n equals Y * ). To prove that σ equals gc(Y * ), it is enough to demonstrate that during search no CPA that can lead to a solution of lower cost than σ is discarded. In AFB BJ + , the CPAs are discarded only in three places (line 19, procedure ExtendCPA(), line 38, procedure ProcessLB(), and line 33, procedure Backtrack()). In all cases above, we are ensured that the lower bound of the discarded CPAs exceeds U B j . Thus, they cannot lead to a solution with a cost smaller than U B j . Now, since σ ≤ U B j when discarding those CPAs, we are ensured that they have a cost larger than σ. Thus, σ is the cost of the optimal solution Y * and then Y n equals Y * . Therefore, when AFB BJ + terminates, v * j is the assignment of x j on Y * . Then, the lemma is proved. This also completes the correctness proof of the AFB BJ + algorithm.

Corollary 1. AFB BJ + is sound, complete, and terminates.

In this section we experimentally compare AFB BJ + to AFB BJ [START_REF] Gershman | Asynchronous Forward Bounding for Distributed COPs[END_REF], BnB-Adopt + [START_REF] Gutierrez | Saving redundant messages in bnb-adopt[END_REF], and BnB-Adopt-DP2 + (BnB-Adopt + combined with DP2 value ordering heuristic [START_REF] Ali | Preprocessing techniques for accelerating the dcop algorithm adopt[END_REF]). Algorithms are evaluated on four commonly used benchmarks: binary random Max-DisCSPs, binary random DCOPs, meeting scheduling and sensor networks. All experiments were performed on the DisChoco 2.0 platform 7 [START_REF] Wahbi | DisChoco 2: A Platform for Distributed Constraint Reasoning[END_REF], in which agents are simulated by Java threads that communicate only through message passing. We evaluate the performance of the algorithms by communication load and computation effort. Communication load is measured by the total number of exchanged messages among agents during algorithm execution (#msg) [START_REF] Lynch | Distributed Algorithms[END_REF]. Computation effort is measured by the number of non-concurrent constraint checks (#ncccs) [START_REF] Zivan | Message delay and DisCSP search algorithms[END_REF]. #ncccs is the metric used in distributed constraint solving to simulate the computation time.

We simulate two scenarios of communication: fast communication (where message delay is null), and slow communication with uniform random message delay, where the delay costs between 0 and 100 #ncccs for each message. On slow communication, the trends are similar to those observed for fast communication, so the results are not reported here.

The benchmark settings

Uniform binary random Max-DisCSPs are characterized by n, d, p 1 , p 2 , where n is the number of agents/variables, d is the number of values per variable, p 1 is the network connectivity defined as the ratio of existing binary constraints, and p 2 is the constraint tightness defined as the ratio of forbidden value pairs (with a cost of 1). We solved instances of two classes of constraint graphs: sparse graphs 10, 10, .4, p 2 and dense graphs 10, 10, .7, p 2 . We varied the tightness from 0.6 to 0.9 by steps of 0.1 and from 0.9 to 0.98 by steps of 0.02. For each pair of fixed density and tightness (p 1 , p 2 ) we report average over 50 instances.

Binary random DCOPs are characterized by n, d, p 1 , where n, d and p 1 are as in Max-DisCSPs [START_REF] Gutierrez | Saving redundant messages in bnb-adopt[END_REF]. For each value combination a cost is selected randomly from the set {0, . . . , 100}. For each p 1 = 0.4, . . . , 0.8, we have generated 50 instances in the class n = 10, d = 10, p 1 .

The meeting scheduling consists of a set of agents, each having a personal private calendar and a set of meetings each taking place in a specified location. The meeting scheduling is encoded as follows. Variables/agents represent meetings. Each meeting/variable has as domain the time slots possible for it. There are constraints between meetings that share participants. We present here 4 cases each with different hierarchical scenarios [START_REF] Yin | USC dcop repository[END_REF].

The sensor network problem consists of a set of sensors that track a set of mobiles. Each mobile must be tracked by 3 sensors. Each sensor can track at most one mobile.

The sensor network problems are encoded as follows. Variables/agents represent mobiles. The possible values of a variable/mobile are all combinations of three sensors that are able to track it. There are constraints between adjacent mobiles. Details are given in [START_REF] Ali | Preprocessing techniques for accelerating the dcop algorithm adopt[END_REF][START_REF] Jung | Argumentation as Distributed Constraint Satisfaction: Applications and Results[END_REF][START_REF] Béjar | Sensor networks and distributed csp: communication, computation and complexity[END_REF]. We present here 4 cases with different topology scenarios [START_REF] Yin | USC dcop repository[END_REF]. 7 http://dischoco.sourceforge.net/ 0 .6 0 0 .7 0 0 .8 0 0 .9 0 0 .9 2 0 .9 4 0 .9 6 0 .9 8 0 .6 0 0 .7 0 0 .8 0 0 .9 0 0 .9 2 0 .9 4 0 .9 6 0 .9 8 0 .6 0 0 .7 0 0 .8 0 0 .9 0 0 .9 2 0 .9 4 0 .9 6 0 .9 8 0 .6 0 0 .7 0 0 .8 0 0 .9 0 0 .9 2 0 .9 4 0 .9 6 0 .9 8

Results & Discussion

The results on instances of the first set of experiments (Max-DisCSPs) are illustrated in Fig. 3. In terms of computational effort (Figs. 3a and3c), AFB BJ + improves the AFB algorithms and performs faster than both BnB-Adopt + algorithms. The factor of improvements is 5 for sparse graphs and 7 for dense graphs. Concerning communication load (Figs. 3b and3d), AFB BJ + requires few messages compared to others algorithms. AFB BJ + improves AFB BJ by a factor of 20 (resp. 15) in sparse (resp. dense) instances. In dense instances, AFB BJ + outperforms BnB-Adopt-DP2 + by a large scale. BnB-Adopt + and BnB-Adopt-DP2 + are the less efficient algorithms for solving Max-DisCSPs, and their performance dramatically deteriorates on dense Max- DisCSP problems. The DP2 heuristic improves the performance of BnB-Adopt + . This improvement is clearer in the sparse problems than in dense ones.

For binary random DCOPs, the results are presented in Table 1. Both versions of BnB-Adopt + dramatically deteriorate compared to algorithms performing assignments sequentially. Again, the DP2 heuristic improves the performance of BnB-Adopt + . AFB BJ + improves the speed-up of AFB BJ by a factor of 6 in dense instances. Regarding the #msg, the factor of improvement is 13.

Table 2 presents the results on meeting scheduling problems. Comparing AFB BJ + to AFB BJ, the obtained results show that AFB BJ + reduces the number of #ncccs by a factor of 10 and the number of required messages by a factor of 50 in all classes. AFB BJ + outperforms BnB-Adopt + by a large factor on both considered measures. However, BnB-Adopt-DP2 + benefits from its preprocessing step and performs faster than AFB BJ + .

For sensor networks, the results are presented in Table 3. Again, AFB BJ + improves the performance of AFB BJ by a large scale. to AFB BJ, AFB BJ + reduces the #ncccs by a factor of 15 and the number of messages by a factor of 50. BnB-Adopt + performs almost the same #ncccs and the same number of messages as AFB BJ + . BnB-Adopt-DP2 + outperforms all other algorithms since it needs very few messages and #ncccs to resolves sensor networks instances.

Looking at all results together, we come to the straightforward conclusion that AFB BJ + performs very well compared to its forward bounding counterparts. The reason for that amounts mainly to refined lower bounds and their use as value ordering heuristic. This guides the search first to promising assignments. The large gap in communication load can be explained by the fact that when an AFB BJ + agent has the token to assign, it sends at-most one request to each lower agent, whereas other AFB algorithms need one message for each lower agent for each possible assignment. In addition, AFB BJ + stores and maintains valid lower bounds to avoid redundant messages.

Both versions of BnB-Adopt + perform very poorly when solving Max-DisCSPs and random DCOPs. One possible reason is that in both algorithms, agents have a strongly asynchronous assignments policy. However, for structured problems, BnB-Adopt-DP2 + has performance close to AFB BJ + . On some highly structured problems (sensor networks), it performs well. When we checked these instances, we found them very sparse with very few constraints. The constraint tree structure used in BnB-Adopt + combined with the very informed DP2 heuristic allows agents, in such very sparse instances, to initialize their lower bounds of values to a cost close to that of the solution.

Our experiments show that AFB BJ + needs less messages than other algorithms. However, AFB BJ + messages can be longer than those sent by other algorithms. The largest messages in AFB BJ + (lb messages) are in O(nd). To see the practical impact of these larger messages, we computed the total number of bytes exchanged by all algorithms. 8 AFB BJ + is improved by BnB-Adopt-DP2 + by a factor up to 2 on meeting scheduling and 47 sensor networks. Except for these two cases, AFB BJ + improves other algorithms in all benchmarks by a factor up to 94 (instead of factor 73 for #msg) for AFB BJ, 144 (instead of 236 for #msg) for BnB-Adopt + , and 80 (instead of 203 for #msg) for BnB-Adopt-DP2 + .

In all our experiments, the longest message sent by AFB BJ + was of size 366 bytes. The minimum datagram size that we are guaranteed to send without fragmentation of a message (in one physical message) is 568 bytes for IPv4 and 1,272 bytes for IPv6 when using either TCP or UDP [START_REF] Bessiere | Agile Asynchronous Backtracking for Distributed Constraint Satisfaction Problems[END_REF]. Thus, counting the number of exchanged messages is equivalent to counting the number of physical messages.

Conclusion

We have proposed AFB BJ + , a revisited version of the AFB BJ algorithm in which we refine the computations of lower bounds by future agents. In AFB BJ + , the lower bounds are computed for the whole domain of the last assigned agent providing him with a very informed value ordering heuristic. Our experiments show that AFB BJ + improves the current state of the art in terms of runtime and number of exchanged messages on different distributed problems. The present work is a step forward in order to address real world applications in multi-agent coordination. Several directions need to be explored in the AFB family. A promising direction is that of variable ordering heuristics. Another direction will be to try to maintain consistencies stronger than forward bounding.

Fig. 2 :

 2 Fig. 2: The AFB BJ + algorithm running on agent Aj.

  + BnB-Adopt-DP2 + (a) #ncccs where p1 = .4

  + BnB-Adopt-DP2 + (b) #msg where p1 = .4

  +BnB-Adopt-DP2 + (c) #ncccs where p1 = .7

7 Fig. 3 :

 73 Fig. 3: Total number of messages sent and #ncccs performed on Max-DisCSP problems in logarithmic scale.

  and lb k (Y )[] that stores the lower bounds received procedure AFB-BJ + () 01. U Bj ← +∞; v * j ← empty; Y ← {}; GC[1..j -1] ← [0, . . . , 0]; 02. mustSendF B ← T rue; 03. foreach ( A k Aj ) do 35. broadcastMsg : stp U Bj ; 36. end ← true; procedure ProcessLB(msg) 37. lb k (Y j ) ← msg.lb ; /* A k is the sender of msg */ 38. if ( lb(Y j ) ≥ U Bj ) then ExtendCPA() ;

	04.	foreach ( vj Dj ) do lb k (Y )[0][vj] ← min v k ∈D k	{c jk (vj, v k )} ;
	05. if ( Aj = A1 ) then ExtendCPA() ;
	06. while ( ¬end ) do
	07.	msg ← getMsg();
	08.	if ( msg.U B < U Bj ) then U Bj ← msg.U B; v * j ← vj ;
	09.	if ( msg.Y is stronger than Y ) then
	10.	Y ← msg.Y ; GC ← msg.GC ;
	11.	clear irrelevant lower bounds ;
	12.	switch ( msg.type ) do
	13.
	23.	U Bj ← gc(Y ) ;	/* Aj = An */
	24.	v * j ← vj ;
	25.	Y ← Y j-1 ;
	26.	ExtendCPA() ;
	27.	else
	28.
	29.	if ( mustSendF B ) then
	30.
	33.	if ( lb(Y )[i-1] < U Bj ) then
	34.

ok? : mustSendF B ← true ; ExtendCPA() ;

14.

back : Y ← Y j-1 ; ExtendCPA() ; 15. fb? : sendMsg : lb lbj(Y i )[], msg.Y to Ai ; /* Ai is msg sender */ 16. lb : ProcessLB(msg); 17. stp : end ← true; procedure ExtendCPA() 18. vj ← arg min v j ∈D j {lb(Y ∪ (xj, v j ))} ; /* Eq. 12 */ 19. if ( lb(Y ∪(xi, vi)) ≥ U Bj ) then Backtrack() ; 20. else 21. Y ← {Y ∪(xj, vj)}; tj ← tj + 1; 22. if ( var(Y ) = X ) then sendMsg : ok? Y, GC, U Bj to Aj+1 ; foreach ( A k Aj ) do sendMsg : fb? Y, GC, U Bj to A k ; 31. mustSendF B ← f alse ; procedure Backtrack() 32. for ( i ← j-1 dowTo 1 ) do sendMsg : back Y i , U Bj to Ai ; return;

Table 1 :

 1 Total number of messages sent and #ncccs performed on binary random DCOPs where costs are randomly selected from 0 to 100.

			#ncccs × 10 3		#msg × 10 3
	p1	0.4 0.5	0.6	0.7	0.8	0.4 0.5		0.6	0.7	0.8
	AFB BJ +	31	77	148	299	554	3	7	14	27	48
	AFB BJ	122 308	654 1,601 3,442	54 111	186	379	658
	BnB-Adopt +	617 3,193 15,436 61,938 98,684 102 419 1,552 5,289 6,312
	BnB-Adopt-DP2 + 180 958 5,266 25,869 75,414	34 151	636 2,549 5,864

Table 2 :

 2 Total number of messages sent and #ncccs performed on Meeting Scheduling.

			#ncccs			#msg
	cases	A	B	C	D	A	B	C	D
	AFB BJ +	4,987	6,536	2,789	2,206	373	871	536	582
	AFB BJ	30,332 101,206 15,841 32,364 7,944 32,262 9,441 17,443
	BnB-Adopt +	272,490 63,352 51,134 30,030 15,507 10,472 8,717 8,278
	BnB-Adopt-DP2 +	5,371	4,224	2,165	1,647	636	749	511	485

Table 3 :

 3 Total number of messages sent and #ncccs performed on Sensor Network.

			#ncccs			#msg
	cases	A	B	C	D	A	B	C	D
	AFB BJ +	5,599	6,182	2,395	4,869 2,043 1,999	325 1,430
	AFB BJ	167,862 190,423 12,084 33,988 127,544 145,421 7,853 33,280
	BnB-Adopt +	4,052	6,337	6,561	8,982	876 1,215 1,198 2,072
	BnB-Adopt-DP2 +	992	1,046	982	1,278	195	238	176	323

The names of the message types here are closer to that used in our approach and then different from that used in the original AFB paper[START_REF] Gershman | Asynchronous Forward-Bounding for Distributed Constraints Optimization[END_REF].

If xi and xj are not neighbors, a simple array is sufficient since the lower bound is the same for all values in Di. Moreover, xj does not known Di.

In our implementation the fb? messages are sent under certain conditions to avoid redundancy, see Section 3.3.

In our implementation we do not perform any message compression.


L Felipe Benites

Nicole Poulton

Karine Labadie

Mike Sieracki

Nigel Grimsley

Gwenael Piganeau

Single Cell Ecogenomics reveals mating types of individual cells and ssDNA viral infections in the Smallest Photosynthetic Eukaryotes

Keywords: Picoeukaryotes, single amplified genome, mating type, cross-contamination, virus, Mamiellophyceae, Tara-Oceans

Planktonic photosynthetic organisms of the class Mamiellophyceae include the smallest eukaryotes (<2 µm), are globally distributed and form the basis of coastal marine ecosystems. Eight complete fully annotated 13 to 22 Mb genomes from three genera, Ostreococcus, Bathycoccus and Micromonas, are available from previously isolated clonal cultured strains and provide an ideal resource to explore the scope and challenges of analysing Single Cell Amplified Genomes (SAGs) isolated from a natural environment. We assembled data from 12 SAGs sampled during the Tara Oceans expedition to gain biological insights about their in-situ ecology, which might be lost by isolation and strain culture. Although the assembled nuclear genomes were incomplete, they were large enough to infer the mating types of four Ostreococcus SAGs. The systematic occurrence of sequences from the mitochondria and chloroplast, representing less than 3% of the total cell's DNA, intimates that SAGs provide suitable substrates for detection of non-target sequences, such as those of virions. Analysis of the non-Mamiellophyceae assemblies, following filtering out crosscontaminations during the sequencing process, revealed two novel 1.6 and 1.8 kb circular DNA viruses, and the presence of specific Bacterial and Oomycete sequences suggests that these organisms might co-occur with the Mamiellales.

Introduction

Planktonic photosynthetic eukaryotes of the class Mamiellophyceae (Chlorophyta) include the smallest eukaryotes, with bacterial-sized cells of less than 2 µm diameter.

Environmental surveys based on the sequence of the highly conserved 18S ribosomal gene (18S) eukaryotic barcode [START_REF] Piganeau | How and why DNA barcodes underestimate the diversity of microbial eukaryotes[END_REF] revealed their worldwide distribution (2)(3)(4) [START_REF] Tragin | Green microalgae in marine coastal waters: The Ocean Sampling Day (OSD) dataset[END_REF]. Their ubiquity, high abundance and turnover suggest they sustain the marine ecosystems in many coastal areas [START_REF] Not | A single species, Micromonas pusilla (Prasinophyceae), dominates the eukaryotic picoplankton in the Western English Channel[END_REF] [START_REF] Worden | Assessing the dynamics and ecology of marine picophytoplankton: The importance of the eukaryotic component[END_REF]. Many flagellate species in this class were described last century before the advent of molecular biology techniques and their first application to these picoalgae [START_REF] Simon | Fluorescent in situ hybridization with rRNA-targeted oligonucleotide probes to identify small phytoplankton by flow cytometry[END_REF]; Monomastix 1912 [START_REF] Scherffel | Zwei neue, trichocystenartige Bildungen führende Flagellaten[END_REF], Micromonas pusilla 1951 [START_REF] Knightjones | Chromulina-Pusilla Butcher, a Dominant Member of the Ultraplankton[END_REF], and Mamiella gilva 1964 [START_REF]Moestrup øjvind. Further studies on Nephroselmis and its allies (Prasinophyceae). II. Mamiella gen. nov., Mamiellaceae fam. nov., Mamiellales ord. nov[END_REF], followed by Crustomastix and Dolichomastix 2005 [START_REF] Manton | Dolichomastix (Prasinophyceae) from arctic Canada, Alaska and South Africa: a new genus of flagellates with scaly flagella[END_REF]. A non-flagellate member was described 1990 as Bathycoccus prasinos [START_REF] Eikrem | The Ultrastructure of Bathycoccus Gen-Nov and Bathycoccus-Prasinos Sp-Nov, a Nonmotile Picoplanktonic Alga (Chlorophyta, Prasinophyceae) from the Mediterranean and Atlantic[END_REF] and the non-scaled coccoid Ostreococcus tauri was described in 1995 [START_REF] Courties | A new marine picoeukaryote: Ostreococcus tauri gen. et sp. nov. (Chlorophyta, Prasinophyceae)[END_REF]. A thorough phylogenetic analysis of the nuclear and plastid encoded rDNA operon led Marin and Melkonian to define the class Mamiellophyceae, divided into three orders Mamiellales, Dolichomastigales and Monomastigales [START_REF] Marin | Molecular phylogeny and classification of the Mamiellophyceae class. nov. (Chlorophyta) based on sequence comparisons of the nuclear-and plastid-encoded rRNA operons[END_REF]. While twenty-two species are presently described in AlgaeBase [START_REF] Guiry | AlgaeBase. World-wide electronic publication[END_REF], recent environmental sequencing suggests that many more species have not been yet isolated [START_REF] Tragin | Green microalgae in marine coastal waters: The Ocean Sampling Day (OSD) dataset[END_REF] [START_REF] Tragin | Novel diversity within marine Mamiellophyceae (Chlorophyta) unveiled by metabarcoding[END_REF] so that the number of species within this class remains an open question. The relative ease of culture of strains of the picoalgae Bathycoccus, Micromonas, and Ostreococcus in Keller's and related media [START_REF] Keller | MEDIA FOR THE CULTURE OF OCEANIC ULTRAPHYTOPLANKTON1,2[END_REF] fostered their development as new models for cell biology [START_REF] Ral | Starch division and partitioning. A mechanism for granule propagation and maintenance in the picophytoplanktonic green alga Ostreococcus tauri(1[w])[END_REF] [START_REF] O'neill | Circadian rhythms persist without transcription in a eukaryote[END_REF] and environmental studies [START_REF] Sanchez-Ferandin | A new, sensitive marine microalgal recombinant biosensor using luminescence monitoring for toxicity testing of antifouling biocides[END_REF] [START_REF] Demory | Picoeukaryotes of the Micromonas genus: sentinels of a warming ocean[END_REF]. Complete and annotated nuclear genomes have been obtained for eight species to date: O. tauri RCC4221 [START_REF] Blanc-Mathieu | An improved genome of the model marine alga Ostreococcus tauri unfolds by assessing Illumina de novo assemblies[END_REF], O. lucimarinus [START_REF] Palenik | The tiny eukaryote Ostreococcus provides genomic insights into the paradox of plankton speciation[END_REF], O. sp. RCC809, O. mediterraneus [START_REF] Yau | Virus-host coexistence in phytoplankton through the genomic lens[END_REF], Micromonas pusilla and M. commoda [START_REF] Van Baren | Evidencebased green algal genomics reveals marine diversity and ancestral characteristics of land plants[END_REF], B. prasinos RCC1105 [START_REF] Moreau | Gene functionalities and genome structure in Bathycoccus prasinos reflect cellular specializations at the base of the green lineage[END_REF], and O. tauri RCC1115 [START_REF] Blanc-Mathieu | Population genomics of picophytoplankton unveils novel chromosome hypervariability[END_REF].

Despite the progress enabled by lab experiments and sequence analyses, the sexual life cycles of these haploid picoalgae, as well as their interactions with bacteria [START_REF] Abby | Bacteria in Ostreococcus tauri cultures -friends, foes or hitchhikers?[END_REF] and viruses [START_REF] Yau | Molecular ecology of Mamiellales and their viruses in the marine environment[END_REF] in situ remain enigmatic. Single amplified genomes (SAGs) assemblies, produced by multiple displacement amplification (MDA) after extraction from sorted single cells have been previously reported to contain foreign DNA [START_REF] Yoon | Single-cell genomics reveals organismal interactions in uncultivated marine protists[END_REF] [START_REF] Bhattacharya | Identification of a Marine Cyanophage in a Protist Single-cell Metagenome Assembly[END_REF]. They may be viewed conceptually as metagenomes, provided they preserve a molecular signature of the ecological context of the cell. SAGs of cells directly sampled from the environment may not only foster knowledge on diversity in microorganisms [START_REF] Vannier | Survey of the green picoalga Bathycoccus genomes in the global ocean[END_REF] [START_REF] Strassert | Single cell genomics of uncultured marine alveolates shows paraphyly of basal dinoflagellates[END_REF], but also open opportunities for the identification of all kinds of novel biological associations, from predation to parasitism [START_REF] Yoon | Single-cell genomics reveals organismal interactions in uncultivated marine protists[END_REF] [START_REF] Galindo | Combined cultivation and single-cell approaches to the phylogenomics of nucleariid amoebae, close relatives of Fungi[END_REF].

Here, we investigated 12 Mamiellales' SAGs sampled in the Indian Ocean during the Tara Oceans expedition [START_REF] Alberti | Viral to metazoan marine plankton nucleotide sequences from the Tara Oceans expedition[END_REF] to explore the biological insights this data provides within an ecogenomic framework. Therefore, we analysed the taxonomic affiliations of the assembled sequence data to discuss the power and challenges of using SAGs to retrieve in-situ ecological information [START_REF] Piganeau | How and why DNA barcodes underestimate the diversity of microbial eukaryotes[END_REF] from the target sequence data, e.g. Mamiellales sequences, about the mating type of the cells and (2) from the non-target sequences identified.

Material and Methods

Data

Single-cells were collected during the Tara Oceans expedition [START_REF] Karsenti | A holistic approach to marine eco-systems biology[END_REF] from surface waters in the Indian Ocean at station 39, 41 and 46 [START_REF] Oceans Consortium | Registry of all stations from the Tara Oceans Expedition[END_REF] and cryopreserved using glycine betaine. Flow cytometric cell sorting, single cell lysis and whole genome amplification by Multiple Displacement Amplification (MDA) [START_REF] Dean | Comprehensive human genome amplification using multiple displacement amplification[END_REF] were performed at the Bigelow Single Cell Genomics facility (Boothbay, Maine USA), as previously described [START_REF] Stepanauskas | Matching phylogeny and metabolism in the uncultured marine bacteria, one cell at a time[END_REF] [START_REF] Martinez-Garcia | Unveiling in situ interactions between marine protists and bacteria through single cell sequencing[END_REF] [START_REF] Stepanauskas | Matching phylogeny and metabolism in the uncultured marine bacteria, one cell at a time[END_REF][START_REF] Martinez-Garcia | Unveiling in situ interactions between marine protists and bacteria through single cell sequencing[END_REF][START_REF] Mangot | Accessing to the genomic information of unculturable oceanic picoeukaryotes by combining multiple single cells[END_REF] [START_REF] Martinez-Garcia | Unveiling in situ interactions between marine protists and bacteria through single cell sequencing[END_REF][START_REF] Stepanauskas | Matching phylogeny and metabolism in the uncultured marine bacteria, one cell at a time[END_REF][START_REF] Martinez-Garcia | Unveiling in situ interactions between marine protists and bacteria through single cell sequencing[END_REF] [START_REF] Stepanauskas | Matching phylogeny and metabolism in the uncultured marine bacteria, one cell at a time[END_REF][START_REF] Martinez-Garcia | Unveiling in situ interactions between marine protists and bacteria through single cell sequencing[END_REF][START_REF] Stepanauskas | Matching phylogeny and metabolism in the uncultured marine bacteria, one cell at a time[END_REF]. The resulting SAGs were screened using universal eukaryotic 18S rDNA PCR primers [START_REF] Mangot | Accessing the genomic information of unculturable oceanic picoeukaryotes by combining multiple single cells[END_REF]. The 12 SAGs we analysed here were selected for sequencing based on the 18S identities and were sequenced using Illumina HiSeq technology at the Genoscope [START_REF] Vannier | Survey of the green picoalga Bathycoccus genomes in the global ocean[END_REF]. All SAGs were multiplexed, SAG1 to 6 and SAG7 to 12 were sequenced on two different runs.

The raw Illumina data was processed as previously described [START_REF] Alberti | Viral to metazoan marine plankton nucleotide sequences from the Tara Oceans expedition[END_REF]: adapters and primers used during library construction were removed from the whole reads and both ends were trimmed for low quality nucleotides (quality value < 20). The longest sequence without adapters and with fewest low quality bases was kept. Sequences between the second unknown nucleotide (N) and the end of the read were also trimmed. Reads shorter than 30 nucleotides after trimming were discarded. These trimming steps were achieved using fastx_clean, an internal software based on the FASTX library. The reads and their mates that mapped onto run quality control sequences (Enterobacteria phage PhiX174 genome, NC_001422.1) were removed using SOAP aligner. The total number of paired-end reads from the processed raw reads ranged between 19 and 36 millions with an average of 26 million reads per SAG.

SAG assembly and taxonomic affiliation of contigs

Individual pair-end read files were assembled using SPAdes v3.9.0 [START_REF] Bankevich | SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing[END_REF] and the final assembly statistics were generated with Quast v4.6.3 [START_REF] Gurevich | QUAST: quality assessment tool for genome assemblies[END_REF]. Contigs smaller than 400 bp were discarded from downstream analysis. Genome completeness was assessed with BUSCO v3.0.2 using the 303 eukaryotic single copy orthologs dataset [START_REF] Simao | BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs[END_REF] and compared to the values obtained on the available closely related sequenced genomes (Table 1).

The SAG assemblies were screened for contigs containing the 18S using a blastn homology search using a custom database of 18S extracted from GenBank complemented with 18S genes of Mamiellales extracted from the genome sequence (supplemental Mamiellophyceae. We considered following high order taxonomic ranks: Archaea, Bacteria, Fungi, Mamiellophyceae, Metazoa, Oomycetes, Other (when several taxa were mixed in uneven proportions within the 10 best hits), Other protists, Plants (Streptophyta),

Unassigned contigs and Virus. To estimate the number of reads recruited to each taxonomic rank in each SAG assembly, we aligned the reads against each assembly with BWA (52) and extracted the read coverage of each contig estimated from the SAMtools mpileup and depth suite [START_REF] Li | The Sequence Alignment/Map format and SAMtools[END_REF].

The sum of contig lengths affiliated to a given group was plotted with ggplots in the R environment [START_REF] Core | R: A language and environment for statistical computing[END_REF].

Cross-contamination filtering

To remove contigs derived from cross-contamination during the sequencing step, all SAG assemblies from the same Illumina run were compared against all with blastn (e-value cutoff 1e-5). This led to the inclusion of one SAG from a MAST lineage (TOSGAG23-1, GenBank accession ERX1271190) that was multiplexed with SAGs 7-12. Each contig was flagged as candidate contaminant if its alignment length with a contig from another SAG was equal or higher than 95% of the length of its best hit, with a nucleotide identity higher than 99.5%. The number of reads recruited to each candidate contig was used to define the source of the contamination: the SAG containing the contig with the highest read coverage was defined as the source SAG, and the corresponding contigs from the other SAGs were considered as contaminants and discarded.

Annotation of non-target contigs

Because the virus contigs detected in SAGs 8 and 12 had best hits with circular viruses, the contigs were first circularized with the Geneious program v10.0.5 [START_REF] Kearse | Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data[END_REF], and open reading frames (ORF) were predicted with Glimmer [START_REF] Salzberg | Microbial gene identification using interpolated Markov models[END_REF].

The conserved motifs of the viral replicase gene (57) [START_REF] Krupovic | Networks of evolutionary interactions underlying the polyphyletic origin of ssDNA viruses[END_REF] were manually searched by inspection of the predicted amino acids sequence located in the N-terminal endonuclease domain and in the C-terminal superfamily 3 helicase domain using MAFFT v7.305b (-localpair -maxiterate 1000). The best blastp hits to these ORFs were complemented with sequences retrieved from GenBank and the Ocean Atlas genes database (59) (supplemental table 2). Multiple sequence alignments were performed with MAFFT v7.305b (--localpair -maxiterate 1000) (60), and trimmed with trimAI v1.2 (-automated1 method) [START_REF] Capella-Gutierrez | trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses[END_REF]. For each alignment, protein evolution models were selected using ProtTest v3.4.2 [START_REF] Darriba | jModelTest 2: more models, new heuristics and parallel computing[END_REF] with the Akaike Information Criterion (AIC) (LG+I+G for ORF1 and RTREV+I+G+F for ORF2).

Phylogenies were built by the maximum likelihood method implemented in RAxML v. 8.2.9 [START_REF] Stamatakis | RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies[END_REF], and plotted in FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).

Screening assemblies for information about mating type

In Ostreococcus species, population genomics analyses have revealed two anciently diverged 500 kb genomic regions with supressed recombination, which are candidate mating type regions [START_REF] Blanc-Mathieu | Population genomics of picophytoplankton unveils novel chromosome hypervariability[END_REF]. As in many bipolar mating type systems, orthologous genes on this mating type locus (MT) cluster phylogenetically by mating type, not by species [START_REF] Blanc-Mathieu | Population genomics of picophytoplankton unveils novel chromosome hypervariability[END_REF]. This is as expected if the origin of mating types predates speciation within this genus, as previously observed in the Volvox genus [START_REF] Ferris | Rapid evolution of sex-related genes in Chlamydomonas[END_REF]. We took advantage of this deep divergence between 23 orthologous genes families of the two mating types, gametologs [START_REF] Garcia-Moreno | Rooting a phylogeny with homologous genes on opposite sex chromosomes (gametologs): a case study using avian CHD[END_REF], to screen the SAG assemblies for homologous sequences. Briefly, orthologous gene families were obtained following an 'all-against-all' protein sequence similarity search, performed with blastp (maximum e-value 1e-4) and gene families were delineated using OrthoFinder (version

Biogeographic occurrence of Mamiellales SAG-associated sequences

The non-targeted bacterial sequences (supplemental table 3) were searched against 957 manually curated partial assemblies (metagenome-assembled genomes MAG) reconstructed from 93 TARA Oceans metagenomes [START_REF] Delmont | Nitrogen-fixing populations of Planctomycetes and Proteobacteria are abundant in surface ocean metagenomes[END_REF], including those stations from the North Indian Ocean the Mamiellophyceae where the SAGs were isolated. A match against a MAG was defined as an alignment of 90% of the length of SAG-associated bacterial contig with a nucleotide identity higher than 98%. Nucleotide identity and total alignment length for each SAG associated bacterial contig were computed from blastn alignments [START_REF] Altschul | Basic local alignment search tool[END_REF].

The geographic occurrences and distribution of the non-target virus sequences was obtained by searching for homologous sequences to the two predicted protein-coding genes, a capsid and a replicase, in the OCEAN GENE ATLAS (59) using blastp (e-value 1e-5). Amino-acid sequences for each hit were downloaded from the database (supplemental table 2). This atlas contains a gene catalogue from 243 metagenomes and 441 metatranscriptomes from the Tara Oceans expedition, as well as the metagenomes from the Global Ocean Sequencing (GOS) expedition (68).

Results

Completeness of Target Assemblies

Phylogenetic analysis of 18S data retrieved from the de novo assembly was consistent with prior taxonomic affiliations. Ostreococcus SAG 18S sequences retrieved from the de novo assembly clustered with O. spp RCC809. However, while SAG9, 10, 11 sequences were 100% identical to RCC809 (and RCC143), the 18S sequence of SAG11 had one difference to RCC809. The 18S of the Bathycoccus SAGs were 100% identical to the 18S sequence from the Bathycoccus prasinos RCC1105 reference genome and the 18S found in Micromonas SAGs were 100% identical to the recently described M. bravo species (69) (Figure 1). We could not retrieve the 18S rDNA sequence from the assemblies in SAG4 and SAG6-7.

The SAG target assemblies recruit the overwhelming majority of reads, from a minimum of 82% of reads in SAG8 to 99.96% in SAG9 (Table 2). Target assemblies are fragmented and summed up to a maximum of 4.2 Mb (SAG2), corresponding to a maximum of one fourth of the expected genome size (SAG2, Bathycoccus, SAG10, Ostreococcus). This modest genome representation is also reflected by the BUSCO analysis, with a maximum of 28.7 % of represented genes (SAG9), to 1% in SAG6, where the target sequences are overwhelmingly represented by the chloroplast and mitochondrial sequences (Table 2). As a comparison, the same BUSCO analyses of the available O. tauri, B. prasinos and M. commoda genomes gave 86%, 84 and 89% respectively, suggesting that BUSCO analysis gives a conservative estimation of genome completeness in the Mamiellales. The relatively low target nuclear genome recoveries are in line with previous studies, which showed that the combination of sequences obtained from several SAGs may increase target genome recovery from 20 to 70% [START_REF] Mangot | Accessing the genomic information of unculturable oceanic picoeukaryotes by combining multiple single cells[END_REF].

There was a large variation in the proportion of reads recruited to the mitochondrial (0.05 to 26%) and chloroplast (0.01 to 73%) genomes. Importantly, the assembly lengths were always large enough (> 20 kb, Table 1) to annotate several complete protein-encoding genes. This sequence data is sufficient for the unambiguous identification of contigs from both the mitochondrial and plastid genomes. However, because the percentage of recruited reads also depends on the copy number of the chloroplast and mitochondria genomes within a cell, we first estimated the nuclear to plastid or to mitochondria ratios in Ostreococcus from available data. In O. tauri, this relative copy number can be inferred from the ratio of the nuclear (28) versus organellar genome coverage (71) obtained from the same Illumina run of a haploid culture. Mitochondrial genome copy number was estimated to range between 2 and 5, chloroplast copy number ranged between 3 and 7, using the available Illumina datasets obtained from 13 independent O. tauri cultures (supplemental table 4). The organellar genomes represent less than 0.5% of the nuclear genome sizes in Mamiellophyceae (Table 1), but may reach 3.5% with the maximum organellar genome copy number.

If the proportion of reads reflects the proportion of DNA in the sampled cell, the percent of reads recruited to the organellar genomes may thus be expected to range from 0.01 to 3.5%. This is the case in all Bathycoccus SAGs and two Ostreococcus SAGs (Table 2).

The complete mitochondrial (mtDNA) or chloroplastic (cpDNA) genomes could not be assembled from individual SAG data (Geneious inbuilt assembler with minimum contig overlap 100 bp and a minimum overlap Identity 90%). However, running the assembler on all contigs affiliated to cpDNA or mtDNA within each genera led to large (>30kb) assemblies of 3 mtDNA and 2 cpDNA. Manual assembly guided by the synteny with the available reference genomes (Table 1) led to 2 complete assemblies, and 3 partial assemblies, summarized in Table 3. To compare the mitochondrial genome assembly obtained for Bathycoccus SAGs with the available reference mitochondrial genome of Bathycoccus prasinos RCC1105 ( 27), the amino-acid identity was calculated along a concatenation of 6 orthologous genes outside the inverted repeat region (nad5, nad4, coc3, cox2, atp6 and nad1) as previously described [START_REF] Piganeau | How and why DNA barcodes underestimate the diversity of microbial eukaryotes[END_REF]. Intraspecific diversity within Mamiellophyceae

organelles is yet unknown for most species, but in O. tauri intraspecific amino-acid divergence is below 0.1% [START_REF] Blanc-Mathieu | Organellar inheritance in the green lineage: insights from Ostreococcus tauri[END_REF]. The amino-acid identity over the 1789 aa long alignment was 81%, supporting the previously published conclusion based on the analysis of the combined nuclear genome assembly from SAG1-4, that the Bathycoccus sequenced belong to a novel, yet uncultured, Bathycoccus species [START_REF] Vannier | Survey of the green picoalga Bathycoccus genomes in the global ocean[END_REF].

Mating type inference in Ostreococcus SAGs

Contigs from the target fraction in the SAGs assemblies were searched for genes located within the candidate MT-and MT+ gene families defined in the available reference genomes. Since the sequences of the two mating types have not yet been identified in Bathycoccus and Micromonas, the inference of mating types is not yet possible in these genera. The amino-acid sequence of orthologous genes could be retrieved for 6 out of the 23 gene families (Figure 2 A). This low gene number is not surprising given the low nuclear genome coverage of each SAG (Table 1), but gene family 16 (GF16 -encoding for a diacylglycerol glucosyltransferase) had matches with 3 out of the 4 Ostreococcus SAGs.

Comparison of the amino-acid sequence of these genes with genes from Ostreococcus sp.

RCC809 indicate that the four Ostreococcus SAGs have the same mating type as RCC809, that is MT+ (Figure 2B).

Taxonomic Diversity of Non-target contigs

The non-target assemblies recruited 0.01 to 6% of the reads (Table 1). Taxonomic affiliation of total non-target contigs is provided in Figure 3A and varies greatly between SAGs. For eight SAGs, most non-target reads were recruited to Bacteria, for 2 SAGs, most non-target reads were recruited to Oomycetes, and for 2 SAGs, most reads were recruited to sequences of viral origin. Sequence comparison between non-target assemblies revealed 100% identical sequences shared between SAGs coming from different stations, but sequenced on the same Illumina lane, a suspicious signature suggesting cross-contamination during the amplification and sequencing step [START_REF] Lejzerowicz | Accurate multiplexing and filtering for highthroughput amplicon-sequencing[END_REF]. This led us to define a cross-contamination filter procedure for all SAGs sequenced within the same run. To find out which of the SAGs was the most likely source of the contamination, we used the number of reads recruited to nontarget contigs. The SAG containing the most sequenced reads was considered as the source SAG of the non-target contig. The taxonomic affiliation of non-target contigs following filtering is presented in Figure 3B. Basically, this led to a 2 to 7-fold reduction of the total non-target assembly length depending on the SAG (Table 2). Following cross-contamination filtering, the maximum percent of reads recruited to non-target sequences dropped from 6% to 1% in SAG8, and became almost negligible (0.002%) in SAG10. The taxonomic affiliation of these filtered non-target sequences could be divided into three groups: Bacteria, Oomycetes/Fungi and Viruses, in order of decreasing frequency.

Most non-target reads were recruited to Bacteria for 11 SAGs and to virus for one SAG (SAG12) (Figure 3B). The non-target bacterial contigs had a large phylogenetic spread (Figure 4A) and while the most common phylum belonged to the Proteobacteria, there was no evidence of any specific bacteria -SAG association from this analysis. To gain more information about the origin of the 83 non-target bacterial sequences, we estimated their similarity with 957 manually curated partial assemblies (metagenome-assembled genomes MAGs) reconstructed from 93 TARA Oceans metagenomes [START_REF] Delmont | Nitrogen-fixing populations of Planctomycetes and Proteobacteria are abundant in surface ocean metagenomes[END_REF]. Over one third (32 out of 83) bacterial contigs had more than 98% identity over 90% of their length with a TARA MAG, suggesting indeed that these bacteria belong to the most abundant bacteria sequenced within the TARA Ocean sampling (Figure 4A). However, only 3 bacterial contigs had a hit against a MAG from the North Indian Ocean, while most hits were against MAGs from the Mediterranean Sea and South East Pacific Ocean. Therefore, except for the 3 bacterial contigs assembled from SAG3 and SAG4, the bacteria sampled together with the SAG do not represent the most dominant bacteria at the sampling site.

The cross-contamination filtering removed most of Oomycetes affiliated non-target contigs (Figure 3B), and the source was one uncultured marine Stramenopile (MAST) SAG [START_REF] Mangot | Accessing the genomic information of unculturable oceanic picoeukaryotes by combining multiple single cells[END_REF] sequenced in the same lane as SAG7 to SAG12. The class Oomycetes is also from the Stramenopiles group and form a distant, yet phylogenetically related clade to MASTs, which currently lack representation in reference databases [START_REF] Seeleuthner | Single-cell genomics of multiple uncultured stramenopiles reveals underestimated functional diversity across oceans[END_REF]. Following the crosscontamination filtration step, non-target Oomycete contigs were affiliated to four families (Figure 4B): Pythiaceae (2 contigs), Albuginaceae (3 contigs) and the more abundant Peronosporaceae (30 contigs) and Saprolegniaceae (33 contigs). Given the available genome size estimations in Oomycetes [START_REF] Mcgowan | Comparative Analysis of Oomycete Genome Evolution Using the Oomycete Gene Order Browser (OGOB)[END_REF] or MASTs [START_REF] Mangot | Accessing the genomic information of unculturable oceanic picoeukaryotes by combining multiple single cells[END_REF], which are larger than 20 Mbp, the size of the largest target Mamiellales, Micromonas (Table 1), the quantity of Oomycetes DNA present within the SAGs precludes any evidence of co-sampling, but rather suggests random capture of DNA during the cell sorting process.

The least frequent non-target sequence group was affiliated to viruses. Analysis of the contigs affiliated to viruses in SAG8 (Micromonas) and SAG12 (Ostreococcus) led us to assemble two novel circular viral genomes, which we named "Mamiellales SAG associated circular virus" (MACV). Genome sizes varies between 1586 bp (SAG12 -Ostreococcus Figure 5A) and 1788 bp (in SAG8 -Micromonas), both coding for two ORFs, a replicase gene (ORF1) and putative capsid gene (ORF2). The two genomes are 100% identical, except for a 202 bp deletion in the shorter version, leading to smaller ORF1 sequence. The experimental analysis performed would not discriminate whether these sequences arose from single or double stranded DNA. The read coverage of the SAG8 and SAG12 viruses is 50x and 7893x respectively. Assuming MDA amplification did not induce coverage biases between the virus and nuclear genomes, the same average coverage would be expected for viral and nuclear nucleotides. Compared to the average coverage of the nuclear genome, the SAG8 associated virus had a 10-fold lower coverage, while the SAG12 associated virus had a 5fold higher coverage. This supports the notion that the virus was present in multiple copies (5 copies) within the cytoplasm of SAG12, suggesting an infected state, while there were fewer copies, and may have been under-amplified, in SAG8. The majority of circular DNA viruses were "Circular Replication-associated protein Encoding Single-Stranded DNA"

(CRESS) viruses [START_REF] Zhao | Eukaryotic Circular Rep-Encoding Single-Stranded DNA (CRESS DNA) Viruses: Ubiquitous Viruses With Small Genomes and a Diverse Host Range[END_REF]. They have been found within diverse hosts and marine metagenomes, even though there is not much information regarding their infection strategies. To explore the similarity of the Mamiellales associated circular virus to other CRESS DNA viruses, we performed Maximum likelihood phylogenetic analyses of their predicted capsid (Figure 5C) and predicted replicase (Figure 5D) proteins, together with their retrieved blast hits (e-value cut off 1E-5) both in GenBank and in the TARA-Oceans Gene Atlas. Capsid homologous sequences were found in seven TARA-Oceans metagenomes (Figure 5B) including station 41 in the Indian Ocean, the same station where the Ostreococcus SAG12 was sampled. The phylogenies of both genes suggested that Mamiellales associated circular viruses formed a well-supported monophyletic clade, which included circo-like viruses, metagenomic sequences, and plant ssDNA viruses (Yerba mate virus and Banana bunchy top virus).

Sequences from other ssDNA viruses clades, such as genes encoded by Bacilladnaviruses (from Diatoms), belonged to divergent branches suggesting that Mamiellales associated circular viruses might form a divergent group of a novel ssDNA virus family.

Discussion

SAGs have so far been previously used to increase our knowledge of the uncultivable multitude of marine microorganisms (40) [START_REF] Seeleuthner | Single-cell genomics of multiple uncultured stramenopiles reveals underestimated functional diversity across oceans[END_REF]. The use of SAGs to discover novel species does not necessitate high genome coverage as only a handful of highly conserved genes is informative to delineate a species. However, partial genome coverage may constitute a greater challenge for the use of SAGs in evolutionary studies. To increase the genome coverage of a target species, data from genetically close SAGs may be merged together (33)(45). This approach has been successfully applied by Vannier et al [START_REF] Vannier | Survey of the green picoalga Bathycoccus genomes in the global ocean[END_REF] to the 4

Bathycoccus SAGs and resulted in an assembly estimated to contain 64% of the complete genome. The amino-acid identity between this novel assembly and the reference B. prasinos reported, but the maximum number of strains isolated from the same sampling point is less than 3 [START_REF] Grimsley | Cryptic sex in the smallest eukaryotic marine green alga[END_REF]. The minimal frequency of sexual reproduction has been estimated indirectly in Ostreococcus tauri from the population polymorphism spectrum, 1 sexual meiosis for 100,000 clonal divisions [START_REF] Blanc-Mathieu | Population genomics of picophytoplankton unveils novel chromosome hypervariability[END_REF]. This low rate may be explained either by a low encounters rate of MT+ and MT-strains in their natural environment, or by a low rate of outcrossing. The identification of informative sequences from the MT locus in 3 of the 4 SAGs demonstrates the utility of the single cell genomics approach for estimating the frequencies of the two mating types within a bloom, provided that a sufficient number of cells is analysed. While the presence of identical MTs in 3 SAG is consistent with clonal reproduction within a bloom, the sample size was too small for a proper estimation of the level of clonality within a bloom.

The estimations of the percentages of reads recruited to the mitochondrion and chloroplast genomes in each SAG attest the ability of single cell sequencing to identify novel candidate associations in the smallest eukaryotic microalgae of the order Mamiellales. Indeed, while organellar copy numbers may reach several thousands of copies in humans [START_REF] Miller | Precise determination of mitochondrial DNA copy number in human skeletal and cardiac muscle by a PCR-based assay: lack of change of copy number with age[END_REF] or Arabidopsis (79), bacterial-sized Ostreococcus have a reduced number of organellar genomes, with an average copy number per cell estimated at 4 and 5 for the mtDNA and cpDNA, respectively. As a consequence, the systematic identification of mtDNA and cpDNA sequences from the SAGs provides support to the ability of this strategy to detect genomic sequences of any intracellular symbionts, provided there is no DNA amplification bias as a consequence of differences in genome architecture or cell wall structure between the symbiont and the organellar genome. Alternatively, these sequences might arise from random capture of amplified dissolved DNA in the marine seawater around the isolated picoalgal cells. In a typical sorting run the droplet containing a single cell may contain about 2.8 nL (80), or 2.8 × 10 6 µm3. The sample volume within a drop is around 1/1000 of the total volume, or 2.8 × 10 3 µm3, whereas a single Ostreococcus cell contains about 0.5 µm3 of cytoplasm, giving about 6 × 10 3 fold more seawater than cell volume, ample space for the "phycosphere" [START_REF] Seymour | Zooming in on the phycosphere: the ecological interface for phytoplankton-bacteria relationships[END_REF] or dissolved DNA (82) that may be sequenced with the targeted algal cell.

To distinguish between random capture of dissolved DNA and intracellular symbionts, relative read coverage of the non-target to the target genomes provides a rule of the thumb about the likely origin of the non-target sequences. In the present data, we found that most of the non-target sequence data can be traced back to cross-contamination from another SAG, sequenced in the same Illumina run. Filtering out cross-contaminants, the amount of non-target DNA identified -except for the virus sequence -is not compatible with cosampling of a Mamiellales and a bacteria or Oomycete cell, and rather suggests random capture of dissolved DNA. These associated non-target sequences may nevertheless provide information about the diversity of microorganisms present in the cell's environment. Indeed, Mamiellales are auxotrophic for B vitamins [START_REF] Paerl | Use of planktonderived vitamin B1 precursors, especially thiazole-related precursor, by key marine picoeukaryotic phytoplankton[END_REF] and have been found

to co-occur with diverse bacteria in lab cultures (29)(84), though no obligate specific association has been yet observed. The taxonomic diversity of the candidate associated bacterial contigs retrieved from the SAG data is thus not surprising.

The presence of a novel CRESS-like DNA virus came as a surprise. Single-stranded DNA (ssDNA) viruses with circular genomes, with sizes from ~1700 bases, are the smallest viruses known to infect eukaryotes [START_REF] Rosario | Diverse circovirus-like genome architectures revealed by environmental metagenomics[END_REF]. Several previous metagenomic studies have identified many novel genomes similar to ssDNA circular viruses through data-mining of public marine metagenomes (85)( 86), but with no hints about their putative hosts.

Previously described CRESS viral hosts range from plants (87) (Geminiviridae, Nanoviridae), bleached kelp macroalgae [START_REF] Beattie | Novel ssDNA Viruses Detected in the Virome of Bleached, Habitat-Forming Kelp Ecklonia radiata[END_REF], vertebrates (89) (Circoviridae), fungi, insects (Genomoviridae) [START_REF] Krupovic | Genomoviridae: a new family of widespread single-stranded DNA viruses[END_REF], and crustaceans [START_REF] Rosario | Virus discovery in all three major lineages of terrestrial arthropods highlights the diversity of singlestranded DNA viruses associated with invertebrates[END_REF]. In protists, CRESS viruses were discovered in SAGs targeting Picobiliphytes [START_REF] Yoon | Single-cell genomics reveals organismal interactions in uncultivated marine protists[END_REF] and in Bacillariophyceae diatoms, where they have been isolated and their lytic effect could be demonstrated [START_REF] Shirai | Isolation and characterization of a single-stranded RNA virus infecting the marine planktonic diatom Chaetoceros tenuissimus Meunier[END_REF]. The Mamiellales associated circular viruses described here are among the smallest viral genomes described to date infecting the smallest photosynthetic eukaryotes.

In conclusion, sorting of these 12 Mamiellales SAGs by flow cytometry relied on the choice of appropriate volumes to favour isolation of one cell, and the "phycosphere" around cell had been subjected to turbulent mixing with the cytometer sheath fluid, thus limiting the suitability of this technique for identification of cell-to-cell interactions. Independent ways of finding such associations are necessary, one of which may be to pick off single cells using micromanipulation (see [START_REF] Galindo | Combined cultivation and single-cell approaches to the phylogenomics of nucleariid amoebae, close relatives of Fungi[END_REF] in this issue), but the technical challenge of performing this for picoeukaryotes < 2 µm in size, such those in the Mamiellales, remain very challenging.

Conventional culture-based approaches of strains and viruses are required to demonstrate the outcome of interactions by co-culture [START_REF] Amin | 4;522(7554):98-101. Reference Genome Nuclear (GC%) Mitochondria (bp) (GC%) Chloroplast (bp) (GC%) B. prasinos[END_REF] and remain an essential tool for deciphering both intercellular and virus-host interactions in aquatic environments and interpreting environmental sequencing data. Nevertheless, we showed here that -provided careful filtration of contaminations between dataset is performed -SAGs provide in-situ ecological information that would be lost in conventional isolation and culture processes, which select for the fastest growing cells in laboratory culture medium.

Tables

(

  Mamiellophyceae) were performed as follows. Each contig was aligned against GenBank using DIAMOND (blastx option e-value 1e-5 maximum aligned sequences 100) and the 10 best hits for each contig were kept. The taxonomic affiliation of each contig was retrieved using a custom perl script that fetched the GenBank accession number from the DIAMOND output, and retrieved the taxonomical nodes from the GenBank taxonomy database for each hit (script available upon request). Manual inspection of taxonomic affiliation of contigs containing several open reading frames led to the correction of the taxonomic affiliation for a few contigs affiliated to Archaea on the base of the best blast hit, while most open reading frames along the contig had a best hit against Bacterial sequences. As the complete genome and organellar sequences of Ostreococcus, Bathycoccus, and Micromonas species are present in GenBank, contigs were considered as target if at least one of the 10 best hits belonged to a

RCC1105 genome led the authors to conclude that these 4 3 Ostreococcus

 43 SAGs belonged to a cryptic novel Bathycoccus species, TOSAG39, which is supported here by the phylogenetic divergence of their mitochondrial genome sequences. Although combining SAGs produces a more complete genome, individual information about the environment of the sorted single cells, which could indicate potential biological interactions, or even stages of infections in the case of parasites, may be lost.Although the individual SAGs were estimated to represent a maximum of 27% of the complete genome sequence in Ostreococcus, a search for the presence of 23 gene families located in the mating-type region in Mamiellales enabled us to retrieve orthologous genes for 3 Ostreococcus SAGs. Amino-acid identity between the SAG sequences and the sequences from available genomes indicated that the candidate mating type of the SAGs was MT+, like the candidate mating type of the closest reference genome from strain O. sp RCC809[START_REF] Blanc-Mathieu | Population genomics of picophytoplankton unveils novel chromosome hypervariability[END_REF]. The relative prevalence of the MT+ and MT-in natural Ostreococcus blooms remains an open question. To the best of our knowledge, the cooccurrence of strains of the two mating types from the same sample has not yet been

FIG1-

  FIG1 -Maximum likelihood phylogeny (under GTR model) of 18S rDNA sequences retrieved from the SAG assemblies of Bathycoccus sp. (SAG1-3), Micromonas sp. (SAG5 and 8) and Ostreococcus sp. (SAG9-12). Mamiellales reference genomes are marked with a red star, and SAGs are marked with orange star. Rapid bootstrap support values are indicated in the key (100 replicates), represented by colored circles in branches. The final alignment contained 1750 nucleotides.

FiguresFIG1-

  Figures

FIG2-A 7 An o p h o re c ir c u la r v ir u s 2 Yerb a mate circu lar DNA virus 1 ACR ES S vir us sp 2 D 6 CR 1 C te n o p h o re c ir c u la r v ir u s 4 Po rc in e ci rc ov iru s 2 2 A 4 C 9 O 1 0 12 SAGFIG5-

 72126142249112 FIG2 -Mating type inference from amino-acid identity with MT+ and MT-orthologous gene families (GF) in Ostreococcus SAGs (a) Expected taxonomic affiliation for MT-and MT+ strains (b) Presence/absence matrix of hits to gene families in each SAGs (c) Maximum Likelihood phylogeny (LG+G model) of the SAGs sequence corresponding to GF16 with orthologous genes from Ostreococcus, Micromonas and Bathycoccus species.

table 1

 1

	). Contigs

phylogenetic trees were built under GTR model using RAxML v. 8.2.9 (51) and plotted with FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/). Taxonomic affiliations to separate non-target contigs (foreign DNA) from target

Table 1 : Assembly statistics of 12 SAGs and statistics of reference genomes from the order

 1

	Mamiellales
	Sample -station	Bathycoccus -39		Micromonas -46	Ostreococcus -41
	SAG	1	2	3	4	5	6	7	8	9	10 11 12
	Assembly size (Mb)	3.25 4.42 3.36 3.72 0.89 0.56 0.58 0.53 0.85 3.46 1.54 2.71
	N50 (kb)	13.6 17.4 19.5 20.2 3.17 2.47 9.67 6.91 6.44 37.5 17.8 36.4
	GC%	45.7 47.0 46.6 46.5 40.0 39.3 41.6 39.1 52.2 57.1 52.1 56.6
	Contigs	3157 3129 3461 2913 2864 1894 1043 1788 2914 2479 2207 2364
	Contigs (> 400 bp)	1048 1241 990 991 691 470 272 275 410 737 546 655

Table 2 :

 2 Summary statistics of assemblies per taxonomic affiliation categories, assembly length and percent of reads recruited are provided for Target (Mamiellophyceae) and Non Target assemblies (Bacteria, Eukaryote, Archaebacteria and virus).

	* Tara-Oceans Expedition station number

Table 3 :

 3 Summary statistics of cpDNA and mtDNA manual assemblies from combined SAG sequence data.

				Target sequences			Non Target sequences		BUSCO
			Assembly total length (% of reads recruited)	Total length (% of reads recruited)
	T.-O.*	SAG	Total	Mitochondria	Chloroplast	Total		Filtered
	Station		x10 3 bp (%)	bp	(%)	bp	(%)	bp	(%)	bp	(%)	(%)
	39	SAG1	3,001 (99.89)	40,481	(1.46)	33,423	(0.03)	101,751 (0.09)	29,645 (0.03)	11.9
	39	SAG2	4,166 (99.47)	49,421	(0.05)	21,954	(0.01)	89,862 (0.28)	22,815 (0.07)	17.2
	39	SAG3	3,130 (99.91)	46,597	(0.11)	35,348	(0.05)	90,824 (0.06)	46,751 (0.03)	15.2
	39	SAG4	3,509 (99.96)	54,044	(0.08)	38,960	(0.04)	118,507 (0.02)	32,028 (0.005)	16.8
	46	SAG5	441 (90.32)	54,577	(3.36)	51,112	(65.97)	228,342 (2.15)	122,933 (1.16)	2.0
	46	SAG6	271 (98.26)	37,972	(20.66)	64,962	(48.70)	135,783 (0.26)	33,796 (0.06)	1.0
	46	SAG7	255 (96.12)	34,404	(0.28)	53,356	(73.30)	136,620 (2.02)	26,962 (0.39)	2.3
	46	SAG8	187 (82.95)	44,825	(14.67)	42,081	(0.21)	144,189 (6.09)	31,813 (1.34)	3.0
	41	SAG9	712 (99.86)	43,157	(1.58)	40,737	(0.15)	41,906 (0.08)	6,307 (0.01)	28.7
	41	SAG10	3,180 (99.88)	47,961	(1.11)	93,684	(0.79)	55,576 (0.01)	11,764 (0.002)	22.8
	41	SAG11	1,225 (92.74)	28,820	(19.66)	101,901	(7.98)	142,887 (1.78)	20,493 (0.25)	7.6
	41	SAG12	2,250 (98.85)	43,041	(0.15)	84,203	(10.60)	79,669 (1.11)	11,876 (0.16)	15.5

Acknowledgments

We would like to thank Tom Delmont and Olivier Jaillon from the Genoscope for information about the SAGs, Lilian Caesar for help with the figures, and the Genotoul platform (www.genotoul.fr) for cluster availability for bioinformatics analysis and Genoscope for sequencing. We are grateful to Fabien Burki for making a perl script for extracting taxonomic affiliation from DIAMOND outputs available. Tara Oceans (which includes both the Tara Oceans and Tara Oceans Polar Circle expeditions) would not exist without the leadership of the Tara Expeditions Foundation and the continuous support of 23 institutes (http://oceans.taraexpeditions.org). We further thank the commitment of the following sponsors: CNRS (in particular Groupement de Recherche GDR3280 and the Research Federation for the study of Global Ocean Systems Ecology and Evolution, FR2022/Tara Oceans-GOSEE), European Molecular Biology Laboratory (EMBL), Genoscope/CEA, The French Ministry of Research, and the French Government 'Investissements d'Avenir' programmes OCEANOMICS (ANR-11-BTBR-0008), FRANCE GENOMIQUE (ANR-10-INBS-09-08), MEMO LIFE (ANR-10-LABX-54), and PSL* Research University (ANR-11-IDEX-0001-02), Gordon and Betty Moore Foundation (award #3790) and the US National Science Foundation (awards OCE#1536989 and OCE#1829831) to MBS, as well as the Ohio Supercomputer for computational support. We also thank the support and commitment of agnès b. and Etienne Bourgois, the Prince Albert II de Monaco Foundation, the Veolia Foundation, Region Bretagne, Lorient Agglomeration, Serge Ferrari, Worldcourier, and KAUST. The global sampling effort was enabled by countless scientists and crew who sampled aboard the Tara from 2009-2013, and we thank MERCATOR-CORIOLIS and ACRI-ST for providing daily satellite data during the expeditions. We are also grateful to the countries who graciously granted sampling permissions. The authors declare that all data reported herein are fully and freely available from the date of publication, with no restrictions, and that all of the analyses, publications, and ownership of data are free from legal entanglement or restriction by the various nations whose waters the Tara Oceans expeditions sampled in. This article is contribution number XX of Tara Oceans.

LFB was funded by the EU Horizon 2020 research and innovation programme, under the Marie Skłodowska-Curie grant agreement No. H2020-MSCA-ITN-2015-675752. We would also like to acknowledge three anonymous referees for constructive comments on a previous version of this manuscript.

Supplementary material

Supplementary material table 1. GenBank accession numbers of 18S sequences used for phylogeny.

Supplementary material table 2. GenBank accession numbers for virus sequences used in Fig. 5. 3. Taxonomic affiliation of non-target bacterial sequences for each SAG.

Supplementary material table

Supplementary material table 4. Plastid, mitochondrion and nuclear genome coverage in 12

Ostreococcus tauri strains and estimation of mitochondrial and plastid genome copy number per haploid nuclear genome.

Additional Information

Ethics

Not applicable.

Data Accessibility

SAGs assemblies analysed in this manuscript have been submitted to GenBank under BioProject PRJNA549236 with accession numbers VIAS00000000 (SAG1); VIAT00000000 (SAG2); VIAU00000000 (SAG3); VIAV00000000 (SAG4); VIAW00000000 (SAG5); VIAX00000000 (SAG6); VIAY00000000 (SAG7); VIAZ00000000 (SAG8); VIBA00000000 (SAG9); VIBB00000000 (SAG10); VIBC00000000 (SAG11); VIBD00000000 (SAG12). Mitochondrial assemblies have been submitted to GenBank under accession MN243825 (Os_mt_TOSAG41), MN243826 (Mi_mt_TOSAG46), plastidial assemblies under accession MN243827 (Os_cp_TOSAG41) and MN243828 (Mi_cpTOSAG46).

Authors' Contributions

LFB and GP performed bioinformatics analyses. NP and MS performed single cell sorting and prior species' taxonomic affiliation. KL performed SAG sequencing. GP, LFB and NG designed the study. All authors contributed to manuscript writing and editing.

Competing Interests

We have no competing interests.