Vivien Maisonneuve
email: vivien.maisonneuve@cri.mines-paristech.fr

Olivier Hermant
email: olivier.hermant@cri.mines-paristech.fr

François Irigoin
email: francois.irigoin@cri.mines-paristech.fr

Preservation of Lyapunov-Theoretic Proofs: From Real to Floating-Point Arithmetic

In a paper, Feron presents how Lyapunovtheoretic proofs of stability can be migrated toward computer-readable and verifiable certificates of control software behavior by relying of Floyd's and Hoare's proof system. However, Lyapunov-theoretic proofs are addressed towards exact, real arithmetic and do not accurately represent the behavior of realistic programs run with machine arithmetic. We address the issue of preserving those proofs in presence of rounding errors resulting from the use of floating-point arithmetic: we present an automatic tool, based on a theoretical framework the soundness of which is proved in Coq, that translates Feron's proof invariants on real arithmetic to similar invariants on floating-point numbers, and preserves the proof structure. We show how our methodology allows to verify whether stability invariants still hold for the concrete implementation of the controller.

We study in details the application of our tool to the open-loop system of Feron's paper and show that stability is preserved out of the box. We also translate Feron's proof for the closed-loop system, and discuss the conditions under which the system remains stable.

I. Introduction

Stability constitutes an essential attribute of control systems, especially when human safety is involved, as in medical or aeronautical domains. Modern system developments, such as adaptive control technologies, rely on robust stability and performance criteria as the primary justification for their relevance to safety-critical control applications. Motivated by such applications, there exist many theorems that ensure system stability and performance under various assumptions and in various settings. Lyapunov's stability theory plays a critical role in that regard.

The low-level software implementation of a control law can be inspected by analysis tools available to support the development of safety-critical computer programs. The simplest program analysis technique consists of performing several simulations, sometimes including a software or hardware representation of the controlled system in the loop. However, simulations provide information about only a finite number of system behaviors. More advanced methods include model checking and abstract interpretation, e.g. using Astrée [START_REF] Cousot | The astrée static analyzer[END_REF]. In these methods, inputs are the very computer programs and outputs are certificates of proper program behavior along the chosen criterion. Another possibility is to use theorem-proving techniques, supported by tools such as Coq, Isabelle or PVS [START_REF]The coq proof assistant[END_REF], [START_REF] Paulson | [END_REF], [START_REF] Owre | PVS specification and verification system[END_REF]. These proof assistants can be used to establish properties of programs and more general mathematical constructs. Model checking, abstract interpretation, and theorem-proving tools are all used to verify safety-critical applications.

In [START_REF] Feron | From control systems to control software[END_REF], Feron investigates how control-system domain knowledge and, in particular, Lyapunov-theoretic proofs of stability for the high-level theoretical modelization, can be migrated towards computer-readable and verifiable certificates of control software behavior by relying on Floyd's and Hoare's proof system [START_REF] Peled | Software reliability methods[END_REF], applied to MATLAB pseudocode (see Sections II and III). Feron presents both openloop and closed-loop implementations of his controller, i.e. depending whether the system action on the environment and its feedback are modelized or not. But errors resulting from the use of floating-point arithmetic are not considered.

In this paper, we address this issue by presenting an automatic approach to translate Lyapunov-theoretic stability proof invariants on pseudo-code with real arithmetic, as provided in Feron's article, to similar invariants on machine code that take into account rounding errors introduced by floating-point arithmetic. We use them to verify whether the stability conditions still hold, in which case system stability with floating-point numbers is established. This approach is based on a generic theoretical framework which soundness is proved in Coq, described further in this document, and is implemented as an automatic tool. This document is organized as follows. First, we describe the second-order dynamical system example used by Feron in [START_REF] Feron | From control systems to control software[END_REF], with the corresponding controller. Next, Feron's analysis of the open-loop controller with real numbers is presented. In the next sections, we introduce our generic translation scheme and its implementation as a Python library, and use it thereafter to rewrite Feron's analysis with floating-point arithmetic. Finally, the case of the closed-loop system is handled. The document concludes with a discussion on the generality of this approach.

II. Motivating Example by Feron

We consider the first system described in the article of Feron [START_REF] Feron | From control systems to control software[END_REF]. It is a dynamical system composed of a single mass and a single spring shown in Figure 1.

The position input y of the mass is available for feedback control. The signal y d is the reference signal, that is, the desired position to be followed by the mass. A discrete-time MATLAB implementation of the controller, using real numbers, is provided in [START_REF] Feron | From control systems to control software[END_REF]. The source code is shown below. The process to produce this code from the system modeling is covered in detail in Feron's paper, but in not necessary to understand this document. In this code, the skip statement is a null operation: when it is executed, nothing happens. It is useful as a placeholder when a statement is required syntactically to be surrounded by invariants, but no code needs to be executed.

Apart from the mechanical system state observation y and the desired system output y d , variables in this code are:

• x c = xc1 xc2 ∈ R 2 is the discrete-time controller state; • y c ∈ [-1, 1]
is the bounded output tracking error, i.e. the input (yy d ) is passed through a saturation function to avoid variable overflow in the controller;

• u ∈ R is the mechanical system input, i.e. the action to be performed according to the controller.

Constants A c , B c , C c and D c are the discrete-time controller state, input, output and feedthrough matrices. The commands send and receive are basically I/O: they respectively send and receive data given in the commands first argument through a specific channel given by the commands second argument.

III. Open-Loop Stability Proof with Reals

The stability proof of this system relies on Lyapunov theory. In simple terms, a system is Lyapunov stable if all states x c reachable from an initial starting state belonging to a bounded neighborhood V of an equilibrium point x e remain in V .

Lyapunov theory provides constraints that must be satisfied by such a V . On linear systems, they are equations that can be solved using linear matrix inequalities [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]. Commonly, V is an ellipsoid.

In this case, to prove Lyapunov stability, we need to show that at any time, x c belongs to the set E P chosen by Feron according to Lyapunov's theory:

E P = {x ∈ R 2 | x T • P • x ≤ 1}, P = 10 -3 0.6742 0.0428
0.0428 2.4651 using E P as the stability neighborhood V . This set is the full ellipse shown in Figure 2, centered around 0 and slightly slanted: A stability proof of the controller is provided in [START_REF] Feron | From control systems to control software[END_REF], using Floyd-Hoare program annotation technique [START_REF] Peled | Software reliability methods[END_REF]

x c ∈ E P ⇐⇒ 0.6742x c 2 1 +0.0856x c1 x c2 +2.4651x c 2 2 ≤ 1000.
skip; % ( xc yc ) ∈ E Qµ , Q µ = µP 02×1 01×2 1-µ , µ = 0.9991 10 u = Cc*xc + Dc*yc; % ( xc yc ) ∈ E Qµ , u 2 ≤ (C c D c ) • Q -1 µ • (C c D c ) -1 11 xc = Ac*xc + Bc*yc; % x c ∈ E R , R = (A c B c ) • Q -1 µ • (A c B c ) T -1 , % u 2 ≤ (C c D c ) • Q -1 µ • (C c D c ) -1 12 send(u, 1); % x c ∈ E R 13 receive(y, 2); % x c ∈ E R receive(yd, 3); % x c ∈ E R 15 skip; % x c ∈ E P 16 end
Most of the proof relies on algebraic arguments. For example, the invariant loosening on Line 9:

% x c ∈ E P , y 2 c ≤ 1 9 skip; % ( xc yc ) ∈ E Qµ , Q µ = µP 02×1 01×2 1-µ , µ = 0.9991 means x c ∈ E P ∧ y 2 c ≤ 1 =⇒ ( xc yc ) ∈ E Qµ , with Q µ = µP 02×1 01×2 1-µ and µ = 0.9991.
The correctness of this assertion stems from the fact that, given any value of µ ∈ [0, 1], the domain E Qµ is a solid ellipsoid, centered around 0, and whose intersection with the planes y c = 1 and y c = -1 is equal to E P . Consequently, the solid bounded cylinder 3).

C = {( xc yc ) | x c ∈ E P ∧ y 2 c ≤ 1} is included within E Qµ (Figure
E Qµ E P C x c2 y c 1 -1 x c1 Figure 3. Inclusion of E P within E Qµ
The next two invariants

% ( xc yc ) ∈ E Qµ , Q µ = µP 02×1 01×2 1-µ , µ = 0.9991 10 u = Cc*xc + Dc*yc; % ( xc yc ) ∈ E Qµ , u 2 ≤ (C c D c ) • Q -1 µ • (C c D c ) -1 11 xc = Ac*xc + Bc*yc; % x c ∈ E R , R = (A c B c ) • Q -1 µ • (A c B c ) T -1 , % u 2 ≤ (C c D c ) • Q -1 µ • (C c D c
) -1 also rely on similar algebraic arguments and theorems. Other invariants are trivial. Finally, only the very last loosening

% x c ∈ E R 13 skip; % x c ∈ E P i.
e. E R ⊂ E P , that closes the loop, is not purely algebraic since its validity relies on the numerical parameters A c , B c , C c , D c . This assertion needs to be checked to ensure the correctness of the proof statements. This can be done either numerically, or algebraically for at most two-dimensional systems like this one.

We have checked this proof using Mathematica [8]. Our proof notebook is available online 1 .

IV. Proof Scheme with Limited-Precision Numbers

We would like to check that the stability proof still holds on a realistic controller device, using limitedprecision numbers (for instance, but not necessarily, floating-point numbers). When using limited-precision arithmetic, the values of constants are slightly altered and calculations are likely to produce rounding errors.

In absolute terms, it is impossible to switch from real to limited-precision numbers without affecting the behavior of the controller. Thus, the stability proof cannot be preserved in the general case. On the other hand, the program still can be stable if rounding errors are small enough and the final inclusion E R ⊂ E P leaves them enough room. We study how proof invariants can be tweaked so that they apply to a limited-precision semantic. Our goal is to derive from the proof scheme for real numbers a proof scheme suited for limited-precision arithmetic, whose correctness, although not guaranteed, can be checked as easily as in the original proof.

A. Formalism

In this section, an abstract formalism is defined to handle a program an its proof, with both real and limitedprecision arithmetic. Let X be the (finite) set of variables in the program.

% d i % d ′ = p(d, i)

1) Values:

Let F ⊂ R be the targeted set of limitedprecision representations of real numbers (e.g., floatingpoint numbers). We assume that each real number c ∈ R has a finite-precision representation c ∈ F and that the property (1) is satisfied:

c ∈ F =⇒ c = c. (1)
Throughout the rest of the document, we will use the uniform notation K to refer to either the set R or F.

2) Functions Symbols:

We consider a set of function symbols F R on real numbers. Each real function symbol f ∈ F R is associated to a finite-precision counterpart f ∈ F F . Functions symbols in F R and F F will be simply referred as functions when there will be no ambiguity.

Each function symbol f in F R or F F has an arity n ∈ N and can be evaluated on its domain:

f ∈ F K =⇒ eval(f ) ∈ K n → K. (2)
Functions f and f have the same arity: arity f = arity f .

As for now, no other assumption is made as to the behavior of f compared to f .

Notice that function symbols are not identified to their evaluations. The reason is that functions with similar evaluations on real arithmetic may behave differently on a limited-precision paradigm, e.g. function x → 2 100 + x -2 100 compared to the identity function: x → x.

3) Valuations:

A valuation v on K is a function that maps variables in X to values in K. The set of valuations on K is noted V K . V K ∋ v : X → K. Notice that V F is a subset of V R . Each valuation v in V R is associated to a valuation v in V F , defined as follows: v : x ∈ X -→ v(x). (3)
Using this definition, it can easily be shown using ( 1) that for any valuation

v ∈ V R , v ∈ V F =⇒ v = v. (4)
Let v 1 and v 2 be two valuations in V R . The distance dist between v 1 and v 2 is the valuation that maps a variable in X to the distance in absolute value between its associated values in v 1 and v 2 :

