F Z Boudara

H Bouzekri

Y Benhammadi

P.-H Cocquet

M Rivaletto

L Pécastaing

A Silvestre De Ferron

Y Le Guer

Accurate numerical modeling of convective heat transfer coefficient for a high power PLA-core toroidal electromagnetic coupler subject to cooling

Keywords:

The toroidal electromagnetic coupler has many advantages over the classical transformer. It is widely used in low-power applications. However, there are some limits that hamper its use in modern and high-power applications, such as the cost, the difficulty of winding and, especially, the limited number of studies carried out on the thermal aspect. In this article, we look at the problem of cooling a high-power and very compact toroidal electromagnetic coupler used in smart grid systems. The cooling is modeled thanks to heat transfer coefficients considering the windings as a surface heat source. A mathematical model which is quick and useful for a first approach is developed, and the results obtained show the temperature distribution for convective boundary conditions inside and outside the hollow torus. In addition, a parametric study is carried out. It aims to determine the optimal internal and external convective heat transfer coefficients for which the temperature of the coupler remains below the material's destruction threshold. Additionally, a 3D numerical model taking into account the wires of the toroidal coupler is studied and compared with the mathematical model. The good consistency between the two methods ensures that this simplified thermal modeling can be trusted for the study of thermal management during the design of future electromagnetic couplers.

Introduction

The increase in greenhouse gases observed in recent decades, mainly caused by non-renewable energy sources, has led to a rise in renewable ones such as solar, wind, or fuel cell energies. For their implementation, these sources require mastery of the power they deliver through adequate and intelligent flow management. In particular, the fast charging of autonomous systems will lead to the transfer of electrical power which is far greater than the capacities of the products currently on the market [START_REF] Nicholas | Lessons learned on early electric vehicle fast-charging deployments[END_REF].

With this in mind, compact toroidal electromagnetic couplers with a circular cross-section were developed [START_REF] Boudara | Design and development of high efficiency 150 kw very compact PLA core electromagnetic coupler for highly resonant power transfer technology[END_REF]. An example is presented in Fig. 1. The primary and secondary windings wound over the entire surface of the torus confine the magnetic field, offering the toroidal coupler very reduced electromagnetic radiation. It also has advantages related to size and volume. Thanks to these qualities, it is widely used in low-power applications such as public lighting, or even in the production of current or voltage sources with a very low internal resistance. However, the toroidal electromagnetic coupler is still not Fig. 1. The toroidal electromagnetic coupler. The primary winding and its electrical connection for the four parallel wires are represented by the color red. The secondary winding and its electrical connection are represented by the color blue. The gray part is the PLA cable holder. Right: The 3D printed PLA cable holder.

Table 1

Thermophysical properties of PLA [START_REF] Urbanlabs | Web site URBANLABS, technical Data sheet of the PLA[END_REF].

Property

Value Unit

Solid density 𝜌 1300 kg m -3 Heat capacity 𝐶𝑝 1800 J kg -1 K -1 Thermal conductivity 𝑘 0.13 W m -1 K -1 Thermal diffusivity 𝛼 5.6 10 -8 m 2 s -1 three pieces. The thermophysical properties of PLA (Polylactic Acid) are given in Table 1.

To our knowledge, no plastic-core electromagnetic toroidal coupler studies have been reported in the scientific literature to date, either for low-power or high-power applications. Our work is therefore the first thermal preliminary study of a high-power electromagnetic toroidal coupler.

In this paper, we chose to model natural or forced convection with constant heat transfer coefficients. We are aware that this is a strong approximation. But as we did not initially know which technical solution would be used and how the flows would be directed to the surfaces of the torus to be cooled, we chose this solution. Depending on the direction of the flow approaching the external torus surface, different vortex structures may develop downstream, and a vortex shedding with a certain frequency may appear [START_REF] Leweke | The flow behind rings: bluff body wakes without end effects[END_REF][START_REF] Inoue | An experimental study on a wake of a torus using UVP monitor[END_REF]. These spatially complex and unsteady flows will therefore imply a local dependence of the convective coefficients on the azimuthal and poloidal directions (cf. Fig. 5) which have barely been studied to date [START_REF] Yovanovich | Natural convection from horizontal circular and square toroids and equivalent cylinders[END_REF][START_REF] Moshkin | Numerical study of flow and heat transfer from a torus placed in a uniform flow[END_REF][START_REF] Zeitoun | Numerical modeling of natural ocnvection heat transfer from a horizontal torus[END_REF]. Inside the torus, the convective flows will develop in the well-known form of Dean's counter-rotating helical cells [START_REF] Berger | Flow in curved pipes[END_REF], which develop within flat (torsionfree) channels of constant curvature (here R). They will also induce variations in the azimuthal and poloidal directions of the convective transfer coefficients [START_REF] Vashith | A review on the potential applications of curved geometries in process industry[END_REF][START_REF] Ghobadi | A review of heat transfer and pressure drop correlations for laminar flow in curved circular ducts[END_REF]. Thus, due to the complexity of the internal and external flow situations that can exist depending on the flow and convection regimes (natural or forced), in this work we chose to use constant mean convective exchange coefficients for the external exchange surface ℎ ext and also for the internal exchange surface ℎ int .

To study the thermal management of the torus, we first carried out mathematical thermal modeling of the electromagnetic coupler. For simplicity, the heat generated by the Joule effect from the primary and secondary windings is understood as a dissipated heat flux on the boundary of the PLA holder (see Fig. 1). Consequently, the coupler can be considered only as a hollow torus which is subjected to a parietal heat flux (the Joule losses) and to external and internal convective cooling with fluid inlet and outlet.

As previously mentioned, the heat transfer problem in a hollow torus with convective boundary conditions has never been studied in the literature, either analytically or numerically. The only works we found for heat transfer within a torus are those of Wiesche [START_REF] Wiesche | Transient heat conduction in a torus: theory and application[END_REF] Alassar and Abushoshah [START_REF] Alassar | Heat conduction from donuts[END_REF][START_REF] Alassar | Cooling of a hot torus[END_REF] but they concern other types of boundary conditions.

We relied on numerical simulations only to compute the temperature distribution inside the torus. Since we wished to know whether the temperature inside the PLA holder rises above 60 • C, we computed the optimal internal and external convective heat transfer coefficients for which the PLA holder is not deteriorated. The optimization was performed parametrically, meaning that only one parameter (e.g. internal or external temperature, minor or major radius, etc.) varied, while the others remained fixed.

It should be emphasized that the previous model approximates the windings as a surface source term. We did not consider the complete geometry of the electromagnetic coupler. It was therefore of interest to compute the temperature inside the coupler while also considering the primary and secondary windings. To achieve this, numerical simulations using ANSYS-Fluent were performed. We used the optimal value of the external heat transfer coefficients obtained with the first modeling to compare the two approaches.

It is worth remembering that this article proposes a preliminary and simplified study of the thermal management of toroidal couplers. This study gives an easy and quick idea of the maximum temperature of the system for different levels of power losses as well as for different dimensions of the toroidal coupler. Among the main objectives of the project is the continuous optimization of the coupler; for this, other toroidal couplers which are more compact than the one presented in this article are under development. Studying their thermal feasibility is a radical point in our optimization process.

Our paper is organized as follow: we first introduce the characteristics of the toroidal electromagnetic coupler and the electrical experimental tests to identify the heat flux to be cooled. Next, the thermal modeling of the torus subjected to external convective boundary conditions with a dissipated heat flux is proposed. In the third section, the complete geometry with the winding is taken into account. The thermal and 3D numerical models are then compared and discussed. Finally, the paper ends with some conclusions and future avenues for research.

The system's electrical performances

This section is dedicated to the presentation of the characteristics and the electrical performances of the electromagnetic coupler. The operating regime of the system is also explained in full in order to provide a better understanding of the usefulness of studying the thermal management of the toroidal coupler.

Operating regime

Before developing a cooling system for the electromagnetic coupler, the electrical performances of the system must be checked. For this reason, the only way to test the toroidal coupler at a high-power level without damaging it is to use burst mode. In this regime, the system works only for a very short period and stops for a long time before the next cycle. Consequently, the system does not have enough time to heat up.

Fig. 2 shows the waveforms of the input voltage while the system is working in burst regime. The first part of the signal is made up of several pulses (here, 30). The system operates only during this period (0.25% of the time) and no power is transmitted otherwise. This signal is repeated every 0.65 ms.

By adopting this operating regime, the electrical performance of the toroidal coupler can be evaluated. The dissipated heat flux can also be determined experimentally in order to study analytically and/or numerically the temperature rise at high power as well as the various possible cooling solutions. However, to confirm these results through experimental measurements, it is first mandatory to develop a cooling system and to operate it continuously, as shown in Fig. 3.

Characteristics of the toroidal electromagnetic coupler

The electromagnetic coupler's parameters are presented in Table 2. The cable holder is hollow, so the inner radius of the tube is 𝑟 1 , the outer radius of the tube is 𝑟 2 and 𝑅 is the torus radius.

To ensure the compactness of the electromagnetic coupler, a high operating frequency of 200 kHz is used. Working at such a frequency level requires limiting AC losses due to the skin and proximity effects. Therefore, a Litz wire composed of 2500 strands with 71 μm diameters is used according to the operating frequency. At this frequency, the measured AC resistances 𝑅 𝐴𝐶 of the primary and secondary windings of the electromagnetic coupler are 15 mΩ. Fig. 4 presents the experimental setup. The measured average transmitted power is almost 200 kW with an efficiency of 97.7%. The primary current (I𝐼 1𝑟𝑚𝑠 ) is measured to be 400 A and the secondary current(𝐼 2𝑟𝑚𝑠 ) is 373 A. According to:

𝑃 𝑙𝑜𝑠𝑠𝑒𝑠 = 𝑅 𝐴𝐶 𝐼 2 1𝑟𝑚𝑠 + 𝑅 𝐴𝐶 𝐼 2 2𝑟𝑚𝑠 , Fig. 3.
Experimental waveform of the input voltage in continuous regime. The project's ultimate goal is for the electromagnetic coupler to operate continuously throughout the period, so a cooling system is mandatory to avoid the PLA cable holder melting. the primary losses are 2.4 kW and the secondary losses are 2.1 kW, so total AC losses measured for this operating point come to 4.5 kW. Therefore, the thermal study is carried out in a steady state mode considering a heat flux of 4.5 kW to be cooled. Therefore, the thermal study will be carried out in a steady state mode considering a heat flux of 4.5 kW to be cooled.

Cooling a torus: thermal modeling

The first part of this section presents the geometrical parametrization of a hollow torus by the poloidal coordinates. Next, the heat equation with convective boundary conditions is presented. The third part is dedicated to the parametric optimization study giving the optimal heat transfer coefficients for which the PLA holder is not deteriorated. Finally, the numerical results are discussed, including the temperature distribution inside the toroidal coupler.

Geometrical parametrization

The parametrization of a hollow torus is (cf. Fig. 5) as follow :

⎧ ⎪ ⎨ ⎪ ⎩ 𝑥(𝑟, 𝜃, 𝜙) = (𝑅 + 𝑟 cos(𝜃)) cos(𝜙), 𝑦(𝑟, 𝜃, 𝜙) = (𝑅 + 𝑟 cos(𝜃)) sin(𝜙), 𝑧(𝑟, 𝜃, 𝜙) = 𝑟 sin(𝜃), (1)
where 𝑅 is the radius of the tore, 𝑟 1 ≤ 𝑟 ≤ 𝑟 2 is the radius of the tube, 𝑅∕𝑟 is the aspect ratio, 𝜃 and 𝜙 are the angles in the poloidal and toroidal directions, respectively, with (𝜃, 𝜙) ∈ [0, 2𝜋] × [0, 2𝜋]. For a torus with minor radius 𝑟, its surface is

