Nadia El Mrabet
email: nadia.el-mrabet@emse.fr

Loubna -3840-584x]

Ghammam

Nicolas Meloni
email: nicola.meloni@univ-tln.fr

Emmanuel Fouotsa
email: fouotsa.emmanuel@uniba.cm

New Versions of Miller-loop Secured against Side-Channel Attacks

Keywords: Miller algorithm, Euclidean addition chains, Fibonacci number, optimal Ate pairing, sidechannel attacks

In this paper, we propose two new versions of Miller algorithm in order to secure pairing computations against existing side-channel attacks (SCA). We have chosen to use the co-Z arithmetic on elliptic curves from which we derive two methods for pairing computations: one based on Euclidean addition chains and one based on Zeckendorf representation. We show that our propositions are resistant to existing side-channel attacks against pairing-based cryptography. We consider differential power analysis and fault attacks. The complexities of our solutions are compared with state-of-the-art one. We demonstrate that our new proposed versions are more efficient by 17%.

Introduction

Pairing-based cryptography (PBC) provides several protocols such as short signature protocols and hierarchical encryption [1,[START_REF] Shamir | Identity-based cryptosystems and signature schemes[END_REF], making it a promising tool for the Internet of things (IoT) or cloud computing. In the first case, pairings would be implemented on constrained devices [START_REF] Duquesne | Memory-saving computation of the pairing final exponentiation on BN curves[END_REF] such as microcontrollers and would be subject to invasive and non-invasive Side-Channel Attacks (SCA) [START_REF] Mrabet | A practical differential power analysis attack against the miller algorithm[END_REF][START_REF] Blömer | Improved side channel attacks on pairing based cryptography[END_REF][START_REF] Unterluggauer | Practical attack on bilinear pairings to disclose the secrets of embedded devices[END_REF][START_REF] Mrabet | A survey of fault attacks in pairing based cryptography[END_REF][START_REF] Mrabet | Guide to Pairing-Based Cryptography[END_REF][START_REF] Jauvart | Improving side-channel attacks against pairing-based cryptography[END_REF]. Efficient and secure pairing computation has thus been an active research area [START_REF] Chatterjee | On the Efficiency and Security of Pairing-Based Protocols in the Ttype 1 and Ttype 4 Settings[END_REF][START_REF] Mrabet | Guide to Pairing-Based Cryptography[END_REF][START_REF] Azarderakhsh | Fast software implementations of bilinear pairings[END_REF].

The most efficient pairing computation algorithms are usually based on the Tate model [START_REF] Frey | The tate pairing and the discrete logarithm applied to elliptic curve cryptosystems[END_REF] and rely on the Miller algorithm in order to compute the rational function f s,Q such that Div(f s,Q ) = s(Q)-([s]Q)-(s-1)(O), where P and Q are two points of an elliptic curve defined over a finite field F q and O represents the identity element for the group of rational points of the elliptic curve. This function is computed thanks to the doubleand-add loop called Miller loop [START_REF] Victor | The weil pairing, and its efficient calculation[END_REF] and followed by a final exponentiation.

Such algorithms are known to be natural targets to SCA ( [START_REF] Page | Fault and side channel attacks on pairing based cryptography[END_REF][START_REF] Whelan | Side Channel Analysis of Practical Pairing Implementations: Which Path Is More Secure?[END_REF][START_REF] Unterluggauer | Practical attack on bilinear pairings to disclose the secrets of embedded devices[END_REF][START_REF] Jauvart | Improving side-channel attacks against pairing-based cryptography[END_REF][START_REF] Mrabet | Physical Attacks, chapter 12[END_REF][START_REF] Jauvart | Improving side-channel attacks against pairing-based cryptography[END_REF]). Our contribution is to provide a new approach to Miller loop to replace the usual double-and-add loop by a loop based on addition chains. Such approaches, e.g., based on the double-and-add algorithm, have been proven to be an interesting alternative to standard scalar multiplication methods when security is in balance with speed, especially when combining to co-Z arithmetic [START_REF] Meloni | New point addition formulae for ECC applications[END_REF]. Concretely, we provide two versions of the Miller algorithm, one based on Euclidean addition chain and the other version based on the Zeckendorf representation of the parameter controlling the number of iterations. We show that both versions are resistant to Side-Chanel Attacks. Though the computational cost of these versions of the Miller algorithm is higher than that of the classic Miller algorithm, they are the most efficient among other countermeasures like masking against SCA.

The paper is organized as follows. Section 2 and 3 are devoted to the necessary mathematical background on pairing computations, co-Z arithmetic, and scalar multiplication using Euclidean additions chains (EAC). In Section 4, we present the new versions of the Miller algorithm, their complexities, and we compare them with the standard version. Section 6 is devoted to the state of the art on SCA and in Section 7 we analyze the security of our algorithms against some of those attacks. Finally, in Section 8, we provide detailed comparisons between our counter-measure, and other counter-measures used in the literature.

Notations: In this paper we denote by:

-M e a multiplication in the field F p e .

-S e a squaring in the field F p e .

A multiplication, a squaring and an inversion in F p are respectively denoted by M , S and I.

Background on pairing computations

Let E be an elliptic curve defined over a prime field F p , with p a large prime number. Let r be a large prime divisor of #E(F p ). In practice, (E, p, r) are provided using parametric families [START_REF] Freeman | A taxonomy of pairing-friendly elliptic curves[END_REF]. Let k be the smallest integer such that r divides p k -1, k is called the embedding degree of E relative to r.

Let G 1 = E(F p )[r] be the r-torsion subgroup of E(F p ). Let G 2 = E (F p k/d )[r]∩ Ker(π p -[p]
) where E is the twist of E (if it exists) of degree d, π p represents the Frobenius map over E and [p] is the scalar multiplication by p over E. The subgroup of F p k consisting of r-th roots of unity is denoted by G 3 = µ r . Let s be an integer constructed for the optimal Ate pairing [START_REF] Vercauteren | Optimal pairings[END_REF] and depending on r and p. Then the (optimal) Ate pairing [START_REF] Vercauteren | Optimal pairings[END_REF] is given by

T r : G 1 × G 2 → G 3 ; (P, Q) → f s,Q (P ) (p k -1)/r .
The pairing computation is divided into two main steps [START_REF] Mrabet | Guide to Pairing-Based Cryptography[END_REF]Chap. 3]. First, one has to compute f s,Q (P ) using an iterative algorithm denoted as the Miller algorithm and described in Algorithm 1 below. As any algorithm based on the double-and-add algorithm, the overall cost directly depends on the length and Hamming weight of integer s.

Algorithm 1 : Miller algorithm (P, Q) → f s,Q (P )

Require: P ∈ G1, Q ∈ G2, s = (sn-1, . . . , . . . s0): binary representation of s Ensure: fs,Q(P ) ∈ G3(⊂ F * p k ) 1: f1 ← 1 2: T ← Q 3: for i = n -2 down to 0 do 4: T ← [2]T 5: f1 ← f 2 1 • l1(P )
where l1 is the tangent to the point T .

6:

if ri = 1 then 7:

T ← T + Q 8: f1 ← f1 • l2(P )
where l2 is the line passing through T and P . 9:

end if 10: end for 11: return f1

The second step of computing the Tate pairing and its variants is the final exponentiation and consists of raising the final result of the main loop, f s,Q (P ), to the power of p k -1 r . The computation of this part can be simplified thanks to the k-th cyclotomic polynomial [START_REF] Duquesne | Memory-saving computation of the pairing final exponentiation on BN curves[END_REF].

Co-Z arithmetic on elliptic curves

Given an elliptic curve E over a finite field F q there are many ways to compute the addition of two points. To avoid the computation of inversions in F q when considering affine coordinates, it is standard to represent the points of E in Jacobian coordinates. A point (x, y) given in affine coordinates is represented in Jacobian coordinates by a triplet (X : Y : Z) such that x = X/Z 2 and y = Y /Z 3 , the point at infinity being then represented by (0 : 1 : 0). A typical point addition costs 11 multiplications (M) and 5 squarings (S) in F q . Considerable literature now exists on the various point addition formulae, coordinate systems and curve shapes that can be used in a different context [START_REF] Cohen | Handbook of Elliptic and Hyperelliptic Curve Cryptography[END_REF].

In the particular case of two points sharing the same Z-coordinate, the addition can be performed using only 5M+2S [START_REF] Meloni | New point addition formulae for ECC applications[END_REF]. This operation is usually referred to as co-Z addition or ZADD. Finding two points sharing the same Z-coordinate is very unlikely but, when combined with the right scheme, one can perform a whole scalar multiplication using the co-Z addition. In this work, we considered two of those schemes: one based on Euclidean addition chains (EAC) and another based on the Zeckendorf representation of integers [START_REF] Meloni | New point addition formulae for ECC applications[END_REF].

Euclidean addition chains

A Euclidean addition chain (EAC) computing an integer k is a sequence n = (n 1 , n 2 , . . . , n w ) such that n 1 = 1, n 2 = 2, n 3 = 3, n w = k and ∀ 3 i w -1, if n i = n i-1 + n j ,
for some j < i -1, then we have n i+1 = n i + n i-1 (big step) or n i+1 = n i + n j (small step). Such an addition chain, computing an integer greater than 3, can easily be represented by a binary chain c 4 , c 5 , . . . , c w representing the succession of big steps (noted 0) and small steps (noted 1) starting from (1, 2, 3). For instance, the EAC (1, 2, 3, 4, 7, 10, 13, 23) can be represented as (1, 0, 1, 1, 0). Another example is the Fibonacci sequence (1, 2, 3, 5, 8, 13, 21, . . . ) which is only made of big steps.

Example 31 The addition chain (1, [START_REF] Shamir | Identity-based cryptosystems and signature schemes[END_REF][START_REF] Duquesne | Memory-saving computation of the pairing final exponentiation on BN curves[END_REF][START_REF] Mrabet | A practical differential power analysis attack against the miller algorithm[END_REF][START_REF] Mrabet | A survey of fault attacks in pairing based cryptography[END_REF][START_REF] Azarderakhsh | Fast software implementations of bilinear pairings[END_REF][START_REF] Whelan | Side Channel Analysis of Practical Pairing Implementations: Which Path Is More Secure?[END_REF][START_REF] Meloni | New point addition formulae for ECC applications[END_REF]34,53) is an Euclidean addition chain computing 53. For example, in step 4, we have computed 4 = 3 + 1. Then for step 5, we must add 3 or 1 to 4, that means that from step 4, we can compute 7 = 4 + 3 or 5 = 4 + 1. In our example, we have chosen to compute, 7 = 4 + 3 (which is a big step). For step 6, we have computed 11 = 7 + 4 (big step), etc.