V R ∋ dist(v 1 , v 2 ) = x ∈ X -→ |v 1 (x) -v 2 (x)|. (5)
We say that v 1 is lower than v 2 , and note

v 1 ≤ v 2 , if ∀x ∈ X, v 1 (x) ≤ v 2 (x). (6)
Finally, the sum of the two valuations, noted v 1 + v 2 , is the valuation defined as follows:

v 1 + v 2 : x ∈ X -→ v 1 (x) + v 2 (x). (7
d = {v | v ∈ d}. (8)
As for valuations in F, it can be shown that

d ∈ D F =⇒ d = d. (9)
We also define an operation to overapproximate a domain with respect to a valuation. The extension of domain

d ∈ D R using a positive valuation v ∈ V R , noted d ⊕ v, is defined as follows: d ⊕ v = {v 2 ∈ V R | v 1 ∈ d ∧ dist(v 1 , v 2 ) ≤ v}. ( 10
)
5) Expressions: We consider a very simple language with variables, constants and function calls. The expressions E K of this language are constructed as follows:

E K ∋ e ::= x ∈ X variable | c ∈ K constant | f (e 1 , . . . , e n ) function call with f ∈ F K , arity f = n and ∀i ∈ [1, n], e i ∈ E K
An expression e ∈ K can be evaluated in the environment described by the valuation v ∈ V K . The result is a value in K.

eval(e, v) =        v(x) if e = x ∈ X c if e = c ∈ K eval(f )(c 1 , . . . , c n ) if e = f (e 1 , . . . , e n ) with ∀i ∈ [1, n], c i = eval(e i , v)
where the notation eval is overloaded to handle expressions.

6) Instructions:

An instruction is the affectation of an expression value to a variable. The set of instructions is noted I K :

I K = X × E K We note x := e the instruction (x, e) ∈ I K .
As for expressions, an instruction

(x := e) ∈ I K can be evaluated in an environment v ∈ V K . The result is a valuation in V K . eval(x := e, v) = y ∈ X -→ eval(e, v) if x = y v(y) if x = y
This definition can be extended to evaluate the image of a domain d ∈ D K , the result being also a domain of D K .

eval(x := e, d) = {eval(x := e, v) | v ∈ d} 7) Invariant Propagation: An invariant propagator, or simply propagator, is a function p that takes as input a domain in D K , known as precondition, an instruction in I K , and returns a domain in D K called postcondition.

p : D K × I K -→ D K
A propagator is expected to return a valid postcondition for a given precondition and instruction. Formally, this can be described by the correctness condition:

∀(d, i) ∈ D K × I K , eval(i, d) ⊂ p(d, i). ( 11
)
The set of (correct) propagators is noted P K .

Propagators are an abstract representation of the arguments used to propagate Floyd's and Hoare's invariant en route to prove a program. Unlike these arguments, propagators have to be defined on the whole domain D K × I K , that is, for any precondition and instruction. To make propagators total functions, the universal postcondition can always be returned if the domain is outside the scope of the corresponding argument, for instance when the hypothesis of some Lyapunov theorem is not verified.

B. Translation Scheme

The general translation scheme that is implemented is the following (see Figure 5). Each instruction i in the abstract controller code with real numbers is turned into an corresponding instruction ī on the concrete controller, where real constants and operators are replaced by their limited-precision counterparts. A post-condition d′ on the limited-precision side is computed using the same proof argument p as in the original code, enlarged to include the error term bound v err that may occur in limited-precision arithmetics in ī. It is computed on an intermediate instruction form ĩ which meaning will be explained later. The precondition d might have been previously replaced by d as the postcondition of a preceding instruction. The error term bound v err depends on the operators of the instruction ī and on the information on program variables given by d. The limited-precision arithmetic error is not necessarily bounded even when the variables are, e.g. considering division. In this case, we are not able to compute a suitable d′ and the proof translation attempt fails.

% d i % d ′ = p(d, i) % d ī % d′ ⊃ p( d, ĩ) ⊕ v err
Also, d′ will be used as a precondition in the next instruction. An effort has to be made so that d′ has a similar shape to d ′ , in order to fit with the proof argument p ′ used in this instruction, if possible. Thus, d′ is an overapproximation of p( d, ī) ⊕ v err in the domain of p ′ . In consequence, the automated analysis falls short if such a d′ cannot be found.

C. Translation Steps

In this section, we consider an instruction i = (x := e) in I R , along with a precondition d and a postcondition d ′ in D R . d, i and d ′ are linked through a propagator p in P R :

d ′ = p(d, i). ( 12
)
As outlined in Section IV-B, we study how the instruction i can be translated into an equivalent instruction ī ∈ I F on limited-precision arithmetic, while using information on the precondition d and propagator p of i to compute an interesting postcondition d′ for ī (given a modification d of d). This is a two-step process:

1) First, the program constants c ∈ R are converted into their limited-precision representations c ∈ F while keeping absolute-precision functions, giving an instruction ĩ ∈ I R ; 2) Second, the functions f ∈ F R that appear in the code are changed into their limited-precision counterparts f ∈ F F .

The resulting code relies on limited-precision arithmetic only, which is what is being sought. The translation scheme is described in Figure 6.

% d i % d ′ = p(d, i) % d ĩ % d′ = p( d, ĩ) % d ī % d′ ⊃ p( d, ĩ) ⊕ v err

1) Converting Constants:

We construct an expression ẽ ∈ E R , obtained from e by converting its constants from R to F. Formally, ẽ is defined recursively as follows:

ẽ =    x if e = x ∈ X c if e = c ∈ R f (ẽ 1 , . . . , ẽn ) if e = f (e 1 , . . . , e n ) (13)
We let ĩ = (x := ẽ).

Let d ∈ D F be the precondition corresponding to d after converting both constants and functions in the preceding instructions of the program. Let d′ = p( d, ĩ) ∈ D R . According to Equation 11, d′ is a valid postcondition for precondition d and instruction ĩ.

2) Converting Functions: The next step is to convert real functions f that appear in the expression into their limited-precision counterparts f , as defined in Section IV-A2. The resulting expression ē ∈ E F is obtained with:

ē =    x if e = x ∈ X c if e = c ∈ R f (ē 1 , . . . , ēn ) if e = f (e 1 , . . . , e n ) (14)
Compared to the definition of ẽ (Equation 13), we just have turned functions symbols f into f . As previously, we also let ī = (x := ē).

Unlike ĩ, instruction ī relies on different functions than i (taken from F F instead of F R ): the underlying invariant propagation argument in propagator p is very likely not to be applicable to it. In this case, p( d, ī) would fall back on the universal domain, which of course is valid but prevents us from doing any interesting proof further in the program.

Functional Arithmetic Errors:

To circumvent this issue, we propose when it is possible to enlarge the invariant d′ = p( d, ĩ) found in Section IV-C1 to take into account arithmetic errors introduced when converting functions from F R to F F . Such errors are called functional arithmetic errors, but we will refer to them as arithmetic errors, or simply errors, later in this document.

We define the arithmetic error on an expression ẽ ∈ E R and a valuation v ∈ V F . The result is a value in R.

err(ẽ, v) = | eval(ē, v) -eval(ẽ, v)| ( 15
)
This definition is extended to define the arithmetic error of an instruction ĩ = (x := ẽ) in I R on a valuation v ∈ V F . The result is a valuation in V R defined as follows:

err(x := ẽ, v) = y ∈ X -→ err(ẽ, v) if x = y 0 if x = y .
∀ d′ ∈ D F , d′ ⊕ v err ⊂ d′ =⇒ eval(ī, d) ⊂ d′ (16)
holds. In other words, if d′ is a postcondition for an instruction ĩ on real numbers, on a given precondition domain d, and d′ is a superset of d′ loose enough to support the arithmetic error v err , then d′ is a valid postcondition for ī, the limited-precision counterpart of instruction ĩ.

Proving this relation raises no particular difficulty, but is a bit too long to fit in this document. A proof in Coq, limited to binary functions (which include all the functions that we need in this paper), was implemented and is available online.

V. Open-Loop Stability Proof Scheme with Floating-Point Numbers

We would like to check whether the stability proof in Feron's program still holds on a controller implemented with floating-point numbers. In this section, we use the IEEE Standard for Floating-Point Arithmetic [START_REF][END_REF] encoded on 64 bits2 , as most of today's floating-point units do. In this standard, both addition and multiplication are correctly rounded depending on the active rounding mode, which allows to bound the rounding error depending on the values of operands. The case of alternative numeric representations is discussed in Section VII.

A. Automatic Translation

We have developed a program to automatically perform these computations, following the translation scheme framework described in the previous section. It is a module called LyaFloat, written in Python and built upon Python libraries SymPy (to handle symbolic mathematics) and Mpmath (for arbitrary-precision floating-point arithmetic).

Here is the listing of a script that automatically computes the floating-point output invariant ellipsoid E R (in variable ERbar) corresponding to E R in the original proof, and tests whether E R ⊂ E P . Two cases are possible:

• either E R ⊂ E P , then the program is Lyapunovstable on a floating-point architecture;

• or E R ⊂ E P : as E R was obtained through overapproximations, we cannot conclude about the program behavior.

In this open-loop case, our program LyaFloat is able to check the inclusion. Thus, the stability of the open-loop system with a 64-bit IEEE 754 compliant implementation is formally proven to hold using our proof translation scheme.

Notice that the precision of the floating-point arithmetic can be set in setfloatify. This allows us to check that the controller is still stable on a 32-bit only architecture, or to compute that the minimum required precision is 17 bits. This can be useful if, instead of adapting the controller code and/or its proof to the hardware, we have to select a controller device that suits the proof.

In the remainder of this section, we describe the operations performed by our tool to compute E R . It follows the two-step scheme described in Section IV-C. Alternatively to our program, an implementation of this proof using Mathematica is available online: see Footnote 1.

u = Cc*xc + Dc*yc; % ( xc yc ) ∈ E Qµ , u 2 ≤ (C c D c ) • Q -1 µ • (C c D c ) -1
First of all, matrices C c and D c must be replaced by their floating-point counterparts C c and D c both in the program instruction Line 10 and the ensuing invariant. This invariant relies only on algebraic arguments and does not depend on the values in the matrices, it still holds considering exact arithmetic operations. But this is not sufficient: indeed, this instruction is a sum of matrix multiplications, i.e. a set of additions and multiplications on floating-point numbers that yield rounding errors.

We can notice that entering this instruction, the values of matrices C c , D c and E Qµ are known, and the values of x c and y c are bounded by the precondition 17), we deduce: Here we are able to find algebraic solutions, but this may be impossible with ellipsoids of higher dimension. Still, we would be able to find bounds using numerical methods.

( xc yc ) ∈ E Qµ , that is 0.000673593x c 2 1 + 0.000085523x c1 x c2 + 0.00246288x c 2 2 + 0.9991y 2 c ≤ 1. (17) From (
         |x c1 | ≤ 3 • 10
In floating-point arithmetic, rounding errors created by addition and multiplication operators can be bounded when the operands are known or bounded by (18), provided that overflow, underflow, and denormalized numbers do not occur [START_REF] Goldberg | What every computer scientist should know about floating-point arithmetic[END_REF], [11].

Here, we need to compute

C c x c + D c y c = C c(0,0) x c1 + D c y c
where all values in the right-hand term are known or bounded. Thus, a constant ε can be computed that bounds the absolute rounding error created when computing u. We obtain:

ε = 5.90 • 10 -12
This way, starting from the algebraic result obtained on real numbers

|u| ≤ (C c D c ) • Q -1 µ • (C c D c )
T we can ensure that with floating-point numbers, the inequality holds:

|u| ≤ Ū = (C c D c ) • Q -1 µ • (C c D c ) T + ε
which leads to the invariants:

% ( xc yc ) ∈ E Qµ , Q µ = µP 02×1 01×2 1-µ , µ = 0.9991 10 u = Cc*xc + Dc*yc; % ( xc yc ) ∈ E Qµ , u 2 ≤ Ū 2
2) Invariant on x c : The next instruction, considering constant changes, is:

11 xc = Ac*xc + Bc*yc; % x c ∈ E R , R = (A c B c ) • Q -1 µ • (A c B c ) T -1 ,
where R is defined the same way R is, using floatingpoint terms A c , B c instead of the real-valued counterparts

A c , B c
, and E R is the ellipsoid built upon R. Again, this invariant holds independently of matrices values.

Here, we compute the values affected to

x c = xc 1 xc 2 : A c(0,0) x c1 + A c(0,1) x c2 + y c A c(1,0) x c1 + x c2
Using the same method as above, absolute rounding errors introduced by floating-point operations can be bounded on dimensions x c1 and x c2 by constants ε 1 = 7.42 • 10 -15 and ε 2 = 3.62 • 10 -15 .

These constants must be taken into account in the postcondition. Then the postcondition can be replaced by

% x c ∈ E R
where E R is an ellipse that includes E R plus the rounding error terms (see Figure 7). As mentioned in Section IV-B, replacing the ellipse in the postcondition by another ellipse has the advantage of introducing little change in the stability proof sketch (instead of using a different domain, which would involve using different theorems), which can greatly facilitate tweaking the rest of the proof in longer codes. Formally, E R must satisfy:

∀x c ∈ E R , ∀x c ′ ∈ R 2 , |x c ′ 1 -x c1 | ≤ ε 1 ∧ |x c ′ 2 -x c2 | ≤ ε 2 =⇒ x c ′ ∈ E R (19) E R E R x c1 x c2 ε 1 ε 2 Figure 7.
Relation between E R and E R At the end of the proof scheme, the system is stable with floating-point numbers if and only if the inclusion E R ⊂ E P holds. To succeed, E R should be as narrow as possible with respect to Equation [START_REF]Qmu = mu * P 14 Qmu = Qmu.col_insert(2, zeros(2, 1))[END_REF]. This is not a clear criterion, as several shapes are possible for E R with no clear winner. We propose to define E R as the smallest homothety of E R centered around 0 that satisfies [START_REF]Qmu = mu * P 14 Qmu = Qmu.col_insert(2, zeros(2, 1))[END_REF]. It can be computed rather easily, for any number of dimension; we give details for two dimensions.

Let a, b, c be the coefficients of E R :

E R = {(x c1 , x c2 ) | ax c 2 1 + bx c 2 2 + cx c1 x c2 ≤ 1}.
a, b and c are known, positive values. Then there exists k ≥ 0 such that

E R = {(x c1 , x c2 ) | ax c 2 1 + bx c 2 2 + cx c1 x c2 ≤ k}.
As E R is wider than E R , k ≥ 1. We need a condition on k that guarantees [START_REF]Qmu = mu * P 14 Qmu = Qmu.col_insert(2, zeros(2, 1))[END_REF].

We consider a point (x c1 , x c2 ) located on the border of

E R : ax c 2 1 + bx c 2 2 + cx c1 x c2 = 1. (20)
By construction, for any values e 1 , e 2 ∈ R such that |e 1 | ≤

ε 1 ∧ |e 2 | ≤ ε 2 the relation (x c1 + e 1 , x c2 + e 2 ) ∈ E R must hold, that is to say: a(x c1 + e 1 ) 2 + b(x c2 + e 2 ) 2 + c(x c1 + e 1 )(x c2 + e 2 ) ≤ k.
It develops into

(ax c 2 1 + bx c 2 2 + cx c1 x c2 ) + (2ae 1 + ce 2 )x c1 + (2be 2 + ce 1 )x c2 + (ae 2 1 + be 2 2 + ce 1 e 2 ) ≤ k, that is 1 + (2ae 1 + ce 2 )x c1 + (2be 2 + ce 1 )x c2 + (ae 2 1 + be 2 2 + ce 1 e 2 ) ≤ k due to (20).
Greatest values for the left-hand term are reached with

|e 1 | = ε 1 ∧ |e 2 | = ε 2 ,
depending on the signs of x c1 and x c2 . As the ellipse E R is symmetric around the origin point (0, 0), we can set e 1 = ε 1 , which leaves only two cases to study. Finally, we numerically verify that greatest values of the term are reached when e 2 = ε 2 . This is the only case we detail here.

We can write:

1 + αx c1 + βx c2 + γ ≤ k with values α = (2aε 1 + cε 2 ), β = (2bε 2 + cε 1 ) and γ = (aε 2 1 + bε 2 2 + cε 1 ε 2 ). We know that x c1 and x c2 are bounded, thus so is the term αx c1 + βx c2 : we can compute a minimum bound δ such that αx c1 +βx c2 ≤ δ. So it is sufficient that k satisfies:

k ≥ 1 + γ + δ.
Consequently, the smallest homothety of E R that satisfies ( 19) is obtained with k = 1 + γ + δ; we take it as our definition of E R . The instruction becomes:

11 xc = Ac*xc + Bc*yc; % x c ∈ E R
In our case, starting from the ellipse that gives E R .

E R = {(x c1 , x c2 ) | 0.00269007x c 2 1 + 0.000341414x c1 x c2 + 0.00247323x c
3) End of Proof Scheme: Then, what remains of the stability proof scheme becomes:

% x c ∈ E R 12 send(u, 1); % x c ∈ E R 13 receive(y, 2); % x c ∈ E R 14 receive(yd, 3); % x c ∈ E R 15 skip; % x c ∈ E P 16 end
At this stage, the final assertion E R ⊂ E P holds and is checked successfully.

VI. Closed-Loop Stability Proof Scheme with

Floating-Point Numbers

We now show how the proof of state boundedness of the closed-loop system specifications can be migrated to the level of the controller code and executable model of the system. To be more precise, we exploit the invariance of the ellipsoid E P to develop a proof of proper behavior, that is, stability and variable boundedness, for the computer program that implements the controller as it interacts with the physical system. Unlike the developments related to the open-loop controller, this proof necessarily involves the presence of a model of the physical system. In [START_REF] Feron | From control systems to control software[END_REF], Feron chooses to represent the physical system and the computer program by two concurrent programs. The code for controller dynamics is unchanged, same as in Section II. The pseudo-code to represent the physical system dynamics is shown below. In this scheme, the computer program representation of the physical system is to remain unchanged, since it only exists for modeling purposes and does not correspond to any actual program, whereas the controller code is allowed to evolve to reflect the various stages of its implementation.

Establishing proofs of stability of the closed-loop system at the code level is necessarily tied to understanding the joint behavior of the controller and the plant. The entire state space therefore consists of the direct sum of the state spaces of the controller and the physical system. The approach described in the previous sections is used to document the corresponding system of two processes. One interesting aspect of these processes is their concurrency, which can complicate the structure of the state transitions. However, a close inspection of the programs reveals that the transition structure of the processes does not need to rely on the extensions of Hoare's logic to concurrent programs: one program at a time is running, through the blocking nature of the receive primitive.

Feron's stability proof with real numbers is much longer than for the open-loop system. We do not detail it, the interested reader is referred to [START_REF] Feron | From control systems to control software[END_REF] for full information. To be noticed, the resulting comments are not much more complex than those available from the study of the controller alone. On the good side, as already mentioned, the Hoare formalism is not significantly affected by the concurrent structure of the closed-loop system.

A floating-point representation of the closed-loop system consists of keeping the listing corresponding to the physical system in its original settings, while replacing the controller part with the corresponding floating-point implementation, as we did in Section V. Using similar techniques to the study of the controller alone, proof invariants can be tweaked to take into account constant changes and rounding errors resulting from the use of floatingpoint arithmetic in the controller and real arithmetic in the plant.

Unfortunately, these invariants are not sufficient to show that the stability condition holds at the end of the loop body in the case of 64-bit floating-point numbers. In this case, our tool cannot prove the system stability on a double-precision floating-point architecture: either the system is not stable with the floating-point based controller, and in this case the proof parameters (E P , µ, …) must be chosen more carefully by the controller designer, or the stability holds but we lost it by overapproximating the errors. Tuning the precision of floating-point arithmetic, as described in Section V-A, learns us that the controller stability holds on a quadruple-precision platform (i.e., with 128 bits). Thus, an alternative solution might be to use a microcontroller device with superior precision.

VII. Alternative Limited-Precision Arithmetics

Our general idea is to replace some of the invariants in the original proof scheme by wider ones that include rounding errors, with the hope that the stability condition is strong enough and still holds. This approach is made possible by the fact that the rounding errors introduced by the operations used in the code are bounded on bounded inputs and bounded state variables.

In previous sections, we mostly dealt with the case of a floating-point representation of real numbers. They are not available on all architectures, especially on microcontrollers that are commonly used to implement control systems. In this section, we quickly discuss alternative realnumber representations.

• We can deal with fixed-point arithmetic the same way we do with floating-point, as long as we stand far enough from extremal values that can lead to overflows or extremely large error terms;

• Another way to represent real numbers is to use two integers, a numerator and a denominator. Considering that the input values are exact, the elementary operations do not introduce rounding errors but can easily lead to overflows, e.g. when computing p 1 q 1 + p 2 q 2 = p 1 q 2 + p 2 q 1 q 1 q 2 .

A strategy must be used to prevent overflows by introducing approximations: in this case, the question is to quantify the errors introduced by these approximations.

In our example, we exclusively used additions and multiplications: divisions are not involved in linear control. Still, programs with divisions can also be analyzed, if the numerator can be shown to be far enough from zero: it is a supplementary constraint, but it is reasonable to assume that it should be respected on a realistic control system that uses divisions. Differentiable, periodic functions such as (sin) can be computed with an abacus and an interpolation function, thus with bounded error. In the same way, functions not periodic, but restricted to finite domains, can also be approximated. Other functions, such as tangent or square root, could raise more issues.

VIII. Conclusion

In this paper, we described a theoretical framework to translate proof invariants on code with real arithmetic to similar invariants on limited-precision numbers, while preserving the overall proof structure. We focused on the case of Lyapunov-theoretic proofs of stability with floatingpoint numbers, for which we implemented LyaFloat, a tool that automatically generates correct invariants for floating-point arithmetic from the provided invariants on real arithmetic, attempts to check whether stability holds and computes the required amount of precision needed for this. We used this tool to analyze Feron's motivating proof of stability, thus addressing the issue of arithmetic accuracy when implementing high-level, proved code on a microcontroller.

Many directions are open for future work. A first idea would be to implement the translation program in Coq rather than Python. This would grant an additional layer of safety, at the price of a deeper formalizing of invariant propagators -or even better proving them -and dealing with the quirks of the formalization of floating-point arithmetic, e.g. using the Flocq library [START_REF] Boldo | Flocq[END_REF]. Alternatively, our tool could simply generate invariants in the form of proof terms that could be checked with Coq or another floating-point-compliant proof checker. We also plan to extend the application scope of our tool, enabling it to analyze a wider range of control programs (which means handle more functions and propagation theorems), so that we could apply it to some of those real-life controllers that come with Lyapunov stability proofs [START_REF] Martin | Estimation d'etat pour les drones aeŕiens[END_REF]. Finally, it would be interesting to be able to deal with different arithmetic paradigms, as discussed in Section VII.

Figure 1 .

 1 Figure 1. Mass-spring system

1

  Ac = [0.4990, -0.0500; 0.0100, 1.0000]; 2 Bc = [1; 0]; 3 Cc = [564.48, 0]; 4 Dc = -1280; 5 xc = zeros(2, 1); 6 receive(y, 2); receive(yd, 3); 7 while (1) 8 yc = max(min(y -yd, 1), -1);

