Yannick Vimont
email: yannick.vimont@ema.fr

Sylvain Boussier
email: sylvain.boussier@ema.fr

Michel Vasquez
email: michel.vasquez@ema.fr

Reduced costs propagation in an efficient implicit enumeration for the 01 multidimensional knapsack problem

Keywords: Multidimensional knapsack problem, Implicit enumeration, Variable fixing, Reduced costs, Constraint propagation 1 Introduction

In a previous work we proposed a variable fixing heuristics for the 0-1 Multidimensional knapsack problem (01MDK). This approach uses fractional optima calculated in hyperplanes which contain the binary optimum. This algorithm obtained best lower bounds on the OR-LIBRARY benchmarks. Although it is very attractive in terms of results, this method does not prove the optimality of the solutions found and may fix variables to a non-optimal value. In this paper, we propose an implicit enumeration based on a reduced costs analysis which tends to fix nonbasic variables to their exact values. The combination of two specific constraint propaga-tions based on reduced costs and an efficient enumeration framework enable us to fix variables on the one hand and to prune significantly the search tree on the other hand. Experimentally, our work provides two main contributions: (1) we obtain sev-eral new optimal solutions on hard instances of the OR-LIBRARY and (2) we reduce the bounds of the number of items at the optimum on several harder instances.

follow:

01MDK

Maximize c • x subject to,

A • x ≤ b and x ∈ {0, 1} n ,
where c ∈ N n , A ∈ N m×n and b ∈ N m . The binary components x j of x are decision variables: x j = 1 if the item j is selected, 0 otherwise. c j is the profit associated to the item j and A ij is the cost (in terms of resource i) of the selecting item j . b i is the amount of the resource i.

The designation comes from the original formulation which describes the problem of a hiker who has to select a subset of items from a list to take in his knapsack, so that the selected items are the most useful possible and the total weight does not exceed a given limit. It can be useful in several practical problems such as budgeting problems (investments problems), manufacturing problems (paper cutting patterns, cargo or plane loading) or telecommunication problems. The reader is referred to [START_REF] Fréville | The multidimensional 0-1 knapsack problem: an overview[END_REF] for more details about the different existing methods for solving this problem.

In a previous work, [START_REF] Vasquez | Improved results on the 0-1 multidimensional knapsack problem[END_REF] proposed a variable fixing heuristics for the 01MDK. This algorithm obtained best results on Chu and Beasley instances available on the OR-LIBRARY. 1 Although it is very attractive in terms of results, this method does not prove the optimality of the found solutions and may fix variables to a non optimal value. [START_REF] Oliva | Constraint and linear programming: using reduced costs for solving the zero/one multiple knapsack problem[END_REF] propose an interesting branch & bound method for solving the 01MDK. They use a reduced costs constraint based on the reduced costs of the non basic variables and the slack variables at the optimum of the relaxed 01MDK. Being inspired by this approach, we propose in this paper an implicit enumeration for the 01MDK which tends to fix variables to their optimal value.

At first, we will explain the main principle of our algorithm in Sect. 2 and we will detail the reduced costs based constraint we use in Sect. 3. Then we will describe the implicit enumeration framework in Sect. 4 and the reduced costs propagation embedded in Sect. 5. The computational results made on benchmarks from the OR-LIBRARY are detailed in Sect. 6.

General principle

In this paper, we propose an implicit enumeration which uses a reduced costs constraint to fix non-basic variables and prune nodes of the search tree. The reduced costs constraint expresses the objective function of the LP-relaxation with reduced cost coefficients and bounds this expression with any available lower bound and the optimal objective value.

Several algorithms which use the reduced costs for solving knapsack problems have been published before, the reader is referred to [START_REF] Saunders | A shrinking boundary algorithm for discrete system models[END_REF], [START_REF] Fayard | An algorithm for the solution of the 0-1 knapsack problem[END_REF] and [START_REF] Oliva | Constraint and linear programming: using reduced costs for solving the zero/one multiple knapsack problem[END_REF] for articles on this subject. The originality of the proposed algorithm lies on the fact that the unpromising parts of the search tree are tackled at first. Indeed, rather than keeping the branching variable to its non-basic value, we explore first the tree that corresponds to fixing it to its opposite value. Let us note that the aim of this approach is to prune the search tree as soon as possible. Hence this algorithm should be more appropriated for solving difficult instances of 01MDK than easy ones.

The method exposed here is divided in three main phases and uses several previously published results: (1) the calculation of a starting lower bound, which is carried out by the local search algorithm published in [START_REF] Vasquez | Improved results on the 0-1 multidimensional knapsack problem[END_REF], (2) the decomposition of the problem in several sub-problems based on the observation of [START_REF] Vasquez | An hybrid approach for the 0-1 multidimensional knapsack problem[END_REF] and (3) the enumeration of each sub-problem embedding the specific reduced cost propagation.

The reduced costs constraint

Let us consider the LP-relaxation of the 01MDK stated as follows:

(01MDK) Max z t = UB + j ∈N - cj x j - j ∈N + cj (1 -x j ) + j ∈M ūi s i subject to Āx + Ss = b, x ∈ [0, 1] n , s ≥ 0.
Where s is the vector of the m slack variables, ( x, s) an optimal solution of the linear program 01MDK and UB its value, ( c, ū) the vector of the reduced costs corresponding to the variables (x, s) for the basic solution ( x, s). Ā, S, b are values corresponding to the base associated to ( x, s), so that S is the identity matrix. N -represents the indexes of the non-basic variables, which are equal to their lower bound and N + the indexes of the non-basic variables which are equal to their upper bound. Let us suppose that we know a lower bound LB ∈ N of the original integer linear program 01MDK, then each current solution better than LB has to satisfy the following relation:

UB + j ∈N - cj x j - j ∈N + cj (1 -x j ) + i∈M ūi s i ≥ LB, - j ∈N - cj x j + j ∈N + cj (1 -x j ) - j ∈M ūi s i ≤ UB -LB.
In our case, we do not take into account the slack variables. Indeed, at the optimum, the slack variables are positive and their associated reduced costs negative, thus the sum j ∈M ūi s i is always negative. Hence, we can withdraw the slack variables from the inequality mentioned above without losing sense. Moreover, we also know that at the optimum of the linear program, the non-basic variables with a positive reduced cost are equal to their upper bound ( cj > 0 ⇒ x j ∈ N + ) and the non-basic variables with a negative reduced cost are equal to their lower bound ( cj < 0 ⇒ x j ∈ N -). Hence, if we consider the absolute values of the reduced costs, we obtain the so called "reduced costs constraint":

j ∈N - | cj |x j + j ∈N + | cj |(1 -x j ) ≤ UB -LB.
(1)

Implicit enumeration framework

In this section, we only describe the implicit enumeration framework. The reduced costs propagation will be detailed in Sect. 5. The first phase, which consists in calculating a starting lower bound, is detailed in Sect. 6. Once we obtain the lower bound, the second phase consists in partitioning the problem according to the observation of [START_REF] Vasquez | An hybrid approach for the 0-1 multidimensional knapsack problem[END_REF]: they highlighted the fact that the number of items k can be easily bounded with two values k min and k max . Calculating the (k maxk min + 1) different values of k using this method, we partition the problem in 01MDK(k) problems for k ∈ {k min , . . . , k max } so that:

01MDK(k) ⎧ ⎪ ⎨ ⎪ ⎩ Maximize c • x subject to, A • x ≤ b and x ∈ {0, 1} n , 1 • x = k ∈ N.
The objective of the proposed algorithm is to solve each 01MDK(k) independently and to return the best value found.

For each sub-problem, the enumeration is tackled as classical branch & bound algorithms, with the main difference that the branching is done on the non-basic (i.e. integer) variables and not on the most fractional one as we usually do.

Let us consider S t a partial solution (or node of the search tree) so that t is the depth parameter. A partial solution S t is represented by a vector (x 1 , . . . , x n ) where x i ∈ {0, 1, * } n and x i = * if x i is a free variable. S t is composed of a subset F (S t ) of variables fixed to zero or one and a subset L(S t ) of free variables. F (S t ) is partitioned in two subsets F 0 (S t ) and F 1 (S t ), containing respectively variables fixed to zero and variables fixed to one. Let 01MDK S t be the linear relaxation of the sub-problem at node S t , where decision variables are free variables and let zS t be its optimal value: the evaluation of S t is defined as follows:

UB(S t ) = j ∈F 1 (S t ) c j + zS t .
(2)

At each node S t of the search tree, the upper bound UB(S t ) is computed. If S t is unfeasible or if UB(S t ) is lower or equal to the original lower bound LB, the enumeration backtracks. In the other case, the branching is done on the non-basic and free variable x i at the optimum of the 01MDK S t with the bigger absolute reduced cost value by fixing it to the opposite of its optimal value (the reason why we branch on the variable with the bigger absolute reduced cost value will be explained in Sect. 5.2).

Once we explore the corresponding partial solution, we are able now to fix x i to its optimal value xS t i and to fix the next non-basic and free variable with the bigger absolute reduced cost value x j to the opposite of its optimal value (1 -xS t j ). Indeed, it is useless to explore the partial solution S t+1 in which x i is fixed to its optimal value since, in that case, 01MDK S t+1 provides the same optimal solution as 01MDK S t . The partial solution where all the non-basic variables are fixed to their optimal values and all the basic variables are free, is explored only once at each node. For example, let us suppose that xS t = (0, 1, 0, 0, 0.5, 0.4) is the optimal solution of 01MDK S t and that the variables are sorted by decreasing order of absolute reduced cost value. Then the nodes generated will be the followings:

(1, * , * , * , * , * ) (0, 0, * , * , * , * ) (0, 1, 1, * , * , * ) (0, 1, 0, 1, * , * ) (0, 1, 0, 0, * , * ).

In the case where all non-basic variables are fixed, and all basic variables are free: a Depth First Search algorithm enumerates the remaining basic variables and if a feasible solution is met, the best known solution is updated. If the sum of the variables fixed to one is equal to k (the number of items corresponding to the explored hyperplane 1 • x = k) and the solution is feasible, the best known solution is updated else the enumeration backtracks.

Reduced costs propagation

In the following, we expose two reduced cost propagations embedded to the implicit enumeration framework. Both of them are based on the same reduced costs analysis but diverge in the way they are carried out. The first one takes the local reduced costs constraint into account at each node of the search tree, and the second one is a more global propagation which deals with the reduced costs constraint of all parents of each node.

Local propagation

Let gap(S t ) = UB(S t ) -LB be the gap value at node S t . The reduced costs constraint (1) allows us to fix several non-basic variables either to zero or to one. Indeed, if we insulate a variable x j which is non-basic at the optimum of the 01MDK S t and consider all others non-basic variables at their bounds, i.e. zero for the variables with a negative reduced cost and one for the variables with a positive reduced cost, the constraint (1) gives the two following relations:

( xS t j = 1) ⇒ x j ≥ 1 - gap(S t ) | ct j | , ( 3
)
( xS t j = 0) ⇒ x j ≤ gap(S t ) | ct j | (4)
where xS t j is the value of the non-basic and free variable x j at the optimality of the 01MDK S t and cS t j its corresponding reduced cost. Using relations (3) and (4), we state the following local reduced costs propagation:

• For each variable x j which is non-basic at the optimum of the 01MDK S t , if | cS t j | > gap(S t ): Fix x j to zero (resp. one) if and only if xS t j is equal to zero (resp. one).

