Philippe Marchina

Alain Brousse

Jean Fontaine

Christophe Dano

Claudine Alonso

Claudine Alonso In

In situ measurement of rock compressibility in a heavy oil reservoir

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

The accurate determination of depletion induced compaction is very important in weakly-to non-cemented oil reservoirs such as encountered in the heavy oil fields, as it controls the surface subsidence and, more essential, the production profile and ultimate recovery from the field. Yet, in this class of reservoirs, the coring process generates such a mechanical damage in the recovered rock samples that the representativity of laboratory derived rock mechanical parameters required for the computation of reservoir compaction can be questioned. Therefore, to palliate that state of facts, the feasibility of measuring in-situ the compressibility of the rock has been investigated. In-situ measurement of rock compressibility with a dilatometer has been identified as a viable approach. Although widely used in geotechnical engineering to measure soil and rock in-situ stiffness, the deployment of a dilatometer in an oil well and the interpretation of the results in terms of rock compressibility had never been attempted so far. Both the operational aspects of the test and the interpretation of the results are presented and discussed in this paper.

Generalities

A vast majority of heavy oil reserves are found in shallow, weakly-to non-cemented sand reservoirs. As a result, the cold production of heavy oils may be controlled to a large extent by a compaction drive mechanism which can greatly affect the total recovery. The amount of compaction will mainly depend on the mechanical properties of the reservoir rock and in particular its compressibility, both elastic and plastic. Since the plastic compressibility is many times higher than the elastic one, it is also crucial to be able to assess whether the reservoir will reach a plastic regime during its production life. Owing to the specificity of their behavior the determination of the mechanical properties of weakly cemented sands from laboratory testing are made difficult. The main reason for this is in relation with the poor degree of representativity of the recovered samples due to the severe disturbances generated by the coring operations. As a consequence, the feasibility of measuring the formation compressibility in-situ has been investigated. Two options could be envisaged. The first one would be to equip a well with radioactive markers, shot into the formation at precisely know depth intervals. By logging, with the appropriate tool and at various depletion levels, the changes in the relative positions of those markers, it is possible to obtain a reservoir compaction profile. Although that technique is field-proven [START_REF] Menghini | Compaction Monitoring in the Ekofisk Area Chalk Fields[END_REF] , it is not a predictive method in the sense that it does only yield a posteriori information on the compaction. In the present case, it was crucial to be able to anticipate the reservoir behaviour with depletion, and therefore the option of the radioactive markers was not further considered. Instead, the use of a dilatometer, was identified as being an adequate technique to obtain the required up-front information.

In-situ Measurement of Rock Compressibility in a Heavy Oil Reservoir Philippe Marchina, SPE (Total), Alain Brousse & Jean Fontaine (Me2i), Christophe Dano (Ecole Centrale de Nantes) and Claudine Alonso (SINCOR) However a number of issues had to be addressed, from the theoretical feasibility of a compressibility measurement with a dilatometer through to the deployment of a non-standard oilfield tool in an oil well.

The dilatometer

The dilatometer is a tool designed for measuring the radial deformation of a borehole in response to a pressure increase inside that borehole. It consists of an inflatable membrane fitted with 3 displacement sensors capable of measuring the radial displacement of the membrane and a pressure gauge reading the fluid pressure inside the tool (Fig. 1). The tool is deployed by lowering it in the borehole to be tested and the formation is submitted to various cycles of inflation-deflation during which the pressure and radial deformations are recorded. The curves generated during a test (Fig. 2) are then inverted to derive the sought formation parameters. Essentially used in geotechnical engineering, the prime use of the dilatometer is to derive the elastic modulus of the tested rock mass.

Radial displacement (mm)

Pressure (bar)

Figure 2 Typical dilatometer test output

Feasibility study

Since, in this particular situation, the information sought deals with the plastic range of rock deformation, which goes beyond the usual field of application of the dilatometer, to validate the approach a theoretical feasibility study was undertaken. This has been achieved by means of numerical simulations. In that exercise, a dilatometer test was modeled numerically, using a Finite Element (FE) code, and the impact of the rock properties on the response (in terms of pressure vs. lateral displacement) was investigated.

Model geometry and boundary conditions

The model was set up to simulate the behavior of a dilatometer (1m long, 8 ½" diameter) in a vertical borehole of the same diameter as the tool. In view of the various symmetries of the problem (vertical axis and horizontal plane through the middle of the dilatometer), only the top half of a vertical section has been considered. The boundary conditions, as shown on the diagram (Fig. 3), assumed constant vertical and horizontal stresses and no vertical displacement across the bottom boundary (plane of symmetry for the problem). To simulate the whole operational process, the loading phase was split in two stages: i)

First, an unloading of the horizontal stress, simulating the drilling phase, was applied along the entirety of the borehole wall (left vertical boundary of the model). This unloading corresponds to the difference between the virgin horizontal in-situ stress and the mud pressure applied during the drilling phase. ii)

Secondly, the dilatometer test itself was simulated by an increasing load (stress) applied on the borehole wall over the half-length of the dilatometer (vertical symmetry). Perfectly drained conditions were assumed during that second stage. The various simulations carried out show that the results of a dilatometer test are significantly affected by the compressibility properties of the rock being tested. In the present case, within the framework of a Cam-clay formulation, both the consolidation pressure (P'c 0 -see Fig. 4) and the plastic compressibility (λ) have a definite influence on the simulated radial displacements of the dilatometer membrane as a function of the applied pressure. In addition, when looking at the results of the simulations in the stress spacenot presented here-, it can be seen that the plastification of the material induces a strong increase in the mean stress which is associated with a strongly contractant behavior. These two observations are a clear sign that the use of a dilatometer to evaluate the compressibility of the rock is feasible and thus validate the approach taken. This statement however does not imply, at this stage, that the sought information (plastic compressibility and consolidation pressure, for example) can readily be obtained by mathematical inversion of the results.

Inversion procedure

Having established the theoretical feasibility of the measurement with a dilatometer, the issue of parameter identification from a dilatometer test was addressed. Existing work on in-situ compressibility measurements using a dilatometer, mainly focussed on soft, low permeability sediments [START_REF] Rangeard | Modified Cam-Clay Parameters Identification From Pressuremeter Test[END_REF] could be adapted to typical high permeability, weakly cemented heavy oil reservoir formations. The numerical method developed there to inverse dilatometer tests consists in minimising the difference between the outcome of a simulation and the actual test data. To achieve this, an iterative direct (or forward) modelling scheme is used in combination with a parameter optimisation routine. The flow chart describing the process is shown in Fig. 5. In addition to the above, a preliminary analysis is usually done to fine tune the initial set of parameters used. This warrants that the inversion scheme converges towards a physically acceptable solution. This is achieved by a rough sensitivity analysis until a set of data is obtained that gives a resulting simulation 'close enough' to the experimental data. The approach described above having been primarily developed to analyse dilatometer tests in low permeability (typically argilaceous), shallow sediments, it remained to be seen whether it was applicable to heavy oil reservoirs. For this, the inversion of a synthetic case representing a highly permeable formation saturated with a viscous fluid, was set up with a particular attention paid to the number and the nature of the parameters that could be simultaneously identified from a single dilatometer curve.

The exercise was carried out by generating the results of a synthetic test from a base set of input parameters chosen so as to represent a dilatometer test in a typical heavy oil well.

Then, a perturbation introduced on selected parameters would generate a new set of input data, which was used as initial data set for the inversion scheme. It then could be checked whether the 'true' parameter values could be found again, thus the validity of the inversion. At the same time the sensitivity of non-perturbed parameters on the obtained value was also evaluated.

From the various attempts, it could be demonstrated that the simultaneous inversion for up to two parameters (e.g. consolidation pressure and plastic compressibility) could be achieved with satisfactory accuracy from one single dilatometer curve.

Optimization of the operational procedure

An important corollary of the above work is that it allows optimizing the operational procedure to be followed so as to maximize the probability of obtaining a set of experimental results that can be conclusively inverted. The general dilatometer test sequence is a continuous pressurization of the tool with a series of two or three unloading-reloading cycles at various pressure levels in the test. The sole objective of these cycles is to allow the determination of the elastic parameters of the formation tested. Previous work [START_REF] Rangeard | Permeability effect on pressuremeter test results[END_REF] had shown that the response of the dilatometer is strongly affected by the drainage conditions during the test, and that, unless the pore pressure could be measured, the inversion of the results could only be achieved for situations close enough to either perfectly drained or else undrained conditions. Thus, for the present case, optimizing the sequence amounted to determining the optimum pressurization rate to be applied during the dilatometer test so as to ensure that the tests are done under conditions close enough to perfectly drained. For a given pressurization rate, the pore pressure response during a test will be governed, among other parameters, by the permeability of the formation and the viscosity of the pore fluid. In the present case, using the formation permeability value and the in-situ oil viscosity yielded a 10 -7 ms -1 intrinsic permeability. To account for some degree of variability around those values, it was assume that the actual intrinsic permeability value would be in the 10 -7 and 10 -6 ms -1 range.

From the results it appears that, under the assumed prevailing conditions, a loading rate of 1 bar/min, corresponds, for all practical purposes. At the other end of the scale, and for the same permeability conditions, an undrained test would require a loading rate of 100 000 kPa/min (1000 bar/min!). This rate being unrealistic in practice, the option of performing an undrained test was not further considered.

The operational program

The operational program consisted of three phases: i) the drilling of the well down to the top of a potential test interval ii) the coring of the test interval, so as to ensure maximum borehole wall quality, and iii) the performing of the test itself.

Drilling program

The proposed drilling program included a two-stage operation.

First, an S-shape pilot hole would be drilled so as to reach the expected top sands vertically but with a limited offset (less than 50 ft) to the top of the hole. The objective of that pilot was to identify the potential sand intervals to be tested; the data acquisition program would include basic MWD (Gamma ray & resistivity) to identify the sand intervals and a wireline caliper. The reason for running that additional log was to evaluate the quality of the hole, considering that it had been observed in several instances that some zones can show a natural tendency to caving. Of course such zones, if present, should be avoided in the coring process. Following the logging acquisition, the pilot hole would then be plugged and abandoned. The actual test well would be a vertical sidetrack of the pilot hole, drilled out of the casing shoe down to the top of the first sand to be tested.

Coring program

