Hugo De Oliveira

Vincent Augusto

Baptiste Jouaneton

Ludovic Lamarsalle

Martin Prodel

Xiaolan Xie

Automatic and Explainable Labeling of Medical Event Logs with Autoencoding

Keywords: process mining, event log, healthcare data, patient pathways, autoencoding

Process mining is a suitable method for knowledge extraction from patient pathways. Structured in event logs, medical events are complex, often described using various medical codes. An efficient labeling of these events before applying process mining analysis is challenging. This paper presents an innovative methodology to handle the complexity of events in medical event logs. Based on autoencoding, accurate labels are created by clustering similar events in latent space. Moreover, the explanation of created labels is provided by the decoding of its corresponding events. Tested on synthetic events, the method is able to find hidden clusters on sparse binary data, as well as accurately explain created labels. A case study on real healthcare data is performed. Results confirm the suitability of the method to extract knowledge from complex event logs representing patient pathways.

I. INTRODUCTION

Data analytic regroups an extensive number of methods to investigate data produced in various systems such as industry, software engineering and healthcare. Knowledge extraction from such data is a lever to improve performances, to predict or simply to describe the reality of facts. Among different types of data, event logs are challenging to analyze because of the presence of time, the high variability of events, and the complex relations between events. Thus, the use of widespread data mining algorithm may not be fully straightforward for some applications. A wise preprocessing step to capture meaningful information may be necessary. Describing processes, these data are present in the manufacture industry, in software engineering and in healthcare [START_REF] Maita | A systematic mapping study of process mining[END_REF]. To analyze event logs, a data-driven approach named process mining has been proposed [START_REF] Van Der Aalst | Workflow mining: discovering process models from event logs[END_REF]. Between data mining and process modeling, event logs are impartially used to extra [START_REF] Aalst | Process Mining: Data Science in Action[END_REF].

The French national health insurance database (SNIIRAM) is a non-clinical claim database. Containing healthcare reimbursements of almost all French citizens, the amount of data is colossal. 66 million inhabitants were part of this database in 2015 [START_REF] Tuppin | Value of a national administrative database to guide public decisions: From the système national d'information interrégimes de l'Assurance Maladie (SNIIRAM) to the système national des données de santé (SNDS) in France[END_REF]. Among all reimbursement information contained in the SNIIRAM, patients' hospitalizations are provided. However, no precise medical information such as test results, H. De Oliveira, V. Augusto and X. Xie are with Mines Saint-Étienne, Univ Clermont Auvergne, CNRS, UMR 6158 LIMOS, Centre CIS, F -42023 Saint-Étienne France (e-mails: hdeoliveira@hevaweb.com; augusto@emse.fr; xie@emse.fr).

H. De Oliveira, B. Jouaneton, L. Lamarsalle and M. Prodel are with HEVA, 186 avenue Thiers, F-69465, Lyon, France (e-mails: bjouaneton@hevaweb.com; llamarsalle@hevaweb.com; mprodel@hevaweb.com).

X. Xie is also with the Antai College of Economics and Management, Shanghai Jiao Tong University, China.

imaging reports, or vital signs are available. Nevertheless, such a database is useful to map patient pathways [START_REF] Prodel | Stochastic simulation of clinical pathways from raw health databases[END_REF]- [START_REF] Oliveira | Optimal process mining of timed event logs[END_REF], perform medical data clustering [START_REF] Vandromme | A Scalable Biclustering Method for Heterogeneous Medical Data[END_REF] and prediction tasks [START_REF] Vandromme | Extraction and optimization of classification rules for temporal sequences: Application to hospital data[END_REF], [START_REF] Oliveira | Binary Classification on French Hospital Data: Benchmark of 7 Machine Learning Algorithms[END_REF]. Regarding healthcare processes, the complexity is multifold. Illustrations are, but not limited to, the presence of free text, the granularity of events analyzed, and the occurrence of multiple event simultaneously, leading to multiple codes describing a given event. These codes, representing medical activities of different types, could be numerous and often inherit of hierarchical structures [START_REF] Prodel | Optimal Process Mining for Large and Complex Event Logs[END_REF]. Even if the hierarchy could be useful to simplify codes and reduce the overall cardinality, the choice of the accurate level in order to produce meaningful events is not obvious, depends on the pathology or the health process studied and often requires a clinical expertise. This aspect of complexity is one of the main challenge regarding non-clinical claim database, such as the SNIIRAM.

Therefore, the main contribution of the present paper is a new methodology to analyze the complexity of events and produce meaningful labels. Using autoencoding and clustering, the proposed method creates artificial labels from initial data. These labels are assigned to events, transforming the raw event log by reducing the overall variability of events. The method provides transparency for practitioners by giving an interpretation for each created artificial label. In practice, the contribution consists of a preprocessing methodology to treat this particular complexity of events. As a result, available process mining tools1 can be used starting from event logs obtained via the proposed methodology.

This paper is organized as follows. An overview of related works is given in Section II. Preliminary notations are presented in Section III. Section IV introduces the problem addressed in this paper. The proposed methodology is described in Section V. To validate the method, a design of experiments is presented in Section VI, followed by a case study based on real-life healthcare data in Section VII. Finally, conclusions and perspectives are given in Section VIII.

II. LITERATURE REVIEW

Healthcare data analysis constitutes a large field to test and apply a wide spectrum of analytic methods. Among them, machine learning and more recently deep learning methods have been largely deployed. Electronic Health Records (EHR) have permitted the development of new models and methods, boosting the field publication activity [START_REF] Ravì | Deep Learning for Health Informatics[END_REF], [START_REF] Shickel | Deep EHR: A Survey of Recent Advances in Deep Learning Techniques for Electronic Health Record (EHR) Analysis[END_REF]. Among the tasks addressed by deep learning, supervised learning and concept embedding emerge for a majority of studies in healthcare [START_REF] Xiao | Opportunities and challenges in developing deep learning models using electronic health records data: a systematic review[END_REF]. In spite of the need to develop high performing predictive models, explainability has been highlighted as a key issue for future model developments [START_REF] Xiao | Opportunities and challenges in developing deep learning models using electronic health records data: a systematic review[END_REF], [START_REF] Miotto | Deep learning for healthcare: review, opportunities and challenges[END_REF]. Valuable for practitioners and experts of the medical field, the explanation of predictive results has already been addressed by deep learning studies [START_REF] Choi | Doctor AI: Predicting Clinical Events via Recurrent Neural Networks[END_REF], [START_REF] Choi | RETAIN: An Interpretable Predictive Model for Healthcare using Reverse Time Attention Mechanism[END_REF].

In addition to supervised learning and concept embedding, process mining is also a promising analytic method. By performing process discovery on event log data, process mining algorithms produce graphical and interpretative representations of occurring processes. The field is active and recent publications are numerous [START_REF] Maita | A systematic mapping study of process mining[END_REF], particularly in healthcare [START_REF] Erdogan | Systematic Mapping of Process Mining Studies in Healthcare[END_REF]. Initiatives of the research community to improve practices and knowledge sharing illustrate the field activity. 2 Recent topics addressed are, but not limited to, privacy [START_REF] Pika | Towards Privacy-Preserving Process Mining in Healthcare[END_REF], clinical guideline [START_REF] Gatta | Clinical Guidelines: A Crossroad of Many Research Areas. Challenges and Opportunities in Process Mining for Healthcare[END_REF], and data cleaning [START_REF] Martin | Interactive Data Cleaning for Process Mining: A Case Study of an Outpatient Clinic's Appointment System[END_REF]. Regarding process discovery, an optimal procedure to construct process models from raw data bases has been proposed by Prodel et. al. in 2018 [6]. Applied to healthcare data event logs, the proposed method has been adapted to take into account temporal information during the optimization process [START_REF] Oliveira | Optimal process mining of timed event logs[END_REF].

Patients' pathways analysis based on real-life healthcare data is valuable to represent and understand how patients' care occurs in real-life. Even if deep representation has been applied on clinical pathways [START_REF] Xu | Learning the Representation of Medical Features for Clinical Pathway Analysis[END_REF], process mining use graphical process models as a support for representation. Thus, it makes the method suitable to discover patients' pathways from raw data when the focus lies on interpretation. However, medical data can be complex to analyze, due to the variety of different medical codes used in claims databases (e.g. diagnoses, procedures and drugs). As a result, the labeling of events is a challenging step in data processing. A commune practice is the definition of labels by hand, based on expert knowledge [START_REF] Oliveira | Optimal process mining of timed event logs[END_REF]. The detection of hidden healthcare sub-processes has been proposed, using Hidden Markov Models (HMM) [START_REF] Alharbi | Towards Unsupervised Detection of Process Models in Healthcare[END_REF]. The presented method allows a reduction of complexity by the enrichment of the log with HMM-derived states, reducing complexity and saving experts time. A practical solution proposed by Prodel et al. in 2018 [6] is the creation of labels during the optimization, using the hierarchy of events from one type of codes, such as main diagnosis. However, when multiple codes from multiple coding systems characterize the events of a process, the selection of the right aggregation level and the combination of codes is not, to the best of our knowledge, a treated problem.

Therefore, the main contribution of this paper is the proposal of a general methodology to treat complex events such as multiple medical activities, in order to apply process mining. Following the definition of Lenz and Reichert [START_REF] Lenz | IT support for healthcare processespremises, challenges, perspectives[END_REF], the proposed method identifies the activities of the medical treatment process by analyzing the coded events of the organizational processes. The core of this methodology is based on recent work in representation learning. A widely used method in representation learning is autoencoding [START_REF] Goodfellow | Deep Learning[END_REF]. The general idea behind autoencoding is the learning of a structure which can encode and decode information while minimizing the loss of 2 http://pods4h.com/ information. By compressing the data, a transformation of the input representation is performed.