Thus, if a non-basic variable x j has a reduced cost cS t j > gap(S t ) at node S t , x j is fixed to the value xS t j . Moreover, changing the value of a variable x j implies that the left-hand side of constraint (1) at node S t is augmented by the value cS t j . Then the value gap(S t ) can be interpreted as the amount of reduced cost available to swap the value of the nonbasic variables at node S t . The current partial solution S t is consistent with constraint (1) if and only if j ∈ S t cS t j ≤ gap(S t ), so that S t = {j / x j is fixed to 1 -xS t j }. Supposing we branch on x j at node S t , then all the other variables x i ∈ L(S t ) so that cS t j + cS t i > gap(S t ) cannot be fixed to the opposite of their optimal value. So, at each node of the search tree, we are able to proceed to this second reduced costs propagation:

• Supposing we branch on x j ∈ L(S t ). For each non-basic and free variable x i ∈ L(S t ) i = j , if cS t j + cS t i > gap(S t ): fix x i to zero (resp. one) if xS t i is equal to zero (resp. one) and let it free otherwise.

For example, let us suppose that xS t = (0, 1, 0, 0, 1, 0.5, 0.4) at the optimality of the 01MDK S t and that cS t 1 , cS t 2 > gap(S t ), then we fix x 1 to 0 and x 2 to 1 according to the local reduced costs propagation and we generate the following children:

(0, 1, 1, * , * , * , * , * ) (0, 1, 0, 1, * , * , * ) (0, 1, 0, 0, 0, * , * ) (0, 1, 0, 0, 1, * , * ).
Let us suppose now that c 3 + c 5 > gap(S t ), then according to the second reduced costs propagation exposed above, x 5 will be directly fixed to its optimal value 1 after the branching on x 3 and the first children generated will be (0, 1, 1, * , 1, * , * ) instead of (0, 1, 1, * , * , * , * , * ). This second propagation enables us to reduce the number of simplex phases in the enumeration.

Global propagation

This propagation comes from the observation that each 01MDK S t provides a reduced costs constraint and that each child of S t has to satisfy this constraint. In their approach, [START_REF] Oliva | Constraint and linear programming: using reduced costs for solving the zero/one multiple knapsack problem[END_REF] use the information brought by those constraints: at each node S t , the corresponding reduced costs constraint is added to a constraint programming model and propagated on the variables domains. This process appears to be high c.p.u. time consuming and to penalize the whole computing time. In the proposed approach, the propagation only deals with reduced costs propagation in terms of available gap quantity as in the local propagation.

Indeed, let P (S t ) be the set of parent nodes of S t . Let us suppose that x j is fixed to 1 -xS t j at node S t and that xS p j = xS t j , p ∈ P (S t ), then we can update gap(S p ) by withdrawing | cS p j |. If it appears that gap(S p ) < 0, S t becomes unfeasible. Supposing x j is fixed to the value v at node S t , then for each parent node S p p ∈ P (S t ) so that xS p j = 1v, the value gap(S p ) is reduced by the corresponding reduced cost | cS p j | and if gap(S p ) < 0, the current partial solution S t is unfeasible.

Regarding the reduced cost propagation scheme, it appears clearly that the more the variable we branch on has a high absolute reduced cost value, the more the gap value is thin and the propagation is effective. Hence, at each node, we first branch on the free variables which have the higher reduced cost value in order to make the algorithm converge faster.

Algorithms 1 and 2 respectively detail the recursive enumeration procedure. Algorithm 1 is devoted to the variable fixing and the branching phase, Algorithm 2 to the creation of sub-problems and to the reduced costs propagation.

Algorithm 1 Enumeration

This method enables us to propagate the reduced costs without adding explicitly new constraints to the model at each node. We only consider one reduced cost constraint and the consistency of each sub-problem solution with it. The propagation is carried out by ascending the search tree from the current node to the root. At each node, we check that the reduced costs constraint is consistent with the fixed variables otherwise the current node is pruned.

Let us note that contrary to the classical Branch & Bound method, this propagation proscribes us to update the lower bound when a better feasible solution is met except at the root of the tree. Indeed, if the lower bound of a node of the branch and bound tree is different from the lower bound of its parents, the global propagation will have no meaning since the gap value will change. However, since the root node does not have any parents, it is possible to update the lower bound at this stage.

Algorithm 2 Auxiliary program 6 Computational experiments

Our enumeration procedure has been tested on benchmarks produced by [START_REF] Chu | A genetic algorithm for the multidimensional knapsack problem[END_REF] issued from the OR-LIBRARY.2 These benchmarks were generated using the procedure proposed by [START_REF] Fréville | An efficient preprocessing procedure for the multidimensional 0-1 knapsack problem[END_REF]. This collection contains problems with n = 100, 250 and 500 variables and m = 5, 10 and 30 constraints. Thirty problems were generated for each n-m combination, giving a total of 270 problems. The a ij integer values were generated by using a uniform random generator U(0, 1000) and the right-hand side values were generated in correlation with the a ij so that b i = α n j =1 a ij where α = 1/4, 1/2 and 3/4 is the tightness ratio. The costs of the objective function (c j ) were generated in correlation with the a ij so that c j = m i=1 a ij /m + 500d j where d j is a random generated number between 0 and 1. The whole name of each problem is cbm.n_r, where m is the number of constraints, n the number of variables and r the number of the instance. Each set of 30 instances with the given numbers m and n is divided into 3 series with a different tightness ratio: problems 1 to 9 correspond to the ratio α = 1/4, problems 10 to 19 to the ratio α = 1/2 and problems 20 to 29 to the ratio α = 3/4. Our algorithm has been tested on a P4 3.2 GHz with 1 GB RAM and the results are compared to those produced by the commercial solver CPLEX 9.2.

Lower bound calculation

The lower bounds were calculated separately with two previously published algorithms: (1) RL tabu proposed by [START_REF] Vasquez | An hybrid approach for the 0-1 multidimensional knapsack problem[END_REF], which is a hybrid approach combining linear programming and an efficient tabu search and (2) Fix + RL tabu proposed by [START_REF] Vasquez | Improved results on the 0-1 multidimensional knapsack problem[END_REF], which is the same algorithm embedding a variable fixing heuristics. We computed the lower bounds with RL tabu for the 100 variable instances and with Fix + RL tabu for the 250 and 500 variable instances. The reason why we use both algorithms is that the variable fixing heuristics embedded in Fix + RL tabu makes the algorithm consume more c.p.u. time and that RL tabu performs quite the same results for the smaller instances. It is important to note that we limited the number of iterations for both algorithms in order to have an appropriate ratio time consumption/quality of the solutions. This explains why the lower bounds exposed in this paper are lower than the ones published in [START_REF] Vasquez | Improved results on the 0-1 multidimensional knapsack problem[END_REF]. For the 100 variable instances, RL tabu is limited to 50 000 iterations and Fix + RL tabu is limited to 100 000 iterations for the 250 variable instances and 200 000 iterations for the 500 variable instances.

Computational results

Our enumeration algorithm found the optimum of all the 10 constraint, 250 variable instances (Table 2) even if the c.p.u. time consumption exceed 4 hours. These instances are hard to solve and we expose here new results never obtained before. CPLEX solves 46% of them in less than 4 hours of computational time and run out of memory for the others 54%. Our enumeration algorithm solves 94% of them in the same time limit. Our algorithm proves all the 500 variable, 5 constraint instances in less than 4 hours of computational time. However the optimal solutions were already published by [START_REF] James | Enumeration methods for repeatedly solving multidimensional knapsack sub-problems[END_REF].

Our algorithm consumes little RAM memory when solving hard instances and this enables us to solve instances to the optimum even if the required computational time is greater than 4 hours (the maximum time is around 24 hours for the cb10.250_7 instance). Indeed, for the majority of the 10 constraint, 250 variable instances, CPLEX exceeds the capacity of RAM memory after 4 hours of computational time when our enumeration solves all the instances. Another positive point of our algorithm is that its structure enables to use parallel processors to reduce the overall run time required to solve a given instance. Indeed, we can assign the resolution of each 01MDK(k) with k ∈ [k min , k max ] to each processor. The parallel approach can widely decrease the run time required to solve hard instances. Tables 1, 2 and 3 detail the results obtained on problems with 250 and 500 variables. Bold face highlights the best c.p.u. time required between our algorithm and CPLEX. The description of the columns is:

• instance: the name of the instance; • z * : the value of the lower bound; • t * : the c.p.u. time in seconds required to calculate the lower bound; • z opt : the optimal solution found by the enumeration; • gap(×10 2 ): the gap between the lower bound and the optimal value; • k opt : the number of items in the optimal solution; best order to explore hyperplanes in a limited time, is to start from the possible extreme integer values of [k min , k max ] to the center integer values: k min , k max , k min + 1, k max -1, etc. The objective is to prove that no solution z > z * with k items exists, so we first tackle the "easier" hyperplanes with the smaller value z(k).

We expose on Table 4 the new bounds for the number of items at the optimum so that k min ≤ k opt ≤ k max , boldface highlights the new bounds provided comparatively to those previously known with the method exposed in Sect. 4. The description of the data per column is:

• Pb: the number of the instance.

• LB: the lower bound used. LB is in boldface if the lower bound has been improved and the star (*) indicates that LB is the optimal value. • k min : the lower bound of the number of items at the optimum.

• k max : the upper bound of the number of items at the optimum.

We reduced the bounds of 88% of the instances within a limit of 10 hours of computational time. A hyperplane is cut if our algorithm explores the whole sub-tree of the corresponding value k and does not find any better solution than the best known solution, in the other case, if it finds a better solution than the best known solution, or if it reaches the time limit, it stops the exploration of the current hyperplane.

Conclusion

In this paper we have demonstrated how an implicit enumeration combining constraint and linear programming can be efficient to solve large-scale 01 multidimensional knapsack problems. Our enumeration method includes two specific reduced costs propagations which enable on the one hand to fix variables and on the other hand, to prune nodes of the search tree. Experimentally, our algorithm leads us to good results in terms of computational times. The new optimal solutions never published before and the new bounding of the number of items at the optimum for hard instances constitutes the two main contributions of our approach. Moreover the provided solutions can help researchers in the field of heuristics to evaluate their solutions. For the harder instances which are not solvable in a reasonable c.p.u time, we plan to experiment a distributed version of our algorithm.

