Xiaodong Liu

Julien Réthoré
email: julienrethore@ec-nantes.fr

Marie-Christine Baietto

Philippe Sainsot

Antonius Adrianus Lubrecht

An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using CT images

Keywords:

X-ray tomography techniques give researchers the full access to material inner structures. With such ample information, employing numerical simulation on real material images becomes more and more common. Material behavior, especially, of heterogeneous materials, e.g. polycrystalline, composite, can be observed at the microscopic scale with the computed tomography (CT) techniques. In this work, an efficient strategy is proposed to carry out simulations on large 3D CT images. A proposed matrix free type finite element based MultiGrid method is applied to improve convergence speed and to reduce memory space requirements. Homogenization techniques are used to obtain specific operators to enhance the convergence of the MultiGrid method when large material property variations are present. Hybrid parallel computing is implemented for memory space and computing time reasons. Free edge effects in a laminate composite with more than 16 billion degrees of freedom and crack opening in a cast iron are studied using the proposed strategy.

Introduction

One of the subjects of material science is to understand material behavior and then to improve material performance to respond to industrial requirements. However, to manufacture materials with such good performance, is not an easy job. Many defects are present in materials, e.g. voids, cracks or broken fibers in composite materials. CT images of such complex materials permit one to evaluate the presence of these defects. Nevertheless, the aim of researchers is to understand how these defects affect material performance. Numerical simulations are often used to analyze material performance. The main idea of this work is to carry out the calculation directly on CT images. This permits one to observe how the material structure affects its performance at the microscopic scale. However, these defects make the simulation expensive or even impractical.

One of the most used methods is to homogenize materials, which means instead of focusing on the local defects, the material is supposed to be homogenized at a macroscopic scale. This technique is called homogenization. [START_REF] Hashin | Analysis of properties of fiber composites with anisotropic constituents[END_REF] explains how to predict the macroscopic property of a fiber composite from its microscopic structure. Özdemir et al. (2008) presents the application of a homogenization method for the thermal conduction of heterogeneous materials.

However, despite many investigations in the homogenization, the drawback of this technique is that it can not represent microscopic details. One loses physical phenomena related to the local scale. Thus, CT simulations at the microscopic scale become essential to overcome these drawbacks of homogenization methods.

Finite element methods (FEM) are often used to carry out CT simulations. [START_REF] Lengsfeld | Comparison of geometrybased and ct voxel-based finite element modelling and experimental validation[END_REF] and [START_REF] Bessho | Prediction of strength and strain of the proximal femur by a ct-based finite element method[END_REF] present an early and a recent CT images application on the human bone problem by FEM, e.g. hip fractures of the human femur. [START_REF] Ferrant | 3d image matching using a finite element based elastic deformation model[END_REF], [START_REF] Michailidis | An image-based reconstruction of the 3d geometry of an al open-cell foam and fem modeling of the material response[END_REF] and [START_REF] Proudhon | 3d simulation of short fatigue crack propagation by finite element crystal plasticity and remeshing[END_REF] analyze properties of industrial materials by applying FEM simulations on CT images. Another well-known method is the fast fourier transform (FFT). One of the most typical examples is presented in [START_REF] Moulinec | A fft-based numerical method for computing the mechanical properties of composites from images of their microstructures[END_REF], which introduces a FFT-based numerical method to simulate the mechanical properties of composite materials. The finite difference method (FDM) is also often used to carry out CT simulations. [START_REF] Gu | An efficient multigrid solver for the 3d simulation of composite materials[END_REF] presents the application of the FDM with the MultiGrid method on simulations of 3D composite materials. All these three well-known methods perform well for certain cases. However, for FEM applications on materials with complex structures, e.g irregular grains in polycrystalline materials, the meshing step requires human interventions, which is very time consuming. The work in [START_REF] Nguyen | Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials[END_REF] presents the application of FEM on crack propagation for polycrystalline materials, the meshing step for this kind of structure is very complex. For FFT-based numerical methods, the limitation is that it is not robust for no-periodic problems. [START_REF] Gu | An efficient multigrid solver for the 3d simulation of composite materials[END_REF] is a good application of a FDM, however, it was developed only for small problems, furthermore, the implementation of the boundary condition in FDM is intricate.

Motivated by the limitations of the current methods, we propose an efficient and automatic strategy to perform numerical simulations directly from large scale 3D tomography images. A proposed matrix free type finite element method is used in this work. The idea of matrix free-FEM (MF-FEM) was first introduced in [START_REF] Carey | Element-by-element linear and nonlinear solution schemes[END_REF] in 1980s. It was investigated by [START_REF] Augarde | An element-based displacement preconditioner for linear elasticity problems[END_REF]; [START_REF] Hughes | Large-scale vectorized implicit calculations in solid mechanics on a cray x-mp/48 utilizing ebe preconditioned conjugate gradients[END_REF][START_REF] Van Rietbergen | Computational strategies for iterative solutions of large fem applications employing voxel data[END_REF] during the last three decades. Nowadays, it is often used in parallel computing, see [START_REF] Tezduyar | Parallel finite-element computation of 3d flows[END_REF].The MF-FEM was proposed to save memory space and to apply parallel computing. It is particularly favorable for problems with regular and few element types. It does a good job for problems arising from voxel conversions, where the mesh generation can be automatic. In this case, all the generated elements have exactly the same geometry and orientation. [START_REF] Arbenz | A scalable multi-level preconditioner for matrix-free µ-finite element analysis of human bone structures[END_REF] presents a human bone structures simulation using CT images by the MF-FEM. The main drawback of this technique is that it is difficult to find an efficient preconditioner. An often used preconditioner is the element-by-element (EBE) preconditioning [START_REF] Hughes | Large-scale vectorized implicit calculations in solid mechanics on a cray x-mp/48 utilizing ebe preconditioned conjugate gradients[END_REF]) or the Jacobi preconditioning (diagonal scaling). In this work, a modification is proposed to the typical MF-FEM. We propose to not use any matrix but to compute the pre-conditioner degree of freedom (DoF) by DoF (see in Section 2).

However, the MF-FEM does not sufficiently efficient for complex problems like CT simulations. An often used accelerator is the well-known MultiGrid algorithms, as presented in [START_REF] Arbenz | A scalable multi-level preconditioner for matrix-free µ-finite element analysis of human bone structures[END_REF]. The MF-FEM coupled with the aggregation-based MultiGrid method. However, their objective was to perform numerical simulations on complex structure, i.e. human bone, with homogeneous material. In this work, the material heterogeneity is accounted for. In recent years, [START_REF] Kronbichler | Multigrid for matrix-free high-order finite element computations on graphics processors[END_REF] studied the performance of MultiGrid methods for MF high-order FEM on GPUs.

The MultiGrid (MG) method enhances convergence of the FEM and the FDM, e.g. see [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF] for the efficiency of MG algorithms applied on the FEM. MG methods were first introduced in the 1960s, and then, enriched by [START_REF] Brandt | Multi-level adaptive technique (mlat) for fast numerical solution to boundary value problems[END_REF][START_REF] Brandt | Multi-level adaptive solutions to boundary-value problems[END_REF] in the 1970s. It is often used with the FDM. Work of [START_REF] Venner | Multi-level methods in lubrication[END_REF] developed an efficient and robust FDM based MG solver for the elasto-hydrodynamic lubrication problem. The work presented in [START_REF] Watremetz | 2d thermo-mechanical contact simulations in a functionally graded material: A multigrid-based approach[END_REF] introduced the application of MG methods for a 2D graded problem. [START_REF] Boffy | An efficient 3d model of heterogeneous materials for elastic contact applications using multigrid methods[END_REF]; [START_REF] Boffy | Multigrid solution of the 3d stress field in strongly heterogeneous materials[END_REF] employed MG FDM for small heterogeneous material mechanical problems. The work of [START_REF] Gu | An efficient multigrid solver for the 3d simulation of composite materials[END_REF] investigates the FDM based MG methods for the application of composite material properties. This work applies the MG algorithm to achieve a good convergence speed.

The aim of this work is to use a tomographic image of a real material as input, and to carry out a numerical simulation with linear elastic equations on this image to observe the material behavior under a certain loading. Due to the large size of the problem, parallel computing is used. However, the detailed analysis of the performance, in terms of memory and computational cost, is beyond the scope of the present paper which focuses on the ability of the proposed methodology to deal with real 3D images of heterogeneous materials.

The outline of this paper is: Section 2 reviews the fundamentals of linear elasticity and its discretization by the proposed MF-FEM; Section 3 introduces the strategy in detail, both the MG method with special intergrid operators and parallel computing are applied to improve efficiency; Section 4 presents the validation of proposed method and some applications; Section 5 provides discussions and conclusions.

Problem statement and theory

To solve a mechanical problem, the equations of equilibrium are the basics. The fundamentals of solid mechanic equations and their numerical solutions are presented in this section.

Governing Equations

Assuming a deformed domain Ω as presented in Figure 1, its boundary is denoted by ∂Ω. Its exterior normal is n. The second order stress tensor is described by σ. The body force in Ω, which is due to gravity, magnetism, etc., is f . The acceleration is presented by a. The density is denoted by ρ. The equations of equilibrium can be described as:

∇ • σ + ρ f = ρ a (1)
For the sake of simplicity, inertia effects are neglected in this paper. The equation of equilibrium can be simplified as:

∇ • σ = 0 (2)
Besides the equilibrium equations, the boundary conditions on ∂Ω are essential to solve a mechanical problem. Two often used boundary conditions are Dirichlet and Neumann boundary conditions. The Dirichlet boundary condition applied on ∂Ω D is often presented as prescribed displacements, e.g. u = u 0 on Γ 1 . The Neumann boundary condition applied on ∂Ω N refers to the external force f ex , which can be described as:

