Pierre-Antoine Laharotte

CEREMA Raballand Wilfried

COMBINING USUAL ROADSIDE SURVEY AND CELLULAR PHONE DATA TO PRODUCE ACCURATE TRANSIT AND EXCHANGE OD MATRICES -APPLICATION TO THE CROSS-BORDER GENEVA 4-STEP MODEL

INTRODUCTION

Knowing the mobility of people is a major challenge to fuel traffic or travel models in order to assess town policies. This knowledge of mobility involves traffic surveys to collect the origins and destinations of users. Owing to the cost of surveys and the emergence of electronic toll collection, it is no longer possible to carry out field surveys on French motorways. Faced with these difficulties, the use of data from on-board technologies is a good way to explore mobility, since it can be used to reconstruct traces of vehicles. They have been assessed in order to estimate the demand for travel on the network [McNeill et al (2017), Bonnel et al (2015)] or to better understand local dynamics [Wang et al (2017), McNeill et al (2017), Nabavi et al (2015), Caceres et al (2008)]. Usually, in the field of transport, digital data is shared between:

• data from on-board geolocation systems (Floating Car Data -FCD),

• data from mobile telephony (Floating Mobile Data -FMD) based on a network of antennas (Edge,2G,3G,4G,Bluetooth,WIFI,etc.). FCD data turns out to be very efficient in reconstructing traffic dynamics [Herring et al (2010), Vanajakshi et al (2009)], especially the speed of the vehicles. Nevertheless, the volumes collected by OD (Origin-Destination) from these tracks require significant corrections aiming to offset the focus made on specific population strata. FMD data do not suffer from such a bias in representativeness and can be used to reconstruct OD matrices from antenna to antenna [Caceres et al (2008), Michau et al (2014)]. However, the projection of these OD flows onto specific zoning and fine reproduction of secondary flows still arise some difficulties [Caceres et al (2008), Raballand (2018)]. The ability of these new data sources to be used instead of conventional roadside surveys requires further methodological refinement before they can be implemented fully. In order to guarantee operational implementation, we introduce a hybrid approach, combining usual surveys and digital collections provided in the form of OD matrices. The aim is to build an OD matrix of transit or exchange trips on the perimeter of Geneva's travel model.

METHODOLOGY

General principle of the introduced protocol

The proposed protocol aims at carrying out a traffic survey picking up the transit and exchange flows over a wide targeted area: the Geneva model area. The usual protocol consists of positioning survey stations on each injection point over the perimeter of the targeted area, i.e. any of the road infrastructures allowing the access to the area. Each conventional survey station consists of roadside investigators interviewing users, and automatic vehicle counting to correct the survey bias.

We propose to set up conventional survey stations as far as possible, i.e. mainly on national and local roads. Then we deal with the motorway cases by setting up digital stations (see Fig. 1). Digital stations are synthetic stations, reconstructed from mobile phone data collections and rectified from automatic counts and number plate readings (NPR).

Specific features of data collected from mobile telephony

The collection process of mobile telephony data requires to define:

• a digital detection area to select all the trips to store within the digital OD matrix (see Fig. 2, green area);

• a slack zoning, known as TEL zoning, guaranteeing a good fit between how the area is divided up and the positioning of the mobile telephone antennas; • a minimum period of immobility to characterize the beginning and the end of a trip, in order to set its origin and destination zone. Such process enables to track mobile phone users and retrieve OD matrices for any timestamp on a specific zoning. Nevertheless, the antenna networks were not designed for this purpose, some bias are therefore introduced:

• the digital detection zone does not fit only the targeted motorway area: it equally picks up motorway users as users with other modes (trains) or on nearby roads (secondary network), as represented in Fig. 3; • TEL zoning is dependent on the layout of the antenna network and does not coincide with the administrative zoning, which requires defining flow distribution rules to move from one zone to the other. The final reconstitution of a digital survey station requires the specific features of these digital data collections to be integrated.

2.3

Theoretical concept to achieve a hybrid road survey Assuming the OD matrices resulting from conventional roadside survey as the reference outputs, some adjustments are required on the digital OD matrices. It is indeed necessary to tackle the issue of:

• mixed mode and mixed road category within the data resulting from digital surveys; • zoning differences between digital surveys (TEL zoning) and conventional roadside survey (administrative zoning).

Mode and road category discrimination process applied to OD flows resulting from digital stations

The targeted motorway OD matrix 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑣𝑒ℎ 𝑚𝑜𝑡𝑜𝑟𝑤𝑎𝑦 (𝑜𝑑) results of the subtraction from the raw OD matrix 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 (𝑜𝑑) of:

• a matrix composed of OD flows observed on the surrounding of the motorway, but driven with alternative modes (ex: train): 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 𝑜𝑡ℎ𝑒𝑟 𝑚𝑜𝑑𝑒 (𝑜𝑑);

• a matrix composed of OD flows observed on the surrounding roads of the motorway: 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑣𝑒ℎ 𝑜𝑡ℎ𝑒𝑟 𝑟𝑜𝑎𝑑𝑠 (𝑜𝑑). These flows need to be removed as they are iterated within the global study: they appear as well within the conventional survey as within the digital survey.

a / Assessment of the OD matrix associated to alternative modes

It is easier to collect available road counts about the roads passing through the targeted area than to collect exhaustive and accurate data about alternative modes. It is assumed that the alternative modes matrix 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 𝑜𝑡ℎ𝑒𝑟 𝑚𝑜𝑑𝑒 (𝑜𝑑) results from the difference between the raw matrix 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 (𝑜𝑑) and an assessment of the road matrix 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 𝑟𝑜𝑎𝑑 (𝑜𝑑), as related in Eq.1. The resulting matrix 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 𝑜𝑡ℎ𝑒𝑟 𝑚𝑜𝑑𝑒 (𝑜𝑑) is punctually compared to available alternative modes data.

The assessment of the road matrix 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 𝑟𝑜𝑎𝑑 (𝑜𝑑), expressed in persons, results from a Fratar type algorithm [Fratar (1954)] applied to the raw matrix 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 . The margins are assigned with the combination of roadside counts surrounding the digital detection areas. This relationship is described by Eq. 2, where 𝜆 𝑜𝑐𝑐𝑢𝑝𝑎𝑛𝑐𝑦 is the average vehicle occupancy rate and 𝜎(𝑜𝑑) is an adjustment coefficient. Thanks to the road count, the unit is switched from person to vehicles. In practice the C and 𝛼 𝑘 contributions are directly assessed through the implementation of a simulator. The simulator aims to reproduce the detection process of the digital stations. It applies this detection process to the OD trips from the conventional roadside surveys, previously assigned to the road network. Finally, the simulator is composed of:

• a detection area associated to any digital survey and depicted as a circular area defined by its radius R; • a road assignment process enabling to make assumptions about the possible route for each OD relationship from the conventional roadside surveys. The matrix 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑣𝑒ℎ 𝑜𝑡ℎ𝑒𝑟 𝑟𝑜𝑎𝑑𝑠 results from the intersection between the simulated detection area and the possible routes. The optimum radius R is a parameter to adjust thanks to digital survey results and further available data such as number plate readings.

c / Conversion of the motorway matrix from TEL zoning to administrative zoning

The motorway OD matrix 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑣𝑒ℎ 𝑚𝑜𝑡𝑜𝑟𝑤𝑎𝑦 is built on the TEL zoning. In order to ensure the communication with OD matrices resulting from conventional surveys, it must be allocated on the administrative zoning. It is equivalent to apply an inverse basis change C -1 , described in concrete terms in paragraph 2.3.2.

Principle for converting TEL zoning to administrative zoning

The inconsistency between the coverage topology of the mobile phone network and the administrative division of the area may lead to certain flow distribution biases, especially when a high-density zone is divided into several parts. The transition between the mobile phone collection zoning and the fine zoning of the model is achieved through a breakdown method based on • explanatory data on the generation (population) and attractiveness (employments) specific to the fine zones; • significant divergences of division affecting highly urbanized and/or adjacent areas.

Non-optimal location of the digital stations

As tracking the mobile phone user is not the main purpose of the phone network antennas, the location and orientation of the antennas could not match with the optimal location of the digital survey stations. Some digital stations could have to be located within the perimeter of the study. It involves that some internal flows are detected by the digital station. Most of these flows are properly assigned to internal flows within the digital matrix. As the study only focuses the external flows, they are removed. Nevertheless, to ensure the results and improve the simulator calibration, NPR have been set up on the perimeter of the concerned detection areas.

Partial coverage of the area by the operator's antenna network

The method set out in 2.3 assumes the continuity of coverage of the mobile phone network, but mobile phone operators can only track their users in their respective countries. Carrying out a crossing-border survey, i.e. on a perimeter covering 2 countries, requires performing 2 joint, simultaneous surveys, and therefore using data from two mobile phone operators. As the Geneva model is a crossing border model, a French operator (Bouygues Telecom) and a Swiss operator (Swisscom) have been asked to carry out a digital survey within their respective territories.

To ensure the continuity of the crossing-border flows within the data, digital stations have been set up at the border by both phone operators. A former border survey is available to check the consistency of the results. Any digital station at the border can pick up movements between the two countries. The connection between the 2 digital surveys can be made, after adjustment of the digital OD matrices, by: • applying a probabilistic assignment according to a gravity model,

• or using the distribution from cross-border travels surveys.

Figure 6: Implementation of the full survey on the case of study of Geneva

OPERATIONAL IMPLEMENTATION AND RESULTS

The methodology introduced in section 2 has been implemented on the Geneva's model perimeter to retrieve the exchange and transit flows driving through the area. It requires to set up:

-a conventional roadside survey composed of 19 stations; -a digital survey composed of 6 stations (3 stations at the border, 3 stations on the perimeter) on the French side; -a digital survey composed of 4 stations (3 stations at the border, 1 station on the perimeter) on the Swiss side.

Analysis of the raw digital surveys 'outputs

The preliminary analysis of the digital raw OD matrices 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 provided by the phone operators highlights 2 sources of inconsistency or bias:

• the flows of crossing-border peoples and its distributions provided by the digital surveys do not match the automatic counts, nor the trends given by the former conventional survey at the border, assumed to be a reference (see Table 1); • the spatio-temporal instability of the TEL zoning involves some mistakes in the distribution of the flows over the zoning (see The application of the mode and road type discrimination process to the current case study enables to split the original flow provided by phone operators into 3 clusters : the alternative mode (ex: railway) flow, the motorway flow and the flow driven on local roads (see Table 3). As, in our case study, the alternative mode flow is assumed to mainly result from trains, this flow has been compared, on the proper areas, to a rough estimation of users on the TGV line LYRIA (Paris-Geneva) for peak hours. It results in a difference of MAPE=55% in the direction of Geneva and MAPE=19% in the reverse direction. The gap may be explained by the rough assessment and definition of the alternative mode. The discrimination according to the driven road type is made thanks to a simulation of the digital detection area and the road assignment of the conventional survey flows. The main parameter of the simulator is the radius R of the digital detection area, which has been assessed around 4000m on the station #4 thanks to a complimentary NPR survey.

Building the exchange and transit OD matrix

Once the motorway OD matrices are built on both sides of the border (French and swiss), the junction can be made on the 3 digital stations at the border thanks to 3 specific probabilistic matrices.

CONCLUSIONS AND PERSPECTIVES

This work develops a method of using mobile phone commercial data to analyse the mobility of people in an area. In particular, the approach is appropriate to meet the specific needs of carrying out a traffic survey by

Figure 1 :

 1 Figure 1: Protocol layout

Figure

  Figure 2: mobile phone detection process with digital survey

  (1) 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 𝑜𝑡ℎ𝑒𝑟 𝑚𝑜𝑑𝑒 (𝑜𝑑) = 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 (𝑜𝑑) -𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 𝑟𝑜𝑎𝑑 (𝑜𝑑) (2) 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 𝑟𝑜𝑎𝑑 (𝑜𝑑) = 𝜆 𝑜𝑐𝑐𝑢𝑝𝑎𝑛𝑐𝑦 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑣𝑒ℎ 𝑟𝑜𝑎𝑑 (𝑜𝑑) = 𝜎(𝑜𝑑) 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑝𝑒𝑟𝑠 (𝑜𝑑) b / Identification of the iterated flows within the road matrix The resulting road matrix 𝑀 𝑑𝑖𝑔𝑖𝑡-𝑣𝑒ℎ 𝑟𝑜𝑎𝑑 , expressed in vehicles, is composed of the motorway flows (𝑀 𝑑𝑖𝑔𝑖𝑡-𝑣𝑒ℎ 𝑚𝑜𝑡𝑜𝑟𝑤𝑎𝑦 ) and parasite flows (𝑀 𝑑𝑖𝑔𝑖𝑡-𝑣𝑒ℎ 𝑜𝑡ℎ𝑒𝑟 𝑟𝑜𝑎𝑑𝑠 ), assumed to be already stored within conventional roadside surveys (𝑇 𝑠𝑢𝑟𝑣𝑒𝑦 𝑘), as related in Eq. 3. The relationship between the iterated flows and flows resulting from the roadside surveys is described in Eq. 4 and Figure4, where 𝛼 𝑘 × 𝑇 𝑠𝑢𝑟𝑣𝑒𝑦 𝑘 represent the proportion of flows collected both by the conventional roadside survey station k and a digital station and C is a change of basis matrix to convert the administrative zoning into a TEL zoning.

Figure

  Figure 4: digital detection area outside the model delimitation

  Any probabilistic matrix describes the probability p(o,d) of a crossing-border trip to begin in an zone o and to end in the zone d. The probabilistic matrices result from the combination of a former (2011) crossing-border survey and the current socioeconomic data from INSEE. The use of a former survey enables to refine any gravity model with local and specific trends. The junction is made between the both motorway matrices with respect to the flow conservation assumption. Final aggregated exchange and transit matrix -Light Vehicles -Morning Peak Hour (MPH)3.4 Digital survey versus conventional roadside surveyTo better feature the digital station abilities, one location at the border (station #3) is covered with 2 kind of survey: a conventional survey and a digital survey. The figure7illustrates the differences between the conventional roadside survey results and the outputs of the digital survey before/after adjustments. Only the French side of the survey is shown, it means that it figures out the distribution of the flows crossing the border through the station #3 to reach Switzerland and coming from a French area (inner area IFX or outer area ZFi). Flows are aggregated by origin. It appears that the main trends are well identified by the mobile phone data (main flows on IF1, IF2 and IF3), but the residual flows are significant. These residual flows are due to the instability of the TEL zoning. Furthermore, as mentioned in previous works [Raballand(2018)], the assessment of the lowest flow by mobile phone data is inconsistency.

Figure 7 :

 7 Figure 7: Flow distribution over the French zones for digital station #3 -MPH

Table 2

 2

	). Especially

Table 2 :

 2 Impact of the instability of TEL zoning

3.2 Isolation of the motorway flows

	French station	Total flow	flows assigned to	Flows assigned
	(as origin)	(raw data)	other roads	to railway
	Station #4	13 757	1 864	NA
	Station #5	4 058	1 086	NA
	Station #6	3 090	270	52

Table 3 :

 3 Impact of the process of mode / road type discriminationnumber of people -Morning Peak Period (MPP)

substituting "digital" stations for field survey stations that are not feasible on the motorway. In the future, work to make the method more reliable and expand the use of mobile data is needed to confirm the effectiveness and robustness of this innovative approach.

BIBLIOGRAPHY


Antoine Legouhy

Olivier Commowick

François Rousseau

Christian Barillot

Online atlasing using an iterative centroid

come

Introduction

An atlas is a generic term covering several uses. Overall, it corresponds to an average model of the brain both in terms of shape and intensity. It is therefore a powerful tool to understand brain variability, to compute statistics on populations [START_REF] Whitcher | Statistical group comparison of diffusion tensors via multivariate hypothesis testing[END_REF] or to segment regions [START_REF] Cabezas | A review of atlasbased segmentation for magnetic resonance brain images[END_REF]. Ideally, atlases are made from a large number of subjects by combining their images through registration and transformation composition which can be computationally quite costly.

Among the many available atlas construction methods, Guimond et al. [START_REF] Guimond | Average brain models: A convergence study[END_REF] proposed to register the subjects onto an arbitrary image, followed by an unbiasing and averaging step, leading to an atlas that becomes the new reference image. By iterating the process, the atlas becomes less and less biased by the initial reference image. Joshi et al. and Rohlfing et al. [START_REF] Joshi | Unbiased diffeomorphic atlas construction for computational anatomy[END_REF][START_REF] Rohlfing | The SRI24 Multichannel Atlas of Normal Adult Human Brain Structure[END_REF] among others focused on another atlas construction class, called template free approaches, working with groupwise registration. Although those methods have the advantage of being unbiased by construction, their computational cost can be prohibitive.

The research leading to these results has been supported by the ANR MAIA project, grant ANR-15-CE23-0009 of the French National Research Agency (http: //recherche.imt-atlantique.fr/maia) and La Région Bretagne.

Historically databases used for atlas construction were small (about 50 to 100 subjects at most). However, there is a recent trend for the creation of large databases, out of which great atlases depicting the variability of the population across the life span could be inferred. Among these databases, the human connectome project (HCP, now completed) is made of around 1200 subjects of the same age range [START_REF] Essen | The human connectome project: A data acquisition perspective[END_REF]. Even bigger databases are currently being constituted from more subjects (up to 100,000 in the UK biobank [START_REF] Miller | Multimodal population brain imaging in the UK biobank prospective epidemiological study[END_REF]) or across the life span (HCP lifespan) on which researchers already wish to start studies. This causes a problem as the databases are still being gathered and thus the atlas needs to be updated. With current methods, the atlas needs to be recomputed each time new acquisitions are added which is very costly. New methods are therefore needed to allow online atlasing, i.e. the capacity of incrementing atlases with new images as they come without having to perform again the whole atlas construction, as opposed to the above-mentioned direct atlasing.

We propose a new method that considers atlasing from a different angle. Inspired by the notion of centroid for a set of points in a Euclidean space and the way it can be constructed following an iterative procedure, Cury et al. [START_REF] Cury | Diffeomorphic Iterative Centroid Methods for Template Estimation on Large Datasets[END_REF] derived a diffeomorphic atlas construction method working on surfaces. We propose here a generalization of this iterative centroid method for image atlasing where the centroid (atlas) is updated gradually as new images are incorporated. The online aspect of this technique has several advantages: when a new images arrives, its integration does not necessitate to restart the atlasing process from scratch. Also, it only requires one registration per new image making it computationally cheap. Benefitting from the log-Euclidean framework for diffeomorphisms from Arsigny et al. [START_REF] Arsigny | A Log-Euclidean Framework for Statistics on Diffeomorphisms[END_REF], the method allows the production of a diffeomorphic atlas unbiased up to a linear transformation (rigid or affine at the choosing of the user).

We first introduce the iterative centroid method and its specificities for image atlasing in Section 2. Then in Section 3, we present experiments comparing how the atlases built using direct atlasing (here a modified Guimond et al. method [START_REF] Legouhy | Unbiased longitudinal brain atlas creation using robust linear registration and log-Euclidean framework for diffeomorphisms[END_REF]) and with our method differ from each other as they include more and more subjects. We demonstrate with three different metrics that the iterative centroid and direct atlases diverge only by a small amount in terms of transformations and that their intensity based features do not vary.

Method

Theoretical background

Considering the registration of an image J onto an image I, we denote by J the registered image and T the transformation operating on coordinates x used to resample J onto I: I(x) ∼ J(x) = J • T (x). Assuming that T lives in a connected Riemannian manifold, we consider γ(t) the geodesic between Id and T such that γ(0) = Id and γ(1) = T . In this framework (named large deformation diffeomorphic metric mapping -LDDMM), the deformation T is parameterized by a time varying vector field integrated over time (the geodesic γ).

Diffeomorphisms are widely used in non linear registration because of their interesting properties: differentiability, bijectivity and differentiability of the inverse map. As the LDDMM is a computationally intensive framework, Arsigny et al. [START_REF] Arsigny | A Log-Euclidean Framework for Statistics on Diffeomorphisms[END_REF] presented a log-Euclidean framework to compute statistics on diffeomorphisms parameterized by a time independent i.e. stationary velocity field (SVF). For such diffeomorphisms, operations are defined as their Euclidean counterparts on the SVFs. One of interest is the power of a diffeomorphism defined as T α = γ(α) = exp (α log(T )). In the following, all non linear transformations are assumed to be diffeomorphisms parameterized with SVFs.

Additionally, given 2 SVFs v and w, it is possible to compute the logarithm of their composition while avoiding a very costly logarithm computation. This is made possible by the Baker-Campbell-Hausdorff (BCH) formula stating that in a BCH-Lie group, for v and w small enough, log(exp(v) • exp(w)) can be approximated as a series of Lie brackets:

log(exp(v) • exp(w)) ≈ BCH(v, w) = v + w + 1 2 [v, w] + . . . (1)
Where [v, w](x) = Jac(v)(x).w(x) -Jac(w)(x).v(x). It has been shown in [START_REF] Bossa | Contributions to 3D Diffeomorphic Atlas Estimation: Application to Brain Images[END_REF][START_REF] Vercauteren | Symmetric Log-Domain Diffeomorphic Registration: A Demons-Based Approach[END_REF] that the use of the BCH formula is well suited for diffeomorphisms parameterized as SVFs. To simplify notations, unless specified otherwise, the composition of two transformations associated to two SVFs v and w is performed via the BCH: exp(v) • exp(w) = exp(BCH(v, w)).

Iterative centroid atlas construction

In a Euclidean space Ω, the centroid b n of a set of points {x 1 , . . . ,

x n } is defined as b n = argmin y∈Ω n i=1 y-x i 2 . A direct solution is 1 n n i=1
x i but an iterative approach is also imaginable. Starting from b 1 = x 1 , we can use the recursion:

b k+1 = k k + 1 b k + 1 k + 1 x k+1 (2)
As depicted in [START_REF] Cury | Diffeomorphic Iterative Centroid Methods for Template Estimation on Large Datasets[END_REF], one can extend this formulation to the Riemannian case for surfaces by using the generalization of the notion of straight line through geodesics. The updated centroid b k+1 , in that case, is on the geodesic linking b k and x k+1 , at 1 k+1 times the distance from b k . We define our iterative centroid atlasing, following the same spirit, on images and diffeomorphisms parameterized by SVFs. Thereafter we assume that the registration of an image I onto an atlas A provides two transformations L and T , where L is linear and T is a diffeomorphism, so that A ∼ I •L•T . Generalizing Eq. 2 to a set of images {I 1 , . . . , I n }, we start from A 1 = I 1 and get:

A k+1 = k k + 1 A k • T -1 k+1 k+1 + 1 k + 1 I k+1 • L k+1 • T k k+1 k+1 (3)
The atlasing scheme is illustrated Fig. 1. We apply to A k a transformation situated at 1 k+1 of the distance from the identity to T -1 k+1 along the geodesic. We also apply to I k+1 , in addition to the linear transformation, a transformation situated at k k+1 of the distance from the identity to T k+1 along the geodesic. Also, A k is made from k images so a weight k k+1 is affected to its intensities while I k+1 intensities get a weight 1 k+1 . At the end, those two images are transformed and weighted according to their contribution to the atlas.

. . . Unlike [START_REF] Cury | Diffeomorphic Iterative Centroid Methods for Template Estimation on Large Datasets[END_REF], we are dealing with images and therefore interpolations on a voxel grid to reconstruct a new image after a transformation. This resampling step has to be avoided as much as possible to keep a maximum of details. To this end it is preferable not to use A k (already a resampled image) when constructing A k+1 , but instead to operate on the initial images using transformation compositions such that each image undergoes only one resampling. This is achieved by expressing the sequence (A k ) k∈N using only the initial images {I j , j = 1, . . . , k}. Rearranging Eq. 3, the following formulation emerges:

A k+1 = 1 k + 1   k j=1 I j • L j • Θ j,k • T -1 k+1 k+1 + I k+1 • L k+1 • T k k+1 k+1   (4)
Where θ j,k is the non linear part of the transformation bringing I j onto A k . There are k compositions using BCH to perform at iteration k which correspond to the update of the transformations θ j,k for j ≤ k.

Initialization: One can start at iteration 1 by initiating A 1 = I 1 and θ 1,1 = Id. However, a very interesting point of the method is the possibility to complement an already existing atlas A p constituted of p images as long as it can be written as A p = 1 p p j=1 I j • L j • T j . To do so, we simply assign θ j,p = T j , ∀j ≤ p.

Iterative procedure: At the end of iteration k, we have, for each j ≤ k, transformations θ j,k that map images I j onto the atlas

A k : A k ∼ I j • L j • θ j,k .
To build A k+1 from these, we do the following: 1-we register I k+1 onto A k : A k ∼ I k+1 •L k+1 •T k+1 (one registration) ; 2-for j ≤ k, we compute the updated θ j,k+1 as θ j,k+1 = θ j,k • T

Unbiased atlas construction

So far, we have only assumed that the registration produced a linear and a non linear part. Going in more details, let us now assume that this registration can be written A ∼ I •R•S •D where: R•S is an affine transformation decomposed into a rigid part R and a stretching part S, and D is a diffeomorphism. Herebefore, by taking L = R • S and T = D i.e. L encapsulating all the linear transformations and T only local displacements, the method will lead -similarly to [START_REF] Guimond | Average brain models: A convergence study[END_REF] -to an atlas unbiased up to an affine transformation. On the other hand, if we take L = R and T = S • D i.e. L only considering rigid motion while T encapsulating both growth and local displacements, the method will lead to an atlas unbiased up to a rigid transformation (useful typically for longitudinal studies).

Experiments and results

MRI database

For our experiments, we have used T1 weighted preprocessed images (brain extracted) from 100 subjects randomly chosen among the Human Connectome Project (HCP) database (for more details, see Van Essen et al. [START_REF] Essen | The human connectome project: A data acquisition perspective[END_REF]). The size of the images is 260 × 311 × 260 with a voxel size of 0.7 mm isotropic.

Implementation details

Registration Our method is agnostic to the registration used. We performed registrations in two steps using symmetric block matching algorithms [START_REF] Commowick | Block-matching strategies for rigid registration of multimodal medical images[END_REF][START_REF] Commowick | Automated diffeomorphic registration of anatomical structures with rigid parts: application to dynamic cervical MRI[END_REF] available in Anima1 :

1. An affine registration that outputs a transformation matrix B which is decomposed through polar decomposition into a rigid part R and a stretching part S as depicted in Legouhy et al. [START_REF] Legouhy | Unbiased longitudinal brain atlas creation using robust linear registration and log-Euclidean framework for diffeomorphisms[END_REF]. 2. A diffeomorphic registration that directly outputs the SVF of the diffeomorphism log(D) linking the two images. This allows to take advantage of the log-Euclidean framework without logarithm computation [START_REF] Vercauteren | Symmetric Log-Domain Diffeomorphic Registration: A Demons-Based Approach[END_REF][START_REF] Bossa | Contributions to 3D Diffeomorphic Atlas Estimation: Application to Brain Images[END_REF].

We chose to compute atlases up to a rigid transformation. Therefore, the linear transformation L is rigid: L = R and the non linear transformation is the combination of stretching and local displacements:

T = S • D.
Transformation composition Compositions of two diffeomorphisms are performed on their SVFs using order 2 BCH. To compose a linear transformation L and a SVF, we first apply log(L) to the real coordinates of an image grid to transform it into an SVF and then compose using order 2 BCH.

Comparison between iterative centroid and direct atlasing

In the following experiments, we evaluate the quality of the iterative centroid atlas and its divergence with respect to direct atlasing, here a method developed in Legouhy et al. [START_REF] Legouhy | Unbiased longitudinal brain atlas creation using robust linear registration and log-Euclidean framework for diffeomorphisms[END_REF] available in Anima-Scripts2 based on Guimond et al. [START_REF] Guimond | Average brain models: A convergence study[END_REF].

For both methods, we used the same registration and composition tools with the same parameters. Starting from an atlas built with 50 images using direct atlasing, we added one by one the 50 other images using our online atlasing. We then compared the direct atlases computed at several steps with the online ones using three metrics.

Divergence Both methods output linear and non linear transformations (denoted T j for online atlasing and Tj for direct atlasing) mapping each image onto the atlas. We first propose to evaluate the divergence between the two atlases A and à by an image δ defined from the transformations: δ A, à = This measure has the advantage of being agnostic to edge issues that can happen in intensities comparisons. Fig. 2 presents a qualitative view of the evolution of the differences between the direct and online atlases with respect to the number of subjects. We observe that most of the differences increase with the number of added subjects and occur in the cortical areas while central structures of the brain are spared. Fig. 3 gives quantitative plots over the brain of the divergence measures. It suggests that the divergence tends to grow (from a median of 0.777 mm for 60 images, to 1.095 mm for 100) but at a slow pace that decreases (+0.153 mm between 60 and 70 images but only +0.025 mm between 90 and 100).

Iconic evaluations

The second part of the evaluation assesses the quality of the atlases based on image intensities. We evaluate first the atlas standard deviation across subject images after transformation on the atlas, the transformations being constructed differently upon the atlasing method. This measure allows to quantify how well the images are aligned after atlas construction (the lower the standard deviation, the better): sdis

(A) = 1 n n i=j (I j • L j • T j -A) 2 .
For better interpretation, we normalized sdis by the average of the intensities of the atlases from the two methods to create sdis n .

In addition, we propose a measure of the sharpness of the atlases. Given N (i), a patch around a voxel i of an image A, we define a local measure of the sharpness at voxel i, sha(A(i)), as the standard deviation of the patch normalized by its average: sha(A(i)) = sd N (i) (A)/mean N (i) (A). In our experiments, we chose a patch of diameter 5 voxels. We present in Fig. 4 a quantitative view over voxels inside the brain mask of these two iconic measures with respect to the number of images in the atlas. These results show that both methods output atlases that are very similar in terms of image quality. The sharpness remains equivalent for both methods independently of the number of images added. sdis n increases only between 60 and 70 images and is then stable, but the evolution and values are the same for both atlas construction methods.

Conclusion

We have presented an online atlasing method that allows the incorporation of new images into an existing atlas without having to restart the atlasing process from the beginning. Based on an iterative centroid process, this approach only necessitates one registration per new image. We observed that the divergence between our method and direct atlasing using the same tools is localized in cortical areas and tends to grow but at a slow and decreasing pace. Also, the obtained atlases from both approaches have shown no differences in terms of image quality. The trend being at large, growing databases, the proposed online atlasing method offers an interesting tool to update an atlas at reasonable computational cost as new images arrive. Finally, we derived a construction method up to either an affine or a rigid transformation. This makes the method especially eligible for longitudinal atlasing as a future work, simply by adapting the weight of each image in Eq. 4.

Fig. 1 .

 1 Fig. 1. Shift of the centroid as new images are incorporated.

  ) ; 3-we assign θ k+1,k+1 = T k k+1 k+1 ; and finally 4-we build A k+1 = 1 k+1 k+1 j=1 I j • L j • θ j,k+1 .

3 d=1

 3 T j,d -Tj,d 2 where d indexes the three components in space.

Fig. 2 .Fig. 3 .

 23 Fig. 2. From left to right atlases with 60, 70, 80 and 100 images. Top: intensity differences between iterative centroid and direct atlases. Bottom: atlases from iterative centroid method overlaid with divergence measure δ.

Fig. 4 .

 4 Fig. 4. Boxplots and violin plots of inter-subjects normalized standard deviations sdisn (a) and sharpness sha (b) for atlases (restrained to the union of the masks) with 60, 70, 80, 90 and 100 images. Blue: iterative centroid. Red: direct atlasing.

Anima: Open source software for medical image processing from the Empenn team. https://github.com/Inria-Visages/Anima-Public -RRID:SCR 017017

Anima-Scripts: Open source scripts using Anima software for medical image processing from the Empenn team. https://github.com/Inria-Visages/ Anima-Scripts-Public -RRID:SCR 017072


Lauris Evariste

Antoine Mottier

Laura Lagier

Stéphanie Cadarsi

Maialen Barret

Cyril Sarrieu

Brigitte Soula

Florence Mouchet

• Emmanuel Flahaut

Eric Pinelli

Laury Gauthier

Assessment of graphene oxide ecotoxicity at several trophic levels using aquatic microcosms

Extensive development of new applications using graphene based materials such as graphene oxide (GO) increases its potential release and occurrence into aquatic environments, raising the question of its biological and ecological risks. As standardized single species based assays fail to highlight t(l)(icological pathways implying interactions between organisms, the use of micro/mesocosms appears as a good solution to fill the lack of environmental realism inherent to these tests. In this work, experiments were achieved using microcosm systems to expose a reconstituted food chain to GO at environmentally relevant concentrations (0.05 and 0.1 mgL 1 ). The trophic chain was composed of a consortium of algae and bacteria as primary producers, chironomid larvae as primary consumers and decomposers while larvae of the amphibian Pleurode/es wa/tii constituted the secondary consumers. Monitoring of multiple ecotoxicological and ecological endpoints allowed to observe changes in bacterial communities while no toxic effects were noticed in chironomids. However, chironomids feeding behaviour changed as a consequence of GO contamination, leading to an increase in leaf litter consumption. Genotoxic effects were noticed in Pleurode/es larvae. This study highlights the importance of using such experimental systems to better encompass the ecotoxic potential of GO through the determination of t(l)(icological routes and consequences on ecosystem's functioning.

Carbon based nanopartides (CBNs) possess unique properties (high sunaoe area, electrical and thermal conductivity, mechanical strength and optical transmittance) triggering tremendous scientific expectations for development of a wide range of industrial applica tians [START_REF] Tahir | The application of graphene and its derivatives to energy conversion, storage, and environmental and biosensing devices[END_REF][START_REF] Liu | Graphene and graphene oxide as new nanocarriers for drug delivery applications[END_REF][START_REF] Nguyen | Promising applications of graphene and graphenebased nanostructures[END_REF][START_REF] Perreault | Environmental applications of graphene-based nanomaterials[END_REF][START_REF] Scida | Recent applications of carbon-based nanomaterials in analytical chemistry: critical review[END_REF]. Graphene oxide (GO) is a chemically oxidized fonn of graphene that consists of a single a tom thick two dimensional sheet of carbon atoms containing epoxide, hydroxyl and carboxyl groups. Due toits high oxygen content, GO is a hydrophilic nanomaterial that is stable when dispersed in aqueous media [START_REF] Compton | Graphene oxide, highly reduced graphene oxide, and graphene: versatile building blocks for carbon-based materials[END_REF]. As other engineered nanoparticles (ENPs), GO is likely to be released into the ecosystems at any stage of its life cycle, from the production, the uses and the waste treatments of the substance and may cause environmental and health issues [START_REF] Arvidsson | Review of potential environmental and health risks of the nanomaterial graphene[END_REF][START_REF] Fadeel | Safety assessment of graphene-based materials: focus on human health and the environment[END_REF][START_REF] Keller | Global life cycle releases of engineered nanomaterials[END_REF][START_REF] Mottier | Environmental impact of engineered carbon nanoparticles: from releases to effects on the aquatic biota[END_REF]. Due to its hydrophilic properties associated to high surface area and chemical functions, GO could potentially react with many components of the environment, modifying its ecological risk toward aquatic ecosystems [START_REF] Zhao | Graphene in the aquatic environment: adsorption, dispersion, toxicity and transformation[END_REF].

Numerous toxicological studies of graphene based nana materials (GBNs) were penormed on mammalian cell lines or human related biological models [START_REF] Fadeel | Safety assessment of graphene-based materials: focus on human health and the environment[END_REF][START_REF] Ou | Toxicity of graphene-family nanoparticles: a general review of the origins and mechanisms[END_REF][START_REF] Seabra | Nanotoxicity of graphene and graphene oxide[END_REF][START_REF] Singh | Toxicity of graphene and its nanocomposites to human cell lines-the present scenario[END_REF] while potential effects on environmentally relevant models were not suffiàently investi gated (15,[START_REF] Montagner | Ecotoxicological effects of graphene-based materials[END_REF]. In addition, existing studies used assays performed through exposure of single speàes to CBNs are essential to un derstand toxicological mechanisms, these tests are poorly repre sentative of environmental conditions [START_REF] Bour | Environmentally relevant approaches to assess nanoparticles ecotoxicity: a review[END_REF]. For this purpose, the use of more complex systems such as micro mesocosms was shown to be effiàent to assess the ecotoxic potential of multiple contaminants under environmentally relevant conditions [START_REF] Bour | Environmentally relevant approaches to assess nanoparticles ecotoxicity: a review[END_REF][START_REF] Garaud | Integrated assessment of ceria nanoparticle impacts on the freshwater bivalve Dreissena polymorpha[END_REF][START_REF] Hommen | A microcosm study to support aquatic risk assessment of nickel: community-level effects and comparison with bioavailability-normalized species sensitivity distributions: aquatic microcosm study with Ni[END_REF][START_REF] Müller | Long-term effects of the fungicide pyrimethanil on aquatic primary producers in macrophyte-dominated outdoor mesocosms in two European ecoregions[END_REF]. These experimental systems allow measurement of bath Thus, performing this type of study is becoming crucial for a better understanding of potential environmental issues implying any contaminant and including GO. Previous works mainly focused on the ecotoxicological effects of GO on organisms from lower trophic levels such as bacteria or algae exposed in single culture. Antibacterial activity of GO was identified [24e26] and among multiple GBNs, GO was shown to exert the highest cytotoxicity toward the bacterial model Escherichia coli [START_REF] Liu | Antibacterial activity of graphite, graphite oxide, graphene oxide, and reduced graphene oxide: membrane and oxidative stress[END_REF]. Concerning adverse effects toward algae, it was demonstrated that GO exposure led to Chlorella sp. growth inhibition and produced membrane damages through mechanical actions as well as oxidative stress induction [START_REF] Hu | Interactions between graphene oxide and plant cells: regulation of cell morphology, uptake, organelle damage, oxidative effects and metabolic disorders[END_REF]29]. Studies carried out on organisms from higher trophic levels (i.e. primary and secondary consumers) are still scarce. It was indicated that GO induced behavioral disturbances of crustacean larvae (Amphibalanus amphitrite and Artemia salina) and was strongly accumulated in digestive tract [START_REF] Mesaric | High surface adsorption properties of carbonbased nanomaterials are responsible for mortality, swimming inhibition, and biochemical responses in Artemia salina larvae[END_REF][START_REF] Mesari C | Effects of nano carbon black and single-layer graphene oxide on settlement, survival and swimming behaviour of Amphibalanus amphitrite larvae[END_REF]. In amphibians, similar accumulation of GO was observed in gut of Xenopus laevis tadpoles and growth was inhibited at high concentrations [START_REF] Lagier | Surface area of carbon-based nanoparticles prevails on dispersion for growth inhibition in amphibians[END_REF]. Observation of toxicological effects induced by GO on low trophic level organisms raises the question of possible larger scale conse quences through trophic network and on ecosystem functioning.

The aim of this study was to investigate the toxicity of GO under realistic conditions to go further in the understanding of its impact on organisms occupying multiple trophic levels and to determine potential consequences on ecosystem's functioning. For this pur pose, microcosm systems were used to expose organisms interacting through a reconstituted trophic chain. The biological models retained for the study were chosen for being key species represen tative of a simplified trophic chain from pond ecosystems [START_REF] Rimet | Life-forms, cell-sizes and ecological guilds of diatoms in European rivers[END_REF][START_REF] Adler | Ecological and societal services of aquatic Diptera[END_REF]. Although they have significant ecological functions, number of hu man activities are threatening these ecosystems as well as the biodiversity they host, especially endangering species such as am phibians [START_REF] Oertli | Conservation and monitoring of pond biodiversity: introduction[END_REF]. For this purpose, the benthic diatom Nitzschia palea and a bacterial consortium were settled to form a biofilm providing food for chironomid larvae (Chironomus riparius) as primary con sumers and decomposers. Secondary consumers were constituted by the carnivorous larvae of the amphibian, Pleurodeles waltl, the Spanish newt, that feed on macro invertebrates. As microbial com munities are at the very basis of every ecosystem functioning and due to their implication on global nutrient cycling [START_REF] Barton | Environmental Microbiology and Microbial Ecology[END_REF], changes in microorganism community composition were monitored during the experiment. In terms of toxicological endpoints, survival and growth parameters of the consumers were controlled as well as morpho logical deformities in chironomids and count of micronucleated erythrocytes in Pleurodeles [START_REF] Gauthier | Chapter 15 genotoxicity in urodele Amphibians Pleurodeles waltl and Ambystoma mexicanum (Lissamphibia, Caudata) exposed to freshwater pollutants: a historical view, in: Ecotoxicol. Genotoxicol. Non-tradit[END_REF][START_REF] Salmelin | Inconsistency in the analysis of morphological deformities in chironomidae (Insecta: Diptera) larvae[END_REF]. Analysis of chironomids mouth part deformities constitutes a marker of teratogenesis that is widely used in pollution monitoring as well as micronuclei induction that is a good indicator of genotoxicity, integrating aneugenic and clasto genic effects that is used as endpoint in ISO 21427e1 procedure. Both teratogenic and genotoxic biomarkers are considered as highly predictable ecological endpoints at the population, community and ecosystem level [START_REF] Anderson | Genetic and molecular ecotoxicology: a research framework[END_REF]. The process of leaf litter decomposition ensured by chironomids and bacteria was also investigated, as it constitutes an interesting functional ecological marker integrative of lower trophic level activities [START_REF] Bour | Impact of CeO 2 nanoparticles on the functions of freshwater ecosystems: a microcosm study[END_REF][START_REF] Pascoal | Role of fungi, bacteria, and invertebrates in leaf litter breakdown in a polluted river[END_REF].

Material and method

Graphene oxide nanoparticles

Graphene oxide nanoparticles were produced from twisted ribbon shaped carbon nanofibers (GANF®) processed by Hummer method [START_REF] Lobato | Large-scale conversion of helical-ribbon carbon nanofibers to a variety of graphene-related materials[END_REF] and supplied by Antolin Group. Full characterization of the tested material was already presented in earlier work [START_REF] Evariste | Thermal reduction of graphene oxide mitigates its in vivo genotoxicity toward Xenopus laevis tadpoles[END_REF] and data are summarized in Table 1. Both experiments were performed within few weeks. GO was stored as dry powder in the dark and dispersions were prepared extemporaneously in order to avoid any possible change of material characteristics.

Stability of graphene oxide suspension in exposure media

The stability of the suspension of GO was evaluated using a Turbiscan LAB (Formulation) equipment, at room temperature. Suspensions were prepared at 10 mg L À1 in deionised water (EDI) or Volvic water immediately before analysis. This concentration was assessed as low GO concentrations low could not be detected by TurbiscanTM. The suspension was bath sonicated for 2 min before insertion into the measurement chamber. Sedimentation was monitored for 24 h with 1 scan per minute by measuring during each scan the transmission of the suspension vs the height in the vial. In case of sedimentation, it is thus expected that the transmission will increase at the heights where particles have settled down, while it should decrease where sedimented particles are present.

Microcosm exposure

Nine microcosms (glass tanks 60 Â 30 Â 30 cm, L x W x H) were filled with 6780 g of reconstituted sediment: (96.6% of silica sand, 2.4% of kaolinite and 1% of CaCO 3 ) and 40 L of the commercial natural spring water Volvic® (Ca 2þ : 12mgL À1 , Na þ : 12 mg L À1 Mg 2þ : 8 mg L À1 , K þ : 6 mg L À1 , Cl À : 1 mg L À1 and Si: 32 mg L À1 ) as used in previous microcosm studies [21,43e45]. In order to maintain correct oxygenation, each tank was equipped with a recirculation system. Light was provided by fluorescent tubes (JBL solar ultra, Natur 9000 K LT 24 WT5 HQ) and photope riod was set on 12:12 light/night. Temperature (21.5 ± 0.9 C), conductivity (234.3 ± 9.1 mS cm À1 ), pH (7.9 ± 0.4), redox potential (306.6 ± 17.3 mV) and dissolved oxygen (9.1 ± 1.2 mg L À1 ) were monitored continuously (Ponsel Odeon open X probes kit). Nitro gen products (NO 3 À , NO 2 À and NH 4 þ ) were monitored twice a week using a HI83203 photometer (Hanna instruments) and dissolved organic carbon was measured at each sampling time by infrared detection of CO 2 produced by catalytic oxidation at 680 C (Table S1). Three different exposure conditions were tested: control (nat ural spring water), graphene oxide at the final concentration of 0.05 mg L À1 and graphene oxide at the final concentration of 0.1 mg L À1 . Each condition was tested in triplicate with random assignation of the microcosms. Contamination of the microcosms by graphene oxide nanoparticles was performed sequentially through the addition of 1/12 e of the final concentration over 10 days to reach the final concentration indicated by the condition name (0.05 mg L À1 and 0.1 mg L À1 ) (Fig. 1).

The progress of the experiment and the order in which the different organisms were introduced are shown in Fig. 1. Senescent alder leaves (Alnus glutinosa Gaertn.) (3.2 g dry mass/mesocosm) were added the first week for measurement of decomposition rate and microorganisms (microbial consortium 1 Â 10 6 cells/ml and Nitzschia palea 2.05 Â 10 4 cells/ml) were allowed to settle and develop during 3 weeks without any contamination. This duration was necessary for the development of the primary producers leading to the formation of a biofilm covering the whole sediment [START_REF] Tella | Chronic dosing of a simulated pond ecosystem in indoor aquatic mesocosms: fate and transport of CeO 2 nanoparticles[END_REF] and equilibration of nitrogen products in the microcosms by bacteria. The first samplings of the biofilm and water were per formed to analyse microorganism communities prior to contami nation and addition of primary and secondary consumers (T0). Other organisms, bred in the laboratory facilities, were added sequentially to the microcosms (Fig. 1). 800 Chironomus riparius individuals aged of 72 h and 15 Pleurodeles waltii individuals at stage 53 of development [START_REF] Gallien | Table chronologique du deyeloppement chez Pleurodeles waltlii[END_REF] were added during the fourth week to allow an exposure duration of 13 and 10 days respectively. These durations were chosen in order to ensure that the surviving C. riparius would develop until the stage 4 (required condition for teratogenicity analysis). Moreover, chironomids were added prior to the newts to grow and start to accumulate GO in order to ensure trophic transfer towards the secondary consumers. For Pleurodeles, exposure duration of 10 days is necessary for genotoxic assessment at low dose [START_REF] Jaylet | A new micronucleus test using peripheral blood erythrocytes of the newt Pleurodeles waltl to detect mutagens in fresh-water pollution[END_REF]. A second sampling of the microorganisms was performed when half of the final GO dose was reached (T1) and at the end of the exposure when the full dose was reached (T2). Xenopus laevis (n 40) at stage 50 of development [START_REF] Nieuwkoop | Normal table of Xenopus laevis (Daudin). A systematical and chronological survey of the development from the fertilized egg till the end of metamorphosis[END_REF] were added during the week 5 as a food complement for Pleurodeles larvae. Chironomids, Pleurodeles and alder leaves were also sampled at T2, constituting the end of the experiment. Water and sediment were sampled at 3 different times during the experiment (T0, T1 and T2). At each sampling time, 100 mL of water and 10 mL of sediment were taken. Water samples were filtered on a 5 mm filter (Whatman® Nuclepore™ Track Etched Membranes) in order to recover microorganisms from the water column. Filters and sediment samples were snap frozen in liquid nitrogen and stored at 80 C before processing for analyses. Total DNA was extracted from water samples and sediment subsamples (644.6 ± 37.2 mg) using the Mobio® power soil DNA isolation kit following manufacturer's instructions. Additionally, negative extraction controls were performed to ensure the absence of DNA contamination during the process. The quantity and quality of DNA extracts was analysed using a NanoDrop 2000 UV spectropho tometer (Thermo Scientific).

The V4eV5 region of archeae and bacteria 16S rRNA gene was targeted using PCR1_515F (5 0 GTGYCAGCMGCCGCGGTA 3 0 ) and PCR1_928R (5 0 CCCCGYCAATTCMTTTRAGT 3 0 ) primers set [START_REF] Wang | Conservative fragments in bacterial 16S rRNA genes and primer design for 16S Ribosomal DNA amplicons in metagenomic studies[END_REF]. PCRs were performed on a Gene AmpTM PCR system 9700 ther mocycler (Applied Biosystems, Forster City, CA, United States) in a final volume of 50 mL containing: 37.5 mL of PCR water, 5 mL of 10X PCR buffer, 2 mL of extracted DNA, 2 mL of both primer (10 mM), 1 mL of dNTP (2.5 mM) and 0.5 mL of Taq DNA polymerase (5 U/mL e Sigma Aldrich). The following PCR protocol was applied: 94 C for 120 s, 30 cycles of 94 C for 60 s, 65 C for 40 s, 72 C 30 s and 72 C for 10 min.

Sequencing of amplicons from 16S rRNA genes was performed using Illumina MiSeq technology (2 Â 250 pb) by the Get_PlaGe platform (Genotoul, Toulouse, France). Bioinformatic analysis was performed using FROGS (Find Rapidly Operational Taxonomic Units (OTU) Galaxy Solution) pipeline on Galaxy [START_REF] Escudi | FROGS: find, rapidly, OTUs with Galaxy solution[END_REF]. Briefly, sequences with mismatches in the primers were excluded and PCR primers were trimmed. Reads were clustered into OTUs using the Swarm clustering method [START_REF] Mah E | Swarm: robust and fast clustering method for amplicon-based studies[END_REF]. Chimera were removed and filters were applied to keep OTUs present in at least 3 samples and representing at least 0.005% of all sequences [START_REF] Bokulich | Quality-filtering vastly improves diversity estimates from Illumina amplicon sequencing[END_REF]. 478 OTUs were assigned at different taxonomic levels (from Kingdom to species) using RDP classifier and NCBI Blast þ on Silva 132 database (pintail 80) [START_REF] Quast | The SILVA ribosomal RNA gene database project: improved data processing and web-based tools[END_REF].

Chironomids

At the end of the experiment, mortality, growth and teratoge nicity were assessed in chironomid larvae that were not predated by Pleurodeles. A careful sorting of the sediment allowed to count the surviving individuals which were then measured using the image J software (total length and width of the cephalic capsule).

Fig. 1. Scheme of the time-progress of the microcosm experiment. Bacterial consortium and diatoms were allowed to settle and develop to colonize the sediment and alder leaves during 3 weeks. Prior to the addition of other organisms and GO, a first sampling of the primary compartment was performed (T0). GO was then added sequentially in the microcosms over 2 weeks to reach 0.05 mg L 1 and 0.1 mg L 1 as final concentration. Pleurodeles (P. waltl) and Chironomids (C. riparius) were added subsequently during the fourth week while a second sampling of microbial consortium was performed at the end of that same week (T1). All organisms were sampled at the beginning of the 6th week to perform analysis of the different endpoints (T2). (A colour version of this figure can be viewed online.)

The cephalic capsules were recovered and discolored in potassium hydroxide (15%, 95 C, 15 min), then mounted on microscope slides with quick hardening medium (Sigma). Teratogenecity was assessed by counting the number of deformities in mentum. Ac cording to Salmelin et al. [START_REF] Salmelin | Inconsistency in the analysis of morphological deformities in chironomidae (Insecta: Diptera) larvae[END_REF], only absolute deformities such as missing and extra teeth and Kohn gaps were taken into account.

Leaf litter degradation

Alder leaves degradation was estimated by the total surface of degradation after microcosm exposure. Leaves were scanned at the beginning and at the end of experiment using a flatbed scanner with 600 DPI resolution. Leaves pictures were processed with image J software: a 8 bit transformation and thresholding allowed us to quantify the degraded surfaces.

Amphibians

No endpoint was measured in X. laevis tadpoles that were all predated during the experiment by Pleurodeles larvae. For Pleuro deles larvae, mortality, growth and genotoxicity were studied after 10 days of exposure (n 15 for each endpoint per microcosm). Mortality was checked daily throughout the exposure. Growth was studied by measurements of size (photo and image analysis with Image J software) at the beginning and end of the experiment. Genotoxicity was assessed by enumeration of micronucleated erythrocytes (MNE) in blood smears after fixation (methanol 99.9%) and staining (groat hematoxylin). The number of cells containing one or more micronuclei was counted among a total of 1000 cells.

Statistical analysis

Statistical analysis was performed using Minitab 16 Statistical software. Data of growth (P. waltl and C. riparius), development and survival chironomids (growth, development and survival) as well as leaf litter degradation endpoints were analysed using one way ANOVA when assumptions of normality and homogeneity of vari ance were met. As it was not the case for genotoxicity data, non parametric Kruskal Wallis followed by Dunn's multiple compari son test was performed. For metagenomics analysis, data manip ulations, OTUs counts, alpha diversity indexes and Unifrac Distances calculations as well as multidimensional scaling (MDS) plot were carried out using "Phyloseq" R package [START_REF] Mcmurdie | phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data[END_REF]. Differential abundance of bacterial genera between exposed conditions compared to the control group was examined using "Deseq2" R package [START_REF] Love | Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[END_REF]. PERMANOVA was performed using Adonis function from the "vegan" R package [START_REF] Oksanen | Package 'vegan[END_REF]. The relationships between the biological responses measured were explored by multivariate Principal Component Analysis (PCA) and correlation analyses using "ade4 R" package [START_REF] Dray | The ade4 package: implementing the duality diagram for ecologists[END_REF].

Results

Behaviour of GO in exposure media

The variation of the transmission data (compared to t0) vs time for EDI and Volvic are presented in Fig. 2. While the dispersion is stable over 24 h in EDI medium, a clear homogeneous sedimenta tion is evidenced in the case of the Volvic water after 5 h (Fig. S1).

Effects on lower trophic level

Bacterial communities in water compartment were different from those in sediment compartment regardless of time and GO concentration, as revealed by MDS and PERMANOVA analysis using unweighted UniFrac distances (Fig. S2) (F 32.554; r 2 0.385; p 0.001). Thus, the two compartments were studied separately to analyse the effects occurring on compartment associated bacterial communities. At T0, prior to contamination, community structures in the water were similar between the condition related micro cosms (PERMANOVA: F 1.193; r 2 0.284; p 0.132), and the same was observed in sediments (F 1.143; r 2 0.276; p 0.252).

Effects on microorganisms from the water column

When analysing the alpha diversity, no difference was measured for the observed diversity, Chao1 and Shannon indexes in the water column, during the whole experiment (Fig. S3A). However, exposure duration as well as GO concentration were shown to significantly affect the microbial community structure (PERMANOVA: Contami nant: F 1.657; r 2 0.104; p 0.010; Exposure duration: F 3.428; r 2 0.216; p 0.001). MDS analysis indicated a higher contribution of the exposure duration (first axis) toward bacterial community composition compared to GO concentration (second axis) (Fig. S3B). In the water column at T0, over 95% of overall community was represented by phyla Proteobacteria and Bacteroidetes with a rela tive abundance of 58.8 ± 21.6% and 37.5 ± 22.2% respectively (Fig. 3A). The relative abundance of none of the phyla analysed was shown to be impacted by GO contamination (Table S2). Focusing on families from the two main phyla found in the water column (Fig. S4), almost no effects were observed on relative abundances (Table S2). Analysis of the genera differential abundance in the GO conditions compared to the control group during the overall expo sure led to a significant decrease of 3 genera from the phylum Pro teobacteria (Legionella, Sphingobium and Reyranella) when exposed to 0.05 mg L À1 of GO, while three other genera from this phylum were decreased in the 0.1 mg L À1 exposure condition (Pseudomino bacter, Sphingorhabdus and Devosia) (Fig. S5A).

Effects on microorganisms in sediment

A significant increase in the biofilm observed richness and Chao1 was noticed between T0 and T1 without being influenced by GO concentration and remained stable until T2 while values of the Shannon index were not significantly influenced (Fig. 4A). Exposure duration and GO concentration significantly affected microbial community structure (PERMANOVA: Contaminant: F 2.862; r 2 0.161; p 0.002, Exposure duration: F 4.875; r 2 0.274; p 0.001) with a greater contribution of the exposure duration over the contaminant concentration (Fig. 4B). At T0, prior to contamination by GO, three main phyla constituted the bacterial compartment of the biofilm. Phyla Proteobacteria, Bacteroidetes and Planctomycetes relative abundances were of 29.1 ± 15.1%, 10.6 ± 4.2% and 3.8 ± 2% respectively. The chloroplastic 16S rRNA gene from diatoms was affiliated to Ochrophyta phylum. This phylum represented 48.7 ± 20.4% of the overall biofilm organisms (Fig. 3B). The relative abundance of the two main bacterial phyla and diatoms was not affected by exposure to GO, while the phyla Planctomycetes, Armatimonadetes and WPS 2 significantly increased in the presence of GO at 0.1 mg L À1 compared to the control group (Table S3). The effects occurring at the family scale in the three main bacterial phyla from the sediment were analysed (see supplementary data and Fig. S6). Analysis of the OTUs differ ential abundances in the biofilm indicated that among 517 OTUs, exposure to GO at 0.05 mg L À1 led to significant changes in the relative abundance of 15 taxa, among which 5 were found in a higher abundance and 10 in a lower abundance (Fig. S5B). At 0.1 mg L À1 , 20 taxa were differentially abundant compared to the control with an equal distribution between over and under repre sented taxa. The relative abundance of genera Gemmata, Azospir illum and Flavobacterium were shown to decrease in a similar manner in the two GO conditions compared to the control group while genera Planctopirus and Armatimonas increased in the same order of manitude.

Effects on chironomids

At the end of the experiment, no significant difference was noticed on the survival rate of C. riparius whatever the exposure condition (ANOVA, p 0.976). A mean of 8.02 ± 4.7% of the overall chironomids were remaining for all conditions. Larval growth and determination of development stage measured from length of the organisms and head width measurement respectively, were not significantly affected by GO exposure (ANOVA, p 0.280 and 0.860 respectively). Considering the whole organisms, 60.6% of the larvae reached the stage 4 at the end of the experiment and a mean size of 8.6 ± 2.1 mm was measured. Analysis of teratogenicity in stage 4 larvae indicated no effect of GO exposure (ANOVA, p 0.913).

Effects on leaf litter degradation

At the end of the exposure, a significant increase in leaf degra dation was noticed after exposure to 0.05 mg L À1 of GO compared to the control group while only a trend in increase of degradation was observed at 0.1 mg L À1 (ANOVA p 0.027 followed by Tuckey) (Fig. 5).

Effects on pleurodeles larvae

No mortality in Pleurodeles larvae was observed during the whole experiment and larval growth was not affected by GO exposure (ANOVA, p 0.723). At the end of the experiment, a sig nificant increase in micronucleated erythrocytes was observed in larvae exposed to 0.1 mg L À1 (Kruskal Wallis: p 0.0013 followed by Dunn's test) (Fig. 6).

Principal component analysis

In the principal component analysis of the responses measured at the different trophic levels, PC1 and PC2 explained 38,6% and 24.3% of the total variance respectively (Fig. 7). The PCA showed a better segregation between the controls microcosms and those contaminated with GO along the PC2 axis based on biological re sponses. The parameters related to chironomids development and survival as well as microbial diversity were negatively associated with PC1 while teratogenesis was positively associated with PC1. The leaf litter degradation and genotoxicity in pleurodeles were mainly associated with the conditions of GO exposure along PC2.

Discussion

To our knowledge, this study is the first one to investigate the effect of GO under environmentally relevant conditions using a reconstituted trophic chain through the use of microcosms. Because of the limited available techniques to quantify GO in complex matrix and changes of materials characteristics over time at environmentally realistic concentrations [START_REF] Goodwin | Detection and quantification of Graphene-family nanomaterials in the environment[END_REF], no relevant quantification analysis could be performed in our study. Thus, the hypothesis of material distribution in our system and toxicological mechanisms involved are based only on Turbiscan data and bio logical responses measured.

Behaviour of GO in exposure media

The fate, transport, and bioavailability of GO are mainly deter mined by their behaviour in aquatic ecosystems including aggre gation deposition that can be influenced by multiple factors [START_REF] Ren | Graphene analogues in aquatic environments and porous media: dispersion, aggregation, deposition and transformation[END_REF]. On the one hand, interactions with organic matter of the experi mental system would increase the colloidal stability of the GO [START_REF] Zhao | Graphene in the aquatic environment: adsorption, dispersion, toxicity and transformation[END_REF][START_REF] Castro | Nanoecotoxicity assessment of graphene oxide and its relationship with humic acid: nanoecotoxicity of graphene oxide and humic acid[END_REF][START_REF] Chowdhury | Colloidal properties and stability of graphene oxide nanomaterials in the aquatic environment[END_REF], while on the other hand, the presence of ions in the water such as Ca 2þ or Na þ could adsorb on the negatively charged functional groups of GO, leading to reduced surface charge and reduced dispersion stability, leading to its deposition [START_REF] Chowdhury | Aggregation and stability of reduced graphene oxide: complex roles of divalent cations, pH, and natural organic matter[END_REF]. Ac cording to the obtained results it seems that the interactions with ions of the exposure media occurring under the experimental conditions are responsible for the sedimentation of the GO. Such sedimentation of the material would decrease its bioavailability for organisms from the water column while it will increase for or ganisms living at the sediment water interface.

Effects on microorganisms

Studying the effects of GO on microorganisms is essential as they constitute the first link of the food chain. Literature is fur nished of data demonstrating the antibacterial capacities of GO [24e26,59e61]. However, few studies focused on the effects generated by CBNs on complex bacterial communities from aquatic ecosystems and the data available were mainly obtained on bac terial consortium from soils [START_REF] Du | Graphene oxide regulates the bacterial community and exhibits property changes in soil[END_REF][START_REF] Xiong | Implication of graphene oxide in Cd-contaminated soil: a case study of bacterial communities[END_REF] or activated sludge [START_REF] Ahmed | Investigation of acute effects of graphene oxide on wastewater microbial community: a case study[END_REF][START_REF] Guo | Effect of graphene oxide on the bioactivities of nitrifying and denitrifying bacteria in aerobic granular sludge[END_REF]. In the present study, complex changes were observed in the bacterial compartment over the whole experiment. It is interesting to point out that GO exposure induced only minor differences in bacterial communities in the water column whereas the biofilm commu nities were more impacted. The difference in magnitude of GO impact between water and sediment communities associated to Turbiscan analysis confirm that the GO is more bioavailable for microbial communities from the sediment rather than for free living microorganisms from the water column. The presence of a complex biofilm containing diatoms N.palea and bacteria also contribute to this change of bioavailability. Indeed, these microor ganisms were shown to stabilize CBNs in the biofilm through secretion of extracellular polymeric substances (EPS) [70e72].

Thus, the trapping of GO may decrease its bioavailability for bac teria from the water column while it would increase the occurrence of direct contact with bacteria in the biofilm (i.e. and other benthic organisms). This direct contact was evidenced as the main mech anism of cytoxicity toward bacteria through impairment of cell membrane integrity [START_REF] Liu | Antibacterial activity of graphite, graphite oxide, graphene oxide, and reduced graphene oxide: membrane and oxidative stress[END_REF][START_REF] Combarros | Toxicity of graphene oxide on growth and metabolism of Pseudomonas putida[END_REF][START_REF] Dizaj | Antimicrobial activity of carbon-based nanoparticles[END_REF].

Our data show that the biofilm compartment is more dynamic compared to free living bacteria as indicated by the time effect on diversity indexes from the control group. However, as there are no other organisms in the microcosms at T0, introduction of chiron omids and newts between T0 and T1 would contribute mainly to the biofilm's dynamics through the increase in bacterial alpha di versity via the production of faeces, for example. In spite of this temporal dynamics inherent to microcosm experimental design as well as to this form of bacterial life [START_REF] Sutherland | The biofilm matrix an immobilized but dynamic microbial environment[END_REF], beta diversity measure ment in pelagic as well as in benthic communities evidenced diverging phylogenetic trajectories in presence of GO compared to the control group. In this case, the stronger the contamination pressure is, the more communities are diverging. As the multiple species in biofilms may possess different sensitivity to GO, the decrease of affected taxa make available new ecological niches benefiting to less sensitive and more opportunistic species leading to different equilibrium in microbial communities. Thus, the phyla Planctomycetes, Armatimonadetes and WPS 2 appear tolerant to GO as their relative abundances increased in the exposed condi tions. A similar trend was previously evidenced for Planctomycetes from soil bacterial communities exposed to GO [START_REF] Forstner | Effects of Graphene Oxide and Graphite on Soil Bacterial and Fungal Diversity[END_REF] or silver nanoparticles [START_REF] Wang | Effects of silver nanoparticles on soil microbial communities and bacterial nitrification in suburban vegetable soils[END_REF]. These changes would contribute to the di vergences in the evolution of the biofilm composition between the tested conditions. However, as a dynamic system, it is possible that the direct toxicity of GO will be attenuated with biofilm maturation. Indeed, it was previously observed in bacterial biofilm using E. coli or S. aureus models that after 48 h of exposure, corresponding to the mature phase of biofilm, cell growth inhibition and induction of ROS production by rGO was suppressed, even at concentrations up to 100 mg L À1 [START_REF] Guo | Toxicity and transformation of graphene oxide and reduced graphene oxide in bacteria biofilm[END_REF]. This study also indicated that bacteria were able to oxidize rGO into GO. Transformations of the tested material potentially occurring under our experimental conditions are remaining to be assessed when technological barriers will be overcome.

Effects on chironomid larvae and consequences for leaf litter degradation

No apparent toxicity on growth and teratogenicity could be highlighted on chironomids. The low amount of larvae remaining at the end of the exposure could be only related to predation by Pleurodeles larvae that is more likely to occur within the last days of the experiment when they reach larger sizes [START_REF] Bour | Toxicity of CeO 2 nanoparticles on a freshwater experimental trophic chain: a study in environmentally relevant conditions through the use of mesocosms[END_REF][START_REF] Green | Shorebird predation affects density, biomass, and size distribution of benthic chironomids in salt pans: an exclosure experiment[END_REF]. It is also possible that the apparent no effects on chironomids may result from a preferential predation by pleurodeles larvae of affected chironomids. However, the low sensitivity of chironomids toward CBNs was previously indicated for fullerene [START_REF] Waissi-Leinonen | Toxicity of fullerene (C60) to sediment-dwelling invertebrate Chironomus riparius larvae[END_REF], multi walled carbon nanotubes [START_REF] Martínez-Paz | Effects at molecular level of multi-walled carbon nanotubes (MWCNT) in Chironomus riparius (DIPTERA) aquatic larvae[END_REF] as well as GO [START_REF] Castro | Nanoecotoxicity assessment of graphene oxide and its relationship with humic acid: nanoecotoxicity of graphene oxide and humic acid[END_REF]. In the later study, no effects were observed on chironomids after 9 days of exposure to GO at concentrations up to 100 mg L À1 , which is consistent with the present results and corroborate the predation hypothesis without any pressure of selection.

Chironomids are shredders, deposit feeders that preferentially feed on biofilm rather than non living organic matter [START_REF] Widenfalk | Sediment microbes and biofilms increase the bioavailability of chlorpyrifos in Chironomus riparius (Chironomidae, Diptera), Ecotoxicol[END_REF]. Despite the known implication of bacteria in this ecological process, chi ronomids constitute the main decomposers of organic matter in such time duration experiment [START_REF] Bour | Impact of CeO 2 nanoparticles on the functions of freshwater ecosystems: a microcosm study[END_REF]. Despite the pressure of pre dation on chironomids, the alder leaves decomposition process is maintained. In this study, it appears that this process increases in presence of GO while opposite results were observed in the case of contamination by other nanoparticles [START_REF] Bour | Impact of CeO 2 nanoparticles on the functions of freshwater ecosystems: a microcosm study[END_REF][START_REF] Pradhan | Can metal nanoparticles be a threat to microbial decomposers of plant litter in streams?[END_REF]. It was previously observed through the comparison of deformities in the mouthparts of chironomids sampled in differentially contaminated sites that teratogenicity was associated with lower leaves degradation ac tivities [START_REF] Macdonald | Incidence of mentum deformities in midge larvae (Diptera:chironomidae) from Northern Nova Scotia, Canada[END_REF]. In the present study, no teratogenicity was measured and as PCA results indicated an independence of these variables, we can suggest that the increase of leaves consumption is associated to indirect effects. The presence of GO in the biofilm might influence their feeding behaviour through two possible ways. 1) Changes in bacterial communities and GO trapping might decrease the biofilm palatability. Indeed, it was demonstrated that food quality prevails on contaminant concentration for chironomids food source selec tion [START_REF] De Haas | Habitat selection by chironomid larvae: fast growth requires fast food[END_REF]. Thus, the loss of biofilm palatability and quality would lead to preferential feeding on alder leaves. 2) The presence of GO at low concentration may increase fungus development as it was demonstrated that GBMs could stimulate the growth of rot fungus at "low concentration" without influencing the decomposition ac tivity [START_REF] Xie | Toxicity of graphene oxide to white rot fungus Phanerochaete chrysosporium[END_REF][START_REF] Yang | Influence of reduced graphene oxide on the growth, structure and decomposition activity of whiterot fungus Phanerochaete chrysosporium[END_REF]. Thus, an increase in fungus biomass on leaves would improve their palatability and influence chironomids feeding behaviour [START_REF] Casas | The effect of diet quality on growth and development of recently hatched larvae of Chironomus gr. plumosus[END_REF][START_REF] Gessner | A perspective on leaf litter breakdown in streams[END_REF]. According to the PCA results, it seems that there is a correlation between the diversity of microorganism species from the biofilm and the development of chironomids that is consistent with the first hypothesis. It is interesting to note that in any case, this change of behaviour had no measurable consequence on chironomids development. Modifications of the feeding behaviour could also theoretically contribute to changes observed in bacterial community compositions. Indeed, diet changes could lead to modification of chironomids gut microbiota and modify microbial composition of the faeces rejected into the mesocosms.

These changes in organic matter decomposition could potentially lead to changes in carbon and nutrient cycling, influencing energy flows of the system [START_REF]Methods in Stream Ecology[END_REF] but it remains to be determined whether or not it is occurring under these conditions. Thus, litter decomposition appears as an interesting marker integrating direct and indirect toxicological effects occurring at the lower trophic levels.

Effects on pleurodeles larvae

The results obtained in this experiment indicated that GO is able to induce genotoxic effects in vivo at low concentration under realistic environmental conditions. Genotoxic potential of GO was also observed in other organisms such as mice after repeated or single injection of GO at concentrations ranging from 0.01 to 20 mg kg À1 [START_REF] Liu | Graphene oxide can induce in vitro and in vivo mutagenesis[END_REF][START_REF] El-Yamany | Graphene oxide nanosheets induced genotoxicity and pulmonary injury in mice[END_REF]. This was also observed in the amphibians Xenopus laevis [START_REF] Evariste | Thermal reduction of graphene oxide mitigates its in vivo genotoxicity toward Xenopus laevis tadpoles[END_REF], which is consistent with the present results. In the latter study case, tadpoles were exposed during 12 days to the same commercial GO as the one used in the current microcosm study. The genotoxic response measured in X. laevis after exposure to 0.1 mg L À1 was associated to oxidative stress and pro inflammatory response leading to an increase in circulating micronucleated erythrocytes [START_REF] Evariste | Thermal reduction of graphene oxide mitigates its in vivo genotoxicity toward Xenopus laevis tadpoles[END_REF]. Despite important biological differences between these two amphibian model species, it can be suggested that the genotoxic effects observed in Pleurodeles resulted from similar pathways. Under our experimental condi tions, Pleurodeles larvae were exposed to GO through two path ways simultaneously: direct exposure due to the presence of GO in the media as well as dietary exposure by feeding on contaminated chironomids. However, PCA analysis indicated that the marker of genotoxicity was independent from those associated to other trophic levels, which reinforces the hypothesis of effects associ ated to direct exposure over trophic effects. In the case of a direct exposure, as the contamination was performed sequentially and according to previous study indicating that at least 4 days of exposure to a genotoxic compound are necessary to induce the formation of micronuclei [START_REF] Jaylet | A new micronucleus test using peripheral blood erythrocytes of the newt Pleurodeles waltl to detect mutagens in fresh-water pollution[END_REF], we can suggest that a concentration of 0.8 mg L À1 of GO is able to induce genotoxic effects (concen tration reached 4 days before the sampling). It seems that the predation behaviour of the Pleurodeles was not influenced by GO exposure. This is supported by the high decrease in chironomid larvae as well as the absence of differences within the chirono mids survival between the tested conditions at the end of the microcosm exposure. However, to date, it is not possible to evaluate separately the contribution of each pathway in the gen otoxic response observed.

Conclusion

This microcosm study allowed to identify the effects of a contamination by GO towards multiple interacting organisms and consequences on ecosystem functioning. Analysis indicated that GO could induce toxicological effects in top predators as well as in microorganisms at environmentally relevant concentration [START_REF] Zhang | Molecular mechanisms of developmental toxicity induced by Graphene Oxide at predicted environmental concentrations[END_REF]. Genotoxicity was evidenced in top predators while bacterial com munities were modified after GO exposure. Communities from the sediment were shown to be more impacted than those from the water column due to sedimentation of the GO under more realistic environmental conditions. Even if no toxicity was observed for chironomids, indirect effects were highlighted and led to changes in the decomposition of organic matter in the system.

This experiment allowed to better encompass ecotoxicity of GO under realistic conditions and supports the need to develop the use of more complex systems including trophic chains in ecotoxicity assessment. Indeed, the use of these experimental systems may evidence routes of exposure that are poorly or not estimated in classical standardized tests based on single species assessments. For families from the phylum Bacteroidetes, the relative abundance AKYH767 and NS9 significantly increased under the exposure condition of 0.1 mg/L of GO (Table S2). Relative abundance of family Pseudomonadaceae from phylum Proteobacteria significantly decreased at both concentrations of GO (Table S2) while for the phylum Planctomycetes, bacteria from family Phycisphaeraceae, Rubinisphaeraceae and Schlesneriaceae increased in presence of GO (Table S2).

2. 4 .

 4 Endpoint measurement 2.4.1. Primary compartment: sampling, DNA extraction, PCR and sequencing

Fig. 2 .

 2 Fig. 2. Variation of the transmission vs time for suspensions of GO at 10 mg L 1 in deionised water (blue) or Volvic water (red). (A colour version of this figure can be viewed online.)

Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .

 3456 Fig. 3. Relative abundance of bacterial phyla from the water column (A) and in the biofilm of the sediment (B) in function of GO concentration (0, 0.05 and 0.1 mg L 1 ). (A colour version of this figure can be viewed online.)

Fig. 7 .

 7 Fig. 7. Biplot of the first two components of principal component analysis (PCA) including biological responses measured in different conditions of GO exposure (0, 0.05 or 0.1 mg L 1 ). Axis 1: 38.6%, Axis 2: 24,3%. Pl pleurodeles, C Chironomus, dev development stage, terato teratogenesis m microorganisms, Obs Observed diversity, Shan Shannon index, Beta beta diversity, Leaves leaf litter degradation. (A colour version of this figure can be viewed online.)

Figure S2 :

 S2 Figure S2: Multi-dimensional scaling analysis of microbiome data from the water (blue) and sediment (red) compartments. Ellipses were drawn with 95% confidence.

Figure S3 :

 S3 FigureS3: Effects of exposure to GO on bacterial communities from the water column of the microcosms as revealed by 16S amplicon sequencing. Observed richness, Chao1 and Shannon evenness indexes are compared between the three exposure conditions (0 mg.L -1 , 0.05 mg.L -1 and 0.1 mg.L -1 ). Observed richness, predicted richness and Shannon evenness index are compared between the three exposure conditions (0 mg/L, 0.05 mg/L and 0.1 mg/L) (A). MDS plot of bacterial communities based on unweighted Unifrac distances (B). Samples from different GO concentrations are represented by different colours with 95% confidence ellipses, while different sampling times are represented by forms surrounded by 95% confidence ellipses in grey.

Figure S4 :

 S4 Figure S4: Relative abundance of bacterial families from phyla Bacteroidetes and Proteobacteria from the water column.

Figure S5 :

 S5 Figure S5: Bacterial genera differentially observed between the exposure conditions compared to the control during the whole experiment in the water compartment (A) and in the biofilm (B).The relative abundance of genera Gemmata, Azospirillum and Flavobacterium were shown to decrease in a similar manner in the two GO conditions compared to the control group while genera Planctopirus and Armatimonas increased in the same order of manitude.

Figure S6 :

 S6 Figure S6: Relative abundance of bacterial families from phyla Bacteroidetes, Proteobacteria and Planctomycetes in the biofilm.

  and toxicological endpoints in multiple interacting species or communities, from different trophic levels, after direct exposure of organisms as well as through trophic routes [21e23].

• Corresponding autbor. E-mail address: florence.moucbet@ensat.fr (F. Mouchet} bttps://doi.org/10.1016/j.carbon2019.09.051 ecosystemic

Table 1

 1

	Physico-chemical characteristics of GO used for the experiment. at. %: atomic
	%; Csp2: sp 2 carbon; Sat.: shake-up satellites (p to p* transitions); TEM:
	transmission electron microscope; HRTEM: high resolution TEM; BET: Bru-
	nauer-Emett-Teller.
		Graphene Oxide
	Carbon content	69.0 at. %
	Oxygen content	31.0 at. %
	Csp2 graphene	35.5 at. %
	C OH/C O C	24.7 at. %
	C O	2.5 at. %
	O C O	5.3 at. %
	Sat.	1.4 at. %
	Number of layers (HRTEM)	1 5
	Lateral size (TEM)	0.2 8 mm
	Specific surface area (BET)

Table S1 :

 S1 Dissolved Organic Carbon (DOC) measured in water of microcosms at three sampling times during the experiment.

			Condition
	Sampling	0 mgGO.L -1	0.05 mgGO.L -1	0.1 mgGO.L -1
	T0	6.2 ± 3.6 mgDOC.L -1	5.1 ± 4.4 mgDOC.L -1	6 ± 5.2 mgDOC.L -1
	T1	1.7 ± 0.2 mgDOC.L -1	1.9 ± 0.05 mgDOC.L -1 1.6 ± 0.6 mgDOC.L -1
	T2	2.2 ± 0.06 mgDOC.L -1 2.4 ± 0.3 mgDOC.L -1	1.9 ± 0.2 mgDOC.L -1

Table S2 :

 S2 Table of results obtained from 2-ways ANOVA performed on data from the water column.

			Contaminant		Time of exposure	Contaminant*Time of exposure
		Df	F-	p-value	Df	F-value	p-value	Df	F-value p-value
			value
	Phylum
	Acidobacteria	2	0.800	0.464	2	1.584	0.232	4	0.868	0.501
	Actinobacteria	2	0.357	0.704	2	14.356	<0.001***	4	0.647	0.635
	Armatimonadetes	2	2.423	0.117	2	0.120	0.887	4	0.704	0.599
	Bacteroidetes	2	0.365	0.699	2	2.988	0.075	4	1.745	0.184
	Chloroflexi	2	1.245	0.311	2	2.530	0.107	4	0.877	0.496
	Ochrophyta	2	0.694	0.512	2	1.730	0.205	4	0.549	0.701
	Cyanobacteria	2	0.859	0.440	2	0.544	0.589	4	1.040	0.413
	Dependentiae	2	1.368	0.279	2	0.898	0.424	4	0.905	0.481
	Firmicutes	2	1.731	0.205	2	4.539	0.025*	4	1.040	0.414
	Fusobacteria	2	0.383	0.687	2	12.645	<0.001***	4	1.061	0.404
	Planctomycetes	2	0.870	0.435	2	0.190	0.828	4	0.368	0.828
	Proteobacteria	2	0.193	0.826	2	2.819	0.086	4	2.337	0.094
	Spirochaetes	2	1.680	0.214	2	1.535	0.242	4	1.162	0.360
	Verrucomicrobia	2	1.712	0.208	2	2.105	0.150	4	0.499	0.737
	WPS-2	2	1.5	0.249	2	1.5	0.249	4	1.5	0.244
	Family within
	Bacteroidetes
	Bacteroidaceae	2	0.511	0.608	2	10.927	<0.001***	4	0.184	0.943
	Chitinophagaceae	2	0.583	0.568	2	3.897	0.039*	4	1.373	0.282
	Crocinitomicaceae	2	2.322	0.126	2	3.082	0.070	4	1.780	0.176
	env.OPS17	2	0.555	0.583	2	1.998	0.164	4	0.567	0.689
	Flavobacteriaceae	2	1.526	0.244	2	1.123	0.346	4	0.475	0.753
	NS11	2	0.786	0.470	2	1.473	0.255	4	0.796	0.543
	NS9	2	0.926	0.400	2	1.096	0.355	4	0.798	0.541
	Saprospiraceae	2	0.814	0.458	2	1.326	0.290	4	0.416	0.795
	Sphingobacteriaceae	2	0.097	0.908	2	2.959	0.077	4	0.585	0.677
	Spirosomaceae	2	0.089	0.915	2	2.467	0.113	4	0.132	0.968
	Family within
	Proteobacteria
	Azospirillaceae	2	1.677	0.215	2	0.667	0.525	4	0.630	0.647
	Burkholderiaceae	2	2.265	0.132	2	3.003	0.074	4	0.378	0.821
	Devosiaceae	2	2.308	0.128	2	2.258	0.133	4	1.141	0.368
	Legionellaceae	2	2.612	0.100	2	0.127	0.881	4	2.027	0.133
	Paracaedibacteraceae	2	0.127	0.881	2	2.449	0.114	4	1.425	0.266
	Reyranellaceae	2	3.956	0.037*	2	1.703	0.210	4	0.820	0.528
	Rhizobiaceae	2	1.517	0.246	2	3.577	0.049*	4	1.364	0.285
	Rhizobiales	2	1.214	0.320	2	3.266	0.061	4	0.610	0.660
	Rhodobacteraceae	2	0.252	0.779	2	0.589	0.565	4	1.073	0.398
	Rickettsiaceae	2	1.564	0.236	2	0.793	0.467	4	0.671	0.620
	Sphingomonadaceae	2	2.286	0.102	2	11.325	<0.001***	4	1.324	0.298
	Xanthobacteraceae	2	0.761	0.481	2	8.263	0.0028**	4	0.249	0.910

Table S3 :

 S3 Table of results obtained from 2-ways ANOVA performed on metagenomics data from the biofilm

	Pirellulaceae	2	3.134	0.067	2	4.172	0.032*	4	1.073	0.398
	Rubinisphaeraceae	2	4.871	0.020*	2	0.191	0.827	4	0.054	0.994
	Schlesneriaceae	2	8.785 0.0021**	2	0.323	0.728	4	0.255	0.902
	Tepidisphaeraceae	2	Contaminant 2.517 0.108	2	Time of exposure 0.517 0.605	exposure Contaminant*Time of 4 0.378 0.820
		Df	F-	p-value	Df	F-value	p-value	Df	F-value p-value
			value
	Phylum
	Acidobacteria	2	5.119	0.0173*	2	1.641	0.221	4	3.119	0.041*
	Actinobacteria	2	0.388	0.683	2	1.351	0.284	4	0.085	0.985
	Armatimonadetes	2	8.994	0.002**	2	0.392	0.681	4	0.159	0.955
	Bacteroidetes	2	0.101	0.904	2	0.893	0.426	4	0.258	0.900
	Chloroflexi	2	2.100	0.151	2	5.644	0.012*	4	0.544	0.705
	Ochrophyta	2	0.398	0.677	2	0.312	0.735	4	0.213	0.928
	Cyanobacteria	2	1.564	0.236	2	4.724	0.024*	4	0.584	0.678
	Dependentiae	2	2.469	0.112	2	0.990	0.390	4	0.733	0.581
	Firmicutes	2	0.663	0.527	2	0.399	0.676	4	0.368	0.828
	Fusobacteria	2	4.214	0.031*	2	24.164	<0.001***	4	3.709	0.022*
	Gemmatimonadetes	2	2.624	0.099	2	3.541	0.0505	4	0.832	0.521
	Planctomycetes	2	7.991	0.003**	2	2.682	0.095	4	0.743	0.574
	Proteobacteria	2	0.167	0.847	2	0.340	0.716	4	0.374	0.824
	Spirochaetes	2	0.043	0.957	2	0.565	0.578	4	0.630	0.647
	Verrucomicrobia	2	1.057	0.368	2	1.858	0.184	4	0.806	0.537
	WPS-2	2	4.048	0.035*	2	6.939	0.006**	4	2.741	0.060
	Family in
	Bacteroidetes
	AKYH767	2	6.231	0.008**	2	0.6251	0.546	4	0.6263	0.649
	Chitinophagaceae	2	0.633	0.542	2	1.3487	0.284	4	0.152	0.959
	Crocinitomicaceae	2	3.144	0.067	2	5.140	0.017*	4	0.830	0.523
	env.OPS17	2	2.073	0.154	2	0.139	0.871	4	0.0511	0.994
	Flavobacteriaceae	2	2.003	0.163	2	0.310	0.736	4	1.0736	0.3985
	Microscillaceae	2	1.522	0.244	2	3.381	0.056	4	0.599	0.667
	NS11	2	2.992	0.075	2	0.409	0.670	4	0.950	0.457
	NS9	2	3.556	0.049*	2	0.528	0.598	4	0.509	0.729
	Saprospiraceae	2	0.626	0.545	2	2.392	0.119	4	1.846	0.163
	Spirosomaceae	2	0.299	0.744	2	0.777	0.474	4	1.379	0.280
	Family within
	Proteobacteria
	Burkholderiaceae	2	1.105	0.352	2	0.094	0.909	4	0.545	0.704
	Caulobacteraceae	2	2.296	0.129	2	3.481	0.052	4	0.438	0.779
	Chromobacteriaceae	2	1.403	0.271	2	2.626	0.099	4	0.218	0.924
	Pseudomonadaceae	2	4.549	0.025*	2	2.115	0.149	4	0.647	0.635
	Reyranellaceae	2	1.051	0.370	2	4.368	0.028*	4	0.502	0.734
	Rhizobiaceae	2	0.564	0.578	2	0.122	0.885	4	0.235	0.914
	Rhodobacteraceae	2	2.269	0.132	2	0.399	0.676	4	1.239	0.329
	Rhodocyclaceae	2	0.220	0.804	2	11.724	<0.001***	4	0.233	0.915
	Rickettsiaceae	2	0.494	0.618	2	5.042	0.018*	4	0.093	0.983
	Sphingomonadaceae	2	0.900	0.424	2	0.159	0.853	4	0.236	0.914
	Xanthobacteraceae	2	0.279	0.759	2	4.775	0.021*	4	0.924	0.471
	Family within
	Planctomycetes
	Gemmataceae	2	1.033	0.376	2	2.012	0.162	4	0.177	0.947
	Phycisphaeraceae	2	5.675	0.012*	2	0.936	0.410	4	1.760	0.180

Acknowledgments

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 696656 and No 78519. This research was also sup ported by the French Ministry of National Education, Higher Edu cation and Research. We would like to thank Alpha Diallo for the support during the experiment.

Conflicts of interest

The authors declare no conflict of interest.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.carbon.2019.09.051.

Assessment of graphene oxide ecotoxicity at several trophic levels using aquatic microcosms

Supplementary data

Lauris Evariste 1 , Antoine Mottier 1 , Laura Lagier