Table 1

 1 Results obtained on cb5.250 problems

	Instance	lb		opt					Cplex
		z *	t *	z opt	gap(×10 2 )	k opt	t opt	t total	t Cpx
	cb5.250_0	59312	80	59312	0	73	0	80	33
	cb5.250_1	61472	79	61472	0	74	3	82	156
	cb5.250_2	62130	61	62130	0	76	1	62	11
	cb5.250_3	59453	106	59463	0.016	71	93	199	861
	cb5.250_4	58951	95	58951	0	74	6	101	261
	cb5.250_5	60062	78	60077	0.024	74	28	106	384
	cb5.250_6	60396	98	60414	0.029	76	6	104	171
	cb5.250_7	61449	101	61472	0.037	73	36	137	239
	cb5.250_8	61885	73	61885	0	76	3	76	69
	cb5.250_9	58959	87	58959	0	72	0	87	16
	cb5.250_10	109086	96	109109	0.021	132	15	111	124
	cb5.250_11	109841	97	109841	0	135	2	99	40
	cb5.250_12	108508	72	108508	0	130	3	75	129
	cb5.250_13	109378	98	109383	0.004	134	12	110	157
	cb5.250_14	110718	95	110720	0.001	130	21	116	671
	cb5.250_15	110256	91	110256	0	132	5	96	227
	cb5.250_16	109040	79	109040	0	133	3	82	195
	cb5.250_17	109042	95	109042	0	132	14	109	314
	cb5.250_18	109971	100	109971	0	130	5	105	165
	cb5.250_19	107058	90	107058	0	131	6	96	18
	cb5.250_20	149659	91	149665	0.004	191	5	96	101
	cb5.250_21	155940	92	155944	0.002	190	3	95	62
	cb5.250_22	149334	95	149334	0	191	24	119	141
	cb5.250_23	152130	103	152130	0	193	2	105	94
	cb5.250_24	150353	98	150353	0	191	4	102	90
	cb5.250_25	150045	73	150045	0	191	0	73	9
	cb5.250_26	148607	75	148607	0	192	0	75	6
	cb5.250_27	149772	95	149782	0.006	193	8	103	127
	cb5.250_28	155061	91	155075	0.009	190	1	92	17
	cb5.250_29	154662	95	154668	0.003	191	4	99	121

Actually, the instances have been generated by using the procedure proposed by[START_REF] Fréville | An efficient preprocessing procedure for the multidimensional 0-1 knapsack problem[END_REF].

Available at http://people.brunel.ac.uk/~mastjjb/jeb/info.html.

This property is directly proved by applying the theorem 10.2 of the book Linear Programming by[START_REF] Chvátal | Linear programming[END_REF].

• t opt : the c.p.u. time in seconds required to find the optimal solution;

• t total : the total time required (t * + t opt );

• t Cpx : the c.p.u. time in seconds required to CPLEX 9.2 to find the optimal solution or X if it exceeds 4 hours of computational time;

6.3 Bounds reduction of the number of items at the optimum For the hard instances (cb30.250, cb10.500 and cb30.500), we attempted to reduce the bounds of the number of items at the optimum by experimenting our algorithm within a limit of 10 hours of computational time. Indeed, even if we stop it prema- turely, the enumeration can provide interesting results: (1) It can fix variables to their optimal value in each explored hyperplane and (2) it can explore all the sub-tree corresponding to 1 • x = k and prove that no solution with k number of items better than the lower bound exists. The main interest here is to reduce the bounds of the number of items at the optimum (i.e. to cut some hyperplanes). Hence, we thought about the best way to explore hyperplanes within a given time limit. Considering the fact that the function z(μ) = c • xμ (optimum of 01MDK(μ) where μ ∈ R) is convex, 3 the


Denise Rui

Ying Ong

Andrés Gutiérrez-Rodríguez

Laurence Garczarek

Dominique Marie

Adriana Lopes

Dos Santos

Nested PCR approach for petB gene metabarcoding of marine Synechococcus

Keywords: marine picocyanobacteria, marine Synechococcus, petB, nested PCR, metabarcoding

The molecular diversity of marine picocyanobacterial populations, an important component of phytoplankton communities, is better characterised using high-resolution marker genes compared to the 16S rRNA gene as it has greater sequence divergence to differentiate between closely related picocyanobacteria groups. Although specific ribosomal primers have been developed, another general disadvantage of bacterial ribosome-based diversity analyses is the variable number of rRNA gene copies. To overcome these issues, the single-copy petB gene, encoding the cytochrome b 6 subunit of the cytochrome b 6 f complex, has been used as a high-resolution marker gene to characterise Synechococcus diversity. We have designed new primers targeting petB gene and proposed a nested polymerase chain reaction (PCR) method (Ong_2022) for metabarcoding of marine Synechococcus populations obtained by flow cytometry cell sorting. We evaluated the specificity and sensitivity of Ong_2022 against the standard amplification protocol (Mazard_2012) using filtered seawater samples. Ong_2022 approach was also tested on flow cytometry sorted Synechococcus populations. Samples (filtered and sorted) were obtained in the Southwest Pacific Ocean, from subtropical (ST) and subantarctic (SA) waters masses. The two PCR approaches on filtered recov-1 Ong et al.

ered the same dominant subclades Ia, Ib, IVa and IVb with small differences in their relative abundance across the distinct samples. For example, subclade IVa was the dominant in ST samples with the Mazard_2012 approach, while the same samples processed with Ong_2022 showed similar contribution of subclades IVa and Ib to total community. Ong_2022 approach generally captured a higher genetic diversity of Synechococcus subcluster 5.1 compared to Mazard_2012, while having a lower proportion of incorrectly assigned ASVs. All flow cytometry sorted Synechococcus samples could only be amplified by our nested approach. The taxonomic diversity obtained with our primers on both sample types was in agreement with clade distribution observed by previous studies that applied other marker genes or PCR-free metagenomic approaches in similar environmental conditions.

IMPORTANCE

The petB gene has been proposed as a high-resolution marker gene to access the diversity of marine Synechococcus. A systematic metabarcoding approach based on the petB gene would improve the characterization/assessment of Synechococcus community structure in marine planktonic ecosystems. We have designed and tested specific primers to be applied in a nested polymerase chain reaction (PCR) protocol (Ong_2022) for metabarcoding of the petB gene. Ong_2022 protocol can be applied to samples with low DNA content, such as those obtained by flow cytometry cell sorting, granting the simultaneous assessment of the genetic diversity of Synechococcus populations with cellular properties and activities (e.g., nutrient cell ratios or carbon uptake rates). Our approach will allow future studies using flow cytometry to investigate the link between ecological traits and taxonomic diversity of marine Synechococcus.

INTRODUCTION

Marine photosynthetic bacterioplankton is dominated by two genera of picocyanobacteria: Prochlorococcus and Synechococcus. Synechococcus is a polyphyletic genus comprising both marine and freshwater lineages [START_REF] Dvořák | Synechococcus: 3 Billion Years of Global Dominance[END_REF]. They contribute up to 16% of the net marine primary productivity [START_REF] Flombaum | Present and Future Global Distributions of the Marine Cyanobacteria Prochlorococcus and Synechococcus[END_REF][START_REF] Scanlan | Ecological Genomics of Marine Picocyanobacteria[END_REF]. Changes in the physico-chemical properties of marine waters induced by climate change are predicted to have an impact on the distribution of phytoplankton groups, including Synechococcus (4). Given their major role in the marine carbon cycling and ocean food web [START_REF] Li | Primary Production of Prochlorophytes, Cyanobacteria, and Eucaryotic Ultraphytoplankton: Measurements from Flow Cytometric Sorting[END_REF][START_REF] Guidi | Plankton Networks Driving Carbon Export in the Oligotrophic Ocean[END_REF], an extensive assessment of Synechococcus community genetic diversity at spatial and temporal scale as well as their ecological traits are critical to understand the future evolution of this group under the ongoing global change.

Whilst studies using the 16S rRNA gene have revealed some degree of diversity within Prochlorococcus and Synechococcus [START_REF] Urbach | Multiple Evolutionary Origins of Prochlorophytes within the Cyanobacterial Radiation[END_REF][START_REF] Fuller | Clade-Specific 16S Ribosomal DNA Oligonucleotides Reveal the Predominance of a Single Marine Synechococcus Clade throughout a[END_REF][START_REF] Zwirglmaier | Global Phylogeography of Marine Synechococcus and Prochlorococcus Reveals a Distinct Partitioning of Lineages among Oceanic Biomes[END_REF], the relatively low sequence divergence between closely related clades or species and the variability in copy number of rRNA genes [START_REF] Fuller | Clade-Specific 16S Ribosomal DNA Oligonucleotides Reveal the Predominance of a Single Marine Synechococcus Clade throughout a[END_REF][START_REF] Dufresne | Unraveling the Genomic Mosaic of a Ubiquitous Genus of Marine Cyanobacteria[END_REF] have prompted the scientific community to use high-resolution marker genes to analyse the genetic diversity within these groups of picocyanobacteria (e.g. ntcA [START_REF] Penno | Diversity of Synechococcus and Prochlorococcus Populations Determined from DNA Sequences of the N-regulatory Gene ntcA[END_REF], petB [START_REF] Mazard | Multi-Locus Sequence Analysis, Taxonomic Resolution and Biogeography of Marine Synechococcus[END_REF], ITS [START_REF] Rocap | Resolution of Prochlorococcus and Synechococcus Ecotypes by Using 16S-23S Ribosomal DNA Internal Transcribed Spacer Sequences[END_REF] and rpoC1 [START_REF] Palenik | Cyanobacterial Community Structure as Seen from RNA Polymerase Gene Sequence Analysis[END_REF]). Although most of these gene markers have been amplified from natural Synechococcus populations by standard polymerase chain reaction (PCR) amplification, others like rpoC1 requires a nested PCR approach [START_REF] Mühling | High Resolution Genetic Diversity Studies of Marine Synechococcus Isolates Using rpoC1-based Restriction Fragment Length Polymorphism[END_REF][START_REF] Xia | Latitudinal and Vertical Variation of Synechococcus Assemblage Composition Along 170°W Transect From the South Pacific to the Arctic Ocean[END_REF].

Nested and semi-nested PCR approaches involve two sequential amplification reactions with different pair of primers. These amplification methods are often applied to increase the sensitivity and/or specificity of the reaction. They are also particularly useful for samples with low nucleic acid concentrations [START_REF] Mühling | High Resolution Genetic Diversity Studies of Marine Synechococcus Isolates Using rpoC1-based Restriction Fragment Length Polymorphism[END_REF][START_REF] Erlich | Recent Advances in the Polymerase Chain Reaction[END_REF]. For example, metabarcoding of flow cytometry (FCM) sorted phytoplankton populations typically require a nested or semi-nested PCR amplification approach because sorted cells are in general in very low abundance compared to those collected in filtered samples [START_REF] Ribeiro | Small Eukaryotic Phytoplankton Communities in Tropical Waters off Brazil Are Dominated by Symbioses between Haptophyta and Nitrogen-Fixing Cyanobacteria[END_REF][START_REF] Metz | Diversity of Photosynthetic Picoeukaryotes in Eutrophic Shallow Lakes as Assessed by Combining Flow Cytometry Cell-Sorting and High Throughput Sequencing[END_REF].

