Vincent Monchiet
email: vincent.monchiet@ed.univ-lille1.fr

Cosmin Gruescu
email: ic.gruescu@univ-lille1.fr

Eric Charkaluk
email: eric.charkaluk@univ-lille1.fr

Djimedo Kondo
email: kondo@univ-lille1.fr

Approximate yield criteria for anisotropic metals with prolate or oblate voids

Keywords: Computational solid mechanics, Ductile porous metals, Anisotropy, Hill criterion, Prolate and oblate voids Mécanique des solides numérique, Métaux ductiles poreux, Anisotropie, Critère de Hill, Cavités allongées et aplaties

 which has extended the well-known approach of Gurson (1975), we propose approximate yield criteria for anisotropic plastic voided metals containing non spherical cavities. The plastic anisotropy of the matrix is described by means of Hill's quadratic criterion. The procedure to establish the closed form expression of approximate macroscopic criteria, in which void shape and plastic anisotropic effects are included, is detailed. The new criteria allow us to recover existing results in the cases of spherical and cylindrical voids in an Hill type plastic matrix. Moreover, they agree with previous criteria for non spherical voids in an isotropic plastic matrix. Finally, for validation purposes, we provide, in the general case of non spherical cavities in the anisotropic matrix, a comparison with the numerical exact two field criteria.

Résumé

 qui étendent l'approche bien connue de Gurson, on propose des critères macroscopiques pour des métaux plastiques anisotropes contenant des cavités non sphériques. L'anisotropie plastique de la matrice est décrite à l'aide du critère quadratique de Hill. On détaille la procédure pour établir les critères approchés incluant la forme des cavités ainsi que les effets d'anisotropie plastique. Les nouveaux critères analytiques obtenus permettent de retrouver les résultats existants dans les cas de cavités sphériques ou cylindriques dans une matrice de type Hill. De plus, ils concordent avec ceux qui sont disponibles pour des cavités non sphériques dans un milieu plastique isotrope. Enfin, dans un but de validation, on fournit, dans le cas général de cavités non sphériques dans la matrice anisotrope, une comparaison avec les expressions des critères issues des solutions numériques exactes du problème à deux champs.

Introduction