f ex = σ • n e.g. applying f 0 on Γ 2 .
To solve equation (2), as stated at the beginning of this paper, one has many choices e.g. the FDM, the FEM, the FFT, etc.. The aim of this work is to build an automatic solver using information directly from CT images. As presented in the above section, the matrix free finite element method (MF-FEM) is suitable for the CT simulations. Instead of taking the heavy meshing step, each voxel in the CT image is supposed to be an elementary node. 8-node cuboid Multiplying equation ( 2) by a test function u * and integrating on Ω, the weak form of equation ( 2) reads:

Ω ∇ • σ u * dΩ = 0 (3)
Applying the divergence theorem and integrating by parts, it reads:

Ω σ∇ u * dΩ = ∂Ω σ • n u * dS (4)
which can be described as the equilibrium of internal force f in and external force f ex :

f in = f ex (5) with f in = Ω σ∇ u * dΩ f ex = ∂Ω σ • n u * dS (6)
Applying finite element discretization, the displacements u can be described by:

u c ≈ N i=1 ψ c i ûc i ( 7
)
where i denotes the node number, and N is the number of nodes. The index c represents three directions with c = {1, 2, 3}. The displacement at the elementary node i is û. The component c of the displacement at the node i is described by ûc i . The shape function ψ is the typical shape function of a 8-node cubic element. The test function u * reads:

u * = ψ c j
where j is the node number.

The internal force for node j and component c reads:

f c in = m e 8 g=1 w g σ cm ( u)∇ m ψ c j (8)
with m = {1, 2, 3}. The Gauss integration point number is g, and w g is its weight. 8 Gauss integration points are used in each element, and e is the element number. The relationship between σ and u is implied by σ( u), which can be obtained by the constitutive law (The constitutive law of a linear elastic material is presented in the next part).

To solve equation ( 4), the typical FEM process is to compute the stiffness matrix and to use direct or iterative solvers. However, a typical FEM process is almost impossible for a large scale CT image simulation because of memory space limitations, e.g. for a problem with more than 18 billion DoF. In this work, a modified MF-FEM is used.

The often used matrix free methods are: the power method, the conjugate gradient method, the locally optimal block preconditioned conjugate gradient method and etc.. However, in these methods, small local matrix are presented. We propose a new matrix free type FEM to not reserve any local matrix, but to compute the unknowns node-by-node (DoF-by-Dof) as presented in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF] for the thermal conduction problem.

In this work, it can be described by equation ( 9).

ˆite+1 u = ˆite u + ωS( f ex -f in ) (9)
where ˆite+1 u and ˆite u are the displacements of one node in the current iteration and previous iteration. The relaxation coefficient is denoted by ω, for 0 < ω < 1, it refers to damped Jacobi, i.e. under-relaxation, which is used in this work. The diagonal matrix S is the compliance of the diagonal value of the stiffness matrix at this node. It reads:

S =   1 Ku 0 0 0 1 Kv 0 0 0 1 Kw   (10)
where K u , K v and K w are the diagonal values of the stiffness matrix at each DoF. To obtain S, the stiffness matrix is not necessary, only its diagonnal terms are used.

With the proposed MF-FEM iterative solver, the large scale problem can be solved node by node.

The linear elasticity

In this work, the above strategy is applied to a linear elastic problem.

For a linear elastic problem, the constitutive equation reads:

σ = Ktr(ε)I + 2G ε - 1 3 tr(ε) (11)
where ε is the second order strain tensor, K and G refer to the bulk modulus and the shear modulus, respectively. The unit tensor is denoted by I. Using with the Young's modulus E and the Poisson ratio υ, K and G read:

K = E 3(1 -2υ) , G = E 2(1 + υ)
Strain ε can be described by displacement u:

ε = 1 2 ∇ u + ∇ T u (12)
Its trace being defined as

tr(ε) = tr 1 2 ∇ u + ∇ T u = tr(∇ u)
Combining equation ( 11) and ( 12), equation ( 11) reads:

σ = K - 2G 3 tr(∇ u)I + G(∇ u + ∇ T u) (13)
Applying the finite element discretization with equation ( 7), the stress reads:

σ cm ( u) = i K - 2G 3 k ∇ k ψ k i ûk i δ cm + G (∇ m ψ c i ûc i + ∇ c ψ m i ûm i ) (14)
where k = {1, 2, 3}. δ cm is the Kronecker delta. It reads:

δ cm = 1 c = m 0 c = m
Combining equation ( 8) and ( 14), one obtains the component c of the internal force on node j:

f c in = e 8 g=1 i w g K g - 2G g 3 k ∇ k ψ k i ûk i ∇ c ψ c j + m G g (∇ m ψ c i ûc i + ∇ c ψ m i ûm i ) ∇ m ψ c j ( 15
)
where K g and G g are the material property at a Gauss point. As mentioned before, the material property is assigned to each node from voxel information. An interpolation is used to obtain the material property at each Gauss point. It reads:

K g = 8 i=1 ψ i K i G g = 8 i=1 ψ i G i .
where K u is the first diagonal term of S in equation ( 10). It reads:

K u = e g w g K g - 2G g 3 ∇ c ψ c j ∇ c ψ c j + m G g ∇ m ψ c j ∇ m ψ c j + G g ∇ c ψ c j ∇ c ψ c j ( 16
)
where c = 1. K v and K w are computed with c = 2 and c = 3, respectively. Equation ( 9) can thus be solved DoF by DoF.

MultiGrid

As presented in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF], using only the MF-FEM Jacobi solver, the convergence rate can be very slow after several iterations. The MG method is therefore applied to the iterative solver to improve convergence. The principle of the MG method is that relaxations on the fine grid can eliminate high frequency errors, the low frequency errors can be eliminated by relaxing on the coarse grid. The first step of MG is therefore to construct several levels of grid. In this work, the grid size on level l is two times larger than that on level l + 1, e.g, for a 128 3 -grid problem with 3 levels, the number of points on level is: 128 3 on level 3, 64 3 on level 2 and 32 3 on level 1.

Once different grids are constructed, three important operators are needed to employ a MG scheme: the coarse grid operator, the restriction operator R and the prolongation operator P . As presented in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF], the standard MG scheme can not deal efficiently with problems with large variations of the material properties. Special intergrid operators can be generated as proposed in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF].

The process to obtain the F in and the S indicates that for a representative coarse grid operator, one needs the material property and the shape function on the coarse grid. The shape function of an element on each coarse grid is simple to obtain with the finite element theorem. For the material property on the coarse grid, the standard method is to compute the arithmetic average, i.e Voigt approximation, of the material property on the coarse grid. However, for heterogeneous materials with large property variations, a simple average can not represent the material on the fine grid. Many researchers have tried to apply MG algorithms for problems with large material property variations. The first work to improve MG convergence for problems with large variations is the work of Alcouffe et al. [START_REF] Alcouffe | The multi-grid method for the diffusion equation with strongly discontinuous coefficients[END_REF] for a 2D problem. Hoekema et al. [START_REF] Hoekema | Multigrid solution of the potential field in modeling electrical nerve stimulation[END_REF] extended it to 3D. This work proposes to apply the homogenization method proposed in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF], which applies the homogenization method on MG algorithms to obtain the coarse grid operator for a thermal conduction problem.

The idea is to first compute the Voigt and Reuss approximations of K and G on all the coarse grids recursively as follow:

   Voigt: K HV = N h 1 K h N h , G HV = N h 1 G h N h Reuss: K HR = N h N h 1 1 K h , G HR = N h N h 1 1 G h
where K HV , G HV , K HR and G HR are Voigt and Reuss approximations of the bulk modulus and shear modulus on the coarse grid, The number of nodes on level l is described by N h , which has the same volume as one element on level l -1.

According to the work on homogenization techniques, the effective material property lies between the Voigt and Reuss approximations, which is called Voigt-Reuss (VR) bounds. VR bounds are not the most accurate bounds, but one can compute them recursively from the finest level. The effective material property on the coarse grid can be obtained by the following equation:

     K H = K HV +K HR 2 + √ K HV K HR 2 G H = G HV +G HR 2 + √ G HV G HR 2 (17)
which means the effective K H and G H on the coarse grid is the mean of the arithmetic and geometric average of the VR bounds on this level. According to the comparison of different homogenization methods in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF], the proposed homogenization scheme is robust and efficient.

With the correct material property and correct shape function on the coarse grid, one can finally define the coarse grid operator. For the restriction operator R and the prolongation operator P , the process of the prolongation of corrections can be described with one element on the coarse grid l -1 and its eight elements on the fine grid l as presented in Figure 2. The subject is to bring the corrections, i.e displacement correction e, from level l -1 to level l. The displacement corrections on level l -1 at black points is known. Instead of computing a simple average, one proposes to account for the material property at each node. The bulk modulus K is used in this strategy.

For the black points on level l, one performs an injection from level l -1 to l, e.g. point A1:

e l A1 = e l-1 A1 . O A1 A2 B1 B2 B3 B4 C1 C2 C3 C6 C5 C4 Figure 2: Prolongation scheme
For red points, one computes them from the e at black point of level l, e.g. point B1:

e l B1 = K A1 e l A1 + K A2 e l A2 K A1 + K A2 .
For blue points, one computes them from the e at red points on level l, e.g. point C1:

e l C1 = K B1 e l B1 + K B2 e l B2 + K B3 e l B3 + K B4 e l B4 K B1 + K B2 + K B3 + K B4 .
For yellow points, one computes them from the e at blue points on level l, e.g. point O:

e l O = 1 K C1 + K C2 + K C3 + K C4 + K C5 + K C6 • K C1 e l C1 + K C2 e l C2 + K C3 e l C3 + K C4 e l C4 + K C5 e l C5 + K C6 e l C6 .
For the restriction process, one can use R = P T to obtain the restriction operator.

After constructing several levels and its operators to go up and go down, one has to choose a MG cycle. The best-known MG cycles are V-cycles, W-cycles and FMG cycles. In this work, FMG cycles are performed to have a good initial solution on the finest grid. FMG cycles with 3 levels are illustrated in Figure 3. On level 1, one performs ν 0 relaxations, ν 1 relaxations are performed on each level going up. ν 2 relaxations are performed on each level going down. n cy V-Cycles are used on each level. A linear interpolation of the solution of level l is applied to obtain the initial solution of level l + 1.