Figure 2 .

 2 Figure 2. The stability domain E P

Figure 4 .

 4 Figure 4. Abstract scheme of invariant propagation on the original program

) 4 )

 4 Domains: A domain d on K is a set of valuations on K, i.e. d ⊂ V K . Domains are used to represent Floyd's and Hoare's invariants seen in previous sections. The set of domains on K is noted D K , with D F ⊂ D R . A real domain d ∈ D R can be associated to a limitedprecision domain d ∈ D F defined by:

Figure 5 .

 5 Figure 5. General translation scheme

Figure 6 .

 6 Figure 6. Refined translation scheme

Finally, a valuation

  v err ∈ V F is an upper bound of the arithmetic error of instruction ĩ on a domain d ∈ D F if: ∀v ∈ d, err(ĩ, v) ≤ v err . Now, let us consider an instruction ĩ in I R surrounded by a precondition d and a postcondition d′ = p( d, ĩ) in D R : eval(ĩ, d) ⊂ d′ Let v err ∈ V F be an upper bound of the arithmetic error of instruction ĩ on domain d. Then, the implication

2 2 ≤

 2 1} we get the following values: α = 1.03246•10 -17 , β = 1.84829•10 -17 , γ = 7.17582•10 -32 . Our program finds δ = 5.35754 • 10 -16 ≫ γ and finally k = 1 + 5.35754 • 10 -16

1

  Ap = [1.0000, 0.0100; -0.0100, 1.0000]; 2 Bp = [0.00005; 0.01]; 3 Cp = [1, 0]; 4 while (1) 5 yp = Cp * xp; 6 send(yp, 2); 7 receive(up, 1); 8 xp = Ap * xp + Bp * up; 9 end

: each program instruction comes with an invariant. The program annotated by Feron is reproduced below. 5 xc = zeros(2,1);

	% x c ∈ E P 6 receive(y, 2); receive(yd, 3); % x c ∈ E P 7 while (1) % x c ∈ E P 8 yc = max(min(y -yd, 1), -1); % x c ∈ E P , y 2 c ≤ 1
	9

Mathematica source file is available at: http://www.cri.ensmp.fr/ people/maisonneuve/lyafloat/resources/lyafloat_stability.nb, and the corresponding PDF file at: http://www.cri.ensmp.fr/people/ maisonneuve/lyafloat/resources/lyafloat_stability.pdf.

The procedure that follows would be exactly the same with

32-bit floating-point numbers, only with different numerical results.

B. Converting Constants

The first step when translating Feron's program is to convert real constants into floating-point numbers. The first lines of code in the program:

1 Ac = [0.4990, -0.0500; 0.0100, 1.0000]; Apart from constants, the proof scheme for the first part of the program is unchanged:

C. Converting Functions

The second step of our proof translation scheme is to convert function symbols from real to floating-point. We show how it is done instruction by instruction.

1) Invariant on u:

The next instruction in the original proof scheme is:


Hélène Agogué

Maaike Brink

Julie Dinasquet

Gerhard J Herndl

Major gradients in putatively nitrifying and non-nitrifying Archaea in the deep North Atlantic

come   L'archive ouverte pluridisciplinaire

Planktonic Archaea consisting of the two major groups, Crenarchaeota and Euryarchaeota, might account for about one-third of all prokaryotic cells in the global ocean 7,8 . Crenarchaeota are generally more abundant than Euryarchaeota in marine waters [7][8][9] and mainly composed of the Marine Crenarchaeota Group I (MCGI) and the deep-branching pSL12-like clade 4 . Until recently, only members the bacterial domain, γ-and β-proteobacteria were known to carry out the oxidation of NH 4 + to NO 2 -and γ-, β-, δproteobacteria and Nitrospira species to oxidize the intermediate product NO 2 to NO 3 -1 . Recent genomic studies suggest that mesophilic Crenarchaeota might play a role in the ammonia oxidation as the amoA gene encoding a subunit of one of the key enzymes in the ammonia oxidation, ammonia monooxygenase, has been detected 2,10,11 . The only two mesophilic marine crenarchaeal isolates (i.e.,

Cenarchaeum symbiosum and Nitrosopumilus maritimus) carry also this amoA gene 12,[START_REF] Hallam | Genomic analysis of the uncultivated marine crenarchaeote Cenarchaeum symbiosum[END_REF] . Recent quantitative studies showed that MCGI oxidize ammonia as their energy source and that archaeal amoA genes are, at least, one order of magnitude more abundant than bacterial amoA in various marine as well as in soil environments [3][4][5][START_REF] Leininger | Archaea predominate among ammonia-oxidizing prokaryotes in soils[END_REF][START_REF] Nakagawa | Distribution of the cold-adapted ammonia-oxidizing microorganisms in the deep-ocean of the Northeastern Japan Sea[END_REF] . Moreover, recently, Mincer et al. 4 suggested that the pSL12-like clade may also contain amoA genes. Taken together, there is evidence that at least some of the MCGI are chemoautotrophic ammonia oxidizers although other studies indicate that Crenarcheaota might be heterotrophic as well 8,[START_REF] Ouverney | Marine planktonic Archaea take up amino acids[END_REF][START_REF] Teira | Archaeal uptake of enantiomeric amino acids in the meso-and bathypelagic waters of the North Atlantic[END_REF][START_REF] Kirchman | Standing stocks and activity of Archaea and Bacteria in the western Arctic Ocean[END_REF][START_REF] Ingalls | Quantifying archaeal community autotrophy in the mesopelagic ocean using natural radiocarbon[END_REF] .

All the studies quantifying archaeal amoA gene abundance in the pelagic realm, however, focused on the euphotic and mesopelagic zone where ammonia is generated by remineralisation processes such as in the oxygen minimum zone [3][4][5]11 . Also the only-free living isolate, Nitrosopumilus maritimus, has been obtained from an ammonia-rich tropical aquarium 12 . However, Crenarchaeota are abundant in the bathypelagic (1000-4000 m depth) waters as well [START_REF] Teira | Distribution and activity of Bacteria and Archaea in the deep water masses of the North Atlantic[END_REF][START_REF] Varela | Latitudinal trends of Crenarchaeota and Bacteria in the meso-and bathypelagic water masses of the Eastern North Atlantic[END_REF] where ammonia concentrations are below 5 nM (M. Woodward, personal communication) and only detectable with highly sensitive methods making it unlikely that ammonia is readily available as an energy source in these deep waters.

The contribution of putatively ammonia-oxidizing Archaea (AOA) and Bacteria (AOB) to the total archaeal and bacterial communities was determined by quantitative PCR (Q-PCR) in the main water masses of the Atlantic from 65°N to 5°S covering a depth range from 100 m to 4000 m. The abundance of specific groups of planktonic Crenarchaeota (MCGI and pSL12-like clade) based on 16S rRNA gene fragments was determined and compared with the abundance of AOA and AOB using two specific primer sets targeting archaeal and bacterial amoA genes at 10 stations, comprising different water masses with specific physical and chemical characteristics (for details see Supplementary Fig. 1 and Table 1, 2).

Highest MCGI abundance was found in the mesopelagic layer in almost all the stations and particularly in the oxygen minimum layer (i.e., between 200 and 600 m depth) (Fig. 1a, Supplementary Tables 1 and3). The abundance of the pSL12-like clade was about 1.5 orders of magnitude lower than MCGI (Supplementary Table 3).

In most of the samples (40 out of 55), the contribution of the pSL12-like clade to total crenarchaeal abundance was less than 5% (Supplementary Table 3) but followed a similar vertical distribution pattern as MCGI. The pSL12-like clade was more abundant just below the euphotic zone (i.e., 250 m depth) than in deeper water masses (Supplementary Table 3). A similar trend was observed by Mincer et al. 4 who found more pSL12-related members below the euphotic zone of the north Pacific subtropical gyre. The distribution of the copy numbers of the 16S rRNA genes (MCGI and pSL12-clade) obtained by Q-PCR corresponds to the distribution pattern of Crenarchaeota with depths reported previously using fluorescent in situ hybridization 7,8,[START_REF] Teira | Distribution and activity of Bacteria and Archaea in the deep water masses of the North Atlantic[END_REF][START_REF] Varela | Latitudinal trends of Crenarchaeota and Bacteria in the meso-and bathypelagic water masses of the Eastern North Atlantic[END_REF] . Overall, marine Crenarchaeota seem to be particularly associated with the oxygen minimum zone 4,5,11,[START_REF] Lam | Linking crenarchaeal and bacterial nitrification to anammox in the Black Sea[END_REF] .

Copy numbers of archaeal amoA genes were always highest at the base of the euphotic zone ranging from 2.5 x 10 3 to 25 x 10 3 ml -1 decreasing with depth (Fig. 1b, Supplementary Table 3). While no latitudinal trend in the abundance of the archaeal amoA genes was detectable at the base of the euphotic zone, in deeper waters, archaeal amoA copy numbers decreased by 2-3 orders of magnitude from the subpolar to the subtropical and equatorial Atlantic (Fig. 1b, Supplementary Table 3). The higher abundance of archaeal amoA genes in the meso-and bathypelagic waters in the northern part of the Atlantic coincides with the higher ammonia concentrations there [START_REF] Varela | Latitudinal trends of Crenarchaeota and Bacteria in the meso-and bathypelagic water masses of the Eastern North Atlantic[END_REF] as these deep waters have been recently formed and therefore, are younger than in the southern part of the North Atlantic [START_REF] Rhein | Labrador Sea Water: Pathways, CFC Inventory, and Formation Rates[END_REF] . In contrast to the high number of archaeal amoA copy numbers, β-proteobacterial amoA copy numbers were at least one order of magnitude lower than archaeal amoA gene abundance in the individual depth layers (Supplementary Table 3). Gamma-proteobacterial amoA genes were not quantified in this study since ammonia-oxidizing γ-proteobacteria were not detected in this region of the Atlantic by Q-PCR in a previous study 3 nor in Arctic and North Pacific samples 4,[START_REF] Hollibaugh | Widespread distribution in polar oceans of a 16S rRNA gene sequence with affinity to Nitrosospira-like ammonia-oxidizing bacteria[END_REF] . Taken together, our data confirm previous findings that the crenarchaeal amoA genes are by far more abundant than bacterial amoA in marine systems 3,5,[START_REF] Nakagawa | Distribution of the cold-adapted ammonia-oxidizing microorganisms in the deep-ocean of the Northeastern Japan Sea[END_REF] .

The average ratio of archaeal amoA copy numbers versus crenarchaeal 16S rRNA genes in the subsurface waters (100-150 m depth) was close to 1, particularly in the northern and equatorial North Atlantic, while it was <0.01 in the bathypelagic waters of the subtropical and equatorial region of the North Atlantic (Fig. 1c, Supplementary Fig. 2). Our ratios for the surface waters are similar to those of Mincer et al. 4 reporting a ratio between archaeal amoA and archaeal 16S rRNA genes of about 1 while Wuchter et al. 3 and Beman et al. 5 reported a ratio of 2.8 and 2.5, respectively. Genomic studies on the two mesophilic crenarchaeal isolates obtained thus far revealed that Cenarchaeum symbiosum and as well as Nitrosopumilus maritimus contain only one amoA gene copy per cell [START_REF] Hallam | Genomic analysis of the uncultivated marine crenarchaeote Cenarchaeum symbiosum[END_REF] . In the 100-150 m depth range, the average ratios between archaeal amoA and crenarchaeal 16S rRNA genes were 0.9, 0.4 and 0.7 for the northern, the subtropical and the equatorial region, respectively (Fig. 1c). We conclude therefore, that in subsurface waters most of the Crenarchaeota are putatively capable of ammonia oxidation throughout the North Atlantic. In the northern part of the Atlantic (65°N-30°N), the ratio archaeal amoA to crenarchaeal 16S rRNA genes decreased only by about one order of magnitude from the surface waters to the bathypelagic realm (Fig. 1c). In the subtropical and equatorial regions of the Atlantic (30°N-5°S), however, this ratio decreased by at least two orders of magnitude from the subsurface to the bathypelagic waters (Fig. 1c).