Void nucleation, growth and coalescence are commonly recognized as the basic micromechanisms of the ductile fracture of metals. Since the pioneering work of Gurson [START_REF] Gurson | Continuum theory of ductile rupture by void nucleation and growth: Part I-Yield criterion and flow rules for porous ductile media[END_REF] for a hollow sphere in a von Mises perfect plastic matrix, the modelling of voids growth has been the subject of several works performed in the context of ductile damage mechanics. Gologanu et al. [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric prolate ellipsoidal cavities[END_REF][START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric oblate ellipsoidal cavities[END_REF] (see also [START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF]), Garajeu et al. [START_REF] Garajeu | A micromechanical approach of damage in viscoplastic materials by evolution in size, shape and distribution of voids[END_REF] provide various extensions of the Gurson model by incorporating void shape effects; either prolate and oblate voids are considered in these studies. Other recent extensions of the Gurson model concern the consideration of plastic anisotropy. For instance, [START_REF] Benzerga | Plastic potentials for anisotropic porous solids[END_REF][START_REF] Liao | Approximate yield criteria for anisotropic porous ductile sheet metals[END_REF][START_REF] Wang | An anisotropic Gurson yield criterion for porous ductile sheet metals with planar anisotropy[END_REF] deal with void growth in a metal matrix, obeying the Hill quadratic criterion [START_REF] Hill | A theory of yielding and plastic flow of anisotropic solids[END_REF]. Unfortunately, all these studies are limited to the case of spherical and cylindrical voids. The purpose of the present Note is to provide a closed form approximated expression of the macroscopic yield function of anisotropic metals containing non spherical voids. Practical applications of the new criteria can be found in the domain of metal forming (see [START_REF] Son | Prediction of forming limits for anisotropic sheets containing prolate ellipsoidal voids[END_REF]).

Principle of the estimate of the macroscopic yield criterion

Statement of the problem

Following Gologanu et al. [START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF], consider a spheroidal (axisymmetric) prolate or oblate cavity with semi-axes a 1 (along e 3 ), and b 1 (along e 1 and e 2 ) embedded in a cell which has the shape of a confocal spheroid with the semiaxes a 2 (along e 3 ), and b 2 (along e 1 and e 2 ). a 1 >b 1 corresponds to a prolate cavity while b 1 >a 1 is associated to an oblate one. Let us denote c the focal distance and e 1 the eccentricity defined by: c = a 2 1b 2 1 , e 1 = c/a 1 for a prolate and c = b 2 1a 2 1 , e 1 = c/b 1 for an oblate cavity (see Fig. 1). For this class of geometry, it is convenient to introduce the system of spheroidal coordinates characterized by λ, β, ϕ, defined, in the cylindrical frame (coordinates ρ,ϕ,z), by: ρ = c sinh λ sin β and z = c cosh λ cos β for a prolate; in the case of an oblate, one has ρ = c cosh λ sin β and z = c sinh λ cos β.

The ). The velocity field, v, in the matrix is classically decomposed in a uniform strain rate, A.x, and an inhomogeneous field v E : v = A.x + Bv E . We adopt the trial velocity field used by [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric prolate ellipsoidal cavities[END_REF] and [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric oblate ellipsoidal cavities[END_REF] for nonspherical cavities analysis in a von-Mises matrix. Even though the accuracy expected from this choice is questionable in the anisotropic case, 1 the estimate of the macroscopic plastic dissipation being greater than the exact one, it is expected that the approximate yield surface will be external to the exact one and will probably need to be improved in further work.

By considering the uniform strain rate condition on the cell boundary, v(λ = λ 2 ) = D.x (D being the macroscopic strain rate), A and B can be identified. Indeed, denoting D h 1 the hydrostatic part of D:

A = D -D h X; B = a 2 b 2 2 c 3 D h with: X = 3 2 (1 -α 2 )( e 1 ⊗ e 1 + e 2 ⊗ e 2 ) + 3α 2 e 3 ⊗ e 3 (2)
α 2 = α(e 2 ). α(e), the function of eccentricity e, is defined by: The strain rate field is defined by d = A + Bd E ,forwhichd E is obtained from the inhomogeneous velocity field v E in the spheroidal coordinates. Let us recall that we adopt the velocity fields proposed by [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric prolate ellipsoidal cavities[END_REF] and [START_REF] Gologanu | Approximate models for ductile metals containing non-spherical voids-case of axisymmetric oblate ellipsoidal cavities[END_REF], instead of those provided in [START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF]. It follows that:

α(e) = ⎧ ⎨ ⎩ 1-e 2
d E = 3c 3 (1 -α) 2ab 2 (1 -3e λ ⊗ e λ ) + 3ac 3 (1 -3α) 2b 2 L 2 λ sin 2 (β)( e λ ⊗ e λ -e β ⊗ e β ) + 3c 3 (1 -3α) 2bL 2 λ sin(2β)e λ s ⊗ e β (4)

Principle of the determination of the macroscopic criteria

We consider now that the plastic matrix obeys the Hill quadratic criterion which can be expressed:

F(σ 11 -σ 22 ) 2 + G(σ 22 -σ 33 ) 2 + H(σ 33 -σ 11 ) 2 + 2Lσ 2 23 + 2Mσ 2 13 + 2Nσ 2 12 - 2 3 σ 2 0 0( 5 )
in which F,G,...,N are material constants and σ 0 the yield stress of the matrix. σ denotes the microscopic stress field. This expression, written in the frame of material anisotropy, can be translated as:

F(σ ) = 3 2 σ : M : σ -σ 2 0 0 ( 6 )
M is a general 2 fourth order anisotropic tensor whose components are expressed in the frame of the spheroid ( e 1 ,e 2 ,e 3 ). When the voids orientation and the direction of plastic anisotropy coincide, M is orthotropic in nature. In any case, the plastic incompressibility condition in the matrix reads M iikl = 0. Let us recall now the definition of the dissipation:

Π(D) = σ 0 |Ω| Ω d eq dV (7)
in which Ω denotes the considered domain whose volume is |Ω|=4πa 2 b 2 2 /3; d eq is the microscopic equivalent plastic strain rate which is written:

d 2 eq = 2 3 d : H : d = A 2 eq + 4 3 A : H : d E + d E eq 2 (8)
1 Indeed, the isotropic expansion field is no longer the exact field under hydrostatic loading. However, note that, by adopting this trial velocity fields, we follow the same kind of consideration as in [START_REF] Benzerga | Plastic potentials for anisotropic porous solids[END_REF] for spherical voids and as in [START_REF] Liao | Approximate yield criteria for anisotropic porous ductile sheet metals[END_REF] for cylindrical ones in the context of plastic anisotropy. 2 Such generality allows us to consider arbitrarily orientated voids with respect to the plastic anisotropy direction. Note also that the (isotropic)

von Mises criterion of the matrix is obtained from ( 6) by putting M = K = I -J where I is the fourth order symmetric identity tensor and

J = 1 3 1 ⊗ 1.
H is the fourth order tensor defined by H : M = K; it also satisfies the matrix incompressibility condition, H iikl = 0.

Expressions of H ij kl as function of the components of M are given in Appendix A. The yield surface related to the macroscopic dissipation is classically deduced from the two following relations:

3Σ h = ∂Π ∂D h = ∂Π ∂A : (1 -X) + ∂Π ∂B a 2 b 2 2 c 3 and Σ = ∂Π ∂D = ∂Π ∂A : ∂A ∂D + ∂Π ∂B ∂B ∂D = ∂Π ∂A (9)
which, by combination, leads to:

3Σ * h = Σ : X = ∂Π ∂B a 2 b 2 2 c 3 ; and Σ = ∂Π ∂A (10)
D and Σ are the deviatoric part of the macroscopic strain rate and the macroscopic stress respectively. Due to the difficulty of integrating the microscopic dissipation on Ω (Eq. ( 7)), the purpose of the next section is to provide some approximate expressions of the macroscopic dissipation and then of the yield criteria, either for prolate and oblate voids in anisotropic plastic matrix.

Macroscopic yield criteria for an anisotropic porous material

Approximate expression of the macroscopic dissipation

In order to establish an approximate expression of the macroscopic criteria by using [START_REF] Son | Prediction of forming limits for anisotropic sheets containing prolate ellipsoidal voids[END_REF], an analytical expression of the Π(D) is due. In the spheroidal frame, one has:

Π(D) = 1 |Ω| Ω d eq dV = 3c 3 4πa 2 b 2 2 λ=λ 2 λ=λ 1 β=π β=0 ϕ=2π ϕ=0 d eq bL 2 λ sin β dλ dβ dϕ (11)
The first step of the approximation consists of replacing d eq by a mean value along each confocal spheroid to the cavity. It follows that (11) can be put in the form:

Π(D) = c 3 a 2 b 2 2 λ=λ 2 λ=λ 1 d 2 eq E 1/2 b 2a 2 + b 2 dλ ( 12
)
The proposal made for (12) consists in replacing d 2 eq by the following mean value on each confocal spheroid:

d 2 eq E = 3 4π(2a 2 + b 2 ) β=π β=0 ϕ=2π ϕ=0 d 2 eq L 2 λ sin(β) dβ dϕ (13)
This is motivated by the fact that, contrarily to the isotropic case, d eq is function of ϕ in the context of plastic anisotropy. Using (8), d 2 eq E reads:

d 2 eq E = A 2 eq + 4 3 A : H : d E E + d E eq 2 E ( 14
)
with

d E E = 3 4π(2a 2 + b 2 ) β=π β=0 ϕ=2π ϕ=0 d E L 2 λ sin(β) dβ dϕ = 3Q(e) 2 Q ( 15
)
where Q is a second order tensor and Q(e) is a function of eccentricity e. With α = α(e) given by ( 3):

Q =- 1 2 (e 1 ⊗ e 1 + e 2 ⊗ e 2 ) + e 3 ⊗ e 3 ; Q(e) = ⎧ ⎪ ⎨ ⎪ ⎩ 2 3α-1+e 2 (1-α) (1-e 2 )(3-e 2 ) e 3 (prolate) 2 3α-1-2e 2 α √ 1-e 2 (3-2e 2 ) e 3 (oblate) (16)
Themeanvalueofd E eq 2 on each confocal spheroid reads:

d E eq 2 E = P(e) (17)
For prolate cavities, one has:

P(e)= 3 2 e 4 (1 -e 2 ) 2 (3 -e 2 ) 3e 2 α(1 -α) -1 + 3α h 1 + 9 1 -e 2 + αe 2 e 2 (1 -α) -1 + 3α h 2 + 4(1 -3α) 1 -e 2 h 3 (18)
and for oblate cavities:

P(e)= 3 2 e 4 3 -2e 2 1 -3α + 2e 2 h 1 + 9 1 -αe 2 1 -e 2 1 -3α + 2αe 2 h 2 -4(1 -3α)h 3 ( 19
)
where parameters h i are related to the anisotropy coefficients (using Voigt's notations) by:

h 1 = H 11 + H 22 + H 66 - 1 2 H 33 ,h 2 = H 33 ,h 3 = H 44 + H 55 (20)
To summarize, d 2 eq E , given by ( 14), takes the form:

d 2 eq E = A 2 eq + 2Q(e)BA : H : Q + P(e)B 2 (21)
The determination of the macroscopic dissipation requires now the integration of ( 12) with ( 21) over the vari-

able λ ∈[ λ 1 ,λ 2 ].
Once again, we need to approximate the result. Let us introduce the following notations:

µ s = 2h 1 +3h 2 +2h 3 10
; µ c = 3h 1 8 and µ p = 3h 3 + 27π 2 32 h 2 . In the case of prolate cavities, it is convenient to apply the following change of variable x = c 3 /(ab 2 ). It is then easy to check that P(e)≃ 4µ 2 s x 2 for x → 0 and P(e)≃ 4µ 2 c x 2 for x →+∞. It appears that P(e) can be considered approximately 'proportional' to x 2 . For oblate cavities, it is easy to check that P(e)≃ 4µ 2 s x 2 for x → 0 and P(e)→ 4µ 2 p for x →+∞. The following change of variable is then introduced: y = µ p x x+µ p . P(e) is then approximately 'proportional' to y 2 .

As a unified notation, let us introduce u such that u = x for prolate cavities and u = y for oblate cavities. The most important approximation is made now, in order to replace (12) by an integral of the type √ F 2 + G 2 u 2 du/u 2 . Following [START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF], P(e) and Q(e) are respectively replaced by p 2 (e)u 2 and q(e)u 2 . Since p(e) exhibits low variation according to e, it is replaced by a mean value which we have to be evaluate later. According to [START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF], q(e) is also replaced by a constant. Although this approximation is not completely justified, q(e) →+∞when e → 0, it was observed in the isotropic case that it gives particularly good results for the approximate criteria. Using the above approximations, d 2 eq becomes:

d 2 eq = A 2 eq - q 2 u 2 p 2 (A : H : Q) 2 + pB + q p A : H : Q 2 u 2 (22)
Since the effect of the term q 2 u 2 p 2 is low, it can also be replaced by a constant, denoted by convenience 2r 3p 2 .T h e macroscopic dissipation is then expressed as:

Π(D) = σ 0 x 2 u 2 u 1 Ã2 + B 2 u 2 1/2 du u 2 = σ 0 x 2 B arcsinh u B Ã - Ã2 + u 2 B 2 u u 2 u 1 ( 23
)
for which is used the following change of variable: Ã2 = A 2 eq -2r 3p 2 (A : H : Q) 2 and B = pB + q p (A : H : Q).Itis recalled that u = x for a prolate cavity and u = y for an oblate one.

Determination of the yield function

We have to determine now the macroscopic yield surface by reporting (23) in [START_REF] Son | Prediction of forming limits for anisotropic sheets containing prolate ellipsoidal voids[END_REF], using the following change of variables (A,B ⇒ Ã, B):

3Σ * h = p x 2 ∂Π ∂ B ; Σ - 3qx 2 p 2 Σ * h H : Q = ∂Π ∂ Ã ∂ Ã ∂A ( 24
)
In order to ease the expression of the yield surface, it is useful to put Ã2 in the form Ã2 = 2 3 A :

H * : A, with H * = H -r p 2 (H : Q) ⊗ (H : Q).
Let us denote by Σ * eq , an equivalent stress defined by:

(Σ * eq ) 2 = 3 2 Σ - 3qx 2 p 2 Σ * h H : Q : M * : Σ - 3qx 2 p 2 Σ * h H : Q = ∂Π ∂ Ã 2 with: M * = M + r p 2 -rQ : H : Q Q ⊗ Q (25)
such that M * : H * = K. For prolate cavities, u 2 = x 2 and u 1 = x 1 = x 2 /f . For oblate cavities, let us introduce g defined by g = x 2 /µ p ; this implies x 1 = µ p g/f , y 2 = µ p g/(g + 1) and y 1 = µ p g/(g + f).

The yield surface, deduced from (24) for prolate and oblate voids, takes the general form:

F macro = Σ * eq σ 0 2 + 2(1 + g)(f + g)cosh 3Σ * h pσ 0 -(1 + g) 2 -(f + g) 2 (26)
with g conventionally taken as zero (g = 0) in the case of prolate voids. In order to make explicit (26) which is the most important result of the study, a more detailed expression of Σ * eq is required, using (25):

Σ * 2 eq = Σ 2 eq + k 1 (Σ * h ) 2 + k 2 Σ : Q 2 -2k 3 Σ * h Σ : Q; Σ * h = 1 3 Σ : X (27)
where X is given by (2). k 1 , k 2 and k 3 depend on p, q, r and Q : H : Q = 9h 2 /4as:

k 1 = 243 2p 2 κ 2 h 2 4p 2 -9rh 2 ; k 2 = 6r 4p 2 -9rh 2 ; k 3 = 18κ 4p 2 -9rh 2 ; (28)
with κ = qx 2 . The macroscopic yield criterion is completely defined by ( 26), ( 27) and (28) in which p, r and κ have to be determined; this is the subject of the following two subsections.

Calculation of the coefficient p

We now come to the determination of the formulae which will give the coefficient p. An efficient method used by [START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF] in the isotropic case consists in considering the exact value in the particular case where A eq = 0. For the plastic anisotropic matrix case, p(e) is also function of h 1 , h 2 and h 3 and the above procedure cannot be used. We propose therefore to identify p as follows:

p 2 = 1 u 2 -u 1 u 2 u 1 p 2 (e) du (29)
With the notation α 1 = α(e 1 ) and α 2 = α(e 2 ), this leads, for prolate voids (using u = x), to:

p 2 = 3 4 1 1 -f 1 -3α 1 e 2 1 -f 1 -3α 2 e 2 2 (h 1 + 9h 2 -4h 3 ) + 18 α 1 (1 -α 1 ) -fα 2 (1 -α 2 ) h 2 + 1 -α 1 -f(1 -α 2 ) (h 1 -9h 2 + 4h 3 ) (30)
and for oblate voids (using u = y):

p 2 = 3 4 (1 + g)(f + g) f(1 -f) f 1 -3α 2 e 2 2 - 1 -3α 1 e 2 1 (h 1 + 9h 2 -4h 3 ) + 2(1 -f)h 1 + 18 fα 2 2 -α 2 1 h 2 + 4(f α 2 -α 1 )(h 1 -2h 3 ) (31)

Calculation of the coefficient κ and r

The procedure used consists of identifying κ and r with the two fields exact solution in the particular case B → 0. The identification method, given in Appendix B, leads to:3

κ = 2(1 + g)(f + g)(α 2 -α 1 ) (1 -f) ; r = 6(1 + g)(f + g)(α 2 -α 1 ) 2 (1 -f) 2 (32)
for which, we recall that g = 0 for prolate voids. It is interesting to notice that in the expression of κ = qx 2 , q could be also deduced from the same procedure used for the determination of p.

Validation

A first step of validation of the results is performed by considering special cases. Indeed, for the hollow sphere which is the particular case of a prolate or an oblate corresponding to e 1 = e 2 = 0, this implies that κ = r = g = 0 and then Σ * eq = Σ eq and Σ * h = Σ h and p 2 = 2 5 (2h ).This is exactly the result already established by [START_REF] Benzerga | Plastic potentials for anisotropic porous solids[END_REF].

The case of cylindrical voids corresponds to the following limiting case of prolate voids: e 1 = e 2 = 1; then κ = r = 0, Σ * eq = Σ eq , Σ * h = (Σ xx + Σ yy )/2 and p 2 = 3h 1 2 = 3 2 (H 11 + H 22 -1 2 H 33 + H 66 ). We also mention that for an isotropic matrix (h 1 = h 3 = 2 and h 2 = 2/3), p = 2 for spherical voids and p = √ 3for cylindrical voids which is in agreement with Gurson's results. We now propose a comparison between the approximate

Appendix B. Determination of parameters κ and r

Let us consider the determination of κ and r. The exact point corresponding to B → 0 is defined by:

Σ * h σ 0 = 1 3x 2 lim B→0 ∂Π ∂D = 1 x 2 1 |Ω| Ω 2 3A eq A : H : d E dV = 2 3 A : H : Q A eq (α 2 -α 1 ) Σ σ 0 = lim B→0 ∂Π ∂D = 2 3 A : H A eq (1 -f) (B.1)
Combining these two relations, it can be shown that:

Σ * h = α 2 -α 1 1 -f Σ : Q; Σ eq = σ 0 (1 -f) (B.2)
Applying B → 0 in the approximate solution (24) yields: Since Σ * h is assumed to be low, we propose to develop the hyperbolic terms in (B.4) at the second order. The identification of κ and r is made by replacing in (B.4), Σ * h and Σ eq by their exact expressions given by (B.2). The two equations (B.4) then become:

⎧ ⎪ ⎨ ⎪ ⎩ 3Σ * h = p
⎧ ⎪ ⎨ ⎪ ⎩ 1-f α 2 -α 1 - 2x 2 2 u 2 u 1 κ -9h 2 2κp 2 3κ 2 2 - rx 2 2 u 2 u 1 = 0 r(1-f) 2 6(α 2 -α 1 ) 2 -κ 1-f α 2 -α 1 + x 2 2 u 2 u 1 + 9h 2
κ = 2x 2 2 (α 2 -α 1 ) u 2 u 1 (1 -f) ; r = 6x 2 2 (α 2 -α 1 ) 2 u 2 u 1 (1 -f) 2 (B.6)

  iso-λ surface defines confocal spheroids, with semi axes a = c cosh(λ), b = c sinh(λ) and eccentricity e = c/a, for a prolate. An oblate spheroid is associated to semi axes a = c sinh(λ), b = c cosh(λ) and eccentricity c/b. The unit vectors of the new base are: e λ = 1 L λ a sin(β)e ρ + b cos(β)e 3 ; e β = 1 L λ b cos(β)e ρa sin(β)e 3 ; e ϕ = e ϕ (1) with: L λ = a 2 sin 2 (β) + b 2 cos 2 (β), ϕ ∈[0, 2π ], β ∈[0,π] and e ρ = cos(β)e 1 + sin(β)e 2 . The porosity is defined by: f = a 1 b 2 1 /(a 2 b 2 2

Fig. 1 .Fig. 1 .

 11 Fig. 1. The geometry of the cell considered: (left) prolate cavity; (right) oblate cavity. Fig. 1. Géométrie de la cellule considérée : (à gauche) cavité allongée ; (à droite) cavité aplatie.

e 3 (

 3 arctanh(e)e) (prolate)

Fig. 2 .

 2 Fig. 2. Comparison between isotropic and anisotropic yield locus for: (a) prolate cavity, a 1 /b 1 = 5, f = 0.1, (b) oblate cavity b 1 /a 1 = 5, f = 0.1.The results (Eq. (26)) are compared to the two field exact ones (full line).

Fig. 2 .

 2 Fig. 2. Comparaison des seuils pour le cas isotrope et pour le cas anisotrope : (a) cavité allongée, a1/b1 = 5, f = 0,1, (b) cavité aplatie b1/a1 = 5, f = 0,1. Les résultats (Éq. (26)) sont comparés aux résultats exacts de la méthode à deux champs (trait plein).

x 2 ∂Π∂A

 2 B = pσ 0 arcsinh(u 2 χ)-arcsinh(u 1 χ) Σ -3κ p 2 Σ * h H : Q =-.Using the following identity: A :H * : Q = (1 -9rh 2 4p2 )A : H : Q and combining relations (B.

5 )

 5 must be verified whatever the value of h 2 , four equations can be deduced from (B.5) and have to fulfilled by κ and r. The solution reads:

  1 + 3h 2 + 2h 3 ) = 4 5 (H 11 + H 22 + H 33 + H 44 + H 55 + H 66

The solutions are very closed to those obtained in[START_REF] Gologanu | Recent extensions of Gurson's model for porous ductile metals[END_REF], and can be deduced by developing hyperbolic terms which appear in their expressions, to the second order.

two fields criteria (26) (full line) and the numerical exact two field criteria (discrete points) on Fig 2 . The yield locus for prolate (a 1 /b 1 = 5) and oblate (b 1 /a 1 = 5) cavities are represented on Figs. 2(a) and 2(b) respectively, for either the particular case of the isotropic von Mises matrix and for the anisotropic Hill matrix. Anisotropy coefficients used for the numerical applications are: M 11 = 0.733, M 22 = 0.57, M 33 = 0.499, M 44 = 3.669, M 55 = 1.141, M 66 = 2.2. The data considered come from [START_REF] Benzerga | Plastic potentials for anisotropic porous solids[END_REF]. The results clearly show the effects of the plastic anisotropy of the matrix. However, a complete validation of the new results (Eq. ( 26)) requires a larger number of velocity fields and probably finite elements computations. This can be the purpose of a forthcoming study.

Appendix A. Expression of components of H

For completeness, let us first indicate the relations between the components of H and the ones of M. Let us first note that in the local frame of the ellipsoidal void ( e 1 ,e 2 ,e 3 ) the components of the tensors are obtained from usual transformation rules:

, where H 0 and M 0 are the anisotropic fourth order tensors defined in the frame related to the orthotropic directions of the matrix. The relations between the nonzero components of H 0 and those of M 0 (directly obtained from ( 5)) are (using Voigt's notation):


Vincent Monchiet
email: vincent.monchiet@ed.univ-lille1.fr

Eric Charkaluk
email: eric.charkaluk@univ-lille1.fr

Djimedo Kondo
email: kondo@univ-lille1.fr

Plasticity-damage based micromechanical modelling in high cycle fatigue

Keywords: Fatigue, Micromechanics, Unlimited endurance, PSB Fatigue, Micromécanique, Endurance illimitée, BGP

de documents 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

Structures and mechanical components submitted to cyclic loadings can fail after a finite number of cycles, even if the macroscopic response remains elastic. Fracture is the result of fatigue macrocracks nucleation and growth. Classically, in uniaxial experiments, a fatigue limit can be distinguished, i.e., there is a maximum value of the stress amplitude under which no fatigue crack can be observed whatever the applied number of cycles. The effect of a mean stress superimposed to a periodic loading has been brought to the fore by Gerber and Goodman [START_REF] Suresh | Fatigue of Materials[END_REF] who postulated constant life diagrams as a relationship between the amplitude and the mean part of the stress.

In the case of a multiaxial stress state, several fatigue criteria have been proposed in the last few decades, using mainly macroscopic stress invariant approaches [START_REF] Crossland | Effect of large hydrostatic pressures on the torsional fatigue strength of an alloy steel[END_REF][START_REF] Sines | Behavior of metals under complex static and alternating stresses[END_REF]. It has been earlier recognized that, under fatigue loadings, metals fracture is the result of the cyclic plastic strain [START_REF] Orowan | Theory of the fatigue of metals[END_REF] and in high cycle fatigue (HCF), at the microscale, some grains (in very low proportion), favorably oriented, undergo plasticity. This local plasticity is however confined and has no significant effect on the macroscopic response which remains elastic. Following Orowan [START_REF] Orowan | Theory of the fatigue of metals[END_REF], Dang Van [5] and then Papadopoulos [START_REF] Papadopoulos | Fatigue polycyclique des métaux : une nouvelle approche[END_REF] (see also [START_REF]High Cycle Metal Fatigue, From Theory to Applications[END_REF]) formalized this multiscale framework. The resulting criterion is based on the consideration that a sufficient condition for non nucleation of fatigue cracks is to ensure elastic shakedown at all the scales. Nevertheless, plastic activity alone is not able to account for the role of pressure on the fatigue criterion. In order to take into account the mean stress effect, Papadopoulos [START_REF] Papadopoulos | Fatigue polycyclique des métaux : une nouvelle approche[END_REF] postulated a 'generalized multiaxial fatigue criterion' (which corresponds to the elastic shakedown limit) depending on the hydrostatic pressure.

In fact, the phenomenon of grain level plasticity is accompanied with the development of localized damage. The explicit incorporation of damage in HCF modelling has been the subject of several recent studies among which [START_REF] Lemaitre | A two scale damage concept applied to fatigue[END_REF] and [START_REF] Flacelière | Contribution à la modélisation du dommage en fatigue multiaxiale d'un acier C36-confrontation à l'expérience[END_REF]. Although these works allow to account for the role of the pressure in fatigue, the principal criticism is that the concept of a scalar variable d used to describe damage cannot be explicitly related to a particular physical micromechanism. Moreover, in the two scale model proposed by [START_REF] Lemaitre | A two scale damage concept applied to fatigue[END_REF], no plasticity occurs at the grain scale below the macroscopic fatigue limit, and this does not correspond to experimental observations [START_REF] Cugy | Microplasticity and temperature increase in low carbon steels[END_REF].

Based on experimental results obtained in the context of faced-centered-cubic structures (FCC) [START_REF] Essmann | Annihilation of dislocations during tensile and cyclic deformation and limits of dislocation densities[END_REF][START_REF] Essmann | A model of extrusions and intrusions in fatigued metals. I. Point-defects production and the growth of extrusions[END_REF], and following the multiscale approach proposed by Dang Van [START_REF] Van | Sur la résistance à la fatigue des métaux[END_REF] (see also [START_REF]High Cycle Metal Fatigue, From Theory to Applications[END_REF]) in HCF, we develop a micromechanical approach of HCF criterion in which plasticity and damage are coupled at the grain scale (Section 2). Considering such coupled micro-mechanisms, the present note proposes a fatigue criterion in which an explicit dependence on the pressure and on the mean stress is shown (Section 3). Finally, through an application to affine cyclic loadings, we investigate the role of the homogenization scheme and the local hardening rule in the qualitative predictions of the approach (Section 4).

Cyclic plasticity and damage at grain level

Under low level of macroscopic loading as in the context of HCF, in the case of FCC structures, the plastic behavior is generally characterized by the activation of a predominant slip system and more precisely by the formation of strain localization band which are the potential sites for microcracks nucleation. Such bands, called persistent slip bands (PSB), have been amply studied by [START_REF] Essmann | Annihilation of dislocations during tensile and cyclic deformation and limits of dislocation densities[END_REF] and [START_REF] Essmann | A model of extrusions and intrusions in fatigued metals. I. Point-defects production and the growth of extrusions[END_REF]; in particular, it is shown that dislocations concentration, which is at the origin of the damage process, is particularly important in those bands.

Monocrystal plasticity

The predominant activated slip system is characterized by a plane, defined by a unit normal vector n and a slip direction m. Let us introduce σ and ε respectively the microscopic (i.e., at the grain level) stress and strain fields. As classically, an additive decomposition of the total strain ε into elastic strain, ε e , and plastic strain, ε p ,a tt h e microscopic scale is adopted: ε = ε e + ε p . The monocrystal plastic behavior is described by the plastic criterion f corresponding to Schmid's law, written in the classical form:

f =|τ -x|-τ 0 -R (1)
τ represents the shear stress acting on the slip plane: τ = σ : , being the second order symmetric tensor =

n s ⊗ m = 1 2 (n ⊗ m + m ⊗ n); τ 0 is the critical shear stress.
x and R are the hardening variables defining respectively:

• A kinematic linear rule: x = cγ p , where γ p is such that: ε p = γ p .

• An isotropic linear or non linear rule:

R(p) = R 0 p R(p) = R s 1 -exp(-rp) (2)
p is the cumulated plastic slip such that: p = t 0 |γ p (t ′ )| dt ′ ; R 0 , R s and r are three parameters.

Microdamage within localization bands

The strain localization into the PSBs is also accompanied by a dislocations annihilation mechanism. Transmission electron microscopy studies [START_REF] Essmann | A model of extrusions and intrusions in fatigued metals. I. Point-defects production and the growth of extrusions[END_REF] reveal that there is a critical spacing between dislocations at which their annihilation is favored. Dislocations will annihilate to form a vacancy or an interstitial point defect. According to Essmann et al. [START_REF] Essmann | A model of extrusions and intrusions in fatigued metals. I. Point-defects production and the growth of extrusions[END_REF], vacancies production is the dominant points defects generation process. The transition from vacancy production to the formation of microvoids or microcracks along the PSBs is not yet well understood. Therefore, we assume that the formation of these microvoids is the result of the agglomeration and the growth of vacancies formed by the dislocation annihilation process. It follows that damage along PSBs is the result of both two mechanisms: vacancies production and voids growth which is the result of the combined effect of the slip-like plastic activity and pressure.

A phenomenological model for points defects production by dislocations annihilation has been already proposed by Essmann et al. [START_REF] Essmann | Annihilation of dislocations during tensile and cyclic deformation and limits of dislocation densities[END_REF]. For this purpose, introducing two parameters ρ 0 and k a , let us first define the dislocation density, ρ p , and the annihilation density, ρ a :

ρ p = ρ 0 k a p; ρ a = ρ 0 k a p -1 + exp(-k a p) (3)
Noting that dislocations annihilation is not an irreversible mechanism (points defects could be mutually annihilated), a points defects density law ρd was introduced:

ρd =ρ a -k e ρ d ṗ (4)
Following [START_REF] Essmann | A model of extrusions and intrusions in fatigued metals. I. Point-defects production and the growth of extrusions[END_REF], it is assumed that all points defects are of vacancy-type. Denoting by V 0 the volume occupied by a vacancy, integrating (4) and taking into account (3), the porosity associated to the mechanism of vacancies production, η a , is given by:

η a (p) = V 0 ρ d = A 0 k a k e (k a -k e ) k a 1 -exp(-k e p) -k e 1 -exp(-k a p) (5)
with k e a constant and A 0 = V 0 ρ 0 . For very large value of p, η a reaches a limit η s a = A 0 k a /k e : there is an equilibrium between points defects annihilation and production. Let us introduce the total porosity at the grain scale:

η = η a + η g (6)
where η g is the porosity associated to the mechanism of voids growth as previously mentioned.

As voids growth induces volume change, the plastic strain at the grain scale can be decomposed in the following form, 1 being the first order identity tensor:

ε p = γ p + ε p h 1 ( 7
)
where the volumetric plastic strain ε p h , due to voids growth and deduced from the mass balance equation, is given by:

η g = 1 -exp -3ε p h (8)
As a first approximation, voids growth is described by means of a modified Rice and Tracey-type model [START_REF] Rice | On the ductile enlargement of voids in triaxial stress fields[END_REF]. For this purpose, we assume spherical voids and for simplicity we replace the crystal by an equivalent von Mises material. The proposed modification leads to (see Appendix A.7):

εp h = 1 2 √ 3 η sinh √ 3 2 σ h τ 0 ṗ (9)
where σ h = tr σ /3.

A local approach of fatigue crack nucleation at the PSB-matrix interface consists then in introducing a critical value of the porosity, η c :

η a + η g = η c (10)
The damage-dependent criterion [START_REF] Cugy | Microplasticity and temperature increase in low carbon steels[END_REF] is the threshold below which no crack nucleation could be observed. From ( 5), ( 8) and ( 9), it is readily seen that this criterion depends explicitly on the hydrostatic pressure σ h , on the cumulated plastic strain p and consequently on the local shear stress τ .

Micro-macro based determination of the fatigue criterion

Expression of the local fatigue criterion

For the determination of the fatigue criterion, following [START_REF] Van | Sur la résistance à la fatigue des métaux[END_REF], we assume that a necessary condition for non nucleation of cracks is that the monocrystal shakedowns elastically. For [START_REF] Papadopoulos | Fatigue polycyclique des métaux : une nouvelle approche[END_REF] (see [START_REF]High Cycle Metal Fatigue, From Theory to Applications[END_REF]), this is also a sufficient condition; then, elastic shakedown limit coincides with the fatigue limit. In order to take into account the dependence on the hydrostatic pressure, a 'generalized multiaxial fatigue limit' has been postulated by [START_REF] Papadopoulos | Fatigue polycyclique des métaux : une nouvelle approche[END_REF].

In our HCF approach, elastic shakedown is a necessary but not a sufficient condition for non nucleation of cracks: in the elastic shakedown regime, a crack may nucleate if a critical porosity along the PSB-matrix interface is attained (see condition [START_REF] Cugy | Microplasticity and temperature increase in low carbon steels[END_REF]). Therefore, as the grain shakedowns elastically, p is bounded and so the porosity is very low; the following approximations can then be made:

(i) η g ≪ 1, then mass balance equation ( 8) reads: η g ≃ 3ε p h . (ii) η a is low compared to its saturation value η s a ; then, point defect annihilation process is neglected (k e = 0) particularly in Eq. ( 4).

Taking into account these approximations in ( 5) and ( 8), the local fatigue criterion, deduced from [START_REF] Cugy | Microplasticity and temperature increase in low carbon steels[END_REF], reads:

A 0 k a p -1 + exp(-k a p) + 3ε p h <η c (11)
Since p and ε p h are defined at the grain scale, in order to establish a macroscopic expression of the fatigue criterion, a non linear homogenization approach is due.

The macroscopic fatigue criterion

Let us first recall that the macroscopic behavior remains elastic. Therefore, the macroscopic stress tensor, , and the macroscopic elastic strain tensor, E, are given by: = C : E where C is the stiffness tensor of the matrix. For simplicity, C is assumed isotropic: C = 3kJ + 2µK, k and µ being respectively the bulk and shear modulus of the matrix, J = 1 3 1 ⊗ 1 and K = I -J with I the fourth order symmetric identity tensor. Moreover, as a first approximation, it is assumed that the elastic behavior of the continuum and of its constituent grains are the same and reads σ = C : ε e .

As a first approach, Lin-Taylor or Sachs models could be considered:

σ = -C : ε p (Lin-Taylor) σ = (Sachs) (12)
A more elaborated, but simple, scheme which can be used to describe crystalline plasticity is the self consistent model (SCM) of Kröner [START_REF] Kröner | Zur plastichen verformung des vielkristalls[END_REF] (see also [START_REF] Bornert | Homogénéisation en mécanique des matériaux, Tome 2, Comportements non linéaires et problèmes ouverts[END_REF]). In the particular case of high cycle fatigue, the plastic strain being confined into a few number of grains embedded in the elastic matrix, the overall behavior is close to the one of polycrystals near the yield limit. For such situation, the Kröner SCM is suitable and gives:

σ = -C : (I -P : C) : ε p ( 13
)
where P is the fourth order Hill tensor; considering a spherical inclusion, P reads:

P = a 3k J + b 2µ K with a = 3k 3k + 4µ and b = 6 5 k + 2µ 3k + 4µ (14)
The interaction laws, given by ( 12) and by ( 13) with ( 14) can be put in the general following form:

σ = -C * : ε p (15)
with C * = 0 for Sachs's model, C * = C for Lin-Taylor's model and C * = C : (I -P : C) for Kröner scheme. 1It follows that the hydrostatic pressure σ h and the shear stress τ , which enter into the local criterion, [START_REF] Cugy | Microplasticity and temperature increase in low carbon steels[END_REF] are given by:

τ = : -2µ * γ p ; σ h = Σ h -3k * ε p h (16) with µ * = µ(1 -b), k * = k(1 -a).
Lin-Taylor estimate is obviously recovered for a = b = 0 and the Sachs model for a = b = 1.

As previously shown, the translation at the macroscopic level of the fatigue criterion [START_REF] Essmann | Annihilation of dislocations during tensile and cyclic deformation and limits of dislocation densities[END_REF] requires the expression of p and ε p h as functions of the macroscopic fields. p is related to τ by means of the Schmid's law (1) and then to the macroscopic stress by the interaction law [START_REF] Cano | Critère de fatigue polycyclique pour des matériaux anisotropes: application aux monocristaux[END_REF]. ε p h is obtained by integrating [START_REF] Flacelière | Contribution à la modélisation du dommage en fatigue multiaxiale d'un acier C36-confrontation à l'expérience[END_REF] in which the expression of σ h ,given by ( 16), is used:

εp h = 1 2 √ 3 η sinh √ 3 2 Σ h -3k * ε p h τ 0 ṗ (17)
In a general loading case, p and ε p h can be determined by using a numerical procedure. For an illustration purpose, we consider now the particular case of affine loadings.

Applications to affine loadings

The macroscopic stress tensor, (t), associated to an affine loading, is decomposed into an alternative and a mean part: (t) = a sin(ωt) + m .

When the grain shakes down elastically, the mean shear stress and the amplitude of the shear stress are given, using [START_REF] Suresh | Fatigue of Materials[END_REF], by: τ m = cγ p m and τ a = τ 0 + R(p). The interaction law [START_REF] Cano | Critère de fatigue polycyclique pour des matériaux anisotropes: application aux monocristaux[END_REF], for the affine loading, reads:

m : = (c + µ * )γ p m ; a : = τ 0 + R(p) (18)
Consider a first situation of an alternated torsion with a mean torsion: T(t)= T a sin(ωt) + T m . Using (18), p and γ p m , are given, in the favorably oriented grain, by:

γ p m = T m c + µ * ; ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ p =
T aτ 0 R 0 (For the linear isotropic hardening rule)

p = 1 r ln R s R s + τ 0 -T a
(For the non-linear isotropic hardening rule)

Since Σ h = 0, one has ε p h = 0 and η g = 0. Consequently, the fatigue criterion [START_REF] Essmann | Annihilation of dislocations during tensile and cyclic deformation and limits of dislocation densities[END_REF] is found to depend only on the variable p; then a mean torsion has no effect on the fatigue limit which is consistent with experimental data provided by [START_REF] Sines | Behavior of metals under complex static and alternating stresses[END_REF]. Let us consider now a second case corresponding to an alternated torsion loading with a mean pressure ( m = Σ h,m 1). For Sachs model, a closed form expression of the fatigue criterion is found:

2 √ 3 A 0 k a H m ( √ 3H m + 2k a ) 2k a exp √ 3 2 H m p -1 + √ 3H m exp{-k a p}-1 <η c (20)
with

H m = sinh{ √ 3 2 Σ h,m τ 0 }.
For Lin-Taylor model, numerical resolution of Eq. ( 17) shows that, for moderate values of the pressure, ε p h reaches a saturation state corresponding to εp h = 0, therefore, σ h,m = Σ h,m -3kε p h = 0 and exp(-k a p) can be neglected in [START_REF] Essmann | Annihilation of dislocations during tensile and cyclic deformation and limits of dislocation densities[END_REF]. An analytical expression of the criterion is then obtained:

A 0 k a p + Σ h,m k <A 0 + η c (21)
It is worth noticing that, for a linear isotropic hardening assumption (see Eq. ( 19)), the fatigue criterion for the Lin-Taylor scheme, reads:

T a τ 0 + α Σ h,m τ 0 <β; with α = R 0 A 0 k a k ; β = R 0 A 0 k a τ 0 η c + A 0 1 + τ 0 k a R 0 ( 22
)
This expression is in agreement with the fatigue criterion proposed by [START_REF] Van | Sur la résistance à la fatigue des métaux[END_REF]. βτ 0 corresponds to the fatigue limit in alternated torsion. Fig. 2 shows representation of the fatigue criterion for an alternated torsion with a mean pressure and for the different homogenization schemes. Figs. 2(a) and (b) correspond respectively to a linear and a non linear isotropic hardening. For the two hardening rules, the three homogenization schemes lead to different predictions, the Kröner model being between the two others. Besides the influence of the homogenization scheme and of the hardening rule, the results show a remarkable effect of the pressure on the fatigue criterion is shown. However, for negative pressure, the results of Sachs approach seem to be physically inconsistent. All those qualitative results have now to be confirmed by a comparison with experimental existing data.

Conclusion

A multiscale approach for the determination of a High Cycle Fatigue criterion has been proposed. It allows us to take into account, not only plasticity activity in some grains, but also damage due to microvoid growth along the PSB-matrix interface. The proposed criterion shows a significant effect of the mean pressure which is due to the considered damage mechanism. It predicts also that a mean shear stress has no effect on the fatigue limit; this qualitative prediction is in agreement with experimental data.

Fig. 1 .

 1 Fig. 1. Representative Volume Element (RVE).

Fig. 2 .

 2 Fig. 2. Application of the fatigue criterion in particular case of an alternated torsion loading with a mean pressure, T a /τ 0 versus Σ h, m /τ 0 : (a) linear isotropic hardening; (b) non-linear isotropic hardening. Parameters used are k = 70 000 MPa, ν = 0.33, A 0 = 2 × 10 -4 , k a = 20, R 0 = 350 MPa, R s = 420 MPa, r = 1.5, τ 0 = 70 Mpa, η c = 0.003.

An extension to the anisotropic behavior can be made, following a recent work of[START_REF] Cano | Critère de fatigue polycyclique pour des matériaux anisotropes: application aux monocristaux[END_REF].

Appendix A. A modified Rice and Tracey void growth model

The objective of this appendix is to derive a new approximation of a void growth model for arbitrary value of triaxiality. For this purpose, let us consider a spherical cavity (radius R), in an infinite perfectly plastic medium which is subjected to a uniform remote strain rate field ε∞ . Let us denote by r, θ , φ, the spherical coordinates and e r the radial unit vector. The strain rate field is decomposed in two terms, the remote field and a strain rate field (related to multiple slip), εD , which allows us to take into account void expansion: ε = ε∞ + εD , with ε∞ =γ . In order to apply the Rice and Tracey [START_REF] Rice | On the ductile enlargement of voids in triaxial stress fields[END_REF] result, we propose to replace the monocrystal by an equivalent von Mises material; this is justified by the activation of multiple slip systems around the cavity. Moreover, we consider a spherically symmetric strain void expansion field:

Therefore, the void growth law reads:

where, σ 0 being the yield stress of the equivalent von Mises material, such that σ 0 ε∞ eq = τ 0 ṗ, i.e., σ 0 = √ 3τ 0 , D is solution of (see [START_REF] Rice | On the ductile enlargement of voids in triaxial stress fields[END_REF]):

In (A.2), z = R 3 r 3 and m = 1 ε∞ eq ε∞ : e r ⊗ e r (-1/2 <m<1/2). The inner integration is over the unit sphere with dS = sin(φ) dφ dθ . Since (A.2) is not analytically integrable, Rice and Tracey proposed an approximation for high value of D, which is valid for large value of | σ h σ 0 |:

We propose here a new approximation for any value of σ h σ 0 . For this purpose, we develop (A.2), around m = 0, in a first time:

Using then the following relations: 1 4π S dS = 1, 1 4π S m dS = 0, 1 4π S m 2 dS = 1 5 , the integration of (A.2) gives:

Since the second term of the left hand of (A.5) is bounded (0

) and low compared to arcsinh(2D), it is neglected in the present study. Then D is given by:

The plastic strain related to the cavity growth, is given by: εp h = 1 3 Ω Ω , in which Ω is the volume of the monocrystal. Considering cavities with an average volume ω (ω = 4π 3 R 3 ), in an incompressible matrix, the mass balance equation, ω = Ω, reads (using (A.1)):


Xueting Sun
email: xuetingsun507@126.com

Jing Tao

Alain Roques
email: alain.roques@inrae.fr

Youqing Luo
email: youqingluo@126.com

Invasion History of Sirex noctilio Based on COI Sequence: The First Six Years in China

Keywords: Sirex noctilio, invasive species, genetic structure, haplotype network, mitochondrial marker

Sirex noctilio F. (Hymenoptera: Siricidae: Siricinae), a new invasive species in China, is a significant international forestry pest which, transported via logs and related wood packing materials, has led to environmental damage and substantial economic loss in many countries around the world. It was first detected in China in 2013, and since then infestations have been found in 18 additional sites. Using a 322 bp fragment of the mitochondrial barcode gene COI, we studied the genetic diversity and structure of S. noctilio populations in both native and invaded ranges, with a specific focus in China. Twelve haplotypes were found across the native and invaded distribution of the pest, of which three were dominant; among these there were only one or two mutational steps between each pair of haplotypes. No obvious genetic structure was found other than in Chinese populations. China has a unique and dominant haplotype not found elsewhere, and compared with the rest of the world, the genetic structure of Chinese populations suggested a multiple invasion scenario.

Introduction

Biological invasions are increasing exponentially as globalization progresses, without any sign of saturation, especially in insects [START_REF] Hulme | transport and trouble: Managing invasive species pathways in an era of globalization[END_REF][START_REF] Seebens | No saturation in the accumulation of alien species worldwide[END_REF]. At present, phytophagous species associated with woody plants constitute the majority of these insect invaders [START_REF] Roques | Temporal and interspecific variation in rates of spread for insect species invading Europe during the last 200 years[END_REF]. The problem warrants increased attention, as the invasion of alien species is likely to impact native biodiversity, and can also cause huge economic losses [START_REF] Simberloff | Nonindigenous Species: A Global Threat to Biodiversity and Stability; Nature and Human Society: The Quest for a Sustainable World[END_REF][START_REF] Clavero | Invasive species are a leading cause of animal extinctions[END_REF][START_REF] Everett | Patterns and pathways of biological invasions[END_REF][START_REF] Pimentel | Economic and environmental threats of alien plant, animal, and microbe invasions[END_REF]. There is a statistically significant correlation between the number of major invasive alien species in China and Chinese trade imports from the year 1980 to 2015 (Pearson correlation r = 0.915, p < 0.01, N = 35). The current primary route for alien forest species invading China is through human activity. This is mainly unintentional introduction; transported seedlings and their propagation materials carry 46% of human activity-generated invaders, with packaging materials or wood accounting for 21% (National Bureau of Statistics, http://www.stats.gov.cn/tjsj/tjgb/ ndtjgb/) [START_REF] Zhao | The identification of risk sources of forestry invasive alien species in China and their controlling countermeasures[END_REF]. Some major invasive species, such as the North American bark beetle Dendroctonus valens, fall webworm Hyphantria cunea and loblolly pine mealybug Oracella acuta, have already caused the destruction of forestry resources in China, negatively influencing species diversity and ecological security, and resulting in major economic loss [START_REF] Lu | Red Turpentine Beetle Dendroctonus valens LeConte[END_REF][START_REF] Cao | Rapid genetic structuring of populations of the invasive fall webworm in relation to spatial expansion and control campaigns[END_REF][START_REF] Chen | Potential distribution of the invasive loblolly pine mealybug, Oracella acuta (Hemiptera: Pseudococcidae), in Asia under future climate change scenarios[END_REF]. In order to monitor and manage the invasion of a pest, including possible biocontrol by natural enemies as well as the prevention of further introductions, it is essential to figure out its invasion history. However, observational records are usually difficult to get, and often biased, particularly when considering a recently-arrived invasive species. From this perspective, genetic analyses comparing populations in the native range with those in the invaded areas constitute an important tool for reconstructing the invasion history.

Molecular markers, such as mitochondrial, ribosomal, chloroplast, and microsatellite DNA, and others, have become increasingly popular for investigating phylogenetic relationships between insect populations. By using mitochondrial genes and microsatellite loci, [START_REF] Wei | Genetic structure and demographic history reveal migration of the diamondback moth Plutella xylostella (Lepidoptera: Plutellidae) from the southern to northern regions of China[END_REF] showed that the invasive diamondback moth, Plutella xylostella, moves within China from south to north, with rare effective travel in the reverse direction [START_REF] Wei | Genetic structure and demographic history reveal migration of the diamondback moth Plutella xylostella (Lepidoptera: Plutellidae) from the southern to northern regions of China[END_REF]. Similar patterns are found for the oriental fruit moth, Grapholita molesta, which originates in southern China and has dispersed to the north [START_REF] Wei | Population genetic structure and approximate B ayesian computation analyses reveal the southern origin and northward dispersal of the oriental fruit moth G rapholita molesta (L epidoptera: T ortricidae) in its native range[END_REF]. In the case of the western flower thrips, Frankliniella occidentalis, Cai (2017) recently determined the 13-year history of its invasion of China, revealing multiple introductions and further human-assisted spread, using mitochondrial genes and 18 microsatellite loci. In order to infer the invasion routes of the Asian long-horned beetle (Anoplophora glabripennis) in Europe, Javal (2019) used a 485 bp fragment of the mitochondrial barcode gene COI, which revealed multiple introductions from China and a bridgehead effect [START_REF] Javal | Complex invasion history of the Asian long-horned beetle: Fifteen years after first detection in Europe[END_REF]. The COI gene (mtDNA cytochrome oxidase I) has been of widespread use in tracing invasions [START_REF] Downie | Locating the sources of an invasive pest, grape phylloxera, using a mitochondrial DNA gene genealogy[END_REF][START_REF] Yang | Invasion genetics of the western flower thrips in China: Evidence for genetic bottleneck, hybridization and bridgehead effect[END_REF][START_REF] Auger-Rozenberg | Inferences on population history of a seed chalcid wasp: Invasion success despite a severe founder effect from an unexpected source population[END_REF][START_REF] Lippens | Genetic structure and invasion history of the house mouse (Mus musculus domesticus) in Senegal, West Africa: A legacy of colonial and contemporary times[END_REF].

Sirex noctilio, a woodwasp native to Eurasia and North Africa, is a significant forestry pest, subject to quarantine regulations around the world [START_REF] Spradbery | Aspects of the ecology of siricid woodwasps (Hymenoptera: Siricidae) in Europe, North Africa and Turkey with special reference to the biological control of Sirex noctilio F. in Australia[END_REF]. Its larvae tunnel the trunks of conifer trees and infect the wood with a damaging fungus, ultimately resulting in tree death. The first record of S. noctilio outside its native Eurasian range was in the southern hemisphere during the 1900s [START_REF] Boissin | Retracing the routes of introduction of invasive species: The case of the Sirex noctilio woodwasp[END_REF]. Following introduction to New Zealand, it spread rapidly through the stands of Pinus radiata planted there, killing trees across approximately 120,000 ha of the north-central region and thereby causing huge economic and ecological loss. During 1987-1989, an outbreak in Australia and Tasmania killed more than 5 million P. radiata trees [START_REF] Haugen | National Strategy for Control of Sirex noctilio in Australia[END_REF]; it was first found there 35 years ago [START_REF] Neumann | The sirex wood wasp in Australian radiata pine plantations[END_REF]. Several areas of South America have also been attacked, including Uruguay in 1980, Argentina in 1985, Brazil in 1988, and Chile in 2000. The mortality rate of pine trees was high, at 60% and 80% in Argentina and Uruguay, respectively. In Brazil, a 350, 000 ha pine plantation was destroyed. In 1994, it was reported in South Africa and in 2006 in the USA and Canada.

In China, S. noctilio was first detected in 2013 in the northeast, at Daqing, Heilongjiang Province [START_REF] Li | Detection and identification of the invasive Sirex noctilio (Hymenoptera: Siricidae) fungal symbiont, Amylostereum areolatum (Russulales: Amylostereacea), in China and the stimulating effect of insect venom on laccase production by A. areolatum YQL03[END_REF][START_REF] Sun | Identification of Sirex noctilio (Hymenoptera: Siricidae) using a species-specific cytochrome C oxidase subunit I PCR assay[END_REF]. Six subsequent years of field investigations uncovered 18 new occurrences, but these were followed by a post-invasion latency period, or time lag. Therefore, the first place where alien forest pests are found does not necessarily represent their point of invasion. However, the relationship between the initial discovery site and the invasion site is typically the closest and most direct. It is generally the case that the site of initial discovery contains the invading colony [START_REF] Zhao | The identification of risk sources of forestry invasive alien species in China and their controlling countermeasures[END_REF]. Moreover, when the invasion site's latitude and habitat are similar to those of the pest's origin, the possibility of invasion and colonization is increased [START_REF] Zhao | The identification of risk sources of forestry invasive alien species in China and their controlling countermeasures[END_REF]. S. noctilio is found in China at similar latitudes to its native sites; moreover, its native regions in Russia are near the Chinese border. Based on climatic conditions and host distribution, domestic and foreign experts analyzed the potential woodwasp distribution area in China. Their results indicated that a band stretching from Yunnan to Heilongjiang represents a highly suitable area for its establishment [START_REF] Dai | Pest Risk Analysis in China and the Economic Impacts Assessment of the Trade Measures for Sirex noctilio[END_REF]. Its potential host tree species are widely distributed in China; hence, there is a high potential invasion risk into additional areas.

In order to investigate the origins and pathways of invasion, Bernard (2012) used microsatellite and sequence data collected from five continents. Two major sources were identified, although the invasion history was quite complex, with most of the populations containing admixtures of independently introduced European insects that subsequently spread through the invaded regions. Similarly, Bittner (2017) investigated the genetic diversity of 924 S. noctilio in nine populations from New York and Pennsylvania (US), Ontario (Canada), and Queensland (Australia), showing that multiple introductions occurred in northeastern North America. However, a more comprehensive understanding requires the collection of more samples from its native Europe and Asia.

In this study, the genetic structure and invasion routes of S. noctilio populations were surveyed using the mtDNA cytochrome oxidase 1 (COI) gene. Our objective was to investigate the genetic structure and genetic diversity of S. noctilio populations around the world, especially in relation to Chinese infestations. We compared patterns in invaded continents and native regions, analyzing for the first time 778 individuals from 6 continents. In particular, samples from original European populations and recent Chinese infestations were included. Overall, this study aims to elucidate the invasive routes of this pest as well as provide data enabling the prevention of the harmful consequences of S. noctilio spread.

Materials and Methods

Sample Collection and Selection

S. noctilio samples were collected from native European, non-Chinese invaders, and Chinese populations, comprising 44, 127, and 607 specimens, respectively. A summary of sampling sites and common parameters of the genetic diversity of S. noctilio in this study is provided in Table 1. GenBank/DRYAD provided 23 of the 44 native European specimens; others were dry or anhydrous ethanol-preserved individuals borrowed from museums or donated by local institutions. Non-Chinese invader populations were represented by 22, 39, 50 and 16 GenBank/DRYAD samples from Africa, North America, South America, and Oceania, respectively. All Genbank/DRYAD samples are retrieved from DRYAD: https://datadryad.org/stash/dataset/doi:10.5061/dryad.37mm8; from Genebank: https: //www.ncbi.nlm.nih.gov/nuccore/?term=sirex+noctilio. Finally, we obtained a total of 607 specimens from 19 Chinese populations (Table 1 and Figure 1), donated by various local forestry bureaus or obtained from the standard quarantine facility of Beijing Forestry University, which contains a large number of locally collected infested stems; these specimens were identified using taxonomic literature on Sirex woodwasps before DNA extraction [START_REF] Schiff | Hymenoptera: Symphyta: Siricoidea) of the western hemisphere[END_REF]. As the Chinese source areas were infested at different times and to different degrees, the number of specimens and sampling year vary. Dried museum samples were handled outside our laboratory, as described below; specimens already in 99.7% pure anhydrous alcohol remained so; all remaining specimens were frozen at -80 • C.

Insects 2020, 11, x FOR PEER REVIEW 3 of 14 the first time 778 individuals from 6 continents. In particular, samples from original European populations and recent Chinese infestations were included. Overall, this study aims to elucidate the invasive routes of this pest as well as provide data enabling the prevention of the harmful consequences of S. noctilio spread.

Materials and Methods

Sample Collection and Selection

S. noctilio samples were collected from native European, non-Chinese invaders, and Chinese populations, comprising 44, 127, and 607 specimens, respectively. A summary of sampling sites and common parameters of the genetic diversity of S. noctilio in this study is provided in Table 1. GenBank/DRYAD provided 23 of the 44 native European specimens; others were dry or anhydrous ethanol-preserved individuals borrowed from museums or donated by local institutions. Non-Chinese invader populations were represented by 22, 39, 50 and 16 GenBank/DRYAD samples from Africa, North America, South America, and Oceania, respectively. All Genbank/DRYAD samples are retrieved from DRYAD: https://datadryad.org/stash/dataset/doi:10.5061/dryad.37mm8; from Genebank: https://www.ncbi.nlm.nih.gov/nuccore/?term=sirex+noctilio. Finally, we obtained a total of 607 specimens from 19 Chinese populations (Table 1 and Figure 1), donated by various local forestry bureaus or obtained from the standard quarantine facility of Beijing Forestry University, which contains a large number of locally collected infested stems; these specimens were identified using taxonomic literature on Sirex woodwasps before DNA extraction [START_REF] Schiff | Hymenoptera: Symphyta: Siricoidea) of the western hemisphere[END_REF]. As the Chinese source areas were infested at different times and to different degrees, the number of specimens and sampling year vary. Dried museum samples were handled outside our laboratory, as described below; specimens already in 99.7% pure anhydrous alcohol remained so; all remaining specimens were frozen at -80 °C. As there is an obvious sample size bias towards the six Chinese populations with high sample numbers (population codes DM, QQ, HG, FJ, JBT, CT), we used only 20 of the sequenced individuals, which included all the haplotypes found in each population. We selected individuals from each haplotype based on the selection ratio and the number of samples in each haplotype. The final total of 364 samples for analysis consisted of 44, 127, and 193 individuals from Europe, non-Chinese invaded areas, and China, respectively.

DNA Sequence Analysis

For the native European samples not supplied as sequences directly from GenBank, DNA extraction, PCR amplification and sequencing of the barcode fragment were performed for 15 specimens at the Laboratory of Forest Zoology URZF, INRAE (Orléans, Paris, France). DNA extracts were prepared from one adult hind leg or thorax muscle using NucleoSpin tissue XS kits (Macherey Nagel, Dylan, Germany) according to the manufacturer's protocol. The COI barcoding fragment was amplified via PCR using the forward primer LCO1490 (5 -GGT CAA CAA ATC ATA AAG ATA TTG G-3') and the reverse primer HC02198, (5'-TAA ACT TCA GGG TGA CCA AAA AAT CA-3') (Vrijenhoek, Rijswijk, Netherlands, 1994). PCR reagents were used and reactions performed following Sun [START_REF] Sun | Identification of Sirex noctilio (Hymenoptera: Siricidae) using a species-specific cytochrome C oxidase subunit I PCR assay[END_REF]. PCR products were purified using the QIAquick PCR Purification Kit (Qiagen, Hilden, Germany). Sequencing was performed using the Sanger method with an ABI Prism BigDye Terminator v3.1 cycle sequencing kit (Thermo-Fisher, Waltham, MA, USA). The reagents and procedure used were: 1 µL 5 × bigdye buffer, 1 µL bigdye, 1 µL LCO/HCO, diluted to 3.3 µM, 2 µL DNA; 30 cycles of 10 s at 96 • C, 5 s at 50 • C, 3 min at 60 • C. The remaining 6 specimens were borrowed from museums and had to remain intact; for these the same COI fragment was amplified at the Canadian Centre for DNA Barcoding (CCDB-Biodiversity Institute of Ontario, University of Guelph) using a slightly different sequencing set (C-LepFolF/C-LepFolR), following the standard high-throughput protocol [START_REF] De Waard | Methods in Molecular Biology: Environmental Genetics[END_REF].

DNA was extracted from the leg or thorax muscle of each Chinese specimen using the Multisource Genomic DNA Miniprep Kit (Axygen, San Francisco, CA, USA), following the manufacturer's protocol to prevent cross-contamination, all tools were sterilized by flame or 75% pure anhydrous ethanol. DNA was eluted using 70 µL of elution buffer and stored at -20 • C until use. Species-specific COI amplification primers were synthesized by a commercial company (SinoGenoMax, Beijing, China). PCR primers LCO1490 and HC02198 were used as above. Reactions were performed in a total volume of 25 µL, containing 12.5 µL of 2 × GoTaq ® Green Master Mix (Promega, Madison, WI, USA), 1 µL of each primer, 1.5 µL of DNA template, and 9 µL of ultrapure water. The reaction procedure was as follows: initial denaturation for 2 min at 94 • C, followed by 35 cycles of denaturation at 94 • C for 30 s, annealing at 45 • C for 30 s, and extension at 68 • C for 2 min, with a final extension at 68 • C for 10 min. PCR products (3 µL) were analyzed by electrophoresis on a 1.5% (w/v) agarose gel (1 × TAE), alongside a DNA marker (D2000, Takara, Kyoto, Japan). After electrophoresis at 130 V for 20 min, the PCR products were visualized using ethidium bromide and ultraviolet light. Finally, products clearly visible after electrophoresis were sent to a commercial company (SinoGenoMax, China) for sequencing in both directions.

All electropherograms were checked manually in CodonCode Aligner v3.7.1. (CodonCode Corporation, Centerville, MA, USA) to assess their quality. The alignment was done independently using Clustal W implemented in MEGA v6.06 with default parameters, and no pseudogenes or stop codons were detected in the sequences [START_REF] Thompson | Improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice[END_REF][START_REF] Tamura | MEGA6: Molecular evolutionary genetics analysis version 6.0[END_REF]. All COI sequences were trimmed to the same length (322 bp) for final alignment.

Details on the collecting data for each specimen, as well as photographs, sequence records, trace files, and primer sequences used for PCR amplification, together with GenBank accession numbers, are available through the dataset in BOLD (www.boldsystems.org).

Insects 2020, 11, 111 7 of 14

Statistical Analysis

Common parameters of genetic diversity, including the number of haplotypes, haplotype, and nucleotide diversity were calculated using DNAsp v6 [START_REF] Librado | DnaSP v5: A software for comprehensive analysis of DNA polymorphism data[END_REF]. The haplotype network was built in PopART [START_REF] Leigh | POPART: Full-feature software for haplotype network construction[END_REF][START_REF] Bandelt | Median-joining networks for inferring intraspecific phylogenies[END_REF][START_REF] Clement | Estimating gene genealogies[END_REF]. Haplotype distribution and frequency were projected on a map using ArcGIS 9.3 (ESRI, Redlands, CA, USA).

Analyses of molecular variance (AMOVA) were performed on native European populations and Chinese populations to measure the partitioning of genetic variation between populations and groups of populations in Arlequin 3.5.2.2 [START_REF] Excoffier | Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data[END_REF][START_REF] Excoffier | version 3.0): An integrated software package for population genetics data analysis[END_REF]. The first AMOVA was performed on four clusters of Chinese populations in the northeast of China defined by geographic proximity: western (mountain area), middle (city area), which was divided into two parts, and eastern (city and mountain areas). A second AMOVA was performed on two clusters of Chinese populations grouped by latitude, north and south, in order to test the trend of spread.

Results

There were 11 variable sites among all insects sampled. Three dominant haplotypes, H1, H2, and H9, were present among the total of 12 haplotypes; the haplotype diversity was 0.4207. Only one or two mutation sites existed across the different haplotypes (Figure 2). There were two mutation sites between H1 and H12. H1 and H9 were haplotypes that were made up of populations from different countries: H1 was widely located on all continents, whereas H9 was not found in China and Africa (Figures 2 and3). H2, a major Chinese haplotype along with H1, was found only in Chinese samples (Figure 2). H3, H4, H5, H6, H7, and H8 were unique to Chinese populations. H8, the haplotype found near the national border with Russia (Figure 4) is worth noting. H10 was unique to South America, as were H11 and H12 to Europe.

Insects 2020, 11, x FOR PEER REVIEW 7 of 14 [START_REF] Leigh | POPART: Full-feature software for haplotype network construction[END_REF][START_REF] Bandelt | Median-joining networks for inferring intraspecific phylogenies[END_REF][START_REF] Clement | Estimating gene genealogies[END_REF]. Haplotype distribution and frequency were projected on a map using ArcGIS 9.3 (ESRI, Redlands, CA, USA). Analyses of molecular variance (AMOVA) were performed on native European populations and Chinese populations to measure the partitioning of genetic variation between populations and groups of populations in Arlequin 3.5.2.2 [START_REF] Excoffier | Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data[END_REF][START_REF] Excoffier | version 3.0): An integrated software package for population genetics data analysis[END_REF]. The first AMOVA was performed on four clusters of Chinese populations in the northeast of China defined by geographic proximity: western (mountain area), middle (city area), which was divided into two parts, and eastern (city and mountain areas). A second AMOVA was performed on two clusters of Chinese populations grouped by latitude, north and south, in order to test the trend of spread.

Results

There were 11 variable sites among all insects sampled. Three dominant haplotypes, H1, H2, and H9, were present among the total of 12 haplotypes; the haplotype diversity was 0.4207. Only one or two mutation sites existed across the different haplotypes (Figure 2). There were two mutation sites between H1 and H12. H1 and H9 were haplotypes that were made up of populations from different countries: H1 was widely located on all continents, whereas H9 was not found in China and Africa (Figures 2 and3). H2, a major Chinese haplotype along with H1, was found only in Chinese samples (Figure 2). H3, H4, H5, H6, H7, and H8 were unique to Chinese populations. H8, the haplotype found near the national border with Russia (Figure 4) is worth noting. H10 was unique to South America, as were H11 and H12 to Europe. AMOVA analysis of the genetic structure of Chinese populations (Table 2) showed that when grouped by geographic proximity (Figure 5), most of the variation (38.64%, p < 0.01) was between geographic clusters, with a relatively small part of the genetic variance (9.15%, p < 0.01) between populations within geographic clusters. When sites were grouped by latitude (Figure 6), the variation (25.58%, p < 0.01) between clusters was slightly higher than that (22.46%, p < 0.01) between populations within clusters. Within this study, no clear genetic structure emerged within native European populations. AMOVA analysis of the genetic structure of Chinese populations (Table 2) showed that when grouped by geographic proximity (Figure 5), most of the variation (38.64%, p < 0.01) was between geographic clusters, with a relatively small part of the genetic variance (9.15%, p < 0.01) between populations within geographic clusters. When sites were grouped by latitude (Figure 6), the variation (25.58%, p < 0.01) between clusters was slightly higher than that (22.46%, p < 0.01) between populations within clusters. Within this study, no clear genetic structure emerged within native European populations. The Chinese sampling sites are divided into four groups based on geographical proximity: west for one group, middle for two groups and east for one group. The middle and eastern groups, where there are cities, contained 4 haplotypes, whereas the western group, which is mainly mountainous, had 3 haplotypes. The first site discovered in China, DM, contained 40% H1 and a similar proportion, 55%, of H2. There is no H1 in the western group, although the proportion of H1 in the two middle area groups is quite high as is also found in the eastern group. There is a decreasing trend from west to east in the proportion of H2. The representative populations on this axis were HDQ, DM, and HG with 100%, 55%, and 20% H2, respectively. Among the other six unique haplotypes in China, H8 is located in MZL, on the border between China and Russia. Four of the remaining five haplotypes are distributed in the populations from which large samples were taken and are located in city areas. The sixth, HHEJ, is located in the area from which P. sylvestris var. mongolica seedlings are taken for distribution and planting elsewhere. H3 is present in 5% of DM individuals, and so does H4, H5, H6, H7, these haplotypes occurred in 10%, 5%, 7%, and 25% of sampled insects of HG, FJ, YS and HHEJ, respectively. The latitude 46.204 • N was used to group samples into the categories north and south. Haplotype diversity in the northern area is obviously higher than that in the south. The northern cluster has seven haplotypes, H1, H2, H3, H4, H5, H7, and H8, from nine populations, compared to two haplotypes (H1 and H6) from ten populations in southern China.

Discussion

Host Unity of S. Noctilio in China

The host range of S. noctilio is rather large, although it is known that the pest is harmful to various kinds of conifers around the world, especially Pinus, Picea, Abies, Larix, and Pseudotsuga [START_REF] Spradbery | Aspects of the ecology of siricid woodwasps (Hymenoptera: Siricidae) in Europe, North Africa and Turkey with special reference to the biological control of Sirex noctilio F. in Australia[END_REF]. Both natural forests and pine plantations can suffer damage, especially when they are overstocked and stressed. In China, the origin of P. sylvestris var. mongolica is in the area of Honghuaerji in Inner Mongolia and north of the Huma River in Heilongjiang [START_REF] Kang | Natural distribution of Pinus sylvestris var. mongolica on sandy land and its cultivation as an exotic species[END_REF]. Field investigation found no S. noctilio in the wild distribution areas of P. sylvestris var. mongolica in either mountain or sandy ecological zones (Figure 1), other than two larvae found in the western, sandy area. This occurrence, and the unique haplotype, probably came from the considerable imports of Russian timber required for the construction of Honghuaerji Scenic Area.

Currently, the main infestation zones in China are those where P. sylvestris var. mongolica has been introduced. Two typical examples are Hegang (Heilongjiang) and Jinbaotun (Inner Mongolia). The former is a mixed forest of Pinus koraiensis, L. olgensis, Picea koraiensis and P. sylvestris var. mongolica. The latter contains Pinus tabuliformis and P. sylvestris var. mongolica. However, fieldwork has indicated that the only infested conifer species in China to date is P. sylvestris var. mongolica, even when other potential hosts are present in the same forest [START_REF] Xu | Volatile profiles of three tree species in the northeastern China and associated effects on Sirex noctilio activity[END_REF]. As is characteristic of invasive pests, S. noctilio appears to be harmless in the natural distribution area of its favorite host tree species.

Patterns of Multiple Invasion and Spread in China

There are currently three more frequent haplotypes, H1, H2, and H9, which have a worldwide distribution. It is thus unsurprising that they are in the center of the haplotype network and closely connected with other well-established haplotypes [START_REF] Freeland | Molecular Ecology[END_REF]. Interestingly, H1, despite being globally the most prevalent haplotype, is absent in western China, although strongly represented in the city regions of central and eastern China. This constitutes evidence that genetic communication is human activity-derived in both invasion and spread scenarios. Given that the western component is the primary seedling source of P. sylvestris var. mongolica, trees from this area will be planted elsewhere. We speculate that H1 is an invasion haplotype carried by humans that have entered the central and eastern urban areas of China. As it existed on all continents, its source is unclear. Based on our data, H2, H7 and H8 were found only in China, and given that they are unlikely to have originated there, we are unable to identify their origin. H2 is the main haplotype of the western component, with a distribution ratio that gradually decreases from west to east; this is in line with the characteristics of spread from a point of origin to other regions. This, as noted above, may represent an invasion originating with Russian wood imported for construction at Honghuaerji and spread eastward by planting. Severely endangered woodlands in Dumeng, Hegang, Fujin, and other eastern areas may have been invaded in this way. H7 and H8 are also unique to the western region; H8 occurs on the border with Russia.

Grouped by latitude, significantly more haplotypes exist in northern Chinese S. noctilio populations than in those to the south. Damage to woodlands in the northeastern region occurs discretely, in multiple forest farms, parks, and other wooded land. This damage was concentrated, as is characteristic of the middle to late stages of an invasion. In contrast, in newly-invaded southern areas, such as Jinbaotun, damaged forests were localized and individual affected trees were distributed discretely in the affected area, indicating the initial stage of introduction. Although the 1130 km between Honghuaerji and Jinbaotun is shorter than the 1400 km between Honghuaerji and the affected area in the northeast, as the haplotype H1 found in the southern region is not present in Honghuaerji, the south regional populations may originate from northern China.

S. noctilio has a relatively strong natural flight ability. In the southern hemisphere, the maximum spread rate of an invading population of S. noctilio is 82 km per year [START_REF] Lantschner | Temperature explains variable spread rates of the invasive woodwasp Sirex noctilio in the Southern Hemisphere[END_REF]. In China, the diffusion rate of the pest is 74.57 km/year [START_REF] Liu | Effects of age (in days), gender and nematode infection on the flight capacity of Sirex noctilio[END_REF]. The straight-line distances between Dumeng, the first site detected, and the other two severely endangered areas, Hegang and Jinbaotun, are 440 km and 390 km, respectively, which would take about 6 years to achieve at 74.57 km/year. Yet the woodwasp was already present in 2014 and 2015. The first recorded occurrence sites in all the affected areas are distributed along with the main traffic networks (Figure 1), which thus appear to have increased the impact caused by human activities, such as wood packaging and log transportation.

Focus for Further Research

When analyzing biological invasions, knowledge of the genetic structure of populations from the area of origin is vital in determining the source. Our target species is native to Eurasia and North Africa, but it does not cause serious harm there and has consequently received less attention, leading to difficulties in collecting samples. In particular, the lack of samples from the Russian region caused sampling breakpoints between the European and northeast Chinese populations. The observation that H2 occurred only in China may be due to the large sample size from some previously unsampled sites from this area. More reliable results may be obtained from sampling efforts focused on this region and on western China. Finally, microsatellite marker use involving cytochrome C oxidase subunit I (COI) could prove effective, showing its utility as a method for determining invasion pathways and spread.

Conclusions

In this study, a 322 bp fragment of the mitochondrial barcode gene COI was used to analysis the genetic diversity and structure of S. noctilio populations in both native and invaded ranges, with a specific focus in China. Twelve haplotypes were found across the native and invaded distribution of the pest, of which three were dominant; among these there were only one or two mutational steps between each pair of haplotypes. No obvious genetic structure was found other than in Chinese populations. China has a unique and dominant haplotype not found elsewhere, and compared with the rest of the world, the genetic structure of Chinese populations suggested a multiple invasion scenario.

Figure 1 .

 1 Figure 1. Sampled Chinese S. noctilio population locations (indicated by 2-4 letter combinations) and year of the first detection; natural distribution area of the host. The triangular marker for HHEJ indicates the Honghuaerji Scenic Area. The pentagram marker for DM indicates the place where S. noctilio was first detected in China.

Figure 1 .

 1 Figure 1. Sampled Chinese S. noctilio population locations (indicated by 2-4 letter combinations) and year of the first detection; natural distribution area of the host. The triangular marker for HHEJ indicates the Honghuaerji Scenic Area. The pentagram marker for DM indicates the place where S. noctilio was first detected in China.

Figure 2 .

 2 Figure 2. S. noctilio haplotype network based on 322 bp sequences, 364 samples around the world. Short solid line(s) between haplotypes show(s) the number of mutation sites.

Figure 2 .

 2 Figure 2. S. noctilio haplotype network based on 322 bp sequences, 364 samples around the world. Short solid line(s) between haplotypes show(s) the number of mutation sites.

Figure 2 .

 2 Figure 2. S. noctilio haplotype network based on 322 bp sequences, 364 samples around the world. Short solid line(s) between haplotypes show(s) the number of mutation sites.

Figure 3 .

 3 Figure 3. Venn diagram of S. noctilio haplotypes by continent. Numbers and colors correspond to the haplotype network shown in Figure 2.

Figure 3 .

 3 Figure 3. Venn diagram of S. noctilio haplotypes by continent. Numbers and colors correspond to the haplotype network shown in Figure 2. Insects 2020, 11, x FOR PEER REVIEW 8 of 14

Figure 4 .

 4 Figure 4. Global S. noctilio haplotype distribution. Colors correspond to the haplotype network shown in Figure 2.

Figure 4 .

 4 Figure 4. Global S. noctilio haplotype distribution. Colors correspond to the haplotype network shown in Figure 2.

Figure 5 .

 5 Figure 5. Chinese S. noctilio populations grouped by geographic proximity (west to east) and their haplotype projections.

Figure 5 .

 5 Figure 5. Chinese S. noctilio populations grouped by geographic proximity (west to east) and their haplotype projections.

Figure 5 .

 5 Figure 5. Chinese S. noctilio populations grouped by geographic proximity (west to east) and their haplotype projections.

Figure 6 .

 6 Figure 6. Chinese S. noctilio populations grouped by latitude (north to south) and their haplotype projections.

Figure 6 .

 6 Figure 6. Chinese S. noctilio populations grouped by latitude (north to south) and their haplotype projections.

Table 1 .

 1 Summary of sampling sites and common parameters of genetic diversity of Sirex noctilio analyzed in this study.

	Locality	Population Code	First Detection	Sampling Time	N(Sequenced Individuals)	N variable Sites	Nht	h	π
	Total				364(778)
	Asia, China				193(607)	7		0.311	0.00114
	Dumeng, Heilongjiang	DM	2013	2016	20(66)	2		0.549	0.00183
	Qiqihaer, Heilongjiang	QQ	2013	2016	20(64)	1		0.185	0.00057
	Hegang, Heilongjiang	HG	2014	2015-16	20(204)	2		0.472	0.00191
	Fujin, Heilongjiang	FJ	2014	2015	20(69)	2		0.344	0.00111
	Hulunbeier, Inner Mongolia	HHEJ, HLBE	2014	2014	8	1		0.4	0.00124
	Keyihe, Inner Mongolia	HDQ	2013	2013	1	0		0	0
	Yushu, Jilin	YS	2015	2015	14	1		0.138	0.00043
	Jinbaotun, Inner Mongolia	JBT	2015	2015-16	20(83)	0		0	0
	Changtu, Liaoning	CT	2016	2016	20(48)	0		0	0
	Changchun, Jilin	CC	2013	2013	1	0		0	0
	Zhanggutai, Liaoning	ZGT	2016	2016	18	0		0	0
	Kangping, Liaoning	KP	2016	2016	18	0		0	0
	Chifeng, Inner Mongolia	CF	2016	2016	2	0		0	0
	Acheng, Heilongjiang	AC	2013	2013	1	0		0	0
	Mudanjiang, Heilongjiang	MDJ	2013	2013	6	0		0	0
		Ningan, Heilongjiang MDJ-DJ		2013
		Suiyang, Heilongjiang MDJ-SY		2013
		Muleng, Heilongjiang MDJ-BMT		2013
		Hailin, Heilongjiang MDJ-CH		2013
		Mudanjiang, Heilongjiang MDJ-SGLK		2013
		Mudanjiang, Heilongjiang MDJ-XS		2013
	Harbin, Heilongjiang	XL	2013	2013	1	0		0	0
	Nongan, Jilin	NAN	2014	2014	1	0		0	0
	Jiamusi, Heilongjiang	HNC	2013	2013	1	0		0	0
	Manzhouli, Inner Mongolia	MZL	2013	2013	1	0		0	0
	Africa				22	0		0	0
	South Africa	SA			1	0		0	0
	Natal Midlands	NM			13	0		0	0
	Western Cape	WC			6	0		0	0
	Eastern Cape	EC			2	0		0	0
	North America				39	1		0.264	0.00082
	Canada Canada-DRYAD	Can Can			27	1		0.352	0.00109
	USA USA-DRYAD	NY, USA USA			12	0 0		0 0	0 0
	South America				50	2		0.563	0.00202

Table 1 .

 1 Cont.Numbers in parentheses are the number of actual sampled individuals from Chinese populations. Others are the sequences of the mitochondrial Cytochrome Oxydase subunit I (COI) gene used in this study based on 322 bp fragment. In bold, the specimens are amplified in our Chinese and French laboratories; Italics are the specimens amplified at the Canadian Centre for DNA Barcoding; others are retrieved from databases. n ht Number of haplotypes, h Haplotype diversity, π Nucleotide diversity.

	Locality	Population Code	First Detection	Sampling Time	N(Sequenced Individuals)	N variable Sites	Nht	h	π
	Argentina (unknow)	Arg
	Arg-DRYAD	Arg, ArgES
	Arg-Arroyo del medio	Arg		2018	25	2		0.568	0.00237
	Arg-Aluminé	Arg		2018
	Arg-Andacollo	Arg		2019
	Uruguay	Uruguay			1	0		0	0
	Chile Chil-DRYAD	Chi, Chil Chi, ChiS			24	1		0.511	0.00159
	Oceania				16	1		0.226	0.0007
	Australia Austr-DRYAD	Austra, Austr AustrES, AustrEB			16	1		0.226	0.0007
	Europe				44	4		0.533	0.00194
	Switzerland Swit-DRYAD	Switz Swit, SW, Swi			15	0		0	0
	Italy	Ita-XS			1	0		0	0
	Poland	Pol-XS			1	0		0	0
	Hungary	Hun, Hun-XS		2018	4,2	0		0	0
	Ger-Berlin	Ger-Ber			2	0		0	0
	Ger-Stuttgart	Ger-Stu			1	0		0	0
	Spain Spain-DRYAD	Spain Spain			3	0		0	0
	Fr-Labenne	Fr6-Lab., Fr-Lab.XS30		2017	1,1	0		0	0
	Fr-Meolans Revel	Fr03-MR			1	0		0	0
	Fr-Villes sur Auzon	Fr01-V.S.A			1	0		0	0
	Fr-Marseille	Fr3-Mar.		2013	1	0		0	0
	Fr-Tarascon	Fr11-Tara		2017	1	0		0	0
	Fr-Brain sur Allones	Fr8-B.S.A., Fr-B.S.A.XS31		2017	1,1	1		0.667	0.00207
	Fr-Etaules	Fr4-Eta.		2017	1	0		0	0
	Fr-Lavercantière	Fr18-Lav.		2011	1	0		0	0
	Fr-Orléans	Fr21-Orle, Fr87-Orle		2017.2	2	1		0.667	0.00207
	Fr-St Alban des Hurtiè	Fr16-St.Al		2017	1	0		0	0
	Fr-Roissy CDG	Fr33-CDG		2017	1	0		0	0
	Fr-Tartigny	Fr12-Tart		2016	1	0		0	0

Table 2 .

 2 AMOVA analysis of genetic variation.

	Source of Variation	d.f.	Sum of Squares	Variance Components	Percentage of Variation	P
	AMOVA-Chinese populations grouped by geographic proximity (West-East)
	Among groups	3	11.561	0.08286	38.64 ***	0.006
	Among populations within groups	15	4.256	0.01961	9.15 ***	0.003
	Within populations	174	19.479	0.11195	52.21 ***	0.000
		AMOVA-Chinese populations grouped by latitude (North-South)
	Among groups	1	6.262	0.05510	25.58 ***	0.004
	Among populations within groups	17	9.555	0.04838	22.46 ***	0.000
	Within populations	174	19.479	0.11195	51.97 ***	0.000

*** p < 0.01; ** p < 0.05; * p < 0.1.

Table 2 .

 2 AMOVA analysis of genetic variation.

	Source of Variation	d.f.	Sum of Squares	Variance Components	Percentage of Variation	P
	AMOVA-Chinese populations grouped by geographic proximity (West-East)
	Among groups	3	11.561	0.08286	38.64 ***	0.006
	Among populations within groups	15	4.256	0.01961	9.15 ***	0.003
	Within populations	174	19.479	0.11195	52.21 ***	0.000
	AMOVA-Chinese populations grouped by latitude (North-South)
	Among groups	1	6.262	0.05510	25.58 ***	0.004
	Among populations within groups	17	9.555	0.04838	22.46 ***	0.000
	Within populations Insects 2020, 11, x FOR PEER REVIEW	174	19.479	0.11195	51.97 ***	0.000	9 of 14

*** p < 0.01; ** p < 0.05; * p < 0.1.

Acknowledgments:

We also thank Marie-Anne Auger-Rozenberg INRAE URZF Orléans, and Thierry Noblecourt, Forest Health Department (DSF) in France, for helping collecting the Sirex populations abroad, especially in native regions in Europe; Natalia Kirichenko (Institute of Forest Siberian Branch of Russian Academy of Sciences, Russia), Emmanuelle Magnoux, Béatrice Courtial and Claudine Courtin (URZF, INRAE, France) for assistance in the lab of URZF; Marie-Anne Auger-Rozenberg and Audrey Bras (URZF, INRAE, France) for advice on data analysis. We thank all workers of the forestry departments in Heilongjiang, Liaoning, Jilin provinces and Inner Mongolia Autonomous Region for their assistance with fieldwork. This study was supported based on the Memorandum of Understanding between The China Scholarship Council and Agreenium (a French public establishment for scientific cooperation bringing together the main agricultural and veterinary research organizations and higher education establishments in France).

Funding: This research was funded by Beijing's Science and Technology Planning Project, grant number Z191100008519004.

Author Contributions: Conceptualization, X.S., J.T., A.R. and Y.L.; methodology, X.S., J.T. and A.R.; software, X.S. and J.T.; validation, Y.L.; formal analysis, X.S.; investigation, X.S. and A.R.; data curation, X.S.; writing-Original draft preparation, X.S.; writing-Review and editing, X.S.; visualization, X.S.; supervision, A.R. and Y.L.; project administration, A.R. and Y.L.; funding acquisition, Y.L. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:

The authors declare no conflict of interest.