Thus, before formally introducing the problem and the proposed methodology, preliminary notations are presented in the following.

III. PRELIMINARIES

A formal description of the data involved is provided in the following, including events, traces and event logs.

Definition 1: Event. Each event denoted e is a couple (c, t, a) where:

• c is the related case ID of the event, with the id function id(e) = c returning the case ID of event e; • t ∈ T with T = N or R corresponds to the event time also called time-stamp, with the event time function time(e) = t returning the time-stamp of event e; • a is a nonempty set called activity set, each element a i ∈ a being an activity. Definition 2: Trace. A trace is a sequence of events denoted as σ = e 1 , e 2 , . . . , e m with m ∈ N * such that time(e k ) < time(e k+1 ) and ∀e, e ∈ σ, id(e) = id(e ). The size of the trace |σ| is defined as the number of events in σ.

Definition 3: Event log. An event log is a set of traces denoted as L = {σ 1 , σ 2 , . . . , σ n } with n ∈ N * . The size of the event log |L| is defined as the number of traces in L. Its length len(L), defined as len(L) = σ∈L |σ| gives the total number of events composing L.

According to the previous definitions, an event log is a group of traces, each trace being a succession of events characterized by activities occurring at a given time-stamp, composing an activity set.

Definition 4: Log activity set. Let L be an event log. The log activity set of L is the nonempty set A defined as:

• A = {a i |∃σ ∈ L, ∃e = (c, t, a) ∈ σ, a i ∈ a} • ∀σ ∈ L, ∀e = (c, t, a) ∈ σ, ∀a i ∈ a, ∃!A j ⊂ A|a i ∈ A j .
Thus, the log activity set is composed of activities encountered in L, with every single activity of event log L belonging to a unique subset A j of A.

Proposition 1:

A = j A j and j A j = ∅ (1)
Proposition 2:

card(A) = j card(A j ) (2)
Definition 5: Activity vector. Let e = (c, t, a) be an event, the activity vector x of e is defined such that x ∈ E X = {0, 1} d , E X being the activity vector space, with d = card(A). A mapping function vect() is also introduced, defined as:

vect : A → E X a → vect(a) = x with its inverse vect -1 : E X → A.
The previously defined functions allow a mapping between an activity set and its corresponding activity vector. The activity vector x of an event is the equivalent representation of an activity set using 1-of-k coding from all possible activities from A to the activity vector space E X . Therefore, a given event e could be defined as e = (c, t, a) or e = (c, t, x) knowing vect and vect -1 without any meaning loss regarding the log activity set. Definition 6: Activity matrix. The activity matrix of an event log L is defined as

M X = (x i ) i∈ 1,len(L) with dim(M X ) = len(L) × d.
The activity matrix M X of an event log L gives a binary representation of activity sets, which is a common representation in machine learning.

Example 1: Let us define an event log L = {σ 1 } with one trace σ 1 = e 1 , e 2 having two events such that e 1 = (c 1 , t 1 , a 1 ) and e 2 = (c 2 , t 2 , a 2 ) with:

• A = A 1 ∪ A 2 with A 1 = {A 1 , B 1 } and A 2 = {A 2 , B 2 , C 2 }; • case ids c 1 = c 2 = c; • time-stamps t 1 = 0, t 2 = 10; • activities a 1 = {A 1 , B 2 }, a 2 = {A 1 , B 1 , A 2 , C 2 }.
According to the previous definitions:

• d = card(A 1 ) + card(A 2 ) = 5; • x 1 = (1, 0, 0, 1, 0) and x 2 = (1, 1, 1, 0, 1). Example 2:
Table I presents a short event log related to patient pathways analysis using data as found in claims database. Each row is a hospitalization event. Events with the same case ID are ordered by increasing time stamp and represent a trace, which is a patient's hospitalization history. Medical activities are of three different categories: diagnosis, medical procedures and drugs, coded using standard notations ICD-10 (International Classification of Diseases 10 th revision), CCAM (Classification Commune des Actes Médicaux) and UCD (Unité Commune de Dispensation), respectively. According to previous definitions, A = A diag ∪ A med.proc. ∪ A drugs , with A diag = {Z511; Z5101}, A med.proc. = {ZZNL053} and A drugs = {261771; 9261110}. Moreover, d = card(A diag ) + card(A med.proc. ) + card(A drugs ) = 5. Representation using activity set and its equivalent using activity vector are presented in Table Ia andTable Ib, respectively. Thus, event log notations as well as activity vector and matrix formalization have been introduced. The following definitions introduce the problem setting.

IV. PROBLEM DEFINITION Definition 7: Label function and set. Given an event log L, a label function λ is a function such that:

λ : A → L (3)
with L being the label set.

The label function maps activities of an event log L to L, a set of possible labels for each activity vector.

Definition 8: Explaining function. Given an event log L and its label set L, an explaining function η on L is a function:

η : L → A (4)
The explaining function allows a mapping of a label l to interpretative elements a from the log activity set : a ⊂ A.

Definition 9: Activity clustering problem. Lets L be an event log with its log activity set A and its activity matrix M X , with dim(M X ) = len(L) × d. The activity clustering problem on L is defined as the search of the triple (L, λ, η) such that:

• L is a label set, L = ∅; • K min ≤ card(L) ≤ K max , with K min , K max ∈ N * ; • K max d; • λ is a label function; • η is an explaining function.
The main objective here is to find an accurate triple (L, λ, η) for the considered event log L. This problem can be seen as a clustering problem, with L being the set of cluster labels and λ the clustering function. The particularity here is for the input data M X to be sparse and of high dimensionality d in terms of features. The label set L should be finite, its cardinality (number of elements) being reasonable. A high cardinality of L will induce difficulties in process mining readability, also encouraging overfitting regarding A. This consideration leads to the proposed upper bound condition K max d. But a small cardinality for the label set will lead to a lack of information and medical meaning, motivating the lower bound K min . Moreover, as explainability is an essential constraint when dealing with medical pathways analysis, elements of L should keep a medical meaning, justifying the search of the explain function η. Given these considerations, it is assumed in this paper that there exist some relevant clusters hidden in medical activities of event log L. This assumption suggests that some combinations of elements of A which could characterize well a sufficient number of events e could be found.

In the following, a methodology based on autoencoding is proposed to find a relevant label set and a label function from a given event log.

V. PROPOSED METHODOLOGY

A. Overview

To solve the activity clustering problem for an event log L and find an accurate triple (L, λ, η), an autoencoding method is proposed. The idea is to transform space data from E X to a latent space E Z of reduced dimensionality, where similar elements are close to each other. This transformation is done using an autoencoding architecture, with an encoder f : E X → E Z and a decoder g : E Z → E X . In the latent space, because the dimensionality is reduced, applying clustering methods is meaningful. Latent clustering is proposed to assign for each element z a label through a function h : E Z → L. Thus, cluster labels defined through clustering in latent space will constitute the label set L defined in Definition 7:

λ : A vect(a) ----→ E X f (x) ---→ E Z h(z) ---→ L (5)
To construct the explaining function η, the function h -1 : L → E Z first returns for each cluster label l, the set of corresponding vectors in latent space Z l . Then, each activity z ∈ Z l can be decoded from E Z to E X using the decoding function g learned during the autoencoder training phase. This results in a set of activity vectors X l , from which an average vector Xl is computed. Finally, Xl is interpreted as a set of activities using vect -1 function:

η : L h -1 (l) ----→ E Z g(Z l ) ---→ E X vect -1 ( Xl ) -------→ A (6)
Consequently, the proposed method is composed of three steps: (1) autoencoder training, (2) latent space clustering, and (3) clusters' related activities decoding. These steps are presented in more detail in the following sections.

B. Autoencoder training

To perform clustering on sparse, binary, high-dimensional activity vectors, the data is transformed into a new space where variables are continuous and the dimensionality is lower. This transformation is performed using an autoencoder, trained on an activity matrix. In this paper, three methods for autoencoder training are investigated.

1) Autoencoder: An autoencoder (AE) is composed of two functions f and g, named the encoding and decoding functions, respectively. The encoder transforms a vector x from the input space into a new vector z from the latent space, its dimensionality being drastically reduced. The decoder takes the vector z from the latent space and decode it back to the input space, resulting in a new vector x . The training of a classic autoencoder is done by minimizing the reconstruction error, usually the binary cross entropy loss function. The dimensionality reduction allows a concentration of information while keeping the useful information in latent space for the reconstruction of input data. In this paper, encoder and decoder are constructed using symmetric feed-forward, fully connected neural networks.

2) Denoising autoencoder: A denoising autoencoder (DAE) is constructed using the same architecture as for an AE. A noisy vector x is created from x, which is encoded and then decoded. The loss function remains unchanged, while the goal of training is to be robust against artificially added noise and to keep useful information to decode data without noise.

3) Variational autoencoder: A variational autoencoder (VAE) is a particular autoencoder where the learned variables are parameters of a distribution. The encoder f is an inference network q(x|z) and the decoder g a generative network p(z|x). The reparameterization trick makes the training of the network possible using gradient descent optimization. The loss function here is the inverse of the expected lower bound ELBO defined as ELBO = E q(z|x) log p(x,z) q(z|x) . In practice, the single sample estimate log p(x|z) + log p(z) -log q(z|x) with z sampled from the inference network is optimized [START_REF] Kingma | Auto-Encoding Variational Bayes[END_REF].

C. Latent space clustering

Once an autoencoder is trained using activity matrix, a representation of every activity vector in latent space E Z can be obtained. In E Z , observations are characterized by a reasonable number of continuous features. As a result, applying clustering in this space is meaningful. K-means algorithm is used to learn the function h : E Z → L and create K clusters of similar observations in latent space. The parameter K corresponds to the final number of labels in L : K = card(L), respecting K min ≤ K ≤ K max (Definition 9). To find such an accurate value of K, one possible criterion could be to maximize the mean silhouette score [START_REF] Rousseeuw | Silhouettes: A graphical aid to the interpretation and validation of cluster analysis[END_REF], defined as:

S = 1 len(L) z∈M Z b z -a z max(a z , b z ) (7) with
• a z : the mean distance between z and all other points of the same cluster; • b z : the mean distance between z and all other points of the nearest next cluster.

D. Clusters' activities decoding

Once the latent space clustering has constructed a function h : E Z → L, the label function λ : A → L is fully defined. A label is assigned to each activity set of event log L, which can be used in the final process model as a node label. To construct the explaining function η : L → A, a methodology based on the decoding of each cluster's activities is proposed. The hypothesis formulated here is that for each cluster defined in latent space E Z , averaging the output of its activities decoding gives an overview of the cluster's meaning. In practice, the vect -1 function needs a threshold to convert an activity vector into decoded activities. The analysis of activity vectors can lead to a judicious choice for the threshold, as presented in Remark 1, Section VI.

To summarize this section, the presented methodology allows for the analysis of activity sets from event logs and defines a set of labels that can be used for process mining. Each event of the event log has a label defined through the label function, λ, and every label is interpreted by the explaining function η. A summary of preliminary notations is proposed in Table II.

VI. DESIGN OF EXPERIMENTS

In this section, a design of experiment on synthetic data is presented. The objectives of such an experiment are multiple:

(1) to verify the accuracy of the methodology in identifying hidden patterns (clusters) in event logs; [START_REF] Van Der Aalst | Workflow mining: discovering process models from event logs[END_REF] to demonstrate the veracity of the decoder in explaining clusters' labels; (3) to compare performances of autoencoding methods with direct clustering on sparse data; and (4) to benchmark autoencoding methods with one another.

A. Data description

The input data used in this paper are event logs, composed of events for which accurate labels are searched, according to the problem defined in Definition 9. Thus, labels hidden in the data were represented using groups of activities, among which vectors will be generated. Such a construction was motivated by data experts observations regarding non-clinical claims data: similar hospitalizations are often described by codes from a same "set", which can be approximated using clinicians knowledge. As an example, a stay for a given operation can be fill in the database using one code representing the operation, another for medical imaging procedures, some codes related to particular diagnosis related to complications of the operation or to the medical condition of the patient. Thus, there exist similarities between two stays for the same operation, because some codes are issued from a same "set", but they are rarely identical.

Considering these remarks, synthetic data were generated representing an activity matrix as defined in Definition 6, where each row represented an activity vector and each column a 1-of-k representation of activities. The number of different labels hidden in the data (i.e. the number of clusters to find) is referred to as κ. For each label k ∈ 1, κ , N k vectors were generated such that

N k = N (µ N , µ N 5 ), µ N = 2500
κ . The number of characteristic activities for each label k is referred to as

M C k ∈ N * , constructed such that M C k = N (µ C , µ C 5 ) with µ C = α × β × µ N .
The number of all different activities involved is referred to as M . For each activity vector of a given label k, a number M a = α × µ N ∈ N * of activities was randomly chosen to construct it. The number M a corresponds to the number of activities randomly chosen among characteristic activities. For these activities, the corresponding attribute value was set to 1, keeping 0 otherwise. An overlapping ratio γ is also introduced, representing the quantity of activities of a label shared with the nearest one. Moreover, a number N noisy = 250 of noisy events was generated. These noisy events were composed of activities among all possible activities in the constructed data set, with no particular pattern related to a given hidden label k. As a result, the total number of traces (i.e. the number of rows) is N = N noisy + κ k=1 N k . The parameters that have been chosen in the design of experiments lead to an approximated sparsity S between 0.9 and 0.99. An example of synthetic data generated is presented in Figure 1.

B. Experiment description

The autoencoder methods AE, DAE, and VAE, presented in Section V-B, were implemented and compared in terms of performances. Neural networks used are feed-forward networks, composed of four fully connected layers of size 10 × d latent , 5 × d latent and d latent , the latter being the dimensionality of the latent space defined as d latent = 8. For DAE, noise was defined as randomly selecting elements in vectors and changing their values (from 0 to 1 or 0 to 1). 1% of noise is added in every vector. For VAE, d latent couple of parameters of Gaussian distributions were learned as latent variables. The inverse single sample Monte Carlo estimate of the ELBO was used as the loss function to minimize during training. For each parameter combination of the data, autoencoder training was done using a symmetric architecture between the encoder and the decoder. Dropout and L-2 regularization were used for each layer to prevent the training from overfitting. The chosen optimizer for training was Adam, with mini-batch of size 32. The total number of epochs was set to 1000. Of the overall data, 80% was used for training, while 20% was kept to evaluate validation error for early stopping (with a patience of 25 iterations). After autoencoder training, K-mean clustering in latent space was applied using all training and validation data, K being fixed by maximizing mean silhouette score for K ∈ K min , K max with K min = 2 and K max = 15.

Performances were evaluated regarding clustering on the one hand, and explainability on the second hand. To evaluate clustering performances, an automatic procedure constructed a confusion matrix between hidden labels and found clusters, maximizing the accuracy (the sum of the diagonal) by permuting columns only (found clusters) to align proposed clusters with potential corresponding ones in hidden ones. The accuracy of resulting confusion matrix was then computed, defined as the ability of a method to accurately assign the right label to each event. To evaluate the ability of the method to explain found clusters, the explaining F -score F η is introduced. Let k ∈ 1, κ = K true be the label of the cluster hidden in data, and c : k → c(k) ∈ K pred being a function returning the corresponding cluster label according to previously described confusion matrix optimization. The average of decoded elements from the cluster c(k) is computed, and its activity set k∈K pred card(a pred c(k) ) the explaining precision. These expressions are analogous to classical binary classification metrics. A high explaining recall means a high ability of the explaining function to find corresponding activities with hidden activities of the identified labels. Furthermore, a high explaining precision corresponds to a decoding that keeps interesting activities without being too general. Ideally, each discovered label corresponds to a hidden one. This is not necessarily the case, as the number of discovered and hidden clusters could be different. Fewer predicted than hidden clusters will impact the explaining recall, and more predicted than true clusters will impact explaining precision.

Remark 1: As mentioned in Section V-D, the threshold will have significant impact in the previously defined metrics. In the design of experiments conducted here, an automatic approach was used. For a cluster l, a list of all decoded values from the average decoding vector Xl was constructed. All elements of this list were ordered in descending order. By differentiation of this curve another list of values was obtained. The minimum value of the resulting curve was used to automatically define a judicious threshold for keeping activities in the explaining set of the related cluster.

For every combination of parameters, 10 data sets were constructed. Columns (activities) were shuffled, right before the shuffling of rows (events). The proposed method was applied using the previously described autoencoders (AE, DAE and VAE) as well as a direct K-mean clustering without autoencoding step, tested as a baseline (referred to as BASIC). Performances were analyzed through mean and standard deviation of clustering and explaining metrics. All algorithms and experiments were conducted using Python 3.7 and Tensorflow 1.14. A schematic description of experiments is presented in Figure 2.

C. Results

Results 3 are summarized in Table III. A total of 24 experiments of increasing difficulty were conducted. For each combination of parameters, the accuracy and the F η score of the tested methods are presented. Results show that the autoencoding methods outperform the direct clustering in 3 Detailed results regarding the design of experiments of Section VI, as well as the case study described in Section VII can be found as supplementary materials on the following website: https://artemis-emse-laparo.hevaweb.com/ sparse high-dimensional space (except regarding experiments 15 and 21 where DAE performances are inferior). Therefore, autoencoding plays an important role in data transformation for the proposed methodology. Furthermore, the results highlight that VAE always outperforms the other methods regarding both accuracy and F η . Event though the standard deviation increases for difficult experiments, VAE shows a lower variation compared to the other methods. As a conclusion, results motivate the choice of VAE as part of the proposed method to obtain accurate clusters and explain them.

VII. CASE STUDY

After proving the accuracy of the method on synthetic data, this section is dedicated to demonstrate the relevance of the method on real healthcare data. In that purpose, a medical case study is presented, were process mining was deployed to extract knowledge about patients' hospital pathways.

A. Overview

A laparotomy is an abdominal surgery consisting of a large incision of the abdomen, sometime necessary to investigate abdominal pain. Incisional hernia (IH) is one of the possible complications following laparotomy. These complication consists in a protrusion of the tissues of the abdomen through the abdominal muscle. The repair on an IH is a common surgery, which can lead to chronic pain and decreased quality of life. Colorectal surgeries, bariatric surgeries and abdominal aortic aneurysm are laparotomy surgeries that may lead to IH [START_REF] Phan | Clinical Pathway Analysis Using Process Mining and Discrete-Event Simulation: an Application to Incisional Hernia[END_REF]. In this case study, we focused on patients developing an IH after a laparotomy operation. The objective was to apply the previously defined methodology to label raw medical event logs before applying process mining. A process mining study using manual labeling was also performed, to illustrate the relevance of the methodology to automatically define interesting labels.

B. Methods

The data 4 5 0.10 2 0.25 0.54 0.09 0.68 0.16 0.70 0.12 0.88 0.01 0.69 0.13 0.84 0.25 0.90 0.20 0.94 0.12 10 5 0.10 5 0.00 0.39 0.07 0.57 0.16 0.64 0.15 0.82 0.02 0.40 0.08 0.68 0.26 0.77 0.25 1.00 0.00 11 5 0.10 5 0.10 0.36 0.05 0.60 0.18 0.69 0.13 0.83 0.01 0.36 0.07 0.72 0.25 0.90 0.25 0.97 0.10 12 5 0.10 5 0.25 0.45 0.09 0.64 0.17 0.63 0.17 0.84 0.02 0.41 0.16 0.79 0.27 0.80 0.28 1.00 0.00 13 10 0.05 2 0.00 0.50 0.06 0.61 0.18 0.62 0.15 0.87 0.04 0.62 0.10 0.72 0.29 0.84 0.24 0.98 0.06 14 10 0.05 2 0.10 0.54 0.06 0.62 0.22 0.66 0.14 0.88 0.02 0.71 0.09 0.73 0.37 0.88 0.25 0.97 0.10 followed by an IH within 5 years after the operation were selected. This resulted in a total number of 7, 906 patients included in the study, for which all hospitalization information was extracted. Each patient's hospitalization was transformed into a trace of his ordered medical activities. Thus, the activity set was structured as follows:

A = A MD A AD A MP A D A TAD (8)
where:

• A MD is the set of main diagnoses, the reasons of the hospitalization (MD, using ICD-10 coding system); • A AD is the set of additional diagnoses (AD, using ICD-10 coding system); • A MP is the set of medical procedures (MP, using French CCAM coding system); • A D is the set of delivered drugs (D, using French UCD coding system linked to ATC -Anatomical Therapeutic Chemical -classes); • A TAD is the set of drugs under temporary authorization for delivery in French hospitals (using French LPP -Liste des Produits et Préstations -coding system). Moreover, for each activity code, hierarchical knowledge (codes of upper levels in the hierarchy) was added as part of the corresponding activity set. This procedure enables relations between activity codes of a same group during autoencoding. It also enriches the explainability of clusters, by providing hierarchical knowledge and setting the level of precision in coding depending on clusters, as shown in the following results.

Stays related to dialyses or chemotherapy, which are known to appears very frequently, have been filtered. Codes appearing less than 50 times were also filtered, resulting in keeping 95.0% of codes while decreasing the size of the log activity

C. Results

Results obtained by automatic labeling (Figure 3a) were compared with a process model of the same dimensionality, constructed with the same process mining procedure but start- ing from an event log with manually defined labels according to authors' prior knowledge about the pathology (Figure 3b). The process models read from left to right, the size of nodes and edges being, for each process model, proportional to the number of patients represented. Explanation results regarding clusters obtained and visualized in Figure 3a are detailed in Table IV. According to Figure 3, similarities are observed between the process models. Most frequent medical procedure codes related to laparotomy are subcategories of the "Therapeutic acts on the digestive system" code in the hierarchy. Thus, the related cluster (label 7) appears at the beginning of the pathways. Label 2, which contains codes related to IH, appears in the following of the pathway, as well as label 14 (combining codes related to laparotomy and IH). Also, stays related to diagnostic procedures (labels 6 and 12) and more precisely to an endoscopy of the alimentary canal (label 10) occupy a significant place in the process analyzed. These stays were not considered during manual labeling but were pointed out by the automatic labeling procedure based on raw data. They may be related to patient's medical control or investigation regarding suspicion of complication after operation. This example illustrates that other interesting information can be extracted from raw data with minimal initial input from the user. However, by comparing the replayability score (which gives a quantitative fitness measure of the resulting process models), a gap was observed between automatic (42, 7%) and manual (77.4%) labeling. The main explanation may be provided by the first laparotomy node (representing 3, 802 vs 7, 906 patients, respectively). In practice, 549 different medical procedure codes, from different chapters of the hierarchy, were selected by medical experts to identify laparotomy in the database. Even if most of the codes are gathered in label 7, remaining laparotomy stays were grouped in other less frequent clusters, which do not appear in the final process model because of the size constraint in optimization. Thus, even if a quantitative replayability gap remains between the two presented methods, the qualitative interpretation resulting from the pathway analysis remains similar, as most of the interesting events were pointed out. Moreover, the final process model and the explanation of clusters furnish an interesting base for discussion with medical experts.

VIII. CONCLUSION In this paper, a methodology to handle the complexity of event logs regarding activities was presented. Based on autoencoding, artificial labels to characterize these events are created, which can be used to apply process mining. Explainability of each label is possible through decoding, which allows the practical application of this method in fields like healthcare where transparency is essential. A design of experiments was presented, designed to mimic non-clinical claims databases regarding authors' knowledge. The ability of the method to both create relevant clusters and explain them accurately was demonstrated. In particular, the Variational Auto-Encoder shows better performances than others tested autoencoders, motivating the use of such learning methods for further applications. Finally, a case study has been presented, illustrating the potential of the methodology when applied on real healthcare data. The presented method sounds promising as a preprocessing solution for process mining, to handle the complexity of medical activities in non-clinical claims databases and of other similar databases.

Further work will focus on the deployment of the method in new case study, to experience the method and generalize its utilization. Particularly, a focus on providing interactive tools to explore the results and facilitate discussions with clinical experts will be made. As the complexity of medical event logs depends on the database, other complexities could be considered in further studies. As an example, the integration of free text with structured medical information to create pertinent medical labels could be interesting to consider in future works. As noted by Helm et al [START_REF] Helm | Towards the Use of Standardized Terms in Clinical Case Studies for Process Mining in Healthcare[END_REF], a lack of sufficient coding in existing case studies remains. Therefore, addressing such complexities in the method will be interesting for practical uses on case studies with insufficiently coded data. Also, the present methodology is a preprocessing step, applied before deploying a process discovery algorithm (in our example, based on optimization). An interesting subject could be the fusion of both steps, by integrating the labeling step directly during the optimization procedure of creating the final process model. Also, as the K-mean algorithm was used to perform clustering in the latent space, future work will focus on testing other clustering algorithms. On a larger scope, future works will focus on the use of the proposed methodology to perform supervised learning on complex event logs. The proposal of a transparent classification algorithm is of interest, particularly for patients pathways data. Furthermore, the bridging of process mining and deep learning is an interesting research track. The use of recent advancements in deep learning for process analysis and prediction seems promising, in particular if process mining is used as an interface between model learning and human understanding.

Fig. 1 :

 1 Fig. 1: Example of generated data. Here, κ = 5, α = 0.05, β = 2 and γ = 0.1, leading to S = 0.9 with each white pixel being a value of 1 and black pixels being a value of 0.

  a pred c(k) is compared to characteristic activities of the corresponding label a true k . To quantitatively analyze the decoding performances, the explaining F -score is defined such as: F η = 2×Rη×Pη Rη+Pη with R η = k∈K true a true k a pred c(k) k∈K true card(a true k ) the explaining recall and P η = k∈K pred a true k a pred c(k)

Fig. 2 :

 2 Fig. 2: Schematic representation of the design of experiments.

  Automatic labeling. (b) Manual labeling.

Fig. 3 :

 3 Fig. 3: Comparison of process models obtained using automatic labeling (left) vs manual labeling (right).

TABLE I :

 I Example of an event log of patient pathways.

TABLE II :

 II Notations summary.

	name	notation	name	notation
	activity	ai	log activity set	A
	activity set	a	activity vector	x
	time-stamp	t	activity matrix	M X
	event	e	activity vector space	E X
	trace	σ	encoder	f
	event log	L	decoder	g
	label set	L	latent vector	z
	label function	λ	latent matrix	M Z
	explaining function	η	latent space	E Z

  were extracted from the SNIIRAM database. All anonimized patients with a first laparotomy operation in 2010,

							ACCURACY		Fη
			Exp.			BASIC	AE	DAE	VAE	BASIC	AE	DAE	VAE
		κ	α	β	γ	AVG STD AVG STD AVG STD AVG STD AVG STD AVG STD AVG STD AVG STD
	1	5 0.05 2 0.00 0.55 0.06 0.66 0.16 0.71 0.05 0.90 0.01 0.71 0.13 0.80 0.27 0.95 0.11 1.00 0.00
	2	5 0.05 2 0.10 0.55 0.09 0.66 0.16 0.71 0.08 0.89 0.01 0.68 0.17 0.82 0.32 0.95 0.14 1.00 0.00
	3	5 0.05 2 0.25 0.56 0.07 0.68 0.16 0.72 0.04 0.89 0.03 0.71 0.17 0.82 0.25 0.92 0.12 0.97 0.09
	4	5 0.05 5 0.00 0.38 0.06 0.59 0.16 0.61 0.18 0.85 0.03 0.37 0.14 0.73 0.27 0.77 0.31 1.00 0.00
	5	5 0.05 5 0.10 0.37 0.04 0.68 0.14 0.57 0.18 0.86 0.03 0.37 0.11 0.83 0.26 0.67 0.37 1.00 0.00
	6	5 0.05 5 0.25 0.43 0.06 0.59 0.19 0.63 0.15 0.85 0.01 0.47 0.13 0.67 0.35 0.80 0.29 1.00 0.00
	7	5 0.10 2 0.00 0.51 0.07 0.69 0.13 0.70 0.07 0.88 0.01 0.65 0.13 0.86 0.23 0.91 0.15 1.00 0.00
	8	5 0.10 2 0.10 0.52 0.07 0.63 0.16 0.70 0.08 0.89 0.01 0.63 0.15 0.80 0.27 0.93 0.16 1.00 0.00
	9

TABLE III :

 III Clustering and explainability performance, measured by accuracy and F η , respectively. For all parameters combinations, average and standard deviation over 10 replications are presented. Best values are highlighted in bold.

TABLE IV

 IV

	: Explanation of clusters appearing on process
	model: relevant decoded activities with corresponding activity
	set and level in the hierarchy.
	set by 85.7%. The final event log constructed for the study
	contained, for 7906 traces (patients), 57533 events (stays) and
	2228 unique activity codes. The previously defined method-
	ology was conducted on the resulting activity matrix (of size
	57533 × 2228), using VAE as autoencoder. The number of
	clusters used was K = 15, which qualitatively appears as
	a wise trade-off between explainability of clusters and final
	process model readability. The process mining framework
	used was the one proposed in [8], designed for application
	to medical event logs. The maximum number of nodes, edges
	and positions for process model optimization was fixed to 15,
	25 and 5, respectively.