With the above strategy, the construction of a matrix free FEM iterative solver based homogenized MG algorithms is finished. One can start performing numerical simulations. However, for a large image, e.g. with more than 18 billion DoF, the computational time is too long using a standard computer. Hybrid MPI-OpenMP programming is therefore applied to achieve good parallel computational performance. For more information about the process and efficiency of parallel computing on the proposed strategy, see [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF].

Validation and applications

In this section, the proposed strategy is validated using a spherical inclusion case. The effective elastic modulus of a spherical inclusion is computed both by the analytical homogenization

Solution validation

To validate the proposed strategy, the typical spherical inclusion case is used. The cubic domain Ω is filled with a spherical inhomogeneity and a homogeneous matrix as presented in figure 4.

The sphere radius is a quarter of the cube size L. The elastic modulus of the sphere is C s , and C m denotes the elastic modulus of matrix. The ratio between these two materials is r e = Em Es , where E m and E s are the Young modulus of matrix and inclusion e.g. Figure 4 presents the case of E m = 500 GPa and E s = 1 GPa. The Poisson ratio equals to 0.3 for both materials. The objective is to compute the effective modulus C of the domain Ω for different r e . For the sake of simplicity, E s equals to 1 GPa.

Figure 4: Young modulus of inclusion and matrix

Eshelby seems to have already solve the spherical inclusion. However, Eshelby analyzed it with an infinite matrix with is not for this case, and it is not possible to have an infinite domain. One of the best-know methods is the Mori-Tanaka (MT) homogenization method, see [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] and [START_REF] Benveniste | A new approach to the application of mori-tanaka's theory in composite materials[END_REF] for details. The effective modulus obtained by MT is:

C M T = C m + V s (C s -C m )A where A = T [(V m -V s )I + V s T ] -1 with T = [I + S s C -1 m (C s -C m )]
-1 where I is the unit tensor, S s is the Eshelby tensor (see [START_REF] Gu | An efficient multigrid solver for the 3d simulation of composite materials[END_REF] for details and Appendix 5 for its formulations), V s denotes the volume fraction of the sphere part, C M T is the effective elastic modulus obtained by MT homogenization.

For the computational homogenization with the proposed MG strategy, the homogeneous displacement boundary condition is applied to obtain the effective modulus C M G . Domain Ω is discretized by 256×256×256 elements. The under-relaxation coefficient equals to 0.5 for the Jacobi solver. The coarsest grid has 4×4×4 elements, thus there are 7 levels. The FMG cycles use n cy = 6, ν 0 = 50, ν 1 = 4 and ν 2 = 2. The comparison between C M T and C M G is qualified using the following equation:

e M T -M G = 1 n t nt i=1 C M T i -C M G i C M T i 2 where C M T i is each term of C M T . C M G i is each term of C M G .
n t is the number of components in elastic modulus, which equals to 36.

The elastic modulus for r e = 10 obtained by the MT and the MG is presented as the following: From these two elastic modulus, we find that the effective material is isotropic, which confirms to this case. Table 1 presents the first item of the effective modulus obtained by the MT and the MG method for different r e . e M T -M G is also presented in Table 1. From this table, a good agreement can be found between the MT and the MG. There is a small difference (less than 0.5%) between these two methods. The proposed method is therefore validated. The convergence check is also carried out by computing the error L 2 norm on level l -1 using the converged solution on l during FMG cycles. Figure 5 presents the error L 2 norm in function of the element size for the spherical case r e = 10. It demonstrates the linear relationship between the error L 2 norm and the element size h in logarithmic scale. This is as expected for FE based multi-grid. The accuracy increases by a factor of 3 (4 for homogeneous materials) when the element size decreases by a factor of 2, which confirms the ability of the proposed multi-grid strategy to converge even for the large material property variation cases. where x is the coordinate of the X direction. Some strain raisers are found at the interface of matrix and inhomogeneity. The strain inside the sphere is smaller than in the matrix which is typical for the soft inclusions.

C M T =        
x y z

Efficiency analysis

The efficiency of the proposed strategy is also analyzed for this mechanical problem. The principle is to perform the numerical simulation for r e = 10 of the spherical inclusion case with both the FMG scheme and the single level Jacobi iterative solver. The total cost of FMG cycles and the V-Cycle is described by the following equation ( 19) according to [START_REF] Venner | Multi-level methods in lubrication[END_REF].

W F M G ≤ n cy 1 -(H/h) -d W V -Cycle ≤ n cy (ν 1 + ν 2 ) (1 -(H/h) -d ) 2 W U (19)
where W U is the cost of one relaxation on the finest grid. W F M G and W V -Cycle are the cost of FMG cycles and the V-Cycle. H/h = 2 in this work which presents the ratio between the grid size of coarse l -1 and fine grid l. d is the problem dimension which is d = 3 in this work. There is also the cost of transform between grids, but problem solved in this work are too large, the relaxation time is much more expensive than others routines. One accounts only the cost of relaxations on all grids.

Figure 4.2 implies that with 2000 relaxations, i.e. a cost of 2000 W U , the single level Jacobi iterative solver still can not achieve the initial residual of the FMG scheme. The FMG achieves a residual of 10 -5 in with 36 relaxations on the finest grid, i.e. a cost of 6(4+2) (1-2 -3 ) 2 W U ≈ 47W U . The FEM convergence slows down after only few relaxations. However, the FMG scheme remains its good convergence. The red line in Figure 7(b) demonstrates the residual at the finest level after each V-Cycle of the FMG scheme, where a good convergence can be found.

As mentioned above, the proposed method uses only the diagonal values of the stiffness matrix. For a N nodes mechanical problem ( discretized using 8-node cubic elements), the memory requirement for the sparse global stiffness matrix is bounded by (3 × 27) × (3 × N ) × 8 bytes (assuming 8 bytes floatting point numbers) which gives 1944N bytes. This estimate assumes that each node has 27 neighbors. However, there are also surface and edge nodes for which the number of neighbors is lower which makes the proposed estimate an upper bound of the actual memory cost. For the proposed method, the memory requirement for the diagonal of the stiffness matrix is 3N × 8 = 24N bytes, which is 81 times cheaper than the global stiffness matrix. For a l levels MG algorithm, the total size is:

24N × 1 + 1 8 + 1 8 2 + • • • + 1 8 l < 27.43N bytes.
It demonstrates the memory requirement of the global sparse stiffness matrix is about 71 times more expensive than the proposed method. e.g. The largest problem in this work has more than 6 billion nodes, which means the size of the global sparse stiffness matrix is more than 10 TB. The proposed method needs only about 0.15 TB.

Laminated material simulations

In this subsection, a simulation on a problem with more than 18 billion DoF is presented. The numerical result is then compared to the experimental result obtained by the digital image correlation (DIC).

The CT image of a laminated material is used in this case. This CT image has 700×1300×1700 voxels with 4.5 µm/pixel. This laminated structure consists of four layers with a fiber orientation of 15 • , -15 • , -15 • and 15 • , respectively. The fiber is the E-glass fiber with a Young's modulus of 80.0 GPa. The matrix is s M9 epoxy with a Young's modulus of 3.2 GPa. The Poisson ratio of these two materials equals to 0.22. Figure 8 illustrates a section view of an interface of layers, one All these defects, which can only be seen by CT, have an impact on the material mechanical behavior under certain loading. More information about this image and this material can be found in [START_REF] Lecomte-Grosbras | Three-dimensional investigation of free-edge effects in laminate composites using x-ray tomography and digital volume correlation[END_REF].

The main subject of this subsection is to analyze free edge effects in the laminated structure by numerical simulation. A qualitative comparison between the numerical results and the experimental results is also presented. The free edge effect was firstly presented by [START_REF] Pipes | Interlaminar stresses in composite laminates under uniform axial extension[END_REF], who found the strain concentrations around free edges and the ply interface. [START_REF] Lecomte-Grosbras | Interlaminar shear strain measurement on angle-ply laminate free edge using digital image correlation[END_REF] illustrates the free edge effects by the DIC experimental method of an unidirectional carbon fiber reinforced plastic laminated structure. For more information about free edge effects, see [START_REF] Pipes | Interlaminar stresses in composite laminates under uniform axial extension[END_REF]. In this paper, instead of carrying out DIC experiments, the CT image of the laminated material is used directly to employ numerical simulations and understand the free edge effect. To carry out numerical simulation, one proposes to take a part of this image which refers to the region of interest (ROI). The ROI is constituted of 577×1153×1153 voxels as presented in Figure 10. To have more voxels in the fiber section, one performs a sub-sampling, i.e. linear interpolation, on this ROI. Figure 11 shows that we have about ten pixels per fiber diameter. The final input image, i.e. domain Ω, has 1153×2305×2305 voxels which means we have more than 6 billion elements, i.e. 18 billion DoF. The coordinates of the center of Ω is (0,0,0) with a size of L×2L×2L.

The material property in each element is indispensable to perform numerical simulations. However, CT images contain the gray level (GL) which varies from 0 to 255. One has to choose a strategy to obtain material properties from CT images. The often used method is the thresholding method. As presented in Figure 12(a), the interface between the E-glass fiber and M9 epoxy matrix is not extraordinarily sharp. It is difficult to distinguish between these two phases (matrix and fiber). From the histogram of this laminated composite in Figure 12(b), the interface is not obvious neither. Instead of applying the thresholding method, one proposes to apply a continuous property.

α = 1 -e -|GL-160.5|

20

• sign (GL -160.5) + 1 where α is the material property variation, e.g. K,G. The parameters in this equation are obtained by the trial and error method to ensure a good agreement between the original image and material property distribution. This equation can be described by Figure 4.3. The boundary conditions are given by the following equation:

u = {0, 0, -0.01L}, on z = -L u = {0, 0, 0.01L}, on z = L (20)
where u z is the displacement in the Z direction.