Finding such chains is simple: choose an integer g co-prime with k and then apply the subtractive form of Euclid's algorithm to those numbers.

-3 = 1 3 -1 = 2 2 -1 = 1 1 -1 = 0
The Euclidean addition chain computing 53 is (1, 2, [START_REF] Duquesne | Memory-saving computation of the pairing final exponentiation on BN curves[END_REF][START_REF] Mrabet | A practical differential power analysis attack against the miller algorithm[END_REF][START_REF] Mrabet | A survey of fault attacks in pairing based cryptography[END_REF][START_REF] Azarderakhsh | Fast software implementations of bilinear pairings[END_REF][START_REF] Whelan | Side Channel Analysis of Practical Pairing Implementations: Which Path Is More Secure?[END_REF][START_REF] Meloni | New point addition formulae for ECC applications[END_REF]34,53). We find this chain by reading the first number of each line of the subtractive form of Euclid's algorithm.

It has been shown that those chains can be an efficient alternative to the classical binary chain used in most scalar multiplication algorithms when side-channel attack resistance is an issue [START_REF] Meloni | New point addition formulae for ECC applications[END_REF][START_REF] Herbaut | Random euclidean addition chain generation and its application to point multiplication[END_REF]. Algorithm 2 [START_REF] Meloni | New point addition formulae for ECC applications[END_REF] describes how to perform a scalar multiplication using an EAC.

The main drawback of this approach is that, although it is easy to find an EAC computing a given integer k, finding a short one seems to be a challenging problem [START_REF] Knuth | Analysis of the subtractive algorithm for greater common divisors[END_REF] making it difficult to use EAC for scalar multiplication directly. However, in the case of PBC, the Miller algorithm uses a fixed scalar which means that it is worth investigating finding an efficient EAC computing the parameter s. However, this method will become less efficient when s grows in size as finding a short EAC will become harder.

Algorithm 2 [START_REF] Meloni | New point addition formulae for ECC applications[END_REF] Calculation of [k]P where k is presented as an Euclidean Additions Chain.

Require: P ∈ E, c = (c4, . . . , . . . cw) an EAC computing k Ensure: [k]P ∈ E (U1, U2) ← ([2]P, P ) for i = 4 to w do if ci = 0 then (U1, U2) ← (U1 + U2, U1) else (U1, U2) ← (U1 + U2, U2) end if end for U1 ← U1 + U2 return U1

Zeckendorf Representation

Let k be an integer and (F i ) i≥0 the Fibonacci sequence, a classical result states that k can be uniquely written in the form k [START_REF] Zeckendof | Représentations des nombres naturels par une somme de nombre de Fibonacci ou de nombres de Lucas[END_REF]. An integer k written in this form is said to be in Zeckendorf representation and will be denoted as k = (d l-1 , . . . , d 2 ) Z . Such a representation is easy to compute as it can be obtained using a greedy algorithm [START_REF] Meloni | New point addition formulae for ECC applications[END_REF].

= l i=2 d i F i , with d i ∈ {0, 1} and d i d i+1 = 0
As mentioned in [START_REF] Zeckendof | Représentations des nombres naturels par une somme de nombre de Fibonacci ou de nombres de Lucas[END_REF] and [START_REF] Meloni | New point addition formulae for ECC applications[END_REF], the disadvantage of the Zeckendorf representation is that it requires 44% more digits than the binary method requires. However, contrary to EAC, it is easy to bound its length [START_REF] Zeckendof | Représentations des nombres naturels par une somme de nombre de Fibonacci ou de nombres de Lucas[END_REF]. Given the Zeckendorf representation of k, one can perform a scalar multiplication using the following Algorithm 3 [START_REF] Meloni | New point addition formulae for ECC applications[END_REF].

Algorithm 3 [START_REF] Meloni | New point addition formulae for ECC applications[END_REF] Computing [k]P using Zeckendorf representation

Require: P ∈ E, k = (dn, . . . , . . . d2)Z Zeckendorf representation of k with dn = 1 Ensure: [k]P ∈ E 1: (U1, U2) ← (P, P ) 2: for i = n -1 to 2 do 3: if di = 1 then 4:
(U1, U2) ← (U1 + P, U2) 5:

end if 6:

(U1, U2) ← (U1 + U2, U1) 7: end for 8: return U1

Co-Z approach to the Miller Algorithm

In this section, we will present two variants of the Miller algorithm. We should remember that the basic idea of Miller's algorithm is to consider the binary representation of the Miller loop length parameter s and apply the double-and-add algorithm. The main idea of our work is to consider different representations of s to take advantage of co-Z formulae adapted for pairing computations.

Miller-Euclide

The idea is to adapt the scalar multiplication scheme described in Algorithm 2 to Miller loop and take advantage of the co-Z coordinates to limit the computational cost. Our approach, referred to as Miller-Euclide, is shown in Algorithm 4.

Algorithm 4 : Miller-Euclide, The computation of f s,Q (P ) using an Euclidean addition chain.

Require: P ∈ G1, Q ∈ G2, c = (c4, . . . , . . . cw) the euclidean addition chain computing s, Ensure: fs,Q(P ) 1: (T1, T2) ← ([2]Q, Q) 2: (f1, f2) ← (lQ,Q(P ), 1) 3: for i = 4 to w do 4: if ci = 0 then 5: (f1, f2) ← (f1 × f2 × T 1 ,T 2 (P ), f1) 6: (T1, T2) ← (T1 + T2, T1) 7: else 8: (f1, f2) ← (f1 × f2 × T 1 ,T 2 (P ), f2) 9:
(T1, T2) ← (T1 + T2, T2) 10:

end if 11: end for 12: f1 ← f1 × f2 × T 1 +T 2 ,T 1 (P ) 13: return f1
First, let us notice that the initialization consists of one doubling and the evaluation of the tangent line passing through Q. It has been shown [START_REF] Meloni | New point addition formulae for ECC applications[END_REF] that the doubling is compatible with the co-Z coordinates without additional cost so that the initialization can be performed in 2 multiplications and 5 squarings for the doubling. We add to this cost 3 multiplications and one squaring for the evaluation of the tangent. The details of this computation are presented in Algorithm 5.

The main loop of the new Miller-Euclide algorithm is only made of addition steps, one for each bit of the EAC, so that the computational cost of our method is simply linked to the length of the EAC representing s.

Algorithm 5 Initialization step of Algorithm 4 lines 1 and 2

Require: P = (xP , yP ) and Q = (XQ, YQ, ZQ)

Ensure: T1 = [2]Q, T2 = Q sharing the same z-coordinate Z , f1 = Q,Q(P ), XP = xP Z 2 and YP = yP Z 3 1: A ← X 2 Q (S k/d ) 2: B ← Y 2 Q (S k/d ) 3: C ← B 2 (S k/d ) 4: D ← 2((X1 + B) 2 -A -C) (S k/d ) 5: E ← 3A 6: F ← E 2 (S k/d ) 7: X2Q ← F -2D 8: Y2Q ← E(D -X2Q) -8C (M k/d ) 9: Z2Q ← 2YQZQ (M k/d ) 10: XP ← xP Z 2 2Q (S k/d + M k/d ) 11: YP ← yP Z 3 2Q (2M k/d ) 12: f1, f2 ← Yp + Xp -3A 2 + 16C, 1 13: T 2 = (X2, Y2) ← D, 8C 14: T 1 = (X1, Y1) ← X2Q, Y2Q
In order to evaluate the cost of our method we must evaluate the cost of a full addition step. Let T 1 = (X 1 , Y 1 , Z) and T 2 = (X 2 , Y 2 , Z) be two points with the same z-coordinate. From [START_REF] Meloni | New point addition formulae for ECC applications[END_REF] we have that the coordinates (X 3 , Y 3 , Z 3 ) of the point T 3 = T 1 + T 2 can be computed using the following equations:

A = (X 2 -X 1 ) 2 , B = X 1 A, C = X 2 A, D = (Y 2 -Y 1 ) 2 X 3 = D -B -C , Y 3 = (Y 2 -Y 1 )(B -X 3 ) -Y 1 (C -B) , Z 3 = Z(X 2 -X 1 ) .
This computation requires 5 multiplications and 2 squarings over F p k/d . On top of that the quantities 3 computed during the addition can be seen as the x and y-coordinates of the point (

X 1 A = X 1 (X 2 -X 1 ) 2 and Y 1 (C -B) = Y 1 (X 2 -X 1 )
X 1 (X 2 -X 1 ) 2 , Y 1 (X 2 -X 1 ) 3 , Z(X 2 -X 1 )) ∼ (X 1 , Y 1 , Z).
Thus it is possible to add P 1 and P 1 + P 2 with the same formulae during the next iteration of the main loop.

Let us now consider the value of l T1,T2 (P ) given by the equation of the line joining T 1 and T 2 evaluated on P = (x P , y P ). The standard equation of such line in affine coordinates is given by

y p = y 1 -y 2 x 1 -x 2 x P + y 2 x 1 -y 1 x 2 x 1 -x 2 , (1)
To convert this equation to co-Z coordinate, we just need to apply the transformation

(x i , y i ) → Xi Z 2 , Yi Z 3
which gives us