Subsequently, a 20ft section of sand would be cored, rather than drilled, to maximize the chances that the resulting borehole is of perfect quality. The diameter of the hole in the test interval has to allow sufficient clearance between the tool and the borehole wall so that it is possible to run the tool, but that clearance has to be minimized so as to maximize the 'useful' range of the displacement sensors. The basis for design was a nominal 4" hole (101.6 mm), which would leave a 5.8 mm clearance with a 90 mm diameter dilatometer. Nevertheless, it was decided to design a 3"7/8 OD coring tool (98.425 mm) in anticipation of some degree of wash-out. With no tool available with that diameter, the entire coring assembly had to be designed and fabricated.

Dilatometer testing program

Once the test hole cored and logged with a caliper, the dilatometer would be rigged up and lowered in place on drillpipe. The testing itself would consist in pressurizing the dilatometer at the required rate of 1 bar/min until either the maximum pressure available or the maximum allowable tool deformation has been reached. A number of unloadingreloading cycles would be performed at various pressure stages during the test to evaluate the elastic properties of the tested interval.

The operations

The drilling and logging of the pilot hole allowed to identify potential test intervals (Fig. 6). The test well was drilled down to the top of the identified interval which was cored with the specially designed coring tool. The coring operation went smoothly and, despite the inability to recover any core material, the hole was perfectly smooth (Fig. 7). The dilatometer was made-up and the whole system (tool + pressure system + data acquisition) successfully tested at surface by inserting the tool in a steel pipe of adequate diameter (a drill pipe rack in this instance) and applying a slight pressure to check the proper functioning of the equipment. The dilatometer was then rigged up and connected to the drill pipe. The tool was run in at moderate speed to allow two to three straps per stand of drill pipe. The experimental set-up did not allow a continuous monitoring of the electric communication with the tool, but checks were done at regular intervals during the running-in. The running speed was lowered to a minimum (estimated at 3-4 ft/min) to enter the test hole. Entry into the small diameter hole was confirmed by a very little slack-off of the hook load, thus confirming the correct depth correlation. The actual test went smoothly as per the program. The raw test data are shown in Fig. 8. The results show that the readings given by the three sensors are slightly different, which could be attributed to a marginal formation anisotropy. Also, it can be noticed during the first unloading-reloading cycle, performed at a similar pressure rate as the test i.e. 1 bar/min, that the recorded deformation continues even with a decreasing tool pressure. This feature indicates some degree of time dependent behavior of the reservoir. Subsequent cycles were done at a faster pressurization rate and do not show that feature.

Depth (ft)

GAMMA RAY Resistivity

Potential test interval

Analysis of the results

The raw data as shown in Fig. 8 were pre-processed to allow the inversion to be carried out. First, in the absence of a strong anisotropy in the data, the average of the three readings was used. In a second step, the 'useful' part of the curve (Fig. 9) was determined (i.e. discarding of the initial part of the curve until the contact between the tool and borehole wall is firmly established). And of course, in the inversion process the portions of the signal associated with the unloading-reloading cycle are not considered. Then, the analysis of the unloading-reloading cycles allowed to deduce the value of the elastic parameters. The determination of the elastic parameters is done using the relation between the shear modulus G and the slope of the (linearized) unloading-reloading cycles:

( ) Tab. 1 : Young modulus values calculated from the test data A series of forward simulations of test #1 have been run with the objective of defining an acceptable range of variation for the parameters related to the compressibility of the sand.

0 2 1 2 R u P E G ∆ ∆ = - = ν ....
In those runs, the following values have been assumed: -Slope of the critical state line M = 1,16 (determined from laboratory testing); -Plastic compressibility λ = 0,1 -Elastic compressibility κ = 0,01 -Poisson ratio ν = 0,3 -initial void ratio = 0,43

In addition, the test conditions were supposed perfectly drained. Having set those values, a simulation was run with the derived Young modulus value, including that of the initial slope. However, because these values seemed rather low, an arbitrary value of 1000 MPa was also considered. From the results (Fig. 10) it appears that none of the values give a satisfactory match, but the initial part of the curve is better described with the lower value of 250 MPa. Although this value does not agree with the interpretation of the unloading-reloading cycles, it will nevertheless be used in the subsequent simulations. Then, the impact of the slope M of the critical state line was assessed without more success (Fig. 11). Finally the combined influence of the consolidation pressure and compressibilities (characterized by β, the difference between the plastic and elastic compressibilities) was considered. It appeared that only the hypothesis of a strong over-consolidation allowed matching the experimental data (Fig. 12). Having identified a set of parameters that gives a satisfactory simulation of the operational results, the following step should have been, in principle, the application of the inversion algorithm. However, at a significant level of overconsolidation, as identified in the preliminary study, it appears that the compressibility no longer impacts the response of the dilatometer. This can be clearly seen on Fig. 12. Therefore, after fine-tuning of the determination of the adequate consolidation pressure, it was not deemed useful to launch the full inversion scheme.

Discussion

On the dilatometer results

The results obtained here indicate a strong over-consolidation of the tested interval. This result can be interpreted in several ways. A light cementation, providing some degree of cohesion in the sand, could be causing a behavior such as recorded.

Another hypothesis for such an over-consolidation could be that a significant erosion, geological or glacial, has taken place, in which case it would imply that, at some stage in its history, the reservoir has been buried at a greater depth than the current one.

In the framework of a Cam-Clay formulation, a rapid estimation of the thickness of the eroded column can even be made under simple assumptions, which are probably somewhat restrictive, but nevertheless adequate to obtain an order of magnitude. These assumptions are:

1 No major tectonic event occurred during the whole erosion process, i.e. the total horizontal to vertical stress ratio has remained constant in time, 2 The pore pressure gradient has remained hydrostatic, 3 The current estimate of the overburden stress gradient at the current reservoir level can be extended down to the maximum historical burial depth of the reservoir determined from this calculation. Only a geological study focussed on a possible erosion would confirm and quantify the above conclusions. However, matching indications of over-consolidation of these formations was found in the analysis of logs in several wells in the same field in the form of a systematic trend in the sonic data across the shale intervals. That trend, indicating higher velocities than the normal trend can be interpreted as an indication of over-compaction, i.e. over consolidation.

Implications for the reservoir

If such a reservoir is produced by depletion, the formation undergoes compaction as the pore pressure decreases. Since the formation is over-consolidated, that compaction will take place in the elastic domain -meaning a small rate of deformation-until a state of stress is reached that corresponds to the current consolidation pressure. Beyond that point, plastic compaction -i.e. a larger rate of deformation-occurs, leading to an increased compaction drive mechanism. Assuming that the reservoir compaction takes place in the vertical direction only -i.e. oedometric deformation-, it is possible to determine the state of stress, and also the associated pore pressure, at which the transition between elastic and plastic deformation will take place. In practical terms, this indicates that compaction drive, if any, will only become significant past that threshold depletion. Depending on the current state of stress and the Over consolidation Ratio, this threshold may only be reached very late in the field life, thus meaning that the compaction drive, hence the associated additional recovery, might be negligible. The above remark will be particularly pertinent for strongly over-consolidated reservoir rocks. With this remark in mind, the impossibility to determine, with the dilatometer, the plastic compressibility, and therefore the importance of the additional recovery from compaction drive, becomes a side issue.

On the inversion scheme

In the above study, two major assumptions have been made : i) the formation properties are homogeneous over the domain investigated by the dilatometer and ii) the behavior of the formation is not time dependent (i.e. any creep phenomenon has been disregarded). It may be that none of those two assumptions are strictly valid.

Homogeneity of the mechanical properties over the investigated domain

Considering the first hypothesis, it amounts to implicitly deny the presence of a filter cake of significant thickness and compressibility and assume that, despite the minute operational clearance available, the positioning of the tool can be achieved without significant mechanical damage being caused to the borehole wall. In the present case, in view of the low in-situ mobility of the pore fluid, it can be assumed that no massive filtration has taken place during the drilling phase and, consequently, that no thick filter cake has been deposited. In addition, given the low tool clearance (less that 3 mm), it is unlikely that any filter cake left after the tool is in place, will be significant enough to affect the test results. As to the presence of a near wellbore mechanical damage, the minimal clearance available when maneuvering to set the dilatometer in place, can be of consequence. Assuming as a fact that a mechanically 'disturbed' zone is generated, the question is then to integrate it into the inversion process which poses the problem of identifying the physical parameters both for the damaged and the intact formation. The resolution of the induced complication would require further work.

Time dependant behavior of the formation

The time dependant behavior of the formations tested has been evidenced during the initial unloading-reloading cycles of the test and had to be performed at pressure rates higher than the rest of the test to shadow the time dependant behavior of the formation. This phenomenon is not taken into account in the inversion procedure. To remedy this, it would be necessary to incorporate into the inversion algorithm the possibility to model a visco-plastic behavior, whereas only elasto-plasticity is currently considered. There are no direct means of evaluating the impact of such a time dependant behavior on the data recorded during the tests. However, disregarding the presence of a time dependant component in the signal would lead to an underestimation of the consolidation pressure, which, in any case, would push the occurrence of plastic compaction even later in the life of the field. Nevertheless, the issue of time dependent behavior ought to be incorporated in the inversion scheme.

Conclusions

A pilot operation to measure in-situ the compressibility of a reservoir rock using a dilatometer has been performed in an oil field. To our knowledge, it is the first time this type of measurement has been attempted.

Operational aspects

To achieve the success of the operation, the skill and commitment of all the parties involved in the planning, design and execution of this non-conventional operation have to be emphasized as a major component of this accomplishment. The highlights of the operations were:

An original drilling program:

An 'S'-shape pilot hole was drilled and logged to identify potential test intervals. The actual test hole could be drilled as a side-track to the pilot hole. The drilling was stopped at the roof of each test interval, which was cored to ensure maximum borehole quality.

An audacious coring operation:

Despite severe operating conditions, given the nature of the formation, the non-standard tool size and the constraint for a high quality borehole, the coring operation resulted in a test hole of extremely high quality.

Tool maneuvering in difficult conditions:

The running of the dilatometer in and out of the test hole required extreme care from the drilling crew given the exceptionally small clearance available which, given the low strength of the formation, made any swab and surge effect potentially dangerous for the integrity of the test hole.

Deployment a civil engineering tool in an oilfield situation:

Because of the dissimilarities between the operational situations in geotechnical and oilfield engineering, a number of technical issues have had to be solved to achieve the test objectives. More important, the drilling crew and dilatometer personnel showed much reciprocal curiosity and provided a great deal of mutual support, which resulted in a seamless operation.

Analytical aspects

The key conclusion of the analysis of the results of the dilatometer data is that the tested sands exhibit some degree of over-consolidation. The major consequence of that overconsolidation could be that plastic compaction of the reservoir may only occur after a significant amount of depletion. It must be noted that the actual value of the plastic compressibility could not be determined. The reason for that is that, at the interpreted degree of over-consolidation, the results of the dilatometer are insensitive to the compressibility of the formation. The conclusion about the occurrence of an over-consolidation is corroborated by the findings of another on-going, independent study on the compaction of the filed shales. However, that later study does not allow yet to confirm the degree of over-consolidation obtained here. A number of recommendations can be formulated in the light of our current experience with the interpretation of dilatometer tests for in-situ measurement of formation compressibility in low strength reservoir rocks: -The time dependant behavior of the formation sands should be taken into account in the inversion methodology. Note that this may impact the operational procedure. -The model used to invert the dilatometer test data assumes homogeneous mechanical properties over the domain of investigation of the dilatometer, thus disregards the possible impact of a thick, compressible filter cake as well as that of a locally damaged formation (e.g. during the drilling or dilatometer operations). Finally, the operation presented here has been performed in a 'logistically simple' situation (onshore, shallow well). Would such an endeavor be envisaged in more 'demanding' environments involving weakly cemented reservoir rocks (typically deep off-shore), a significant engineering effort in the dilatometer equipment ought to be considered.

Figure 1 :

 1 Figure 1 : Schematics of a dilatometer

Fig. 3 :

 3 Fig. 3 : Model geometry and boundary conditions for feasibility study

Fig. 4

 4 Fig. 4 Impact of consolidation pressure on dilatometer response

Fig. 5 :

 5 Fig. 5 : Flow chart of the inversion procedure

Fig. 6 :Fig. 8 :

 68 Fig. 6 : MWD data in the pilot hole for dilatometer interval selection

Fig. 9 :

 9 Fig. 9 : 'Useful' part of the dilatometer results

Fig. 10 :

 10 Fig. 10 : Impact of Young modulus on direct test simulations

Fig. 11 :

 11 Fig. 11 : Impact of the value of M on direct test simulations

Fig. 12 :

 12 Fig. 12 : Impact of the value of the consolidation pressure on direct test simulations

st cycle 2 nd cycle 3 rd cycle Final unloading

  .................................... (A.1) MPa was found, except for the initial part of the curve, from which a value of 220 MPa was deduced. Although the results were fairly consistent from one cycle to the other, the values found were unexpectedly low, considering the depth of the tested interval.

	Where u is the radial displacement, and
	R 0 the initial wellbore diameter
	Using each of the unloading cycles (including the final one)
	the following values of G are derived. The corresponding
	values of the Young modulus, assuming a 0.3 Poisson ratio,
	are calculated (Tab. 1). A Young modulus value ranging from
	640-720 Initial slope 1 Shear modulus 84	247	247	271	229
	(MPa)
	Young Modulus	219	642	642	705	595
	(MPa)

Acknowledgement

The authors would like to acknowledge permission to publish the paper from Total and SINCOR management. The help and enthusiasm of the many people involved in this project has been a key element in the successful achievement of this operation.


Zined Marzouk
email: zinedmarzouk@hotmail.com

Didier Aurelle

Anne Chenuil

Cryptic lineages and high population genetic structure in the exploited marine snail Hexaplex trunculus (Gastropoda: Muricidae)

Keywords:

Although the banded murex Hexaplex trunculus (Linnaeus, 1758) has no dispersal stage, it is widely distributed in a relatively broad range of habitats. These features make it a particularly suitable model to reconstruct the history of the fragmentation of its geographical range. We investigated its genetic structure from the eastern Mediterranean to the Atlantic coast, by sequencing a nuclear intron (i29) and a mitochondrial marker (cox1). We found strong genetic differentiation between all population pairs, congruent with the absence of a dispersing larval stage. A deep phylogeographical break separated two parapatric lineages, a western lineage (I) and an eastern one (II). The two lineages were separated by a vicariance event dated to the Pliocene by both markers, analysed independently. They co-occur in southern Italy and in the Siculo-Tunisian straight, where some individuals display recombined genotypes (lineage I for one marker, lineage II for the other), suggesting that introgression occurred in sympatric populations. We were unable to determine whether the vicariance was across the STS or located further east across the Peloponnese Arc, but the presence of lineage I in the northern Adriatic is more parsimoniously explained by a more eastern vicariance. Lineage I displayed a stronger signal of demographic expansion than lineage II, and its expansion was estimated to be more recent. This result, which has been reported in other marine species, suggests less drastic conditions for the eastern (and possibly also the central) Mediterranean benthos than for western populations during past climatic oscillations. The parapatric distribution still observed today suggests that human exploitation of this snail, which dates back to Antiquity, did not result in efficient introductions among basins, although present-day introductions were recently reported in the Bay of Biscay.

INTRODUCTION