There are 1153×2305×2305 elementary nodes on the finest level and 9×18×18 elements on the coarsest level with a total of 8 levels. The under-relaxation coefficient equals to 0.5 for the Jacobi solver. The parameters of the FMG cycles are: n cy = 9, ν 0 = 500, ν 1 = 4 and ν 2 = 2.

Figure 13 illustrates the Young modulus and the shear strain XZ on surface y = -L. This figure demonstrates that the shear strain field mimics the Young modulus distribution. Shear strain concentrations are found on the two interfaces. This is the so called free edge effect in laminated materials. In the area with few fibers with too much matrix, they also aches strain concentrations which can lead to preventive damage in industrial applications. The same phenomena are found by the DIC experimental method in [START_REF] Lecomte-Grosbras | Interlaminar shear strain measurement on angle-ply laminate free edge using digital image correlation[END_REF]. The 2009) is a carbon fiber reinforced plastic laminated structure which is similar to the material used in this work. A good correlation between the numerical and the experimental results is found.

Figure 14 reveals the displacement in the Z direction on the two opposite surfaces y = -L and y = L. Equally, the free edge effects can be found on the interfaces, which leads to the displacement variations on the interfaces. This displacement variation can also be found in [START_REF] Lecomte-Grosbras | Interlaminar shear strain measurement on angle-ply laminate free edge using digital image correlation[END_REF], the displacement curves illustrated in Figure 15, have the same tendency both for the numerical results and the DIC results. Another phenomenon that we can observe is that the displacement variations on y = -L is almost anti-symmetric to y = L. This is because the fiber orientations on these two opposite surfaces are anti-symmetric. The Z displacements along the X axis in the center of these two surfaces are illustrated in Figure 15 where a clear anti-symmetry can be found.

Cast iron applications

In this work, we perform the proposed strategy on a nodular graphite cast iron CT image. The image used in this paper is the image obtained by [START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF] using X-Ray tomography. This image consists of 340×340×512 voxels with a 5.06 µm/pixel size. For more information about this CT image, please refer to [START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF]. A ROI with 257×257×257 voxels is taken from this CT image. As illustrated in Figure 16, many carbon nodules with a random distribution, can be found at the microscopic scale of the CT image.

The objective is to perform the linear elastic simulation on the ROI, with a prescribed rectangular crack to see how carbon nodules affect the crack opening. Assuming Ω is the ROI. The center of the Ω is the origin of axis. The size of Ω is L. 256 3 8-node cubic elements are used to discretize Ω. The rectangular crack is presented in Figure 17. The width of the crack is the size of the cube L. Its length is L 3 . The crack is constructed by setting material property as 0 one three layers of nodes in Z direction. The prescribed boundary conditions are given by Here, the thresholding method, i.e. Figure 4.4, is used to obtain the material property due to an obvious boundary between iron and carbon nodular. To confirm this choice, the volume fraction of carbon nodular after the thresholding is equal to the value provided by the manufacturing company.

The material properties are given in Table 2. Table 2 implies that the crack is defined by setting the material property to 0 at the crack nodes. The crack domain is considered as nodes with a material property 0, which means the computational domain is the entire cubic Ω. To perform this numerical simulation, a specific treatment is required to deal with the infinite material property jump. This treatment involves in the following steps.

• Step A: The MF-FEM iterative solver for the relaxation step.

• Step B: The injection of solutions from the fine grid to the coarse grid.

• Step C: The restriction of residuals from the fine grid to the coarse grid.

• Step D: The prolongation of the corrections from the coarse grid to the fine grid. For "Step A": during the relaxation processing, the displacements on nodes without material are not updated.

For "Step B": as mentioned above, the Reuss approximation is

K HR = N h N h 1 1 K h , G HR = N h N h 1 1 G h
To avoid 0 in the denominator, one defines: if there is K h = 0 and G h = 0, then K HR = 0 and G HR = 0. For "Step C" and "D": the restriction and the prolongation process are done only when the material property on this node and on all of its nearest neighbor nodes is not 0 at the fine grid. If 0 appears in the denominator when computing the restriction matrix, one replaces it with 10 -6 .

For the grids, one has 256×256×256 elements on the finest level and 4×4×4 elements on the coarsest level with 7 levels. The under-relaxation coefficient is taken to 0.5 for the Jacobi solver. The parameters of the FMG cycles are: n cy = 4, ν 0 = 100, ν 1 = 8 and ν 2 = 4. Since the displacement and the crack thickness are too small for the visualization, i.e. 1%, in the following figures, the amplitude of mesh deformation is multiplied by 20.. Figure 19 illustrates the Young's modulus and the strain field in the cast iron. It presents the strain concentrations on the crack front. The largest strain can be found in carbon nodules on the crack front. The largest strain is 10 times larger than the prescribed strain at the macroscopic scale, i.e. 1%. Other strain concentrations are located on carbon nodules, because carbon nodules are 10 times softer than iron. Another simulation is carried out to compare the crack opening in a homogeneous material and in a heterogeneous material. The principle is to replace all the carbon nodules in the CT image by iron, which means the simulation is carried with a prescribed crack in the iron. The same boundary conditions are applied on the simulation of the crack opening in iron. Figure 20 shows the strain field in the homogeneous iron and heterogeneous cast iron. The "butterfly" strain field in the iron presents typical strain concentrations in the homogeneous materials after the crack opening. Compared to the strain concentrations in the iron, the strain concentrations in the cast iron are not only in the vicinity of the crack front but also in the carbon nodules over the entire volume. The material heterogeneity spreads strain concentrations in a large volume. The "butterfly" strain field is destroyed by carbon nodules. Figure 21 illustrates the strain concentrations on the crack front in the iron and in the cast iron, respectively. The largest strain is located on the carbon nodules on the crack front in the cast iron, but in iron, a uniform strain concentration can be found on the entire crack front. This should lead to non rectilinear crack front (if crack propagation is performed) and strong interactions between the crack and the material microstructure.

Conclusions

In this paper, several strategies are proposed to achieve a CT simulation with more than 16 billion of degree of freedom. The application of the matrix free finite element method permits us to carry out such large simulation with the minimum memory space. The MultiGrid method improves convergence compared to a single level Jacobi iterative solver. The use of the proposed

Iron Cast iron

Iron

Cast iron homogenization techniques improve the convergence stability of the MultiGrid method. The applications presented in this paper confirms the CT image based simulations permit one to analyze the materials behavior at the microscopic scale. The good correlation between the DIC experiments and the CT simulations for the laminated structures confirms the strategy proposed in this paper. The free edge effect in the laminated materials proves the importance to carry out numerical simulations at the microscopic scale. The analysis of the crack opening application in the cast iron demonstrates how the carbon nodules affect the strain concentrations in the materials. It implies the importance to account for the carbon nodules in the numerical simulation and also in the application of materials in the industrial domain. Since the largest strain concentrations are located in the soft carbon nodules, the damage starts firstly in this areas. Both these applications imply the importance to carry out simulations at the microscopic scale by using the CT images. Performing numerical simulations with real material microstructures permits ones to observe many phenomena that one can never observe at the macroscopic scale or the theoretical microstructure.

Figure 1 :

 1 Figure 1: An elastic body Ω and its boundary

Figure 3

 3 Figure 3: FMG scheme

  Figure 5: Convergence check

Figure 6 :

 6 Figure 6: Strain field of spherical inclusion case.

Figure 7 :

 7 Figure 7: Convergence of the Jacobi solver (a) and FMG scheme (b) on a 256 3 nodes problem

Figure 8 :

 8 Figure 8: Sections of material

Figure 9 :ZFigure 10 :

 910 Figure9: Gray level and material property variations. "0" represents epoxy, "1" represents fibers

Figure 12 :

 12 Figure 12: Laminated composite image and its histogram

Figure 13 :

 13 Figure13: Young modulus and strain on the surface y = -L; the white dashed lines are the interfaces of layers. Strain o is the macroscopic strain, which is 1%.

Figure 14 :Figure 15 :

 1415 Figure 14: Displacement u z on the surface y = -L and on the surface y = L; the while dashed lines are two interfaces of layers, the green line is the z = 0 line which is the center of the surface

Figure 16 :Figure 17 :

 1617 Figure 16: The GL of the ROI of the CT image of cast iron

Figure 18 :

 18 Figure 18: Thresholding method to obtain material property from GL in image

Figure 19 :

 19 Figure 19: The Young's modulus and the strain zz in cast iron. Strain o = 1%, which is the macroscopic strain. The amplitude of mesh deformation is multiplied by 20.

Figure 20 :

 20 Figure 20: Strain zz in iron and in cast iron. Strain o = 1%, which is the macroscopic strain. The amplitude of mesh deformation is multiplied by 20.

Figure 21 :

 21 Figure 21: Strain concentrations on the crack font in the cast iron and in the iron. Strain o = 1%, which is the macroscopic strain. The amplitude of mesh deformation is multiplied by 20.

Table 1 :

 1 C 1111 obtained byMT and MG

	r e	C M T 1111 /GPa C M G 1111 /GPa e M T -M G
	1	1.346	1.346	0
	10 11.934	11.986	0.31%
	50 58.220	58.557	0.40%
	100 116.037	116.730	0.42%
	200 231.6670	233.068	0.43%
	500 578.549	582.073	0.43%

Table 2 :

 2 Material properties in cast iron

	Component	Young's modu-lus/GPa	Poisson ratio
	Iron	210	0.3
	Carbon granules 21	0.2
	Crack	0	0

Acknowledgments

We gratefully thank for the support of Région Pays de la Loire and Nantes Métropole through a grant Connect Talent IDS.

Appendix

The following Eshelby's tensor is for the spherical inclusion of a linear elastic problem in macroscopically isotropic materials.

where ν is the Poisson's ratio of the inclusion.


Nathalie Vega

Claudie Pinteur

Gaël Buffelan

Emmanuelle Loizon

Hubert Vidal

Danielle Naville

Dr Brigitte Le Magueresse-Battistoni

Exposure to pollutants altered glucocorticoid signaling and clock gene expression in female mice. Evidence of tissue-and sex-specificity