y p = Y 1 -Y 2 Z(X 1 -X 2 ) x P + Y 2 X 1 -Y 1 X 2 Z 3 (X 1 -X 2 ) ⇔ y p (X 1 -X 2 )Z 3 = (Y 1 -Y 2 )Z 2 x p + Y 2 X 1 -Y 1 X 2 ⇔ y p (X 1 -X 2 )Z 3 = (Y 1 -Y 2 )(x p Z 2 -X 1 ) -Y 1 (X 1 -X 2 ) ⇔ y p Z 3 = (Y 1 -Y 2 )(x p Z 2 -X 1 ) -Y 1 where X 1 = X 1 (X 1 -X 2 ) 2 , Y 1 = Y 1 (X 1 -X 2 ) 3 and Z = Z(X 1 -X 2 )
. Now we remark that all those values have been computed during the point addition as well as the values of (X 1 -X 2 ) 2 and (X 1 -X 2 ) 3 . If we suppose that the values of y p Z 3 and x p Z 2 were stored from the previous iteration of the main loop, the evaluation of the line joining the points T 1 and T 2 required exactly 3 multiplications over F p k/d . Moreover, the value of Z = Z(X 1 -X 2 ) itself is not needed during the whole process, so that we can spare one multiplication per iteration, for a total cost of 7 multiplications and 2 squarings over F p k/d . Finally, we have to compute f 1 × f 2 × T1,T2 (P ) which requires one sparse multiplication between elements of F p k/d and F p k and one full multiplication over F p k .

Algorithm 6 sums up the computations required to perform one addition step with their respective costs. Therefore, each addition step of Algorithm 4 requires

M k + M k,k/d + 7M k/d + 2S k/d operations. At the end the total cost of the algorithm is 5M k/d + 6S k/d + (n -2)(M k + M k,k/d + 7M k/d + 2S k/d ).

Algorithm 6 Computing one addition step of Algorithm 4

Require: T1, T2 with same z-coordinate Z, Xp = xpZ 2 , Yp = yP Z 3 and f1, f2

∈ F p k Ensure: Z = Z(X1 -X2), T3 = T1 + T2, T1 = (X3, Y3, Z ), T2 = (X1Z 2 , Y1Z 3 , Z ), Xp = xpZ 2 , Yp = yP Z 3 , f1 = f1 × f2 × T 1 ,T 2 (P ) 1: A ← (X1 -X2) 2 (S k/d ) 2: B ← X1A (M k/d ) 3: C ← X2A (M k/d ) 4: D ← (Y1 -Y2) 2 (S k/d ) 5: E ← C -B 6: F ← Y1E 7: X3 ← D -B -C 8: Y3 ← (Y2 -Y1)(B -X3) -F (2M k/d ) 9: XP ← XP A (M k/d ) 10: YP ← YP E (M k/d ) 11: L ← YP -(Y1 -Y2)(XP -X3) -Y3 (M k/d ) 12: f3 ← f1L (M k,k/d ) 13: f3 ← f1f2 (M k ) 14: f1, f2 ← f3, f1 15: X2, Y2 ← B, F 16: X1, Y1 ← X3, Y3

Miller-Fibonacci

Let us now adapt Algorithm 3 to the Miller loop to present our new Algorithm 7 which computes f s,Q (P ) using the Zeckendorf representation of s. (f1, f2) ← (f1 × lT 1 ,T 2 (P ), f2) 6:

(T1, T2) ← (T1 + Q, T2) 7:

end if 8:

(f1, f2) ← (f1 × f2 × lT 1 ,T 2 (P ), f1) 9:

(T1, T2) ← (T1 + T2, T1) 10: end for 11: return f1

First, let us note that the first computation will be a doubling instead of an addition, independently of the value of s n-1 , but for the sake of simplicity, we only use the additive notation. For this particular step, we use Algorithm 5.

Each step of this algorithm is similar to Miller-Euclide. If the current digit is a 0, we use the procedure described in Algorithm 6. If the current digit is a 1 then we have to perform an additional sparse multiplication and then the evaluation of the line passing through two points of the elliptic curve. To do so, we must first upgrade the coordinates of T 1 and Q so that they share the same z-coordinate as well as the values of X p and Y p (this is easily done in 7 multiplications and 2 squarings) and then use the formulae given in Algorithm 6. On top of that, we have to keep track of the z-coordinate of T 1 during the whole process which adds another multiplication. In the end the cost when the

d i = 0 is M k +M k,k/d +8M k/d +2S k/d plus 7M k/d +2S k/d if d i = 1. The total cost of Algorithm 7 is then 5M k/d +6S k/d +(n-2)(M k +M k,k/d +7M k/d +2S k/d )+hz(s)(7M k/d +2S k/d )
where hz(s) is the hamming weight of the Zeckendorf representation of s.

Comparison

In this section, we compare three versions of the Miller algorithm that allow us to compute the rational function f s,Q evaluated at P : Miller's algorithm, Miller-Euclide, and Miller-Fibonacci. For the sake of simplicity, we restraint our comparisons to the computation of the Optimal Ate pairing over BN curves [START_REF] Pereira Geovandro | A family of implementation-friendly BN elliptic curves[END_REF] and BLS 12 curves [START_REF] Paulo | Constructing elliptic curves with prescribed embedding degrees[END_REF].

After Barbulescu and Duquesne results presented in [START_REF] Barbulescu | Updating key size estimations for pairings[END_REF], it is recommended to use BLS12 for computing Optimal Ate pairing for the 128-bit security level instead of BN curves. However, BN curves are still considered in practice for several schemes [START_REF]Black-box wallets: Fast anonymous two-way payments for constrained devices[END_REF][START_REF] Canard | Blindids: Market-compliant and privacy-friendly intrusion detection system over encrypted traffic[END_REF]. That's why, in this paper, we consider the two curves in our computation.

Both of these curves have the embedding degree k = 12, therefore, the arithmetic in their extension tower has the same complexity (since we don't consider additions in our cost evaluation). The tower extension used for pairing computation on both curves is given by: The field F p 12 is built using the following extension tower.

-F p 2 = F p [α]/(α 2 + 1) -F p 6 = F p 2 [β]/ β 3 -(α + 1)) -F p 12 = F p 6 [γ]/ γ 2 -β
In Table 1, we also present the cost of each operation needed for computing Miller algorithm using projective coordinates [START_REF] Sylvain Duquesne | Choosing and generating parameters for low level pairing implementation on BN curves[END_REF] and also the cost of the operations needed for our new versions of Miller's algorithm. Note that these two curves, BN and BLS12, have a twist of degree d = 6 and the twist isomorphism is given by the following map.

Ψ : E F p 2 → E(F p 12 ) (x Q , y Q ) → (γ 2 x Q , γ 3 y Q )

BN curves

A Barreto-Naehrig (BN) curve [START_REF] Pereira Geovandro | A family of implementation-friendly BN elliptic curves[END_REF] is an elliptic curve E : y 2 = x 3 + b defined over a finite prime field F p such that b ∈ F p , -E(F p ) has prime order n = #E(F p ), p = p(u) = 36u 4 + 36u 3 + 24u 2 + 6u + 1, and n = n(u) = 36u 4 + 36u 3 + 18u 2 + 6u + 1 for some u ∈ Z.

We choose the following elliptic curve which is efficient in practice [START_REF] Sylvain Duquesne | Choosing and generating parameters for low level pairing implementation on BN curves[END_REF]:

E(F) : y 2 = x 3 + 2.
The twisted elliptic curve is defined by the equation:

E (F p 2 ) : y 2 = x 3 + (1 -i).
Recall that the Optimal Ate pairing over BN curves is the following map:

E(F p )[r] × Ψ (E (F p 2 ))[r] -→ F p 12 (P, Q) -→ = (f 6u+2,Q (P )l [6u+2]Q,πp(Q) (P )l [6u+2]Q,π 2 p (Q) (P )) p 12 -1 r
Using the costs presented in Table 1, we obtain the comparison between the different versions of Miller's algorithm to compute f 6u+2,Q (P ), in terms of operations in F p . This comparison is presented in Table 2. We recall that the pairing parameter in our case is s = 6u + 2 with u = -2 114 + 2 101 -2 14 -1. In this case, it is possible to find a EAC of length 168 (using the subtractive euclidean algorithm with parameter s with s = 174601120802286487002979030408532519 when we consider u = 47085091321987640844053613709579372). Moreover, the Zeckendorf representation of s = 6u + 2 = F ibo(169) + F ibo(61) + F ibo(99) with a chosen parameter u = 15533701296897238962071235758772159, has a length of 168 and its hamming weight is 3. Thus, we have: From this Table 2, we can notice that the classical Miller algorithm is more efficient than our two new versions. However, it is not secure when considering DPA attack. We show in the next section that our method protects against this attack.

BLS 12 Curves

In 2002, Barreto, Lynn and Scott presented in [START_REF] Paulo | Constructing elliptic curves with prescribed embedding degrees[END_REF] a method to generate pairing-friendly elliptic curves over a prime field F p with embedding degree k = 12. BLS12 are defined over F p by the following equation:

E : y 2 = x 3 + b
and by a parameter u ∈ Z such that:

   p = (u -1) 2 (u 4 -u 2 + 1)/3 + u r = u 4 -u 2 + 1 t = u + 1 (2)
where t is the trace of the Frobenius map on the curve. The parameter u is chosen such that p and r are prime and have the sizes corresponding to the desired security level. For the 128-bit security level and as recommended in [START_REF] Barbulescu | Updating key size estimations for pairings[END_REF], p and r are of at least 461 and 308 and the proposed parameter u = -2 77 + 2 50 + 2 33 . The choosed BLS12 elliptic curve is defined over F by E(F) : y 2 = x 3 + 4 which admits a twist of degree d = 6 and given by E (F p 2 ) :

y 2 = x 3 + 4(1 -i).
The Optimal Ate pairing on BLS12 curves is defined by the following map.