The marine environment is considered to be fertile for cryptic speciation (i.e. species that are not distinguishable morphologically), suggesting that marine biodiversity has been underestimated at the species level [START_REF] Knowlton | Sibling species in the sea[END_REF][START_REF] Féral | How useful are the genetic markers in attempts to understand and manage marine biodiversity?[END_REF]. Several genetic studies of marine organisms in the Mediterranean Sea and the north-eastern Atlantic Ocean have revealed the existence of cryptic species complexes, highlighting the evolutionary processes that shaped the genetic diversity within these complexes (e.g. [START_REF] Boissin | Defining reproductively isolated units in a cryptic and syntopic species complex using mitochondrial and nuclear markers: the brooding brittle star, Amphipholis squamata (Ophiuroidea)[END_REF][START_REF] Calvo | Hidden Mediterranean biodiversity: molecular evidence for a cryptic species complex within the reef building vermetid gastropod Dendropoma petraeum (Mollusca: Caenogastropoda)[END_REF][START_REF] Krebes | Glacial-driven vicariance in the amphipod Gammarus duebeni[END_REF]Boissin, Hoareau & Berrebi, 2011a;[START_REF] Modica | DNA-barcoding of sympatric species of ectoparasitic gastropods of the genus Cerithiopsis (Mollusca: Gastropoda: Cerithiopsidae) from Croatia[END_REF][START_REF] Durand | Mitochondrial phylogeny of grey mullets (Acanthopterygii: Mugilidae) suggests high proportion of cryptic species[END_REF].

The biogeography of the Mediterranean sea has been shaped by a complex geological history since the closure of the Tethys (during the Miocene), passing through the Messinian Salinity Crisis, the Pliocene inundation from the Atlantic Ocean and Pleistocene glacial episodes [START_REF] Patarnello | Pillars of Hercules: is the Atlantic-Mediterranean transition a phylogeographical break?[END_REF][START_REF] Sabelli | The making of the Mediterranean Molluscan biodiversity[END_REF]and references therein). Traces of past Tethys distribution ranges are reflected in the present distributions of several taxa, either in the distinct species of a genus (e.g. Posidonia, [START_REF] Aires | Evolutionary history of the seagrass genus Posidonia[END_REF] or between cryptic species of a given nominal species [START_REF] Egea | Morphological and genetic analyses reveal a cryptic species complex in the echinoid Echinocardium cordatum and rule out a stabilizing selection explanation[END_REF]. Independently of past events, numerous studies report a role of contemporary physical barriers in shaping the distribution of the genetic diversity of marine organisms in the Mediterranean Sea and the north-eastern Atlantic Ocean [START_REF] Duran | Genetic diversity and population structure of the commercially harvested sea urchin Paracentrotus lividus (Echinodermata: Echinoidea)[END_REF][START_REF] Patarnello | Pillars of Hercules: is the Atlantic-Mediterranean transition a phylogeographical break?[END_REF][START_REF] Maggs | Evaluating signatures of glacial refugia for North Atlantic benthic marine taxa[END_REF][START_REF] Zitari-Chatti | Mitochondrial DNA variation in the caramote prawn Penaeus (Melicertus) kerathurus across a transition zone in the Mediterranean Sea[END_REF][START_REF] Maltagliati | Phylogeography and genetic structure of the edible sea urchin Paracentrotus lividus (Echinodermata[END_REF][START_REF] Tarnowska | Complex genetic population structure of the bivalve Cerastoderma glaucum in a highly fragmented lagoon habitat[END_REF]Boissin, Stöhr & Chenuil, 2011b;[START_REF] Fadhlaoui-Zid | Genetic structure of Octopus vulgaris (Cephalopoda, Octopodidae) in the central Mediterranean Sea inferred from the mitochondrial COIII gene[END_REF]. Such barriers are created by distinct hydrological and ecological conditions that may be sufficient to reduce gene flow even on a relatively small geographical scale within seemingly continuous populations [START_REF] Hellberg | Gene flow and isolation among populations of marine animals[END_REF]. These barriers to dispersal often coincide with contact zones between differentiated lineages, which can be explained by the dynamics of hybrid zones causing their displacement and stabilization in ecotones or lower carrying capacity areas [START_REF] Barton | Dynamics of hybrid zones[END_REF].

Species life-history traits such as dispersal ability may have a major influence on the level and structuring of genetic diversity [START_REF] Weber | Influence of the larval phase on connectivity: strong differences in the genetic structure of brooders and broadcasters in the Ophioderma longicauda species complex[END_REF], but see also [START_REF] Boissin | Contrasting population genetic structures in Amphipholis squamata, a complex of brooding and self-reproducing sister species sharing lifehistory traits[END_REF]. A strong dispersal potential enhances genetic homogenization among populations throughout the distribution range of a species, whereas restricted dispersal (absence of larval phase or nektobenthic larvae, and sessile or sedentary adults) may lead to strong population differentiation and increased speciation rates [START_REF] Kelly | Genetic structure among 50 species of the north eastern Pacific rocky intertidal community[END_REF]Boissin et al., 2011b). Oceanographic barriers do not affect all marine species in the same way, and some species do not display deep phylogeographical breaks [START_REF] Patarnello | Pillars of Hercules: is the Atlantic-Mediterranean transition a phylogeographical break?[END_REF][START_REF] Aurelle | Phylogeography of the red coral (Corallium rubrum): inferences on the evolutionary history of a temperate gorgonian[END_REF]. The link between life-history traits and the impact of these barriers [START_REF] Patarnello | Pillars of Hercules: is the Atlantic-Mediterranean transition a phylogeographical break?[END_REF], as well as the correlation between potential dispersal and observed levels of genetic differentiation, are not straightforward [START_REF] Selkoe | Marine connectivity: a new look at pelagic larval duration and genetic metrics of dispersal[END_REF]. It is often difficult to distinguish the impact of habitat fragmentation from that of low dispersal ability, as both contribute to increased genetic differentiation (e.g. [START_REF] Tarnowska | Comparative phylogeography of two sister (congeneric) species of cardiid bivalve: strong influence of habitat, life history and post-glacial history[END_REF]. Species displaying narrow habitats and high dispersal are relatively rare, like low dispersers with wide niches.

The muricid gastropod, Hexaplex trunculus, is a common inhabitant of the Mediterranean Sea and adjacent Atlantic coasts (Madeira and Canary archipelagos, and the Moroccan and Portuguese coasts). This species occurs in intertidal and infralittoral zones between 1-and 100-m depth. It inhabits both hard and soft substrates, from rocky to sandy-muddy shores [START_REF] Chiavarini | Butyltins concentration levels and imposex occurrence in snails from the Sicilian coasts (Italy)[END_REF]. Hexaplex trunculus is a gonochoric species with collective spawning, which produces intracapsular development larvae (no planktonic larval stage). It is also highly restricted in mobility at the adult stage, which potentially limits its colonization ability [START_REF] Vasconcelos | Spawning of Hexaplex (Trunculariopsis) trunculus (Gastropoda: Muricidae) in the laboratory: description of spawning[END_REF]. This species has been harvested by humans since prehistoric times [START_REF] Fernandez-Lopez De Pablo | El Collado shell midden and the exploitation patterns of littoral resources during the Mesolithic in the Eastern Iberian Peninsula[END_REF]. It had high commercial value since Antiquity as a dye provider [START_REF] Valenzuela | An ancient fishery of Banded dye-murex (Hexaplex trunculus): zooarchaeological evidence from the Roman city of Pollentia (Mallorca, Western Mediterranean)[END_REF], and more recently for human consumption. Recent studies found it a potential source of gelatin for industry [START_REF] Zarai | Process for extracting gelatin from marine snail (Hexaplex trunculus): chemical composition and functional properties[END_REF], and it is a sensitive bioindicator of organic pollution via imposex [START_REF] Anastasiou | Imposex and organotin compounds in ports of the Mediterranean and the Atlantic: is the story over?[END_REF]. Given its wide geographical distribution, relatively low habitat fragmentation (unlike intertidal gastropods) and its highly restricted dispersal ability, we considered H. trunculus to be an ideal candidate for studying physical and evolutionary factors inducing population genetic structure in the Mediterranean Sea and the north-eastern Atlantic Ocean.

To date, few studies have been carried out on the population genetics of H. trunculus. [START_REF] González-Tizón | Genetic diversity in fishery-exploited populations of the banded murex (Hexaplex trunculus) from the southern Iberian Peninsula[END_REF] compared the genetic diversity of the harvested population from Ria Formosa (southern Portugal) with two adjacent populations in order to provide source populations for stock enhancement. [START_REF] Villamor | Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates[END_REF] characterized Italian populations from the Ligurian, Tyrrhenian, Ionian and Adriatic Seas, encompassing putative biogeographical barriers for the cox1 marker, and reported strong genetic differentiation across the Siculo-Tunisian strait, with a potential gene flow from Ligurian-Tyrrhenian towards Ionian-Adriatic populations. However, these studies were based on a single mitochondrial marker and/or were geographically too restricted to address the genetic variability of H. trunculus across its distribution range or to identify major genetic divisions. Recently, a morphometric study, based on shell characters (such as shell height, aperture length etc.) from 17 locations throughout the whole distribution area of the species strongly suggested the presence of four main groups corresponding to different geographical areas: the Atlantic; the Alboran Sea; the western Mediterranean (France, northern Tunisia and the Siculo-Tunisian strait); and an extensive eastern Mediterranean group grouping eastern Tunisian, Adriatic and Greek populations [START_REF] Marzouk | Morphometric variation of fishery-exploited Muricidae (Hexaplex trunculus) in the Mediterranean Sea and the Northeastern Atlantic Ocean: Implications for stock identification[END_REF]. This structure was evidenced by multidimensional scaling (MDS), hierarchical clustering analysis and canonical discriminant analysis. All morphometric indices were significantly differentiated between groups. Since morphology may be influenced by environment, through developmental plasticity, genetic markers were necessary to establish the validity of this phylogeographical pattern and, in the case of long-term divergence, to estimate its age. Such data are also crucial to optimize the management of this exploited resource.

In the present study, we investigated the genetic structure of H. trunculus at 19 localities, including those of the morphological study [START_REF] Marzouk | Morphometric variation of fishery-exploited Muricidae (Hexaplex trunculus) in the Mediterranean Sea and the Northeastern Atlantic Ocean: Implications for stock identification[END_REF], for a mitochondrial marker (cox1) and a nuclear intron (i29). In addition, we included data from the six locations characterized by [START_REF] Villamor | Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates[END_REF] and one cox1 sequence of a specimen from Haifa (Levantine Sea) retrieved in the BOLD-system database [START_REF] Ratnasingham | A DNA-based registry for all animal species: the barcode index number (BIN) system[END_REF]. On the basis of extensive sampling coverage, we investigated (1) the main phylogeographical splits and the possibility of distinct species and secondary contacts within the nominal species H. trunculus, (2) its genetic structure at a smaller geographical scale and (3) its demographic history.

MATERIAL AND METHODS

Sample collection

Samples were collected from 19 localities across most of the distribution range of this species (Fig. 1; Table 1). All specimens were collected between 0.5-and 30-m depth from January 2009 to July 2013. Collected samples were consistent with the taxonomic description of the species by [START_REF] Houart | A review of the recent Mediterranean and Northeastern Atlantic species of Muricidae[END_REF]. Following capture, samples were transported to the laboratory, and a small portion of foot tissue was cut and preserved in absolute ethanol at -20 °C until further processing. A hundred published cox1 sequences (GenBank accessions: KF297407-KF297433) from six Italian populations and a single specimen from 1 for details), Letters (A→F) are the six populations studied by [START_REF] Villamor | Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates[END_REF]. Pie charts show the distribution of the major haplogroups for cox1 (A 1 : present study, A 2 : [START_REF] Villamor | Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates[END_REF] and for i29 (B). For the cox1 sequences of [START_REF] Villamor | Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates[END_REF], it was not possible to distinguish haplogroups Ia and Ib due to shorter sequence length. ATO, Atlantic Ocean; ALS, Alboran Sea; WMB, Western Mediterranean Basin; STS, Siculo-Tunisian Strait; EMB, Eastern Mediterranean Basin; ADS, Adriatic Sea; AGS, Aegean Sea. Haifa (Levantine Sea) (BOLD record ID: 3443789) were included in our analysis (Fig. 1).

Dna extraction, amplification anD Sequencing

Total genomic DNA was extracted with Chelex resin (Bio-Rad) solution at 10% [START_REF] Walsh | Chelex 100 as a medium for simple extraction of DNA for PCR-based typing from forensic material[END_REF] and preserved at -20 °C. PCR was used to amplify (1) 710 bp of the mitochondrial Cytochrome Oxydase I gene (cox1), using the primers LCO1490 and HCO2198 from [START_REF] Folmer | DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates[END_REF] and ( 2) 581 bp of the EPIC (Exon Primed Intron Crossing) marker i29 that was successfully amplified in gastropods and bivalves in a recent study [START_REF] Gérard | PCR survey of 50 introns in animals: cross-amplification of homologous EPIC loci in eight non-bilaterian, protostome and deuterostome phyla[END_REF].

For i29, we used the primers HTi29F1 and HTi29R2, which were designed from sequence data of the fragment, amplified with universal primers developed by [START_REF] Chenuil | An efficient method to find potentially universal population genetic markers, applied to metazoans[END_REF]. The sequences of the primers are HTi29F1: 5′-AGAGACAGACACACTGTCTGAC-3′ and HTi29R2: 5′-TCCGTGGTCTCCCACACACAAT-3′. The cox1 marker exhibits levels of sequence variation suitable for population analysis of numerous marine invertebrates [START_REF] Reuschel | Phylogeny and geographic differentiation of Atlanto-Mediterranean species of the genus Xantho (Crustacea: Brachyura: Xanthidae) based on genetic and morphometric analyses[END_REF][START_REF] Pérez-Losada | Testing hypotheses of population structuring in the Northeast Atlantic Ocean and Mediterranean Sea using the common cuttlefish Sepia officinalis[END_REF][START_REF] Calvo | Hidden Mediterranean biodiversity: molecular evidence for a cryptic species complex within the reef building vermetid gastropod Dendropoma petraeum (Mollusca: Caenogastropoda)[END_REF][START_REF] Zitari-Chatti | Mitochondrial DNA variation in the caramote prawn Penaeus (Melicertus) kerathurus across a transition zone in the Mediterranean Sea[END_REF][START_REF] Penant | Planktonic larvae do not ensure gene flow in the edible sea urchin Paracentrotus lividus[END_REF]; but see [START_REF] Calderón | Evaluation of the utility of COI and ITS markers as tools for population genetic studies of temperate gorgonians[END_REF]. The nuclear intron i29, like many other EPIC markers, is conserved in position across many metazoans, yet is potentially highly polymorphic [START_REF] Chenuil | An efficient method to find potentially universal population genetic markers, applied to metazoans[END_REF]. PCR were conducted in a Mastercycler Epgradient S from Eppendorf using the following programs. (1) For cox1, initial denaturation 3 min at 94 °C, 35 cycles, denaturation 1 min at 94 °C, primer annealing for 1 min at 42 °C, extension for 1 min at 72 °C and final elongation for 5 min at 72 °C PCR reactions were carried out in 25-μL reaction volumes, with 1× buffer, 2 mM MgCl 2 , 0.2 mM dNTP, 2 μM of each primer, 1 U Taq polymerase (Promega) and 2 μL of DNA. (2) For i29, the PCR programme was similar as for cox1, but the annealing temperature was 50 °C and the final elongation step lasted 8 min. PCR products were sequenced on ABI automated sequencers at the LGC Genomics Firm using the LCO1490 or the HTi29F1 primers. For each marker, a first set of 24 individuals was also sequenced with the opposite primer. On the basis of the results (see below), we decided to use a single primer for sequencing.

phylogenetic analySeS anD Dating

Sequences were aligned using ClustalW [START_REF] Thompson | CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice[END_REF] implemented in Bioedit [START_REF] Hall | Bioedit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/ NT[END_REF]. The DNA sequences were deposited in the GenBank database with the following accession numbers: KJ135221-KJ135286 and KT934232-KT934250 for cox1 and KT934104-KT934231 for i29.

We used the NETWORK software [START_REF] Bandelt | Median-joining networks for inferring intraspecific phylogenies[END_REF] to reconstruct haplotype networks. We used SplitsTree4 [START_REF] Huson | Application of phylogenetic networks in evolutionary studies[END_REF] to perform the Phi-test of the hypothesis that recombination occurred in i29 and to compare several methods of network reconstruction. Phylogenetic reconstructions were performed using (1) maximum likelihood algorithms in Phyml online at the ATGC Montpellier bioinformatics platform [START_REF] Guindon | PhyML Online a web server for fast maximum likelihood based phylogenetic inference[END_REF] and (2) Bayesian inference (BI) in BEAST [START_REF] Drummond | BEAST: Bayesian evolutionary analysis by sampling trees[END_REF]. All analyses were performed using the HKY + I + G model [START_REF] Hasegawa | Dating of the human-ape splitting by a molecular clock of mitochondrial DNA[END_REF], as determined by the JModelTest [START_REF] Posada | jModelTest: phylogenetic model averaging[END_REF]. Sequences of Bolinus brandaris (GenBank accession numbers: DQ280020 and KU821678) and Nucella lapillus (GenBank accession number: EU391582) were used as outgroups. A hundred bootstrap replicates were performed for Phyml. BI reconstruction implemented in BEAST v1.7.2 [START_REF] Drummond | BEAST: Bayesian evolutionary analysis by sampling trees[END_REF] was used under a strict clock and Yule (pure birth) process as a model of speciation to estimate TMRCA (time to the most recent common ancestor) and their credibility intervals (95% HPD). To calibrate the cox1 analyses, we used a fossil calibration of 40 Mya for the divergence between N. lapillus and the clade grouping Hexaplex/Bolinus, following Merle, Garrigues & Pointier (2011 and pers. commun.). Then, the divergence time between Hexaplex and Bolinus deduced from the BEAST phylogenetic analyses on cox1 was used to calibrate the i29 analyses.

Posterior distributions of parameters were approximated using five independent MCMC analyses of 20 million generations with 10% burn-in. Samples from the two runs, which yielded similar results, were combined, and convergence of the chains was checked using the program Tracer v1.6 [START_REF] Drummond | BEAST: Bayesian evolutionary analysis by sampling trees[END_REF]. Results were analysed and visualized with TreeAnnotator v1.7.2, using the maximum clade credibility tree and mean node heights, and with FigTree v1.4 [START_REF] Rambaut | FigTree version 1.3.1[END_REF].

genetic DiverSity anD Structure

Nucleotide and haplotype diversities (π, Hd) were calculated using the DNAsp software [START_REF] Librado | DnaSP v5: a software for comprehensive analysis of DNA polymorphism data[END_REF]. Differentiations between population pairs were estimated with pairwise F ST and Φ ST statistics. ARLEQUIN was used [START_REF] Excoffier | Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows[END_REF] to compute F ST and Φ ST , and the null hypothesis of no differentiation was tested by performing 10 000 permutations. The [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] correction for multiple tests was applied. Kimura two-parameters (K2P) mean genetic distances between lineages were computed, for both markers, in MEGA5 [START_REF] Tamura | MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods[END_REF].

mDS baSeD on genetic DiStanceS between populationS

Ordinal MDS analyses, based on Φ ST pairwise genetic distances between populations, were performed using XLSTAT (AddinSoft, 2007) to visualize relationships among populations with the two distinct markers. Populations were coloured according to their morphological group assignment, which was spatially interpolated for the two populations without morphological data [START_REF] Marzouk | Morphometric variation of fishery-exploited Muricidae (Hexaplex trunculus) in the Mediterranean Sea and the Northeastern Atlantic Ocean: Implications for stock identification[END_REF], Corsica and Tabarka.

multilocuS analySeS of inDiviDualS in areaS of

Sympatry among mitochonDrial haplogroupS

Since very divergent haplogroups were found in cox1 (see Results), we tested the hypothesis that those haplogroups corresponded to separate biological species [START_REF] Mayr | Systematic and the origin of species[END_REF]. Under this hypothesis, when in sympatry, these species should not exchange genes efficiently (hybrids should be unviable or sterile). Thus, genetic differentiation for the nuclear marker should be observed between mitochondrial groups within populations. We thus built, for each population wherein distinct haplogroups coexisted, contingency tables splitting individuals based on their cox1 and i29 lineages. We then tested the hypothesis that i29 lineages were independent of cox1 lineages using Fisher's exact tests.

Demographic hiStorieS

Mismatch distributions of pairwise nucleotide differences [START_REF] Rogers | Population growth makes waves in the distribution of pairwise genetic differences[END_REF] in the main geographic regions (ATO, ALS, WMB, STS, EMB, ADS and AGS) were built using DNAsp. For neutral polymorphisms, when population sizes are constant, these distributions are generally chaotic, whereas in case of past demographic expansion, a unimodal curve is observed. Mismatch distributions were also performed separately for the main lineages found for each marker. We used ARLEQUIN [START_REF] Schneider | Estimation of past demographic parameters from the distribution of pairwise differences when the mutation rates vary among sites: application to human mitochondrial DNA[END_REF] to compute the tau estimator and its confidence interval for the cox1 marker. The tau statistic (τ), under a model of demographic expansion, estimates the expansion's starting time in mutational units, using the formula of [START_REF] Li | Distribution of nucleotide differences between two randomly chosen cistrons in a finite population[END_REF]: τ = 2 μt, where μ is the mutation rate for the whole sequence and t is the time in generation units. Mismatch distribution tests were reported to be conservative [START_REF] Ramos-Onsins | Statistical properties of new neutrality tests against population growth[END_REF], thus we also computed Tajima's D tests [START_REF] Tajima | Statistical testing for the neutral mutation hypothesis by DNA polymorphism[END_REF]), Fu's Fs [START_REF] Fu | Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection[END_REF] and Ramos-Onsins and Rozas's R 2 tests [START_REF] Ramos-Onsins | Statistical properties of new neutrality tests against population growth[END_REF], implemented in DNAsp, for the cox1 marker. Significance levels for these tests were calculated using 1000 coalescent simulations. When these parameter values depart from zero, the null hypothesis of demographic stability and selective neutrality of the cox1 marker is rejected.

RESULTS genetic DiverSity

A total of 343 H. trunculus cox1 sequences of 658 bp were obtained. Among them, 85 different haplotypes were identified; 51 sites were variable and 39 were parsimony informative. For the nuclear marker i29, 339 specimens were successfully sequenced for 500 bp, displaying 128 distinct haplotypes, 58 polymorphic sites, of which 45 were parsimony informative. Only one seventh of superimposed peaks were obtained for i29, after Sanger amplicon sequencing, probably due to individual heterozygosity at this marker. Most individuals displayed no ambiguity for their i29 sequence. Some i29 haplotypes were thus probably not detected. For this reason, data analyses did not use individual diploid genotypes at i29. Polymorphism parameters are reported in Table 1 for each marker and each population.

haplotype networkS, phylogenetic analySeS anD DatationS

Haplotype networks and phylogenetic trees evidenced strong phylogeographical structures for both markers. The cox1 network showed that the haplotypes were distributed into two main haplogroups (Fig. 2). Using cox1 sequences of two outgroups, B. brandaris and N. lapillus, phylogenetic analyses recovered two lineages perfectly matching the two haplogroups (I and II). These main lineages were supported by bootstrap values higher than 80% and posterior probabilities higher than 0.97 (Fig. 3). The haplotype network showed that the haplogroup I could be divided into two groups of haplotypes: Ia, which was found in the Atlantic, Alboran, Marseilles, Corsica, Tunisia (except M. Bourguiba) and the Adriatic, and Ib, which was only found in Corsica, Tunisia and the Adriatic. The first lineage or haplogroup (I) displayed a star-like pattern (Fig. 2). The second haplogroup (II) was mainly composed of unique haplotypes from Aegean and Siculo-Tunisian samples and a few haplotypes shared between these two regions (Figs 1,2). The two cox1 lineages (I cox1 and II cox1 ) were separated by seven mutations corresponding to a divergence of at least 1.06% among haplotypes between those lineages. The mean K2P genetic distance was estimated at 0.024 between cox1 lineages and 0.133 between B. brandaris and the two cox1 lineages. The divergence time between the two lineages (I cox1 and II cox1 ) was estimated at 5.18 Mya (95% HPD 3.5-7.3), corresponding to the Pliocene. The estimated cox1 substitution rate mean value was 0.00295 substitutions per Mya (SD = 0.000418) corresponding to the divergence rate of 0.6% per Mya. This was different from the cox1 divergence rates estimated by [START_REF] Lessios | The great American schism: divergence of marine organisms after the rise of the Central American Isthmus[END_REF] and [START_REF] Miura | Molecular phylogenetics reveals differential divergence of coastal snails separated by the Isthmus of Panama[END_REF] (2.4-3.3% per Mya), based on a marine gastropod in the Isthmus of Panama. Applying this rate provided a divergence time of about 1 Mya for the two H. trunculus lineages. When we integrated the cox1 sequences published by [START_REF] Villamor | Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates[END_REF] and by BOLD-system database [START_REF] Ratnasingham | A DNA-based registry for all animal species: the barcode index number (BIN) system[END_REF] in our data set, the distinction between haplogroups Ia and Ib was lost because Villamor's fragment was smaller than ours (we truncated our sequence file to the fragment size of [START_REF] Villamor | Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates[END_REF]. We identified the two main haplogroups: I cox1 that occurred in the six Italian localities (A, B, C, D, E and F) and II cox1 occurring in the Tyrrhenian (C), Ionian (D and E), Adriatic (F) and Levantine (Haifa) Seas (Supporting Information, Appendix S1; Fig. 1).

The i29 network showed two lineages (I and II) and three haplogroups (I-1, I-2 and II), which appeared less divergent (as % of sequence differences) than the two mitochondrial ones (Fig. 2). The divergence value among i29 lineages is 0.4%. The mean K2P genetic distance was estimated at 0.019 between i29 lineages and 0.135 between B. brandaris and i29 lineages.

As illustrated by Figure 2, the phylogeographical patterns are very similar between cox1 and i29, in particular: (1) Aegean populations display a single haplogroup and (2) this haplogroup is found nowhere else except in the two populations of the Siculo-Tunisian strait (STS). A particular feature for i29 is the presence of a group of haplotypes (I-2 i29 ), which is only found in the STS. Using phylogenetic trees for this marker, we recovered, with low node supports, the distinction between two main lineages (I i29 and II i29 ) and their divergence time was estimated at 5.21 Mya (95% HPD 3.4-7.5) (Fig. 4), very similar to the divergence between I cox1 and II cox1 . Haplogroups I-1 i29 and I-2 i29 formed distinct clades in phylogenetic trees but were poorly supported.

In Tunisia, where the two main lineages are in contact, we investigated whether haplogroup I-2 i29 , which is restricted to the Siculo-Tunisian strait, may have been generated by a recombination event between haplogroups I-1 i29 and II i29 , which are both present in Tunisian populations. For this, we compared the three haplogroups and checked in the sequence alignment whether one of them may have resulted from a Colours indicate the origin of the individuals displaying each haplotype, as in the legend: to limit the use of colours, spatially close populations were combined according to the regions defined in Table 1, except for the four populations from the WMB region because (1) they may be very distant and (2) two very close localities (northern Tunisia) do not fall into the same haplogroup (for cox1). Black circles represent unsampled haplotypes. recombination between the other two: no pair of sites supported this hypothesis, thus recombination events cannot explain the peculiar geographical distribution of i29 haplogroups. The Phi-tests, both in the total data set and in the area of sympatry, were not significant, supporting the hypothesis that no recombination occurred in this DNA region.

population Differentiation

For cox1 and i29, all pairwise comparisons between populations based on F ST , and nearly all of those based on Φ ST , remained significant after correction for multiple tests [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF] (Supporting Information, Appendix S2). F ST values were high for both markers (about 0.301 for cox1 and 0.286 for i29), even between pairs of adjacent populations [e.g. F ST (cox1) > 0.6 between the two Atlantic sites and F ST (cox1) > 0.3 between adjacent Moroccan populations]. Overall, Φ ST was about 0.419 and 0.438 for cox1 and i29, respectively. mDS illustrated that genetic differentiation patterns were broadly compatible with morphological grouping (which were highly congruent with geography), although some differences were visible (Supporting Information, Appendix S3). The genetic difference between the Atlantic (with or without the Alboran Sea) and the Mediterranean Sea (excluding or not the Alboran Sea) was not salient, in contrast with the differentiation of Aegean localities vs. all other localities, which was striking for cox1, and very clear for i29. The central position of samples from the STS in the MDS reflected the fact that they are a mixture of both cryptic lineages (species), so they are intermediate between Aegean Sea and other populations for both markers.

multilocuS teStS for cox1 anD i29 haplogroupS aSSociation in inDiviDualS

We built contingency tables to test the independence of the two genetic markers' major splits. This test was done within populations (to avoid the effect of genetic structure due to restricted dispersal) where the haplogroups I and II coexisted (Haouaria and Kelibia). In those tables, we distributed the individuals according to their belonging to major lineages for the cox1 marker (raws) or i29 (columns) (Supporting Information, Appendix S4). Some heterogeneous or 'recombined' individuals displayed the eastern haplogroup (named II) for one marker and the western haplogroup (named I) for the other marker. The probability of obtaining a more skewed distribution (onesided Fisher's exact test) was 0.031* for Haouaria and 0.049* for Kelibia, but the Fisher's exact tests were not significant (probability of obtaining a more or equally skewed distribution, respectively, 0.521 and 0.237). We thus cannot reject the null hypothesis that genes are randomly associated between these lineages in this area.