Metabarcoding is a suitable technique to assess the taxonomic diversity of flow cytometry sorted populations with high sensitivity (e.g. trace concentrations of DNA can be PCR amplified and sequenced) and of many samples, including Synechococcus. Although alternative approaches of metagenome or whole genome sequencing could similarly obtain the taxonomic diversity without metabarcoding limitations of primer and amplification biases, these approaches require a much larger sequencing depth and sample processing time which would increase the cost and time needed [START_REF] Santoferrara | Perspectives from Ten Years of Protist Studies by High-Throughput Metabarcoding[END_REF]. Whole genome sequencing would instead be more advantageous when applied to functional diversity and population genetics [START_REF] Malmstrom | Ecology of Uncultured Prochlorococcus Clades Revealed through Single-Cell Genomics and Biogeographic Analysis[END_REF][START_REF] Santoro | Contributions of Single-Cell Genomics to Our Understanding of Planktonic Marine Archaea[END_REF]. Flow cytometric sorted Synechococcus cells have provided important information about their metabolism. These sorted populations have been used to measure cellular properties such as nutrient stoichiometric and isotopic composition [START_REF] Lomas | Varying Influence of Phytoplankton Biodiversity and Stoichiometric Plasticity on Bulk Particulate Stoichiometry across Ocean Basins[END_REF][START_REF] Hansman | Measuring the in Situ Carbon Isotopic Composition of Distinct Marine Plankton Populations Sorted by Flow Cytometry[END_REF] and cellular activity such as glucose [START_REF] Duhamel | Mixotrophic Metabolism by Natural Communities of Unicellular Cyanobacteria in the Western Tropical South Pacific Ocean[END_REF], phosphorus [START_REF] Casey | Phytoplankton Taxon-Specific Orthophosphate (Pi) and ATP Utilization in the Western Subtropical North Atlantic[END_REF], nitrogen [START_REF] Berthelot | NanoSIMS Single Cell Analyses Reveal the Contrasting Nitrogen Sources for Small Phytoplankton[END_REF][START_REF] Aldunate | Nitrogen Assimilation in Picocyanobacteria Inhabiting the Oxygen-Deficient Waters of the Eastern Tropical North and South Pacific[END_REF] and CO 2 uptake rates [START_REF] Li | Primary Production of Prochlorophytes, Cyanobacteria, and Eucaryotic Ultraphytoplankton: Measurements from Flow Cytometric Sorting[END_REF][START_REF] Duhamel | Mixotrophic Metabolism by Natural Communities of Unicellular Cyanobacteria in the Western Tropical South Pacific Ocean[END_REF][START_REF] Berthelot | NanoSIMS Single Cell Analyses Reveal the Contrasting Nitrogen Sources for Small Phytoplankton[END_REF][START_REF] Rii | Diversity and Productivity of Photosynthetic Picoeukaryotes in Biogeochemically Distinct Regions of the South East Pacific Ocean[END_REF][START_REF] Jardillier | Significant CO2 Fixation by Small Prymnesiophytes in the Subtropical and Tropical Northeast Atlantic Ocean[END_REF][START_REF] Hartmann | Efficient CO2 Fixation by Surface Prochlorococcus in the Atlantic Ocean[END_REF].

However, only few of these studies combined quantitative cell measurements with a fine taxonomic identification of the sorted populations by molecular methods (clone library sequencing and fluorescence in-situ hybridisation only) [START_REF] Rii | Diversity and Productivity of Photosynthetic Picoeukaryotes in Biogeochemically Distinct Regions of the South East Pacific Ocean[END_REF][START_REF] Jardillier | Significant CO2 Fixation by Small Prymnesiophytes in the Subtropical and Tropical Northeast Atlantic Ocean[END_REF][START_REF] Hartmann | Efficient CO2 Fixation by Surface Prochlorococcus in the Atlantic Ocean[END_REF].

The petB gene, encoding the cytochrome b 6 subunit of the cytochrome b 6 f complex, possesses several features that make this locus a good candidate for a gene marker. The petB is a single-copy gene, therefore reducing amplification bias [START_REF] Mazard | Multi-Locus Sequence Analysis, Taxonomic Resolution and Biogeography of Marine Synechococcus[END_REF]. It is highly conserved in length and sequence allowing the obtained reads to be easily aligned [START_REF] Mazard | Multi-Locus Sequence Analysis, Taxonomic Resolution and Biogeography of Marine Synechococcus[END_REF] and has a comprehensive reference nucleotide sequence database that encompass most of the genetic diversity identified within marine Prochlorococcus and Synechococcus [START_REF] Farrant | Delineating Ecologically Significant Taxonomic Units from Global Patterns of Marine Picocyanobacteria[END_REF]. Few specific primers targeting petB gene amplification are currently available in the literature. Mazard et al. (12) developed the original primer set, petB-F and petB-R. Ohnemus et al. [START_REF] Ohnemus | Silicon Content of Individual Cells of Synechococcus from the North Atlantic Ocean[END_REF] modified the original pair by increasing their degeneracy (up to 32 for both forward and reverse primers) in order to expand the diversity coverage of marine picocyanobacteria groups. The annealing region of these primers set is very similar and the degeneracy too high, precluding their use in nested or semi-nested PCR (polymerase chain reaction) amplification approaches.

With the objective of establishing a new nested PCR methodology to survey the diversity of flow cytometry sorted Synechococcus populations, we have designed and tested a new pair of primers (petB-50F and petB-634R) targeting the petB gene locus.

Our approach to assess the diversity of marine Synechococcus populations by metabarcoding of petB includes a nested-PCR amplification (Ong_2022) which combines the newly developed pair (petB-50F and petB-634R) with the original primer set (petB-F and petB-R) established by Mazard et al. (12). We have evaluated the sensitivity and specificity of our nested PCR protocol against the standard PCR protocol with petB-F and petB-R (Mazard_2012) on Synechococcus populations recovered from filtered seawater samples and sorted by flow cytometry from the same water parcels as the filtered samples.

RESULTS

Primer design and PCR amplification. The design of the new primers followed three main criteria. First, the primers must contain none or a very low number of mismatches (≥ 3) across all Synechococcus reference sequences present in the database in order to avoid bias in the amplification. Second, the resulting amplicon sequence must contain enough nucleotide variation to differentiate between the clades and subclades and therefore provide good taxonomic resolution. Last, the amplicon length ideally should be less than 550 bp long, so sequencing could be performed with Illumina 2x300 bp chemistry.

While the first two benchmarks were fulfilled with the new primers set petb-50F and petB-634R, the final amplicon length criteria could not be achieved. In Ong_2022 protocol, the new primers combined with petB-F and petB-R from Mazard et al. [START_REF] Mazard | Multi-Locus Sequence Analysis, Taxonomic Resolution and Biogeography of Marine Synechococcus[END_REF] (petB-634R during the first round and petB-50F during the second round of amplification) produced an amplicon of 571 bp length (from 50-618 position of the petB gene Synechococcus sp. WH8109 reference, accession number CP00688) while the standard Mazard_2012 PCR generates a fragment of 597bp (Table 1). There was no overlap between the forward and reverse reads obtained with both standard and nested PCR.

To solve this problem, during sequence processing, we have merged our reads by adding 10 degenerate base 'N at their 3 prime end. Amplicons without overlap have been successfully processed with DADA2 pipeline using this bioinformatic manoeuvre, which does not significantly influence k-mer based taxonomic classification such as the naive Bayesian classifier implemented in DADA2 [START_REF] Needham | Pronounced Daily Succession of Phytoplankton, Archaea and Bacteria Following a Spring Bloom[END_REF][START_REF] Jeraldo | IM-TORNADO: A Tool for Comparison of 16S Reads from Paired-End Libraries[END_REF]. Moreover, concatenated sequences (i.e. joined without overlap) for longer amplicons were shown to improve taxonomic recovery and classification over merged sequences [START_REF] Dacey | Concatenation of Paired-End Reads Improves Taxonomic Classification of Amplicons for Profiling Microbial Communities[END_REF].

A total of 71 DNA samples extracted from filtered seawater and 18 flow cytometry sorted samples were obtained from subtropical (ST1 and ST2) and subantarctic waters (SA1, SA2 and SA3) along the Chatham rise (Figure S1, Table S1). Among the DNA extracted samples, an amplification product was obtained for all environmental samples using Ong_2022 nested PCR and one sample could not be amplified with the standard Mazard_2012. Hence, from this point, the comparison between Mazard_2012 and Ong_2022 approaches will only include the 70 samples successfully amplified by both methods (Table S1). The results from sorted samples were only obtained with Ong_2022 PCR since none of the sorted samples were amplified with success by standard PCR reaction.

Following sequence processing, the total number of reads retained in the final dataset was about 49% and 23% of the initial number of reads generated with Mazard_2012 and Ong_2022 respectively (Table 2). The chimera detection step removed 66% of reads obtained with Ong_2022 and 35% for Mazard_2012 after merging forward and reverse reads (Table 2). The resulting mean number of reads obtained with the standard PCR was almost 3 times higher than the mean number of reads obtained with the nested PCR. However, the number of amplicon sequence variants (ASVs) generated with the nested PCR was higher (6307) when compared with Mazard_2012 standard reaction (4716) (Table 2).

Community composition from filtered samples. Four major Synechococcus subclades that were present in almost all samples were Ia, Ib, IVa and IVb (Figure 1). While these subclades represented nearly 99% of the total number of reads obtained by both approaches, the percentage of reads per sample for each of the major subclades differed between the amplification methods: Mazard_2012 produced a higher proportion of reads from subclades IVa and IVb while Ong_2022 amplified a higher proportion of subclade Ia and Ib sequences (Figure 2A). Yet, patterns of community composition such as the higher relative abundance of clade IVa in ST compared to SA waters were captured by both approaches. Among the seven minor clades and subclades present in both amplifications (CRD1, EnvA, EnvB, UC-A, Ic, II-WPC2, V), the proportion of reads of clades CRD1, EnvA, EnvB and V were higher in Ong_2022 compared to Mazard_2012 (Figure 2B, Table S2). Six clades and subclades with low contribution to the total community were only amplified by either Ong_2022 (WPC1, IIh, IIe, VIb) or Mazard_2012 (IIb, VIc) PCR method (Figure 2B). Paired samples Wilcoxon test indicated that the median percent of reads for each clade and subclade in Mazard_2012 was significantly different from Ong_2022, except for clade V present at very low proportion (Table S2). The mean percent of reads for subclades Ia and Ib with Ong_2022 were higher than with Mazard_2012 (8.51% and 46.73% vs. 3.47% and 31.04%, respectively). The mean percent of reads for subclades IVa and IVb in Ong_2022 were lower compared to Mazard_2012 (11.96% and 31.40% vs. 20.39% and 44.04%, respectively) (Table S2). Among clades and subclades with minor contribution of reads to the total community, the mean percentage of reads of clade UC-A and subclades Ic and II-WPC2 were higher in Mazard_2012 compared to Ong_2022 (Figure 2B, Table S2). Clade UC-A was detected in all samples across both amplification types, but the average percentage of reads was higher in Mazard_2012 (0.63%) compared to Ong_2022 (0.13%). Subclades II-WPC2 and Ic were present in higher proportion and detected in more samples with Mazard_2012 compared to Ong_2022. For subclade II-WPC2, the average percentage of reads in Mazard_2012 was 0.26% and found in 67 samples, compared to Ong_2022, 0.01% in 27 samples. The average percentage of reads of subclade Ic in Mazard_2012 was 0.12% and found in 51 samples, compared to Ong_2022 of 0.03% and found in 20 samples. The average percentage of reads of clades CRD1, EnvA, EnvB, and V were higher in Ong_2022 (0.93%, 0.18%, 0.11% and 0.0006%, respectively) compared to Mazard _2012 (0.02%, 0.001%, 0.03% and 0.0003%, respectively) (Figure 2B, Table S2). Clades CRD1, EnvA and V were present in more samples in Ong_2022 (35, 21 and 3 samples, respectively) compared to Mazard_2012 (24, 5 and 1 samples, respectively). Clade EnvB was detected in the same 23 samples for both amplification types. Clade WPC1 and subclades IIh, IIe and VIb were present only in Ong_2022, although they were found in low proportion (less than 0.21% of reads in a sample) and only in 1 to 4 samples each. Subclades VIc and IIb were present only in Mazard_2012, but in very low proportion (less than 0.10 % of reads in a sample) and only in 1 to 2 samples (Table S2).