𝑆 𝑡𝑜𝑟𝑢𝑠 = 4𝜋 2 𝑅𝑟, (2)
and its volume is

𝑉 𝑡𝑜𝑟𝑢𝑠 = 2𝜋 2 𝑅𝑟 2 . ( 3
)
For a hollow torus with inner and outer radius 𝑟 1 and 𝑟 2 , its volume can be computed from (3) and reads

𝑉 ℎ𝑜𝑙𝑙𝑜𝑤-𝑡𝑜𝑟𝑢𝑠 = 2𝜋 2 𝑅 ( 𝑟 2 2 -𝑟 2 1
) .

(4)

Heat equation with convective boundary conditions

The transient heat conduction equation to solve is:

𝛥𝑇 = 1 𝛼 𝜕𝑇 𝜕𝑡 ( 5
)
where 𝑇 is the temperature, 𝛼 is the thermal diffusivity of the material, and 𝑡 is the time. We first rewrite the heat equation in the (𝑟, 𝜃, 𝜙) coordinates. The metric tensor is defined as 𝑔 = 𝐽 𝑇 𝐽 where 𝐽 is the Jacobian matrix of the mapping (𝑟, 𝜃, 𝜙) ↦ (𝑥, 𝑦, 𝑧). Some computations give

𝑔 = ⎛ ⎜ ⎜ ⎝ 1 0 0 0 𝑟 2 0 0 0 (𝑅 + 𝑟 cos(𝜃)) 2 ⎞ ⎟ ⎟ ⎠ , 𝑔 -1 = ⎛ ⎜ ⎜ ⎜ ⎝ 1 0 0 0 1 𝑟 2 0 0 0 1 (𝑅+𝑟 cos(𝜃)) 2 ⎞ ⎟ ⎟ ⎟ ⎠ .
As a result, denoting (𝑔 𝑖𝑗 ) 𝑖𝑗 = 𝑔 -1 , the Laplace operator in the moving frame reads (we use the summation over repeated indices and refer to [START_REF] Fecko | Differential Geometry and Lie Groups for Physicists[END_REF] for the differential operators associated with a metric tensor)

Δ𝑓 = 1 √ 𝑔 𝜕 𝑗 (√ 𝑔𝑔 𝑗𝑘 𝜕 𝑘 𝑓 ) = 1 √ 𝑔 ( 𝜕 𝑟 (√ 𝑔𝜕 𝑟 𝑓 ) + 𝜕 𝜃 ( √ 𝑔 𝑟 2 𝜕 𝜃 𝑓 ) + 𝜕 𝜙 ( √ 𝑔 (𝑅 + 𝑟 cos(𝜃)) 2 𝜕 𝜙 𝑓 ))
,

where √ 𝑔 = 𝑟(𝑅 + 𝑟 cos(𝜃)) is the square root of the determinant of the metric tensor. The transient heat equation ( 5) is then transformed into (𝑟, 𝜃, 𝜙) coordinates as follows

Δ T = 1 𝛼 𝜕 T 𝜕𝑡 , (6)
where T (𝑡, 𝑟, 𝜃, 𝜙) = 𝑇 (𝑡, 𝑥, 𝑦, 𝑧) is the temperature in the (𝑟, 𝜃, 𝜙) coordinates.

We consider a hollow torus of respectively inner radius 𝑟 1 and outer radius 𝑟 2 of the toroidal tube. The torus is to be air or watercooled, by convection (natural or forced) inside the torus (ℎ int ) and air-cooled outside the torus (ℎ ext ) (see Fig. 5). As a first approach, we will consider the constant convective heat exchange coefficients as explained above. The temperatures of the fluid inside and outside the torus are respectively 𝑇 int and ext . We recall that the Litz wires are considered as a surface heat source. The heat flux 𝑄 dissipated by the Litz wires is considered to be homogeneously distributed at the outer surface of the toroidal tube 𝑆 ext , where 𝑆 ext = 4𝜋 2 𝑟 2 𝑅 (see [START_REF] Boudara | Design and development of high efficiency 150 kw very compact PLA core electromagnetic coupler for highly resonant power transfer technology[END_REF]). Taking into account the convective heat flux with the Fourier law, the boundary conditions finally reads as follows:

𝑟 = 𝑟 2 ∶ -𝑘 𝜕 T 𝜕𝑟 | | | | |𝑟=𝑟 2 = - 𝑄 𝑆 ext + ℎ ext ( T | 𝑟=𝑟 2 -𝑇 ext ) , ( 7
)
𝑟 = 𝑟 1 ∶ 𝑘 𝜕 T 𝜕𝑟 | | | | |𝑟=𝑟 1 = ℎ int ( T | 𝑟=𝑟 1 -𝑇 int ) ,
where 𝑘 is the thermal conductivity of the material. We now give the dimensionless form of the previous heat equation. To do so, we introduce the aspect ratio of the torus 𝜂 and the dimensionless quantities

𝜂 = 𝑅 𝑏 , 𝑡 * = 𝑡 𝑡 0 = 𝛼𝑡 𝑏 2 , 𝑟 * = 𝑟 𝑏 ,
where 𝑡 0 and 𝑏 are arbitrary. Using the notation 𝑓 * (𝑡 * , 𝑟 * , 𝜃, 𝜙) = 𝑓 (𝑡∕𝑡 0 , 𝑟∕𝑏, 𝜃, 𝜙), we get

𝜕𝑓 * 𝜕𝑟 * = 𝑏 𝜕𝑓 𝜕𝑟 , 𝜕𝑓 * 𝜕𝑡 * = 𝑏 2 𝛼 𝜕𝑓 𝜕𝑡 .
The dimensionless Laplace operator now reads

Δ * 𝑓 * = 1 𝑏 2 ( 𝜕 2 𝑓 * 𝜕𝑟 * 2 + 1 (𝑟 * ) 2 𝜕 2 𝑓 * 𝜕𝜃 + 1 (𝜂 + 𝑟 * cos(𝜃)) 2 𝜕 2 𝑓 * 𝜕𝜙 2 ) + 1 𝑏 2 ( 2𝑟 * cos(𝜃) + 𝜂 𝑟 * (𝜂 + 𝑟 * cos(𝜃)) 𝜕𝑓 * 𝜕𝑟 * - 1 (𝑟 * ) 2 𝑟 * sin(𝜃) (𝜂 + 𝑟 * cos(𝜃)) 𝜕𝑓 * 𝜕𝜃
) .

We also introduce the notation Δ * 𝑓 * = 1 𝑏 2  * [𝑓 * ] and the following dimensionless temperature

𝑇 * (𝑡 * , 𝑟 * , 𝜃, 𝜙) = 𝑇 ( 𝑡 𝑡 0 , 𝑟 𝑏 , 𝜃, 𝜙 ) -𝑇 ∞ 𝑇 ref ,
where 𝑇 ∞ and 𝑇 ref are arbitrary. The dimensionless heat equation is then

𝜕𝑇 * 𝜕𝑡 * =  * [𝑇 * ].
Denoting 𝑟 * 1 = 𝑟 1 ∕𝑏 and 𝑟 * 2 = 𝑟 2 ∕𝑏, the boundary conditions are recast into

𝑟 * = 𝑟 * 2 ∶ - 𝜕𝑇 * 𝜕𝑟 * | | | |𝑟 * =𝑟 * 2 = -𝑄 * + 𝐵𝑖 ext ( 𝑇 * | 𝑟 * =𝑟 * 2 + 𝑇 ∞ -𝑇 ext 𝑇 ref ) , 𝑟 * = 𝑟 * 1 ∶ 𝜕𝑇 * 𝜕𝑟 * | | | |𝑟 * =𝑟 * 1 = 𝐵𝑖 int ( 𝑇 * | 𝑟 * =𝑟 * 1 + 𝑇 ∞ -𝑇 int 𝑇 ref ) ,
where

𝑄 * = 𝑄 𝑏 𝑘𝑇 ref 𝑆 ext
is the dimensionless surface heat source and

𝐵𝑖 int = ℎ int 𝑏 𝑘 , 𝐵𝑖 ext = ℎ ext 𝑏 𝑘
are the inner and outer Biot numbers. The latter, for example, characterizes the ratio of the heat transfer resistances, conductive inside the body and convective at the external boundary of the torus. Above 𝑏 is a characteristic length defined by the ratio of the volume to the surface of the body at the surface in question, namely 𝑟 = 𝑟 2 :

𝑏 = 1 2𝜋 𝑟 2 2 -𝑟 2 1 𝑟 2 . ( 8
)
We recall that the temperature inside the torus does not depend on the spatial variable (i.e. lumped analysis) if 𝐵𝑖 < 0.1. In our case, this yields ℎ < 0.1 𝑘 𝑏 = 0.664, which indicates that the lumped analysis does not apply to our torus when natural convection or forced convection is involved. We thus rely on numerical simulations to get the temperature inside the torus.

Numerical method

Recall that the electromagnetic coupler is intended to operate continuously. The coupler reaches its maximal temperature in the steadystate regime. Since this is the most difficult case related to cooling, we are now working in the steady regime. We can consider that the heat flux that emanates from the torus wall is transferred to the cooling fluid that flows inside the torus. Thus, if we consider the assumption of a steady flow, along a differential element 𝑑𝑥 of the inner torus channel, the fluid will undergo a temperature variation 𝑑𝑇 𝑚 which can be calculated by knowing the heat flow transferred to the wall 𝑑𝑞: 𝑑𝑞 = 𝑚 𝐶 𝑝 𝑑𝑇 𝑚 = 𝑞 ′′ 𝑠 𝑃 𝑑𝑥 where 𝑚 is the constant flow rate of the fluid, 𝐶 𝑝 its specific heat, 𝑃 is the inner surface perimeter of the torus and 𝑞 ′′ 𝑠 is the inner surface heat flux density.

Since no explicit solution can be found in the setting considered, numerical simulation is then used to compute the temperature inside the torus slice (a differential element 𝑑𝑥). The conduction problem is solved using finite elements and Free-Fem ++ [START_REF] Hecht | New development in freefem++[END_REF]. Let  be a quasi-uniform triangular mesh of [0, 1] × [0, 2𝜋] and

X 𝜏 = { 𝜓 𝜏 is continuous | 𝜓 𝜏 | 𝑇 ∈ P 1 and 𝜃 ↦ 𝜓 𝜏 (𝑟, 𝜃) 2𝜋-periodic } .
Let 𝜏 𝑇 be the diameter of its inscribed circle. The meshsize 𝜏 is then defined as 𝜏 = max 𝑇 ∈ 𝜏 𝑇 . The variational formulation of ( 6)-( 7) is obtained by multiplying (6) with a test function √ 𝑔𝜑, integrating by part and then using the boundary condition [START_REF] Moshkin | Numerical study of flow and heat transfer from a torus placed in a uniform flow[END_REF] (see [17, Chapter IX]). The discretized variational formulation of ( 6)-( 7) then reads:

Find 𝑇 𝜏 ∈ X 𝜏 satisfying ∫ 1 𝑟=0 ∫ 2𝜋 𝜃=0 ( 𝜕 𝑟 𝑇 𝜏 𝜕 𝑟 𝜑 𝜏 + 1 𝑟 2 𝜕 𝜃 𝑇 𝜏 𝜕 𝜃 𝜑 𝜏 ) √ 𝑔 𝑑𝑟 𝑑𝜃 + ∫ 𝜃=2𝜋 𝜃=0 ℎ ext 𝑘 (√ 𝑔𝑇 𝜏 𝜑 𝜏 ) | | |𝑟=𝑟 2 + ℎ int 𝑘 (√ 𝑔𝑇 𝜏 𝜑 𝜏 ) | | |𝑟=𝑟 1 𝑑𝜃 (9) -∫ 2𝜋 𝜃=0 (√ 𝑔𝜑 𝜏 ) | | |𝑟=𝑟 2 ( 𝑄 𝑘𝑆 ext + ℎ ext 𝑘 𝑇 ext ) + (√ 𝑔𝜑 𝜏 ) | | |𝑟=𝑟 1 ℎ int 𝑘 𝑇 int 𝑑𝜃 = ∫ 𝑟 2 𝑟=𝑟 1 ∫ 𝜃 √ 𝑔𝐹 𝜑 𝜏 𝑑𝑟 𝑑𝜃, ∀𝜑 𝜏 ∈ X 𝜏 .
where 𝐹 is a volume source term set to zero in the following numerical results.

