Lina Aliouat
email: lina.aliouat@univ-fcomte.fr

Hakim Mabed
email: hmabed@femto-st.fr

Julien Bourgeois
email: julien.bourgeois@femto-st.fr

Layer Protocol for Traffic Regulation in Ultra-Dense Nanonetwork Traffic Regulation in Ultra-Dense Nanonetwork

Keywords: Terahertz Nanonetwork, Ultra Dense nanonetwork, directional antenna, MAC Layer

The nano terahertz networks represent one of the promising areas in the field of wireless telecommunications. Technological advances in miniaturization of antennas and terahertz communications have paved the way for new network applications such as the body network, the programmable material and multi-core processors. Some of these applications require the concentration of a very large number of tiny nodes in a limited space. In this ultra-dense context and in the absence of centralized access control units, we propose to implement a distributed strategy of spatial and temporal traffic regulation to guard against the risks of congestion, interference and energy over-consumption. In this paper, we propose a protocol for optimizing terahertz radio links using beam steering antenna, distributed time division technique and sleep mode in order to reduce the flow of redundant traffic over the network, smooth the volume of communications exchanged over time, and preserves the lifetime of the nodes.

Introduction

New network applications have emerged in recent years driven by major advances in the miniaturization of electronic devices and radio antennas. In ultra dense nanonetworks, a very large number of radio devices are confined in a small space. In this context, the Terahertz frequency band has the double advantage of combining a high bandwidth with a low energy and coverage range. The use of short-range terahertz communications finds many applications in the field of massively multi-core computer architectures [START_REF] Abadal | Grapheneenabled wireless communication for massive multicore architectures[END_REF] and programmable material [START_REF] Bourgeois | Programmable matter as a cyber-physical conjugation[END_REF].

Due to the limited computation and energy capabilities of the nanonetwork nodes (sub-millimeter scale), the multiple access protocols must meet simplicity and scalability requirements. Access to the channel must be done with a reduction of the number of control messages and without resorting to centralized entities. To this end, several innovative techniques have been proposed such as: PHLAME [START_REF] Jornet | Sol-Pareta: PHLAME: A Physical Layer Aware MAC protocol for Electromagnetic nanonetworks[END_REF], ASRH-TSOOK [START_REF] Mabed | Enhanced spread in time on-off keying technique for dense Terahertz nanonetworks[END_REF], HLMAC [START_REF] Jianling | High-throughput low-delay MAC protocol for TeraHertz ultra-high data-rate wireless networks[END_REF], DRIH-MAC [START_REF] Mohrehkesh | DRIH-MAC: A Distributed Receiver-Initiated Harvesting-Aware MAC for Nanonetworks[END_REF], etc. However, in view of the extreme density of the network, the classical multiple access protocols are not sufficient to spread the traffic load and to control the multi-hop flows on the network. Indeed, given the density of the network, the message broadcast causes numerous feedback loops that saturate the system (broadcast storms).

The traffic regulation protocol is seen as a 2.5 networking layer that allows to extend the access control layer missions with some routing considerations. The traffic regulation amounts to defining a logical topology of the network starting from the physical topology where any two nodes can communicate when they are within range of each other. The logical topology designates a subset of neighboring nodes which can communicate in a predetermined direction and at predetermined time. Formally, the logical topology represents a directed sub-graph of the physical topology. The logical topology is said robust when the directed sub-graph is strongly connected. The connectivity degree of the sub-graph could be used as a measurement of the logical topology robustness. A directed graph is said k-connected if it remains connected whenever fewer than k nodes are removed.

Traffic regulation protocols for ad hoc networks have been widely studied in the literature [START_REF] Jacquet | Optimized Link State Routing Protocol for Ad Hoc Networks[END_REF][START_REF] Rodolakis | Multicast Overlay Spanning Tree Protocol for Ad Hoc Networks[END_REF][START_REF] Radhakrishnan | Protocol for dynamic ad hoc networks using distributed spanning trees[END_REF]. One of the best known is the Optimized Link State Routing Protocol (OLSR) [START_REF] Jacquet | Optimized Link State Routing Protocol for Ad Hoc Networks[END_REF]. However the adaptation of this protocol in the case of an ultra-dense network is complicated. OLSR is based on the exchange of neighboring lists between nodes. Due to the density of the nanonetwork, those lists are heavy to transmit and difficult to store or to process. Other protocols aim to define a spanning tree over the network's nodes [START_REF] Rodolakis | Quality-of-Service Multicast Overlay Spanning Tree Algorithms for Wireless Ad Hoc Networks[END_REF]. The logical topology of the network represents then a tree where the nodes close to the root concentrate more traffic then nodes near the leafs. Therefore these protocols are adapted when the physical structure of the network involves different types of nodes : simple nodes and super-nodes like in Wireless Body Sensor Nanonetwork architecture [START_REF] Afsana | An Energy Conserving Routing Scheme for Wireless Body Sensor Nanonetwork Communication[END_REF]. In addition, such protocols present only one path to link every two nodes, which makes the logical topology unreliable.

By conclusion, few works from literature deal with traffic regulation in dense homogeneous ad hoc networks.

Our proposed protocol has three main objectives. First, traffic load evolution presents peaks that lead to congestion phenomena. The principle of communication by appointments allows to spread the traffic and to schedule communications in time. Secondly, dense networks present various path to transmit data from one node to another. In terms of routing (layer 3), this implies a greater complexity of choice. The risk of local congestion on the network is therefore higher and more redundancy is expected (multiple receptions of the same message by the same node). The use of electronically steerable antennas improves the control of the transmitted radio signals (interference). Moreover, each node selects the subset of neighboring nodes with which it can communicate directly. Finally, communication by appointment allows nodes to plan their waking and sleeping periods. In addition, the energy consumed by communications is reduced because the emission power is channeled in specific directions at given periods. Finally, only a subset of covered nodes are selected as sources or successors nodes and messages of not selected sources are ignored (not delivered to the layer 3).

In this paper, we propose an original procedure for the Layer 2.5 networking protocol that takes into account the terahertz frequencies particularities in a dense context. The idea is to extract a logical topology from a dense homogeneous nanonetwork that allows to both reduce the number of direct communication links and maintain the robustness against temporal nodes unavailability. This procedure exploits the available antenna steering techniques to schedule over time and space the data transmission. Unlike the other approaches, our method does not impose any conditions on the physical platform and presents, according to our knowledge, the first 2.5 networking layer protocol adapted for terahertz nanonetworks.

Traffic regulation problem modeling

Let W be a nano wireless network composed of N nodes. Each node in the network has a reconfigurable directional antenna that can be steered dynamically to cover a particular direction. Let T ch be the time needed to change the orientation of an antenna. Let G(X, A) the connected graph describing the physical topology of the network with X the set of nodes (|X|= N ) and A the communication links between the nodes. (x, y) ∈ A means that there is a particular configuration of the x and y antennas that makes the two nodes communicate directly.

Traffic regulation constraints

The traffic control problem consists of calculating a directed sub-graph G (X, E) with (x, y) ∈ E → (x, y) ∈ A where the following conditions are satisfied:

-G is strongly connected: given two nodes in the graph, there is a way to route the data from one node to the other in the two directions.

(1) ∀x, y ∈ X 2 , ∃ a path from x to y -G is robust: whatever the node, there are enough ways, p, to receive the data from the other nodes and enough means, s, for the node to broadcast its own data. The values p and s denote the desired level of robustness represented by the number of predecessors and successors of each node. Choosing a large value of p and s allows a higher level of robustness that derives from the reliability level of the nodes. When nodes are prone to a high risk of outages or if the energetic capacity of nodes makes it regularly in charging phase, then a high value of p and s is more suitable. (2) T c = N S × T s For a given period S i , a node x is either in listening, transmitting or sleeping mode. For each listening period S i of x, there is one and only one node y ∈ X such that (y, x) ∈ E, and S xy = S i , which means that there is only one listened node at a time. If the period S i is a transmission period of x, then there is at least one node y such that directed edge (x, y) ∈ E and S xy = S i . Sleeping mode corresponds to periods of time where the node x has no directed edge (x, y) ∈ E or (y, x) ∈ E with S xy = S i , which means that the node is not listening and not transmitting.

Given the asynchronous nature of the network, the index of a period S i of a given node has only a local signification. The figure 2 When two edges (x, y1) and (x, y2) with the same tail have the same index value on x, the two incoming nodes y1 and y2 are then served by the same multicast stream. By the way, the edges (x1, y) and (x2, y) with the same head can not have the same period index on y in order to avoid interference. For all the edges (x, y) ∈ G , the associated listening period on y must be equal to the duration of the transmitting time on x in order to maximize the sleep periods.

Along with the respect of traffic regulation constraints, in particular the connectivity and the robustness of the sub-graph G , the traffic control algo-rithm must take care to maximize the useful listening and transmission times, T i , (T i ⊂ S i ) as well as the sleep periods of the nodes. A transmission period of a node x is said to be useful when throughout all its duration, all nodes y, (x, y) ∈ E are at listening mode. A listening period of a node y is said useful when throughout its duration, the listened node (there is only one) is in transmission phase to y.

Traffic regulation in Terahertz network

In DAMC modulation technique [START_REF] Han | Distanceaware multicarrier (DAMC) modulation in Terahertz band communication[END_REF], the terahertz frequency band is mainly subdivided into three frequency windows which are allocated according to the transmitter-to-receiver distance which is either short, mid or long. To take into account this particularity, the traffic regulation protocol associates with each active period a coverage range: short, mid or long. Then, only nodes with a distance in the selected range are considered to be successors. At a given transmitting time slot, the successors of a given node are all in the same range, allowing to use the same optimal frequency window to serve them (See Figure 3). Fig. 3 At every period, the transmitting node selects its successors according their distance from the node in order to optimize the used frequency window. The colored arcs represent the orientation and coverage of the node at a given period.