Most clades and subclades had a higher number of ASVs with Ong_2022 compared to Mazard_2012 (Figure 2C). Subclades with a higher mean percent of reads also had a higher number of ASVs. The exceptions were clade UC-A, and subclades Ic and II-Ong et al.

WPC2 which had a higher number of ASVs relative to their percent of reads and a higher number of ASVs in Mazard_2012 (Table S2).

The relative abundance of the major subclades followed a similar pattern across amplification methods and cycles, with notable exception for subclades IVa and IVb in ST1 and SA2 cycles, respectively. Subclade IVa was the dominant in ST1 samples with the Mazard_2012 approach, while the same samples processed with Ong_2022 showed similar contribution of subclades IVa and Ib to total community (Figure 1, Table S3). In cycles ST2 and SA1, the proportions of subclade IVa decreases and subclade IVb increases in both methods. In SA2 and SA3, subclade IVb becomes dominant only in Mazard_2012 samples (Figure 1, Table S3). Clades CRD1, EnvA and EnvB were found in higher proportion in subantarctic oceanic cycles SA2 and SA3 compared to more coastal SA1 and oceanic subtropical cycles for both amplifications (Figure 1, Table S3).

The mean percent of reads for each subclade in cycles ST1, ST2 and SA1 was less than 0.004% and 0.06% for Mazard_2012 and Ong_2022, respectively. The mean percent of reads increased in cycles SA2 and SA3 to 0.06% and 2.8% for CRD1, 0.003% and 0.58% for EnvA, and 0.09% and 0.35% for EnvB, for Mazard_2012 and Ong_2022, respectively.

ASV sequences similarity to reference database. Within each clade and subclade, we compared the nucleotide similarity of the ASVs obtained with Mazard_2012 and Ong_2022 and with the reference sequences in the database (32) (Figure 3 and Table S4). In general, similarity between the ASVs within each clade and subclade obtained with both PCR approaches was lower when compared with the reference database, showing that both primers sets were able to capture a higher diversity than the one available in the database (Figure 3). The only exceptions were clades CRD1, EnvA and EnvB, and subclade VIc (Figure 3 and Table S4). Amidst the clades and subclades amplified by both PCR approaches, similarity within the ASVs obtained with Ong_2022 (from 85.9% to 94.7%) were lower than with Mazard_2012 (from 92.3% to 99.8%) for all the subclades with two exceptions, clade UC-A (87.7% and 30.5% respectively) and subclade II-WPC2 (85.9% and 25.9%, respectively) (Figure 3 and Table S4), indicating that Ong_2022 captured a higher microdiversity than Mazard_2012.

To confirm the taxonomic assignment of the ASVs, the nucleotide sequences of ASVs and references were translated into protein sequences. For each clade and sub-clade, we compared the translated ASVs from each PCR approaches (Mazard_2012 and Ong_2022) against the translated reference sequences. ASVs for which the protein sequence had 95% pairwise identity against translated reference sequences were regarded as correctly assigned. Overall, Ong_2022 had a higher proportion of total ASVs that were correctly assigned compared to Mazard_2012 (98.2% vs. 84.2%; Table S4).

The percentage of correctly assigned ASVs in each clade and subclade was comparable between Mazard_2012 (from 96.2% to 100%) and Ong_2022 (from 90.5% to 100%) except for clade UC-A (no correctly assigned ASVs) and subclade II-WPC2 (2.4%) obtained with Mazard_2012 (Table S4). This result shows that several ASVs in Mazard_2012 assigned to UC-A and II-WPC2 did not represent petB gene sequences, hence the very low percentage of nucleotide similarity among the ASVs for these lineages as mentioned previously. These ASVs assigned to UC-A and II-WPC2 from Mazard_2012 amplification were further examined by BLAST against GenBank nucleotide database. More than half of these ASVs had no sequence match in the database, indicating that they could have resulted from spurious amplification or the current database did not have a corresponding reference sequence.

Community composition from flow cytometry sorted Synechococcus popula-

tions. We tested our nested PCR with Synechococcus cells sorted by flow cytometry from 18 water samples collected in five cycles (Table S1). The 998 ASVs obtained from sorted cells with Ong_2022 were assigned to fewer clades and subclades belonging from subcluster (SC) 5.1 compared to ASVs obtained from filtered samples (Figure 4,

S5).

Similar to the results obtained with filtered samples, Ia, Ib, IVa and IVb were the major subclades present in all the samples, representing an average of 99.98% of total reads (Figure 4 and Table S5). Clades UC-A, CRD1 and EnvB and subclade II-WPC2 had low contribution to the total community (Figure 4 and Table S5). The mean percent of reads for subclades Ia and Ib (3.84% and 31.33% respectively) from sorted cells were lower than those obtained from filtered samples (8.51% and 46.73%, Table S2 andS5).

In contrast, for subclades IVa and IVb, the mean percent of reads were higher (25.47% and 39.35% respectively) in sorted cells than filtered samples (11.96% and 31.40%, Table S2 andS5).

Among minor clades and subclades, clade UC-A was found in 11 out of 18 sorted samples with an average proportion reads of 0.01% in contrast with filtered samples where this clade represented 0.13% of reads and was present in all samples. Clades CRD1 and EnvB and subclade II-WPC2 were present in 1 or 2 sorted samples, while in filtered samples they were present in one-third to half of the samples. Clades CRD1 and EnvB, and subclade II-WPC2 had a lower average percent of reads (0.003%, 0.002% and 0.002%) compared to Ong_2022 (0.93%, 0.11% and 0.01%). The 7 minor clades and subclades (EnvA, WPC1, V, IIh, IIe, Ic and VIb) that were found in filtered samples were not present in sorted samples.

The average number of ASVs per sample obtained from sorted Synechococcus cells was almost 2.5 times lower compared to filtered samples (82 ± 35 ASVs and 217 ± 63 ASVs, respectively). Each clade and subclade also had a lower number of ASVs.

The relative abundance of major subclades followed a similar pattern between sorted Synechococcus and filtered seawater samples, with the exceptions of subclades IVa and IVb (Figure 1 and4). Subclade IVa and subclade IVb were dominant in ST and SA sorted Synechococcus populations, respectively. This pattern was also observed with Mazard_2012 approach on filtered samples but differ of Ong_2022 where the overall contribution of subclades IVa and IVb was similar to that of subclades Ia and Ib (Figure 1).

DISCUSSION

We have established a nested PCR protocol to assess the diversity of marine Synechococcus based on metabarcoding analysis of the petB gene marker. Our protocol combines a pair of newly developed primers (petB-50F and petB-634R) and the original primer set (petB-F and petB-R) established by Mazard et al. [START_REF] Mazard | Multi-Locus Sequence Analysis, Taxonomic Resolution and Biogeography of Marine Synechococcus[END_REF]. We have compared the sensitivity and specificity of the new nested amplification protocol (Ong_2022) against the original standard PCR (Mazard_2012) on DNA extracted from filtered seawater samples and determined the performance of the proposed nested amplification protocol on flow cytometry sorted Synechococcus populations. Whilst both standard and nested PCR worked in nearly all filtered seawater samples, sorted Synechococcus samples could not be amplified with Mazard_2012 standard PCR likely due to low ini-tial concentration of Synechococcus DNA. Nested PCR is often employed when DNA concentration is too low to be amplified by a standard PCR [START_REF] Yu | Nested PCR Biases in Interpreting Microbial Community Structure in 16S rRNA Gene Sequence Datasets[END_REF][START_REF] Haro | Evaluation of Established Methods for DNA Extraction and Primer Pairs Targeting 16S rRNA Gene for Bacterial Microbiota Profiling of Olive Xylem Sap[END_REF], thus Ong_2022 approach could be an alternative method when the initial concentration of Synechococcus cells is expected to be low, especially in sorted samples.

None of the ASVs obtained by both methods on filtered samples have been assigned to Prochlorococcus, although Prochlorococcus cells were detected through flow cytometry (from 0.3 x 10 3 to 90.6 x 10 3 cells ml -1 ) in subantarctic water samples (SA1, SA2 and SA3) (data not shown). Similar to our new primers, Mazard et al. [START_REF] Mazard | Multi-Locus Sequence Analysis, Taxonomic Resolution and Biogeography of Marine Synechococcus[END_REF] petB-F and petB-R are also biased toward Synechococcus [START_REF] Mazard | Multi-Locus Sequence Analysis, Taxonomic Resolution and Biogeography of Marine Synechococcus[END_REF]. However, during their study, Prochlorococcus petB sequences were recovered at low frequency from samples where Prochorococcus cell abundances were at least 45 times higher compared to Synechococcus. In our study, Synechococcus cells were always in equivalent or higher abundance In addition, Synechococcus cells from SC 5.2 and 5.3 were probably absent in our samples. SC 5.2 representatives are mainly found in estuaries and river-influenced coastal waters [START_REF] Jing | Phylogenetic Diversity of Marine Synechococcus Spp. in the Sea of Okhotsk[END_REF][START_REF] Xia | Comparison of the Seasonal Variations of Synechococcus Assemblage Structures in Estuarine Waters and Coastal Waters of Hong Kong[END_REF][START_REF] Xia | Synechococcus Assemblages across the Salinity Gradient in a Salt Wedge Estuary[END_REF][START_REF] Doré | Global Phylogeography of Marine Synechococcus in Coastal Areas Reveals Strikingly Different Communities than in the Open Ocean[END_REF] and although SC 5.3 contain both freshwater and marine representatives [START_REF] Farrant | Delineating Ecologically Significant Taxonomic Units from Global Patterns of Marine Picocyanobacteria[END_REF][START_REF] Cabello-Yeves | Ecological and Genomic Features of Two Widespread Freshwater Picocyanobacteria[END_REF], the latter have been sporadically detected at high abundance only in some regions such as the Red Sea, Mediterranean Sea and northwestern Atlantic Ocean [START_REF] Farrant | Delineating Ecologically Significant Taxonomic Units from Global Patterns of Marine Picocyanobacteria[END_REF][START_REF] Huang | Novel Lineages of Prochlorococcus and Synechococcus in the Global Oceans[END_REF][START_REF] Sohm | Co-Occurring Synechococcus Ecotypes Occupy Four Major Oceanic Regimes Defined by Temperature, Macronutrients and Iron[END_REF].

One of the main concerns with nested PCR is the generation of chimeric sequences that are not true members of the community [START_REF] Yu | Nested PCR Biases in Interpreting Microbial Community Structure in 16S rRNA Gene Sequence Datasets[END_REF]. Indeed, based on chimera removal analysis, our nested PCR approach had twice the proportion of chimeras detected than the standard amplification. However, Ong_2022 had an overall higher proportion of correctly assigned ASVs compared to Mazard_2012, indicating that the nested PCR was more specific toward Synechococcus. While the use of higher number of cycles might increase the chimera generation and sensitivity (e.g., Ong_2022 generally increased the proportion of reads of low abundant clades and subclades) in nested approaches, Ong et al.