Such as ProM, Disco, PM4Py or bupaR.

Access to an extraction of the SNIIRAM database was provided by the French CNIL under the agreement number DR-2019-147.


Céline Garcia
email: celine.garcia@ifremer.fr

Aurélie Mesnil

Delphine Tourbiez

Mirna Moussa

Christine Dubreuil

Amélie Gonçalves De Sa

Bruno Chollet

Yoann Godfrin

Subsp Francensis

1 Nov

Lionel Degrémont

Delphine Serpin

Marie-Agnès Travers

Description of Vibrio aestuarianus subsp. cardii subsp. nov. pathogenic to the edible cockles, Cerastoderma edule in France and establishment of Vibrio aestuarianus subsp. aestuarianus subsp. nov. and Vibrio aestuarianus subsp

Keywords: nov., pathogenic to the edible cockles Cerastoderma edule nov. and Vibrio aestuarianus subsp. francensis subsp. nov Shellfish, Bivalvia, mortality, Vibrio aestuarianus, cockle, bacterial pathogen 1.5 Repositories MK307685-307690, MK307692, MK307693, MK307695, MK307696

HAL is

INTRODUCTION

The edible cockle Cerastoderma edule is one of the most common bivalve species widely distributed along the European coasts, and is commercially important in the UK, Netherlands, Spain and France. For many years, mortalities of cockles and dramatic declines in their populations have been reported [START_REF] Longshaw | A review of the infectious agents, parasites, pathogens and commensals of European cockles (Cerastoderma edule and C. glaucum)[END_REF]. Predation, diseases, climatic events, pollution, poor recruitment and over-fishing have all been suggested as driving factors [START_REF] Ducrotoy | A comparison of the population dynamics of the cockle (Cerastoderma edule, L.) in north-western Europe[END_REF]. Cockles are host to a wide variety of viruses, bacteria, fungi and protistans. Many of these microorganisms are normally innocuous but can become problematic if prevailing environmental conditions are poor. Some parasites have been isolated during mass mortality events [For review, 1]: Protistan Marteilia cochillia in Spain [START_REF] Carrasco | Marteilia cochillia sp nov., a new Marteilia species affecting the edible cockle Cerastoderma edule in European waters[END_REF], gregarines Nematopsis scheideri and Nematopsis incognito in Portugal [START_REF] Azevedo | Fine structure of the apicomplexa oocyst of Nematopsis sp. of two marine bivalve mollusks[END_REF][START_REF] Ordas | Histopathological study of abnormal mortalities of cockle (Cerastoderma edule) in Carril (Galicia, NW Spain)[END_REF], Amoeba [START_REF] Azevedo | Some observations on an amoeba, destructive agent of the gills of Cerastoderma edule (Mollusca, Bivalvia)[END_REF], Haplosporidia [START_REF] Azevedo | Description of Haplosporidium edule n. sp. (Phylum Haplosporidia), a parasite of Cerastoderma edule (Mollusca, Bivalvia) with complex spore ornamentation[END_REF], and a wide range of Digenea (trematodes) [START_REF] De Montaudouin | Digenean trematode species in the cockle Cerastoderma edule: identification key and distribution along the north-eastern Atlantic shoreline[END_REF][START_REF] Magalhães | Spatio-temporal variation of trematode parasites community in Cerastoderma edule cockles from Ria de Aveiro (Portugal)[END_REF]. Cockles also host emerging parasites, such as Perkinsus chesapeaki, newly described in Europe [START_REF] Carrasco | Perkinsus chesapeaki observed in a new host, the European common edible cockle Cerastoderma edule, in the Spanish Mediterranean coast[END_REF]. They can also suffer from a transmissible disease known as disseminated neoplasia [START_REF] Metzger | Widespread transmission of independent cancer lineages within multiple bivalve species[END_REF][START_REF] Díaz | Long-term epidemiological study of disseminated neoplasia of cockles in Galicia (NW Spain): temporal patterns at individual and population levels, influence of environmental and cockle-based factors and lethality[END_REF]. Some studies have suggested that different factors, such as genetic alterations, viruses, retrotransposons and contaminants [START_REF] Díaz | Long-term epidemiological study of disseminated neoplasia of cockles in Galicia (NW Spain): temporal patterns at individual and population levels, influence of environmental and cockle-based factors and lethality[END_REF], can act in concert [START_REF] Barber | Neoplastic diseases of commercially important marine bivalves[END_REF].

However, bacteria have rarely been isolated associated with mortalities, apart from mycoplasma-like organisms in Spain [START_REF] Azevedo | Occurence of an anusual branchial mycoplasma-like infection in cockle Cerastoderma edule (Mollusca, Bivalvia)[END_REF].

Vibrio are the most abundant cultivable bacteria found in the marine environment. Free, associated with plankton or with marine animals, bacteria belonging to this genus are classified in clades [START_REF] Sawabe | Inferring the evolutionary history of vibrios by means of multilocus sequence analysis[END_REF][START_REF] Sawabe | Updating the Vibrio clades defined by multilocus sequence phylogeny: proposal of eight new clades, and the description of Vibrio tritonius sp[END_REF][START_REF] Thompson | Phylogeny and Molecular Identification of Vibrios on the Basis of Multilocus Sequence Analysis[END_REF][START_REF] Thompson | Biodiversity of Vibrios[END_REF] clustering more than 100 different species [START_REF] Romalde | New Vibrio species associated to molluscan microbiota: a review[END_REF]. Some of these species have virulent strains implicated in fish and shellfish diseases either alone, in bacterial consortia [START_REF] Goudenege | A single regulatory gene is sufficient to alter Vibrio aestuarianus pathogenicity in oysters[END_REF][START_REF] Lemire | Populations, not clones, are the unit of vibrio pathogenesis in naturally infected oysters[END_REF][START_REF] Lasa | Dynamics of the Pacific oyster pathobiota during mortality episodes in Europe assessed by 16S rRNA gene profiling and a new target enrichment nextgeneration sequencing strategy[END_REF] or in concert with other pathogens [START_REF] De Lorgeril | Immune-suppression by OsHV-1 viral infection causes fatal bacteraemia in Pacific oysters[END_REF]. For bivalves, the main pathotypes belong to the Splendidus clade (V. tasmaniensis, V. crassostreae, V. splendidus) or the species V. tapetis, V. tubiashii, V. bivalvicida, V. europaeus or V. aestuarianus [START_REF] Romalde | New Vibrio species associated to molluscan microbiota: a review[END_REF][START_REF] Travers | Bacterial diseases in marine bivalves[END_REF].

Vibrio aestuarianus was first isolated in the 1980s in the USA, from seawater and several shellfish including oysters, clams and crabs presenting no signs of disease [START_REF] Tison | Vibrio-Aestuarianus -a New Species from Estuarine Waters and Shellfish[END_REF]. This species was isolated for the first time in France in 2001, during mortality events of Pacific oyster, Crassostrea gigas. Garnier et al. [START_REF] Garnier | Molecular and phenotypic characterization of Vibrio aestuarianus subsp francensis subsp nov., a pathogen of the oyster Crassostrea gigas[END_REF] compared the French strains with the American one and found that the French strains had some different biochemical and pathogenic characters. These authors concluded that the French strains could constitute a new subspecies within V. aestuarianus and proposed to name the American strain V. aestuarianus subsp. aestuarianus and the French strains V. aestuarianus subsp. francensis, although these two subspecies were never officially validated.

The present study describes the isolation and classification of bacteria from cockle mortality events that occurred in France in 2012 and 2015. Strains were described and sequenced that allowed us to reproduce and describe this new cockle disease. These strains were compared with the French 'oyster' strains [START_REF] Garnier | Molecular and phenotypic characterization of Vibrio aestuarianus subsp francensis subsp nov., a pathogen of the oyster Crassostrea gigas[END_REF] and with the type strain of Vibrio aestuarianus (LMG 7909T) [START_REF] Tison | Vibrio-Aestuarianus -a New Species from Estuarine Waters and Shellfish[END_REF], and led us to propose two novel subspecies of Vibrio aestuarianus.

Isolation and Ecology

In 2012, mass mortality events of Cerastoderma edule were reported in two wild beds on the English Channel Coast: Baie des Veys and Baie de Somme, the latter of which is the largest cockle bed in France. In 2015, another mortality event occurred in both these locations as well as in Binic, also on the English Channel coast (Figure S1). Fishers reported these mortality events to the competent authority and national network for surveillance and monitoring of mollusc health, REPAMO (REseau de PAthologie des MOllusques), in charge of recording mortality cases and analysing representative samples associated with them. In 2012 and 2015, 15 to 30 adult cockles were collected from the different wild beds. Each collected individual was examined for macroscopic characteristics and processed using either histological or bacteriological analysis in addition to molecular analyses (Table S1). The mortality rate in 2012 was estimated at 80-89% according to site. For each site, the mortality rate was lower in 2015 than in 2012 ranging from 33 to 49% (Table S1). In all cases, cockle mortality was uniform, massive and sudden, concerning mainly adults. No other species were affected by these mortality events except in Baie des Veys where a number of stocks of Crassostrea gigas spat and juveniles also suffered mortalities. The cockle mortality occurred during summer, and lasted from 1 to 2 months, according to fishermen. No specific environmental events were reported during the summers of 2012 or 2015, except a high level of macroalgae at Binic. The sampled individuals were mainly moribund, lying on the bottom and characterized by very slow valve closing valves and limited intrapalleal fluid. No specific macroscopic signs such as pustules, conchyolin deposits or necrosis and degradation of the hinge could be observed.

For the detection and quantification of bacterial pathogens, pieces of gill, muscle and mantle of five moribund cockles per sample were individually crushed in 150 μl sterile seawater, diluted 100fold and 10,000-fold and spread onto Marine agar (Conda, Spain) incubated for 48-72 h at 20°C.