Hence, AOA are apparently only abundant throughout the water column in the northern North Atlantic with its newly formed and therefore, young deep waters [START_REF] Rhein | Labrador Sea Water: Pathways, CFC Inventory, and Formation Rates[END_REF][START_REF] Smethie | Rates of North Atlantic Deep Water formation calculated from chlorofluorocarbon inventories[END_REF][START_REF] Van Aken | The hydrography of the mid-latitude northeast Atlantic Ocean: I: The deep water masses[END_REF] .

The major fraction of the abundant Crenarchaeota in the subtropical and tropical

North Atlantic deep waters most likely use energy sources other than ammonia.

To assess the diversity of AOA in the meso-and bathypelagic waters of the North Atlantic, archaeal amoA clone libraries were established for two mesopelagic samples where Q-PCR analyses revealed relatively high numbers of AOA, and for one bathypelagic sample where the AOA abundance was very low (Supplementary Table 3). The bathypelagic clone library exhibited the lowest observed and estimated operational taxonomic units (OTUs) richness (Supplementary Fig. 3 and Table 4). The richness levels of the two mesopelagic clone libraries were comparable or higher than those observed for archaeal amoA clone libraries from other marine pelagic environments 11 (Supplementary Table 4). In a DNA-based phylogenetic tree, the obtained sequences fell into two distinct groups: cluster A or "shallow" cluster and cluster B or "deep" cluster 5,11,[START_REF] Hallam | Pathways of carbon assimilation and ammonia oxidation suggested by environmental genomic analyses of marine Crenarchaeota[END_REF] (Fig. 2). Based on phylogenetic analyses of PCR clone libraries, individual regions of the mesopelagic North Atlantic harbour specific archaeal amoA sequences (i.e., North Atlantic and equatorial subclusters, Fig. 2).

More than half of the sequences of the clone library from the mesopelagic North Atlantic formed a coherent cluster closely related to sequences obtained from the same region 3 . Similarly, specific sequences recovered from the bathypelagic clone library fell into the bathypelagic subcluster 1 within cluster A and into the bathypelagic subcluster 2 of the phylogenetically distinct cluster B (Fig. 2). About subcluster harbouring amoA sequences from a diverse range of marine habitats [3][4][5]11,[START_REF] Lam | Linking crenarchaeal and bacterial nitrification to anammox in the Black Sea[END_REF] (Fig. 2). Taken together, our phylogenetic analyses indicate latitudinal and vertical differences among archaeal amoA sequences with generally rather low diversity in equatorial bathypelagic waters (Supplementary Fig. 3) coinciding with the rather low abundance of AOA. There is, however, still the possibility that deep ocean Crenarchaeota harbour specific versions of amoA genes not targeted by the thus far available primer sets.

An archaeal 16S rRNA clone library was also established for one mesopelagic sample and 3 bathypelagic samples (for details see Supplementary information). In a DNA-based phylogenetic tree, the obtained crenarchaeal sequences fell into the MCGI and pSL12-like groups (Supplementary Fig. 4). The phylogenetic analyses of the archaeal 16S rRNA genes strengthen also the idea that specific water masses harbour specific ecotypes and that the presence of OTU 2, closely related to Nitrosopumilus maritimus, especially in the TRANSAT sample might explain the high abundance of amoA genes in this specific water mass (for details see

Supplementary information).

Ammonia oxidizing Archaea are generally chemoautotrophic using acetylcoenzyme A carboxylase to fix inorganic carbon [START_REF] Berg | A 3-hydroxypropionate/4hydroxybutyrate autotrophic carbon dioxide assimilation pathway in archaea[END_REF] . Hence, if AOA are the dominating nitrifiers as suggested by the much higher archaeal than bacterial amoA copy numbers, one would expect a tight relation between archaeal amoA copy numbers and bulk CO 2 fixation rates. Relating the dark CO 2 fixation rates obtained at these specific stations and depths by Herndl et al. 8 and Reinthaler and Herndl (unpublished) to the archaeal amoA copy numbers determined in this study, it appears that about 50% of the variation in CO 2 fixation rates can be explained by archaeal amoA abundance (Fig. 3). This provides further independent evidence that the vast majority of the bathypelagic Crenarchaeota are not chemoautotrophic nitrifiers but most likely use organic substrate as carbon and energy source as indicated also by microautoradiography combined with fluorescence in situ hybridization and genomic analyses 8,[START_REF] Ouverney | Marine planktonic Archaea take up amino acids[END_REF][START_REF] Teira | Archaeal uptake of enantiomeric amino acids in the meso-and bathypelagic waters of the North Atlantic[END_REF][START_REF] Kirchman | Standing stocks and activity of Archaea and Bacteria in the western Arctic Ocean[END_REF][START_REF] Varela | Latitudinal trends of Crenarchaeota and Bacteria in the meso-and bathypelagic water masses of the Eastern North Atlantic[END_REF][START_REF] Hallam | Pathways of carbon assimilation and ammonia oxidation suggested by environmental genomic analyses of marine Crenarchaeota[END_REF] . This conclusion is further supported by recent comparative metagenomic analyses of uncultivated Archaea of various meso-and bathypelagic oceanic regions by Martin-Cuadrado et al. [START_REF] Martin-Cuadrado | Hindsight in the relative abundance, metabolic potential and genome dynamics of uncultivated marine archaea from comparative metagenomic analyses of bathypelagic plankton of different oceanic regions[END_REF] . In none of the nine fosmid libraries of MCGI from 500 m to 3000 m depth, amoA genes were detected [START_REF] Martin-Cuadrado | Hindsight in the relative abundance, metabolic potential and genome dynamics of uncultivated marine archaea from comparative metagenomic analyses of bathypelagic plankton of different oceanic regions[END_REF] .

In summary, it appears that only in the northern part of the North Atlantic mesophilic Archaea are predominantly oxidizing ammonia as an energy source throughout the water column including the bathypelagic realm. The abundance of putatively ammonia oxidizing Crenarchaeota gradually decreases in the meso-and bathypelagic waters from the north to the equator in the Atlantic coinciding with an increasing age of these water masses and concomitantly, a decrease in the diversity of archaeal amoA sequences. In the subtropical and equatorial bathypelagic waters, AOA contribute less than 1 % to the crenarchaeal community. Hence, only a minor fraction of the bathypelagic marine Crenarchaeota is putatively oxidizing ammonia in the temperate and subtropical North Atlantic. Consequently, a heterotrophic life mode of bathypelagic Crenarchaeota seems to become increasingly important from north to south in the Atlantic's interior.

Methods summary

Sampling was conducted aboard RV Pelagia during two cruises (TRANSAT-1 and ARCHIMEDES-2) following the North Atlantic Deep Water from 65°N to 5°S in the eastern basin of the Atlantic Ocean. Specific water masses were identified and sampled based on their physical and chemical characteristics.

Using quantitative PCR, the copy numbers of MCGI 16S rRNA, pSL12-like 16S rRNA, archaeal amoA and β-proteobacterial amoA genes were determined using specific primers and SYBRGreen. The resulting mean disintegrations per minute (DPM) of the samples were converted into organic carbon produced over time and corrected for the natural dissolved inorganic carbon (DIC) concentration. Water from 65°N to 5°S in the eastern basin of the Atlantic Ocean. Specific water masses were identified and sampled based on their physical and chemical characteristics (Supplementary Table 1). At each station, 6 depths were sampled with 10 l NOEX samplers mounted on a CTD (conductivity-temperature-depth) frame. Ten stations were used for quantitative PCR analysis (Supplementary Fig. 1).

DNA extraction of marine samples. Extraction of total DNA was performed using

UltraClean soil DNA isolation kit and Mega soil DNA isolation (Mobio) for TRANSAT-1 and ARCHIMEDES-2 samples, respectively.

Quantitative PCR. The copy numbers of MCGI 16S rRNA, pSL12-like 16S rRNA, archaeal amoA and β-proteobacterial amoA genes were determined in triplicate on the non-diluted samples and in two different dilutions per sample using specific primers (Supplementary Table 2) and SYBRGreen. The 4 real-time approaches yielded highly reproducible standard curves with DNA of the plasmid 088exp4, the fosmid HF770_041/11, Nitrosopumilus maritimus and a mix of 4 β-proteobacteria, for MCGI, pSL12-like, archaeal amoA and β-proteobacterial amoA standards, respectively. All reactions were performed in duplicate for the standards and in triplicate for the samples and an average value was calculated for each sample. Two different methods were used to determine the presence of non-specific amplified products. No unspecific products were found in any of the analysis (See Supplementary Information for details).

Archaeal amoA clones libraries. Archaeal amoA clone libraries were generated from two stations: Station A33 at 2 distinct depths (749 m and 2502 m) and Station T28 at 597 m depth by PCR amplification using specific archaeal amoA primers (Supplementary Fig. 1, Table 2). The obtained sequences were compared with sequences from public database and analyzed using the ARB phylogenetic package (See Supplementary information for details). 2). The obtained sequences were compared with sequences from public database and analyzed using the ARB phylogenetic package (see Supplementary information for details).

Archaeal CO 2 fixation. Dark CO 2 fixation was measured via the incorporation of Archaeal amoA (copies.ml -1 )

CO 2 fixation (m- 3 .day -1 )
2750 m and 4000 m depth, respectively. For each depth, ten liters of seawater were filtered through a 0.22 µm Sterivex filter GP unit (Millipore) and 1.8 ml of lysis buffer (40 mM EDTA, 50 mM Tris-HCl, 0.75 M sucrose) was added and the filters subsequently stored at -80°C until analysis.

Molecular analysis

DNA extraction and preparation of the Q-PCR standards. Archaeal 16S rRNA gene fragments were amplified from the plasmid 088exp4 (from the archaeal clones library) and the fosmid HF770_041|11 (obtained from T. Mincer 1 ) as 16S rDNA, Marine Crenarchaeota Group I (MCGI) and pSL12-like standard, respectively, using the 16S specific primers (Parch519f and ARC915R) 2 . The archaeal amoA gene of the standard was amplified from DNA of Nitrosopumilus maritimus (obtained from M. Könneke 3 ), using specific archaeal amoA primers (Supplementary Table 2). The β-proteobacterial amoA gene for the standard was amplified using a mix (25% each) of Nitrosomonas europaea, N. eutropha, N. marina and Nitrospira briensis using specific β-proteobacterial amoA primers (Supplementary Table 2).

Each amplification was performed under the following conditions: 4 min initial denaturation; 35 cycles at, 94°C for 30 s, specific annealing temperature of the primer set for 40 s (Supplementary Table S2), 72°C for 2 min, 80°C for 25 s using for 1 U of Pico Maxx high fidelity DNA polymerase (Stratagene), 10x Pico Maxx PCR buffer, 0.25 mM of each dNTP, 8 µg of BSA, 0.2 µM of primers, 3mM of MgCl 2 and ultra pure sterile water (Sigma).