Distributed algorithm of traffic regulation

The algorithm 1 represents our traffic regulation protocol. The design of the algorithm aims to satisfy two main constraints: a reduced computation requirement and limited messages exchange. When a node wants to join the nanonetwork, it defines for each slot a covering range (short, mid or long). Then, the node alternates between two modes. In the first mode, the node listens to the channel and switches over the set of configurations C = {C 1 ..C N S } with a frequency of 1/T s . In the second mode, the node launches invitations in different directions looking for successor nodes and changes its configuration with a frequency of 1/T c . The use of two different reconfiguration speeds aims to prevent the hidden node problem.

In successor search mode, the node keeps the same antenna configuration during all a cycle T c of N S slots T s . The value of N S depends on several parameters such as the reconfiguration delay T ch , the cycle duration and the number of needed sources and successors. In our tests, we have set N S = 14. At each period T I , the node launches invitations. After each reception of an acceptance, the source updates its useful transmission period, its coverage list and sets the period mode to 'transmission' mode. Nodes in source search mode, listen for any invitations. Based on the target range of the source, the strength of the received signal and the remaining listening time, the node chooses to accept the invitation or not. In case of acceptance, the node registers the useful listening period. Once the number of necessary source nodes is reached (parameter p), the node stops using the source search mode.

To avoid a slow start when few nodes are active, the number of listening and announcement cycles is limited to M . After every M ordinary operating cycles, a node with not enough sources (resp. successors) listens (resp. resends invitations) during an entire cycle. These two procedures allow nodes that start very early or late compared to the other nodes to complete their lists of prefixed links.

The disconnection of the graph G (X, E) is avoided thanks to a long-term procedure which provides that each node, at a very important time interval, listens in all directions whether neighboring nodes belong to other connected components. For this purpose, the nodes of the same connected component share an identifier of the component which corresponds to the smallest MAC identifier of the nodes belonging to the component. When a node detects a neighboring node with a different component identifier, a connection procedure of the two nodes is started.

Tests and results

To study the impact of the traffic regulation algorithm over dense nanonetworks, we have established several test scenarios, which are differentiated by the number of nodes, the spatial dimension of the network, the range of the radio signal and the density of the network. We first begin by evaluating the impact of traffic regulation on network performance in terms of interference. A first indicator of the impact on interference is the comparison of the number of edges in the graph G(X, A) and the number of edges in the graph G (X, E). When |E|<< |A|, the average number of signals arriving on each node is reduced considerably. Interference reduction also benefits from time division access mode, beam control antenna and selectivity of listened sources. All these factors make it possible to reduce the risks of massive arrival of communications at the same time on the same node. For each scenario, on the left, the graph G(X, A) shows the physical topology of the network and on the right the graph G (X, E) is obtained by the traffic control algorithm 1. The traffic regulation is carried out with the maximum number of sources equal to p = 4 and the maximum number of successors equal to s = 10. For Scenario 1, the traffic regulation algorithm reduces the density of the graph from 1202 edges to 744 edges, i.e. a reduction of 38%. For the Scenario 3, the traffic control algorithm reduces the graph from 30891 edges to 3997 edges, which corresponds to a reduction of 87%.

We also assess the impact of traffic regulation on network data broadcasting. To this end, we have selected three evaluation criteria: The total number of receptions including redundancies, the maximum number of receptions of the same message on one node and the time for a total broadcasting of the message. According to the first two criteria, the traffic regulation allows to improve the behavior of the network when a node broadcasts a message over the network. Regarding the total number of messages received by the nodes, an approach without traffic regulation will generate for the scenarios 1, 2 and 3 respectively 2404, 15082 and 61618 messages. The number of receptions with traffic control decreases to 744, 1997 and 3997 messages respectively. The maximum number of reception of the same message varies in the approach without traffic regulation between 24, 45 and 96 for the three scenarios whereas with traffic control, it remains 4 for all three scenarios. This corresponds to the value of the maximum number of sources: parameter p.

The number of hops allowing all nodes to receive the message is also a crucial factor for the performance of the traffic control algorithm. The simulations show that the total diffusion of a message depends on the position of the source node in the network. We studied 4 different placement of the source nodes for the Scenario 1, shown in Figure 4, circled in red and labelled A, B, C and D. Broadcasting without traffic regulation from node (A) requires 9 hops while 14 hops are required with traffic control. Similarly, the application of traffic control requires respectively 26, 18 and 22 hops to broadcast a message from the nodes (B), (C) and (D) instead of 12, 12 and 9 jumps without regulation traffic. For Scenario 1 which is the least dense, the number of hops for the total broadcast doubles with traffic regulation, a negligible additional cost in return for reducing the total number of exchanged messages.

Conclusion

In this paper, we have proposed a novel distributed protocol for optimizing the logical topology of homogeneous, ultra-dense terahertz nanonetworks. The objective of this optimization is to reduce the amount of interferences due to the concentration of numerous nodes transmitting in all directions at random instants. The logical topology we proposed fixes for each node the moments during which it can send data in a given direction to given selected destination nodes. Each node ignores messages sent by nodes that do not belong to the predefined subset of neighbors, which reduces the energy consumption.

Furthermore, the reduction of direct communication links between nodes does not impact data broadcasting coverage. Simulations show that the diminution of the broadcasting speed is limited compared to the gain in terms of transmission redundancy and generated interferences.

Finally, the traffic control protocol is adapted to the DAMC protocol used for terahertz communication. The successors of a given node in a specific direction are in the same range of distance from the transmitter. Therefore, all the nodes covered in the same direction are served by the same frequency channel. The distance between nodes and their successors varies according to the transmission direction leading to a better use of the terahertz band.

Fig. 1

 1 Fig. 1 Traffic regulation problem: from physical to logical topology

  shows an example of traffic regulation involving 5 nodes. The node (A) has two active periods S 1 in transmission and S 2 in reception. The period S 1 = S AC = [0.3 -0.4] covers 1/10 of the cycle time T c between instants 0.3 × T c and 0.4 × T c . During this period, the node (A) covers the node (C) which is listening during its period S 3 = S CA = [0 -0.1] as well as the node (E) which listens during its period S 1 = S EA = [0.7 -0.8].

Fig. 2

 2 Fig.2Example of TDMA synchronization: each edge (x, y) is indexed by the index of the transmission period on x (S i ) and the index of the reception period on y (S j ). Each period S[tb -te] is designated by its beginning and end time (the time is given relatively to the concerned node). The colored arcs display the antenna orientation at the corresponding period.

Fig. 4

 4 Fig. 4 Traffic regulation for 3 scenarios: 200 nodes, 500 nodes and 1000 nodes. Tests are implemented using Microsoft Excel VBA.

Figure 4

 4 Figure 4 presents three simulations of traffic regulation using different number of nodes 200, 500 and 1000 in the same area using Microsoft Excel VBA. For each scenario, on the left, the graph G(X, A) shows the physical topology of the network and on the right the graph G (X, E) is obtained by the traffic control algorithm 1. The traffic regulation is carried out with the maximum number of sources equal to p = 4 and the maximum number of successors equal to s = 10. For Scenario 1, the traffic regulation algorithm reduces the density of the graph from 1202 edges to 744 edges, i.e. a reduction of 38%. For the Scenario 3, the traffic control algorithm reduces the graph from 30891 edges to 3997 edges, which corresponds to a reduction of 87%.We also assess the impact of traffic regulation on network data broadcasting. To this end, we have selected three evaluation criteria: The total number of receptions including redundancies, the maximum number of receptions of the same message on one node and the time for a total broadcasting of the message. According to the first two criteria, the traffic regulation allows to improve the behavior of the network when a node broadcasts a message over the network. Regarding the total number of messages received by the nodes, an approach without traffic regulation will generate for the scenarios 1, 2 and 3 respectively 2404, 15082 and 61618 messages. The number of receptions with traffic control decreases to 744, 1997 and 3997 messages respectively. The maximum number of reception of the same message varies in the approach without traffic regulation between 24, 45 and 96 for the three scenarios whereas with traffic control, it remains 4 for all three scenarios. This corresponds to the value of the maximum number of sources: parameter p.The number of hops allowing all nodes to receive the message is also a crucial factor for the performance of the traffic control algorithm. The simulations show that the total diffusion of a message depends on the position of the source node in the network. We studied 4 different placement of the source nodes for the Scenario 1, shown in Figure4, circled in red and labelled A, B, C and D. Broadcasting without traffic regulation from node (A) requires 9 hops while 14 hops are required with traffic control. Similarly, the application of traffic control requires respectively 26, 18 and 22 hops to broadcast a message from the nodes (B), (C) and (D) instead of 12, 12 and 9 jumps without regulation traffic. For Scenario 1 which is the least dense, the number of hops for the total broadcast doubles with traffic regulation, a negligible additional cost in return for reducing the total number of exchanged messages.

Require: Ts, Tc, t0, N S, P = {C 1 ..C N S }, T = ∅, nbcycles = 0, t0 = now(), nbsec = 0, nbsrc = 0 1: for i ∈ {1 to NS} do 2:

rng i = rand(1..


CEA Marcoule Murielle Bertrand

Sophie Lalleman

L.-Z Mojica-Rodriguez

Christian Sorel

Philippe Moisy

F Auger

E Plasari

H Muhr

Esat

Luz-Adriana Mojica-Rodriguez

Applied thermodynamics in nuclear industry modelling of precipitation processes

HAL is

PRECIPITATION PROCESS MODELLING Precipitation Modelling :

Coupling Chemistry/Turbulence

  ) ( , ) , ( ) , ( ) , ( 1 L N R V t L Q t L Q L t L G t t L V V i i o           

Empirical approaches

Bromley method or SIT specific interaction theory Exp. data (Rard and Miller, 1979) Davies and Jones' model (Söhnel and Garside, 1992) Bromley's model (Bromley, 1973)

                      ν ν 1 ν eq , X ν eq , M ν X ν M z z z z a a a a S O .xH X M X ν M ν 2 ν ν O H - z z 2          

Good compromise J:

       3 3 3 3 3 3 NO Nd NO Nd ) Nd(NO δ δ B B B        Cl Nd Cl Nd NdCl δ δ B B B

  Calculation of the driving force = supersaturation S during the precipitation process In concentrated solutions involving salts and ions of high valence  Coefficient activity models:Physico-chemical approachesPitzer equations



  very low solubility  widely used in the Ln and An separation, analytical chemistry and treatment of high-level liquid wastes -in nitric acid medium Bromley method for multicomponent systems Based on the contributions B M1X1 of every electrolytes present in solution  Neodymium individual contributions Tabulated by Bromley in 1973 After  New experimental binary data published later for neodymium nitrate and chloride  The Bromley model not in good agreement  Determination of new values for neodymium contributions (B + ,δ + ) using tow salts:

Nd

  Very good agreement between exp. data and the Bromley model combined with the contributions of HNO 3 and HCl New contributions of Nd 3+ in the Bromley model 7 [1] Rard J.A. and Miller D.G, 1979. J. Chem. Eng. Data 24, 348-353 [2] Spedding et al., 1976. J. chem. Eng. Data 21, 341-360. 1 st STEP: Validation of the Bromley individual contributions for NO 3 Experimental data (Charrin, 1999) Bromley model with B = 0.0776 (Bromley, 1973) Experimental data (Hamer and Wu, 1972) Bromley model with B = 0.1433 (Bromley, 1973) (Rard and Miller, 1979) Bromley's model with the new B value (Bromley, 1973) (Spedding et al., 1976) Bromley's model with the new B value (Bromley, 1973) Very good agreement between exp. data and the Bromley model with the new contributions of Nd(NO 3 ) 3 et NdCl 3 These new contributions  application to ternary solutions Nd(NO 3 ) 3 /HNO 3 /H 2 O  Measurements of water activity a W in ternary mixtures Using the concept of simple solutions and Mikulin equation*: a W  activity coefficients 𝛾 𝑁𝑑(Precipitation in sulphuric media UO 2 SO 4 / H 2 SO 4 / Na 2 SO 4 / H 2 O 2 / H 2 O Bromley method not adapted: sulphuric acid partially dissociated + no data for hydrogen peroxide  Specific Interaction Theory* ** Based on specific interaction coefficients between ions of opposite charge  1 st step: Validation of the specific interaction coefficients of ions: 2 𝑂 ↔ 𝑈𝑂 2 (𝑂 2 ) • 4𝐻 2 𝑂 *T. Vercouter, B. Amekraz, C. Moulin, E. Giffaut and P. Vitorge, Inorg. Chem., 2005, 44, 7570-7581Bromley L.A., 1973. AIChE J. 19, 313-320 **R. M. Cigala, F. Crea, C. De Stefano, G. Lando, D. Milea and S. Sammartano, Geochim. Cosmochim. Acta, 2012, 87, 1between experimental pH and SIT  2 nd step: determination of the specific interaction coefficients of other ions from experimental data of uranium peroxide solubility (ICP-AES and ICP-MS) And by application of the SIT hypothesis: Commutative function (i,j) = (j,i) Independent from I (i,j) = 0 for ions of same charge For unknown coefficients  analogies between ions of same charge ε (kg.mol -12 ± 0,06 → -0,18 -0,12 ± 0,06 → -0,12 UO 2 (SO 4 Comparison between experimental and calculated solubilities over a wide range of operating conditions Good agreement higher deviations at low sulphate concentrations  very low solubilities  high analytical uncertainties Conclusion and outlooks 14 These methods can be applied to many other systems provided -the binary data of 2 salts are known for the Bromley model -solubility measurements can be performed for SIT  For precipitation process modelling  development of thermodynamic models to calculate the supersaturation taking into account the nonideality of the solutions  On the basis of simple empirical methods  Advantage  easy to use and accurate  Application  Experiments New individual contributions of the neodymium ion for the Bromley model New specific interaction coefficients for SIT  Comparison between experiments and models  a good agreement Scientific team

French Atomic Energy Commission


Marwa Haddad

Jean-Marc Nicod

Marie-Cécile Péra

Hydrogen infrastructure: data-center supply-refueling station synergy

Nowadays, data-center are available in a lot of cities in different countries around the world. Their continuous growth in size and complexity due to the increasing demand for storage, networking and computation has become a world wide energetic problem. Taking into account their huge consumption of electricity and the amount of carbon emission in the air, researchers work on making data-centers as green as possible using renewable energy sources. On the other hand, electrical cars, more and more common as they have the potential to provide a high autonomy and to reduce significantly the emmissions that contribute to climate change. These green cars have however an inconvenient which is the availability of hydrogen station to fill their tank. Thus, this paper presents first a sizing of a hybrid system to power a data-center of 500 kW composed of photovoltaic panels, winding turbines, batteries and fuel cells. Consequently, this paper presents an option for this hydrogen storage to become easily a secure hydrogen station to supply fuel cell cars. Our contribution consists then in finding the accurate sizing (wind turbines, solar panels as primary sources, hydrogen as a long term storage and batteries as a back up) for both data centers and fuel cell vehicles.

I. INTRODUCTION

Nowadays, everybody accesses daily to public or personalized services on Internet. Thus, we can reach our bank, health insurance accounts, plan an agenda in a very transparent way and link it to the emails, we are able to share music, online photos and to communicate directly by using social networks. This space, in which we transparently navigate, is named "the Cloud". The Cloud then offers remote access to all of these spaces from different resources such as phones, tablets or computers. Thereby, the loss of a phone no longer rhymes with the loss of the address book.

Getting the possibility to connect to our email address and to check all the rest of data means that these emails are memorized somewhere on one hard disk accessible at any time into a real data-center. Indeed, there is diverse resources that receive our requests and provide us with the content of the data that we want to access. This supposes complex network connections to route the information and deal with conflict problems. These devices which were mainly deployed at the level of companies' resources a few years ago, are now clustered together in what we call a data-center.

The vast growth of internet and the increase of the number of online services have highlighted the cost imposed by the Cloud in general and by data-centers in particular at the energy consumption level. Indeed each of our stored e-mails runs permanently on one hard disk. Moreover any request to any online service, search engine, data consultation needs a particular processing. Thus using the Internet consumes electricity. According to the Environment and Energy Control Agency, sending an e-mail with a file as an attachment requires the cost of operating a high-powered bulb for one hour. Every day, 10 billion emails are sent around the world, which corresponds on average to 50 GW h, means the electrical production of 15 nuclear power stations operating for one hour [START_REF] Tison | Internet : la pollution cachée[END_REF]. Consequently, the global electricity part dedicated to their consumption has reached unprecedented levels. According to the estimates of the year 2010, data-centers consumed roughly 1.5% of the total electricity used in the world within the half attributed to cooling [START_REF] Koomey | Growth in data center electricity use 2005 to 2010[END_REF]. This electricity consumption influences on climate change as most of the electricity is produced from fossil fuels. A 2008 study has estimated that data-centers around the world emitted 116 million tons of carbon, a little more than what Nigeria emitted in the same year [START_REF] Mankoff | Some computer science issues in creating a sustainable world[END_REF]. The cost of electricity is therefore a significant burden on data-center operators, which has led to many green energy initiatives.

One of these initiatives is to allow a data-center either to generate its own renewable energy or to use power directly from the nearest renewable power plant. Several recent initiatives are proposed in that end. For example, McGraw-Hill has recently completed a 14 M W solar field for its data-center [START_REF] Miller | Huge solar array will support nj data center[END_REF] and Apple build a 40 M W solar field for its North Carolina data-center [START_REF]Apple facilities, environmental footprint report[END_REF]. Nevertheless, many researchers investigate in the same course and try to combine renewable energy sources to supply their data-centers: the European project named "GreenDataNet" combines solar panels and batteries, has designed and validated a new smart energy management system allowing urban data-centers to radically improve their energy performance [START_REF]Green data net project[END_REF]; the "Renew it" project aims at developing a simulation tool to evaluate the energy performance of different technical solutions integrating diverse renewable energy sources as wind turbines and photovoltaics in different European climate regions [START_REF]Renewit project[END_REF].

Moreover, the DATAZERO project 1 , funded by the French National Research Agency (ANR), project of which we are engaged, aims at building an innovative adaptive data-center concept that optimizes both IT and power loads. Relying on the latest technologies in terms of renewable power sources, electrical equipment, servers and cloud management systems, it aims at drastically reducing the carbon impact of the energy consumption of tomorrow's data-centers. Thus, we argued that an option design for green data-centers is to supply them with both solar/wind energy sources and fuel cells. Batteries are also essential to assure a backup function. The major research challenge with solar and wind power is that they are intermittent sources, so they are not always available. For example, the wind turbines work following the wind speed. In case the minimal wind speed is less than 15 km/h, a turbine cannot be activated to produce electricity and in case it's more than 90 km/h, a device voluntarily stops the turbine for security reasons [START_REF]Pourquoi les éoliennes sont-elles parfois arrêtées alors que le vent souffle?[END_REF]. Then, the sizing of the storage of the data-center represents a big challenge for that kind of data-center. The second challenge consists in sizing the hydrogen so it can supply both data-centers as a back-up generator in time of need and fuel cell cars as a fuel cell station.

In fact, in 2012, the number of data-center worldwide was estimated at 509,147 [START_REF] Glanz | Power, Pollution and the Internet[END_REF] and it will continue growing until it peaks at 8.6 million in 2017 with a total area of more than 186 million square meters almost the equivalent of more than 25,100 football fields, according to the recent report published by IDC [10]. Thus, these data-centers could effectively be an hydrogen station for fuel cell cars with the appropriate sizing.

Thus in this paper, we will explain the methodology used to dimension the renewable energy sources used within the DATAZERO project, to power a data center of 500kW. In particular, we are designing an first approach to make a sizing of this hybrid system. Moreover, it explains how datacenters could be a supplier for these FCV in the future if the dimensioning is more or less over sized.

II. DATA-CENTER

A. Presentation

A data-center is a set of resources that houses all manner of IT equipment, including servers, storage subsystems, networking switches, routers and firewalls, as well as the cabling and physical racks used to organize and interconnect the IT equipment. Diverse requests of different natures are addressed to a data-center ranging from consultation of information to the triggering of treatment more or less important as the daily treatments or SaaS type calculations. That's why these datacenters has to maintain a high level of security.

B. Power distribution model

The data-center energy distribution system is the equipment responsible for supplying electrical energy to the load of the system (computer and mechanical equipment). Taking in consideration the micro-breaks or longer outages that the public grid may have, it can cause a dis-function of a data-center. Thus, it is essential to ensure adequate power supply. This equipment must therefore meet the required energy quality and security of supply. Figure 1 shows a diagram of a power infrastructure of a standard data-center.

C. Implementation of renewable energy on data-centers

Data-centers are important energy consumers, consequently Renewable Energy Sources (RES) can play an important role in this industry. Researches on RES and their applications to different fields has attracted a great interest for many researchers around the world. To give arguments of the integration of renewable energy sources into real systems, Dalton et al. conducted technical and financial analysis of the RES-only and grid/RES hybrid power supply configurations for a large hotel connected to the grid [START_REF] Dalton | Feasibility analysis of renewable energy supply options for a grid-connected large hotel[END_REF]. Recently, researchers are more interested in using RES to provide some of the overall energy consumption of data-centers. However, while renewable energies are available intermittently, the computation demands for data-centers has to be satisfied even when green energy is not available.

Other companies used fuel cells as a building's primary source of power, using the local grid as back-up power if needed. In addition, the excess heat generated by fuel cells can be used to provide cooling for the IT rooms. For example, EBAY operates 6 M W of fuel cells powering its South Jordan, Utah data-center [START_REF] Curtin | The business case for fuel cells reliability, resiliency & savings[END_REF] and Apple has installed 10 M W of fuel cells coupled with solar panels to power its Maiden, North Carolina iCloud data-center [START_REF] Association | Fuel cells and data center[END_REF]. The rub is that fuel cells are expensive, depending on local markets. They can produce electricity at a cost higher than grid power, which does not make economical sense for a data-center. In this situation, data-centers use brown energy drawn from the grid and emit carbon in the environment.

III. MODEL OF HYBRID SYSTEM COMPONENTS

The hybrid power system composed of solar panels (PV), wind turbines (WT), fuel cells (FC), batteries, electrolyzer, hydrogen tank and other energy conversion devices is a bit complicated. In order to analyze and size correctly the introduced hybrid system, it is necessary to understand the model of each component. As a result, we can obtain the energy relationships between various components.

A. Photovoltaic panel

The relation between the irradiation data and output power of the PV panels could be described as the following equation (1):

P pv = I × A pv × η pv (1)
where, I is the hourly solar irradiance in kW/m 2 , A pv is the area of the PV panel in m 2 and η pv is PV panels' efficiency. For our study, we choose the performance of an IBC PolySol 265 CS4 (265 W ) solar panel as provided in Table I.

B. Wind Turbine Generator

The model for the output power of one wind turbine generator follows the power curve shown in Figure 2. So, the turbine starts generating power at the cut-in wind speed V ci . Then, the generated output power increases with the increases of wind speed from the cut-in V ci to the rated wind speed V r . When the wind speed varies between the rated wind and the cut-out wind speed V co , which is the maximum wind speed value at which the turbine can correctly work, the turbine produces a constant or rated power. Once the wind speed goes beyond the cut-out speed, the turbine stops generating for safety reasons. The mathematical model of this behavior is as follows the equation ( 2):

P w =          0 if v(t) ≤ v ci or v(t) ≥ v co P r v 3 (t) -v 3 ci v 3 r -v 3 ci if v ci < v(t) < v r P r if v r < v(t) < v co (2)
where v(t) is the wind speed (m.s -1 ) at any time t (s). The total power output of the wind turbine generator at an hour t is given by the equation (3):

P wt = P w × A w × η w (3)
where A w is the total swept area by the blades and η w is the efficiency of the wind turbine. Different wind turbines have diverse power output and performance curves. Therefore, the modeling equation of a wind system is strongly influenced by the power curve of the wind turbine that is used. For our study, we choose the parameters given in Table II. Above all, the energy generated in one hour is given as follows:

P tot = P pv + P wt (4)

C. Hydrogen System

The hydrogen system is composed of electrolyzer, fuel cell and hydrogen tank. Taking into account the limitation of the batteries, the over production of electrical energy (when batteries are full), the electrolyzer could transform the difference to dihydrogen molecules and store them into tanks so there is no waste. This process is explained by the following chemical reaction ( 5):

H 2 O → H 2 + 1 2 O 2 (5)
In order to raise density of stored energy, a compressor may pressurize electrolyzer's output up to 200 bar [START_REF] Garcia | A wind-diesel system with hydrogen storage: Joint optimisation of design and dispatch[END_REF]. Then, the power from electrolyzer to hydrogen is expressed as follows:

P tank = P diff × η el (6)
where η el is the electrolyzer efficiency which assumed as a constant. P tank is the transferred power from the electrolyzer to the tank and P diff is the transferred power from generated power to the electrolyzer.

To realize the same reaction in the other sense, producing electrical energy from hydrogen, we can use a fuel cell. The chemical reactions at each cathode are the following:

• Anode reaction: H 2 → 2H + + 2e - • Cathode reaction: 1 2 O 2 + 2H + + 2e -→ H 2 O
Under given working conditions, 1 mole of hydrogen is equivalent to 2.10 3 kg equal to 241.8 kJ. To be more specific: 1 kg hydrogen = 120.910 3 kJ = 33.6 kW h

The converted energy can be calculated at any time t as follows:

P tank = P diff × η el × η f c (7)
where η f c is the efficiency of the fuel cell.

IV. SIZING STRATEGY

In this paper, the sizing of the hybrid system operates as the following strategy:

1) Calculation of the needed energy demand during 1 year of the data-center and for a portion of 10 hybrid cars. 2) Calculation of the energy produced from 1 wind turbine during 1 year . 3) Calculation of the energy produced from 1 photovoltaic panel during 1 year. 4) Calculation of number of WT and PV that we need. 5) Calculation of the energy difference between generated and needed. 6) Conversion of the energy difference to a volume of hydrogen. 7) Calculation of number of tank. We propose in this paper to dimension the infrastructure as an iterative process. At the beginning, we start by using the maximum of n wind turbines and completing the difference needed with photovoltaic panels. Then, we use (n -1) wind turbines and complete panels till we are able to reach the power demand with only using photovoltaic modules.