Colonies were enumerated and those bacterial colonies found in abundance were isolated and tested by duplex real time PCR for the detection of V. aestuarianus and V. Splendidus-related. [START_REF] Saulnier | Development of a duplex Taqman realtime PCR assay for rapid identification of Vibrio splendidus-related and V. aestuarianus strains from bacterial cultures[END_REF].

Pure cultures of these bacterial strains were conserved frozen at -80°C in Zobell broth with 15% glycerol. These bacteriological analyses revealed that all analysed samples yielded numerous bacterial colonies. The species V. aestuarianus was detected in all samples after 48h at 20°C (Table S2). Bacteria of the Splendidus clade were also isolated in some samples but these bacteria were not the dominant strains in the Petri dishes. Moreover, V. aestuarianus DNA was detected by molecular analyses [START_REF] Saulnier | Real-time PCR assay for rapid detection and quantification of Vibrio aestuarianus in oyster and seawater: A useful tool for epidemiologic studies[END_REF] in all samples, at a generally high detection frequency. DNA of OsHV-1 (Ostreid Herpevirus type 1) known to contribute to mortality events in Crassostrea gigas [START_REF] De Lorgeril | Immune-suppression by OsHV-1 viral infection causes fatal bacteraemia in Pacific oysters[END_REF][START_REF] Garcia | Ostreid herpesvirus 1 detection and relationship with Crassostrea gigas spat mortality in France between 1998 and 2006[END_REF], was not detected in any of the samples (Table S2).

Fifteen to twenty individuals per sample were also processed for histological analyses. After 48 h in Davidson's fixative, tissues were maintained in 70% ethanol until histological preparation by dehydration and embedding in paraffin following standard procedures. Sections of 2-3 µm thickness were stained with haematoxylin and eosin. Histological analyses revealed the presence of other organisms at low infection intensity (a few parasites per individual), despite a high prevalence in some cases such as gregarine spores or rickettsia-like organisms (Table S2).

16S RNA phylogeny

A total of 59 bacterial strains were isolated during periods when the REPAMO surveillance network reported high levels of cockle mortalities. After isolation, the predominant bacteria appeared small (1-2 mm), translucid, regular and white to cream-colored on marine agar after 48 h at 20°C (Figure S2). For genotyping, DNA was extracted by boiling in distilled water [START_REF] Saulnier | A large-scale epidemiological study to identify bacteria pathogenic to Pacific oyster Crassostrea gigas and correlation between virulence and metalloprotease-like activity[END_REF] and the ldh (lactate dehydrogenase) gene was amplified and sequenced according to methods previously described [START_REF] Goudenege | A single regulatory gene is sufficient to alter Vibrio aestuarianus pathogenicity in oysters[END_REF], as were the 16S rRNA and gyrB genes [START_REF] Thompson | Phylogeny and Molecular Identification of Vibrios on the Basis of Multilocus Sequence Analysis[END_REF] (Table S4). Sequences were aligned and phylogenetic trees built using Mega X [START_REF] Kumar | Molecular Evolutionary Genetics Analysis across Computing Platforms[END_REF]. The trees were drawn using the Neighbor-Joining method with the Kimura two-parameter model [START_REF] Saitou | The Neighbor-Joining Method: a new method for reconstructing phylogenetic trees[END_REF] and maximum parsimony method [START_REF] Eck | Atlas of protein sequence and structure[END_REF]. A total of 17 isolates (isolated in 2012 and 2015 from cockle mortality events) were identified as belonging to Vibrio aestuarianus by qPCR (Table S3) and gene sequencing (16S rDNA and gyrB, Table S4). They showed a high degree of similarity with 16S rDNA reference sequences of Vibrio aestuarianus (99% with the type strain LMG 7909T). Phylogenetic analyses based on the ldh gene sequence of these strains led to their affiliation to the species V. aestuarianus, in a group close but apart from the oyster strains and type strain. Moreover, oyster strains formed an individual branch clearly separated from the type strain of V. aestuarianus (Figure 1) [START_REF] Garnier | Molecular and phenotypic characterization of Vibrio aestuarianus subsp francensis subsp nov., a pathogen of the oyster Crassostrea gigas[END_REF].

This classification was confirmed by DNA-DNA hybridization (DDH) performed by the BCCM TM /LMG laboratory at Gent University. After genomic DNA extraction according to a modified Wilson protocol [START_REF] Wilson | Preparation of genomic DNA from bacteria[END_REF], hybridizations were performed in microplates at 40°C [START_REF] Ezaki | Fluorometric Deoxyribonucleic Acid-Deoxyribonucleic Acid Hybridization in Microdilution Wells as an Alternative to Membrane Filter Hybridization in which Radioisotopes Are Used To Determine Genetic Relatedness among Bacterial Strains[END_REF] and DNA relatedness values (%) based on means of a minimum of six hybridizations. DDH of two strains isolated from diseased cockles in 2012 and 2015 and one strain isolated from diseased oysters were performed with Vibrio aestuarianus, V. ordalii and V. salilacus reference strains. With more than 70% hybridization to the V. aestuarianus type strain, we confirm the affiliation of the cockle and oyster strains to the species V. aestuarianus. Indeed, hybridizations were 70, 71 and 79% with the V. aestuarianus type strain for strains 12/122 3T3, 15/061 4T2 and 02/041, respectively, and were 80 and 80% with the V. aestuarianus oyster strain for cockle strains 12/122 3T3 and 15/061 4T2, respectively (Table 2). Garnier et al. generally found similar results for hybridizations between the V. aestuarianus type strain and different oyster strains: their hybridizations were between 54 and 87% [START_REF] Garnier | Molecular and phenotypic characterization of Vibrio aestuarianus subsp francensis subsp nov., a pathogen of the oyster Crassostrea gigas[END_REF].

The 3 representative strains for the 3 subspecies were further studied by whole cell matrixassisted laser desorption ionization-time of flight mass spectrometry profiling. MALDI-TOF mass spectrum measurements and processing were carried out as described by Mougin et al [START_REF] Mougin | Rapid Identification of Vibrio Species of the Harveyi Clade Using MALDI-TOF MS Profiling With Main Spectral Profile Database Implemented With an In-House Database: Luvibase[END_REF]. The analyses were performed using Flex Analysis (version 3.4) and ClinProTool (version 3.0) softwares (Bruker Daltonics). Figure 2 shows the comparison of the MALDI-TOF spectra. First, all 3 strains were identified with Bruker database as Vibrio aestuarianus with correct scores. Moreover, significant differences observed between spectra (12 per strains) with prominent mass ions at 3877.4 and 7321.7 , characteristic for V. aestuarianus francensis subsp. nov type strain and V. aestuarianus cardii subsp. nov type strain, respectively. The dendrogram based on whole-cell MALDI-TOF mass spectra as well as principal component analysis demonstrated a clear separation of the three strains.

Comparison in other isolates is now needed to confirm the peak combination allowing discrimination [START_REF] Rahmani | MALDI-TOF MS as a promising tool to assess potential virulence of Vibrio tapetis isolates[END_REF].

Genome Features

The existence of three distinct and coherent groups was also supported by results from the comparison of the complete genomes of four cockle isolates, genome of the V. aestuarianus type strain, and the five genomes already available for different V. aestuarianus oyster strains. Briefly, whole genomes of cockle isolates were sequenced by BIOFIDAL (Biofidal, France) after DNA isolation using NucleoSpin microbial DNA (Macherey-Nagel). DNA Libraries were made from 1 ng of genomic DNA using a Nextera XT DNA kit (Illumina). Paired-end isolate sequencing (2x300) was performed on a MiSeq system (Illumina). Genomes were assembled de novo using Spades (3.11.0 version, [START_REF] Bankevich | SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing[END_REF]), with read error corrections and 20x coverage cut-off. Assembly quality was checked using QUAST [START_REF] Gurevich | QUAST: quality assessment tool for genome assemblies[END_REF]. Assemblies were added in MAGE [START_REF] Vallenet | MicroScope: an integrated platform for the annotation and exploration of microbial gene functions through genomic, pangenomic and metabolic comparative analysis[END_REF], excluding contigs with lengths below 1000 bp. Oyster and cockle Vibrio aestuarianus strain draft genomes presented typical features for the species. They ranged between 4.02 and 4.51Mb, with a relatively close number of CDS, genes coding for transfer RNA (between 91 and 102 for cockle isolates and between 84 and 85 for oyster isolates), and G+C contents (42.68 to 42.89% for cockle isolates and 42.52 to 42.83% for oyster isolates) (Table 3). 16S sequences and whole genome comparisons were performed by calculating average nucleotide identity (ANI) between strains (JSpeciesWS, [START_REF] Richter | JSpeciesWS: a web server for prokaryotic species circumscription based on pairwise genome comparison[END_REF]). A distance matrix was calculated inversing ANI results based on whole genomes. Similarity percentage calculated between 16S sequences extracted from whole genomes ranged from 99.09 to 100% (Data not shown). All values were above the 98.7 % cut-off generally accepted for differentiating two species [START_REF] Kim | Towards a taxonomic coherence between average nucleotide identity and 16S rRNA gene sequence similarity for species demarcation of prokaryotes[END_REF]. Moreover, ANI values calculated with whole genomes, indicated an intra-clade similarity of 99.21% to 99.99% (mean value of 99.69%) for V. aestuarianus oyster strains, and an intra-clade similarity ranging from 97.88% to 99.99% (mean value of 98.47%) for cockle isolates. The inter-clade similarity ranged from 96.83% and 97.8% depending on the isolates (Table S5) and confirmed that all strains belong to the same species (ANI>95%) but can be classified into three different subspecies.