Parametric optimization in the steady regime

First, we show that cooling the torus so that the PLA is not deteriorated is mandatory. For this purpose, the external and internal heat transfer coefficients are chosen equal to be equal 25 W m -2 K -1 . Fig. 6 presents the temperature rise at the level of the tube and not the torus, it represents a cross-section of the coupler's loop. For the sake of simplicity, only one half of the tube is presented as there is a symmetry between the top and the bottom of the tube and also because the convective heat exchange coefficients are considered constant. From this figure, it can be noticed that the temperature is lower at the inner radius of the loop and reaches its maximal value on the outer surface, namely where the flux dissipated by the wires is involved.

Concerning the temperature rise with natural convection, the maximum reached temperature reached is 758 K (485.55 • C). It is obvious that at such a level of heat losses, the temperature rise will be huge and the electromagnetic coupler will need to be cooled down.

Fig. 6 shows the evolution of the temperature at different places. The red curve is when 𝜃 = 0, which is at the extremity of the loop, as seen in Fig. 6. The green curve is for 𝜃 = 𝜋∕3 and the blue one is for 𝜃 = 2𝜋∕3. This clearly shows that the temperature distribution is not the same at all points of the loop since, in poloidal coordinates, there is no symmetry according to 𝜃. We emphasize that this is the main difference between resolving the system in cylindrical coordinates and in poloidal ones.

In Fig. 6, we can see that the electromagnetic coupler needs larger ℎ to avoid destruction of the material. Next, we numerically compute the minimal value for the external convective heat transfer coefficient ℎ ext for which the PLA does not melt. Let ℎ min ext be the minimal value of the external heat transfer coefficient such that the maximal temperature of the torus does not exceed 60 • C. The latter is thus defined as

ℎ min ext = min { ℎ ext such that max 𝑇 < 60 • C } ,
where T is the temperature computed from the heat equation ( 6)- [START_REF] Moshkin | Numerical study of flow and heat transfer from a torus placed in a uniform flow[END_REF]. Before going further into the parametric optimization, it is worth noting that ( 6)-( 7) can have closed-form solutions when

ℎ int = 0.
Indeed, the only solution to ( 6)-( 7) is the constant temperature 𝑇 , satisfying

0 = - 𝑄 𝑆 ext + ℎ ext ( 𝑇 -𝑇 ext ) , Table 3
Investigated configurations for cooling a torus.

𝑇 ext [K] 𝑇 int [K] 𝑟 1 [mm] 𝑟 2 [mm] 𝑅 [mm] ℎ int [W m -2 K -1 ]
Case which gives

𝑇 = 𝑄 𝑆 ext ℎ ext + 𝑇 ext .
The value of ℎ min ext is thus

ℎ min ext = 𝑄 𝑆 ext (60 -𝑇 ext ) . ( 10
)
We now compute the value of ℎ min ext when ℎ int ≠ 0 in the settings given in Table 3.

All the following numerical experiments have been done with the same meshsize 𝜏 ≈ 0.0141 mm and the numerical results can be found in Fig. 7 for Cases 1-2 and Fig. 8 for Cases 3-4.

From our numerical results, we may draw the following conclusions. First, the value of ℎ min ext increases with 𝑄 and with the exterior and the interior temperatures (cf. Fig. 7 -Left). In addition, as seen in Fig. 7 (Right), ℎ min ext is not affected by the value of the interior heat transfer coefficient, which indicates that it is not worth trying to cool the torus from the inside. We expect that this result to stem from two mains factors. First, the thermal conductivity of the PLA material is small, hence this material is a thermal insulator. In addition, the power is dissipated on the outer surface of the torus. In fact, all the cooling work is done by the external heat transfer coefficient and increasing the internal heat transfer coefficient hardly affects the cooling.

Cooling the torus with winding

In this section, we perform numerical simulations that take into account the wires to verify the mathematical model. The steady-state thermal analysis of the electromagnetic coupler was carried out using the commercial software ANSYS-Fluent. The geometry simulated includes the hollow PLA cable holder and the primary and secondary Litz wire windings. The diameter of a single strand of Litz wire in highfrequency applications is very small (a few tens of micrometers). Taking into account the depth of the skin effect requires substantial meshing and is time-consuming. Therefore, the method proposed by [START_REF] Etemadrezaei | Equivalent complex permeability and conductivity of litz wire in wireless power transfer systems[END_REF] to simulate the Litz wire is used. It consists in homogenizing the whole Litz bundle by finding its complex quantities. Each strand is replaced by an equivalent represented by complex permeability and conductivity and then the whole bundle is homogenized with its equivalent complex quantities. The software takes as heat sources a primary power loss of 2.3 kW and a secondary power loss of 2.2 kW. In order to simplify these numerical simulations, the following assumptions are made:

• The ambient temperature is constant and is equal to 20 • C. It is also the initial temperature chosen for the simulation. • There is no influence of the temperature rise on the thermal properties of materials (PLA and copper). • Heat losses are constant and uniformly distributed in each turn of the Litz wires, which implies that the heat flux density is no longer homogeneous, as assumed in the previous study, but now depends on 𝜃 and 𝜙.

Firstly, it is mandatory to determine the temperature rise in the system with natural convection. The cable holder and the windings have been meshed by tetrahedral volume elements generated automatically, then refined, adapted and optimized in order to obtain the best mesh which gives a good compromise between precision of results and calculation time. The mesh has 165,658 nodes and 80,163 elements, and a time-accurate thermal analysis is used for the solution procedure.

To obtain the following results, the Euler-explicit method was used with a step time of 0.5 s. Additionally, 8 h of CPU-time were needed using a PC-server with 8 processors (3.6 GHz processor and 64 GB RAM).

The temperature distribution is shown in Fig. 10 for a convective heat transfer coefficient of 25 W m -2 K -1 . The highest temperature in the electromagnetic coupler is the temperature of the winding, which is equal to 471.37 • C. From the mathematical model presented in Section 3.3, the temperature rise with natural convection was found to be 484.55 • C which means, that the difference between the two methods is only 2.8%.

The temperature reached is much higher than the limit temperature supported by the PLA material. As indicated in Fig. 6, it is therefore mandatory to introduce an efficient cooling system in order to avoid exceeding the operating temperature limit. Therefore, a convective heat transfer coefficient needs to be determined to consider the effect of the fluid on the solid. Generally, this coefficient can be provided through experimental tests. In our case, this coefficient is determined by using the approach in Section 3. According to the dimensions of the coupler and for heat losses of 4.5 kW, the calculated value for the external heat transfer coefficient is 305 W m -2 K -1 for the primary and secondary windings and the PLA cable holder.

Fig. 11 represents the temperature evolution of the coupler as a function of Joule losses for two values of the heat transfer coefficient. The red curve is the temperature rise obtained by using the mathematical model and the green curve is the result of the full numerical modeling. The results are very consistent: the difference between the two methods varies between 1.25% for low losses and 9.5% maximum for significant ones.

To enable the continuous functioning of the coupler without developing a cooling solution, maximum losses should not exceed 300 W, which means that only 1.3 kW can be transmitted from the primary to the secondary windings. This is at odds with the goal and ultimate application of the system. For this reason, the experimental test under the normal condition at a high-power level could not be conducted, mainly in order not to destroy the electromagnetic coupler and the power converters located upstream and downstream of the coupler. Using a cooling system (with h = 305 W m -2 K -1 ) will enable the coupler to work at high power, up to 200 kW, without exceeding the PLA temperature limit.

Fig. 12 represents the temperature field for the electromagnetic coupler. The maximum temperature reached in the primary and secondary windings is 57.36 • C, and in the PLA cable holder 38.08 • C. The heating of the coils is more intense than of the holder given their high thermal conductivity. In addition, the thermal limits of the different system components are respected.

Fig. 13 represents the temperature rise of the coupler's loop with forced convection using the simplified model from Section 3. For the same convective heat transfer coefficient, the maximum temperature reached is 332 𝐾 (59.95 • C). Although this value is equal to the limit temperature supported by the PLA, it does not present any danger of destroying the system.

Conclusion and perspectives

This paper presented two numerical approaches in order to study the cooling of a high-power toroidal electromagnetic coupler. The first methodology adopted estimates an approximate temperature distribution of the system by solving a 2D heat equation in poloidal coordinates and numerically optimizing the heat transfer coefficient without taking the Litz windings into account. The second approach considers the complex geometry of the coupler with windings and relies on 3D numerical modeling.

Both methods give results that are in good agreement. The maximum difference between the two methods used is about 9.5% only for the chosen parameters. This verifies the model presented in this paper and shows that it can be trusted when designing an electromagnetic coupler, with the consequence that complex, expensive CAD and CFD software, with time-consuming calculations, is not really required.

A CFD-type modeling will, however, be necessary when the choice of internal and/or external cooling solutions for the toroidal electromagnetic coupler is made. This numerical study allowed us to obtain the local convective exchange coefficients under cooling conditions. Furthermore, once the cooling solution is developed and the electromagnetic coupler is able to operate under normal conditions, experimental tests to evaluate the temperature rise with a high-power level will be performed, in order to validate the two developed methods.

Moreover, the parametric study showed that is very difficult to cool down the electromagnetic coupler using PLA material due to its low thermal conductivity and temperature limit. We believe that exchanging the PLA material for another, higher thermal conductive material (such as TC Poly) will allow the coupler to be cooled from inside, which will be a more effective and innovative solution.

Appendix. Code verification for the thermal modeling

We proceed here with a validation of the code used to compute the value of ℎ min ext To do this, we consider the following manufactured solution One can then check that 𝑇 exact is the exact solution to ( 6)-( 7) with 𝐹 = -(2𝑟 cos(𝜃) + 𝑅)𝑎 (𝑟(𝑅 + 𝑟 cos(𝜃)))

𝑇
as source term. The 𝐿 2 error defined as

√ ∫ 1 𝑟=0 ∫ 2𝜋 𝜃=0 ( T -𝑇 exact ) 2 √ 𝑔 𝑑𝑟 𝑑𝜃,
can be found in Table 5. These results show the convergence of the numerical solution toward the exact one, thus validating the code. The numerical simulations in Section 3.3 were all done using Free-Fem ++ software [START_REF] Hecht | New development in freefem++[END_REF] available at www.freefem.org. The code can be found at https://osur-devspot.univ-reunion.fr/phcocque/conduction_ torus.

Fig. 2 .

 2 Fig. 2. Experimental waveform of the input voltage in burst regime.

Fig. 4 .

 4 Fig. 4. The experimental setup. The DC source supplies the system with DC voltage, then the power converter (developed by our industrial partner EFFITECH) converts the DC voltage to AC voltage at the input of the electromagnetic coupler. The resonant capacitors are used to compensate for the inductive effect of the primary side and enhance the power capability transfer. The current and voltage probes are used to measure the input and output current and voltage, enabling us to determine the average input and output power. The scope is used to show the waveforms of the input and output current and voltage signals. A resistive load of 1 Ω is used at the output of the electromagnetic coupler.

Fig. 5 .

 5 Fig. 5. Parametrization of a torus by the coordinates (𝑟, 𝜃, 𝜙).