neutrality teStS anD miSmatch DiStributionS

Mismatch distributions within regions were generally unimodal for both markers, suggesting population expansions (or past selective sweeps), except for the STS where the coexistence of distinct cox1 lineages generated clear-cut bimodality (one mode for intra-lineage, one for inter-lineage comparisons). The fact that bimodality was not observed for the i29 marker in the same region may be explained by the lower genetic divergence between lineages than for cox1. The mismatch distribution did not appear unimodal (nor unambiguously bimodal) for the Aegean Sea in i29 (Supporting Information, Appendices S5, S6; Fig. 2), a result compatible with a more stable demography. Numerous significant negative values of Fs and Tajima's D (but R 2 values were often not significant) in cox1 (Table 2) rule out the hypothesis of constant effective sizes and neutral polymorphism and suggest either selective sweeps or past demographic expansion. Since past population expansion is likely, it is of interest to examine tau parameters and corresponding expansion times (Table 2). They suggested a more recent and/or more extensive expansion in the lineages I cox1 and I-1 i29 compared to lineages II (much more star-like haplotype networks). The most recent expansion was estimated for the Adriatic Sea (about 88 thousand years ago).

DISCUSSION

populationS are genetically DifferentiateD even at Small geographical Scale Our analyses based on extensive geographical sampling evidenced strong genetic differentiation among populations of H. trunculus at different geographical scales. At small geographical distance, differentiation was revealed by high and significant F ST values among all studied population pairs, even between the geographically closest ones (about 60 km between Moroccan populations and about 20-40 km between Tunisian populations). Since the habitat of H. trunculus is not particularly fragmented, reduced dispersal ability is the most likely cause of such a genetic structure. The benthic habit of the adults [START_REF] Chiavarini | Butyltins concentration levels and imposex occurrence in snails from the Sicilian coasts (Italy)[END_REF], the benthic eggs hatching directly to nektobenthic juveniles and restricted mobility of both juveniles and adults [START_REF] Vasconcelos | Spawning of Hexaplex (Trunculariopsis) trunculus (Gastropoda: Muricidae) in the laboratory: description of spawning[END_REF] strongly limit dispersal distances. Several studies of species without a dispersing larva (e.g. [START_REF] Boissin | Defining reproductively isolated units in a cryptic and syntopic species complex using mitochondrial and nuclear markers: the brooding brittle star, Amphipholis squamata (Ophiuroidea)[END_REF] reported genetic structure at a small geographical scale. Studies comparing pairs of closely related species with contrasting larval dispersal [START_REF] Tarnowska | Comparative phylogeography of two sister (congeneric) species of cardiid bivalve: strong influence of habitat, life history and post-glacial history[END_REF][START_REF] Weber | Influence of the larval phase on connectivity: strong differences in the genetic structure of brooders and broadcasters in the Ophioderma longicauda species complex[END_REF], as well as reviews of the literature (e.g. [START_REF] Selkoe | Marine connectivity: a new look at pelagic larval duration and genetic metrics of dispersal[END_REF], confirmed the influence of the initial life stage for benthic invertebrates.

a major phylogeographical break acroSS the StS or acroSS the peloponneSe arc A distinction between the Aegean and other populations was reported in other marine species and was associated with the phylogeographical barrier located at the straits of the Peloponnese arc. The so-called Peloponnese break opposes the western, Adriatic and central basins (including the Ionian Sea) to the Aegean and Levantine basins [START_REF] Bahri-Sfar | Fragmentation of sea bass populations in the western and eastern Mediterranean as revealed by microsatellite polymorphism[END_REF][START_REF] Nikula | Phylogeography of Cerastoderma glaucum (Bivalvia: Cardiidae) across Europe: a major break in the Eastern Mediterranean[END_REF][START_REF] Marzouk | Genetic variability of the banded Murex (Hexaplex trunculus) revealed by ND2 and ITS2 sequences[END_REF]). An Aegean-Black Sea-Caspian Sea mitochondrial clade was also found in the bivalve Cerastoderma glaucum, which lives in a fragmented environment, the larvae of which have a short planktonic duration [START_REF] Tarnowska | Complex genetic population structure of the bivalve Cerastoderma glaucum in a highly fragmented lagoon habitat[END_REF]. In a brittle-star species complex, this region also seems to represent a distribution limit for brooders [START_REF] Weber | Genetic data, reproduction season and reproductive strategy data support the existence of biological species in Ophioderma longicauda[END_REF][START_REF] Weber | Influence of the larval phase on connectivity: strong differences in the genetic structure of brooders and broadcasters in the Ophioderma longicauda species complex[END_REF]. By contrast, and unlike numerous species (e.g. [START_REF] Borsa | Genetic structure of the flounders Platichthys flesus and P. Stellatus at different geographic scales[END_REF][START_REF] Patarnello | Pillars of Hercules: is the Atlantic-Mediterranean transition a phylogeographical break?[END_REF], there is no deep genetic differentiation between Atlantic and Mediterranean populations in H. trunculus, and several haplotypes are shared across the strait of Gibraltar and across the Almeria-Oran front. Phylogeographical breaks between the western and the eastern (including central) Mediterranean basins (located at the Siculo-Tunisian sill), and isolating the Adriatic Sea, have been reported for many marine invertebrates (e.g. [START_REF] Borsa | Genetic structure of the flounders Platichthys flesus and P. Stellatus at different geographic scales[END_REF][START_REF] Pérez-Losada | Testing hypotheses of population structuring in the Northeast Atlantic Ocean and Mediterranean Sea using the common cuttlefish Sepia officinalis[END_REF][START_REF] Zitari-Chatti | Genetic variation and population structure of the caramote prawn Penaeus kerathurus (Forskäl) from the eastern and western Mediterranean coasts in Tunisia[END_REF], 2009;[START_REF] Fadhlaoui-Zid | Genetic structure of Octopus vulgaris (Cephalopoda, Octopodidae) in the central Mediterranean Sea inferred from the mitochondrial COIII gene[END_REF][START_REF] Villamor | Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates[END_REF]Deli et al., 2015a;Deli, Said & Chatti, 2015b). In H. trunculus, the Adriatic Sea was not associated with any deep divergence. The parapatric distributions of the two main lineages strongly suggest that they arose by vicariance, i.e. their divergence is due to long-term separation without any genetic exchanges of their distribution ranges. Since lineage II also appears in southern Italy and in the STS, two scenarios appear equally likely.

(1) Vicariance initially separated Aegean populations (possibly, together with Levantine populations) from all the other Mediterranean populations, subsequently followed by migration from the Aegean Sea to the west following the main currents, along the northern shores of the Mediterranean (the general Mediterranean coastal circulation is clockwise). If only Aegean populations were initially isolated by vicariance, they may subsequently have expanded towards the Levantine basin and westward towards southern Italy. (2) Vicariance separated two geographical areas across the STS, and was followed by colonization following the main clockwise Mediterranean currents, explaining the presence of lineage I along eastern Tunisian coasts and lineage II in southern Italy. Morphological data [START_REF] Marzouk | Morphometric variation of fishery-exploited Muricidae (Hexaplex trunculus) in the Mediterranean Sea and the Northeastern Atlantic Ocean: Implications for stock identification[END_REF], reflecting the average of many genes, and possibly the influence of environment, grouped eastern Tunisia and Adriatic with Aegean populations, and STS with western Mediterranean populations. They thus support the STS vicariance rather than the Peloponnese vicariance hypothesis. However, the presence of lineage I in the northern Adriatic Sea is not easily explained by an initial STS vicariance. This region was above sea level until very recently (after the last glacial maximum). Humans, large fishes or birds may have transported lineage I specimens to the Adriatic Sea [START_REF] Green | Recent advances in the study of long-distance dispersal of aquatic invertebrates via birds[END_REF].

Even though vicariance led to cryptic speciation in H. trunculus, hybridizations are now possible, as shown by recombined genotypes found in STS populations. Our data do not establish whether or not lineages I and II evolved genomic incompatibilities, and thus whether they can be considered as distinct species or subspecies. The non-significant exact tests (for association between lineages of different markers) could be the consequence of small sample sizes or of the analysis of genomic regions unlinked to genomic incompatibilities.

hiStorical iSolation between the two parapatric lineageS, poSSibly cryptic SpecieS

For the two studied markers (cox1 and i29), the major split was estimated to have occurred about 5 Mya, which coincides with the Pliocene. This divergence time may be an overestimate because it was based on an older calibration point (40 Mya) [START_REF] Ho | Time dependency of molecular rate estimates and systematic overestimation of recent divergence times[END_REF]. Using the mutation rate from gastropods separated by the Panama Isthmus (3.3 Mya) (phylogenetically distant from Hexaplex), led to a much more recent divergence, of 1 Mya. In [START_REF] Ho | Time dependency of molecular rate estimates and systematic overestimation of recent divergence times[END_REF], the drastic changes of inferred rates with calibration times occurred between 0 and 2 Mya, and calibration points older than 2 Mya for mitochondrial DNA led to similar divergence rates. This suggests that time dependency of divergence rate is not responsible for the differences we found because the two calibration points leading to the conflicting divergence time estimations were older than 2 Mya. We thus favour a change in evolutionary rates among gastropods, and for this reason, we trust the times estimated from Muricinae fossils more than that of Panama one. The similar dates of divergence of the two lineages for the two independent markers lend additional support to the hypothesis that this major split corresponds to vicariance and not to selective sweeps (which can produce similar haplotype divergence patterns), because such sweeps are unlikely to occur simultaneously for distinct markers. This divergence estimate of ~5 Mya is compatible with the end of the Messinian salinity crisis, during which the Mediterranean basin was split into several sub-basins with varying salinity, and may have led to allopatric divergence (eventually speciation) [START_REF] Sabelli | The making of the Mediterranean Molluscan biodiversity[END_REF][START_REF] Calvo | Surviving the Messinian Salinity Crisis? Divergence patterns in the genus Dendropoma (Gastropoda: Vermetidae) in the Mediterranean Sea[END_REF]. Lineages with distinct or allopatric geographic distributions were reported in other Mediterranean gastropods with a short larval phase, as in the muricid species complex Ocinebrina edwardsii [START_REF] Barco | Molecular data reveal cryptic lineages within the northeastern Atlantic and Mediterranean small mussel drills of the Ocinebrina edwardsii complex (Mollusca: Gastropoda: Muricidae)[END_REF] and the vermetid Dendropoma petraeum [START_REF] Calvo | Hidden Mediterranean biodiversity: molecular evidence for a cryptic species complex within the reef building vermetid gastropod Dendropoma petraeum (Mollusca: Caenogastropoda)[END_REF][START_REF] Calvo | Surviving the Messinian Salinity Crisis? Divergence patterns in the genus Dendropoma (Gastropoda: Vermetidae) in the Mediterranean Sea[END_REF]. As for H. trunculus, their distribution also spanned the eastern Atlantic coasts and their history was probably similar. The distribution of individual mitochondrial lineages in these taxa is even more limited than in H. trunculus. Since the same marker was used (cox1) in all three taxa, their contrasted range sizes may have been caused by contrasted dispersal ability, habitat fragmentation or human introductions in H. trunculus. Since the banded murex has no larval stage, less fragmented habitat (e.g. greater depth range) and human introductions are the most likely explanations (see below).

We could not establish whether the genomes of the two lineages were perfectly compatible or not, that is whether hybrids resulting from crosses between the two lineages had lower fitness due to genomic incompatibilities, owing to restricted sample sizes in the sympatry area and to a limited number of loci. The location of this contact zone, concordant with other marine species' transition zones [START_REF] Borsa | Genetic structure of the flounders Platichthys flesus and P. Stellatus at different geographic scales[END_REF], could be explained by a coupling between genetic incompatibilities (reduced or null fitness of the 'hybrids' or introgressed individuals), and exogenous limits to gene flow (differential survival of the two lineages in the environments separated by the contact zone or hydrodynamic limits) [START_REF] Barton | Dynamics of hybrid zones[END_REF]. It would be of interest to analyse more markers in order to assess the extent of genomic incompatibilities between the lineages.

In conclusion, in most of its distribution area, H. trunculus populations are composed of a single lineage, but there may be two cryptic species, currently hybridizing in the STS and possibly along southern Italian coasts. The absence of the whole western haplogroup I-1 i29 in the STS suggests the historical isolation of these populations, favouring the ancient hybridization event hypothesis (current absence of parental species and F1 individuals). Since this typical haplogroup is not recovered with cox1, despite an earlier coalescence time of the mitochondrial genome, it is also likely that selective events led to the evolution of this i29 haplogroup in the STS. Additional markers and/or samples are needed to test the possible scenarios.

Different Demographic hiStorieS in the two lineageS

The neutrality tests, expansion parameter estimations and haplotype network shapes suggested that either demographic expansion had occurred or selective sweeps had affected the genetic markers. The globally similar shapes of the haplotype networks from the nuclear and the mitochondrial genes support the demographic rather than the selective explanation. In effect, a selective sweep would have occurred in lineage I (or occurred more recently in this lineage than lineage II) for both markers cox1 and i29, to account for the data, which is a less parsimonious and less likely explanation than demographic expansion [START_REF] Hoareau | Late glacial demographic expansion motivates a clock overhaul for population genetics[END_REF]. Such haplotype network patterns were often related to expansion after periods of small effective population size caused by Pleistocene climatic fluctuations [START_REF] Patarnello | Pillars of Hercules: is the Atlantic-Mediterranean transition a phylogeographical break?[END_REF][START_REF] Crandall | Expansion dating: calibrating molecular clocks in marine species from expansions onto the Sunda Shelf following the Last Glacial Maximum[END_REF]. Using the cox1 mutation rate estimated from Muricinae fossils in the present study, demographic expansion onset times range from c. 74 to 481 kya. Such estimations are affected differently according to the sampling scheme [START_REF] Städler | The impact of sampling schemes on the site frequency spectrum in nonequilibrium subdivided populations[END_REF]: in the case of expansion, the Tajima and Fu estimators are overestimated (leading to non-significant tests, even in case of true expansion). The results of these tests, interpreted as expansion signals, are therefore expected to be conservative. However, the time estimates of the expansion, in addition to their large confidence interval, may be strongly biased. The robust result is that the two lineages display a very distinct distribution of genetic diversity (illustrated by the different shapes of their haplotype networks), but it is difficult to separate the effect of demographic history from the effect of population structure in a quantitative way. The shapes of the haplotype networks and the expansion time estimates strongly suggest that the two lineages underwent distinct histories. During the Pleistocene period, H. trunculus populations may have undergone several size fluctuations. Populations from the Adriatic Sea are estimated to have the most recent demographic expansion (Table 2), in agreement with knowledge on sea levels and past climates that establish that nearly the whole of the Adriatic was emerged during the last glacial maximum, in contrast to other basins. However the expansion onset time estimated is too ancient to be explained directly by the recolonization of the Adriatic basin after its emersion. The low levels of nucleotide diversity of H. trunculus from Venice, Marseille, Alboran and Atlantic populations, as compared to those of the eastern Mediterranean, suggest that populations from the western part of the distribution range experienced severe genetic drift due to founder effects and/or bottlenecks. However, the haplotypes from the west are definitely not derived from eastern lineage haplotypes (e.g. [START_REF] Tarnowska | Complex genetic population structure of the bivalve Cerastoderma glaucum in a highly fragmented lagoon habitat[END_REF][START_REF] Ragionieri | Population genetics, gene flow, and biogeographical boundaries of Carcinus aestuarii (Crustacea: Brachyura: Carcinidae) along the European Mediterranean coast[END_REF]and references therein), as shown by rooted phylogenetic trees for both markers, so one cannot conclude that the eastern basin played the role of a refugium for the distribution of the whole species. A similar contrast between eastern Mediterranean lineages or populations (displaying demographic stability) and those of the western Mediterranean ones (with signs of demographic expansion) has been reported in other marine invertebrates, including gastropods and bivalves [START_REF] Tarnowska | Complex genetic population structure of the bivalve Cerastoderma glaucum in a highly fragmented lagoon habitat[END_REF][START_REF] Calvo | Surviving the Messinian Salinity Crisis? Divergence patterns in the genus Dendropoma (Gastropoda: Vermetidae) in the Mediterranean Sea[END_REF].

fiSherieS management anD archaeological implicationS

The abundance of this species of high commercial value was reported to have decreased in some places, probably due to over-harvesting [START_REF] Vasconcelos | The artisanal fishery for muricid gastropods (banded murex and purple dye murex) in the Ria Formosa lagoon (Algarve coast -southern Portugal)[END_REF]. Our results suggest that H. trunculus may be considered as at least two main distinct management units, and the strong genetic differentiation within each basin may imply a need for an even smaller management scale. This approach is rarely followed, although it would be appropriate for other harvested Mediterranean species (e.g. [START_REF] Ledoux | Genetic survey of shallow populations of the Mediterranean red coral [Corallium rubrum (Linnaeus, 1758)]: new insights into evolutionary processes shaping nuclear diversity and implications for conservation[END_REF][START_REF] Penant | Planktonic larvae do not ensure gene flow in the edible sea urchin Paracentrotus lividus[END_REF]. Tunisia and southern Italy are the only places in the world where the two lineages were found together in the field. This is a source of strength, because there is more genetic diversity available for breeding and for adaptation to disturbances, but also a weak point, because (1) each lineage is at its range limit and (2) there are at least two stocks which probably cannot replenish one another. Our findings are potentially important for archaeology because this species provided a highly valuable dye, comparable in price with gold and precious metals, from Antiquity (more specifically since the time of the Minoan civilization) [START_REF] Guckelsberger | Purple murex dye in antiquity[END_REF]. The strong population structure and parapatric distribution of the two lineages suggest that living specimens were not translocated by humans or that if they were, they did not settle or reproduce efficiently, between the eastern and western Mediterranean (with the possible exception of Adriatic populations (see above). Interestingly, [START_REF] Merle | First record of Trunculariopsis trunculus (Linnaeus, 1758) (Mollusca: Gastropoda) from the basin of Arcachon (Gironde): the second exotic muricid on the French Atlantic coasts[END_REF] reported the presence of exotic H. trunculus specimens in the basin of Arcachon, more than 1100 km from their natural range [START_REF] Houart | A review of the recent Mediterranean and Northeastern Atlantic species of Muricidae[END_REF], and after more than 2 Myr of absence in the regional fossil record. Determining the area of origin of this exotic population should be possible with our two molecular markers, given the highly significant genetic differentiation they reveal between population pairs. (CNR-ISMAR, Venice-Italy), Andrés Martínez-Lage (Universidad da Coruña, Spain), Zeid MARZOUK, Najeh BOUIDA and Marwene BOUIDA (Tunisia) for collaboration in the collection of samples. We acknowledge indirect support from the SCBM (Service Commun Biologie Moléculaire) of the IMBE. Didier Merle is thanked for his explanations of the complex fossil record situation. Abigail Cahill and Michael Paul are thanked for improving the final version of the manuscript and correcting the English. We thank the anonymous reviewers for their helpful comments. This work was supported by the Ministry of higher Education and

Figure 1 .

 1 Figure 1. Map showing the sampling sites of Hexaplex trunculus. Numbers (1→19) are populations sampled in the present study (see Table1for details), Letters (A→F) are the six populations studied by[START_REF] Villamor | Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates[END_REF]. Pie charts show the distribution of the major haplogroups for cox1 (A 1 : present study, A 2 :[START_REF] Villamor | Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates[END_REF] and for i29 (B). For the cox1 sequences of[START_REF] Villamor | Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates[END_REF], it was not possible to distinguish haplogroups Ia and Ib due to shorter sequence length. ATO, Atlantic Ocean; ALS, Alboran Sea; WMB, Western Mediterranean Basin; STS, Siculo-Tunisian Strait; EMB, Eastern Mediterranean Basin; ADS, Adriatic Sea; AGS, Aegean Sea.

Figure 2 .

 2 Figure 2. Haplotype networks for cox1 (A) and i29 (B). Each circle represents an haplotype. The size of the circle is proportional to the number of individuals displaying each haplotype. Numbers of mutations (when more than one mutation separates the sequences) between haplotypes are indicated near branches. Colours indicate the origin of the individuals displaying each haplotype, as in the legend: to limit the use of colours, spatially close populations were combined according to the regions defined in Table1, except for the four populations from the WMB region because (1) they may be very distant and (2) two very close localities (northern Tunisia) do not fall into the same haplogroup (for cox1). Black circles represent unsampled haplotypes.

Figure 3 .

 3 Figure 3. Bayesian reconstruction of phylogenetic relationships among cox1 haplotypes. Lineage I is symbolized by circles; Lineage II by squares, followed by numbers that indicate individual haplotypes. Yellow squares indicate haplotypes occurring in Siculo-Tunisian strait; red squares indicates occurrence of haplotype in both Aegean Sea and Siculo-Tunisian strait. Numbers above the branches show support values (Bayesian posterior probability for BI and bootstraps greater than 50% for maximum likelihood reconstructions, respectively). Numbers below the branches show the TMRCA (Mya) with their confidence intervals.

Figure 4 .

 4 Figure 4. Bayesian reconstruction of phylogenetic relationships among i29 haplotypes. Lineage I is symbolized by circles, blue circles indicate haplogroup I-1 and green circles indicate haplogroup I-2. Lineage II is symbolized by squares; red squares indicate haplotypes occurring in both Aegean Sea and Siculo-Tunisian strait. Numbers above the branches show support values (Bayesian posterior probability for BI and bootstraps greater than 50% for maximum likelihood reconstructions, respectively). -Bootstraps < 50%. Numbers below the branches show the TMRCA (Mya) with their confidence intervals.

Table 1 .

 1 Diversity measures for the populations of Hexaplex trunculus studied

	cox1
	Code
	Geographical origin

Table 2 .

 2 Neutrality tests, demographic parameter (τ) and expansion times (Mya) with their confidence intervals for cox1 marker

	Geographical origin	n	D	Fs	R 2	τ	ET
	Algarve	24	-0.89	0.92	0.08*	2.56 (0.359-3.50)	0.324 (0.045-0.443)
	Marseille	26	-1.10	-2.25**	0.08*	2.07 (0.346-4.09)	0.262 (0.043-0.517)
	Corse	24	0.24	-1.72**	0.14	3.23 (0.057-6.90)	0.408 (0.007-0.873)
	Tabarka	24	-1.85*	-2.85**	0.09*	0.00 (0.00-0.00)	Constant
	Oropos Attiki	20	-0.04	-7.25***	0.12	2.99 (1.89-4.09)	0.378 (0.239-0.517)
	ATO	35	-1.75*	-2.86**	0.08	2.39 (0.40-3.50)	0.302 (0.050-0.443)
	ALS	44	0.43	-2.02*	0.13	1.43 (0.63-2.49)	0.181 (0.079-0.315)
	WMB	92	-1.68	-15.50***	0.04	1.70 (0.95-2.63)	0.215 (0.120-0.332)
	STS	58	-0.31	-6.06*	0.11	14.5 (4.5-17.46)	1.835 (0.569-2.210)
	EMB	57	-1.05	-9.54***	0.06	2.40 (0.89-5.32)	0.303 (0.112-0.673)
	ADS	10	-0.69	-0.59*	0.17	0.70 (0.00-1.56)	0.088 (0.000-0.197)
	AGS	47	0.13	-15.53***	0.11	3.73 (1.82-5.13)	0.474 (0.057-0.162)
	Lineage I	283	-1.84**	-64.79***	0.03*	0.59 (0.00-0.97)	0.074 (0.000-0.122)
	Lineage II	60	0.21	-14.76***	0.08	3.80 (0.94-5.82)	0.481 (0.118-0.736)
	Geographical region (as in Fig. 1), sample size (n), Tajima's neutrality test (D), Fu's test of neutrality (Fs), Ramos-Onsins and Rozas statistic (R 2 ), ET expansion times using the cox1 mutation rate estimated from Muricinae fossils in the present study (0.6% per Mya).
	*P < 0.05, **P < 0.01, ***P < 0.001.

Z. MARZOUK ET AL. © 2017 The Linnean Society of London, Biological Journal of the Linnean Society, 2017, XX, 1-18

Z. MARZOUK ET AL. © 2017 The Linnean Society of London, Biological Journal of the Linnean Society, 2017, XX, 1-18

© 2017 The Linnean Society of London, Biological Journal of the Linnean Society, 2017, XX, 1-18

0.969 0.0092 Sample size (n), haplotypes and their frequencies, number of haplotypes (Nh), haplotype diversity (Hd) and nucleotide diversity (π).6 Z. MARZOUK ET AL.© 2017 The Linnean Society of London, Biological Journal of the Linnean Society, 2017, XX, 1-18

Z. MARZOUK ET AL. © 2017 The Linnean Society of London, Biological Journal of the Linnean Society, 2017, XX, 1-18

Z. MARZOUK ET AL. © 2017 The Linnean Society of London, Biological Journal of the Linnean Society, 2017, XX, 1-18

Scientific Research of the Republic of Tunisia.

ACKNOWLEDGEMENTS

The authors express their gratitude to Benlahcen RIJAL LEBLAD (INRH, Tanger-Morocco), Dimitris POURSANIDIS and Dr. Drosos KOUTSOUBAS (HCMR, Crete-Greece), Fréderic Zuberer (Pythéas Institute, Marseille-France), Davide TAGLIAPIETRA

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article at the publisher's web-site: Appendix S1. Haplotype networks for cox1 including sequences from [START_REF] Villamor | Genetic structuring across marine biogeographic boundaries in rocky shore invertebrates[END_REF] and the BOLDsystem database [START_REF] Ratnasingham | A DNA-based registry for all animal species: the barcode index number (BIN) system[END_REF]. Appendix S2. F ST (below the diagonal) and Φ ST (above the diagonal) between different populations. Not significantly different values after the multiple correction tests are written in bold. Appendix S3. Non-parametric multidimensional scaling analyses (MDS) based on Φ ST distances between populations for each marker. Appendix S4. Contingency tables comparing main splits for cox1 and i29 markers, in STS populations. Appendix S5. Pairwise mismatch distribution for cox1 marker. Appendix S6. Pairwise mismatch distribution for i29 marker. Appendix S7. GenBank accession numbers of all haplotypes for both markers.

View publication stats View publication stats