Keywords: Pollutant mixture, mouse, Adipose tissue, Glucocorticoid signaling, Clockrelated genes Highlights

tetrachlorodibenzo-p-dioxin (20 pg/kg body weight/day [bw/d]), polychlorobiphenyl 153 (200 ng/kg bw/d), di-[2-ethylhexyl]-phthalate (500µg/kg bw/d) and bisphenol A (40µg/kg bw/d).

In female mice fed a standard diet (ST), we observed a decrease in plasma levels of leptin as well as a reduced expression of corticoid receptors Nr3c1 and Nr3c2, of leptin and of various canonical genes related to the circadian clock machinery in visceral (VAT) but not subcutaneous (SAT) adipose tissue. However, Nr3c1 and Nr3c2 mRNA levels did not change in high-fat-fed females exposed to pollutants. In ST-fed males, pollutants caused the same decrease of Nr3c1 mRNA levels in VAT observed in ST-fed females but levels of Nr3c2 and other clock-related genes found to be down-regulated in female VAT were enhanced in male SAT and not affected in male VAT. The expression of corticoid receptors was not affected in the livers of both sexes in response to pollutants. In summary, exposure to a mixture of pollutants at doses lower than the no-observed adverse effect levels (NoAELs) resulted in sexdependent glucocorticoid signaling disturbances and clock-related gene expression modifications in the adipose tissue of ST-fed mice.

Introduction

Metabolic-related diseases, such as diabetes and obesity, are prominent public health challenges of the 21 st century. In the last decades, prevalence of these chronic diseases reached epidemic proportions worldwide. Metabolic diseases have multifactorial causes.

Genetic predisposition, high calorie intake, poor physical activity and also environmental pollutants are among the etiological factors identified. The capacity of some pollutants to interfere with hormonal action, also known as endocrine disruptors has been well documented [START_REF] Carpenter | Environmental contaminants as risk factors for developing diabetes[END_REF][START_REF] Diamanti-Kandarakis | Endocrine-disrupting chemicals: an Endocrine Society scientific statement[END_REF], Thayer et al., 2012[START_REF] Ibrahim | Chronic consumption of farmed salmon containing persistent organic pollutants causes insulin resistance and obesity in mice[END_REF][START_REF] Ruzzin | Persistent organic pollutant exposure leads to insulin resistance syndrome[END_REF][START_REF] Alonso-Magdalena | Endocrine disruptors in the etiology of type 2 diabetes mellitus[END_REF][START_REF] Heindel | Environmental Obesogens: Mechanisms and Controversies[END_REF], Le Magueresse-Battistoni et al., 2015). Endocrine disruptors have been shown to target the estrogen, androgen, thyroid, and steroidogenesis pathways and others (e.g., through binding to estrogen-related receptor, peroxisome proliferator-activated receptor, xenobiotic receptors) [START_REF] Diamanti-Kandarakis | Endocrine-disrupting chemicals: an Endocrine Society scientific statement[END_REF][START_REF] Reif | Endocrine profiling and prioritization of environmental chemicals using ToxCast data[END_REF][START_REF] Delfosse | A structural perspective on nuclear receptors as targets of environmental compounds[END_REF], depending on the repertoire of co-expressed nuclear receptors and hormonal environment [START_REF] Casals-Casas | Endocrine disruptors: from endocrine to metabolic disruption[END_REF]Desvergne, 2011, Le Magueresse-Battistoni et al., 2017). Therefore, chemicals targeting energy homeostasis could have extensive metabolic implications. Indeed, energy homeostasis is a complex process orchestrated by multiple hormones adjusting energy intake to energy expenditure including sex hormones [START_REF] Clegg | The role of estrogens in control of energy balance and glucose homeostasis[END_REF] and glucocorticoids (Vegiopoulos and Herzig, 2007).

Over the course of a lifetime, humans are exposed to a plethora of chemicals that can interact in vivo by modifying each other's effect. In addition, some chemicals such as bisphenol A (BPA) were shown to cause adverse effects in experimental studies (Vandenberg et al., 2012[START_REF] Angle | Metabolic disruption in male mice due to fetal exposure to low but not high doses of bisphenol A (BPA): evidence for effects on body weight, food intake, adipocytes, leptin, adiponectin, insulin and glucose regulation[END_REF] even at doses lower than those considered to be safe in humans [START_REF] Koch | Human body burdens of chemicals used in plastic manufacture[END_REF]Calafat, 2009, Schafer and[START_REF] Schafer | Persistent toxic chemicals in the US food supply[END_REF]. Therefore, the cocktail effect resulting from exposure to complex mixtures makes it difficult to predict the health consequences (Le Magueresse-Battistoni et al., 2018b[START_REF] Kortenkamp | Low dose mixture effects of endocrine disrupters: implications for risk assessment and epidemiology[END_REF], Heindel et al., 2017b). Interestingly, we have recently implemented a model in which mice were exposed through diet to a mixture of pollutants made of BPA, di-[2-ethylhexyl]-phthalate (DEHP), polychlorinated biphenyl (PCB) 153 and 2,3,7,. These chemicals are part of a large group of archetypal endocrine disruptors present in food [START_REF] Casals-Casas | Endocrine disruptors: from endocrine to metabolic disruption[END_REF], Heindel et al., 2017a, Vandenberg et al., 2012) with different characteristics regarding half-lives, endocrine disrupting activities and activation of nuclear receptors/transcription factors (Le Magueresse-Battistoni et al., 2018a).

BPA is probably the most emblematic chemical since humans are almost continuously exposed to it, with estimated daily doses of 0.4 to 5 µ/kg body weight (Vandenberg et al., 2012). BPA has large industrialized applications, particularly in plastics and resins, and its global production exceeds millions of tons per year (Vandenberg et al., 2007). Initially described as pro-estrogenic because of its binding to estrogen receptors (ERs), the modes of action of BPA may involve non-estrogenic receptors including the peroxisome proliferator activator receptor (PPAR)γ, the glucocorticoid receptor (GR) or the thyroid receptor, all intricately linked to metabolic health and body weight (Le Magueresse-Battistoni et al., 2018a, Zoeller et al., 2005[START_REF] Casals-Casas | Endocrine disruptors: from endocrine to metabolic disruption[END_REF]. Disruptive effects on glucose and lipid metabolism or insulin resistance may also result from mechanisms linked to hormonal metabolism, oxidative stress, inflammation, mitochondrial dysfunction and epigenetic changes (Vandenberg et al., 2012, Vom Saal et al., 2012[START_REF] Ma | The adverse health effects of bisphenol A and related toxicity mechanisms[END_REF], Nadal et al., 2017). Phthalates are other short-lived chemicals that are used primarily to soften plastics.

They are of great concerns for human health with an average exposure of 0.5 to 25 µ/kg body weight per day for DEHP (Vandenberg et al., 2012). As with BPA, the diet is the main source of contamination and results from the chemical leaching of food and drink containers and packaging [START_REF] Koch | Human body burdens of chemicals used in plastic manufacture[END_REF]Calafat, 2009, Vandenberg et al., 2007). Experimental studies have shown that DEHP exposure can cause adverse effects on energy balance and glucose homeostasis through binding to PPARs [START_REF] Feige | The endocrine disruptor monoethyl-hexyl-phthalate is a selective peroxisome proliferatoractivated receptor gamma modulator that promotes adipogenesis[END_REF]. Phthalates can also bind to xenobiotic receptors but not to androgen receptors although they exert anti-androgenic effects [START_REF] Schug | Endocrine disrupting chemicals and disease susceptibility[END_REF][START_REF] Casals-Casas | Endocrine disruptors: from endocrine to metabolic disruption[END_REF], Le Magueresse-Battistoni et al., 2017).

Dioxins and polychlorinated biphenyls (PCBs) are lipophilic and belong to the group of persistent organic pollutants (POP). Epidemiological and experimental studies demonstrated that exposure to POPs could lead to detrimental metabolic effects including obesity and type 2 diabetes (Wahlang et al., 2013[START_REF] Merrill | Toxicological function of adipose tissue: focus on persistent organic pollutants[END_REF][START_REF] Gauthier | The metabolically healthy but obese phenotype is associated with lower plasma levels of persistent organic pollutants as compared to the metabolically abnormal obese phenotype[END_REF][START_REF] Ibrahim | Chronic consumption of farmed salmon containing persistent organic pollutants causes insulin resistance and obesity in mice[END_REF], van Esterik et al., 2015). Dioxins result from incomplete combustion, which occurs during industrial processes, forest fires or volcanic eruptions. PCBs had a large number of industrial applications because of their low-flammability and high-thermal conductivity properties. Today, the production of PCBs is banned and industrial emissions of dioxins are limited. However, they still contaminate food by bio-amplification in the food chain and are concentrated in fatty foods due to their low degradability [START_REF] Guo | Persistent Organic Pollutants in Food: Contamination Sources, Health Effects and Detection Methods[END_REF]. It was reported that average daily exposure to TCDD (the most toxic dioxin congener that includes the coplanar PCBs) was below the tolerable daily intake (TDI) dose but it was up to 80 ng/kg body weight per day for non-coplanar PCBs thus exceeding the toxicological reference values for certain groups of consumers [START_REF] Baars | Dioxins, dioxin-like PCBs and non-dioxin-like PCBs in foodstuffs: occurrence and dietary intake in The Netherlands[END_REF], EFSA, 2005). Interestingly, PCB153 a non-coplanar PCB is considered as a marker for the total PCB exposure [START_REF] Hagmar | Intraindividual variations and time trends 1991-2001 in human serum levels of PCB, DDE and hexachlorobenzene[END_REF].