the improved specificity of Ong_2022 PCR is derived from the binding of two separate sets of oligonucleotides to the same target template [START_REF] Green | Nested Polymerase Chain Reaction (PCR)[END_REF]. Additionally, Ong_2022 PCR captured a higher genetic diversity of Synechococcus SC 5.1 compared to Mazard_2012. Within common subclades amplified by both Mazard_2012 and Ong_2022, almost all subclades had a higher number of ASVs in the latter approach. These ASVs had lower nucleotide similarity within each subclades indicating that our approach was able to capture a higher microdiversity within Synechococcus subclades. This directly contrasts with previous studies, which reported that nested PCR approaches resulted in a lower number of lineages and operational taxonomic units compared to a standard PCR [START_REF] Yu | Nested PCR Biases in Interpreting Microbial Community Structure in 16S rRNA Gene Sequence Datasets[END_REF][START_REF] Haro | Evaluation of Established Methods for DNA Extraction and Primer Pairs Targeting 16S rRNA Gene for Bacterial Microbiota Profiling of Olive Xylem Sap[END_REF].

The number and proportion of reads of minor clades and subclades decreased in sorted samples (11 vs. 4 clades and subclades, 1% vs. 0.2% of reads for filtered and sorted samples, respectively). Metz et al. [START_REF] Metz | Diversity of Photosynthetic Picoeukaryotes in Eutrophic Shallow Lakes as Assessed by Combining Flow Cytometry Cell-Sorting and High Throughput Sequencing[END_REF] investigated the diversity of photosynthetic picoeukaryotes in lakes using a similar approach. They found that the richness of filtered samples was on average lower than in the sorted samples. One major difference between Metz et al. [START_REF] Metz | Diversity of Photosynthetic Picoeukaryotes in Eutrophic Shallow Lakes as Assessed by Combining Flow Cytometry Cell-Sorting and High Throughput Sequencing[END_REF] sorting protocol and Ong_2022 is the number of cells sorted. In Metz et al. [START_REF] Metz | Diversity of Photosynthetic Picoeukaryotes in Eutrophic Shallow Lakes as Assessed by Combining Flow Cytometry Cell-Sorting and High Throughput Sequencing[END_REF], the number of cells sorted was 22 to 50 times higher than in Ong_2022. Since microbial plankton communities contain a large number of taxa that are present in low abundance [START_REF] Sogin | Microbial Diversity in the Deep Sea and the Underexplored "Rare Biosphere[END_REF], a larger sample size would probably increase the number of taxa detected and improve the representation of rare taxa with Ong_2022 approach on sorted cells [START_REF] Cermeño | Sampling the Limits of Species Richness in Marine Phytoplankton Communities[END_REF]. However, the trade-off between flow cytometry sorting time and the detection of rare microbial groups needs to be considered when designing the study. Cell sorting by flow cytometry is a time-consuming process, and increasing the number of cells sorted for numerous natural marine samples might not be feasible.

The taxonomic diversity obtained with our primers agreed with the clade distribution observed in previous studies conducted with other marker genes and PCR free methods in geographical regions with similar environmental conditions. For example, subclades Ia, Ib, IVa and IVb were the major lineages observed in the datasets obtained from filtered samples (∼99%) and sorted samples (99.8%) with nested and standard PCR. This agrees with previous reports, in which clades I and IV most often co-occur in the open ocean at latitudes above 30°N and below 30°S [START_REF] Zwirglmaier | Global Phylogeography of Marine Synechococcus and Prochlorococcus Reveals a Distinct Partitioning of Lineages among Oceanic Biomes[END_REF][START_REF] Mazard | Multi-Locus Sequence Analysis, Taxonomic Resolution and Biogeography of Marine Synechococcus[END_REF][START_REF] Farrant | Delineating Ecologically Significant Taxonomic Units from Global Patterns of Marine Picocyanobacteria[END_REF][START_REF] Xia | Latitudinal and Vertical Variation of Synechococcus Assemblage Composition Along 170°W Transect From the South Pacific to the Arctic Ocean[END_REF][START_REF] Sohm | Co-Occurring Synechococcus Ecotypes Occupy Four Major Oceanic Regimes Defined by Temperature, Macronutrients and Iron[END_REF][START_REF] Kent | Parallel Phylogeography of Prochlorococcus and Synechococcus[END_REF]. However, the relative abundance of clades I and IV varied between filtered samples amplified with Mazard_2012 and Ong_2022 as well as in sorted samples. Filtered samples amplified with Mazard_2012 and sorted samples amplified with Ong_2022 followed a similar pattern, whereby clade IV dominated in subtropical and subantarctic waters.

In contrast, the proportion of clades I and IV remained similar across sampling locations when filtered samples were amplified with Ong_2022. Prior studies conducted in surrounding areas report different proportions of clade I to IV. Both clades were detected in similar proportions along the coast of New Zealand at overlapping latitudes (45°S) by PCR independent method [START_REF] Doré | Global Phylogeography of Marine Synechococcus in Coastal Areas Reveals Strikingly Different Communities than in the Open Ocean[END_REF]. Clade IV has been observed in higher proportion compared to clade I in the western and eastern South Pacific Ocean at slightly lower latitudes with warmer conditions (around 30°S compared to 42°S in this study) [START_REF] Zwirglmaier | Global Phylogeography of Marine Synechococcus and Prochlorococcus Reveals a Distinct Partitioning of Lineages among Oceanic Biomes[END_REF][START_REF] Farrant | Delineating Ecologically Significant Taxonomic Units from Global Patterns of Marine Picocyanobacteria[END_REF][START_REF] Sohm | Co-Occurring Synechococcus Ecotypes Occupy Four Major Oceanic Regimes Defined by Temperature, Macronutrients and Iron[END_REF]. Hence, we cannot ascertain which method is more closely reflecting the natural abundance of the dominant clades.

We observed that relative abundance of clades CRD1, EnvA and EnvB increased in high nutrient low chlorophyll (HNLC) subantarctic waters where iron is typically the proximate limiting nutrient (filtered samples amplified with both approaches) [START_REF] Boyd | Role of Iron, Light, and Silicate in Controlling Algal Biomass in Subantarctic Waters SE of New Zealand[END_REF][START_REF] Longhurst | Ecological Geography of the Sea[END_REF].

Clade CRD1 was reported to co-occur with clades I and IV in the southeastern Pacific Ocean and other cold and high latitude HNLC conditions [START_REF] Farrant | Delineating Ecologically Significant Taxonomic Units from Global Patterns of Marine Picocyanobacteria[END_REF][START_REF] Xia | Latitudinal and Vertical Variation of Synechococcus Assemblage Composition Along 170°W Transect From the South Pacific to the Arctic Ocean[END_REF][START_REF] Sohm | Co-Occurring Synechococcus Ecotypes Occupy Four Major Oceanic Regimes Defined by Temperature, Macronutrients and Iron[END_REF]. More recent reports suggested that EnvA and EnvB also correspond to lineages found in low iron conditions [START_REF] Farrant | Delineating Ecologically Significant Taxonomic Units from Global Patterns of Marine Picocyanobacteria[END_REF][START_REF] Sohm | Co-Occurring Synechococcus Ecotypes Occupy Four Major Oceanic Regimes Defined by Temperature, Macronutrients and Iron[END_REF][START_REF] Caputi | Community-Level Responses to Iron Availability in Open Ocean Plankton Ecosystems[END_REF][START_REF] Ahlgren | Genomic Mosaicism Underlies the Adaptation of Marine Synechococcus Ecotypes to Distinct Oceanic Iron Niches[END_REF]. The Ong_2022 approach recorded subclades IIh and IIe in cycle ST1 at very low abundance. Although members from clade II are mostly warm thermal types that dominate in nutrient depleted tropical areas (3), subclades IIe and IIh were also reported at low relative abundances in cooler waters (14.1-17.5 °C) where iron was not limiting for phytoplankton growth [START_REF] Farrant | Delineating Ecologically Significant Taxonomic Units from Global Patterns of Marine Picocyanobacteria[END_REF]. Surface waters at ST1 were at a similar temperature at the time of sampling (13 °C), which is likely within the temperature range where subclades IIe and IIh were detected. Subclades WPC1, VIb, VIc and V were detected in low concentrations, but their precise ecological niche and biogeographical distribution have not yet been established [START_REF] Farrant | Delineating Ecologically Significant Taxonomic Units from Global Patterns of Marine Picocyanobacteria[END_REF].

In summary, we have designed a nested PCR protocol with two new primer sequences to amplify the petB gene from marine Synechococcus. Results obtained from Ong et al.

filtered seawater samples suggest that the new nested PCR approach is specific of Synechococcus and captured a wider genetic diversity, especially for rare groups. The nested PCR protocol successfully amplified Synechococcus cells from flow cytometry sorted samples and recovered a composition for the dominant subclades similar to that of filtered samples. Our amplification protocol will improve our understanding Synechococcus communities by allowing the determination of Synechococcus diversity directly associated with quantitative measurements (e.g. nutrient cell ratios or carbon update rates) obtained on sorted samples.

MATERIALS AND METHODS

Sampling Location. Samples were collected during the Salp Particle expOrt and Oceanic Production (SalpPOOP) cruise (TAN1810, 21st October -21st November 2018) conducted near the Chatham Rise (Aotearoa -New Zealand) on board the R/V Tangaroa. The Chatham Rise is a dynamic region where northward-moving subtropical (ST) water masses mix with southward moving subantarctic (SA) water masses to form the Subtropical Convergence Zone [START_REF] Heath | Observations on Chatham Rise Currents[END_REF][START_REF] Chiswell | Acoustic Doppler Current Profiler Measurements over the Chatham Rise[END_REF][START_REF] Stramma | The South Pacific Current[END_REF][START_REF] Lüskow | Biology of Salpa Thompsoni at the Chatham Rise, New Zealand: Demography, Growth, and Diel Vertical Migration[END_REF]. ST waters are warm, saline and macronutrient-depleted while SA waters are cool, less saline and high-nitrate, low chlorophyll, low-silicate (HNLC-LSi) waters where low iron and silicate are the primary limiting factors for phytoplankton growth and productivity [START_REF] Boyd | Role of Iron, Light, and Silicate in Controlling Algal Biomass in Subantarctic Waters SE of New Zealand[END_REF][START_REF] Hall | Structure and Dynamics of the Pelagic Microbial Food Web of the Subtropical Convergence Region East of New Zealand[END_REF]. We adopted Lagrangian sampling strategy using a satellite-tracked 15 m drogue drifted array that allowed us to track and sample the same water parcel multiple times during 3 to 6 days (hereafter referred to as 'cycles') [START_REF] Landry | Lagrangian Studies of Phytoplankton Growth and Grazing Relationships in a Coastal Upwelling Ecosystem off Southern California[END_REF]. Five cycles were conducted, two in subtropical waters (ST1 and ST2) and three in subantarctic waters (SA1, SA2 and SA3) (Figure S1). In each cycle, water column profiles of temperature, salinity, dissolved oxygen, photosynthetically active radiation (PAR) and fluorescence were measured 4-6 times daily by Seabird (SBE 911plus) Conductivity -Temperature -Depth (CTD) with a PAR sensor. The cycles were identified as subtropical or subantarctic based on physical measurements during sampling following Lüskow et al. [START_REF] Lüskow | Biology of Salpa Thompsoni at the Chatham Rise, New Zealand: Demography, Growth, and Diel Vertical Migration[END_REF] definition. Samples used for picocyanobacterial community analysis presented in this study were collected using a CTD rosette sampler equipped with 10 L Niskin bottles from six depths (5-70 m). For DNA analysis of filtered seawater samples, 1.5-2 L seawater from all depths was filtered through a 0.22 µm pore-size Sterivex filter using a peristaltic pump (n samples = 71) (Table S1). For DNA analysis of flow cytometry sorted samples, 1.5 mL samples were collected from the surface mixed layer (SUR, 12 m) and the deep chlorophyll maximum (DCM, 40 m or 70 m) (Table S1), where they were preserved with a solution of DMSO and 1/100 pluronic acid mix at 10% final concentration, and flash-frozen in liquid nitrogen (n samples = 21). All samples were stored at -80 °C until laboratory processing. 61) with default parameters and used to locate the biding region attending the criteria described. The selected forward primer petB-50F is located at positions 50-66 with reference to the first nucleotide of the petB gene from Synechococcus sp. WH8109 (accession number CP00688) (Table 1). The selected reverse primer petB-634R is located at positions 615-634 with reference to the first nucleotide of the petB gene from Synechococcus sp. WH8109 (accession number CP00688) (Table 1). Primer petB-50F had no mismatches in the binding region for Synechococcus SC 5.1 reference sequences, 2 to 3 mismatches for SC 5.2 and 5.3, and excluded