Fig. 6 .

 6 Fig. 6. The temperature rise of the coupler's loop as function of r (mm) with 𝑄 = 5000 W, ℎ ext = ℎ int = 25 W m -2 K -1 with 𝑇 ext = 20 • C. Left: Temperature field. The center of the torus is located on the left. Right: Cut along some angles, 𝜃 = 0, 𝜃 = 𝜋∕3 and 𝜃 = 2𝜋∕3.

Fig. 7 .

 7 Fig. 7. Value of ℎ min ext for the numerical experiments from Table 3 for different values of the heat flux 𝑄 (W) dissipated by the Litz wires. Left: Case 1, Right: Case 2.

Fig. 11 .

 11 Fig. 11. Temperature's evolution of the electromagnetic coupler as a function of Joule losses and for natural convection with h = 25 W m -2 K -1 and forced convection with h = 305 W m -2 K -1 .

Fig. 12 .

 12 Fig. 12. The temperature rise of the coupler with forced convection (ℎ ext = 305 W m -2 K -1 ) with 𝑇 ext = 20 • C.

Fig. 13 .

 13 Fig. 13. The temperature rise of the coupler's loop as function of r (mm) with forced convection (ℎ int = 25, ℎ ext = 305 W m -2 K -1 ) with 𝑇 ext = 𝑇 int = 20 • C.

Table 2

 2 Parameters of the electromagnetic coupler[START_REF] Boudara | Design and development of high efficiency 150 kw very compact PLA core electromagnetic coupler for highly resonant power transfer technology[END_REF].

	Parameter	Value	Unit
	𝑅	124	mm
	𝑟 2	52.25	mm
	𝑟 1	26.12	mm
	Frequency	200	kHz
	𝑅 𝐴𝐶	15	mΩ

Table 5

 5 𝐿 2 error between the numerical and the explicit solution for varying meshsize.

	Meshsize 𝜏	0.1414	0.02828	3.06204 × 10 -5
	𝐿 2 -error	0.148168	3.06204 × 10 -5	4.82811 × 10 -7

  exact (𝑟, 𝜃) = 𝑎𝑟 + 𝑏,with 𝑎 = -ℎ int (𝑆 ext 𝑇 ext ℎ ext -𝑆 ext 𝑇 int ℎ ext + 𝑄) 𝑆 ext (ℎ ext ℎ int 𝑟 1 -ℎ ext ℎ int 𝑟 2 -ℎ ext 𝑘 -ℎ int 𝑘) 𝑏 = 𝑆 ext 𝑇 ext ℎ ext ℎ int 𝑟 1 -𝑆 ext 𝑇 int ℎ ext ℎ int 𝑟 2 -𝑆 ext 𝑇 ext ℎ ext 𝑘 -𝑆 ext 𝑇 int ℎ int 𝑘 + 𝑄ℎ int 𝑟 1 -𝑄𝑘 𝑆 ext (ℎ ext ℎ int 𝑟 1 -ℎ ext ℎ int 𝑟 2 -ℎ ext 𝑘 -ℎ int 𝑘) .

Acknowledgments