TCDD and co-planar PCBs can bind to the aryl hydrocarbon receptor (AHR) [START_REF] Barouki | The aryl hydrocarbon receptor system[END_REF][START_REF] Bock | Aryl hydrocarbon receptor (AHR): From selected human target genes and crosstalk with transcription factors to multiple AHR functions[END_REF], Van den Berg et al., 1998) while non-coplanar PCBs cannot bind to AHR although they can interact with its signaling mechanisms. They may as well activate the pregnane X receptor (PXR) and the constitutive androstane receptor (CAR) [START_REF] Chen | Non-coplanar and coplanar polychlorinated biphenyls potentiate genotoxicity of aflatoxin B1 in a human hepatocyte line by enhancing CYP1A2 and CYP3A4 expression[END_REF]. Dioxins and PCBs can also interact with the estrogen and thyroid signaling pathways, respectively [START_REF] Diamanti-Kandarakis | Endocrine-disrupting chemicals: an Endocrine Society scientific statement[END_REF][START_REF] Casals-Casas | Endocrine disruptors: from endocrine to metabolic disruption[END_REF][START_REF] Pavek | Pregnane X Receptor (PXR)-Mediated Gene Repression and Cross-Talk of PXR with Other Nuclear Receptors via Coactivator Interactions[END_REF].

Using the exposure model to a mixture of BPA, DEHP, PCB153 and TCDD, we previously demonstrated that the metabolic adverse impact in mice was strongly sex-dependent altering estrogen signaling in the liver of exposed females and hepatic cholesterol metabolism in males [START_REF] Naville | Low-dose food contaminants trigger sex-specific, hepatic metabolic changes in the progeny of obese mice[END_REF]. In addition, the mixture of pollutants affected other metabolic tissues including the adipose tissues and involved other nuclear receptors than the estrogen and xenobiotic receptors [START_REF] Julien | Evidence for estrogeno-mimetic effects of a mixture of low-dose pollutants in a model of ovariectomized mice[END_REF], Julien et al., 2019, Naville et al., 2015, Naville et al., 2019). The purpose of the current study was to address the implication of glucocorticoid signaling in mediating the effects of the pollutants in female mice. We also looked for sexual dimorphism benefiting from male samples taken in a previous study [START_REF] Naville | Chronic exposure to a pollutant mixture at low doses led to tissue-specific metabolic alterations in male mice fed standard and high-fat highsucrose diet[END_REF].

Glucocorticoids (GC) are steroid hormones with strong anti-inflammatory properties. They are produced by the adrenal glands under the negative feed-back involving the hypothalamicpituitary-adrenal (HPA) axis. They exert vital roles on behavior, immunity and metabolism through binding to corticoid receptors in target cells including the liver and adipocytes (Vegiopoulos and Herzig, 2007). Recent reports highlighted that xenobiotics could bind to GR [START_REF] Gulliver | Xenobiotics and the Glucocorticoid Receptor[END_REF], and that BPA and phthalates had a weak GR agonist activity in adipocytes [START_REF] Gulliver | Xenobiotics and the Glucocorticoid Receptor[END_REF], Le Magueresse-Battistoni et al., 2017[START_REF] Casals-Casas | Endocrine disruptors: from endocrine to metabolic disruption[END_REF][START_REF] Sargis | Environmental endocrine disruptors promote adipogenesis in the 3T3-L1 cell line through glucocorticoid receptor activation[END_REF].

Materials and Methods

2.1.Animals, diet and experimental design

All procedures were performed with the approval of the Regional Committee of Ethics for Animal Experiments. Four-week (wk) old C57Bl/6J female and male mice (Envigo; Gannat, France) were randomly housed two per polypropylene cage at 21°C with normal light/dark cycle and free access to water and standard chow diet (from Genestil, Royaucourt, France) for a one-week acclimatization period. Mice were next randomly assigned to four different groups (n=7-8 / group). Two groups were fed a standard diet (ST; 10% Fat and 35.6% maltodextrin) and two groups were fed a high-fat high-sucrose diet (HF; 39.4% Fat and 16.6% maltodextrin +16.6% sucrose, Envigo). The STp and HFp groups correspond to the ST or HF diet, containing respectively a mixture of 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD, CAS N° 1746-01-6 from LGC Promochem, Molschein, France), polychlorinated biphenyl (PCB) 153 (CAS N° 35065-27-1), bisphenol A (BPA, CAS N° 80-05-7) and di-[2ethylhexyl]-phthalate (DEHP, CAS N°117-81-7) (all three from Sigma-Aldrich). The ST0 and HF0 groups correspond to the ST or HF diet containing the vehicle, only. To prepare the contaminated food, stock solutions of pollutants dissolved in DMSO were diluted in corn oil and they were mixed into diet (316 µl DMSO in 60 ml corn oil/kg diet) using a foodprocessor to ensure good distribution in the food. Specifically, mice were exposed to 20 pg TCDD /kg bw/day, 200 ng PCB153/kg bw/day, 500µg DEHP/kg bw/day and 40µg BPA/kg bw/day (Supplementary Table 1). These doses were at least 10 times lower than the NoAEL (No-observed Adverse Effect Level) dose for each pollutant when taken individually (DEHP, BPA) or representative congener (TCDD, PCB153) [START_REF] Beausoleil | Regulatory identification of BPA as an endocrine disruptor: Context and methodology[END_REF], van Leeuwen et al., 2000, WHO, 2003). To ensure that animals ingested the correct amount of polluted food/kg bw/day and were fed ad libitum, the calculations of intake of pollutants in the diet were based on 1 g of contaminated food /17 g of body weight / day, and pollutant-free food was provided to the animals at libitum. Food was given on Mondays, Wednesdays and Fridays. Body weight was recorded weekly and the amount of polluted food distributed was adjusted, accordingly over a 15 weeks period. It should be noted that the male samples used in the present study have previously been analyzed for lipid metabolism in the liver and jejunum, and adipogenesis and inflammation of the fat pads [START_REF] Naville | Chronic exposure to a pollutant mixture at low doses led to tissue-specific metabolic alterations in male mice fed standard and high-fat highsucrose diet[END_REF], but not glucocorticoid signaling.

2.2.Tissue collection and blood and plasma analyses

Before sampling, mice were fasted six-hour and weighed. Blood was collected by retro-orbital sampling and mice were euthanized by cervical dislocation. Adipose tissues and liver were quickly removed, weighed, frozen in liquid nitrogen and stored at -80°C. Blood glucose concentrations (Accucheck Performa glucometer; Roche Diabetes Care France, Meylan), and plasma levels of insulin and corticosterone (Mouse Ultrasensitive ELISA, Eurobio, Courtaboeuf, France), leptin and adiponectin (ELISA kits from Crystal Chem Europe, Zaandam, Netherlands), triglycerides (Biolabo, Maizy, France) and free fatty acids, total cholesterol and cholesteryl esters (all kits from Sigma-Aldrich) were measured.

RNA extraction and real-time PCR analyses

Total RNA was extracted from subcutaneous (SAT) and visceral (VAT) adipose tissue using TRI Reagent (Applied Biosystems, Courtaboeuf, France). RNA integrity was determined with the Agilent 2100 Bioanalyzer and RNA 6000 Nano Kit (Agilent Technologies, Massy, France). Samples were reverse-transcribed using the Prime Script RT Reagent kit (Takara Bio Europe SAS, Saint-Germain-en-Laye, France). We performed real time PCR using a set of specific primers specifically designed to encompass an intron (Supplemental Table 2). Data were normalized using Tbp encoding TATA-Box Binding Protein as described for adipose tissues and beta-glucuronidase (Gusb) for liver samples [START_REF] Naville | Direct and indirect impact of 2,3,7,8tetrachlorodibenzo-p-dioxin (TCDD) on adult mouse Leydig cells: an in vitro study[END_REF], Naville et al., 2019).

2.4.Western blotting analyses

Proteins prepared from adipose tissues were separated by SDS-10% polyacrylamide gel electrophoresis and transferred to a polyvinylidene difluoride membrane. Immunoblotting was performed as previously described [START_REF] Naville | Low-dose food contaminants trigger sex-specific, hepatic metabolic changes in the progeny of obese mice[END_REF][START_REF] Longin | MT1-MMP in rat testicular development and the control of Sertoli cell proMMP-2 activation[END_REF] using a rabbit polyclonal anti-glucocorticoid receptor (dilution 1/5000, Ab3580; Abcam France) or β-Actin (dilution1/10,000, A5060, Sigma). Blots were revealed using Luminata Classico Western HRP substrate (Milllipore, Molsheim, France), detected using the ChemiDoc XRS+ imaging system (Biorad), and analyzed with Image Lab software (Bio-Rad). Data were normalized relative to β-Actin.

2.5.Statistical analyses

Results are presented as means ± S.E.M. Graphpad Prism 8.01 software was used for statistical analyses. The overall impact of the pollutant exposure and of the diet was determined using the two-way Anova analysis followed by the Bonferroni's multiple comparisons post-hoc test to determine significant differences between relevant groups. In experiments only assessing the impact of pollutant exposure, differences were tested using the non-parametric Mann-Whitney test. Differences between means were considered significant at p-value < 0.05.

Results

Metabolic phenotyping of female mice