Amplification products were run on an agarose gel (1%), stained with SYBRGold ® (Invitrogen), bands were isolated and purified using QuickClean 5M gel extraction kit (GenScript). Purified products were quantified using a Nanodrop ® spectrophotometer and the copy numbers of the 16S rRNA genes of each QPCR standard were calculated directly from the concentration of the purified DNA. Ten-fold serial dilutions ranging from 10 7 to 10 2 copies of a known copy number of the standard were used per QPCR reaction in duplicate to generate an external quantification standard. Accumulation of newly amplified double stranded gene products was followed online as the increase of fluorescence due to the binding of the fluorescent dyes SYBRGreen ® and fluorescein. Specificity for QPCR reaction was tested on agarose gel electrophoresis and with a melting curve analysis (60°C-94°C) in order to identify unspecific PCR products such as primer dimers or fragments with unexpected fragment lengths. PCR efficiencies and correlation coefficients for standard curves were as follows: for the MCGI 16S rDNA assay, 87.9-105.8% and r² = 0.991-0.999, for the pSL12-like 16S rDNA assay, 94.7-110.5% and r² = 0.949-0.999, for the archaeal amoA assay, 86.5-97% and r² = 0.983-0.999 and for the βproteobacterial amoA assay, 77.8-88.3% and r² = 0.989-0.999. Each gene fragment was detected using a standard for the specific quantification of MCGI 16S rRNA, pSL12-like 16S rRNA, archaeal amoA and β-proteobacterial amoA genes and primer combinations and annealing temperature as listed in Table 3. Thermocycling was performed as follows: initial denaturation at 95°C for 4 min; amplification: 41 cycles at, 95°C for 30 s, primer annealing temperature for 40 s, and extension at 72°C for 30 s, 80°C for 25 s, with a plate read between each cycle; melting curve 60 -94°C with a read every 0.5°C held for 1 s between each read.

DNA extraction and quantitative

AmoA gene amplification, sequencing and phylogenetic analysis of archaeal amoA sequences. Archaeal amoA sequences were amplified with the primers Arch-amoA-for and Arch-amoA-rev (Supplementary Table 2) as described in the Q-PCR section and 35 cycles of amplification. AOA amoA PCR were purified (QuickClean 5M PCR purification Kit, GenScript) and cloned into the pCR4-TOPO ® vector for sequencing (Invitrogen) according to manufacturers instructions. Clones were checked for the right insert by PCR and sequencing was performed (Eurofins MWG GmbH, Germany). Sequence data were compiled using ARB software 5 and aligned with complete length sequences of closest relatives obtained from NCBI database. Using ARB, the phylogenetic tree was first generated with the aligned, almost complete length sequences of closest relatives from the NCBI data base using the neighbor-joining method. Then the short-aligned sequences were added to the tree using the ARB parsimony option.

To compare the archaeal amoA-based richness within each clone library, rarefaction analysis and Chao1 non parametric richness estimations were performed by using DOTUR 6 .

Operational taxonomic units (OTUs) were defined as sequence group in which sequences differed by ≤ 2% (Supplementary Fig. 3 and Table 4). 3). The 16S rDNA cloning was performed as described above. Seventy-seven, 90, 97 and 80 archaeal clones from the stations T9 at 2071 m, A3 at 249 m, A3 at 2750 m, and A33 at 7155 m, respectively, were sequenced. Among these clones, 44 (57%), 50 (56%), 83 (86%) and 32 (40%) clones from the stations T9 at 2071 m, A3 at 249 m, A3 at 2750 m, and A33 at 7155 m, respectively, were affiliated to the crenarchaeal group. A DOTUR analysis was performed on the 219 crenarchaeal sequences, 19 different OTUs were defined as sequence groups in which sequences differed by ≤ 3%.

Archaeal

The sequences were compared with sequences from public database and analyzed using the ARB phylogenetic package (Supplementary Fig. 4).

Primers specificity.

16S rRNA genes primers. To test whether groups not targeted by our MCGI and pSL12 primers might be responsible for the observed Q-PCR pattern, archaeal 16S rRNA gene fragments were amplified, cloned, screened and sequenced (see above). Using the ARB phylogenetic package 5 , the specificity of the primers was tested on our MCGI and pSL12clade affiliated sequences from the four archaeal clone libraries. The primers used in our quantitative PCR approach for detecting MCGI and the pSL12-like cluster targeted all the obtained sequences (data not shown). The specificity of the primers used by Mincer et al. 1 was also tested in silico on our MCGI sequences. These primers targeted also all the MCGI sequences (data not shown).

AmoA genes primers. To test whether groups of archaeal amoA sequences not targeted by the archaeal amoA primers used in this study 8 amplified by these primers is at the end of the amoA gene and is not or not entirely recovered by the primers used by Wuchter et al. 8 and Francis et al. 9 .

Using the primers of Francis et al., we recovered 38 amoA sequences with a length of 595 bp affiliating exclusively to cluster B (Supplementary Fig. 5). Half of these sequences are closely related to the bathypelagic subcluster 2. The sequences belonging to cluster A recovered in our North Atlantic clone library are not detected by this set of primers. Using this set of primers for Q-PCR analysis in North Atlantic waters would underestimate the number of AOA. Hence, the archaeal amoA primers used in our study provide a better coverage of the entire group of marine archaeal amoA sequences (clusters A and B, Fig. 2), however, do not recover some subclusters within cluster B which are covered by the primers of Francis et al. 9 (Supplementary Fig. 5). Recently, Beman et al. 10 designed two new forward primers, Arch-amoAFA and Arch-amoAFB, able to amplify specifically cluster A and B of AOA, respectively. The long and closest neighbor sequences of our North Atlantic sequences (cluster A and B) are recovered in silico by these primers (Fig. 2: EF414272, EU049823, EU049824, EU049846, AB373359, EU049842, AM901414 and AB305361). Thus, we conclude that our sequences would also be recovered by the set of primers of Beman et al. 10 .

Archaeal

  amoA clone libraries were generated from two stations: Station A33 at 2 distinct depths (749 m and 2502 m) and Station T28 at 597 m depth by PCR amplification using specific archaeal amoA primers. Archaeal 16S rRNA gene clone libraries were generated from 3 stations: Station T9 at 2071 m, Station A3 at 2 distinct depths (249 m and 2750 m) and one additional station in the Romanche Fracture Zone (-0.21°N, -18.43°E) near the station A33 (1°N, -20.60°E) at 7155 m depth by PCR amplification using specific archaeal primers. Dark archaeal CO 2 fixation was measured via the incorporation of [ 14 C]-bicarbonate.

Figure 1 │

 1 Figure 1 │ Latitudinal and depth distribution of crenarchaeal 16S rRNA and archaeal amoA copy numbers in the eastern North Atlantic. a, b. Abundance of MCGI (a) and AOA (b) in the subsurface (n=9), mesopelagic (Meso) (n=24) and bathypelagic (Bathy) (n=22) waters. c. Ratio of archaeal amoA copies to crenarchaeal (MCGI + pSL12-like clade) 16S rRNA copies per ml of marine water (n=55). Bars indicate standard error.

Figure 2 │

 2 Figure 2 │ Phylogenetic tree of archaeal amoA sequences recovered from the North Atlantic. Neighbor-joining bootstrap tree of long amoA-like nucleotide sequences to which we added our partial sequences (217 bp, EU795424-EU795460) and the sequences of Wuchter et al. 3 (217 bp, DQ784528-DQ784537) using the ARB parsimony tool (NA: North Atlantic). One representative of sequence group ≥99% identical is shown (in color); the number in front of the respective symbol indicates the number of clones. Database sequences are shown in black. Bootstrap values (>50%) are indicated at branch point.

Figure 3 │

 3 Figure 3 │ Archaeal amoA copy numbers vs. CO 2 fixation. Relationship between archaeal amoA gene copy numbers and CO 2 fixation across all the stations and depths (except station T9) (n=31). Regression equation: y=0.0063x 0.5905 .