The authors are grateful to the French Ministry for Defense (Direction Générale de l'Armement) for its financial support for this work (RAPID CEMPER N • 172906093).

The authors would also like to thank members of the French company EFFITECH for their collaboration.

In Fig. 8, we notice that ℎ min ext decreases as the minor and major radius of the electromagnetic coupler increases, which is obvious because the larger the exchange surface, the easier the cooling. In both cases, when the radius (𝑟 2 or 𝑅) is doubled, the value of the convective exchange coefficient ℎ min ext required to maintain the PLA electromagnetic coupler at an acceptable temperature is halved (from a value of 400 W m -2 K -1 to a value close to 200 W m -2 K -1 ).

We conclude this parametric optimization by computing the percentage error between the numerically optimized ℎ min ext and the explicit one from [START_REF] Vashith | A review on the potential applications of curved geometries in process industry[END_REF], namely where the max is taken on each varying parameter according to the case considered. The results are shown in Table 4. Therefore, although the temperature distributions inside the torus for ℎ int = 0 (adiabatic wall) and ℎ int > 0 are very different, the numerical values of ℎ min ext are actually very closed to the explicit value (10) except, for Cases 2-4.

It is worth noting that the explicit formula [START_REF] Vashith | A review on the potential applications of curved geometries in process industry[END_REF] does not depend on 𝑘, 𝑟 1 nor ℎ int and we thus now study the influence of these two remaining parameters on ℎ min ext . Our numerical results are found with

with varying ℎ int from 25 to 100 W m -2 K -1 .

The numerical results can be found in Fig. 9. Both values diverge from each other if the inner heat transfer coefficient and the conductivity are larger and, in such cases, a smaller ℎ min ext is obtained for the torus cooling.

From our numerical investigations we can conclude that formula (10) can be used to get ℎ min ext if the inner heat transfer coefficient is small enough, and that in the other cases the choice of cooling by the interior of the torus makes it possible to reduce the value of the external heat transfer coefficient in a non-negligible way.

The heat transfer exchange coefficients for free or forced convection (internal and external) chosen correspond to realistic cases that can be obtained from correlations of the evolution of the Nusselt number for a flow of air outside the torus or for a flow of air or liquid inside the torus (case of a curved channel with a circular section) [START_REF] Vashith | A review on the potential applications of curved geometries in process industry[END_REF].

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Lucile Blouquy

Claire Mottet

Jérôme Olivares

Christophe Plantamp

Myriam Siegwart

Benoit Barrès
email: benoit.barres@anses.fr

Je ´ro ˆme Olivares

How varying parameters impact insecticide resistance bioassay: An example on the worldwide invasive pest Drosophila suzukii

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

Introduction

Bioassays are considered as the gold standard for pesticide resistance testing. They can detect new resistances and can assess resistance in an integrative manner, whether the underlying mechanism is known or unknown [START_REF] R4p | Trends and challenges in pesticide resistance detection[END_REF]. With the growing ambition to employ sustainable practices of pest control, monitoring pesticide resistance is becoming the key to implementing efficient, integrated pest management strategies. Accurate and reliable phenotypic data, in addition to molecular data, are also essential to address numerous theoretical evolutionary questions linked to pesticide resistance (see [START_REF] Hawkins | The evolutionary origins of pesticide resistance[END_REF] for examples with insecticides). As in most experiments involving living material, the output of pesticide resistance bioassays can be affected by many parameters, depending on the type of bioassay used and the study species.

When dealing with a new species/active substance pair, dose-response analyses are the most frequently used approach. The standard procedure consists in exposing insects to a range of doses of the active substance of interest and in counting the number of dead individuals for each dose. Using regression methods, it is then possible to determine various parameters of interest, and in particular the median lethal dose (LD 50 ), the dose that kills 50% of the exposed individuals. Results of bioassays may be affected by experimental conditions and design [START_REF] Owusu | Comparability between insecticide resistance bioassays for mosquito vectors: time to review current methodology?[END_REF] and also by the biological characteristics of the insect species [START_REF] Hamby | Integrating circadian activity and gene expression profiles to predict chronotoxicity of Drosophila suzukii response to insecticides[END_REF][START_REF] Cotter | Macronutrient balance mediates trade-offs between immune function and life history traits[END_REF][START_REF] Kulma | Effects of age and larval nutrition on phenotypic expression of insecticide-resistance in Anopheles mosquitoes[END_REF][START_REF] Chouaibou | Increase in susceptibility to insecticides with aging of wild Anopheles gambiae mosquitoes from Co ˆte d'Ivoire[END_REF][START_REF] Ong | Larval Age and Nutrition Affect the Susceptibility of Culex quinquefasciatus (Diptera: Culicidae) to Temephos[END_REF]. For each pest, it is necessary to develop, or at least to adapt, a standard experimental protocol to obtain reliable and repeatable results that will allow comparisons of pesticide susceptibility over space and time.

In this study, we used the spotted wing drosophila Drosophila suzukii (Matsumura) for a case study. It is a harmful pest that severely affects fruit production and causes severe economic losses [START_REF] Goodhue | Spotted wing drosophila infestation of California strawberries and raspberries: economic analysis of potential revenue losses and control costs[END_REF][START_REF] Lee | In Focus: Spotted wing drosophila, Drosophila suzukii, across perspectives[END_REF][START_REF] Digiacomo | Economic Impact of Spotted Wing Drosophila (Diptera: Drosophilidae) Yield Loss on Minnesota Raspberry Farms: A Grower Survey[END_REF] by altering berries and other soft-skinned fruits. D. suzukii, native to Asia, is also known as a worldwide invasive species [START_REF] Santos | Global potential distribution of Drosophila suzukii (Diptera, Drosophilidae)[END_REF][START_REF] Fraimout | Accounting for intraspecific variation to quantify niche dynamics along the invasion routes of Drosophila suzukii[END_REF] and was first introduced in North America and Europe in 2008 before spreading throughout the globe [START_REF] Cini | A review of the invasion of Drosophila suzukii in Europe and a draft research agenda for integrated pest management[END_REF][START_REF] Fraimout | Deciphering the routes of invasion of Drosophila suzukii by means of ABC random forest[END_REF][START_REF] Rota-Stabelli | Distinct genotypes and phenotypes in European and American strains of Drosophila suzukii: implications for biology and management of an invasive organism[END_REF][START_REF] Asplen | Invasion biology of spotted wing Drosophila (Drosophila suzukii): a global perspective and future priorities[END_REF][START_REF] Nagy | Distribution of the spotted-wing drosophila (Drosophila suzukii) in the north-eastern part of the Carpathian lowlands[END_REF]. The main control strategies for D. suzukii rely on chemical insecticides [START_REF] Haye | Current SWD IPM tactics and their practical implementation in fruit crops across different regions around the world[END_REF][START_REF] Mermer | Canopy spray deposition and related mortality impacts of commonly used insecticides on Drosophila suzukii Matsumura (Diptera: Drosophilidae) populations in blueberry[END_REF] and marginally on alternative control methods [START_REF] Cattel | Back and forth Wolbachia transfers reveal efficient strains to control spotted wing drosophila populations[END_REF][START_REF] Knoll | Seasonal and regional presence of hymenopteran parasitoids of Drosophila in Switzerland and their ability to parasitize the invasive Drosophila suzukii[END_REF][START_REF] Kruger | Effects of Irradiation Dose on Sterility Induction and Quality Parameters of Drosophila suzukii (Diptera: Drosophilidae)[END_REF][START_REF] Schetelig | Environmentally sustainable pest control options for Drosophila suzukii[END_REF][START_REF] Cahenzli | Screening of 25 different natural crop protection products against Drosophila suzukii[END_REF][START_REF] Fanning | Efficacy of biopesticides on spotted wing drosophila, Drosophila suzukii Matsumura in fall red raspberries[END_REF][START_REF] Alnajjar | Behavioral and preventative management of Drosophila suzukii Matsumura (Diptera: Drosophilidae) in Maine wild blueberry (Vaccinium angustifolium Aiton) through attract and kill trapping and insect exclusion-netting[END_REF][START_REF] Dam | Natural compounds for controlling Drosophila suzukii. A review[END_REF][START_REF] Larson | Controlled-Release Dispenser and Dry Trap Developments for Drosophila suzukii Detection[END_REF][START_REF] Scho ¨neberg | Pruning of small fruit crops can affect habitat suitability for Drosophila suzukii[END_REF][START_REF] Finetti | Modulation of Drosophila suzukii type 1 tyramine receptor (DsTAR1) by monoterpenes: a potential new target for next generation biopesticides[END_REF][START_REF] De Souza | Chemical composition of essential oils of selected species of Piper and their insecticidal activity against Drosophila suzukii and Trichopria anastrephae[END_REF][START_REF] Foye | A Rare, Recently Discovered Nematode, Oscheius onirici (Rhabditida: Rhabditidae), Kills Drosophila suzukii (Diptera: Drosophilidae) Within Fruit[END_REF][START_REF] Nikolouli | Combining sterile and incompatible insect techniques for the population suppression of Drosophila suzukii[END_REF][START_REF] Sial | Evaluation of organic insecticides for management of spotted-wing drosophila (Drosophila suzukii) in berry crops[END_REF]. Its dispersal and reproductive abilities make it a species with high evolutionary potential [START_REF] Asplen | Invasion biology of spotted wing Drosophila (Drosophila suzukii): a global perspective and future priorities[END_REF][START_REF] Tochen | Temperature-related development and population parameters for Drosophila suzukii (Diptera: Drosophilidae) on cherry and blueberry[END_REF], and it may be prone to evolve insecticide resistance [START_REF] Gress | Identification and risk assessment of spinosad resistance in a California population of Drosophila suzukii[END_REF]. Moreover, other species of the Drosophila genus have already demonstrated their ability to resist insecticides [START_REF] Perry | The biology of insecticidal activity and resistance[END_REF][START_REF] Wilson | Resistance of Drosophila to toxins[END_REF][START_REF] Wilson | Drosophila: sentinels of environmental toxicants[END_REF][START_REF] Rane | Detoxification genes differ between cactus-, fruit-, and flower-feeding Drosophila[END_REF]. All these reasons call for a careful and large-scale plan to monitor insecticide resistance in this species, but such monitoring plans require reliable and standardized methods. To date, studies have evaluated pesticide efficiency or resistance in D. suzukii populations using various protocols [START_REF] Hamby | Integrating circadian activity and gene expression profiles to predict chronotoxicity of Drosophila suzukii response to insecticides[END_REF][START_REF] Mermer | Canopy spray deposition and related mortality impacts of commonly used insecticides on Drosophila suzukii Matsumura (Diptera: Drosophilidae) populations in blueberry[END_REF][START_REF] Gress | Identification and risk assessment of spinosad resistance in a California population of Drosophila suzukii[END_REF][START_REF] Andreazza | Toxicity to and egg-laying avoidance of Drosophila suzukii (Diptera: Drosophilidae) caused by an old alternative inorganic insecticide preparation[END_REF][START_REF] Beers | Developing Drosophila suzukii management programs for sweet cherry in the western United States[END_REF][START_REF] Bruck | Laboratory and field comparisons of insecticides to reduce infestation of Drosophila suzukii in berry crops[END_REF][START_REF] Cuthbertson | Preliminary Screening of Potential Control Products against Drosophila suzukii[END_REF][START_REF] Mishra | High throughput sequencing reveals Drosophila suzukii responses to insecticides[END_REF][START_REF] Shawer | Laboratory and field trials to identify effective chemical control strategies for integrated management of Drosophila suzukii in European cherry orchards[END_REF][START_REF] Smirle | Laboratory studies of insecticide efficacy and resistance in Drosophila suzukii (Matsumura) (Diptera: Drosophilidae) populations from British Columbia, Canada[END_REF][START_REF] Wise | Curative activity contributes to control of spotted-wing drosophila (Diptera: Drosophilidae) and blueberry maggot (Diptera: Tephritidae) in highbush blueberry[END_REF][START_REF] Van Timmeren | Baseline susceptibility of spotted wing Drosophila (Drosophila suzukii) to four key insecticide classes[END_REF][START_REF] Pavlova | Laboratory Bioassays with Three Different Substrates to Test the Efficacy of Insecticides against Various Stages of Drosophila suzukii (Diptera: Drosophilidae)[END_REF][START_REF] Van Timmeren | Development of a rapid assessment method for detecting insecticide resistance in spotted wing Drosophila (Drosophila suzukii Matsumura)[END_REF][START_REF] Valtierra-De-Luis | Quantification of dose-mortality responses in adult Diptera: Validation using Ceratitis capitata and Drosophila suzukii responses to spinosad[END_REF], without extensive investigations on the reliability of the methods used.

Because most insecticides targeting D. suzukii belong to the group of contact insecticides, we chose to use tarsal-contact bioassays. This type of bioassay is easy to set up and relatively inexpensive, facilitating its use worldwide. Organophosphates and pyrethroids have been the most effective and the most used insecticides in the field in France [42-44, 47, 49]; we therefore focused our experiments on two active substances of these families commonly used in Europe against D. suzukii: the organophosphate phosmet and the pyrethroid lambda-cyhalothrin. The formulated products act via contact, ingestion and inhalation for phosmet and via contact and ingestion for lambda-cyhalothrin.

The main aim of this study was to illustrate and to investigate how various biological and technical parameters can affect the accuracy and reliability of a pesticide resistance bioassay.

The key biological parameters we explored were (i) sex, (ii) age of the imago (adult stage) and (iii) genetic diversity of the tested population. Two methodological parameters were also tested: (iv) number of individuals used to evaluate the level of insecticide resistance and (v) duration of exposure to the active substance. One interesting and valuable output of this work is a reliable, accurate and standardized tarsal-contact bioassay to test insecticide resistance.

Our work highlights the essential factors to control in the development of reliable bioassays on any other pest species.

Materials and methods

Sample collection and preparation

The Ste-Foy population used in the experiments was originally established from approximately 20 D. suzukii females, collected from southeastern France in an urban area on raspberries (Sainte-Foy-les-Lyon, France) in May 2012 by Roland Allemand (Biometry and Evolutionary Biology Laboratory, CNRS-University of Lyon). The population was then mass-reared and maintained in standard drosophila vials (Ø 25 mm x h 95 mm) containing ~10 ml of a standard food media composed of 10 g of agar diluted in 1 l of water, 60 g of glucose, 30 g of saccharose, 80 g of malted yeast, 20 g of yeast extract, 20 g of peptones, 0.5 g of magnesium sulfate, 0.5 g of calcium chloride and 1 g of nipagin, pre-diluted in 10 ml of ethanol. They were maintained in a climatic chamber at 23±1˚C, under a relative humidity of approximately 70% and a 16:8 h light-dark cycle. Temperature was strictly controlled because it can influence D. suzukii susceptibility to insecticides [START_REF] Saeed | Temperature Alters the Response to Insecticides in Drosophila suzukii (Diptera: Drosophilidae)[END_REF]. Then, 25-30 males and females per vial were left to mate and oviposit for four to seven days before being removed from the rearing vials. After two weeks, the newly emerged adults (imagoes) were isolated and transferred to a new vial to start a new generation. To maintain the Ste-Foy population, each generation consisted of approximately 20 vials and the emerged adults of all the vials were mixed together to maximize genetic variability before distributing them into new vials.

A low genetic diversity population, SF-IsoA, was generated from the Ste-Foy population by inbreeding the flies for three generations: each generation, a single brother and virgin sister pair were isolated and left to reproduce in a new vial. The SF-IsoA population was then maintained in the same conditions as the Ste-Foy population.

Tarsal-contact bioassays

All the experiments were conducted using a common framework for the bioassays. Based on this framework, the five parameters of the protocol were modified to test for their effect on the assessed pesticide resistance. The tarsal-contact bioassay consists in exposing fly distal part of the legs (tarsi) to insecticide for a pre-determined time before assessing mortality. We used 20 ml scintillation glass vials (Ø 28 mm x h 61 mm), which have previously been used on D. suzukii with good results for resistance monitoring [START_REF] Smirle | Laboratory studies of insecticide efficacy and resistance in Drosophila suzukii (Matsumura) (Diptera: Drosophilidae) populations from British Columbia, Canada[END_REF][START_REF] Van Timmeren | Development of a rapid assessment method for detecting insecticide resistance in spotted wing Drosophila (Drosophila suzukii Matsumura)[END_REF]. A volume of 500 μl of the insecticide solution at the desired concentration (dissolved in acetone) was deposited on the walls of the vials. The insecticide solution was uniformly distributed with acetone evaporation by rolling the vials at room temperature for 1 h, on a hot-dog roller. Seven doses of the active substance were tested, including a control with acetone only. The range of concentrations differed for the two insecticides and depended on what concentrations allowed a dose-response curve ranging from 0 to 100% mortality. For the organophosphate phosmet (95% of purity), an acetylcholinesterase inhibitor and obtained from Gowan, the chosen range of concentrations was 0, 4.7, 9.4, 18.9, 37.7, 75.5 and 150.9 mg/l. Lambda-cyhalothrin (92.8% of purity) is a pyrethroid that targets the sodium channel and was furnished by Syngenta France SAS. The selected concentration range was adjusted according to the bioassay and consisted in 7 doses chosen from the following concentrations: 0, 0.01, 0.05, 0.1, 0.25, 0.5, 1, 3 and 5 mg/l. Because contact time lasted up to 24 h, 100 μl of an agarose-sucrose mixture (5% sucrose, 8 g/l agarose) was provided in each vial to prevent death by starvation. For a bioassay, one to four vials for the same dose were used, depending on the number of flies tested. The number of D. suzukii flies of a uniform age class (the standard being 24 to 48 hours old in this study) from a single population ranged from 3 to 32 by vial. Males and females were not separated prior to the bioassay to avoid anesthesia (necessary to determine the sex of the flies). The vials were plugged with thin netting held in place by perforated caps to prevent flies from escaping and still allow good ventilation. The vials were maintained in a climatic chamber during the bioassay at 20 ±1˚C, under a relative humidity of approximately 75% and a 16:8 h light-dark cycle. After a certain duration of exposure (the standard being 24 h), the vials were briefly shaken and the numbers of dead, moribund and live flies were assessed. Individuals that could not remain on their legs and those that showed unusual behavior (i.e. uncertain or irregular flight, twitching legs and/or uneven movements) were considered as moribund. Immobile adult flies were considered dead. The sex of the flies for the different categories was determined during this mortality assessment by checking for the presence of black spots on the wings which characterizes males. The susceptibility of a population to the chemical was assessed by calculating the LD 50 value (see Statistical analysis section below).

This protocol was used for conducting five experiments with some variation regarding the D. suzukii individuals or the conditions of the bioassay (see Table 1). For each experiment, several bioassays were performed and each bioassay was done on a different date. Within a bioassay, multiple vials per dose could be used in order to have a limited number of flies per vial (up to 32). For Experiment 1 (24 bioassays), the test aimed to compare the susceptibility to phosmet according to sex at 24 h of exposure on flies from Ste-foy population aged of 24 to 48 h. In Experiment 2 (26 bioassays), we explored the influence of age class on susceptibility of Ste-Foy population to 24h of exposition to phosmet by testing males and females of three different age classes: 0 to 24 h, 24 to 48 h and 72 to 96 h. Experiment 3 (5 bioassays) consisted in observing the impact of the genetic diversity of a population on its resistance to phosmet after 24h of exposure. To do so, the reference population Ste-Foy was compared to the inbred SF-IsoA population and flies were 24 to 48 h old. The respective levels of genetic diversity of the two populations were assessed using molecular tools (see below). Experiment 4 (24 bioassays) investigated the impact of the number of flies tested per bioassay (mean number of flies per dose per bioassay ranging from 6 to 37 for the females and from 9 to 35 for the males, see Table 2) on the accuracy of the LD 50 estimates after 24 h of exposure of 24 to 48 h old flies from Ste-Foy to phosmet. Finally, Experiment 5 (3 bioassays) explored the effect of the duration exposure to pyrethroid insecticide (lambda-cyhalothrin) on the evaluation of susceptibility of Ste-Foy population. Unlike organophosphate insecticides, pyrethroid can potentially induce a knockdown effect (due to several mutations in the sodium channel molecular target), which may cause an erroneous mortality assessment depending on the time of observation after initial exposure to the Assessing the neutral and adaptive genetic diversity of the populations DNA extraction. Adults (aged from 24 to 72 h old) were picked from the rearing stock in order to assess the genetic diversity of Ste-Foy and SF-IsoA populations. Flies were killed directly in 70% ethanol and the total DNA was extracted from the whole body within a week. DNA from 44 individual D. suzukii samples (30 females, 14 males) from each population (Ste-Foy and SF-IsoA) was extracted following a modified method based on Walsh et al. (1991) [START_REF] Walsh | Chelex 100 as a medium for simple extraction of DNA for PCRbased typing from forensic material[END_REF]. A volume of 100 μl of a 10% Chelex 100 solution (Bio-Rad) and 3% of 10 mg/ml proteinase K (Eurobio) was added to each sample. Then, each sample was crushed using 2 mm steel beads on a 1600 MiniG tissue homogenizer (Spex1 SamplePrep) at 1500 strokes/min for 5 sec, in a 96-well format PCR plate. Tissues were digested for 14 h at 56˚C with a Mastercycler thermocycler (Eppendorf) with a final temperature step of 30 min at 98˚C, the supernatant was used as DNA template for PCR reaction.

Microsatellite genotyping. The genetic diversity of the two studied populations was assessed using microsatellite markers. We used a slightly modified nested PCR approach described by Schuelke [START_REF] Schuelke | An economic method for the fluorescent labeling of PCR fragments[END_REF] coupled with a selection of 13 microsatellite markers developed by Fraimout et al. [START_REF] Fraimout | New set of microsatellite markers for the spotted-wing Drosophila suzukii (Diptera: Drosophilidae): A promising molecular tool for inferring the invasion history of this major insect pest[END_REF] (S1 Table ). A forward specific primer was conjugated with a 5'-GTTGTAAAACGACGGCCAGT-3' M13-tail at its 5' end, the same labeled universal M13-tail was used for fluorescence detection. We used a 1:10 ratio with 1 unit tailed forward specific primer to 10 units labeled M13 tails and 10 specific reverse primers. PCR amplifications were carried out on a Mastercycler thermocycler (Eppendorf) in a 12 μl reaction volume containing 1X GoTaq1 Flexi Buffer, 1.5 mM MgCl 2 , 0.1 mg/ml bovine serum albumin (BSA), 200 μM of each dNTP, 0.4 μM of the labelled M13 tail, 0.4 μM of the specific reverse primer, 0.04 μM of M13-tailed specific forward primer, 1 unit of GoTaq1 Flexi DNA Polymerase (Promega) and 2 μl of DNA template. The PCR conditions were: 3 min at 95˚C followed by 30 cycles at 95˚C for 30 sec, 57˚C for 45 sec, 72˚C for 45 sec and 10 cycles at 95˚C for 30 sec, 54˚C for 30 sec, 72˚C for 45 sec with a final extension step at 72˚C for 20 min.

Labeled PCR products were pool-plexed (up to 7 loci) by using 2-8 μl of each PCR in 50 μl of H 2 O, according to their 5' end-labeled dyes (dilution 1:8 for Tamra, 1:16 for Hex and Atto-565 and 1:32 for 6-Fam Dyes). A volume of 2 μl of the diluted PCR mixture, 7.8 μl of HiDi formamide, and 0.2 μl GeneScan™ 600 LIZ1 size standard (Applied Biosystems) were injected in an ABI 3730xl DNA Analyzer (Applied Biosystems) using POP7 polymer. These genotyping runs were analyzed using GeneMapper1 V4.1 Analysis Software (Applied Biosystems).

Checking for the kdr mutation. Knockdown resistance to pyrethroids is known to be induced by some mutations in its molecular target, the voltage-gated sodium channel. To assess the presence of the main kdr mutations that lead to the L1014F substitution [START_REF] Williamson | Identification of mutations in the houseflypara-type sodium channel gene associated with knockdown resistance (kdr) to pyrethroid insecticides[END_REF], we adapted a PCR-RFLP protocol from Franck et al. [START_REF] Franck | Genetic architecture in codling moth populations: comparison between microsatellite and insecticide resistance markers[END_REF]. Twenty adults (aged from 24 to 72 h old) were randomly picked from the rearing stock of Ste-Foy population as well as 20 flies from the same population that survived an exposition of 24 h at 0.25 mg/l of lambda-cyhalothrin (from Experiment 5) and their DNA was extracted as described above. A 371 bp PCR fragment was amplified on a Mastercycler thermocycler (Eppendorf) in a 12 μl reaction volume containing 1X GoTaq1 Flexi Buffer, 1.5 mM MgCl 2 , 0.1 mg/ml BSA, 200 μM of each dNTP, 0.4 μM of each primer (CKDR1 and Cgd2; see Table 3), 1 unit of GoTaq1 Flexi DNA Polymerase (Promega) and 2 μl of DNA template. The PCR conditions were: 3 min at 95˚C followed by 35 cycles at 95˚C for 30 sec, 56˚C for 45 sec, 72˚C for 45 sec with a final extension step at 72˚C for 20 min. The PCR products (5 μl) were digested at 37˚C for 16 h with 2 units of MluCI endonuclease and 1X of NEB Buffer (New England Biolabs) in 10 μl of reaction volume. The digested products (4 μl) were separated on a 2% agarose gel at 100 V for 30 min using the RunOneTM Electrophoresis System (Embi Technology). The size of the fragments was estimated by comparison with a 100 bp DNA ladder (Promega). The L1014F mutation reveals a restriction site recognized by the MluCI endonuclease, which cuts the 371 bp PCR fragment into two fragments of respectively 123 bp and 248 bp. Susceptible genotypes remain undigested.

Primers Ds-SKdr-F and Ds-Seq-R were designed specifically for this study. Primers CKDR1 and Cgd2, designed in our laboratory and used in previous studies [START_REF] Brun-Barale | Involvement of a sodium channel mutation in pyrethroid resistance in Cydia pomonella L, and development of a diagnostic test[END_REF][START_REF] Franck | Multiple origins of the sodium channel kdr mutations in codling moth populations[END_REF], were used as is given their perfect homologies.

Voltage-gated sodium channel gene sequencing. To verify the absence of the L1014F substitution and other secondary substitutions previously described in other insect species (i.e. M918T, L925I, T929I, L932F, C933A, I936V, G943A, Q945R, I1011M/V, N1013S and V1016G) [START_REF] Soderlund | The molecular biology of knockdown resistance to pyrethroid insecticides[END_REF][START_REF] Dong | Molecular biology of insect sodium channels and pyrethroid resistance[END_REF][START_REF] Hopkins | The Helicoverpa zea (Boddie) (Lepidoptera: Noctuidae) voltage-gated sodium channel and mutations associated with pyrethroid resistance in field-collected adult males[END_REF][START_REF] Usherwood | Mutations in DIIS5 and the DIIS4-S5 linker of Drosophila melanogaster sodium channel define binding domains for pyrethroids and DDT[END_REF], we relied on the complete genome of D. suzukii (GenBank assembly accession: GCA_000472105.1) focusing on scaffold023, positions 597.317 to 698.969, which correspond to the trans-membrane segments 4 to 6 of the domain II region of the voltage-gated sodium channel (S1 Fig) . A 1564 bp PCR fragment was amplified for two individuals from Ste-Foy population (that survived an exposition of 24 h at 0.25 mg/l of lambda-cyhalothrin) on a Mastercycler thermocycler (Eppendorf) in a 30 μl reaction volume, containing: 1X GoTaq1 Flexi Buffer, 1.5 mM MgCl 2 , 0.1 mg/ml BSA, 200 μM of each dNTP, 0.4 μM of each primer (Ds-SKdr-F and Cgd2) (see Table 1), 1 unit of GoTaq1 Flexi DNA Polymerase (Promega) and 2 μl of DNA template. The PCR conditions were: 3 min at 95˚C followed by 35 cycles at 95˚C for 30 sec, 54˚C for 45 sec, 72˚C for 2 min with a final extension step at 72˚C for 20 min. The PCR products (20 μl) were sequenced (Eurofins Genomics) using primers Ds-SKdr-F, Ds-Seq-R and CKDR1 (see Table 3). DNA sequences were manually aligned and analyzed using Bioedit software [START_REF] Hall | BioEdit: a user-friendly biological sequence alignments editor and analysis program for Windows 95/98/NT[END_REF].

Statistical analysis

Moribund and dead individuals were combined and considered as dead. Prior to data analyses, we kept only bioassays with a mortality rate in the control less than 15%. Survival data were analyzed in a non-linear regression framework using the 'drc' package [START_REF] Ritz | Dose-Response Analysis Using R[END_REF] in R [START_REF] Team | R: a language and environment for statistical computing[END_REF]. For all experiments, we assumed a binomial distribution of errors. For Experiments 1, 2 and 3, survival data were fitted to the three-parameter log-normal model, which is equivalent to the classic probit model [START_REF] Finney | Probit analysis[END_REF] with an additional parameter that takes into account the 'natural' mortality rate observed in the control of each categories. Due to the very low mortality rate for the different categories of Experiment 4 and Experiment 5 (<10%) as well as the sometimes limited number of individuals tested, survival data were fitted to the two-parameter log-normal model equivalent of the classic probit model [START_REF] Finney | Probit analysis[END_REF]. For Experiments 1, 2 3 and 5, the results of the bioassays were homogeneous and were therefore pooled together. The results of the bioassays included in Experiment 4 were not pooled because the purpose of this experiment was to assess the effect of the mean number of individuals in the bioassay on the LD 50 estimates. Models were fitted for the different categories separately within each experiment. The 'drc' package allows the estimation of LD 50 and the associated 95% CI and standard errors. Difference in LD 50 between males and females (Experiment 1) was tested directly using the 'compParm' function implemented in the 'drc' package. For the other experiments, the datasets were split into male and female subsets prior to the analysis. Pairwise comparisons of LD 50 between the different age classes, the populations with different levels of genetic diversity and the different times of exposure were performed for each sex separately for Experiments 2, 3 and 5, respectively, using the 'compParm' function. LD 50 estimated on each bioassay of Experiment 4 separately were compared with the LD 50 estimated on the bioassays pooled by sex based on the overlapping of the 95% CI and the respective values. Additionally, for Experiment 3, two genetic diversity indices were computed for the Ste-Foy and the SF-IsoA populations.

The number of alleles (N a ) and gene diversity (H e ) [START_REF] Nei | Molecular evolutionary genetics[END_REF] were estimated for both populations and for all loci with GENEPOP V4.2.2 [START_REF] Rousset | genepop'007: a complete re-implementation of the genepop software for Windows and Linux[END_REF]. The dataset and the code used for the different analyses, as well as for the production of the figures, have been deposited in an online repository (doi: 10.5281/zenodo.2842939).

Results

Effect of biological parameters on LD 50 estimates: Sex, individual fly age and population genetic diversity

The influence of sex on insecticide susceptibility. A total of 2745 females and 2341 males were tested in 24 bioassays. The LD 50 values for phosmet were estimated at 39.6 (95% CI: 37.2-42.0) and 19.8 (95% CI: 18.6-21.1) mg/l for females and males, respectively (Fig 1). The two-fold higher resistance to phosmet of females compared with males was highly significant (relative female:male potency = 1.99, t-value = 11.3, p-value < 0.001).

The influence of fly age on insecticide susceptibility. A significant effect was observed between the 0-24 h and 24-48 h age classes for males (relative potency 0-24 h:24-48 h = 0.78, t-value = -3.09, p-value = 0.002). No significant differences were observed between the other different age classes for male flies (relative potency 0- Unlike males, age seemed to involve a clear change in insecticide susceptibility for females (Fig 2). The two older age classes showed similar LD 50 values with 40.3 mg/l (95% CI: 36.9-43.8) for the 24-48 h females and 42.9 mg/l (95% CI: 37.3-48.5) for the 72-96 h age class (relative potency 24-48 h:72-96 h = 0.94, t-value = -0.80, p-value = 0.421). The youngest females were twice as susceptible to phosmet, with a LD 50 value of 18.6 mg/l (95% CI: 16.2-20.9), compared with the two other age classes (relative potency 0-24 h:24-48 h = 0.46, t-value = -14.9, p-value < 0.001 and relative potency 0-24 h:72-96 h = 0.43, tvalue = -14.0, p-value < 0.001).

The influence of genetic diversity on insecticide susceptibility. Microsatellite genotyping confirmed that the genetic diversity of the SF-IsoA population was reduced compared with the Ste-Foy population at 13 loci (Fig 3 andS2 Table). The mean number of alleles (per microsatellite marker) and the mean gene diversity was higher in Ste-Foy (N a = 3.31; H e = 0.595) than in SF-IsoA (N a = 1.46; H e = 0.179), confirming that Ste-Foy and SF-IsoA populations are appropriate biological materials to test the influence of genetic diversity on resistance.

The LD 50 values estimated by dose-response curve analyses were very similar with overlapping 95% confidence intervals for females (31.

Effect of the experimental setup on LD 50 values: Number of tested individuals and exposure duration

The influence of the number of tested individuals on insecticide susceptibility. We analyzed the variation in the estimated LD 50 values according to the number of tested individuals per bioassay per dose of phosmet (from a minimum mean number per dose of 9.4 and 10 flies for males and females, respectively, to a maximum mean number per dose of 34.9 and 37.1 for males and females, respectively). Our results clearly highlight the importance of the number of tested flies on the accuracy of the estimation of LD 50 (Fig 5). For both females and males, an insufficient number of tested individuals led to imprecise and haphazard evaluations of LD 50 (Fig 5). The standard errors associated with the LD 50 values estimated with a limited number of individuals tended to be higher when fewer individuals were included in the experiment. In contrast, when the mean number of individuals tested per dose exceeded a threshold of approximatively 30, the LD 50 assessed for the various bioassays were similar and close to the value estimated with every tested individual pooled together as one unique total bioassay according to sex. Males with 29, The influence of the duration of insecticide exposure on D. suzukii susceptibility. Exposure duration to lambda-cyhalothrin had a significant effect on the LD 50 estimates in the study population, both for females and males. The three bioassays were very similar: only 2 out of 30 LD50 comparisons between bioassays at each time of exposure were significantly different. For this reason, subsequent analyses were performed on pooled data from the bioassays. During the first phase of the experiment (from 1 to 5 h duration), we observed a decrease in the LD 50 estimate. In the second phase, the LD 50 estimate was stable (from 20 to 24 h duration) (Fig 6A). The duration of exposure before scoring had a significant effect on the LD 50 estimate (e.g. relative potency 5h:24h = 1.55, t-value = 2.53, p-value = 0.011) (Fig 6B). There was a clear significant break between 5 h and 20 h of exposure time. In particular, the estimation of LD 50 for 1, 2, 3 or 4 h exposure was less precise than after longer exposure times, with larger standard errors and 95% confidence intervals. The LD 50 estimates for 20 to 24 h of exposure were similar (relative potency close to 1) and not significantly different (Fig 6B ).

Detection of putative mutations in the voltage-gated sodium channel. The genotyping of 20 random individuals from Ste-Foy and 20 individuals from Ste-Foy that survived for 24 h at 0.25 mg/l of lambda-cyhalothrin did not reveal the presence of mutant L1014F kdr alleles. Therefore, the results observed in our lambda-cyhalothrin assays were not biased by this mutation. The partial sequencing of the gene encoding for the voltage-gated sodium channel of two flies from the Ste-Foy population that survived at 0.25 mg/l of lambda-cyhalothrin and exhibited a kdr-like phenotype did not reveal the presence of any of the 12 main resistance mutations described in the literature (GenBank accession nos. MK645039-MK645040).

Discussion

Our study explored a selection of critical parameters to control in pest bioassays. From an applied point of view, this study defined a useful reproducible protocol for a reliable tarsalcontact bioassay to evaluate insecticide resistance in D. suzukii.

We found that D. suzukii does not respond in the same way to insecticide according to sex and age. In most studies on insecticide resistance of D. suzukii, tests are done on females only or on a mix of males and females without taking sex into account. A few experiments have reported that males have greater insecticide susceptibility than females [START_REF] Bruck | Laboratory and field comparisons of insecticides to reduce infestation of Drosophila suzukii in berry crops[END_REF][START_REF] Smirle | Laboratory studies of insecticide efficacy and resistance in Drosophila suzukii (Matsumura) (Diptera: Drosophilidae) populations from British Columbia, Canada[END_REF][START_REF] Van Timmeren | Development of a rapid assessment method for detecting insecticide resistance in spotted wing Drosophila (Drosophila suzukii Matsumura)[END_REF]. This phenomenon has also been observed in various insect species and specifically in Diptera for various insecticides, including pyrethroids and organophosphates [START_REF] Carrière | Sex, and the Evolution of Pesticide Resistance[END_REF][START_REF] Koren | Biochemical Properties Characterizing the Development of Tolerance to Malathion in Ceratitis capitata Wiedemann (Diptera: Tephritidae)[END_REF][START_REF] Mcdonald | Linkage, Expression, and Distribution of the Pyrethroid Resistance Gene in the Horn Fly (Diptera: Muscidae)[END_REF][START_REF] Pruett | Evaluation of horn flies (Diptera: Muscidae) from a pyrethroid susceptible colony for general and permethrin esterase activities[END_REF]. We observed the same phenomenon in our experiment with males being twice as susceptible as females. The origin of this difference has already been discussed in other studies, being sometimes attributed to the difference in size between sexes due to sexual dimorphism [START_REF] Pluthero | Genetic differences in malathion avoidance and resistance in Drosophila melanogaster[END_REF][START_REF] Riedl | Monitoring susceptibility to azinphosmethyl in field populations of the codling moth (Lepidoptera: Tortricidae) with pheromone traps[END_REF][START_REF] Knight | Use of sex pheromone traps to monitor azinphosmethyl resistance in tufted apple bud moth (Lepidoptera: Tortricidae)[END_REF], although detoxification enzymes also seem to be involved in these differences [START_REF] Navarro-Roldan | Comparative Effect of Three Neurotoxic Insecticides With Different Modes of Action on Adult Males and Females of Three Tortricid Moth Pests[END_REF][START_REF] Shearer | Sex-related response to organophosphorus and carbamate insecticides in adult Oriental fruit moth, Grapholita molesta[END_REF]. Our results also showed an influence of age on susceptibility, but only for females. Susceptibility to pesticide can vary according to the development stage of the insect [START_REF] Reyes | Resistance monitoring in codling moth: a need for standardization[END_REF]. Smirle et al. [START_REF] Smirle | Laboratory studies of insecticide efficacy and resistance in Drosophila suzukii (Matsumura) (Diptera: Drosophilidae) populations from British Columbia, Canada[END_REF] reported that malathion (an organophosphate) induces more toxicity on 5-8-day-old D. suzukii adults than on 2-dayold individuals. Most experiments use flies at least 5 days old, which is time-consuming in terms of sample preparation. No experiment has indicated insecticide susceptibility in D. suzukii flies younger than 24 h. In our experiment, newly hatched females (0-24 h) showed significantly higher susceptibility to insecticides than older age classes. This period corresponds to the sexual maturation stage in this synovigenic species [START_REF] Revadi | Sexual Behavior of Drosophila suzukii[END_REF][START_REF] Kanzawa | Studies on Drosophila suzukii mats[END_REF]. During the ovarian maturation period, metabolic detoxification pathways and energy investments are different compared with the rest of the adult life [START_REF] Aguila | The role of larval fat cells in adult Drosophila melanogaster[END_REF][START_REF] Wheeler | The role of nutrition in oogenesis[END_REF]. Testing the effects of mating status on insecticide resistance might be an interesting avenue of research. Mating occurs mainly after the female refractory period of 1 to 2 days after emergence [START_REF] Revadi | Sexual Behavior of Drosophila suzukii[END_REF][START_REF] Greenspan | Courtship in drosophila[END_REF] and generates strong changes in female physiology and behavior in Drosophila [START_REF] Lee | Sexual dimorphism in nutrient intake and life span is mediated by mating in Drosophila melanogaster[END_REF][START_REF] Avila | Insect seminal fluid proteins: identification and function[END_REF]. On the contrary, young D. melanogaster flies have been described as more resistant to several abiotic stresses [START_REF] Colinet | Age-related Decline of Abiotic Stress Tolerance in Young Drosophila melanogaster Adults[END_REF]. Our finding is of interest for establishing a reliable bioassay method, but may also have an applied interest in targeting treatments on the early life stage of young females.

Inbred populations with low genetic diversity may display lower resistance to pathogenic organisms [START_REF] Spielman | Does inbreeding and loss of genetic diversity decrease disease resistance?[END_REF] and may be more sensitive to environmental stress [START_REF] Bulsma | Does inbreeding affect the extinction risk of small populations?: predictions from Drosophila[END_REF]. By creating a low genetic diversity population SF-IsoA, we were able to test the hypothesis that reduced genetic diversity may (i) reduce the variability of insecticide response estimates (i.e. LD 50 values) within the population and increase repeatability, and (ii) genetically fix the strain and obtain a stable reference population through generations. We therefore expected a different insecticide response in the SF-IsoA compared with the Ste-Foy population. However, no differences in susceptibility were detected in our comparison between original and inbred populations. Our finding is thus likely an illustration of the fact that the Ste-Foy population is initially devoid of resistance alleles and therefore allelic paucity does not change this status.

Regarding the effect of the experimental setup on insecticide susceptibility, the number of tested flies per bioassay and per dose and the duration of exposure to the tested insecticide are essential parameters to control. As suggested in previous studies, a range from 10 to 25 insects per dose is necessary [START_REF] Miller | Bioassays for monitoring insecticide resistance[END_REF][START_REF] Brogdon | Guideline for evaluating insecticide resistance in vectors using the CDC bottle bioassay[END_REF]. Our results showed that the number of tested insects is the most important factor among those explored in this work. According to our data, a reliable mean number of tested insects per dose is a minimum of 29 males and 34 females. We therefore recommend about 30 individuals per dose for the tests to be accurate and reproducible. For most of the existing protocols on D. suzukii, the duration of exposure to the chemical is usually about 24 h [START_REF] Bruck | Laboratory and field comparisons of insecticides to reduce infestation of Drosophila suzukii in berry crops[END_REF][START_REF] Smirle | Laboratory studies of insecticide efficacy and resistance in Drosophila suzukii (Matsumura) (Diptera: Drosophilidae) populations from British Columbia, Canada[END_REF][START_REF] Van Timmeren | Baseline susceptibility of spotted wing Drosophila (Drosophila suzukii) to four key insecticide classes[END_REF][START_REF] Pavlova | Laboratory Bioassays with Three Different Substrates to Test the Efficacy of Insecticides against Various Stages of Drosophila suzukii (Diptera: Drosophilidae)[END_REF] and can reach 48 or 72 h [START_REF] Hamby | Integrating circadian activity and gene expression profiles to predict chronotoxicity of Drosophila suzukii response to insecticides[END_REF][START_REF] Beers | Developing Drosophila suzukii management programs for sweet cherry in the western United States[END_REF][START_REF] Cuthbertson | Preliminary Screening of Potential Control Products against Drosophila suzukii[END_REF][START_REF] Shawer | Laboratory and field trials to identify effective chemical control strategies for integrated management of Drosophila suzukii in European cherry orchards[END_REF]. A study on the mortality induced by ingestion of spinosad on this species emphasized the importance of time elapsed (5 days) before assessing mortality, suggesting that exposure duration is crucial in different types of bioassays [START_REF] Valtierra-De-Luis | Quantification of dose-mortality responses in adult Diptera: Validation using Ceratitis capitata and Drosophila suzukii responses to spinosad[END_REF]. Bioassays can be extremely time consuming, making it tempting to minimize exposure duration to reduce the total testing time. However, for pyrethroids, several species have shown a knockdown phenotype that can induce a transient state of apparent mortality in the insects a few minutes after exposure to the insecticide. In a study of mosquitoes exposed to lambda-cyhalothrin, the time to 50% knockdown (KT 50 ) ranged from 14 (Mansonia africana) to 128 min (Anopheles gambiae) with KT 95 values of 42 and 277 min, respectively [START_REF] Agossa | Impact of insecticide resistance on the effectiveness of pyrethroid-based malaria vectors control tools in Benin: Decreased toxicity and repellent effect[END_REF]. In another study, the KT 50 for D. melanogaster flies ranged from 11 to 40 min, depending on the pyrethroid and the strain [START_REF] Chen | DSC1 channel-dependent developmental regulation of pyrethroid susceptibility in Drosophila melanogaster[END_REF]. Therefore, assessing the mortality too early after exposure may distort results, while a late, albeit reliable, assessment [START_REF] Shaw | Implications of sub-lethal rates of insecticides and daily time of application on Drosophila suzukii lifecycle[END_REF] may result in lost time. We observed a few individuals with a kdr-like phenotype in our experiment. They showed temporary paralysis followed by complete recovery of their locomotive capacities. Additional research did not reveal any of the main resistance mutations or the kdr mutation at the expected locus (1014) of the para gene. Similar results have been reported from a study on D. suzukii adults exposed to spinosad [START_REF] Van Timmeren | Development of a rapid assessment method for detecting insecticide resistance in spotted wing Drosophila (Drosophila suzukii Matsumura)[END_REF]-to which resistance is emerging in the US [START_REF] Gress | Identification and risk assessment of spinosad resistance in a California population of Drosophila suzukii[END_REF]. This study's results showed highly variable LD 50 values, calculated from the large number of moribund flies due to an insufficient exposure time (6 h). To avoid potential misinterpretation due to a knockdown effect, a rapid detoxification mechanism (potentially present in susceptible insects) or a slow mode of action of an active substance, we recommend assessing insect mortality at least 24 h after insecticide exposure (or more depending on the mode of action of the insecticide). Because the phenotypic effect of the insecticide is estimated by visual assessment, it is of utmost importance to establish precise scoring criteria that can be standardized across experiments and laboratory workers, in particular to differentiate moribund individuals from the dead and live individuals.

As illustrated by the low to moderate levels of resistance of D. suzukii populations that have been described in the US [START_REF] Gress | Identification and risk assessment of spinosad resistance in a California population of Drosophila suzukii[END_REF], monitoring insecticide resistance of this pest is important to be able to react quickly and adapt control management. We hereby propose a validated and operational method for a reliable worldwide monitoring of D. suzukii resistance to insecticides, and highlighted several parameters that are essential to control for in the design of resistance assessment bioassays on other species as well.

  24 h:72-96 h = 0.93, t-value = -0.64, pvalue = 0.523 and relative potency 24-48 h:72-96 h = 1.19, t-value = 1.51, p-value = 0.132) (Fig 2). The estimates of the LD 50 values were very similar for the three age classes with 15.2 mg/l (95% CI: 12.9-17.5), 19.5 mg/l (95% CI: 17.7-21.2) and 16.4 mg/l (95% CI: 13.4-19.4), from youngest to oldest.

Fig 1 .

 1 Fig 1. Effect of sex on LD 50 in a D. suzukii population. Dose-response curves of 24 to 48 h old male (light gray triangles, shading and dashed line) and female (dark gray circles, shading and solid line) adults of D. suzukii after 24 h of tarsal exposure to phosmet. The 95% confidence intervals were derived from the dose-response model. https://doi.org/10.1371/journal.pone.0247756.g001

  9 mg/l (95% CI: 15.2-48.6) and 41.8 mg/l (95% CI: 35.7-48.0) for SF-IsoA and Ste-Foy, respectively) and males (22.4 mg/l (95% CI: 21.3-23.6) and 17.7 mg/l (95% CI: 13.2-22.3) for SF-IsoA and Ste-Foy, respectively) (Fig 4). Intra-

Fig 2 .

 2 Fig 2. Effect of age class on LD 50 in a D. suzukii population. Results obtained after 24 h tarsal exposure to phosmet for females (upper panel) and males (lower panel). Red, green and blue represent the 0-24 h, 24-48 h and 72-96 h age classes, respectively. Dose-response curves and 95% confidence intervals were derived from the dose-response model. https://doi.org/10.1371/journal.pone.0247756.g002

Fig 3 .Fig 4 .

 34 Fig 3. Diversity indices of the Ste Foy and SF-IsoA populations. Radar plot of (A) gene diversity (H e ) and (B) the number of observed alleles (N a ) at 13 microsatellite loci for the Ste-Foy (red) and SF-IsoA populations (blue). https://doi.org/10.1371/journal.pone.0247756.g003

  30 and 35 individuals per dose led to a similar estimation of LD 50 values with overlapping 95% confidence intervals: 21.1 mg/l (95% CI: 17.7-24.6), 22.4 mg/l (95% CI: 19.3-25.5) and 21.8 mg/l (95% CI: 18.1-25.4), respectively. Similarly, females gave LD 50 estimates of 38.3 mg/l (95% CI: 30.8-45.8), 46.1 mg/l (95% CI: 38.0-54.1) and 36.0 mg/l (95% CI: 30.2-41.8), for respectively 34, 36 and 37 individuals per dose.

Fig 5 .

 5 Fig 5. Effect of the mean number of individuals tested per dose on LD 50 in a D. suzukii population. Distribution of LD 50 values and associated 95% confidence intervals after 24 h of tarsal exposure of 24 to 48 hour-old adults to phosmet, according to the number of individuals tested for females (left panel) and males (right panel). The continuous red line indicates the LD 50 value computed on pooled individuals as one unique bioassay and the dashed lines show the corresponding 95% confidence interval. https://doi.org/10.1371/journal.pone.0247756.g005

Fig 6 .

 6 Fig 6. Effect of duration of exposure on LD 50 estimates in a D. suzukii population. The bioassay consisted of a tarsal exposure to lambda-cyhalothrin and was conducted on 24 to 48 hour-old D. suzukii adults from the Ste-Foy population. Panel A shows the change in the dose-response curves for female individuals. Early scoring (after 1, 2, 3, 4 and 5 h) curves are shown in black and late scoring (after 20, 21, 22, 23 and 24 h) curves are shown in red. The black arrow shows the variation in LD 50 between 1 and 5 h. The half-matrix in Panel B shows the significance levels of relative potencies between different exposure times for the female individuals (red: p< 0.001; orange: 0.001 �p <0.01; yellow: 0.01 � p <0.05; gray: ns). https://doi.org/10.1371/journal.pone.0247756.g006

Table 1 . Parameters tested in the five experiments. Colonne1 Tested parameters Insecticide Population Age of the flies Exposure duration to insecticide Number of bioassays Experiment 1

 1

		Sex	Phosmet	Ste Foy	24-48h	24h	24
	Experiment	Age	Phosmet	Ste Foy	0-24h; 24-48h; 72-	24h	26
	2				96h
	Experiment	Genetic diversity	Phosmet	Ste Foy;	24-48h	24h	5
	3			SF-IsoA
	Experiment	Insect number per dose	Phosmet	Ste Foy	24-48h	24h	24
	4
	Experiment	Duration of insecticide	Lambda-	Ste Foy	24-48h
	5	exposure	cyhalothrin

1; 2; 3; 4; 5; 20; 21; 22; 23; 24h 3 https://doi.org/10.1371/journal.pone.0247756.t001 chemical. That is why we conducted the bioassays with lambda-cyhalothrin in this Experiment 5. Mortality was assessed repeatedly on the same vials at 10 different times: 1, 2, 3, 4, 5, 20, 21, 22, 23 and 24 h after initial exposure.

Table 2 .

 2 Characteristics of the study's bioassays. For each of the bioassays, the date of performance, the inclusion in the different experiments, the mean number of flies per dose (and per sex) and the total number of individuals used are indicated.

	Mean

number of flies per dose Date of the bioassay Experiment 1 Experiment 2 Experiment 3 Experiment 4 Experiment 5 ♀ ♂ ♀♂ Total number of flies per bioassay

	25/05/2016	X	X	X	20 14	34
	01/06/2016	X	X	X	23 20	43
	08/06/2016	X	X	X	15 16	31
	13/06/2016		X		9	7	16
	14/06/2016		X		22 21	42
	16/06/2016		X		21 15	36
	20/06/2016		X		9	7	16
	22/06/2016	X	X	X	22 21	44
	23/06/2016	X	X	X	32 29	61
	27/06/2016		X		11	7	18
	28/06/2016		X		13 13	26
	30/06/2016	X	X	X	18 15	33
	04/07/2016		X		13 12	25
	05/07/2016		X		12	9	22
	07/07/2016	X	X	X	18 14	32
	11/07/2016		X		11	8	19
	12/07/2016		X		6	4	10
	21/07/2016	X	X	X	13	9	22
	25/07/2016		X		13	8	22
	01/09/2016	X	X	X	13 10	23
	08/09/2016	X	X	X	13 14	27
	15/09/2016	X	X	X	20 17	36
	22/09/2016	X	X	X	23 18	41
	13/10/2016	X	X	X	10 10	20
	20/10/2016	X	X	X	10 10	20
	27/10/2016	X	X	X	12 11	24
	01/03/2017	X	X	X	17 17	34
	15/03/2017	X	X	X	18 15	33
	23/03/2017	X	X	X	14 15	29
	29/03/2017	X		X	18 12	30
	05/04/2017	X	X	X	21 16	37
	06/04/2017	X	X	X	24 22	47
	04/05/2017	X		X	37 29	66
	11/05/2017	X		X	36 30	66
	17/05/2017	X		X	34 35	69
	09/09/2020			X	37 24	61
	30/09/2020			X	23 13	35
	01/10/2020			X	38 28	66
	♀-Female flies
	♂-male flies of Ste-Foy or SF-IsoA population of Drosophila suzukii.
	https://doi.org/10.1371/journal.pone.0247756.t002

Table 3 . Primers used to amplify and sequence the trans-membrane segments 4 to 6 of the domain II region of the voltage-gated sodium channel of Drosophila suzukii.

 3

	PRIMER NAME	SEQUENCE (5'! 3')
	Ds-SKdr-F	TGGCCAACACTTAATTTACTC
	Ds-seq-R	CAAGAAGAAGGGAATGCAC
	CKDR1	CACAGCTTCATGATCGTGTTC
	Cgd2	GCAAGGCTAAGAAAAGGTTAAG
	https://doi.org/10.1371/journal.pone.0247756.t003

PLOS ONE | https://doi.org/10.1371/journal.pone.0247756March 5, 2021

Acknowledgments

We are grateful to the Biometry and Evolutionary Biology Laboratory (LBBE, University of Lyon) for providing the Ste-Foy population. We also thank Christine Brazier, Elorri Segura and Ce ´lia Goutebroze for their assistance with the lab work.

BB received funding from ANSES via the tax on the sales of plant protection products. As stipulated in the French Act on the Future of Agriculture of 13 October 2014, the proceeds of this tax are allocated to ANSES to finance research on establishing a system for monitoring the adverse effects of plant protection products ('phytopharmacovigilance', PPV). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Supporting information

S1 Table. Thirteen microsatellites markers used in