We explored the metabolic alterations of female mice in response to pollutant exposure including the analysis of the systemic levels of corticosterone. No changes in body weight, blood glucose, plasma insulin and corticosterone levels and no dyslipidemia were observed between control and pollutant-exposed mice. As expected, female mice fed a HF diet showed enhanced body weight, blood glucose, plasma insulin levels, and dyslipidemia as compared to ST-fed mice. Plasma corticosterone levels were in the control range (Fig. 1). receptors (Swedenborg et al., 2009[START_REF] Lee | Glucocorticoid Receptor and Adipocyte Biology[END_REF]. In addition to analyzing these receptors by RT-qPCR, we also studied Hsd11b1 which encodes the 11β-hydroxysteroid dehydrogenase1 (11βHSD1). Indeed in the adipose tissue, 11βHSD1 catalyzes the conversion of the inactive cortisone into the active corticosterone thus locally regulating availability of GCs [START_REF] Sandeep | Pathophysiology of modulation of local glucocorticoid levels by 11beta-hydroxysteroid dehydrogenases[END_REF]. Two-way Anova analysis indicated that the "Pollutant" factor had significant effects on all the genes studied, excepted for Esr1 and Gper1 encoding ERα and GPR30, respectively. The "Diet" factor had significant effects on all the genes studied, but not on Nr3c2 encoding MR. A trend was observed for Esr1 (Fig. 2, insert). Post Hoc analysis determined that pollutant exposure significantly down-regulated Nr3c1 encoding GR, Nr3c2, Hsd11b1and Esr2 (encoding ERβ) in ST-fed but not in HF-fed nutritional conditions. A trend was observed for Ahr (Fig. 2A). Western-blotting analysis of VAT proteins confirmed decreased protein levels of GR in females exposed to the mixture of pollutants (Fig. 2B).

Expression of various

Expression of various nuclear receptor/transcription factors and of 11βHSD1 in

the subcutaneous adipose tissue of female mice fed a ST or a HF diet and exposed to the mixture of pollutants Contrary to the results observed in VAT, exposure to pollutants in SAT did not affect the levels of expression of the glucocorticoid and estrogen receptors, of 11βHSD1 and of AHR (Fig. 3). In addition, the mRNA levels of Nr3c1 and Nr3c2 remained constant in the SAT of females fed a HF diet (Fig. 3). Two-way Anova analysis indicated that the "Diet" factor had significant effects on Hsd11b1 and the estrogen receptors, and a trend was observed for Ahr (Fig. 3, insert). Specifically, as in VAT (Fig. 2), the expression levels of Hsd11b1 and the estrogen receptors Esr1 and Esr2 were significantly decreased and the mRNA levels of Gper1 were significantly increased in the SAT (Fig. 3).

Low-dosed pollutant exposure targeted leptin gene expression in females

To decipher mechanisms underlying changes in glucocorticoid signaling observed in the adipose tissues of standard-fed females, plasma leptin was measured. Indeed, leptin is an adipokine involved in satiety mechanisms which gene is regulated by glucocorticoids (Zhang and Chua, 2017). We observed significant changes in leptin plasma levels (-55%; p < 0.05) together with a small but significant increase in food intake (Fig. 4), despite no changes in the weight of SAT and VAT (Fig. 1). Gene expression analysis showed a significant reduction in Lep mRNA levels in VAT but not in SAT. Adiponectin is another hormone produced by adipocytes. Exposure to pollutants had no effects on adiponectin levels measured in plasma or on the mRNA levels of the adipose tissue depots (Fig. 4). To determine if other markers of adipogenesis displayed altered expression following pollutant exposure, we quantified the mRNA of Nr1c3 encoding PPAR and other target genes including Fabp4 encoding fatty acid binding protein 4 (FABP4) and Cidec encoding fat protein specific 27 (FSP27). Expression mRNA levels were similar to controls in both adipose tissues (Supplementary Fig. 1). We also quantified inflammatory cytokines in VAT and SAT. Pollutants did not affect Tnfa, Ccl2, Ccl5, and Il1b mRNA levels in SAT and VAT (Supplementary Fig. 2).

Pollutant exposure caused changes in the expression levels of glucocorticoid-induced genes including circadian clock genes in the adipose tissues

It is well demontrated that leptin displays a circadian cycling pattern [START_REF] Fu | The molecular clock mediates leptin-regulated bone formation[END_REF] and that adrenalectomy results in changes in leptin secretion (Yilmaz et al., 2002) and in circadian gene expression in VAT (Sotak et al., 2016). Therefore, we quantified a glucocorticoid-and clock controlled gene namely Tsc22d3 encoding the glucocorticoid induced leucine zipper (GILZ) protein [START_REF] Ayyar | Quantitative tissue-specific dynamics of in vivo GILZ mRNA expression and regulation by endogenous and exogenous glucocorticoids[END_REF] and also various core clock genes including Clock, Arntl, Per1, Per2 and Nr1d1 encoding CLOCK, BMAL1, PER proteins and NR1D1, respectively [START_REF] Reppert | Coordination of circadian timing in mammals[END_REF]. Quantification of these genes by RT-qPCR revealed a significant down-regulation of Per1 and Per2 and a trend for Tsc22d3 (p=0.07) in VAT but not in SAT from pollutant-exposed females. Clock and Arntl mRNA levels remained in the control range (Fig. 5A-B).

Analysis of the response to pollutant mixture in male mice

To determine if changes in glucocorticoid signaling also occurred in males in response to the pollutant mixture in ST-fed conditions, we screened all the genes investigated in females, including Nr3c1, Nr3c2, Hsd11b1 as well as Tsc22d3, Clock, Arntl, Per1, Per2 and Nr1d1 in the VAT and the SAT of males fed a standard diet and exposed to the mixture of pollutants. In the VAT of ST-fed males, we observed a significant decrease of Nr3c1 mRNA levels and of Arntl while in the SAT, several genes were significantly up-regulated including Nr3c2, Tsc22d3, Per1, Per2 and Nr1d1 (Fig. 6). Finally, plasma corticosterone levels, food intake, and expression of Lep and Adipoq were analyzed in males (Fig. 7). Plasma corticosterone levels did not change. Contrary to females, food intake and Lep mRNA levels did not change in males but Adipoq mRNA levels were enhanced in the SAT (Fig. 7). The mRNA levels of Nr3c1, Nr3c2 and Hsd11b1 were not affected by exposure to pollutants in the liver, in both sexes (Supplementary Fig. 3).

Discussion

Glucocorticoids exert a strong effect on lipid metabolism via binding to the steroid receptor GR in the adipose tissues and the liver [START_REF] Lee | Glucocorticoid Receptor and Adipocyte Biology[END_REF], Woods et al., 2015) in which 11βHSD1 has predominantly reductase activities allowing local amplification of glucocorticoid action in the absence of changes in plasma corticosterone levels [START_REF] Lee | The glucocorticoid receptor, not the mineralocorticoid receptor, plays the dominant role in adipogenesis and adipokine production in human adipocytes[END_REF]. Although weakly expressed, MR mediates as well metabolic effects related to glucocorticoid binding rather than aldosterone binding in liver and ATs [START_REF] Kuhn | The mineralocorticoid receptor: a new player controlling energy homeostasis[END_REF]. Based on findings presented herein showing coordinated and reduced expression of Nr3c1, Nr3c2 and Hsd11b1 mRNA levels and decreased GR protein level in VAT of females fed a standard diet; we concluded that pollutant exposure caused a defect in glucocorticoid signaling. Importantly, these effects exhibited tissue-specificity because they were not observed in female SAT or in liver of both sexes. They also exhibited sex-specificity. Indeed, while the levels of Nr3c1 were down-regulated in VAT as shown in ST-females, the levels of Hsd11b1 did not change in VAT and the levels of Nr3c2 were up-regulated in the SAT of STmales exposed to pollutants, unlike females.