Prochlorococcus (≥ 3 mismatches) and other marine cyanobacteria (e.g., Richelia sp.).

Primer petB-634R had mainly no mismatches in the binding region for Synechococcus S1). As the community was similar at each depth, only samples that had 10,000 cells sorted were used for analysis and for filtered samples, all depths were combined. Analysis was performed with R Phyloseq package Version 1.36.0 (66) and Geneious Prime Version 2019.2.3 [START_REF] Kearse | Geneious Basic: An Integrated[END_REF]. All sequence processing and analysis scripts along with the reformatted database can be found in https: //github.com/deniseong/marine-Synechococcus-metaB.

FIGURES

List of Figures

( 1 .

 1 2 to 45 times) than Procholorococcus which might explain the absence of ASVs assigned to the latter.ASVs assigned to Synechococcus SC 5.2 and 5.3 were also not detected. Mazard et al.[START_REF] Mazard | Multi-Locus Sequence Analysis, Taxonomic Resolution and Biogeography of Marine Synechococcus[END_REF] (petB-F and petB-R) and our primers (petB-50F and petB-634R) are biased toward SC 5.1, which is represented in the database by the majority of the sequences[START_REF] Farrant | Delineating Ecologically Significant Taxonomic Units from Global Patterns of Marine Picocyanobacteria[END_REF].