E(F p )[r] × Ψ (E (F p 2 ))[r] -→ F p 12 (P, Q) -→ = ((f u,Q (P )) p 12 -1 r .
Using the costs presented in Table 1, we obtain the comparison between the different versions of the Miller algorithm to compute f u,Q (P ), in terms of operations in F p . This comparison is presented in Table 3. We recall that the pairing parameter in our case is u = -2 77 + 2 50 + 2 30 . In this case, it is possible to find a EAC of length 112 (using the subtractive euclidean algorithm with parameter u = 204035480723636378792533 and s = 126100861998131272870737). Note that u is chosen such that p and r are primes. Moreover, the Zeckendorf representation of u = 70492524767089125860497 = F ibo(3) + F ibo(7) + F ibo(24) + F ibo(111) has length 111 and its hamming weight is 4 (u is chosen such that we obtain an efficient Zeckendorf representation where p and r are primes. Thus, we have:

From this Table 3, we can notice that the classical Miller algorithm is more efficient than our two new versions. However, it is not secure when considering DPA attack. We show in the next section that our method protects against this attack. 6 SCA against Miller's algorithm

The SCA re the perturbation analysis that are invasive and observation analysis that are non-invasive. Noninvasive attacks against pairing, such as DPA, CPA, template attack, have been widely studied [START_REF] Page | Fault and side channel attacks on pairing based cryptography[END_REF][START_REF] Whelan | Side Channel Analysis of Practical Pairing Implementations: Which Path Is More Secure?[END_REF][START_REF] Mrabet | A practical differential power analysis attack against the miller algorithm[END_REF][START_REF] Ghosh | Security of prime field pairing cryptoprocessor against differential power attack[END_REF][START_REF] Unterluggauer | Practical attack on bilinear pairings to disclose the secrets of embedded devices[END_REF][START_REF] Pan | A correlation power analysis attack against Tate pairing on FPGA[END_REF]34,[START_REF] Jauvart | Improving side-channel attacks against pairing-based cryptography[END_REF][START_REF] Jauvart | Improving side-channel attacks against pairing-based cryptography[END_REF], so as the invasive attacks (basically fault attacks) [START_REF] Mrabet | A survey of fault attacks in pairing based cryptography[END_REF][START_REF] Mrabet | Guide to Pairing-Based Cryptography[END_REF]. In this section, we briefly recall the non-invasive attacks against pairing, then, we will illustrate the fact that our new algorithms seem to be resistant to existing attacks.

Non-invasive attacks

Non-invasive attacks can use for instance the power consumption or electromagnetic emission of a device. Those attacks are also called observation analysis, they include Differential Power Analysis, Correlation Power Analysis, Template attacks and can be performed in vertical or horizontal mode. We will describe the resistance of our algorithms by using the DPA attack model to resume the non invasive attacks. Those attacks can allow an attacker to compute the intermediate values within cryptographic computations through statistical analysis of data collected from multiple cryptographic operations. In pairing-based cryptography, the secret is one of the inputs, either the point P either the point Q. The public data is the other point. The vulnerabilities of a pairing implemented with a classical Miller algorithm have been illustrated in several articles [35,36,

. The scheme is as follow:

1. First choose the secret (P or Q).

2. Then find during the algorithm an operation that involves both the secret data and a data that can be manipulated freely during the attack. This operation is the target of the attack. In pairing-based cryptography, the freely chosen data are the coordinates of the public point. 3. For a vertical attack: the operation is executed at a specific time t 0 , the statistical analysis is performed on a collection of traces for the same secret point and hundred to millions of public point. 4. For a horizontal attack: the same secret data is used more than once during a single algorithm execution, the statistical analysis is performed using one public entry at different times t i .

Existing Countermeasure

To avoid SCA in pairing-based cryptography, several counter-measures have been proposed in the literature:

1. The counter-measure proposed by Öztürk et al. in [START_REF] Öztürk | Tate pairing with strong fault resiliency[END_REF]. The principle of this counter-measure is to avoid any perturbation against Miller's algorithm using resilient counters. The disadvantage of this counter-measure is that it is a proof of concept implementation and doesn't have a theoretical basis for measuring its security. 2. The second counter-measure is proposed by Gosh et al. in [START_REF] Ghosh | Fault attack, countermeasures on pairing based cryptography[END_REF]. It consists of implementing a modified version of the Miller algorithm. Unfortunately, this method introduces an additional calculation which is expensive and which, moreover, does not really improve security for the Miller algorithm. 3. The use of the homogeneity property of projective coordinates (or Jacobian coordinates) of the point P or of the point Q, is a counter-measure against this attack on Miller's algorithm [START_REF] Page | Fault and Side-Channel Attacks on Pairing Based Cryptography[END_REF]. This counter-measure causes a modification of the result of the Miller algorithm: thus, we obtain a multiple of the final result of Miller algorithm. However, this is not a problem because this multiple being in the subfield of F p k , will be eliminated by the final exponentiation. The advantage of this counter-measure is that it is not expensive. However, unfortunately, El Mrabet et al. have shown in [START_REF] Lashermes | Practical validation of several fault attacks against the miller algorithm[END_REF] that this is not enough to protect Miller's algorithm. 4. Another efficient counter-measure is to choose two integers a et b such that a × b = 1 mod r [START_REF] Page | Fault and Side-Channel Attacks on Pairing Based Cryptography[END_REF]. Then, the idea is to compute the pairings between the points [a]P and [b]Q instead of computing the pairing between P and Q. This operation is possible thanks to the bilinearity of pairings, we have:

e([a]P, [b]Q) = e(P, Q) ab .
Also, thanks to the final exponentiation, the exponent, ab will be canceled. Finally, we obtain:

e(P, Q) ab×( p k -1 r ) = e(P, Q) ( p k -1 r ) .
5. Finally, we cite the most considered counter-measure: masking [START_REF] Page | Fault and Side-Channel Attacks on Pairing Based Cryptography[END_REF]. When the secret is the point P , this counter-measure consists of choosing a random point R of the elliptic curve E such that e(P, R) is defined, then computing the pairing between P and Q passing through the point R. Using the property of bilinearity, we obtain: e(P, Q) = e(P, Q + R) × e(P, -R).

Therefore, this counter-measure consists of computing two pairings instead of one and then performing their product.

The last two counter-measures are the most considered in practice. In a theoretical study of the various existing counter-measures, El Mrabet et al. suggested in [START_REF] Mrabet | A survey of fault attacks in pairing based cryptography[END_REF] to consider the counter-measure by masking to protect the Miller algorithm against SCA attacks. The counter-measure based on the bilinearity of pairings implies the computation of two scalar multiplications, one on E(F p ), and one on E(F p k ), both of them are scalar multiplication by an integer of order r. Make these two computations secure against side channel attack is a classical problem in elliptic curve cryptography. In [START_REF] Benoit Feix | Alexandre Side-channel analysis on blinded regular scalar multiplications[END_REF], Feix at al. demonstrate that even a secure scalar multiplication can be sensitive to side channel attack. The use of that counter-measure in pairing-based cryptography would imply to add two other sensitive computations with their counter-measure [START_REF] Meloni | New point addition formulae for ECC applications[END_REF][START_REF] Herbaut | Random euclidean addition chain generation and its application to point multiplication[END_REF][START_REF] Fournaris | Fault and Power Analysis Attack Protection Techniques for Standardized Public Key Cryptosystems[END_REF]. On the other side the masking counter-measure, imply one addition over E and in the worst case two pairing computations shorter than the scalar multiplication. As a consequence, we consider the masking counter-measure in the sequel of the article.

Resistance of our methods against non-invasive attacks

Since Miller-Fibonacci algorithm is based on the operation used in the Miller-Euclide so that proving that it is secure against a non invasive would show that the Miller-Euclide is also secured.

The Miller-Fibonacci algorithm is based on writing the pairing parameter s as a sequence of Fibonacci numbers. In our new algorithm, we have only additions steps. Therefore we only have to evaluate the line passing through the points T 1 and T 2 (or T 1 and Q) of respective coordinates (X 1 , Y 1 , Z) and (X 2 , Y 2 , Z) evaluated on the public point P . Therefore, the only equation involving known data that the non invasive attack can use is the following one:

l T1,T2 (P ) = y P Z 3 (X 1 -X 2 ) 3 -Y 1 (X 1 -X 2 ) 3 -(Y 1 -Y 2 ) x P Z 2 (X 1 -X 2 ) 2 -X 1 (X 1 -X 2 ) 2
The DPA attack model in pairing-based cryptography [START_REF] Page | Fault and Side-Channel Attacks on Pairing Based Cryptography[END_REF][START_REF] Mrabet | A practical differential power analysis attack against the miller algorithm[END_REF], is efficient when we have an operation that involved a known data and one of the temporary variables X 1 , X 2 , Y 1 , Y 2 or Z. If the secret is the point Q, the known values are x p and y p . With a DPA we can find Z 2 (X 1 -X 2 ) 2 using the multiplication x P Z 2 (X 1 -X 2 ) 2 and Z 3 (X 1 -X 2 ) 3 using y P Z 3 (X 1 -X 2 ) 3 . This gives us access to the value Z(X 1 -X 2 ). Recall that our algorithm is described using the co-Z arithmetic, as a consequence the value of Z is the same for the points at each execution step of the Miller-Fibonacci algorithm. We analyze now which information can be found given two hypotheses on Z, either Z Q = 1 or Z Q = 1. If Z Q = 1, i.e. Z = 1, then using a vertical DPA attack we can find the values (X 1 -X 2 ), (X 1 -X 2 ) 2 , and (X 1 -X 2 ) 3 . Then using another DPA we can find either Y 1 using Y 1 (X 1 -X 2 ) 3 or X 1 using X 1 (X 1 -X 2 ) 2 . With Z = 1 and X 1 or Y 1 and the elliptic curve equation, we can find Y 1 or X 1 . The attack would be successful. But if Z Q = 1, we cannot perform the attack. Indeed, we made the hypothesis that Q is secret, then Z Q is also secret and different from 1. We do not know either Z, neither X 1 and X 2 as they are derived from Q. We can make hypotheses on the value of Z and then on (X 1 -X 2 ), but we cannot try all the possible values of Z in F p . As a consequence, as long as Z Q = 1, our Miller-Fibonacci algorithm, based on presenting the parameter pairing s with the Zeckendorf representation, is natively protected against a DPA attack. Indeed, once we eliminate the leakage of information during the operations y P Z 3 (X 1 -X 2 ) 3 and x P Z 2 (X 1 -X 2 ) 2 , the other registers depend only on the secret point and consequently are protected against any vertical DPA analysis. In fact, our algorithm is also secure against horizontal attacks. Indeed, in horizontal attacks, the hypothesis is that a same entry is used more than once during one single execution of the algorithm. Using the co-Z arithmetic, we use the same Z coordinates during one step of the Miller algorithm. It is not sufficient to perform an horizontal SCA.

The same analysis can be done for Miller-Euclide which is also protected against SCA as long as Z = 1.

Comparison

In the Section 7, we have proved that our new algorithms resist against DPA attacks. In this section, we compare their respective complexities with a classical Miller algorithm protected again DPA attack. We choose to compare these two counter-measures in the context of computing Optimal Ate pairing on BN and BLS 12 curves for the reasons explained in the Introduction Remark 81 Our methods for computing the Miller loop using Euclidean addition chains or Fibonacci Sequences are available for computing all pairings and in any elliptic curve.

We have to evaluate f 6u+2,Q (P ) using Miller algorithm where 6u + 2 is an integer of 114 bits and Hamming weight 4. We present the following table the cost of our counter-measure by comparing it with the counter-measure by masking. 4. Comparison between the costs of the counter-measures Table 4 shows that our counter-measures are more efficient than the counter-measure proposed in the literature. The Miller-Euclide and the Miller-Fibonacci algorithm are, respectively, 16% and 17% faster than the Miller standard algorithm combined with the masking.

Algorithm

Conclusion

In this paper, we have proposed two new versions of Miller's algorithm: Miller-Euclide and Miller-Fibonacci. The complexity of these algorithms depends on the Miller loop parameter. Therefore, we showed how to choose a short representation for the parameter s which makes our algorithms as efficient in terms of complexity as possible. We proved also in this paper that these two algorithms resist against Differential Power Analysis attacks.

We compare the complexity of our new algorithm Miller-Fibonacci with the most efficient counter-measure of the classical Miller algorithm. We proved that our proposal Miller-Fibonacci algorithm is more efficient than the counter-measure by masking about 10.7% in the worst case and in the best case about 20%. // We only investigate the vertical power analysis, we leave as an open problem the analysis of our algorithms against horizontal power analysis. Indeed, this kind of analysis needs implementations and practical attacks.

Example 32

 32 We have presented 53 as an Euclidean addition chain. Let g = 34, gcd(34, 53) = 1. Let apply now the subtractive form of Euclid's algorithm:

Algorithm 7 :

 7 Miller-Fibonacci: The computation of f s,Q (P ) using Fibonacci sequences Require: P ∈ G1, Q ∈ G2, (dn, . . . , . . . d2) The Zeckendorf representation of s, Ensure: fu,Q(P ) 1: (T1, T2) ← (Q, Q) 2: (f1, f2) ← (1, 1) 3: for i = n -1 down to 2 do 4: if di = 1 then 5:

Table 1 .

 1 Cost of necessary operations for Miller loop

	Operation	Cost in Fp
	Multiplication in F p 12	54 M
	Sparse Multiplication in F p 12	39 M
	Squaring in F p 12	36 M
	Multiplication in F p 2	3 M
	Squaring in F p 2	2 M
	Addition step for Miller classic ( Projective coordinates)	80 M
	Doubling step for Miller classic ( Projective coordinates) 100 M
	Initialization step for Miller-Euclide/Fibonacci	27 M
	Addition step for Miller-Euclide	118 M
	Addition step for Miller-Fibonacci(for ui = 0)	121 M
	Addition step for Miller-Fibonacci (for ui = 1)	146 M

Table 2 .

 2 Cost of each Miller algorithm for BN curves

	Method	Complexity in
	Miller classic Algo. 1	12068 M
	Miller-Euclide Algo. 4	19991 M	+65, 5%
	Miller-Fibonacci Algo. 7	20683 M	+71, 3%

Table 3 .

 3 Cost of each Miller algorithm for BLS12

	Method	Complexity in Fp
	Miller classic Algo. 1	7438 M
	Miller-Euclide Algo. 4	13147 M	+76, 8%
	Miller-Fibonacci Algo. 7	13936 M	+87, 3%


Loubna Abou

El Qassim

Jaime Alonso

Ke Zhao

Sandrine Le Guillou

Jorge Diez

Fernando Vicente

Manuel Fernández-Sanjurjo

Eduardo Iglesias-Gutiérrez

Leluo Guan

Luis J Royo
email: ljroyo@serida.org

Loubna Abou El Qass

Jorge Di

Differences in the microRNAs Levels of Raw Milk from Dairy Cattle Raised under Extensive or Intensive Production Systems

Keywords: milk, miRNA, dairy production systems, biomarker

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Consumers' growing concern about food characteristics has contributed to the creation of a new concept of quality, particularly for animal products. This concept includes traditional attributes related to nutritional value, flavor, aroma, and color, together with new indicators related to ethical aspects, such as animal welfare and environmental impact of the production system [START_REF] Luykx | An Overview of Analytical Methods for Determining the Geographical Origin of Food Products[END_REF]. Consumers assume that products from cows raised under pasture and/or grazing are more natural and better for meeting animal welfare demands than those from the cows raised under the cereal-rich diets typical of intensive production systems. Intensive systems are considered less sustainable than production based on pastures, particularly because they have a larger ecological footprint and because they divert cereals from human consumption. These considerations highlight the need to compare whether agri-food products obtained from one production system or the other may have advantages, so that consumers can make fully informed purchasing decisions [START_REF] Valenti | Infrared Spectroscopic Methods for the Discrimination of Cows' Milk According to the Feeding System, Cow Breed and Altitude of the Dairy Farm[END_REF].

Milk is a complex secretory product and a source of nutrients for children and adults. Milk also contains numerous biomolecules, including microRNAs (miRNAs), a group of non-coding RNAs that bind to specific regions within messenger RNA to regulate gene expression post-transcriptionally [START_REF] He | MicroRNAs: Small RNAs with a Big Role in Gene Regulation[END_REF]. In vitro models [START_REF] Thorning | Milk and Dairy Products: Good or Bad for Human Health? An Assessment of the Totality of Scientific Evidence[END_REF] have suggested that exogenous miRNAs, such as bovine milk miRNAs, may influence the health of humans due to their resistance to stomach digestion, and that they may affect cellular function. Nevertheless, the potential bioactivity of dietary miRNAs is still under study, and several aspects remain poorly understood, including the dose needed to produce biological effects [START_REF] De Las Hazas | Dietary bovine milk miRNAs transported in extracellular vesicles are partially stable during GI digestion, are bioavailable and reach target tissues but need a minimum dose to impact on gene expression[END_REF]. The miRNA level of milk depends on genetic factors, such as the animal's breed and physiology [START_REF] Billa | Deep RNA-Seq Reveals MiRNome Differences in Mammary Tissue of Lactating Holstein and Montbéliarde Cows[END_REF][START_REF] Guillou | Characterization of Holstein and Normande Whole Milk MiRNomes Highlights Breed Specificities[END_REF], as well as environmental factors, such as the animal's lifestyle, pathological state, and diet [START_REF] Miretti | MicroRNAs as Biomarkers for Animal Health and Welfare in Livestock[END_REF][START_REF] Tzelos | MiRNAs in Milk Can Be Used towards Early Prediction of Mammary Gland Inflammation in Cattle[END_REF][START_REF] Muroya | Grazing Affects Exosomal Circulating MicroRNAs in Cattle[END_REF][START_REF] Muroya | Differences in Circulating MicroRNAs between Grazing and Grain-Fed Wagyu Cattle Are Associated with Altered Expression of Intramuscular MicroRNA, the Potential Target PTEN, and Lipogenic Genes[END_REF]. Furthermore, the miRNA level depends on the milk fractions, milk fat, whey, and cells [START_REF] Li | Comparative Analysis of the MiRNome of Bovine Milk Fat, Whey and Cells[END_REF]. The miRNAs found in milk are highly resistant to acidic pH and to enzymes such as RNases [START_REF] Zhang | Exogenous Plant MIR168a Specifically Targets Mammalian LDLRAP1: Evidence of Cross-Kingdom Regulation by MicroRNA[END_REF], which implies that they are very stable and resistant to industrial treatments [START_REF] Zhang | Exogenous Plant MIR168a Specifically Targets Mammalian LDLRAP1: Evidence of Cross-Kingdom Regulation by MicroRNA[END_REF][START_REF] Chen | Identification and Characterization of MicroRNAs in Raw Milk during Different Periods of Lactation, Commercial Fluid, and Powdered Milk Products[END_REF][START_REF] Izumi | Bovine Milk Contains MicroRNA and Messenger RNA That Are Stable under Degradative Conditions[END_REF]. For these reasons, miRNAs have been proposed as biomarkers of quality control in dairy products, and they have been used to control for fraud in the labeling of milk powder [START_REF] Chen | Identification and Characterization of MicroRNAs in Raw Milk during Different Periods of Lactation, Commercial Fluid, and Powdered Milk Products[END_REF].

Changes in diet have been shown to affect the expression of genes in the mammary gland [START_REF] Bauman | Nutrigenomics, Rumen-Derived Bioactive Fatty Acids, and the Regulation of Milk Fat Synthesis[END_REF][START_REF] Ollier | Whole Intact Rapeseeds or Sunflower Oil in High-Forage or High-Concentrate Diets Affects Milk Yield, Milk Composition, and Mammary Gene Expression Profile in Goats[END_REF]. Considering that miRNAs are essential regulators of gene expression, miRNA profiles in milk may also change with diet. However, few studies have investigated the effect of diet on miRNA profiles in bovine milk [START_REF] Quan | Replacement of Forage Fiber with Non-Forage Fiber Sources in Dairy Cow Diets Changes Milk Extracellular Vesicle-MiRNA Expression[END_REF], such as through comparisons of the quite different diets in intensive or extensive production systems. Therefore, our objective was to evaluate differences in miRNA levels of bovine milk produced under extensive or intensive systems. Characterizing production systems as intensive or extensive is not always straightforward, given the lack of regulatory definitions and the complexity of the factors involved [START_REF] Marin | Sistema de Explotación Ganadera: Notas En Torno a Su Concepto[END_REF]. In the present study, we defined intensive systems as those without grazing and with high amounts of maize silage and concentrate in the diet. We defined extensive systems as those with grazing and a diet with a low amount of concentrate and no maize silage. Our work was carried out in Asturias, north of Spain, where the edapho-climatic conditions permit the coexistence of different milk production systems, varying from intensive systems, with animals fed indoors, to extensive pastured-based systems. This results in a fairly wide range of production systems covering all management possibilities [START_REF] Ruiz | Definición y Caracterización de la Extensividad en las Explotaciones Ganaderas en España. Situación de la Ganadería Extensiva en España (I)[END_REF]. Furthermore, we also explored the use of miRNA markers for the characterization and traceability of agri-food products. We used bioinformatics to predict what target genes may be regulated by miRNAs whose levels differ between the two production systems, in an effort to examine how certain miRNAs may reveal the functional properties of agri-food products.

Materials and Methods

Sample Collection and Preparation

In order to maximize potential differences in milk composition between intensive and extensive production systems, for sequencing, we sampled six farms at the end of one spring season (May 2016). On three farms (extensive production), animals grazed for at least for 12 h each day and ate a diet based on fresh grass with a small amount of concentrates. On three farms (intensive production), animals did not graze, instead eating a diet based on conserved feed and concentrates. Farm characteristics and details of the diets are described in Table 1. All animals on all farms were Holstein cows. For validation by RT-qPCR, twenty farms were sampled during autumn 2017 and spring 2018. Ten dairy farms applied an extensive system in which animals grazed and consumed a diet based on fresh grass, without maize silage and with a low amount of concentrates (4-8 kg/animal/day), while the other ten applied an intensive system in which animals were housed and ate a diet consisting of a high dry matter intake of maize silage (16-30 kg/animal/day) and a high amount of concentrate (>10 kg/animal/day) (Table 2). From each type of farm, bulk tank milk was sampled after an even number of milkings (to avoid differences due to afternoon and morning milk composition). Samples were maintained at 4 • C and immediately transported to the laboratory for processing. Milk was well mixed, and 50 mL of each sample was centrifuged at 1900× g for 20 min. The fat in the upper phase was transferred to a new 50-mL RNase-free tube. Then, 7.5 mL of QIAzol lysis reagent (Qiagen, Crawley, UK) was added, and the emulsion was vigorously mixed until the fat was well dispersed. The pellet (cellular fraction) from the initial centrifugation was washed twice with phosphate-buffered saline, then homogenized with 1 mL of QIAzol lysis reagent. All samples were stored at -80 • C until RNA extraction.

RNA Isolation

For each sample, total RNA was extracted from 2 mL of milk fat in QIAzol lysis reagent, and from 1 mL of cell pellet resuspended in QIAzol lysis reagent. RNA was extracted using the miRVana miRNA isolation kit (Applied Biosystems, Foster City, CA, USA) following the manufacturer's instructions, then stored at -80 • C. The concentration and integrity of RNA (RIN: RNA Integrity Number) for sequencing was further determined on an Agilent 2100 Bioanalyzer using an RNA 6000 Pico kit (both from Agilent Technologies, Santa Clara, CA, USA).

Search for miRNAs Candidates from Sequencing 2.3.1. RNA Sequencing

In order to identify miRNAs differing between the two production systems, 12 libraries corresponding to the fat or cellular fractions of milk from six samples of bulk tank milk (three per production system) were prepared and sequenced using the Illumina platform (Illumina, San Diego, CA, USA) as 50 bp single reads. The raw sequence data were processed for quality control, and low quality reads were removed from raw data using CASAVA 1.8 based on chastity. Then, adaptors were trimmed, and sequences with read lengths between 15 to 40 nt were mapped to the bovine genome (bostau 7) and the miRNA database (miRBase, release version 21) in order to identify the known miRNAs. Expression levels of each miRNA were estimated based on the frequency of reads, and results were normalized to the number of reads per million (RPM) using the following formula: RPM = (specific miRNA reads number/total mapped miRNA reads per library) × 10 6 .

Identification of Reference miRNAs for qPCR Normalization

In order to select the miRNAs to be used as candidates for normalizers in RT-qPCR, we chose those miRNAs with more stable expression among samples for each milk fraction, that is, miRNAs with the smallest coefficient of variation (CV = standard deviation/mean).

Identification of miRNAs Whose Levels Differed between Production Systems

In order to identify those miRNAs whose levels differ between production systems, the results from miRNA sequencing were analyzed using three statistical tests. One test was the ratio of the difference between the means to the sum of the standard deviations between the two production systems: value = |(mean1-mean2)|/(standard deviation1 + standard deviation2) (strict Cohen's d). This ratio indicates how many sums of deviations fit between the means; the higher the ratio is, the greater the difference between the means [START_REF] Cohen | Statistical Power Analysis for the Behavioral Sciences[END_REF]. Another test was Student's t test, for which lower t values indicated greater differences between the means for the two production systems [START_REF] Student | The Probable Error of a Mean[END_REF]. The third test was the absolute value of the correlation coefficient; the higher this value was, the greater the difference between the two systems [START_REF] Pearson | Note on Regression and Inheritance in the Case of Two Parents[END_REF]. Afterward, to choose the miRNAs whose levels differed between production systems, they were ranked first according to each statistical test, and then according to the average of the classification of the three tests.

Validation of Candidate miRNAs Using RT-qPCR 2.4.1. RT-qPCR Analysis

A subset of the miRNAs which were determined to differ significantly between the two production systems identified in Section 3.2 were validated using quantitative RT-PCR and milk samples from twenty dairy farms featuring extensive or intensive production systems (Table 2).

Total RNA was used for cDNA synthesis using the TaqMan Advanced miRNA cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA), and the resulting cDNA was stored at -20 °C until use. Levels of miRNAs were determined using quantitative RT-PCR (TaqMan Advanced miRNA Assays; ThermoFisher Scientific, Waltham, MA, USA) in a StepOne thermocycler (Applied Biosystems, Foster City, CA, USA). The final reaction solution contained 10 µL of 2 × TaqMan Fast Advanced Master mix (ThermoFisher Scientific, Waltham, MA, USA), 1 µL of 20× TaqMan Advanced miRNA Assay (ThermoFisher Scientific, Waltham, MA, USA), 4 µL of RNase free water, and 5 µL of cDNA (diluted 1:10). The thermocycling program was set at 95 • C for 20 s, followed by 40 cycles at 95 • C for 1 s and 60 • C for 20 s. All PCR reactions were performed in duplicate, and a maximum of 0.5 threshold cycles were permitted between duplicates.

Selection of Stable Reference miRNAs. GeNorm Analysis

Normalization is an essential component of a reliable qPCR assay. geNorm [START_REF] Vandesompele | Accurate Normalization of Real-Time Quantitative RT-PCR Data by Geometric Averaging of Multiple Internal Control Genes[END_REF] is one of the most popular algorithms to find stable reference genes from a set of tested candidate reference genes in a given experimental condition. We used geNorm to find the optimal number and choice of reference genes for normalization, using miRNAs identified in Section 3.2 in 22 tank milk samples, representing the experimental variation of the dairy production systems existing in the area of study [START_REF] Abou El Qassim | Variaciones de Los Perfiles de MicroARN en la Leche Cruda de Vaca Según el Sistema de Alimentación[END_REF].

miRNAs Levels Normalization and Estimation

Levels of miRNA were normalized based on the geometric mean of the selected reference miRNAs selected by geNorm, estimated using QBase+ 3.1 software, Biogazelle (Gent, Belgium) [START_REF] Hellemans | QBase Relative Quantification Framework and Software for Management and Automated Analysis of Real-Time Quantitative PCR Data[END_REF] and expressed in base log10. Unless otherwise noted, results were reported as mean ± standard deviation. Mean miRNAs levels between extensive and intensive production systems were compared using Student's t test in R-Commander 2.7-1. Significance was defined as p < 0.05.

Prediction and Functional Analysis of Genes Targeted by miRNAs

The Target Scan 7.2 bioinformatics tool [START_REF] Agarwal | Predicting Effective MicroRNA Target Sites in Mammalian MRNAs[END_REF] was used to predict bovine mRNA targets of candidate miRNAs. A pathway analysis of targeted genes was performed using Panther bioinformatics tool version 16.0 (http://www.pantherdb.org/ accessed on 8 June 2021) based on Gene Ontology classification [START_REF] Mi | PANTHER Version 16: A Revised Family Classification, Tree-Based Classification Tool, Enhancer Regions and Extensive API[END_REF].

Results

miRNAs Levels in Fat and Cellular Fractions of Milk

The total mean RNA concentration was 192.4 ± 37.4 ng/µL in milk fat and 30.6 ± 15.8 ng/µL (mean ± SD) in milk cells. The RIN value was 2.6 ± 0.15 for RNA from milk fat, and 6.5 ± 0.3 for RNA from milk cells.

The six libraries from the cellular fraction of milk yielded a mean of 24,017,290.17 reads, significantly more than the 6,964,122.33 reads from the six libraries from the fat fraction (p = 0.004, Table 3). Almost half of the reads came from small RNAs, of which the most abundant were transfer RNAs (tRNAs) in the cellular fraction and non-coding RNAs (ncRNAs) in the fat fraction. Significant differences were found between the fat and cellular fractions for a percentage of all small RNAs, except for small nuclear RNAs (snRNAs) and miRNAs. Ribosomal RNAs (rRNAs), small nucleolar RNAs (snoRNAs), and ncRNAs were more abundant in the fat fraction than in the cellular fraction, while the converse was true for tRNAs. The levels of miRNAs in the fat fraction differed significantly between extensive and intensive systems (p = 0.040, Table 4). Intensive production was associated with higher miRNAs levels. We identified 518 known miRNAs in the cellular fraction of milk and 477 in the fat fraction. Most of these miRNAs (454) were present in both fractions (Table S1).

Validation of miRNAs Whose Levels Differed between Intensive and Extensive Production

The first 10 miRNAs in the fat fraction whose levels differed between the two production systems are ranked in Table 5, and those in the cellular fraction are ranked in Table 6. Levels of the first five miRNAs from each fraction were validated using quantitative RT-PCR and milk samples from another 20 farms. The following miRNAs in the fat fraction were subjected to validation: bta-miR-215, bta-miR-369-3p, bta-miR-6520, bta-miR-7863, and bta-miR-133a. The following miRNAs in the cellular fraction were subjected to validation: bta-miR-574, bta-miR-3432a, bta-miR-2285e, bta-miR-197, and bta-miR-2284y. Six miRNAs in each milk fraction were chosen as candidates for normalization (Table 7). Following GeNorm, normalization of levels in the fat fraction was optimal using the geometric mean of bta-miR-151-3p and bta-miR-30a-5p. In the case of cellular fraction, GeNorm recommended the use of the geometric mean of the three most stable miRNAs (bta-miR-103, bta-miR-107, bta-miR-28).

Based on analysis of the normalized miRNAs levels (Figures 1 and2), the only miRNA in the fat or cellular fraction that differed significantly between intensive and extensive production systems was bta-miR-215, in the fat fraction. This miRNA was significantly more abundant in the fat fraction of milk from intensive production (p = 0.030). The miRNAs bta-miR-2284y and bta-miR-2285e in the cellular fraction showed a trend towards lower levels in milk from extensive production, but the differences did not reach statistical significance.

Putative target gene and pathway analyses.

Using Target Scan, 143 potential target genes were identified for bta-miR-215, which was abundant in intensive compared to extensive dairy production (Table S2). Among these targets, one gene was particularly involved in lipid metabolism and energy metabolism: the fatty acid binding protein 3 (Fabp3).

The pathway analysis of bta-miR-215 target genes allowed for the identification of 41 associated biological pathways. The gonadotropin-releasing hormone receptor pathway and the Transforming Growth Factor β (TGF-β) signaling pathway (Figure 3) were highlighted, as they include the most bta-miR-215 target genes, i.e., Activin Receptor Type-2B (Acvr2b), Activin A Receptor, Type 2A (Acvr2a), Mitogen-Activated Protein Kinase 1 (Mapk1), and bone morphogenetic protein receptor type II (Bmpr2).

Putative Target Gene and Pathway Analyses

Using Target Scan, 143 potential target genes were identified for bta-miR-215, which was abundant in intensive compared to extensive dairy production (Table S2). Among these targets, one gene was particularly involved in lipid metabolism and energy metabolism: the fatty acid binding protein 3 (Fabp3).

The pathway analysis of bta-miR-215 target genes allowed for the identification of 41 associated biological pathways. The gonadotropin-releasing hormone receptor pathway and the Transforming Growth Factor β (TGF-β) signaling pathway (Figure 3) were highlighted, as they include the most bta-miR-215 target genes, i.e., Activin Receptor Type-2B (Acvr2b), Activin A Receptor, Type 2A (Acvr2a), Mitogen-Activated Protein Kinase 1 (Mapk1), and bone morphogenetic protein receptor type II (Bmpr2).

Discussion

The objective of this study was to identify a set of miRNAs in milk whose levels differed with the production system, in order to (1) determine whether miRNAs profiling can be used to authenticate milk from a given production system, and (2) investigate whether the production system can influence the functional properties of agri-food products such as bovine milk. To address these questions, we sampled milk from dairy farms in Asturias in northern Spain that applied intensive or extensive production practices, which differ significantly in numerous characteristics that can alter milk quality, including feeding management, animal density, and access to exercise through grazing [START_REF] Vicente | Capacity of Milk Composition to Identify the Feeding System Used to Feed Dairy Cows[END_REF].

The origin of miRNAs in milk is controversial, and most miRNAs in milk are not found in blood [START_REF] Alsaweed | Human Milk MiRNAs Primarily Originate from the Mammary Gland Resulting in Unique MiRNA Profiles of Fractionated Milk[END_REF]. In addition, the miRNA profile in milk differs across the fractions of fat, whey, and cells [START_REF] Li | Comparative Analysis of the MiRNome of Bovine Milk Fat, Whey and Cells[END_REF]. In the present work, we decided not to study the whey fraction because its miRNA content is lower, its miRNA profile is highly similar to that of milk fat [START_REF] Li | Comparative Analysis of the MiRNome of Bovine Milk Fat, Whey and Cells[END_REF], and it is not exploited in certain dairy industry practices, such as cheese production.

We successfully isolated total RNA from fat and cellular fractions of bovine milk, especially considering that the samples came from milk tanks on commercial dairy farms. The RIN was low for RNA from milk fat, which likely reflects its abundant content of low-molecular-weight RNA, which is different from ribosomal RNA [START_REF] Li | Comparative Analysis of the MiRNome of Bovine Milk Fat, Whey and Cells[END_REF].

RNA sequencing analysis confirmed different RNA profiles in the fat and cellular fractions of bovine milk (Table 3). Interestingly, we did not identify significant differences in miRNAs levels between the two fractions, although we did obtain what appears to be the first evidence that miRNAs in the fat fraction differ between milk from intensive or extensive production systems. Given that milk has been proposed as a major epigenetic modulator of the gene expression of the milk recipient [START_REF] Melnik | Milk Exosomes and MicroRNAs: Potential Epigenetic Regulators[END_REF], and its modulation can be dependent on the amount of miRNAs [START_REF] Pozo-Acebo | Bovine Milk-Derived Exosomes as a Drug Delivery Vehicle for Mirna-Based Therapy[END_REF], our results imply that the animal production system can influence the functional properties of agri-food products of animal origin.

Sequencing results did not allow for the identification of miRNAs that were specific to a given production system in either the fat or cellular fraction. Therefore, we focused on miRNAs whose levels differed between the two systems, and we validated a subset of the promising miRNAs using quantitative RT-PCR. The only miRNA that we validated to differ significantly between the two production systems was bta-miR-215 in the fat fraction (Table 5). This miRNA was upregulated in milk from intensive production. At the moment, we can only speculate as to why intensive production might upregulate this miRNA. One possibility is that its upregulation somehow compensates for poor efficiency of feed conversion into milk: dairy cows with medium potential show lower conversion efficiency during indoor feeding than during grazing in pastures [START_REF] Knaus | Perspectives on Pasture versus Indoor Feeding of Dairy Cows[END_REF], and Angus cows, less efficient at feed conversion, show upregulation of bta-miR-215 [START_REF] Al-Husseini | Characterization and Profiling of Liver MicroRNAs by RNA-Sequencing in Cattle Divergently Selected for Residual Feed Intake[END_REF]. Consistent with this possibility is that heat and oxidative stress, which can easily occur in intensive production systems, upregulate bta-miR-215 in the serum of Holstein cows [START_REF] Zheng | Identification and Bioinformatics Analysis of MicroRNAs Associated with Stress and Immune Response in Serum of Heat-Stressed and Normal Holstein Cows[END_REF].

Among the 143 genes identified by bioinformatic tools, Fabp3 stands out for its involvement in fatty acid transport and activation in the bovine mammary gland [START_REF] Zhang | Abundance of Solute Carrier Family 27 Member 6 (: SLC27A6) in the Bovine Mammary Gland Alters Fatty Acid Metabolism[END_REF]. Vargas-Bello-Perez et al. (2020) [START_REF] Vargas-Bello-Pérez | Effects of Dietary Vegetable Oils on Mammary Lipid-Related Genes in Holstein Dairy Cows[END_REF] demonstrated that when the diet of Holstein dairy cows is supplemented with hydrogenated vegetable oil, Fabp3 is downregulated in milk somatic cells. Likewise, when fed a high-concentrate diet, the fatty acid transporter Fabp3 is inhibited in mammary glands [START_REF] Ma | A High-Concentrate Diet Induces an Inflammatory Response and Oxidative Stress and Depresses Milk Fat Synthesis in the Mammary Gland of Dairy Cows[END_REF]. In our study, bta-miR-215 increased in milk in response to intensive farming conditions, which are generally associated with high-energy diets, which is consistent with its likely involvement in the downregulation of Fabp3. A direct relationship between bta-miR-215 and Fabp3 has, indeed, been validated in bone marrow mesenchymal stem cells [START_REF] Zhou | Exosomes Isolated from the MiR-215-Modified Bone Marrow Mesenchymal Stem Cells Protect H2O2-Induced Rat Myoblasts via the MiR-215/FABP3 Pathway[END_REF].

Apart from Fabp3, another two genes are of particular interest: Mapk1 and Bmpr2. Mapk1 is involved in the regulation of milk protein synthesis [START_REF] Lu | Proteomic and Functional Analyses Reveal MAPK1 Regulates Milk Protein Synthesis[END_REF], and Bmpr2 is associated with the glucose metabolism and insulin response. In fact, when Bmpr2 is altered, it likely blunts glucose response and lipid uptake [START_REF] Hemnes | BMPR2 Dysfunction Impairs Insulin Signaling and Glucose Homeostasis in Cardiomyocytes[END_REF]. Altogether, they likely have a negative impact on protein and lactose synthesis in animals bred under intensive dairy systems.

We also identified the miRNAs bta-mir-2284y and bta-mir-2285e in the cellular fraction of milk whose levels tended to differ between intensive and extensive production systems. Further analysis showed that these differences became significant when we compared animals on diets containing more or fewer than 10 kg of concentrates per day (data not shown). These results suggest that our failure to detect significant differences in levels of bta-mir-2284y and bta-mi-2285e in the overall analysis may reflect milk sampling from commercial farms, as well as the difficulty of defining intensive and extensive production. Further studies should examine these two miRNAs as additional potential biomarkers for authentication and functional analysis of agri-foods.

Finally, our study presents several limitations. The study was carried out in commercial farms, where managing conditions and diets are quite different even within the same group. Validation in controlled conditions at experimental farms, which reduce internal variance of the experimental groups, should help to identify miRNAs as putative biomarkers for dairy production systems. In vitro experiments will help to elucidate the functional features of milk produced under different dairy systems.

Conclusions

We investigated differences in miRNA profiles of raw cow tank milk from commercial farms which applied either intensive or extensive production systems. We identified bta-miRNA-215 in the fat fraction of milk as a possible biomarker of milk from intensive production systems. Our results imply that the type of production system can influence miRNA levels, and, therefore, functional properties of bovine milk, as well as, potentially, other agri-food products of animal origin.

Supplementary Materials:

The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/vetsci9120661/s1, Table S1: Sequencing RPM; Table S2: Selection of Predicted bta-miR-215 Targets from Target Scan database.

Figure 1 .

 1 Figure 1. Average levels of miRNAs estimated by qRT-PCR in fat fractions of milk from extensive and intensive farms. On the X axis, production systems are represented; on the Y axis, average levels of each miRNA are represented. Levels of bta-miR-369-3p showed a skewed distribution, so they are shown using a box-and-whisker plot. The errors bars in the other plots indicate standard deviation.

Figure 1 .

 1 Figure 1. Average levels of miRNAs estimated by qRT-PCR in fat fractions of milk from extensive and intensive farms. On the X axis, production systems are represented; on the Y axis, average levels of each miRNA are represented. Levels of bta-miR-369-3p showed a skewed distribution, so they are shown using a box-and-whisker plot. The errors bars in the other plots indicate standard deviation.

Figure 2 .

 2 Figure 2. Average levels of miRNAs estimated by qRT-PCR in cellular fractions of milk from extensive and intensive farms. On the X axis, production systems are represented; on the Y axis, average levels of each miRNA are represented. The errors bars indicate standard deviation.

Figure 2 .

 2 Figure 2. Average levels of miRNAs estimated by qRT-PCR in cellular fractions of milk from extensive and intensive farms. On the X axis, production systems are represented; on the Y axis, average levels of each miRNA are represented. The errors bars indicate standard deviation.

Figure 3 .

 3 Figure 3. Functional Gene Ontology classification of bta-miR-215 pathway analysis of targeted genes.

Figure 3 .

 3 Figure 3. Functional Gene Ontology classification of bta-miR-215 pathway analysis of targeted genes.

Table 1 .

 1 Characteristics and diet of cows on the dairy farms where milk was sampled for RNA sequencing.

	Production System	Number of Cows	Milk Production (L/Day/Cow)	Grazing (h/Day)	Grass Silage (kg F */Day/Cow)	Maize Silage (kg F */Day/cow)	Hay (kg F */Day/Cow)	Concentrate (kg F */Day/Cow)
			30	0	17	10.0	2	10
	Intensive (3 farms)		28	0	10	15.0	10	11
			29	0	16	20.0	2	12
			21	20	10	0.0	2	7
	Extensive (3 farms)		31	>12	14	0.0	6	6
			29	18	15	0.0	3	6
			* kg F, kg of fresh matter.

Table 2 .

 2 Characteristics and diet from twenty dairy farms featuring extensive or intensive production systems, where milk was sampled for quantitative real-time PCR validation.

Production System Number of Cows Milk Production (L/Day Cow) Grazing (h/Day) Grass Silage (kg F */Day/Cow) Maize Silage (kg F */Day/Cow) Hay (kg F */Day/Cow) Concentrate (kg F */Day/Cow)

		124	37.4	0	8.0	30.0	0.8	11.5
		116	37.0	0	10.0	30.0	0.0	10.5
			29.0	0	16.0	20.0	3.0	12.0
		240	36.0	0	10.0	16.0	2.5	12.0
	Intensive	250	38.0	0	12.0	30.0	0.9	12.3
	(10 farms)		27.0	0	14.0	28.0	0.0	10.5
		110	30.0	0	16.0	20.0	0.0	11.0
			28.0	0	15.0	20.0	2.5	12.0
		118	36.0	0	16.0	22.0	0.0	10.5
		124	37.0	0	11.0	20.0	4.5	12.0
			21.0	6	10.0	0.0	6.8	6.5
			21.0	20	12.0	0.0	0.0	7.0
			26.2	22	0.0	0.0	6.0	4.0
		8	18.8	21	0.0	0.0	4.6	4.1
	Extensive		19.5	20	14.0	0.0	0.0	4.7
	(10 farms)		29.0	18	33.0	0.0	0.7	6.2
			27.0	18	0.0	0.0	4.0	8.0
		7	20.0	22	6.0	0.0	0.0	5.0
			23.0	21	0.0	0.0	0.0	5.0
			24.0	20	0.0	0.0	0.0	6.0

* kg F, kg of fresh matter.

Table 3 .

 3 Read mapping statistics.

			Total	Small RNA				Small RNAs (%)
			Reads	Reads	Total	rRNA	snRNA	snoRNA	tRNA	miRNA	ncRNA
			32,560,190	16,693,358	0.5	12.7	3.3	5.7	45.0	2.2	31.1
	Cellular fraction	Intensive Farms Extensive Farms	26,678,803 23,616,528 19,944,125 21,166,161	13,626,331 11,871,084 9,753,050 9,650,109	0.5 0.5 0.5 0.5	8.3 16.5 13.6 14.9	2.0 3.7 3.1 3.9	4.7 5.5 5.6 6.6	64.3 37.6 51.5 45.4	0.9 2.0 1.4 1.4	19.8 34.8 24.7 27.8
			20,137,934	10,053,989	0.5	15.9	3.9	6.1	36.8	2.3	34.9
			7,620,977	3,867,636	0.5	20.6	3.4	6.1	29.8	1.7	38.4
	Fat fraction	Intensive Farms Extensive Farms	7,209,138 6,180,232 7,056,743 6,682,264	3,139,312 3,115,110 2,599,438 3,279,814	0.4 0.5 0.4 0.5	21.9 26.0 19.4 25.0	3.7 3.8 3.7 4.1	6.6 6.0 8.4 7.1	25.3 15.8 24.8 16.9	1.9 1.6 2.2 2.4	40.6 46.9 41.5 44.5
			7,035,380	3,854,280	0.6	25.9	4.3	7.0	19.2	2.3	41.4

rRNA: ribosomal ribonucleic acid; snRNA: small nuclear RNA; snoRNA: small nucleolar RNA; tRNA: transfer RNA; miRNA: microRNA; ncRNA: non-coding RNA.

Table 4 .

 4 Differences in abundance of small RNA classes depending on milk fraction and production system.

		p-Value Based on Student's t Test
	Small RNA Class	Cellular vs. Fat Fraction *	Cellular Fraction: Extensive vs. Intensive	Fat Fraction: Extensive vs. Intensive
	rRNA	0.004	0.513	0.827
	snRNA	0.196	0.268	0.184
	snoRNA	0.024	0.127	0.050
	tRNA	0.004	0.827	0.513
	miRNA	0.260	0.825	0.040
	Non-coding RNA	0.004	0.827	0.513

* Both intensive and extensive production systems together. rRNA: ribosomal ribonucleic acid; snRNA: small nuclear RNA; snoRNA: small nucleolar RNA; tRNA: transfer RNA; miRNA: microRNA; ncRNA: non-coding RNA. In bold are significant differences.

Table 5 .

 5 Ranking average of the first ten differentially expressed miRNAs in the fat fraction, the three applied tests' values, and their rankings.

	miRNA		Result for		Ranking According to	Average Ranking
		Test 1 a Test 2 b Test 3 c	Test 1	Test 2	Test 3
	bta-miR-215 *	3.180	0.010	0.960	1	1	1	1.0
	bta-miR-369-3p *	1.760	0.020	0.900	3	2	2	2.3
	bta-miR-6520 *	1.360	0.030	0.850	4	3	4	3.7
	bta-miR-7863 *	1.970	0.080	0.860	2	7	3	4.0
	bta-miR-133a *	1.300	0.040	0.840	5	5	5	5.0
	bta-miR-532	1.260	0.040	0.840	6	4	6	5.3
	bta-miR-148a	1.210	0.120	0.780	7	13	7	9.0
	bta-miR-138	1.000	0.070	0.770	22	6	8	12.0
	bta-miR-450a	1.190	0.140	0.760	8	18	10	12.0
	bta-miR-6527	1.010	0.090	0.770	21	8	9	12.7

* These miRNAs were validated using quantitative RT-PCR. a Strict Cohen's d test. b Student's t test. c Correlation coefficient test.

Table 6 .

 6 Ranking average of the first ten differentially expressed miRNAs in the cellular fraction, the three applied tests' values, and their rankings.

	miRNA		Result for		Ranking According to	Average Ranking
		Test 1 a Test 2 b Test 3 c	Test 1	Test 2	Test 3
	bta-miR-574 *	5.770	0.000	0.990	1	1	1	1.0
	bta-miR-3432a *	5.520	0.010	0.980	2	3	2	2.3
	bta-miR-2285e *	2.540	0.010	0.950	5	2	3	3.3
	bta-miR-197 *	1.970	0.010	0.920	6	4	5	5.0
	bta-miR-2284y *	2.750	0.020	0.940	3	8	4	5.0
	bta-miR-219	1.740	0.010	0.910	9	5	7	7.0
	bta-miR-2397-3p	1.770	0.020	0.900	8	7	8	7.7
	bta-miR-2308	2.560	0.050	0.910	4	14	6	8.0
	bta-miR-2419-5p	1.620	0.020	0.890	11	6	9	8.7
	bta-miR-2409	1.790	0.040	0.890	7	12	10	9.7

* These miRNAs were validated using quantitative RT-PCR. a Strict Cohen's d test. b Student's t test. c Correlation coefficient test.

Table 7 .

 7 Ranking of the most stable miRNAs according to the coefficient of variation.

	Milk Fraction	miRNA	Coefficient of Variation
		bta-miR-532	0.060
		bta-miR-151-3p	0.070
		bta-miR-27b	0.090
	Fat	bta-miR-103	0.090
		bta-miR-30a-5p	0.090
		bta-miR-99a-3p	0.090
		bta-miR-103	0.080
		bta-miR-107	0.090
		bta-miR-181a	0.090
	Cellular	bta-miR-28	0.100
		bta-miR-345-3p	0.100
		bta-miR-28342	0.100

Vet. Sci. 2022, 9, x FOR PEER REVIEW 9 of 15

Funding: This research was funded by Principado de Asturias Regional Government (project IDI/2021/000102), co-financed by the European Union through the ERDF (European Regional Development Fund) and Ministry of Science and Innovation (AEI) project PID2021-126010OR-I00. Loubna Abou el qassim was founded by FICYT Severo Ochoa Grant (BP17-49).

Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.

Data Availability Statement: Not applicable.

Acknowledgments:

We want to acknowledge all dairy farmers that contribute to this study.

Author Contributions: Conceptualization, L.J.R., L.A.e.q., and L.G.; methodology, J.A., J.D., L.J.R., L.A.e.q., S.L.G., and K.Z.; formal analysis, L.A.e.q., S.L.G., M.F.-S., E.I.-G., and L.J.R.; investigation, L.A.e.q., L.J.R.; resources, L.A.e.q. and F.V.; data curation, L.J.R., K.Z., and L.A.e.q.; writing-original draft preparation, L.A.e.q. and L.J.R.; writing-review and editing, L.G., K.Z., S.L.G., M.F.-S., E.I.-G., and F.V.; visualization, L.A.e.q. and L.J.R.; supervision, L.J.R. and L.G.; project administration, L.J.R.; funding acquisition, L.J.R., L.G., and F.V. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:

The authors declare no conflict of interest.