Archaeal

  16S rRNA gene clone libraries. Archaeal 16S rRNA gene clone libraries were generated from 3 stations: Station T9 at 2071 m, Station A3 at 2 distinct depths (249 m and 2750 m) and one additional station in the Romanche Fracture Zone (-0.21°N, -18.43°E) near the station A33 (1°N, -20.60°E) at 7155 m depth by PCR amplification using specific archaeal primers (Supplementary Fig. 1, Table

[ 14 CFigure 1 .

 141 Figure 1.

  PCR of the environmental samples. DNA extraction. Total DNA extraction from marine samples was performed using the UltraClean soil DNA isolation kit (Mobio) for the TRANSAT-1 samples and the UltraClean Mega soil DNA isolation kit (Mobio) followed by a DNA concentration with a Centricon unit (Millipore) for ARCHIMEDES-2 samples. Q-PCR analysis. All Q-PCR analyses were performed on an iCycler iQ 5 thermocycler (Biorad) equipped with i-Cycler iQ software (version 3.1, Bio-Rad). The copy numbers of MCGI 16S rDNA, pSL12-like 16S rDNA, archaeal amoA and β-proteobacterial amoA in all environmental samples were determined in triplicate on the non-diluted sample and for two different dilutions of the sample (5 times and 25 times diluted). The reaction mixture (20 µl) contained 1 U of Pico Maxx high fidelity DNA polymerase (Stratagene), 2 µl of 10x Pico Maxx PCR buffer, 0.25 mM of each dNTP, 8 µg of BSA, 0.2 µM of primers, 50 000 times diluted SYBR Green ® (Invitrogen) (optimized concentration), a final concentration of 10 nM of fluorescein, 3mM of MgCl 2 and ultra pure sterile water (Sigma) 4 . All reactions were performed in a 96-well QPCR plates (Bio-Rad) with optical tape (Bio-Rad). One µL of diluted or non-diluted environmental DNA was added to 19 µL of mix in each well.

  16S rRNA gene amplification, sequencing and phylogenetic analysis of archaeal 16S rRNA sequences. Archaeal 16S rRNA gene clone libraries were generated from 3 stations: Station T9 at 2071 m, Station A3 at 2 distinct depths (249 m and 2750 m) and one station (ARCHIMEDES-3 cruise, 0.21°N, -18.43°E) close to the station A33 (1°N, -20.60°E) at 7155 m depth by PCR amplification using a forward primer specific for Archaea and a reverse universal primer (Arch21f, 1492r) 7 (Supplementary Table

  might be responsible for the observed Q-PCR pattern, we established two clone libraries with the primers used by Mincer et al. 1 and Francis et al. 9 . These two clone libraries were established from the same bathypelagic sample as used with our primers (Station A33, 2502 m depth). From the 45 short sequences (only 85 bp) recovered from the clone library with the Mincer et al. 1 set of primers, only 2 sequences fell into cluster A, all the other sequences affiliated with cluster B (data not shown). The relatively short sequences of the amoA gene amplified by these primers do not allow phylogenetic analyses. Mincer and coauthors used this set of primers to determine the abundance of archaeal amoA genes by Q-PCR. The short fragment

  Figure S1.

Table S1 :

 S1 Characteristics of the water masses sampled during the TRANSAT-1 and ARCHIMEDES-2 cruise based on CTD a profiles at the stations (Archaeal amoA cloning and archaeal 16S rRNA gene cloning were performed at stations in bold and italic, respectively). North Atlantic Central Water ; SACW, South Atlantic Central Water ; transitional CW, transitional Central Water ; LSW, Labrador Sea Water ; ISOW, Iceland-Scotland Overflow Water ; SAIW, Subarctic Intermediate Water ; AAIW, Antarctic Intermediate Water ; transitional IW, transitional intermediate water ; UNEADW, Upper North East Atlantic Deep Water ; MNEADW, Middle North East Atlantic Deep Water ; LNEADW, Lower North East Atlantic Deep Water, AABW, Antarctic Bottom Water.Parameters: temp, temperature (°C) measured with the CTD rosette ; salinity, salinity (pss) measured with the CTD rosette, pot temp, potential temperature (°C) ; density, potential density anomaly (kg.m -3 ) ; CTD oxy, oxygen concentration (µmol.kg -1 ) measured with the CTD rosette ; AOU, apparent oxygen utilization (µmol.kg -1 ).

	ARCHIMEDES-2	11	11.98	-20.21	250	transitional CW	12.101	35.207	12.068	26.738	86.5	176.013
	ARCHIMEDES-2	11	11.98	-20.21	100	Subsurface	14.206	35.446	14.192	26.492	73.9	177.046
	ARCHIMEDES-2	19	4.90	-15.19	4001	LNEADW	2.368	34.892	2.023	27.883	239.1	92.827
	Cruise ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-3 TRANSAT-1 ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-2 TRANSAT-1 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 a : CTD, conductivity, temperature, depth Station Latitude (°N) 19 4.90 9 61.63 19 4.90 9 61.63 19 4.90 9 61.63 19 4.90 9 61.63 19 4.90 9 61.63 28 -3.39 20 55.91 28 -3.39 20 55.91 28 -3.39 20 55.91 28 -3.39 20 55.91 28 -3.39 20 55.91 28 -3.39 28 49.73 23 -0.21 28 49.73 33 1.00 28 49.73 33 1.00 28 49.73 33 1.00 28 49.73 33 1.00 3 24.50 33 1.00 3 24.50 33 1.00 3 24.50 48 19.19 3 24.50 48 19.19 3 24.50 48 19.19 5 21.32 48 19.19 5 21.32 48 19.19 5 21.32 48 19.19 5 21.32	Longitude (°E) -15.19 -20.82 -15.19 -20.82 -15.19 -20.82 -15.19 -20.82 -15.19 -20.82 -14.83 -32.21 -14.83 -32.21 -14.83 -32.21 -14.83 -32.21 -14.83 -32.21 -14.83 -26.13 -18.43 -26.13 -20.60 -26.13 -20.60 -26.13 -20.60 -26.13 -20.60 -18.69 -20.60 -18.69 -20.60 -18.69 -25.24 -18.69 -25.24 -18.69 -25.24 -20.87 -25.24 -20.87 -25.24 -20.87 -25.24 -20.87	Depth (m) 2751 2071 900 1633 500 1014 251 498 100 99 2501 2500 1799 1897 901 1147 498 450 252 149 100 3878 7155 2999 2502 1897 1800 597 749 146 498 2750 250 800 99 499 2500 249 1798 99 901 4001 499 2748 249 899 99 501	Water mass MNEADW ISOW AAIW slope water SACW AAIW SACW NACW Subsurface Subsurface MNEADW ISOW UNEADW ISOW AAIW LSW SACW NACW SACW Subsurface Subsurface LNEADW AABW MNEADW MNEADW LSW UNEADW NACW AAIW Subsurface SACW MNEADW SACW SAIW Subsurface NACW MNEADW NACW UNEADW Subsurface SAIW LNEADW NACW MNEADW NACW SAIW Subsurface NACW	Temp 2.787 2.571 5.117 3.390 7.938 5.644 13.183 7.877 16.219 9.100 2.948 2.942 3.781 3.184 4.419 3.757 7.499 6.257 11.901 7.508 16.398 2.707 1.4212 2.873 2.962 3.233 3.900 5.701 5.307 11.214 7.195 2.920 11.987 8.342 16.469 11.824 3.196 15.902 4.091 19.751 6.707 2.368 10.193 2.954 13.918 7.069 18.733 11.371	Salinity 34.926 34.983 34.629 34.952 34.758 35.065 35.314 35.196 35.655 35.228 34.925 34.974 34.971 34.972 34.543 34.915 34.658 35.048 35.139 34.938 35.689 34.928 34.752 34.943 34.934 34.888 34.982 34.946 34.513 35.439 34.637 34.952 35.150 35.188 35.666 35.603 34.963 36.259 35.007 36.898 34.939 34.897 35.245 34.956 35.735 35.045 36.508 35.506	Pot temp 2.566 2.416 5.041 3.263 7.887 5.552 13.148 7.825 16.203 9.089 2.749 2.739 3.635 3.035 4.349 3.669 7.450 6.216 11.868 7.493 16.382 2.359 0.700 2.620 2.762 3.084 3.753 5.649 5.244 11.196 7.147 2.697 11.955 8.256 16.453 11.758 2.993 15.862 3.942 19.732 6.621 2.023 10.133 2.730 13.881 6.980 18.716 11.306	Density 27.865 27.924 27.377 27.822 27.101 27.660 26.607 27.455 26.205 27.283 27.848 27.888 27.801 27.860 27.386 27.753 27.086 27.563 26.723 27.301 26.190 27.884 27.870 27.874 27.855 27.788 27.798 27.555 27.260 27.082 27.112 27.875 26.715 27.383 26.156 27.105 27.857 26.749 27.798 26.282 27.423 27.888 27.121 27.874 26.781 27.457 26.248 27.114	CTD oxy 236.0 292.3 127.1 287.4 86.9 247.8 79.0 265.9 86.2 264.2 239.8 284.0 240.3 284.9 149.8 278.8 108.2 237.4 94.9 264.2 109.6 253.0 222.4 269.8 245.0 286.7 244.8 225.6 142.9 240.7 127.5 230.9 103.3 117.1 117.0 145.1 234.4 183.9 212.3 216.1 115.5 244.3 64.0 234.3 83.5 126.8 137.5 109.4	91.279 AOU 181.609 36.100 201.628 34.200 177.496 56.100 154.580 22.100 85.957 15.800 78.298 41.795 164.233 38.444 183.404 39.707 168.898 61.839 130.293 26.546 121.400 76.054 80.673 56.977 72.817 36.515 164.551 77.905 166.187 26.282 160.023 95.237 122.601 168.165 89.284 118.393 103.796 57.515 181.197 7.355 209.541 87.562 168.521 91.597 90.726 167.168 156.863
	ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 ARCHIMEDES-2 Water masses: NACW, 9 5 21.32 -20.87 249 NACW 13.821 35.763 13.785 26.823 92.3 160.108 5 21.32 -20.87 100 Subsurface 15.642 35.869 15.627 26.503 56.7 186.391 11 11.98 -20.21 3999 LNEADW 2.365 34.894 2.020 27.886 242.2 89.782 11 11.98 -20.21 2751 MNEADW 2.859 34.938 2.636 27.869 237.5 89.244 11 11.98 -20.21 900 transitional IW 5.993 34.778 5.912 27.389 105.0 196.989 11 11.98 -20.21 499 transitional CW 9.680 35.076 9.622 27.076 44.0 232.912

Table S2 -

 S2 Primers used for detection and quantification of archaeal 16S rRNA and amoA genes and bacterial amoA genes

				Fragment	Annealing
	Target	Gene	Name	length	temperature	Sequence (5' to 3')	Reference
				(bp)	(°C)
	Marine Crenarchaeota	16S rRNA 16S rRNA	MCGI-391f MCGI-554r	122	61.0	AAGGTTARTCCGAGTGRTTTC TGACCACTTGAGGTGCTG	8 8
	pSL12-like clade	16S rRNA 16S rRNA	pSL12_750F pSL12_876R	326	65.0	GGTCCRCCAGAACGCGC GTACTCCCCAGGCGGCAA	1 1
	Archaeal amoA	amoA amoA	Arch-amoA-for Arch-amoA-rev	256	58.5	CTGAYTGGGCYTGGACATC TTCTTCTTTGTTGCCCAGTA	8 8
	β-proteobacterial amoA	amoA amoA	amoA-1F amoA-r NEW	490	61.5	GGGGTTTCTACTGGTGGT CCCCTCBGSAAAVCCTTCTTC	11

Table S3 :

 S3 Archaeal 16S rRNA and amoA gene, and β-proteobacterial amoA gene copy numbers in the samples from the Atlantic Ocean MW: marine water ; n.d.: not detected (i.e., no specific product formed upon Q-PCR or below detection limit).

	ARCHIMEDES-2	11	3999	1.44x10 4 (2.78x10 3 )	4.38x10 2 (2.07x10 2 )	1.03x10 1 (1.73x10 0 )	1.41x10 0 (7.43x10 -1 )
	ARCHIMEDES-2	11	2751	2.54x10 4 (1.15x10 3 )	1.92x10 2 (1.04x10 1 )	3.40x10 1 (5.97x10 0 )	2.56x10 0 (1.32x10 0 )
	ARCHIMEDES-2	11	900	MCGI 16S rDNA 1.21x10 5 (4.37x10 4 )	pSL12-like 16S rDNA 2.13x10 2 (6.17x10 1 )	Archaeal amoA gene 3.75x10 1 (1.08x10 1 )	β-proteobacterial amoA gene 1.33x10 1 (-)
	Cruise ARCHIMEDES-2	Station 11	Depth (m) 499	copies ml -1 MW (S.D.) 1.50x10 5 (3.05x10 4 )	copies ml -1 MW (s.d.) 4.95x10 2 (2.84x10 1 )	copies ml -1 MW (s.d.) 1.92x10 2 (3.50x10 1 )	copies ml -1 MW (s.d.) 5.30x10 1 (8.84x10 0 )
	TRANSAT-1 ARCHIMEDES-2	9 11	2071 250	7.84x10 4 (7.50x10 3 ) 4.02x10 5 (7.04x10 4 )	9.33x10 2 (7.52x10 1 ) 2.29x10 3 (3.64x10 2 )	2.36x10 3 (5.06x10 2 ) 2.57x10 3 (6.73x10 2 )	2.21x10 2 (5.13x10 1 ) 1.13x10 2 (3.84x10 1 )
	TRANSAT-1 ARCHIMEDES-2	9 11	1633 100	1.39x10 4 (1.83x10 3 ) 2.35x10 5 (3.53x10 4 )	8.70x10 2 (2.63x10 2 ) 2.87x10 2 (5.75x10 1 )	4.98x10 2 (2.52x10 1 ) 2.53x10 4 (7.42x10 2 )	7.20x10 1 (1.88x10 1 ) 1.87x10 1 (-)
	TRANSAT-1 ARCHIMEDES-2	9 19	1014 4001	1.35x10 5 (4.01x10 4 ) 2.86x10 3 (6.52x10 2 )	3.74x10 3 (1.20x10 3 ) 5.73x10 2 (1.48x10 2 )	2.39x10 3 (3.87x10 2 ) 3.80x10 0 (2.28x10 0 )	6.28x10 1 (1.74x10 1 ) n.d.
	TRANSAT-1 ARCHIMEDES-2	9 19	498 2751	3.99x10 4 (1.68x10 4 ) 3.23x10 4 (6.54x10 3 )	1.81x10 3 (3.07x10 2 ) 2.85x10 2 (1.90x10 1 )	4.68x10 3 (1.23x10 3 ) 6.05x10 0 (4.95x10 -1 )	1.75x10 2 (5.55x10 1 ) 4.33x10 -2 (-)
	TRANSAT-1 ARCHIMEDES-2	9 19	99 900	6.31x10 4 (1.09x10 4 ) 1.60x10 4 (4.98x10 3 )	1.19x10 2 (1.18x10 1 ) 6.75x10 2 (-)	7.08x10 3 (1.98x10 3 ) 1.84x10 1 (4.62x10 0 )	1.99x10 3 (7.19x10 2 ) n.d.
	TRANSAT-1 ARCHIMEDES-2	20 19	2500 500	7.53x10 2 (1.65x10 2 ) 2.85x10 4 (6.88x10 3 )	9.51x10 2 (8.07x10 1 ) 4.80x10 2 (2.09x10 2 )	3.47x10 2 (2.02x10 1 ) 1.41x10 1 (-)	1.97x10 1 (2.90x10 0 ) n.d.
	TRANSAT-1 ARCHIMEDES-2	20 19	1897 251	1.12x10 3 (1.25x10 2 ) 1.27x10 5 (5.84x10 4 )	3.69x10 2 (5.48x10 1 ) 1.57x10 3 (7.26x10 2 )	4.48x10 2 (1.93x10 2 ) 8.15x10 2 (1.15x10 2 )	9.98x10 0 (-) n.d.
	TRANSAT-1 ARCHIMEDES-2	20 19	1147 100	1.36x10 3 (1.71x10 2 ) 7.09x10 4 (3.77x10 4 )	2.82x10 2 (1.07x10 2 ) 1.15x10 2 (1.18x10 1 )	1.38x10 2 (4.28x10 1 ) 3.32x10 4 (1.05x10 4 )	n.d. n.d.
	TRANSAT-1 ARCHIMEDES-2	20 28	450 2501	7.34x10 3 (3.02x10 3 ) 3.31x10 4 (3.71x10 3 )	5.89x10 2 (1.20x10 2 ) 3.67x10 3 (1.46x10 3 )	3.56x10 3 (1.11x10 3 ) 5.54x10 1 (7.27x10 0 )	1.60x10 1 (-) 1.52x10 -1 (-)
	TRANSAT-1 ARCHIMEDES-2	20 28	149 1799	8.87x10 3 (2.18x10 3 ) 5.19x10 2 (9.93x10 1 )	3.78x10 2 (7.95x10 0 ) 2.87x10 3 (5.08x10 2 )	2.23x10 4 (3.90x10 3 ) 6.46x10 0 (1.20x10 0 )	5.93x10 2 (1.14x10 2 ) n.d.
	TRANSAT-1 ARCHIMEDES-2	28 28	3878 901	5.10x10 4 (1.10x10 4 ) 8.54x10 4 (2.30x10 3 )	2.72x10 3 (1.17x10 3 ) 1.12x10 3 (5.30x10 2 )	8.52x10 1 (7.78x10 0 ) 1.50x10 1 (4.22x10 0 )	8.93x10 0 (3.28x10 0 ) 1.27x10 -1 (4.54x10 -2 )
	TRANSAT-1 ARCHIMEDES-2	28 28	2999 498	5.36x10 4 (1.44x10 4 ) 3.83x10 4 (2.07x10 4	3.53x10 3 (1.41x10 3 ) 1.83x10 3 (5.90x10 2 )	1.40x10 2 (8.25x10 0 ) 1.70x10 2 (1.27x10 1 )	3.00x10 0 (4.56x10 -1 ) 4.36x10 -2 (2.09x10 -2 )
	TRANSAT-1 ARCHIMEDES-2	28 28	1897 252	5.00x10 4 (8.02x10 3 ) 4.71x10 5 (8.41x10 4 )	3.18x10 3 (-) 2.08x10 4 (2.25x10 3 )	7.98x10 1 (1.68x10 1 ) 2.39x10 3 (4.87x10 2 )	n.d. 2.34x10 -1 (7.09x10 -2 )
	TRANSAT-1 ARCHIMEDES-2	28 28	597 100	2.13x10 5 (4.69x10 4 ) 1.28x10 4 (4.17x10 3 )	4.52x10 3 (1.53x10 3 ) 1.38x10 2 (9.85x10 0 )	5.52x10 3 (1.07x10 3 ) 1.20x10 4 (8.52x10 2 )	9.78x10 2 (1.61x10 1 ) 4.31x10 -1 (1.67x10 -1 )
	TRANSAT-1 ARCHIMEDES-2	28 33	146 2502	1.38x10 5 (4.49x10 4 ) 5.35x10 4 (1.36x10 4 )	1.45x10 3 (6.67x10 2 ) 2.90x10 3 (6.65x10 2 )	9.52x10 3 (9.83x10 2 ) 6.47x10 1 (1.88x10 1 )	1.33x10 2 (2.70x10 1 ) n.d.
	ARCHIMEDES-2 ARCHIMEDES-2	3 33	2750 1800	2.30x10 4 (2.03x10 3 ) 2.70x10 4 (2.12x10 3 )	1.67x10 2 (1.63x10 1 ) 1.53x10 3 (2.92x10 2 )	7.80x10 1 (2.96x10 1 ) 3.53x10 1 (3.69x10 0 )	5.23x10 -2 (1.54x10 -2 ) n.d.
	ARCHIMEDES-2 ARCHIMEDES-2	3 33	800 749	2.21x10 5 (4.93x10 4 ) 4.51x10 5 (1.88x10 5 )	1.41x10 3 (5.94x10 1 ) 3.61x10 3 (2.02x10 2 )	7.98x10 1 (2.18x10 1 ) 1.29x10 4 (1.91x10 3 )	n.d. 3.95x10 -2 (-)
	ARCHIMEDES-2 ARCHIMEDES-2	3 33	499 498	2.49x10 4 (3.76x10 3 ) 3.24x10 5 (4.54x10 4 )	1.10x10 3 (2.62x10 2 ) 1.27x10 3 (2.9 x10 2 )	1.38x10 2 (1.95x10 1 ) 1.39x10 2 (3.09x10 1 )	n.d. 3.89x10 -2 (-)
	ARCHIMEDES-2 ARCHIMEDES-2	3 33	249 250	1.24x10 4 (9.64x10 2 ) 3.36x10 5 (9.88x10 4 )	4.69x10 3 (9.04x10 2 ) 4.84x10 3 (1.04x10 3 )	7.48x10 2 (2.51x10 2 ) 2.68x10 3 (1.24x10 2 )	n.d. 9.82x10 -2 (4.36x10 -2 )
	ARCHIMEDES-2 ARCHIMEDES-2	3 33	99 99	2.60x10 3 (4.98x10 2 ) 3.10x10 4 (3.67x10 3 )	2.30x10 1 (1.77x10 0 ) 1.56x10 3 (3.36x10 2 )	2.49x10 3 (8.78x10 1 ) 2.50x10 4 (1.20x10 3 )	8.85x10 -2 (4.53x10 -2 ) 3.19x10 -2 (-)
	ARCHIMEDES-2 ARCHIMEDES-2	5 48	4001 2500	1.54x10 4 (2.31x10 3 ) 1.69x10 3 (1.18x10 2 )	5.79x10 1 (-) 2.63x10 2 (8.87x10 1 )	2.04x10 2 (4.90x10 1 ) 6.65x10 0 (7.07x10 -2 )	6.58x10 -1 (5.23x10 -2 ) 1.62x10 -1 (-)
	ARCHIMEDES-2 ARCHIMEDES-2	5 48	2748 1798	4.01x10 3 (4.62x10 2 ) 7.07x10 3 (1.90x10 3 )	6.10x10 1 (9.82x10 0 ) 1.27x10 2 (4.36x10 1 )	1.51x10 1 (2.29x10 0 ) 6.78x10 0 (1.31x10 0 )	6.02x10 -2 (1.10x10 -2 ) n.d.
	ARCHIMEDES-2 ARCHIMEDES-2	5 48	899 901	1.39x10 4 (2.00x10 3 ) 5.04x10 3 (2.14x10 3 )	1.32x10 2 (7.71x10 1 ) 7.18x10 1 (1.06x10 0 )	7.68x10 1 (2.09x10 1 ) 1.73x10 0 (1.16x10 0 )	9.70x10 -1 (7.57x10 -1 ) 1.96x10 -1 (-)
	ARCHIMEDES-2 ARCHIMEDES-2	5 48	501 499	7.16x10 4 (1.68x10 4 ) 2.43x10 5 (7.20x10 4 )	2.20x10 2 (1.20x10 2 ) 7.48x10 2 (3.10x10 2 )	1.59x10 3 (2.76x10 2 ) 5.33x10 2 (1.08x10 2 )	9.28x10 -1 (2.47x10 -2 ) 1.81x10 -1 (4.08x10 -2 )
	ARCHIMEDES-2 ARCHIMEDES-2	5 48	249 249	8.37x10 5 (2.29x10 5 ) 1.83x10 4 (5.20x10 3 )	6.16x10 2 (7.07x10 1 ) 5.93x10 2 (1.88x10 2 )	4.73x10 3 (4.68x10 2 ) 3.24x10 2 (8.53x10 1 )	2.31x10 0 (5.90x10 -1 ) 1.21x10 -1 (-)
	ARCHIMEDES-2 s.d.: standard deviation ;	5	100	1.03x10 5 (3.88x10 4 )	4.85x10 2 (1.54x10 2 )	2.32x10 4 (2.62x10 3 )	5.94x10 1 (6.55x10 0 )
							11

Table S4 -Observed and estimated richness in archaeal amoA gene libraries

 S4

				No. of	No. of
	Station	Depth	Water Mass	clones	OTUs	Chao1 (2%)
				sequenced	(2%)
	TRANSAT-1, station 28	597 m	NACW	45	15	20
	ARCHIMEDES-2, station 33	749 m	AAIW	39	8	18
	ARCHIMEDES-2, station 33	2502 m	MNEADW	39	7	8.5

Acknowledgements

We thank the captain and crew of the RV Pelagia for their help during work at sea. We thank H.M. van Aken for the characterization of the water masses, J.M. Arrieta for collecting the samples during the TRANSAT-1 cruise, B. Abbas, J. van Bleijswijk and H. Witte for technical discussions about the quantitative PCR and phylogenetic analyses and T. Reinthaler for helpful discussions. We also thank A. Hunting, A. M. Schmitz and W. D. Lienhart for help with data processing. This research was supported by a Marie Curie Fellowship of the European Community to H.A. and shiptime was provided through grants of the Earth and Life Science Division of the Dutch Science Foundation (ALW-NWO) (TRANSAT and ARCHIMEDES projects) to G.J.H. The work was carried out within the frame of the 'Networks of Excellence' MarBef and EurOceans supported by the 6th Framework Program of the European Union.

Supplementary information is linked to the online version of the paper at www.nature.com/nature.

Authors contributions

Author information

The archaeal sequences reported in this study have been deposited at GenBank (NCBI) under accession numbers EU650236-EU650270 (Station 3, ARCHIMEDES-2), FJ002858-FJ002876 (Station 23, ARCHIMEDES-3), FJ150794-FJ150834 (Station 9, TRANSAT-1) for 16S rRNA genes and EU795424-EU795460 and EU810209-810235 for amoA genes. Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interest. Correspondence and requests for materials should be addressed to G.J.H.

(herndl@nioz.nl).

Supplementary information

Study sites and sampling.

During the TRANSAT-1 cruise (September 2002), in total 36 stations were occupied with RV Pelagia following a track from 65°N to 35°N in the eastern basin of the North Atlantic.

Quantitative PCR analyses were done on 3 stations: T9, T20 and T28 (Supplementary Fig. 1).

Sampling was performed with 22 10-L NOEX (no oxygen exchange) bottles mounted in a frame also holding sensors for conductivity, temperature, depth, chlorophyll fluorescence, and optical backscattering. The specific water masses were identified based on their temperature and salinity characteristics and their oxygen concentrations (Supplementary Table 1). One subsurface sample per station was taken from the lower euphotic layer at around 100 m depth.

The waters below the surface thermocline (500 m depth) were formed by North Atlantic Central Water (NACW). Below, the 1000-2000 m depth layer was occupied by either Labrador Seawater (LSW; Sts T20, T28) or Antarctic Intermediate Water (AAIW, St T9). At Station 9, below the AAIW, slope water was present. At around 2000 m depth, the Iceland-Scotland Overflow Water (ISOW) was sampled at Stations T9 and T20. At Station T28, Middle Northeast Atlantic Deep Water (MNEADW) and Lower Northeast Atlantic Deep Water (LNEADW) were sampled at around 3000 m and 4000 m depth, respectively. For each depth, one liter of seawater was filtered onto 0.2 µm polycarbonate filter (Millipore) and the filters were subsequently stored at -80°C until further processing in the lab.

During the ARCHIMEDES-2 cruise (November to December 2006), a total of 20 stations was occupied with RV Pelagia following a track from 26°N to 5°S in the eastern (sub)tropical Atlantic. Quantitative PCR analyses were done at 7 stations: A3, A5, A11, A19, A28, A33 and A48 (Supplementary Fig. 1). Sampling was performed as described above for TRANSAT-1 and according to the water masses (Supplementary Table 1). One subsurface sample per station was taken from the lower euphotic layer at around 100 m depth. The waters