A. Calculation of the needed energy

The task of this paper is to find a sizing of system components to power a data-center of 500 kW and to be in the same time a station for hybrid cars. Thus, we supposed that the data-center needs 500 kW each hour during one year. Then the energy needed over one year is calculated in the following equation ( 8):

D1 = P datacenter × 24 × 365 = 4380000 kW h (8)
We assume that a hybrid cars consumes around 0.76 kg of hydrogen to travel 100 km [16]. If we suppose that every cars travels 31 km on average per day [17], their consumption in hydrogen in (kg) for one year is given by C in the following formula:

C = 31 × 0.76 100 × 10 × 365 = 875 kg (9)
This quantity of hydrogen corresponds to an electrical energy value 29433.6 kW h. The total energy needed is then:

D = D1 + C = 4409433.6 kW h (10)

B. Calculation of the energy produced

In calm weather, the idealized daily wind speed distribution can be represented like the typical curve of photovoltaic system. Wind speed is higher during the day (from sunrise to sunset) than during the period after sunset and sunrise. Wind speed follows a sine curve during the day, and it remains constant for periods before sunrise and after sunset [START_REF] Marinelli | Wind and photovoltaic large-scale regional models for hourly production evaluation[END_REF]. Moreover we assume that the production of energy by the primary sources during one year linearly increases when days get longer and linearly decreases when days get shorter (see Figure 3). In this figure n represents the number days (n = 365), max and min are respectively the maximum and minimum production of energy in the year for each day. So 0 or n corresponds to December the 21th (the day with the shortest light) and n/2 to June the 21th (the day with the longest light). During one year, the electrical production varies greatly over season. With making a zoom on the daily scale, it's obvious to differentiate the best production day (E) and the worst one (e) as showed in Figure 4 and where D is the energy demand of the data-center. Considering our model, the overproduction of one day has to offset the underproduction of another day. And using both the wind turbine and the photovoltaic panel models introduced before, we can estimate their electrical production in kW h during that two days. This allows us to evaluate the primary source sizing. To make this sizing, we choose the worst and the best days (0 and n/2). As mention before, the objective fixed is to fulfill the lack of energy of the worst day by the overproduction of the best day and to repeat this operation all over the year by taking into account the storage efficiency η. This problem is summarized by these equations [START_REF] Goiri | Matching renewable energy supply and demand in green datacenters[END_REF]:

Energy profile (kW

Energy (kW h)

(E -D) × η = D -e E × η -D × η = D -e ηE + e = D(1 + η) (11)
We assume the energy produced during the best day E and the energy produced during the worst day e could be written using a proportionality coefficient α: e = αE.

ηE + αE = D(1 + η) (η + α)E = D(1 + η) E = D(1 + η) (η + α) (12)
Then, the sizing of WT and PV is made following the Power demand that should be provided on the best day. Considering the value of E and the energy model for irradiance and wind, we can found the pic of power that the photovoltaic panels and the wind turbines have capture to produce E at the end of the best day. This power is P demand .

C. Sizing of WT and PV

Once the electrical production of one WT is fixed, the program proceed to the calculation of the number of wind turbines as given in the following statement [START_REF] Curtin | The business case for fuel cells reliability, resiliency & savings[END_REF]:

N wt = P demand /P wt (13)
where N wt is the number of wind turbines, P wt is the power that is observed of the WT when P demand is received by the wind conditions. Then, the area of photovoltaics N i pv is given by ( 14) where i is the number of WT are is chosen:

N i pv = P demand -N i wt × P wt P pv (14)

D. Calculation of the energy difference

Taking into account the storage efficiency, the model used in the beginning is not anymore accurate. Actually the real model then representing the system could be defined in Figure 5. where max and min are respectively the maximum and minimum energy produced in a day during the year, D is the power of the data-center and n = 365 is a year.

To conclude, the energy difference E diff consists actually of the area of the triangle above the demand and is calculated by [START_REF] Garcia | A wind-diesel system with hydrogen storage: Joint optimisation of design and dispatch[END_REF]:

E diff = D×24×(1+η) η+α × n 2 2 (15)

E. Conversion of the energy difference to a volume of H 2

To convert, in an easy way, an electrical energy to a mass or a volume of hydrogen, these steps should be followed. In a first time, use the enthalpy of H 2 in its lowest value (when dihydrogen is gas) ( 16

To convert it in kW h, more used in the electrical language:

1 kg = 120.9 3600 = 33.6 kW h

In the end, to transform it in a volume using the unit normolitre:

1 mol ↔ 2.10 -3 kg ↔ 22.4 N l 1 kg = 22.4 210 -3 = 11.210 3 N l = 11.2N m 3 (18)
To conclude, the conversion in the program is executed following this expression (19):

V = E diff × η el × η f c × 11.2 33.6 (19)
where V is the volume of dihydrogen in N m 3 , E diff is the energy difference between generated and demand and η el , η f c are respectively the efficiency of the electrolyzer and the fuel cell.

V. RESULTS

In order to perform the sizing methodology, results are performed using input data profiles. Solar radiations are obtained from PVGis, an online free solar photovoltaic energy calculator for stand alone or connected to the grid PV systems and plants, in Europe and the load from the instruction of the DATAZERO project.

When we use the calculation described before, we obtain only 5 possible cases as detailed in Table III: These results represent a first simple sizing which seems efficient for powering the data-center and a range of 10 hybrid cars per day during a year. As represented in Table III, the volume of hydrogen does not change following the change of the primary sources.

Actually the results of wind Turbines and Photovoltaics are drawn in Figure 6. In this paper, we have proposed a first sizing for the hybrid system powering the new data-center DATAZERO. The sizing was precise for the primary sources composed of photovoltaics and wind turbines. We calculated the energy needed to the data-center and a range of 10 hybrid cars by using an iterative process that calculates the number of wind turbines and PV area needed in each case.

The objective was also to size in an accurate way the storage part composed of batteries and fuel cells. In this paper, the batteries were supposed existing and our focus was on the hydrogen system. Perspectives is now to focus on sizing of fuel cells, electrolyzer needed, batteries, etc.

To conclude, by presenting an over sizing of the primary sources and hydrogen tank, we explain how data-centers all over the world can help in the development of FCV by being their gas stations.

Fig. 1 .

 1 Fig. 1. Electric System Distribution of a data-center (source: [11])

Fig. 2 .

 2 Fig. 2. Ideal wind turbine power output

Fig. 3 .

 3 Fig. 3. Primary energy produced by PV and WT

Fig. 4 .

 4 Fig. 4. Primary energy produced by PV and WT

Fig. 5 .

 5 Fig. 5. Primary energy produced by PV and WT

  ): ∆H = 241.8 kJ/mol 1 mol ↔ 2.10 -3 kg ↔ 241.8 kJ 1 kg = 241.8 210 -3 = 120.9 kJ

Fig. 6 .

 6 Fig. 6. Number of Wind Turbines in function of PV area

http://www.datazero.org

ACKNOWLEDGMENT

This work was supported in part by the ANR DATAZERO (contract "ANR-15-CE25-0012") project and by the Labex ACTION project (contract "ANR-11-LA BX-01-01").


Christophe Guyeux
email: christophe.guyeux@univ-fcomte.fr

Stéphane Chrétien
email: stephane.chretien@univ-fcomte.fr

Nathalie M Côté

Jacques M Bahi
email: jacques.bahi@univ-fcomte.fr

M.-L Côté

Systematic investigations of gene effects on both topologies and supports: an Echinococcus illustration

Keywords: Molecular phylogeny, Echinococcus, Statistical tests

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

Introduction

Echinococcus (Cestoda: Taeniidae) species are parasitic organisms responsible for echinococcosis, a serious disease with fluid-filled cycle lesion that affects humans, livestock, and wildlife. Based on morphology and life cycle data, only E. granulosus (10 genotypes), E. multilocularis, E. vogeli, E. oligarthra, and E. shiquicus have been identified as true Echinococcus species [START_REF] Thompson | Towards a taxonomic revision of the genus Echinococcus[END_REF][START_REF] Nakao | Phylogenetic systematics of the genus echinococcus (cestoda: Taeniidae)[END_REF]. But, despite extensive research, the phylogenetic relationship between species of this genus remains unclear. The large amount of genetics data, recently obtained with complete mitochondrial (mt) genomes, can provide a key to better understand the phylogeny of the Echinococcus genus.

Genes from mt genomes are standard markers for phylogeny, as this genome presents interesting features for such types of analysis. Indeed, these genes are present in almost all animals and in only one copy. The organelle is maternally inherited and the genome does not recombine in most cases [START_REF] Ballard | The incomplete natural history of mitochondria[END_REF]. Moreover, animal cells contain many mitochondria, which also contain many copies of mt genomes. This feature can be useful when only minute amount of biological material is available for molecular analysis. Finally, coding sequences within this genome are well known: twenty-two tRNAs and 2 rRNA are encoded for mitochondrial translation machinery, while in most cases there are 13 genes coding for proteins of the respiratory chain (in some species, some genes have been lost during evolution, such as atp8 within the Cestoda and Trematoda for instance [START_REF] Donald P Mcmanus | Genomics of parasitic flatworms[END_REF]). This is why phylogenetic analyses of animals are frequently performed by using a subset of mitochondrial coding sequences [START_REF] Hardman | Comparison of the phylogenetic performance of neodermatan mitochondrial protein-coding genes[END_REF], while chloroplasts are often considered for plant phylogenies [START_REF] Alkindy | Hybrid genetic algorithm and lasso test approach for inferring well supported phylogenetic trees based on subsets of chloroplastic core genes[END_REF][START_REF] Alkindy | Gene similarity-based approaches for determining core-genes of chloroplasts[END_REF][START_REF] Alkindy | Finding the core-genes of chloroplasts[END_REF]. For instance, in tetrapods and mammals, the best results have been obtained with nad5, nad4, nad2, cytb, and cox [START_REF] Zardoya | Phylogenetic performance of mitochondrial protein-coding genes in resolving relationships among vertebrates[END_REF]. The optimal size and the content of this subset of mt genes have however significant consequences: a smaller subset of sequences necessitates a lower cost for DNA extraction and sequencing, while demanding a lower computation time. Conversely, a larger set may be more representative, leading to a more accurate inference. Additionally, it is sometimes reported that, due to their own origin and history, the evolution of mt genomes may be different from the one of nucleus genomes, which may reflect more the species relationship. And various situations have appeared to be hard to solve by using only mitochondrial coding sequences. For instance, in the Echinococcus case, no combination of such mt genes has been able to produce a well supported tree, and so enlarging the set of sequences to nuclear genes has been considered as a way to estimate the phylogeny.

The problem in finding a relevant subset of sequences to infer a well supported phylogeny is due to the presence of "blurring" sequences that can increase uncertainties in some locations of the tree. At the origin of such blurring effects, explanations include homoplasy (spurious similarity due to convergence or reversion and not to common ancestry), stochastic errors, undetected paralogy, incomplete lineage sorting, horizontal gene transfers, or even hybridization. totally different phylogenetic trees. But, if we except the use of gene trees, there is a lack of tools that allow to reinforce the confidence put in a biomolecular phylogeny by understanding the impact of gene selection on a given phylogeny. In this article, our aim is to present a way to obtain robust taxonomies for phylogenetic problems which are in particular hard to solve, by using a statistical approach (logit regression with potentially dummy variables and LASSO tests) and extensive computations on available genetic data. Our proposed methodology is able to detect blurring genes and to understand their effects on topologies and on bootstrap supports. This method is presented and illustrated by addressing a case that is well known to be difficult, namely the Echinococcus genus. The methodology encompasses a de novo annotation of each available mt genome that is completed with nuclear genes, the multi-alignment of each coding sequence, and the systematic investigations of trees that can be inferred using all possible large subsets of sequences. The computation of 43, 796 trees, corresponding to all the subsets from 13 to 19 genes, leads to the possibility to understand, by using advanced statistical techniques, the effects of each gene on both topologies and supports, and thus to choose the best phylogenetic representation of the considered set of species.

State-of-the-art

To date, 9 complete mitochondrial genomes of Echinococcus have been published, they are listed with their accession numbers in Table 1. These genomes or other nuclear genes have been used recently to update the molecular phylogeny of this genus. Indeed, 4 different phylogenetic trees obtained with biomolecular analyses can be found in the literature (Figure 1). Tree A has been inferred on the full mitogenome data set [START_REF] Nakao | A molecular phylogeny of the genus echinococcus inferred from complete mitochondrial genomes[END_REF] while the three other Echinococcus trees based on nuclear genes have been proposed in [START_REF] Knapp | Phylogenetic relationships within echinococcus and taenia tapeworms (cestoda: Taeniidae): An inference from nuclear protein-coding genes[END_REF]. Trees B and C have been inferred using maximum likelihood analysis (PAUP 4.0 software [START_REF] Swofford | PAUP*: phylogenetic analysis using parsimony[END_REF]) of rpb2 (RNA polymerase II second largest subunit), pepck (phosphoenolpyruvate carboxykinase), and pold (DNA polymerase delta) genes (gDNA and exon data sets respectively, with 5008 and 4726 base pairs resp.). Finally, Tree D is the result of a Bayesian inference using BEAST software [START_REF] Drummond | Beast: Bayesian evolutionary analysis by sampling trees[END_REF] based on the exon data set detailed above.

V.mustelae has been used as the outgroup for the tree B, C, and D. A general topology with low bootstrap values of the Echinococcus phylogeny can be deduced from these trees: E. vogeli and E. oligarthra appear as sister species, and this is also the case for E. canadensis and E. ortleppi ; and for E. Granulosus and E. felidis. Additionally, all these species plus E. equinus seem to constitute a clade sister to both E. vogeli and E. oligarthra. But obviously, these four trees are not well supported, and the following questions remain unsolved: what is the correct position of E. equinus ? Are E. shiquicus and E. multilocularis sister species of E. canadensis and E. ortleppi, or of E. granulosus and E. felidis ? And what is the real relation between E. vogeli and E. oligarthra ? Only a well-supported tree can provide clear responses.

Materials and methods

Alignment and annotations of coding sequences

To answer these questions, we have performed first Bayesian and maximum likelihood analyses on (1) the whole mitogenome, (2) its 12 protein-coding genes (vogeli,(((granulosus,felidis),((equinus,(multilocularis,shiquicus)),(canadensis,ortleppi))),oligarthra)) nad2 (felidis,(granulosus,((((equinus,vogeli),(canadensis,ortleppi)),(multilocularis,shiquicus)),oligarthra))) nad3 (oligarthra,(vogeli,(multilocularis,(((felidis,granulosus),(shiquicus,equinus)),(canadensis,ortleppi))))) nad4 (vogeli,((((granulosus,felidis),(multilocularis,(shiquicus,(canadensis,ortleppi)))),equinus),oligarthra)) nad4l (shiquicus,(((oligarthra,vogeli),multilocularis),((canadensis,ortleppi),(equinus,(felidis,granulosus))))) nad5

(oligarthra,(vogeli,((((felidis,granulosus),(canadensis,ortleppi)),equinus),(shiquicus,multilocularis)))) nad6 (((equinus,(((felidis,(granulosus,oligarthra)),vogeli),(shiquicus,multilocularis))),ortleppi),canadensis) pepck ((oligarthra,vogeli),((equinus,(shiquicus,multilocularis)),(granulosus,((canadensis,ortleppi),felidis)))) pold (oligarthra,(((felidis,granulosus),(((shiquicus,multilocularis),(canadensis,ortleppi)),equinus)),vogeli)) rpb2 ((((oligarthra,vogeli),(equinus,((multilocularis,shiquicus),(felidis,granulosus)))),canadensis),ortleppi) rrnL (oligarthra,((((canadensis,ortleppi),((shiquicus,multilocularis),equinus)),vogeli),(felidis,granulosus))) rrnS (multilocularis,(oligarthra,((((canadensis,ortleppi),(felidis,granulosus)),(vogeli,equinus)),shiquicus)))

merged after alignments, and (3) on all 19 available genes taken alone (leading to 19 gene trees). Theses analyses were realized using nucleotides and translated amino acid sequences. Tools used during these first analyses were:

• Muscle [START_REF] Edgar | Muscle: a multiple sequence alignment method with reduced time and space complexity[END_REF] for complete mitogenome alignments and T-Coffee [START_REF] Notredame | T-Coffee: A novel method for fast and accurate multiple sequence alignment[END_REF] for gene alignments;

• NCBI annotations for the coding sequences in a first step, and then DOGMA [START_REF] Wyman | Automatic annotation of organellar genomes with dogma[END_REF] for further analyses;

• PhyloBayes [START_REF] Lartillot | PhyloBayes 3: a Bayesian software package for phylogenetic reconstruction and molecular dating[END_REF] for Bayesian inference, PhyML [START_REF] Guindon | New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of phyml 3.0[END_REF] and RAxML [START_REF] Stamatakis | Raxml-omp: An efficient program for phylogenetic inference on smps[END_REF] for maximum likelihood.

Contrary to the intuition, for every experiment, a problem with the support (i.e. support was less than 95% recovery) was found in at least one branch in each trees obtained. Obtained gene trees illustrate well the statistical hardness of this phylogeny: in Table 2, we can surprisingly see that 19 different topologies have been obtained when considering each coding sequence taken alone. Shared structures can be found between two given gene trees, but it is hard to deduce a consensus topology from their phylogenies, due to the large variability we obtained. Furthermore, each gene tree contains numerous problematic supports, as can be seen in the depicted trees of the supplementary material (see the appendix).

Based on these preliminary results, we are able to conclude that: (1) using coding sequences is better than using the whole mitogenome, (2) there are numerous inconsistencies in NCBI annotations (cf. the appendix), (3) the alignments obtained using T-Coffee lead to better supported trees than the ones using Muscle, (4) to separate the natural history of the genus from the individual history of each gene seems very difficult, when considering the large variability of gene trees, and (5) to increase the amount of genetic data leads to better supported trees. Having these results in mind, we our proposal consists in a general methodology for estimating a phylogeny by measuring gene effects. Application to the Echinococcus case will demonstrate the potential of our approach.

Methodological approach

To estimate the phylogeny of Echinococcus and to determine which genes do not support this phylogeny, a simple solution consists in considering all the available coding sequences (mt and nucleus ones) shared by these 10 species, and to investigate how the phylogeny is affected when using subsets of these sequences. Doing so will extend the first investigations of Hardman et al. [START_REF] Hardman | Comparison of the phylogenetic performance of neodermatan mitochondrial protein-coding genes[END_REF], who studied the phylogeny of 5 Taeniidae based on the 12 individual mitochondrial coding genes.

Nineteen sequences have been extracted from each of the considered species: 5 nuclear genes, and 12 protein-coding sequences and 2 rDNA genes from the mitochondrial genomes. They are listed below. • Nuclear genes: ezrin-radixin-moesinlike protein (elp), elongation factor 1 alpha (ef1a), RNA polymerase II second largest subunit (rpb2 ), phosphoenolpyruvate carboxykinase (pepck ), DNA polymerase delta (pold ).

• Mitochondrial protein coding sequences: ATP synthase 6 (atp6 ), cytochrome b (cob), cytochrome c oxydase 1 (cox1 ), cox2, cox3, NADH dehydrogenase subunit 1 (nad1 ), nad2, nad3, nad4, nad4l, nad5, nad6 .

• Mitochondrial rDNAs: the large subunit rDNA (rrnL), and the small subunit rDNA (rrnS ).

DOGMA has been used to annotate from scratch each up-to-date complete mitochondrial genome downloaded from NCBI [6] (accession numbers provided in Table 1). Default parameters of DOGMA have been selected, with identity cutoffs for amino acids equal to 60% and 80% for coding genes and rDNA genes. These thresholds have been reduced to 55% and 75% respectively for V.mustelae, due to a problem for the detection of nad6 and rrnL. The e-value was equal to 1e -5, and the number of blast hits to return has been set to 5, see Table 8 of the appendix section for a comparison with annotated genomes from NCBI.

Each of these 19 coding sequences has been aligned separately by using T-Coffee (M-Coffee mode, using 6 cores for multiprocessing). Then 43,796 trees have been constructed, corresponding to all the possible combinations of [START_REF] Lartillot | PhyloBayes 3: a Bayesian software package for phylogenetic reconstruction and molecular dating[END_REF] bust phylogenetic trees and the effects of each gene on topologies and supports. RAxML version 8.0.20 was used for maximum likelihood inference, with 3 distinct models/data partitions with joint branch length optimization at each computation, corresponding to the nuclear genes, the mitochondrial rDNA genes, and the mitochondrial protein-coding sequences [START_REF] Guyeux | The study of unfoldable self-avoiding walks. application to protein structure prediction software[END_REF]. All free model parameters have been estimated by RAxML for the GAMMA model of the rate heterogeneity and the ML estimate of the alpha-parameter. At each time, a maximum of 1000 non-parametric bootstrap inferences were executed, with MRE-based boots stopping criterion, and V.mustelae has been used as outgroup.

for each gene g do multi-align all the 10 sequences of g using T-coffee; end for k=13,..., [START_REF] Nakao | Phylogenetic systematics of the genus echinococcus (cestoda: Taeniidae)[END_REF]

Discussion and results

Results

Only 15 topologies have been obtained during our computations, among the 34,459,425 possible phylogenetic tree topologies with 9 species and 1 outgroup. Only 4 of these 15 topologies have a number of occurrences larger than 300, when considering the 43,796 obtained trees. These trees are depicted in Figure 2, while further information is provided in Table 3: for each topology, the lowest bootstrap of the best tree (that is, the lowest bootstrap of the tree that maximizes the minimum taken over all its bootstraps), the number of trees having this topology, the average minimal bootstrap value, and the list of genes that have been removed to obtain the best tree having this topology. We can remark that Topology 0 corresponds to Tree B of the literature (Figure 1), with more robust bootstraps. Topology 2 is similar to Tree A, while the Trees C and D have not been obtained during our computations. The ordinate provides the number of trees in the topology whose lowest bootstrap is larger than the associated abscissa, for 5 removed genes. B. As in (A), but for 6 discarded genes.

The 4 best topologies represent 98.19% of the whole obtained trees, and they share most of their structure. And, at this stage, a first question raised in Section 1 can be answered: E. shiquicus, E. multilocularis, E. canadensis, and E. ortleppi are within the same clade, which is sister to the clade consisting of E. granulosus and E. felidis. For the correct position of E. equinus, Topologies 0 and 1 are opposed to Topologies 2 and 3, while Topologies 0 and 2 are in contradiction with the two other topologies for the relation between E. vogeli and E. oligarthra.

Several reasons lead us to consider the Topology 0 depicted in Figure 2(a) as the most probable one. Firstly, this is the most frequent topology, representing 53.69% of the trees, while the second one (Topology 2) represents only 24.40%. This topology remains the most frequent even when 1, 2, 4, 5, or 6 genes were removed. Surprisingly, Topology 3 is a bit more frequent than Topology 0 when 3 genes are discarded, see Figure 3(a). Secondly, this topology contains the most robust trees: the smallest bootstraps in Topology 0 are largest compared to the bootstraps of the other topologies. This property is preserved when 1 to 6 genes are removed from the dataset, see Figure 3(c). In particular, the best tree for Topology 0, depicted in Figure 2(a), outperforms the best tree of the other topologies, even when considering the smallest bootstrap or its average. Thirdly, the number of trees with no bootstrap under 80 is equal to 1,135 for Topology 0, and is equal to 29, 54, and 49 for Topologies 1, 2, and 3 respectively. This is illustrated in Figure 3(b) and this property is further investigated in Figure 4, which shows the number of trees with no bootstrap under a given threshold for each topology.

Let us remark now that we analyzed 14 mitochondrial, but only 5 nuclear genes: available data being what they are, our approach may preferably support trees based on mitochondrial genes. Additionally, the possibility of a systematic bias has been reported in General Time Reversible (GTR) maximum likelihood analyses of mitochondrial data. Trees based on specific nuclear genes can be different from the overall species-level phylogenetic tree, which can be due to incomplete lineage sorting. Thus, we have to deal with potential artifacts due to the fact that mitochondrial trees are more represented than nuclear ones. In the next sections, the relationship between genes and tree topologies/supports has been investigated in detail, in order to find genes that behave differently from the other ones (for biological or methodological reasons), and to understand their differences regarding the phylogeny.

The frequency of occurrences of the genes

Using a simple Python script, the evenness of each gene in the 4 most frequent topologies has been evaluated. Since all the 43,796 trees have been constructed using a subset of the 19 available sequences, it is relevant to consider the frequency of occurrence of each of these sequences. Table 4 Some genes, either over or under represented, seem to be associated to particular topologies. More precisely, the following points may be outlined.

1. The 3 least frequent genes of Topology 0 are 3 nuclear genes, while all mitochondrial genes are well ranked in that topology. Conversely, several mitochondrial genes are wrongly placed on the other topologies, while the first positions in Topologies 1, 2, and 3 are taken by nuclear genes. Based on mitochondrial data, Topology 0 seems to provide a better evaluation of the evolution, while pepck, rpb2, and ef1a modify this topology, leading to Topologies 1, 2, and 3. 2. The nuclear gene pepck moves E. equinus in the tree, as it is ranked 19 for Topologies 0 and 1, while it is respectively ranked 1 and 2 for Topologies 2 and 3. However, these claims need to be further investigated by a more rigorous statistical approach, which is the aim of the next sections.

Influence of genes on topology using Dummy logit model

To investigate more soundly the effects of each coding sequence on the species topology, 4 dummy binary choice logit models have been set up for each best topology using scikit-learn [START_REF] Pedregosa | Scikitlearn: Machine learning in Python[END_REF] module of Python language. The reference to the exogenous design is a 19 × 43, 796 array, each row being a vector of 0's and 1's: a 0 in position i of row k means that, in the k-th tree computation, gene number i (in alphabetic order) were discarded, and conversely the gene was conserved if the coefficient is 1. The rows are thus the "observations" while the columns correspond to regressors. The 1-d endogenous response variable is a vector of size 43, 796, with a 1 in position k if and only if the Topology 0 has been produced with the choice of genes corresponding to the row number k in the exogenous design (resp. Topology 1, 2 or 3 in the three other binary choice logit models). The model has then been fitted using the maximum likelihood approach and a Newton-Raphson solver for performing the maximization. Convergence was obtained after 8 iterations, and the results of the Logit regression are summarized in Table 5. The first conclusion that can be drawn from the investigation of the impact of each gene is that, except for the particular case of nad2 and rrnS (see column P > |z|), all the coding sequences bring phylogenetic information. Additionally, all coding sequences reinforcing the choice of Topology 0 as the best topology are mitochondrial genes (see the largest values in coef column). Conversely, nuclear genes tend to participate in the rejection of this tree (lowest values). This trend is minimal for the genes ef1a, elp, and pold, but very important for genes pepck and rpb2. These last 2 genes have been used in [START_REF] Knapp | Phylogenetic relationships within echinococcus and taenia tapeworms (cestoda: Taeniidae): An inference from nuclear protein-coding genes[END_REF] and our results could explain the limited robustness of the Topology 0 in Tree B, and how the position of E. equinus has changed, leading to the Tree C.

A

Another conclusion is that the negative effect of mitochondrial coding sequences on Topology 0 has a low impact. Cob, cox1, nad4, nad5, and nad6 have a strong impact on this topology. All these findings are consistent with the frequency of occurrences of each gene in the choice of Topology 0: nad5 has been included in 22,335 computations leading to this topology (94.99%), while only 14, 657 computations with pepck has led to this topology (62.33%), as described in Table 4.

Further investigations of the role of each sequence show that rpb2 and ef1a are responsible for the inaccurate phylogenetic position of E. vogeli in the Topologies 1 and 3 (this species is considered as the sister species of E. oligarthra in these topologies). Their logit regression coefficients are among the lowest ones in Topologies 0 and 2, and among the largest ones in the two other topologies (see the Tables 9, 10, and 11 of supplementary data, which contain the results of the dummy logit regression test for Topologies 1, 2, and 3). Similarly, the inaccurate phylogenetic position of E. equinus in the Topologies 2 and 3 is mainly due to the gene pepck, whose coefficient is the lowest one in the Topologies 0 and 1, and almost the largest one in the Topologies 2 and 3. These points contribute to explain the results obtained by Knapp et al. in [12].

Regressing the bootstrap on the genes using the LASSO method

To better understand the effects of each gene on each topology, the 4 bootstrap values depicted in Figure 5 in the Topologies 0 and 2, and Bootstrap B corresponds to the basal position of E. oligarthra, considered as a sister species of all the other Echinococcus in Topologies 1 and 3.

With and without each gene, the average bootstrap in A (resp. B, C, and D) of all the trees that present this particular topology has been computed. To have an overall and global view of the effect of each gene on the bootstrap values, average computations have been preferred to minimum ones. Obtained results are presented in the Table 6. In the absence of gene rpb2, an average bootstrap of 75.98 (dropping to 55.93 when including rpb2 ) has been computed, pointing out the influence of genes rpb2 and ef1a on the bootstraps A and B. Bootstraps C and D, for their part, are mainly affected by gene pepck.

To have a more rigorous validation of this type of influence and to determine which genes are associated with bootstraps of the different topologies, a Linear Model trained with L1 prior as regularizer (LASSO) test has been performed 4 times. By doing so, a gene that summarizes the evolution of the considered bootstrap, either by increasing it or by reducing it, has been identified. Then, by reducing the constant that multiplies the L1 term in the LASSO test, supplementary genes appear one by one, ordered by their ability to summarize the bootstrap value.

According to the bootstraps values C and D shown in Figure 5, two additional LASSO tests have been performed to determine the genes that are related to the position of E. equinus in the phylogenetic tree. They have been performed using linear_model class of sklearn [START_REF] Pedregosa | Scikitlearn: Machine learning in Python[END_REF] Python module, with inputted data being a Boolean array with 19 columns (one per gene presence) and a vector of bootstrap. The number of rows in the array corresponds to the number of trees having the topologies with this bootstrap (Topologies 0 and 2 for bootstrap A, and so on). The associated row in the vector provides the bootstrap in the tree obtained with this selection of genes. Obtained results for each bootstrap are provided in Table 7. Using these results in combination with the previous ones, supplementary arguments are provided, that show that pepck is the gene breaking sister relationship between E. oligarthra and E. vogeli, while rpb2 and ef1a tend to move the position of E. equinus. To some extent, it is possible to conclude that pepck has a negative impact on the position of E. equinus while rpb2 has a role in the erroneous position of E. vogeli in Topologies 1 and 3.

Finally, the most likely phylogenetic tree for Echinococcus genus, resulting from the statistical analyses presented in previous sections, appears to be the Topology 0. As stated before, it is the most frequent tree obtained during our investigations of the influence of genes on the phylogeny of Echinococcus, and it is also the most robust tree.

Conclusion

Deep investigation of the molecular phylogeny of the genus Echinococcus has been performed in this paper. 19 coding sequences, taken from both mitochondrial and nuclear genomes, have been considered for maximum likelihood phylogenetic reconstruction. As the trees that we obtained were not as robust that we could expect, combinations from 13 to 19 genes have been further investigated. This analysis has generated 43, 796 trees that represent 15 topologies. The position of E. equinus and the sister relationship between E. oligarthra and E. vogeli led us to the selection of 4 specific topologies. Using the logit model and the LASSO, Topology 0 has appeared as the most probable one. Finally, negative effects of several genes for that particular topology have been emphasized.

In future work, the phylogeny of the genus Taenia [START_REF] Al-Nayyef | Taenia biomolecular phylogeny and the impact of mitochondrial genes on this latter[END_REF] will be explored using a similar approach, as some species seem to have discrepant phylogenetic positions. All the possible combinations of the 14 mitochondrial genes (12 proteincoding genes and 2 ribosomal RNAs) will be considered, leading to the production of 16, 384 phylogenetic trees. Their topologies will be compared, and the influence of each gene on these topologies will be rigorously measured, in order to determine the most probable phylogenetic tree of this species. Finally, the phylogeny of the class Eucestoda will be investigated using a similar approach.

All computations have been performed using the Mésocentre de Calcul de Franche-Comté facilities. Authors would like to thank Arnaud Mouly for helpful discussions, and the staffs of Asahikawa Medical University (Japan) and of other institutes that have kindly submitted their sequences to the GenBank database.

Appendices

Species

  q u i n u s E . o r t l e p p i E . m u l t i l o c u l a r i s E . o l i g a r t h r u s E . g r a n u l o s u

6 0EFigure 2 :

 62 Figure 2: The 4 trees with more than 300 occurrences, when considering 43,796 trees obtained with Algorithm 1 (V. mustelae is an outgroup)

1 :Figure 3 :

 13 Figure 3: Comparison of the 4 best topologies, according to the number of discarded genes. A. The number of trees in each topology, according to the number of discarded genes. B. As in (A), but considering only trees whose supports are larger than 80. C. Minimal support in the best tree of each topology regarding the number of discarded genes. D. Averages of all minimal bootstraps in each tree of each topology.

Figure 4 :

 4 Figure 4: Comparison of the best topologies, according to the number of supported trees. A.The ordinate provides the number of trees in the topology whose lowest bootstrap is larger than the associated abscissa, for 5 removed genes. B. As in (A), but for 6 discarded genes.

Figure 6 : 22 (

 622 Figure 6: Gene trees of the 19 considered genes

Figure 7 : 23 (Figure 8 : 24 (Figure 9 : 25 (Figure 10 :Figure 11 : 27 (Figure 12 :Figure 13 :Figure 14 :Figure 15 :

 72382492510112712131415 Figure 7: Gene trees of the 19 considered genes (cont.)

  Due to such blurring genes, two different subset of sequences may lead to two

	Species	Accession	reference
	Echinococcus canadensis	NC_011121 [17]
	Echinococcus equinus	NC_020374 [20]
	Echinococcus felidis	NC_021144 [16]
	Echinococcus granulosus	NC_008075 [14]
	Echinococcus multilocularis NC_000928 [18]
	Echinococcus oligarthra	NC_009461 [17]
	Echinococcus ortleppi	NC_011122 [17]
	Echinococcus shiquicus	NC_009460 [17]
	Echinococcus vogeli	NC_009462 [17]
	Vesteria mustelae	NC_021143 [16]

Table 1 :

 1 Echinococcus species analyzed in this paper and the accession numbers of mitochondrial genomes. V. mustelae is an outgroup.

Table 2 :

 2 Gene tree topologies

	gene topology (newick)	(vogeli,((felidis,granulosus),(equinus,(multilocularis,((canadensis,ortleppi),(shiquicus,oligarthra)))))) atp6	(multilocularis,((equinus,((felidis,granulosus),(oligarthra,vogeli))),(shiquicus,(canadensis,ortleppi)))) cob	(oligarthra,(vogeli,((((granulosus,felidis),(canadensis,ortleppi)),(shiquicus,multilocularis)),equinus))) cox1	(((((vogeli,(canadensis,ortleppi)),(shiquicus,(multilocularis,oligarthra))),equinus),felidis),granulosus) cox2	(vogeli,(((granulosus,(oligarthra,felidis)),(canadensis,ortleppi)),((shiquicus,multilocularis),equinus))) cox3	(felidis,((canadensis,ortleppi),(((multilocularis,shiquicus),((oligarthra,vogeli),granulosus)),equinus))) ef1a	(oligarthra,(((equinus,(shiquicus,multilocularis)),((canadensis,ortleppi),(felidis,granulosus))),vogeli)) elp	nad1

Table 3 :

 3 Detail of obtained topologies. Lowest bootstraps, number of occurrences, and not included genes are provided.

  summarizes obtained results.

	Topologies	0		1		2		3
		number	rank	number	rank	number	rank	number	rank
	atp6	16805	7	3684	9	7586	14	2932	13
	cob	17074	5	4038	4	7467	15	2881	15
	cox1	20409	2	2884	18	7601	13	1459	19
	cox2	16926	6	3568	11	7609	12	2890	14
	cox3	15310	16	3390	14	9789	2	3454	4
	ef1a	15117	17	4635	2	7765	8	3594	3
	elp	16097	13	3403	13	8695	5	3006	10
	nad1	15622	15	3377	15	8886	4	3442	5
	nad2	16500	11	3756	8	7742	9	3347	6
	nad3	16758	8	3551	12	8654	6	2948	12
	nad4	17222	3	3795	5	6848	18	2813	17
	nad4L	16676	9	3760	7	7694	11	3054	9
	nad5	22335	1	3203	16	7133	16	2776	18
	nad6	17160	4	4299	3	6897	17	2844	16
	pepck	14657	19	2832	19	21336	1	3612	2
	pold	16322	12	3675	10	7698	10	2959	11
	rpb2	15052	18	4803	1	6795	19	12057	1
	rrnL	15890	14	3033	17	8922	3	3133	8
	rrnS	16590	10	3787	6	7849	7	3210	7

Table 4 :

 4 Number of times each gene were present to produce a tree having one of the 4 most relevant topologies.

Table 5 :

 5 Dummy logit regression results for Topology 0

3. Similarly, rpb2 changes the relationship between E. vogeli and E. oligarthra.

  have been investigated with more detail . Bootstrap A corresponds to the sister relationship of E. oligarthra and E. vogeli

					Bootstraps
		A		B		C		D
		without with without with without with without with
	atp6	73.29	73.29	62.07	60.19	81.92	81.61	76.48	74.91
	cob	72.14	75.56	61.22	62.19	81.42	83.64	75.90	77.27
	cox1	72.05	75.46	61.67	60.90	81.35	83.09	76.84	70.18
	cox2	72.87	74.18	61.83	60.77	81.26	83.39	76.20	75.87
	cox3	73.75	72.03	62.31	58.96	81.34	82.70	76.20	75.97
	ef1a	73.31	73.23	57.43	67.35	82.25	80.83	73.09	82.91
	elp	73.57	72.61	61.92	60.36	81.74	82.03	76.27	75.75
	nad1	73.75	72.06	62.38	58.61	82.17	81.24	76.87	74.37
	nad2	73.03	73.93	61.30	62.25	81.99	81.47	76.13	76.14
	nad3	73.11	73.71	61.94	60.38	81.46	82.73	76.45	75.15
	nad4	71.86	75.88	60.22	64.29	81.69	82.59	75.96	77.16
	nad4l	73.24	73.41	61.71	61.27	81.78	81.99	76.02	76.44
	nad5	72.51	74.55	61.73	61.06	81.25	84.29	76.67	71.95
	nad6	71.82	76.01	60.21	63.55	81.84	81.84	75.92	77.41
	pepck	73.53	68.85	61.78	56.12	78.32	82.83	72.83	78.03
	pold	74.29	70.80	62.82	57.93	85.04	73.49	78.24	69.81
	rpb2	75.98	55.93	50.06	61.68	83.07	68.26	66.75	77.18
	rrnl	74.24	70.84	62.72	55.48	83.21	79.36	78.22	70.53
	rrns	73.03	73.93	61.35	62.12	80.53	85.09	74.77	79.58
	total	73.29		61.58		81.84		76.13

Table 6 :

 6 Evolution of average bootstraps A, B, C, and D when considering or not each gene.

Table 7 :

 7 Genes with the largest impact regarding bootstraps A, B, C, and D (important ones are listed first).

Table 8 :

 8 Comparaison between DOGMA and NCBI annotations

		Gene	DOGMA location NCBI location Differences (bp)
	E. canadensis	NAD4	4445-5650	4394-5650	51
	E. equinus	NAD4	4339-5547	4291-5547	48
	E. equinus	COX2	12475-13041	12475-13053	12
	E. felidis	NAD4	4355-5557	4301-5557	54
	E. felidis	COX2	12497-13069	12497-13075	6
	E. granulosus	NAD4	2070-3275	2019-3275	51
	E. granulosus	COX2	10182-10748	10182-10760	12
	E. granulosus	NAD5	11763-13328	11763-13331	3
	E. multilocularis COX2	12602-13168	12602-13180	12
	E. oligarthra	NAD4L 4216-4473	4201-4473	15
	E. oligarthra	NAD4	4488-5696	4437-5696	51
	E. oligarthra	NAD2	6419-7294	6413-7294	6
	E. oligarthra	COX2	12659-13225	12659-13237	12
	E. ortleppi	NAD4	4442-5647	4391-5647	51
	E. ortleppi	COX2	12588-13154	12588-13163	9
	E. shiquicus	NAD4L 4218-4472	4212-4475	9
	E. shiquicus	NAD4	4490-5695	4439-5695	51
	E. shiquicus	COX2	12664-13236	12664-13242	6
	E. vogeli	NAD4	4470-5675	4419-5675	51
	E. vogeli	COX1	9189-10793	9189-10782	12
	E. vogeli	COX2	12623-13189	12623-13204	15
	V. mustelae	NAD5	434-1990	434-1999	9
	V. mustelae	NAD4L 4074-4325	4074-4331	6
	V. mustelae	NAD4	4340-5551	4295-5551	45
	V. mustelae	ATP6	5766-6269	5757-6269	9
	V. mustelae	NAD2	6298-7149	6280-7152	21
	V. mustelae	NAD3	8588-8905	8561-8905	18
	V. mustelae	COX1	9041-10654	9041-10657	3
	V. mustelae	COX2	12466-13038	12466-13041	3