The origin of the sex-dimorphism probably stems from differences in sexual steroids, estrogens and androgens. It is well described that manifestation of metabolic syndrome caused by GC excess depend on interactions between glucocorticoids and sex hormones [START_REF] Alemany | Do the interactions between glucocorticoids and sex hormones regulate the development of the metabolic syndrome?[END_REF][START_REF] Quinn | Sexually dimorphic actions of glucocorticoids: beyond chromosomes and sex hormones[END_REF]. For example, a high-fat feeding resulted in enhanced corticosterone plasma levels in males (not shown) but not in females (the present study). Origin of the sexdimorphism could as well stem from cross-interactions of AHR, glucocorticoids and estrogen signaling as reviewed [START_REF]Gender differences in glucose homeostasis and diabetes[END_REF], Swedenborg et al., 2009[START_REF] Casals-Casas | Endocrine disruptors: from endocrine to metabolic disruption[END_REF][START_REF] Leblanc | A dual mixture of persistent organic pollutants modifies carbohydrate metabolism in the human hepatic cell line HepaRG[END_REF], raising the hypothesis that a mixture of pollutants could locally contribute to the alteration of GC signaling in females without impacting corticosterone plasma levels. Consistently, we detected mRNA changes in Ahr and Esr2 in the VAT of pollutant-exposed ST-females in addition to the decreased expression of the glucocorticoid receptors. These effects were restricted to the VAT in pollutant-exposed females and distinct from those observed in males [(present study and (Naville et al., 2019)],

suggesting that changes in GC but also in estrogen and AHR signaling in response to exposure to the mixture of pollutants may cause the overall metabolic changes in females, described in the present study. Certain effects caused by pollutant exposure overlapped those described following a high-fat feeding as previously reported [START_REF] Labaronne | Low-dose pollutant mixture triggers metabolic disturbances in female mice leading to common and specific features as compared to a high-fat diet[END_REF]. It included decreased of Nr3c1, Hsd11b1, AhR mRNA levels but also of the estrogen receptor Esr2 (a trend was observed for Esr1) and enhanced expression of Gper1. This may explain why exposure to pollutants did not cause changes in the VAT of females fed a high-fat diet.

Noteworthy, Gper1 increased in SAT and VAT in both females (present study) and males [START_REF] Naville | Chronic exposure to a pollutant mixture at low doses led to tissue-specific metabolic alterations in male mice fed standard and high-fat highsucrose diet[END_REF] in response to a high-fat diet indicating that this estrogen receptor was not regulated in a sex-dependent way.

Glucocorticoids regulate almost every aspect of adipose tissue biology metabolism including adipogenesis, adipocyte metabolism, inflammation and adipokine production preferentially targeting visceral adipose tissue as exemplified in Cushing's syndrome with the development of visceral obesity at the expense of the subcutaneous adipose tissue [START_REF] Lee | The glucocorticoid receptor, not the mineralocorticoid receptor, plays the dominant role in adipogenesis and adipokine production in human adipocytes[END_REF]. Importantly, the pollutant-induced effects described herein in female mice occurred in the absence of changes in the weight of SAT and VAT excluding impact of pollutant exposure on adipogenesis. Accordingly, adipogenic markers remained in the control range. To validate glucocorticoid-signaling disturbances in pollutant-exposed females, we analyzed the mRNA levels of one primary GR target gene presenting glucocorticoid response elements (GREs) in the promoter, Tsc22d3 encoding GILZ for glucocorticoid induced leucine zipper [START_REF] Ayroldi | Glucocorticoid-induced leucine zipper (GILZ): a new important mediator of glucocorticoid action[END_REF]. Tsc22d3 mRNA levels decreased in VAT (a trend) but no changes were observed in SAT in females, thus suggesting pollutants interfered with glucocorticoid signaling in VAT, only. In addition, as pollutants may have pro-inflammatory actions [START_REF] Naville | Metabolic outcome of female mice exposed to a mixture of low-dose pollutants in a diet-induced obesity model[END_REF][START_REF] Kim | Inflammatory pathway genes belong to major targets of persistent organic pollutants in adipose cells[END_REF][START_REF] Pardo | Exposure to air pollution interacts with obesogenic nutrition to induce tissue-specific response patterns[END_REF][START_REF] Rebourcet | The effects of an in utero exposure to 2,3,7,8-tetrachloro-dibenzo-p-dioxin on male reproductive function: identification of Ccl5 as a potential marker[END_REF], we measured mRNA levels of several cytokines in VAT (and SAT) but they remained in the control range. Eventually, we demonstrated that leptin plasma levels and Lep mRNA levels significantly decreased in VAT of the pollutant-exposed females (not in SAT) while adiponectin plasma levels and Adipoq mRNA levels in VAT and SAT remained in the control range. Interestingly, males exposed to the mixture of pollutants exhibited enhanced mRNA levels of several adipogenic markers including Adipoq in SAT but not in VAT and no changes in Lep mRNA levels or leptin plasma levels [the present study and [START_REF] Naville | Chronic exposure to a pollutant mixture at low doses led to tissue-specific metabolic alterations in male mice fed standard and high-fat highsucrose diet[END_REF]], again highlighting sex-differences in the response to pollutant-exposure in metabolicallyactive tissues. Further studies will aim at understanding the coordinated regulation of Nr3c2 and Tsc22d3 mRNA levels in the SAT of exposed males together with the enhanced expression of several adipogenic markers including Nr1c3, Fabp4 and Cidec described elsewhere [START_REF] Naville | Chronic exposure to a pollutant mixture at low doses led to tissue-specific metabolic alterations in male mice fed standard and high-fat highsucrose diet[END_REF].

Leptin level was the only significant metabolic phenotype difference measured in females in response to exposure to the mixture of pollutants. Leptin is an adipokine that exerts satiety effects and reduced leptin fits with enhanced food intake as we observed in females but not in males. Leptin is expressed as a function of adipose tissue mass and is more expressed in females than in males (Zhang and Chua, 2017); thus, reduced leptin levels in basal conditions of feeding without gain of weight despite enhanced food intake could indicate enhanced energy expenditure and thermogenesis. Since we did not measure these parameters in the present study, validation of the hypothesis will await further studies. A second hypothesis to explain the reduction of glucocorticoid signaling and leptin expression in females exposed to pollutants could be a possible dysregulated peripheral clock in VAT. There are several clues in the literature to support this hypothesis. Glucocorticoids are secreted in circadian cycles to coordinate the circadian rhythm with energy balance and leptin synthesis (So et al., 2009).

Adrenalectomy diversely affected peripheral circadian clocks including GILZ expression in the adipose tissue (Sotak et al., 2016) and led to decreased serum leptin (Yilmaz et al., 2002).

The powerful GR agonist, dexamethasone led to the disappearance of clock gene rhythmicity in adipose tissues [START_REF] Gomez-Abellan | Glucocorticoids affect 24 h clock genes expression in human adipose tissue explant cultures[END_REF].

The circadian clock regulates metabolism and energy homeostasis in peripheral tissues including the adipose tissues and synchronization with the master clock in the hypothalamic suprachiasmatic nuclei (SCN) would be in part orchestrated by glucocorticoids to keep in step diurnal functions with the environmental cycles of light and nutrition [START_REF] Dickmeis | Glucocorticoids and the circadian clock[END_REF][START_REF] Balsalobre | Resetting of circadian time in peripheral tissues by glucocorticoid signaling[END_REF]. Circadian rhythms are generated by feedback loops in core clock gene expression [START_REF] Reinke | Crosstalk between metabolism and circadian clocks[END_REF]. Among core clock components investigated, we demonstrated herein significant reduction of Per1 and Per2 mRNA levels in VAT but not in SAT in response to pollutant exposure in females, and no effect on Clock and Arntl.

Interestingly, both Per1 and Per2 mRNA levels enhanced in the SAT (not in VAT) of males exposed to the mixture of pollutants. These findings are consistent with the presence of GRE elements in the promoters of Per1 and Per2 genes [START_REF] Koyanagi | Glucocorticoid regulation of 24-hour oscillation in interferon receptor gene expression in mouse liver[END_REF], So et al., 2009).

They also validate changes in GC signaling in the VAT of exposed females and in the SAT of exposed males. Overall, these findings pointed that several peripheral circadian clock-related genes were diversely affected by exposure to a mixture of pollutants in a sex-and tissuedependent way (at least at the mRNA level in the ATs) extending previous reports [START_REF] Labaronne | Low-dose pollutant mixture triggers metabolic disturbances in female mice leading to common and specific features as compared to a high-fat diet[END_REF] [START_REF] Shen | Low-dose PCB126 compromises circadian rhythms associated with disordered glucose and lipid metabolism in mice[END_REF][START_REF] Kopp | Exposure to endocrine disrupting chemicals perturbs lipid metabolism and circadian rhythms[END_REF]. Further studies are needed to explore why VAT in females and SAT in males were differentially regulated in response to the mixture of pollutants in a context of standard diet feeding, and why SAT in females was not affected. It may explain the mild metabolic phenotype described in females (the present study). Indeed, circadian metabolic disturbances may be early events and exposure to pollutants greater than 15 weeks could have resulted in more harmful phenotypes than what was described herein. This hypothesis merits further attention.

Consistent with the definition that endocrine disruptors may influence any hormonal action to mediate their disruptive effects, we also identified that the mixture of pollutants targeted 11βHSD1 in addition to corticoid receptors. These findings extend previous reports that endocrine disrupters can affect enzymes involved in hormone metabolism including the steroidogenic acute regulatory protein (STAR) [START_REF] Guyot | Diethylstilbestrol inhibits the expression of the steroidogenic acute regulatory protein in mouse fetal testis[END_REF], the estrogen sulfotransferase (SULT1E1) [START_REF] Naville | Low-dose food contaminants trigger sex-specific, hepatic metabolic changes in the progeny of obese mice[END_REF], aromatase and others [START_REF] Beausoleil | Regulatory identification of BPA as an endocrine disruptor: Context and methodology[END_REF]) (Zoeller et al., 2012, Nadal et al., 2017). Noteworthy, pollutants used in the mixture were incorporated at a dose lower than their no-observed adverse effect level (NoAEL) and just above their tolerable daily intake (TDI) dose. The TDI is supposedly the dose without adverse effect in humans and it is derived from the NoAELs determined through experimental studies applying corrective factors to take into account inter-and intra-species variations [START_REF] Dorne | Metabolism, variability and risk assessment[END_REF] of 100 to 1000 (Supp. Table 1). Therefore, this study gave further evidences that the risk assessment procedure which is based on NoAELS for single chemicals may not be protective enough. In addition, comparison of metabolic disturbances elicited by the mixture versus pollutants individually will allow determining if the observed effect resulted from the action of one or several pollutants or from a cocktail effect.

Conclusions

We here demonstrated that exposure to a mixture of low-dosed pollutants caused reduced expression of glucocorticoid receptors in the visceral but not in the subcutaneous adipose tissue or in the liver in female mice and that it was associated with decreased plasma levels of leptin and reduced expression of Lep and of various canonical genes related to the circadian clock machinery in visceral adipose tissue. Inasmuch as lowering glucocorticoid signaling following adrenalectomy prevents physiological oscillations of circadian genes and impact adipokine secretion (Sotak et al., 2016, Yilmaz et al., 2002), our data fuel the hypothesis that a mixture of low-dosed chemicals may alter clock genes resulting in metabolic effects at least in the visceral adipose tissue in female mice.

  nuclear receptor/transcription factors and of 11βhydroxysteroid dehydrogenase1 (11βHSD1) in the visceral adipose tissues of female mice fed a ST or a HF diet and exposed to the mixture of pollutants Adipose tissues express several nuclear receptor/transcription factors that pollutants can bind to initiate downstream activation signals among which the estrogen receptors, the aryl hydrocarbon receptor (AhR) but also the glucocorticoid (GR) and the mineralocorticoid (MR)

Figure 1 :

 1 Figure 1: Metabolic phenotyping of female mice fed a standard (ST) or a high-fat (HF) diet

Figure 2 :

 2 Figure 2 : Effect of the mixture of pollutants and of the high-fat high-sucrose (HF) diet on the

Figure 3 :

 3 Figure 3 : Effect of the mixture of pollutants and of the high-fat high-sucrose (HF) diet on the

Figure 4 :

 4 Figure 4 : Effect of the mixture of pollutants on plasma leptin and adiponectin levels, food

Figure 5 :

 5 Figure 5 : Effect of the mixture of pollutants on the expression of Tsc22d3 and the circadian

Figure 6 :

 6 Figure 6: Effect of the mixture of pollutants on the expression of Nr3c1, Nr3c2, Hsd11b1,

Figure 7 :

 7 Figure 7 : Effect of the mixture of pollutants on food intake and corticosterone plasma levels

Acknowledgments

This work was supported by INSERM and University of Lyon1 to INSERM U1060. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. We gratefully acknowledge Sandra Pesenti for the technical help with PCR primers and Patrick Manas for animal care. The authors wish to thank Dr.

Alexandrine Derrien-Colemyn for critically reading the manuscript.

Declaration of competing interest

The authors declare that there are no conflicts of interest.

health protection: a statement of principles from The Endocrine Society. Endocrinology, 153, 4097-110.

Supplementary File

Supplementary Table 1: References doses of the pollutants present in the mixture and doses used in the diets (all per kg and per day) (From Naville et al., 2019).