Comparison in predicted CDS revealed a core genome, shared by all 3 strains, with 2858 CDS and specific genomes containing 611 to 797 specific CDS (figure S3). As oysters and cockles' isolates are mollusc pathogens, predicted virulome was compared in the 3 representative strains. Globally, the 3 strains contain homologous CDS encoding putative virulence effectors, virulence regulators, secretion systems (SST1, SST2, SST4 and SST6) and pili/flagella/chemotaxis-related proteins. But, some predicted virulence factors are missing in the different strains as haemolysin co-regulated protein in LMG7909T, zonular occludens toxin zot and T1SS-2 in 12/122-3T3 and a complete T6SS locus in 02-041. However, their implication in pathogenesis must definitively be demonstrated.

Physiology and Chemotaxonomy

The phenotypes of selected isolates were characterized in different ways. Bacterial growth was observed on Marine Agar (Conda), TCBS medium (AES Chemunex), and CHROMagar™ Vibrio medium (CHROMagar) (Figure S2). The capacity of bacteria to grow in Zobell broth with 0% to 10% NaCl at 20°C, or with 2% NaCl at 4°C, 10°, 37°C and 40°C was determined (Table 1). NaCl, were suspended in SASW to obtain about OD 600 = 0.1. After inoculation, the strips were incubated for 48 h at 20°C and results were determined according to the manufacturer's instructions.

Antibiograms were carried out using disc diffusion assays on

Finally, FAME (Fatty Acids) determination was performed by the BCCM TM /LMG laboratory at Gent University using a standard commercial identification system (MIDI) and gas chromatography.

The 17 cockle strains and two oyster strains were able to grow in Zobell liquid medium at different temperatures from 4°C to 37°C, with a limited growth rate at the extremes. However, no growth was observed at 40°C, with or without 6% salt. Only the Vibrio aestuarianus type strain was able to grow in the presence of 6% salt. All the strains, like the type strain, were shown to be sensitive to the tested antibiotics (Table 1). Cockle and oyster strains appeared under blue light on Chromagar medium and, unlike the type strain, had a limited growth capacity on TCBS medium, where a few yellow colonies appeared (Figure S2). The observations on oyster strains and the type strain were similar to those made by Garnier et al. [START_REF] Garnier | Molecular and phenotypic characterization of Vibrio aestuarianus subsp francensis subsp nov., a pathogen of the oyster Crassostrea gigas[END_REF]. Additional biochemical characterization analyses were performed to test this proposition, using API20E and 50CH strips (Table 4) and FAME analyses (Table 5). The strains isolated from cockles and oysters, like the V. aestuarianus type stain, could form acids from D-Glucose, Sucrose, D-Trehalose, N-acetyl-D-glucosamine and have Betagalactosidase activity and indole production. Cockle strains could be distinguished from oyster strains by their utilization of D-mannitol and arginine and fermentation of starch, D-ribose and glycogen.

They could be distinguished from the type strain by the non-fermentation of potassium gluconate.

Oyster strains cannot use D-mannitol and arginine. However, as fermentation and utilization of substrates appeared variable among the cockle and oyster strains, no character clearly separated the three groups of strains.

Other characters that can also differentiate strains include their fatty acids (Table 5). V. aestuarianus strains can be distinguished from V. ordalii and V. furnisii by their lower level of 3 hydroxy-dodecanoic acid and higher level of iso-C16:0 fatty acid. Between the V. aestuarianus strains, differences were observed for the C12:0-3-OH, iso-C14:0 and C14:0, for which cockle strains appeared intermediate and oyster strains had a high level compared with the type strain. However, the most remarkable difference concerned the C18:1 w7c fatty acid, present at a higher level in the type strain compared with the oyster (low level) or cockle strains (intermediate level). In contrast, oyster and cockle strains showed a high level of the C16:0 fatty acid compared to the type strain.

Finally, the pathogenicity of nine cockle isolates (12/122 1T1, 12/122 2 T3, 12/122 3T3, 12/122 4T3, 15/061 1T1, 15/061 4T2, 15/064 3T2, 15/075 1T2, 15/075 3T2) on reared cockles was investigated by experimental injections, which also included oyster strain 02/041, pathogenic for oysters, and Vibrio splendidus-related 15/064 1T2, isolated from cockles (Figure 1). Cockles were induced to spawn in September 2015 at the Ifremer hatchery in La Tremblade and larvae were grown according to standard protocols used for C. gigas, except that seawater temperature was maintained at 20°C for cockle larvae [START_REF] Degremont | Relative importance of family, site, and field placement timing on survival, growth, and yield of hatchery-produced Pacific oyster spat (Crassostrea gigas)[END_REF]. Settlement took place from 11 to 16 days post fertilization. All progenies were maintained under controlled environmental conditions in our facilities using unheated UV treated (40 mj/cm 2 ) seawater in a flow-through system. All animals were fed ad libitum.

Following an initial screening of 30 of these animals (no detection of Vibrio aestuarianus on 159 screened colonies, no detection of V. aestuarianus or herpesvirus OsHV-1 DNA and no detection of any other pathogens by histology) and in the absence of mortality over the rearing period, we chose to use this biological material for the controlled infections. Experiments were conducted in May and June 2016 when cockles weighed 4.39 g (s.d. 0.80) and measured 23.28 mm (s.d. 1.52).

Bacteria isolated from diseased cockles as well as controls (V. splendidus-related from cockles 15/064 1T2 and V. aestuarianus oyster strain 02/041) were grown at 20°C for 20 h in marine broth with constant shaking at 40 rpm (Rotator SB3; Stuart). After a SASW washing, bacterial concentrations were evaluated spectrophotometrically at an optical density (OD) of 600 nm, adjusted to 0.4 and 0.04 and checked by plating. Cockles were anesthetized for 2 to 4 hours in a magnesium chloride solution (MgCl 2 , Fluka) at a final concentration of 50 g l -1 (1/4: v/v seawater/freshwater) under aeration. Subsequently, 50 μl of bacterial suspension (10 7 CFU to 10 6 per cockle) was injected into the adductor muscle. A group of cockles were injected with SASW to serve as negative controls. After injection, the cockles were transferred to tanks (triplicates of 10 cockles in 2 litres) filled with UV-treated, 1-μm-filtered seawater and maintained under static conditions at 20°C with aeration. Mortality was monitored twice a day and any newly dead cockles were removed from the tanks over a 7-day period. Animals were considered to be dead when the valves did not close following stimulation.

After 24 h, cumulative mortalities ranged from 0 to 83% for cockles injected with V. aestuarianus cockle strains. Mortalities continued over 7 days, with a peak during the first two days. For all moribund cockles, the predominant bacterial isolates were obtained from muscle and haemolymph (205 strains screened from 116 animals, all belonging to V. aestuarianus species using the dnaJ-qPCR diagnosis assay). We also noticed the high levels of bacteria belonging to the Splendidus clade in 31 individuals. No mortality was noted in the control or in cockles injected with V. splendidus-related 15/064 1T2 or V. aestuarianus oyster strain 02/041. Garnier et al. [START_REF] Garnier | Molecular and phenotypic characterization of Vibrio aestuarianus subsp francensis subsp nov., a pathogen of the oyster Crassostrea gigas[END_REF] also tested the pathogenicity of 11 Vibrio aestuarianus oyster strains and V. aestuarianus type strain by injection of 10 8 bacteria per oyster. They observed mortality ranging between 2 and 55% depending on the oyster strain, whereas a low mortality (10%) was observed with the V. aestuarianus type strain. The same observation was made by Goudenège et al. [START_REF] Goudenege | A single regulatory gene is sufficient to alter Vibrio aestuarianus pathogenicity in oysters[END_REF], who classified oyster strains as either highly virulent for oysters (i.e. inducing > 50% mortalities at 10 2 CFU animal-1) or non-virulent (i.e. inducing < 50% mortalities at 10 7 CFU animal-1).

In summary, Vibrio aestuarianus isolates obtained from oysters or cockles present phenotypic and genetic traits distinct from the V. aestuarianus type strain and these differences support the separation of these oyster and cockle strains into two novel subspecies of V. aestuarianus distinct from V. aestuarianus LMG 7909T and other V. aestuarianus strains. Therefore, we propose that Vibrio aestuarianus should be divided into three subspecies, for which the names Vibrio aestuarianus subsp. aestuarianus subsp. nov., Vibrio aestuarianus subsp. francensis subsp. nov. and Vibrio aestuarianus subsp. cardii subps. nov. are proposed, with strains LMG 7909T, 02/041 and 12/122 3T3 as the type strains, respectively.

Protologue Emended description of Vibrio aestuarianus

The characteristics modified or added to those reported in original descriptions by [START_REF] Tison | Vibrio-Aestuarianus -a New Species from Estuarine Waters and Shellfish[END_REF] [START_REF] Tison | Vibrio-Aestuarianus -a New Species from Estuarine Waters and Shellfish[END_REF] are as follows:

Vibrio aestuarianus is a gram-negative, oxidase-positive, motile bacterium. Colonial morphology is variable according to the subspecies. Growth occurs on Marine Agar from 4 to 30°C with an optimum between 20 and 25°C and strains do not grow at 40°C or in the absence of NaCl. Most strains utilize N-acetyl-D-glucosamine, Beta-galactosidase and produce indole. Growth occurs from 4 to 30°C, difficult at 37°C and not at 40°C. Strains can use D-mannitol and do not ferment potassium gluconate. Concerning its fatty acid composition, it presents a high level of C16:0 and an intermediate level of other fatty acids. Some strains, including the type strain 12/122 3T3, are pathogenic to the edible cockle Cerastoderma edule but the type strain is not pathogenic to Pacific oyster Crassostrea gigas.