Flow cytometry cell sorting.

  Synechococcus cells were sorted using FACSAria™ flow cytometer (Becton Dickinson, San Jose, CA), equipped with a laser emitting at 488 nm, 70 mm nozzle and the following filters: 488/10 band pass (BP) for size side scatter (SSC), 576/26 BP for orange fluorescence, and 655 long pass for red fluorescence. Synechococcus populations were detected by their signature in plots of chlorophyll-a red fluorescence versus size scatter (SSC) and discriminated from photosynthetic picoeukaryotes based on the orange fluorescence of phycoerythrin. Tris-HCl 50 mM, pH 8.0 and, NaCl 10 mM buffer filter-sterilized (0.22 µm) was used as sheath liquid. Sheathpressure was set at 70 PSI and nozzle frequency was 90,000 Hz with a deflection voltage of 6,000 V. Cells were sorted in purity mode and collected into Eppendorf tubes containing Tris-EDTA lysis buffer (Tris 10 mM, EDTA 1 mM, and 1.2% Triton, final concentration) as described previously[START_REF] Ribeiro | Small Eukaryotic Phytoplankton Communities in Tropical Waters off Brazil Are Dominated by Symbioses between Haptophyta and Nitrogen-Fixing Cyanobacteria[END_REF]. Sheath fluid samples collected during sorting were processed and analysed as negative controls in all subsequent steps, including sequencing, to test for contamination in the flow sorting process. petB-50F and petB-634R primer design. The criteria for primer design were: (i) 3 or fewer mismatches in binding region, (ii) no nucleotide ambiguity or mismatches at the 5 and 3 prime ends, (iii) low degeneracy, (iv) binding region between 15 and 20 base pairs, (v) annealing temperature difference between the forward and reverse primers not greater than 4 °C, (vi) primer binding region should be absent of insertions and deletions among the different groups and (vii) no hairpins and self-dimerisation within primers sequence. Four hundred and seven Synechococcus nucleotide sequences of the petB gene from Farrant et al. (32) reference database were aligned using Multiple Alignment and Fast Fourier Transform (MAFFT) alignment Version 7.450 (60) in Geneious Prime Version 2019.2.3 (

SC 5 .

 5 1 to 5.3, and excluded Prochlorococcus (≥ 3 mismatches). After the binding region was selected, the sequences were truncated to the length amplified by the new forward primer sequence and petB-R from Mazard et al.[START_REF] Mazard | Multi-Locus Sequence Analysis, Taxonomic Resolution and Biogeography of Marine Synechococcus[END_REF] and used to infer a maximum likelihood phylogenetic tree with FastTree Version 2.1.11 (62) in Geneious Prime with default parameters. The phylogenetic tree was used to determine if the shorter petB amplicon could distinguish between subclades. All alignments and phylogenetic trees can be found in https://github.com/deniseong/marine-Synechococcus-metaB.DNA extraction and PCR. DNA from Sterivex filters was extracted using the Nucleospin Plant II DNA extraction kit from Macherey-Nagel with modification of manufacturer instructions. The modified protocol can be found on https://github.com/ deniseong/marine-Synechococcus-metaB. DNA from flow cytometry sorted Synechococcus cells was obtained using three cycles of flash-freezing and thawing in liquid nitrogen[START_REF] Ribeiro | Small Eukaryotic Phytoplankton Communities in Tropical Waters off Brazil Are Dominated by Symbioses between Haptophyta and Nitrogen-Fixing Cyanobacteria[END_REF]. Mazard_2012 protocol refers to PCR amplification with petB-F (5'-TACGACTGGTTCCAGGAACG-3') and petB-R (5'-GAAGTGCATGAGCATGAA-3') primers[START_REF] Mazard | Multi-Locus Sequence Analysis, Taxonomic Resolution and Biogeography of Marine Synechococcus[END_REF], performed on extracted DNA from filtered seawater. PCR reactions were done using the following reaction: extracted DNA template with an average final concentration of 1.3 ng µL -1 (from 0.007 to 4.11 ng µL -1 ), 12.5 µL KAPA HiFi HotStart ReadyMix 2X, 0.3 µM primers (petB-F and petB-R) final concentration, 0.1 to 0.4 µg µL -1 final concentration of bovine serum albumin (BSA) and H 2 O to 25 µL final volume reaction. Thermal conditions were: 94 °C for 5 min, followed by 30 cycles of 94 °C for 30 s, 55 °C for 30 s, 72 °C for 45 s, and a final cycle of 72 °C for 6 min. All PCRs were performed in duplicates and pooled together. Ong_2022 protocol refers to a nested PCR performed on extracted DNA from either filtered seawater or flow cytometry sorted Synechococcus cells. The first round of PCR amplification was done using the following reaction: extracted DNA template with an average final concentration of 0.54 ng µL -1 (from 0.015 to 2.055 ng µL -1 ) or volume corresponding to approximately 160-400 sorted Synechococcus cells (4 µL maximum), 5 µL KAPA HiFi HotStart ReadyMix 2×, 0.3 µM final concentration of primer petB-F, 0.3 µM final concentration of primer petB-634R (5'-GCTTVCGRATCATCARGAAG-3'), 0.1 µg µL -1 final concentration of BSA (only for filtered samples) and H 2 O for a 10 µL reaction. Thermal conditions were: 94 °C for 5 min, followed by 30 cycles of 94 °C for 30 s, 59 °C for 30 s, 72 °C for 45 s, and a final cycle of 72 °C for 6 min. For the second round of amplification the following conditions were used: 2.5 µL of first round product, 12.5 µL KAPA HiFi HotStart ReadyMix 2X, 0.3 µM final concentration of primer petB-50F (5'-TYCAGGACATYGCTGAY-3'), 0.3 µM final concentration of primer petB-R and H 2 O for a 25 µL reaction. Thermal conditions were: 94 °C for 5 min, followed by 30 cycles of 94 °C for 30 s, 55 °C for 30 s, 72 °C for 45 s, and a final cycle of 72 °C for 6 min. All first and second round PCRS were performed in duplicates and pooled together. Samples were purified, barcoded, and sequenced by the GeT-PlaGe platform of GenoToul (INRAE Auzeville, France) using an Illumina MiSeq platform (2 x 300 bp). petB rRNA gene amplicon analysis. Sequences were processed on RStudio Version 1.4.1717 (63). Primer sequences were removed using Cutadapt Version 3.4 (64). Fastq files were trimmed and quality filtered using function 'filterandtrim' with the DADA2 R package Version 1.12 (65). Based on sequence quality score, reads obtained with Mazard_2012 and Ong_2022 were trimmed and filtered using the following options: truncLen = c(280, 280) and c(250, 280), minLen = c(280, 280) and minLen = c(250, 280) respectively with truncQ = 2, maxEE = c(10, 10) in both datasets. Forward and reverse sequences were dereplicated and grouped. To merge the forward and reverse reads, the 'mergePairs' function on DADA2 was used with the justConcatenate = TRUE option to add 10 degenerate base 'N' between forward and reverse reads. Chimeras were identified and removed using the function 'removeBimeraDenovo'. Taxonomy was obtained with the function 'assignTaxonomy'. ASVs were assigned against petB reference database from Farrant et al. (32) reformatted for use with DADA2. The following taxonomic nomenclature was maintained: three major Ong et al. Synechococcus/Cyanobium lineages called subclusters (SC) 5.1 through 5.3; fourteen clades (I to IX, XVI, XX, UC-A, WPC1, CRD1, EnvA and EnvB) within SC 5.1 and 20 subclades (Ia to Ic, IIa to IIh and II-WPC2, IIIa to IIIc, IVa and IVb, and VIa to VIc). ASVs assigned as Richelia sp. were removed prior to normalising reads. At station SA1_5, sorted samples from two depths (12 m and 40 m) were sorted to two and three different total number of cells (Table

FIG 1 FIG 2 FIG 3 FIG 4

 1234 FIG 1 Synechococcus taxonomic composition at clade and subclade level from filtered samples obtained with Mazard_2012 and Ong_2022 approaches. Samples were grouped by cycles and ordered across a spatial gradient, from subtropical (ST) to subantarctic (SA) cycles. Samples from all depths were combined.

Figure 1 27 Figure 2 28 Figure 3 29 Figure 4

 1272283294 Figure 1 Synechococcus taxonomic composition at clade and subclade level from filtered samples obtained with Mazard_2012 and Ong_2022 approaches. Samples were grouped by cycles and ordered across a spatial gradient, from subtropical (ST) to subantarctic (SA) cycles. Samples from all depths were combined. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

FIG S1 FIG S2

 S1S2 FIG S1 Study area with cycle locations at the Chatham Rise east of Aotearoa-New Zealand. Black dots represent the sampling location within each cycle. Bathymetric lines (m) are indicated. Bathymetry was obtained from Amante & Eakins (67) and plotted with ggOceanMaps (68). Sea surface temperature data (°C) was obtained from MODIS (NASA) and averaged over November 2018. ST and SA refers to subtropical and subantarctic respectively and numbers to the cycle

TABLE 1

 1 Primers sequence and PCR conditions used for amplification. Based on the Synechococcus sp. WH 8109 (CP006882), the petB gene is located between position 440,919 and 441,575 of the genome. The total length of the petB gene is 657bp. The amplicon length includes the primers. The position of the primer in petB gene is indicated in reference to the first nucleotide of the gene. Deg = Degeneracy, T(°C) = Annealing temperature PCR type Round Amplicon length Primer name Direction Sequence Gene position GC% Deg T (°C) Citation

	Mazard_2012 Standard 1	597	petB-F	Forward TACGACTGGTTCCAGGAACG 22-41	55	0	55	Mazard et al. (12)
					petB-R	Reverse GAAGTGCATGAGCATGAA	601-618	44	0		Mazard et al. (12)
	Ong_2022	Nested	1	613	petB-F	Forward TACGACTGGTTCCAGGAACG 22-41	55	0	59	Mazard et al. (12)
					petB-634R	Reverse GCTTVCGRATCATCARGA	615-634	47	12		This study
			2	569	petB-50F	Forward CAGGACATYGCTGAY	50-66	50	8	55	This study
					petB-R	Reverse GAAGTGCATGAGCATGAA	601-618	44	0		Mazard et al. (12)

Table 2

 2 Number of reads and ASVs per sample (Mean ± SD) after each sequence processing step using DADA2 on filtered samples obtained with Mazard_2012 and Ong_2022 approaches. . . . . . . . . . . . . . . 26 Table S1 Location and date of samples. . . . . . . . . . . . . . . . . . . . . . . . . 33

	Table S2	Number of ASVs and percentage of reads at clade and subclade level
		for 70 distinct samples: mean, maximum, minimum and standard de-
		viation, and the total number of ASVs using Mazard_2012 and Ong_2022.
		Significant differences of the percent of reads per clade or subclade
		between the two approaches were tested using a paired sample Wilcoxon
		test (P-value < 0.05). The rows are arranged in descending order ac-
		cording to the mean percentage of reads from Ong_2022. . . . . . . .	35
	Table S3	The mean percentage of reads for each clade or subclade in each cycle,
		using Mazard_2012 and Ong_2022 amplification protocol. The rows
		are arranged in descending order according to the mean percentage
		of reads of each clade or subclade from Ong_2022. . . . . . . . . . . .	36
	Table S4	Percent similarity of nucleotide sequences/ASVs and percent of cor-
		rectly assigned ASVs in each clade or subclade from Mazard_2012 and
		Ong_2022 approaches on filtered samples and reference database. .	37
	Table S5	Percentage of reads at clade and subclade level across 18 samples:
		mean, maximum, minimum and standard deviation, and the total num-
		ber of ASVs of sorted Synechococcus cells. The rows are arranged in
		descending order according to the mean percentage of reads. . . . .	37

TABLE S1 (

 S1 continued)

							Filtered seawater	Sorted Synechococcus
	Date	Latitude Longitude cycle station CTD	Depth (m) Sample number Mazard 2012	Ong 2022	Sample name No. of cells
	15/11/2018 -43.73 -179.82	ST2 2	U9159 20	172	CTD-petB-172	CTD-NpetB-172
	15/11/2018 -43.73 -179.82	ST2 2	U9159 30	173	CTD-petB-173	CTD-NpetB-173
	15/11/2018 -43.73 -179.82	ST2 2	U9159 40	174	CTD-petB-174	CTD-NpetB-174
	15/11/2018 -43.73 -179.82	ST2 2	U9159 50	175	CTD-petB-175	CTD-NpetB-175
	13/11/2018 -43.48 179.94	ST2 3	U9152 12		syn-69	2,000
	13/11/2018 -43.48 179.94	ST2 3	U9152 40		syn-74	2,000

TABLE S2

 S2 Number of ASVs and percentage of reads at clade and subclade level for 70 distinct samples: mean, maximum, minimum and standard deviation, and the total number of ASVs using Mazard_2012 and Ong_2022. Significant differences of the percent of reads per clade or subclade between the two approaches were tested using a paired sample Wilcoxon test (P-value < 0.05). The rows are arranged in descending order according to the mean percentage of reads from Ong_2022.

			Mazard_2012		Ong_2022
				% of reads		% of reads
	Clade Subclade No. of ASVs Mean	Max Min SD	No. of ASVs Mean	Max Min	SD	P-value
	1 I	Ib	899	31.04	53.07 10.16 13.46	1444	46.73	56.32 29.30	4.33	4.07E-12
	2 IV	IVb	1543	44.04	84.40 0.33 28.00	2023	31.40	44.39 0.35	14.99	2.33E-06
	3 IV	IVa	934	20.39	83.48 0.00 27.60	1526	11.96	39.21 0.01	13.45	6.17E-03
	4 I	Ia	249	3.47	8.52 1.32 1.42	518	8.51	17.03 0.00	4.26	3.30E-12
	5 CRD1		11	0.02	0.18 0.00 0.03	92	0.93	5.08 0.00	1.40	5.28E-11
	6 EnvA		5	9.78E-04 0.02 0.00 4.13E-03 30	0.18	1.48 0.00	0.37	1.67E-06
	7 UC-A		662	0.63	1.81 0.01 0.52	565	0.13	0.50 9.98E-03 0.11	1.72E-09
	8 EnvB		8	0.03	0.21 0.00 0.05	15	0.11	1.12 0.00	0.21	2.38E-07
	9 I	Ic	26	0.12	0.79 0.00 0.15	21	0.03	0.26 0.00	0.05	2.57E-09
	10 II	II-WPC2 374	0.26	2.26 0.00 0.38	44	0.01	0.32 0.00	0.04	1.20E-11
	11 II	IIh	0			11	0.01	0.26 0.00	0.05
	12 II	IIe	0			3	4.30E-03 0.25 0.00	0.03
	13 WPC1		0			8	3.01E-03 0.21 0.00	0.03
	14 VI	VIb	0			4	1.00E-03 0.03 0.00	4.86E-03
	15 V		2	3.26E-04 0.02 0.00 2.73E-03 3	5.71E-04 0.02 0.00	2.89E-03 0.87
	16 VI	VIc	2	2.72E-04 0.01 0.00 1.63E-03 0
	17 II	IIb	1	1.47E-03 0.10 0.00 0.01	0

TABLE S4

 S4 Percent similarity of nucleotide sequences/ASVs and percent of correctly assigned ASVs in each clade or subclade from Mazard_2012 and Ong_2022 approaches on filtered samples and reference database.

			No. of	% similarity of nucleotide seq	% of ASVs correctly assigned
	Clade Subclade ref seq Ref seq Mazard_2012 Ong_2022 Mazard_2012 Ong_2022
	I	Ib	55	94.8	92.3	91.9	98.0	99.6
	IV	IVb	20	97.6	94.9	91.9	99.3	96.7
	IV	IVa	21	95.9	93.9	91.9	99.5	98.4
	I	Ia	19	95.6	93.2	91.6	99.2	98.5
	CRD1		23	90.8	94.1	91.6	100	93.5
	EnvA		10	90.9	95.7	94.7	100	100
	UC-A		1		30.5	87.7	0	98.1
	EnvB		13	91.3	96.9	92.7	100	100
	I	Ic	18	94.7	95.0	92.3	96.2	90.5
	II	II-WPC2	1		25.9	85.9	2.4	97.7
	II	IIh	10	92.9		90.8		100
	II	IIe	4	98.1		93.2		100
	WPC1		11	92.8		89.7		100
	VI	VIb	20	97.6		92.5		100
	V		2	99.7	92.7	87.4	100	100
	VI	VIc	4	97.1	99.8		100
	II	IIb	13	93.3	N/A		100

TABLE S5

 S5 Percentage of reads at clade and subclade level across 18 samples: mean, maximum, minimum and standard deviation, and the total number of ASVs of sorted Synechococcus cells. The rows are arranged in descending order according to the mean percentage of reads.

				% of reads
	Clade Subclade No.of ASVs	Mean Max Min	SD
	1 IV	IVb	413	39.35 78.82 0.00 30.39
	2 I	Ib	274	31.33 46.38 12.84 9.82
	3 IV	IVa	244	25.47 78.76 0.00 30.04
	4 I	Ia	37	3.84 16.84 0.00 5.55
	5 UC-A		24 8.11E-03 0.05 0.00 0.01
	6 CRD1		2 3.34E-03 0.05 0.00 0.01
	7 II	II-WPC2	3 2.15E-03 0.04 0.00 0.01
	8 EnvB		1 2.15E-03 0.04 0.00 0.01

± 63

0 0 0 2000

ACKNOWLEDGMENTS

We acknowledge the crew of RV Tangaroa for their efforts in facilitating the sampling throughout the TAN1810 voyage. We thank Florian Humily and Sophie Mazard for providing the primer design of petB-634R. We thank Daniel Vaulot (CNRS) for feedback on data processing and editing comments. We are grateful to Jaret Bilewitch and Debbie Hulston (National Institute of Water and Atmospheric Research, New Zealand) for their laboratory work on the DNA extraction of filtered seawater samples and Debbie Hulston for coordinating the logistics of transporting the DNA extracted samples from New Zealand to Singapore. We thank Priscillia Gourvil (CNRS) for coordinating the logistics of transporting samples from France to Singapore. We are grateful to Karen Selph (University of Hawaii at Manoa, School of Ocean and Earth Science & Technology) for her comments on the abundance of Prochlorococcus cell obtained by flow cytometry in the samples collected during TAN1810. We thank Molly A. Moynihan (Marine Biological Laboratory) for laboratory assistance and guidance, as well as coordinating the logistics of Illumina sequencing. We thank Clarence Sim Wei Hung (NTU) for laboratory assistance. ALS and DO were supported by the Singapore Ministry of Education, Academic Research Fund Tier 1 (RG26/19). LG was supported by the French "Agence Nationale de la Recherche" Programs CINNAMON (ANR-17-CE02-0014-01). AGR was supported by NIWA via the New Zealand Ministry of Business, Innovation and Employment's Strategic Science Investment Funding to the National Coasts Ocean Centre.

Availability of data and materials. Raw sequencing data were deposited in NCBI Sequence Read Archive under the bioProject number PRJNA885274. Source code, DNA extraction protocol as well as all metadata generated are available on GitHub (https: //github.com/deniseong/marine-Synechococcus-metaB). Raw data files from flow cytometry sorting is available on http://flowrepository.org/experiments/1773 (Repository ID: FR-FCM-Z5P8). Step-by-step protocol for DNA extraction and PCR amplification are published on protocols.io at DOI: 10.17504/protocols.io.5qpvorkj7v4o/v1.

Author contributions statement DRYO and ALS conceived the study. AGR organized the SalpPOOP voyage and sampling structure. AGR and ALS developed field protocol and processed the samples during the voyage. ALS and DM performed flow cytometry cell sorting. DRYO developed primers and processed the samples for sequencing. DRYO analysed data and drafted the manuscript. All authors were involved in data interpretation and critically revised the manuscript. Table 1 Primers sequence and PCR conditions used for amplification. Based

List of Tables

SUPPLEMENTAL MATERIAL