The type strain is 12/122 3T3 (LMG 31436 = DSM 109723, Genbank genome accession number JAAKZK000000000) and was isolated from a diseased cockle sampled from a wild bed (Baie de Somme) in France.

  Mueller-Hinton agar (Oxoid) with commercial impregnated discs (Oxoid): Flumequine 30 µg (UB30), Kanamycin 30 µg (K30), Tetracycline 30 µg (TE30), Sulfamethoxazole/trimethoprim 19:1 25 µg (SXT25), Streptomycin 25 µg (S25), Sulfonamides 300 µg (S3), Erythromycin 30 µg (E30), Penicillin 10 UI (6 µg) (P10), Trimethoprim 5 µg (W5), Chloramphenicol 30 µg (C30), Ampicillin 10 µg (AMP10), Gentamicin 30 µg (CN30) or the Vibriostatic agent, Pteridine 150 µg and 10 µg (O129/150 and O129/10). Inhibition was observed after 24-48 h at 22 °C. API 20E, API 20NE and API 50CH (Biomérieux) were used in accordance with manufacturer's recommendations. Fresh colonies, obtained from culture bacteria on TSA (Difco ™ ) adjusted to 2%

  Strains ferment the sucrose, Dmannose, D-fructose, D-glucose, D-galactose and D-trehalose. Strains are susceptible to Ampicillin, Erythromycin, Chloramphenicol, Kanamycin and Tetracycline. The major fatty acids are C12:0 (2.2-4.4%), C12:0; 3-OH (1.3-3.4%), C14:0 (4-8.7%), iso-C16:0 (4.1-5.6%), C16:0 (13.4-20.06%) and C18:1 w7c (11.3-22.6%). The genome size is 3.9-4.5 Mb (Genbank genome accession number JAAZTU000000000), and genomic DNA G+C content is 42.5-42.9%. The type strain is ATCC 35048 (=DSM 19606, CAIM 592, CCUG 28583, CECT 625, CIP 102971, LMG 7909, NBRC 15629, NCMB 2236). The GenBank/EMBL/DDBJ accession number for the genome of the type strain is JAAZTU000000000. Description of Vibrio aestuarianus subsp. aestuarianus subsp. nov. Vibrio aestuarianus subsp. aestuarianus Tison and Seidler [25] (a.es.tu.ri.a'nus. N.L. masc. adj. aestuarianus, pertaining to an estuary). The description is essentially the same as that given above for the species Vibrio aestuarianus, with the following additions. Colonies are opaque, 2-3 mm in diameter on marine agar (MA) after 24h at 20°C. It grows on TCBS, with 6% salt and at 37°C. It has gelatinase and arginine dihydrolase activity and utilizes starch, Dribose, D-maltose, D-cellobiose, glycogen and potassium gluconate. Among the three subspecies, it presents the highest level of C18:1 w7c fatty acid and the lowest level of C12:0; 3-OH and iso-C14:0 fatty acids. The type strain does not induce mortality of Pacific oyster, Crassostrea gigas. The genome size is predicted around 4.2 Mb and genomic DNA G+C content is 42.6% (Genbank genome accession number JAAZTU000000000). The type strain is ATCC 35048 (=DSM 19606, CAIM 592, CCUG 28583, CECT 625, CIP 102971, LMG 7909, NBRC 15629, NCMB 2236). Description of Vibrio aestuarianus subsp. francensis subsp. nov. Vibrio aestuarianus subsp. francensis as proposed by Garnier et al. [26] (fran.cen'sis. N.L. masc. adj. francensis: of, or belonging to, France) It displays typical characteristics for the species V. aestuarianus. Colonies are white and slightly opaque on Marine Agar and are less than 1 mm in diameter after 24 h at 20°C. No strain can grow in the absence of NaCl or in the presence of 6% NaCl (except for strain 02/041) or on TCBS. Growth occurs from 4 to 30°C. Most strains have no arginine dihydrolase and cannot use D-mannitol. Some strains use citrate. It presents the highest level of C16:0 and C14:0 fatty acids and the lowest level of C18:1 w7c fatty acids compared to the two other subspecies. The majority of strains, including the type strain 02/041, are pathogenic to the Pacific oyster, Crassostrea gigas, but the type strain is not pathogenic to edible cockles, Cerastoderma edule. The genome size is estimated around 4.18-4.53 Mb with the draft genomes, and genomic DNA G+C content is about 42.52-42.83 mol%. The type strain is 02/041 (LMG 24517 = CIP 109791, Genbank genome accession number JAAZTV000000000) and was isolated from a diseased oyster maintained in the Argenton experimental hatchery (Ifremer, France). Description of Vibrio aestuarianus subsp. cardii subsp. nov. Vibrio aestuarianus subsp. cardii (car'di.i. N.L. fam. cardii referring to Cardium -latin: cockle) It displays typical characteristics for the species V. aestuarianus. Colonies are around 1 mm in diameter after 24 h at 20°C. No strain can grow in the absence of NaCl and growth was poor on TCBS.

Table 2 .

 2 DDH values measured by comparing Vibrio aestuarianus strains from cockles and oysters with V. aestuarianus type strain (LMG 7909T) and other 1 close species 2 Taxa. 1, V. aestuarianus subsp. cardii subsp. nov. 12/122 3T3; 2, V. aestuarianus subsp. cardii subsp. nov. 15/061 4T2; 3, V. aestuarianus subsp. francensis 3 subsp. nov. 02/041; 4, V. aestuarianus subsp. aestuarianus subsp. nov. LMG 7909T; 5, Vibrio furnisii LMG 7910T; 6, Vibrio ordalii LMG 13544T; 7, Vibrio 4

	5	salilacus JCM 19265T.
	6
			1	2	3	4	5	6	7
		V. aestuarianus subsp. cardii subsp. nov. 12/122 3T3	100
		V. aestuarianus subsp.cardii subsp. nov. 15/061 4T2	84 (7) 100
		V. aestuarianus subsp.francensis subsp. nov. 02/041	80 (18) 80 (11) 100
		V. aestuarianus subsp.aestuarianus subsp. nov. LMG 7909T	70 (2) 71 (5) 79 (22) 100
		V. furnisii LMG 7910T	14 (3) 18 (4) 16 (5) 17 (20) 100
		V. ordalii LMG 13544T	14 (12) 16 (5) 16 (5) 9 (16) 8 (19) 100
		V. salilacus JCM 19265T	15 (4) 16 (4) 12 (2) n.d.	n.d.	n.d. 100

Table 5 .

 5 Fatty acid compositions (as percentage of the total fatty acids present) in Vibrio 18 aestuarianus strains from cockles and oysters, V. aestuarianus type strain and related species.[START_REF] Romalde | New Vibrio species associated to molluscan microbiota: a review[END_REF] Taxa. 1, V. aestuarianus subsp. cardii subsp. nov. 12/122 3T3; 2, V. aestuarianus subsp. francensis 20 subsp. nov. 02/041T; 3, V. aestuarianus subsp. aestuarianus subsp. nov. LMG 7909T; 4, V. furnisii 21

	22	LMG 7910T; 5, V. ordalii LMG 13544T.
	23	nd = no data
	24
			1	2	3	4	5
		C12:0	4.4 (0.4)	4	3.2	6.7	-
		C12:0; 3-OH	2.1 (0.4)	3.4	1.3	6.1	4.5
		iso-C14:0	1.2 (0.1)	2.1	TR	-	-
		C14:0	4.8 (1.7)	8.7	4	4.5	5.9
		iso-C16:0	4.1 (1.4)	5.4	5.6	TR	-
		C16:0	18.6 (1.5)	16.3	13.4	12	19.2
		C18:1 w7c	14.5 (2.7)	11.3	22.6	17.1	23.4
		summed features 2	2.9 (1.6)	2.3	1.1	6.7	1.5
		summed features 3	41.6 (1.8)	34.6	41.8	37.5	39.8

Acknowledgements

We wish to thank the coastal correspondents of the REPAMO network (particularly those at Boulogne sur Mer, Port en Bessin and Dinard) for their contribution to collecting and providing samples and related information. We also thank the hatchery and nursery teams of the Laboratory of Genetics and Pathology at Ifremer La Tremblade and Ifremer Bouin for their assistance with cockle production. We are also grateful to M. Thuillier and M. Le Piouffle from the laboratory LABOCEA, Mme Cauvin and Mme Benoît from the laboratory, LABEO-Manche and M. Treilles from the laboratory QUALYSE for their assistance in MALDI-TOF analyses. ABBREVIATIONS ANI: Average Nucleotide Identity DDH: DNA-DNA hybridization FAME: Fatty Acid Methyl Esters MALDI-TOF MS : Matrix-Assisted Laser Desorption Ionization-Time Of Flight Mass Spectrometry REPAMO: REseau de PAthologie des MOllusques SASW: Sterilized Artificial SeaWater TCBS: Thiosulfate Citrate Bile salts Sucrose agar TSA: Tryptone Soya Agar CDS: coding DNA sequence

Funding information

This work was supported by Ifremer; by DGAL (Direction Générale de l'Alimentation) through the REPAMO network and the National Reference Laboratory for Mollusc Diseases, Ifremer, La Tremblade; by DPMA (Direction des Pêches Maritimes et de l'Aquaculture) through MIRAGE project funding and by the European Commission through the European Union Reference Laboratory, Ifremer, La Tremblade. This work was also partly supported by the EU-funded project VIVALDI ('Preventing and Mitigating Farmed Bivalve Diseases') H2020 program, no. 678589.

AUTHOR STATEMENTS 1.6 Authors and contributors

DT, CD, AGS, BC, YG, DS, CG, AM and MAT performed the experiments. LD produced and conserved the animals. MAT and CG supervised the project. DT, CD, AM, CG and MAT analysed the data. MAT, AM and CG wrote the manuscript that was revised and approved by all the authors.

Conflicts of interest

The authors declare that there are no conflicts of interest.

FIGURES AND TABLES