E Perrier

B Champion

N Diot

Yohan Vuillerme
email: nicolas.vuillerme@agim.eu

Payan

M Bucki

V Luboz

A Perrier

E Champion
email: eric.champion@texisense.com

B Diot
email: b.diot@ids-assistance.com

N Vuillerme

Y Payan

Clinical workflow for personalized foot pressure ulcer prevention

Keywords: Foot pressure ulcer, Soft tissues, Patient-Specific, Finite element method

Clinical workflow for personalized foot pressure ulcer prevention

 [START_REF] Shaw | The pathogenesis of diabetic foot problems. An overview[END_REF]. In addition to causing pain and morbidity, foot lesions in diabetic patients have substantial direct and indirect economic consequences [START_REF] Shearer | The influence of foot posture, support stiffness, heel pad loading and tissue mechanical properties on biomechanical factors associated with a risk of heel ulceration[END_REF] [START_REF] Gordois | The healthcare costs of diabetic peripheral neuropathy in the UK[END_REF]. Diabetic foot ulcers result from multiple pathophysiological mechanisms, including neuropathy, peripheral vascular disease, high foot pressures, foot deformity, and diabetes severity [START_REF] Telfer | What Has Finite Element Analysis Taught Us about Diabetic Foot Disease and Its Management ? A Systematic Review[END_REF]. Several studies [START_REF] Mueller | Etiology, evaluation, and treatment of the neuropathic foot[END_REF] [START_REF] Loerakker | The effects of deformation, ischaemia and reperfusion on the development of muscle damage during prolonged loading[END_REF] recognized at least three mechanisms leading to pressure ulcer: (1) ischemia caused by increased pressure duration even for low induced strains, (2) high internal tissue strains created by increased pressure magnitude, and/or (3) tissue fatigue caused by increased number of periodic pressure loads. Time and strain have an inversely proportional contribution to ulceration [START_REF] Kosiak | Etiology and pathology of ischemic ulcers[END_REF] [START_REF] Loerakker | The effects of deformation, ischaemia and reperfusion on the development of muscle damage during prolonged loading[END_REF] (Van Schie et al., 2006): high strains take a relatively short time (a few minutes) to cause ulceration whereas low strains induce lesion after a longer period (between two and four hours). Short and long term lesion inducing strain thresholds have been characterized by [START_REF] Loerakker | The effects of deformation, ischaemia and reperfusion on the development of muscle damage during prolonged loading[END_REF] in muscle tissues. The obtained values were around 50 % of deformation for short term high strains and 20 % of deformation for long term low strains. This study also showed that fat tissues have large strain variations (although not as large as muscle tissues) and they might suffer from pressure ulcer. The two strain thresholds aforementioned are therefore key values in pressure ulcer prevention.

Daily monitoring by the patient or clinical staff is the main tool to prevent foot pressure ulcers and results in an estimated reduction of foot ulcers and amputations from 50% to 80% (Boulton et al., 2005). Because early stages of ulceration are not always visible, both patient's and staff's vigilances tend to decrease over time. Unfortunately, in the case of diabetic patients, it is precisely when the first ulcers appear that serious complications develop, mainly because of the angiopathy, which severely limits healing.

It is consequently essential to introduce new monitoring tools to promote awareness and as a result, patient's autonomy in everyday life. Measuring pressure loads at the skin surface, all around the foot, and, if possible, estimating the corresponding internal strains could help preventing further ulceration and facilitate wound healing [START_REF] Gefen | The biomechanics of heel ulcers[END_REF]. Measuring interface pressures can be performed with a pressure sensor such as the ones proposed by Novel (http://www.novel.de), Tekscan (http://www.tekscan.com), Vista Medical (http://www.pressuremapping.com) or Texisense (http://www.texisense.com), however several studies established that using pressure measurements at the skin interface is not sufficient to prevent foot pressure ulcers, especially the ones starting deep in the tissues and causing substantial subcutaneous damage underneath intact skin (Linder-Ganz et al., 2008) (Atlas et al., 2009) [START_REF] Gefen | Plantar soft tissue loading under the medial metatarsals in the standing diabetic foot[END_REF]. Indeed surface measurements do not provide enough information as to predict ulcer formation in a reliable way [START_REF] Linder-Ganz | Strains and stresses in sub-dermal tissues of the buttocks are greater in paraplegics than in healthy during sitting[END_REF]. For example, with the same pressure map, a patient with a sharp calcaneus, or a thinner heel pad, could develop a pressure ulcer while another one, with a different morphology, might not. Pressure ulcer risk is consequently highly patient-related and integrates a number of factors such as bones' curvature [START_REF] Luboz | Influence of the Calcaneus Shape on the Risk of Posterior Heel Ulcer Using 3D Patient-Specific Biomechanical Modeling[END_REF], or soft tissue thickness (in skin, fat and muscles) [START_REF] Gefen | The biomechanics of heel ulcers[END_REF]. Monitoring internal strains is currently a consensual criterion to assess the risk of pressure ulcer and has been widely used in previous studies (Linder-Ganz et al., 2008) [START_REF] Oomens | Can Loaded Interface Characteristics Influence Strain Distributions in Muscle Adjacent to Bony Prominences[END_REF]. Nevertheless, measuring internal strains in vivo being impossible, a biomechanical model integrating the behavior of the foot internal soft tissues and bones is needed to assess internal strains and the resulting risk of ulceration. Furthermore because of inter-individual anatomical variability, personalized biomechanical models must be resorted to in order to accurately estimate internal strains and implement an adequate prevention strategy [START_REF] Gefen | The biomechanics of heel ulcers[END_REF] [START_REF] Luboz | Influence of the Calcaneus Shape on the Risk of Posterior Heel Ulcer Using 3D Patient-Specific Biomechanical Modeling[END_REF] [START_REF] Tenenbaum | Effects of foot posture and heel padding devices on soft tissue deformations under the heel in supine position in males: MRI studies[END_REF].

Several studies have demonstrated the use of biomechanical foot models to estimate internal strains. Most of the feasibility studies are limited to a single foot model generated for a specific patient, and seem difficult to extend in an automatic fashion to a wider group of subjects -not to mention -in clinical routine. For example, [START_REF] Ledoux | A quasi-linear, viscoelastic, structural model of the plantar soft tissue with frequency sensitive damping properties[END_REF]) modelled the soft tissues under the foot (skin, fat and muscles) as a Finite Element (FE) mesh with a homogeneous linear elastic material, bones as rigid FE meshes; joints were accounted for as idealized contacts between bones, and around 20 ligaments connecting the mid foot bones were modelled as cables. In another study, [START_REF] Chen | Effects of internal stress concentrations in plantar soft-tissue -preliminary three-dimensional finite element analysis[END_REF] proposed a more detailed FE foot model including almost all foot ligaments and using a large deformation Mooney Rivlin constitutive law for the soft tissue bulk. Even though this model is fairly complete, it lacks computational efficiency and does not distinguish between different tissue types. These drawbacks were addressed in the model that we recently proposed [START_REF] Perrier | Conception and evaluation of a 3D musculoskeletal finite element foot model[END_REF] with foot soft tissues represented as four different Neo Hookean materials for skin, fat, heel pad fat, and muscles respectively. In this model, bones were represented as rigid bodies connected by the most significant ligaments of the foot, modeled as cables. Nevertheless, this last model, just like the two previously cited ones, was generated from a single subject dataset and is consequently only representative of this particular morphology. In this paper, inspired by our previous study on patientspecific modelling of the calcaneus [START_REF] Luboz | Influence of the Calcaneus Shape on the Risk of Posterior Heel Ulcer Using 3D Patient-Specific Biomechanical Modeling[END_REF], we propose to use this complex foot model as an atlas -or generic model -and to generate new patient-specific models by deforming this atlas to fit the patients' specific morphology. The goal is to design a process making it possible to produce patient-specific biomechanical models in the most automated and user-friendly way possible. The proposed modeling technique could be used to study the influence of variability in morphology on pressure ulcer formation. Its further goal is to provide insight at how morphological specificities should be accounted for in the design of medical devices to optimize strain monitoring-based prevention for each individual. The following study has been carried out in a static analysis framework i.e. does not take into account the duration or repetitive mechanisms leading to pressure ulcer but only tissue compression resulting from a static stance.

Methods

Foot model atlas

The shape of the atlas model is based on a single subject (male, 33, healthy) and is presented in details in (Perrier et al. 2014a) [START_REF] Perrier | Conception and evaluation of a 3D musculoskeletal finite element foot model[END_REF]. The contours of the skin, heel fat pad, muscles, and bones were manually segmented from the CT scan of the right foot of this healthy subject. An automatic FE mesh generator (developed by Texisense) was run on the resulting surfaces, and produced a conforming multi-domain FE mesh containing four layers: muscles, fat, heel fat pad, and skin (Figure 1). The meshing algorithm generates as many hexahedrons as possible in the core of the continuum to limit the locking effect observed for tetrahedral elements under quasi-incompressible assumption. Smooth and conforming boundaries between the different internal domains are defined using transition elements such as pyramids, wedges, or tetrahedrons. The meshing procedure led to a FE mesh having 44,220 elements, including 3,610 hexahedrons, 12,062 pyramids, 8,674 wedges, and 19,874 tetrahedrons, for a total of 19,574 nodes.

Finite Element analyses are carried out on the 3D simulation platform ArtiSynth (Lloyd et al., 2012)(www.artisynth.org). Soft tissues (skin, fat, and muscles) are modelled using Neo Hookean materials in order to account for large deformations. Each layer is assigned distinct material properties, drawn from literature (Sopher et al., 2011). Young moduli were set to 200 kPa for the skin, 60 kPa for the muscles, 100 kPa for the heel fat pad, and 30 kPa for the rest of the fat. Assuming these tissues are quasi-incompressible, we set the Poisson ratios to 0.485 for the skin, 0.495 for the muscles, 0.499 for the heel fat pad, and 0.49 for other fatty tissue. Bones, featuring a significantly higher stiffness, are modeled as rigid bodies and their shapes are cut out within the soft tissue continuum, i.e. without any finite element inside the volume. The foot 28 bones, the tibia, and fibula are integrated in the model. Each bone can collide with its neighbors and is connected to them by several ligaments, forming the joints. FE nodes nearby bony surfaces are automatically attached to the neighboring bones, which results in a non-sliding soft tissue-bone interface. The main musculoskeletal structures are modeled by active cables elements within the Artisynth framework.

The numerical foot model is divided in two components. The first component is a musculoskeletal model accounting for rigid body motion within the foot. This model implements anatomical constrains such as contacts between adjacent bones, action of the ligaments, or simulated muscle contractions. The second component is the soft tissue continuum modeled by the FE mesh which is iteratively coupled with the musculoskeletal component and translates boundary conditions and internal rigid body motion into elastic deformation in soft tissues.

Patient data

In order to generate patient-specific models, a description of the patient's morphology is needed and more specifically the contours of bones as well as the external skin surface. These shapes are used by the atlas-to-patient registration procedure in order to compute the anatomical transfer of musculoskeletal data. This three-dimensional information can be provided by various sources such as CT, MR scans, or EOS images (provided a 3D reconstruction of the morphology is performed by the latter bi-planar imaging device). Extracting tissue contours from these modalities involves specific automatic or semi-automatic procedures and is a challenge in itself that lies beyond the scope of this article. In our study, this task has been carried out manually. In the remainder of the article, we describe our workflow assuming that the medical images are already labelled (i.e. segmented), forming so-called "binary images". Each label in a binary image represents a distinct bone or soft tissue, which will be implemented as a modeling domain in subsequent mesh generation steps.

Three patients were included in our study: a 70 year old male (BR), a 67 year old male (FP) and a 55 year old male (FC). Images were acquired while the patients were in dorsal decubitus, legs and feet supported by the table. The exams were performed in the context of a vascular exam using CT angiography. Clinical exam revealed that the feet in all three patients were healthy. The size of the CT volumes is as follows: 220x404x519 for BR, 196x276x575 for FC, and 183x297x580 for FP, and the resolution is 1x1x1 mm 3 in all three datasets. The right foot is modeled for both FC and FP, and the left foot for BR.

Patient-specific model generation

Shapes of patient's bones and skin were recovered from binary images. Bony surfaces were directly used as rigid bodies in the numerical model. Each patient's personalized musculoskeletal model was generated by registering the patient's bones and skin with their counterparts in the atlas. The Mesh-Match-and-Repair (MMRep) registration algorithm is used to perform this task automatically in about a minute (Bucki et al., 2010). Our implementation of the atlas to patient registration procedure is divided into three steps increasingly introducing distortion in the data: (1) a rigid registration that roughly positions the patient data set with respect to the atlas model, (2) an affine deformation that compensates for global scale discrepancies, followed by 3) an elastic registration that accurately fits the bony contours and the skin surface. Once all three deformation functions are combined, the resulting deformation is applied to the atlas dataset to transfer the atlas information (muscles, ligament insertions, fat pad, plantar fascia) into the patient's referential. The procedure -producing the musculoskeletal component of the patient model -is automatic and takes about two minutes.

The assembly of the soft tissue continuum mainly consists of the generation of the FE mesh corresponding to the fat, muscle, fat pad and skin domains. The outlines of these domains however are not present in the binary images in our dataset. Indeed, for the sake of integration of our procedure within a realistic clinical workflow, we assumed that only a basic segmentation could be performed to extract the prominent morphological features that are the bones and skin. In order to overcome this practical hurdle, we again resort to the atlas approach to infer the missing information in the patient data. The volumetric deformation function computed using the MMRep algorithm for the tendon and ligament insertions is applied to the soft tissue domains defined in the atlas that we wish to replicate in the patient's model. The outlines (materialized by triangular surface meshes) of both the muscles and fat pad are deformed and their position is adjusted to fit the patient's bones and skin. Then Texisense mesher algorithm is used to automatically produce a conforming FE mesh of the domains. This step takes about three minutes. Before the mesh generation, a preprocessing step for the selection of the region of interest containing the foot (simple cropping of the image around the skin and bones) is performed manually and requires about 10 minutes of user intervention. The definitions of the muscle and fat pad subdomains in the patient model rely on assumptions formulated in the atlas as well as approximations involved by the registration procedure. However, in the current state of the art, we believe that the atlas paradigm provides the right balance between the efficiency of existing image processing techniques and a level of accuracy required for the targeted biomechanical simulation in the context of pressure ulcer prevention. In the future, should a new segmentation algorithm (or a new imaging modality) appear that would enable an accurate and cost effective segmentation of one or all of these domains, our approach could easily take advantage of it by replacing our registered domain by the actual one yielded by the novel technique.

Once the FE mesh has been produced, a quality control and optimization step is performed to improve the elements that might put FE analysis in jeopardy. These elements are identified using popular FE quality metrics such as Aspect Ratio and Jacobian Ratio and derivations of such. This step is semi-automatic as an informed user needs to supervise this mesh untangling. However, once the parameters set, mesh relaxation is automatic and takes approximately five minutes.

Lastly, since the atlas is a right foot, plain mid-sagittal plane symmetry flips it to accommodate a left foot in the patient.

The whole above described model specialization process takes less than 20 minutes.

Estimation of foot ulcer risk through simulation

To study the influence of foot morphology on the location and magnitude of internal strains, and therefore on the risk of pressure ulcer, a common pressure pattern simulating a static unipodal stance was applied below all three virtual foot soles. The chosen plantar pressure pattern was measured using a commercially available pressure sensor (Zebris platform, http://www.zebris.de/) under the right foot of the subject used to build the atlas (Figure 2). Most of the plantar pressures mainly appear below the heel and the metatarsal heads, with a peak of 14.5 N.cm -2 below the calcaneus. The pressure pattern has been mirrored prior to applying the boundary conditions on the left foot of patient BR. The pressure map was aligned under the foot by fitting the highest pressure peak below the lowest point of the calcaneus for each patient, which implicitly forms the assumption that this bony prominence is the source of peak pressures under the heel. The axis of the foot was given by the vector pointing from this lowest calcaneus point to the lowest point under the second metatarsal head.

During the simulations, the tibia was fixed while the rest of the foot and the fibula were left free to move. The first phase of the simulations allowed the foot to relax under the influence of the multiple tendons, ligaments and muscles of the model which tend to recover their equilibrium length and thus generate pre-stresses in the FE continuum. These musculoskeletal structures were not initially at their resting length because their morphology is derived from the medical image dataset, and consequently from the pre constrained position imposed upon the patient during the medical exam. Once steady state reached, pressure patterns were projected below the foot and the FE nodes at the surface of the foot model were assigned the pressure corresponding to their position in the pressure pattern. These normal pressure values were converted into nodal forces. Plantar pressure was applied gradually. Once 100 % of the pressure was applied, the loading was maintained until equilibrium was reached. The reported strain values are those measured in the final equilibrium configuration of the models.

Cluster analysis and regionalization

In order to assess the risk of pressure ulcer, besides monitoring the maximal internal strains, we introduce a novel paradigm in our Finite Element simulations by considering the volume of the largest group of adjacent elements with nodes exhibiting a VM strain over one of the considered 20 or 50 % thresholds, as suggested by [START_REF] Loerakker | The effects of deformation, ischaemia and reperfusion on the development of muscle damage during prolonged loading[END_REF]for the muscle tissue. We call "clusters" the isolated groups of such adjacent and over-strained elements. Clusters are determined by aggregating all the neighboring elements with strains higher than a given threshold. The external boundary of a cluster is defined as the set of all cluster element faces not shared by another element in the cluster. The shape, volume and hence boundary of clusters depend on each individual's morphology and tissue behavior. A cluster can be heterogeneous in tissue nature which means that fat, muscle, and/or skin elements can share a common cluster. In the absence of results for fatty tissues at the time of this study, we applied the same strain thresholds for the whole soft tissue bulk. Should any future experiments lead to strain threshold values in fat or skin, these parameters could easily be integrated in the cluster definition by merely adjusting each element's inclusion test to the threshold of the tissue it stands for in the model.

One possible interpretation that can be derived from clusters is referred to, herein, as "cluster volume" and is merely the volume (in mm 3 ) of the considered cluster. However, other "indicators"i.e. scalar interpretations of these clustered subsets of elements within the continuum mesh -can be contemplated but were not investigated in this study. Let's just mention a few : the extent of a cluster in a specific direction -e.g. in anisotropic tissues, the amount of blood vessels that it encompasses and is likely to influence, the proportion of different kinds of soft tissues within, or its geometrical correlation with the patient's previous lesion history.

Clusters allow quantitative comparison of tissue suffering levels among individuals while overcoming the lack of reliability of peak VM strains which can locally stem from numerical uncertainty in the FE analysis. The mathematical rationale here is that important jumps in the gradient of the solution (here the displacement field) are a known indicator of local numerical uncertainty in Finite Element analysis [START_REF] Verfürth | A review of a posteriori error estimation techniques for elasticity problems[END_REF] and it is unwise to draw conclusions from these local epiphenomena.

Cluster analysis however introduces a new unknown which is the "minimal lethal cluster volume" i.e. the value of the volume of tissue above which a pressure ulcer may develop following one of the above-defined thresholds i.e. short term or long term lesion. In the absence of physiological definition of such volume, we are at the moment restricted to relative conclusions which can be formulated as: "this patient is more at risk than that patient, or this insole is better suited for this patient than that insole." Indeed, cluster volume might either indicate the lack of relevance of a strain value if it is associated with a negligible volume (yet to be defined) or, on the contrary, show that the strain value is observed at a macroscopic scale and most likely affects a significant volume of tissue, hence leading to a possible lesion. To rephrase our proposition, the assertion "The maximal volume of tissue clusters undergoing a VM strain > X % is Y mm 3 " (which can be formulated in the framework of cluster analysis and interpreted in a comparative study) is more intuitive and numerically more robust than the assertion "The maximal VM strain in the model is Z %."

Cluster localization within the foot volume also provides information on the areas where lesions are prone to appear. To make the interpretation of the results more intuitive and clinically relevant the foot was partitioned into eight key anatomical regions defined as follow (Figure 3): (1) the Achilles tendon, (2) the top of the foot, (3) the heel, (4) the medial foot, (5) the first metatarsus, (6) the four other metatarsi, (7) the hallux, ( 8) and the four other toes. This partitioning makes it possible to correlate cluster volumes with their respective locations within the foot anatomy which in turn relates to corresponding foot functions that might be affected by a potential lesion. A risk level per region can thus be assessed and used to refine a prevention strategy or the design of an orthosis. The foot regions are defined within the atlas and are automatically adjusted to each patient's foot using the MMRep inference procedure described above.

Results

Finite Element meshes for the atlas model and the three studied patients are presented in Figure 4. Morphological differences between the three patients are obvious: various external foot shapes, and various individual bone shapes and orientations. The main difference in terms of external shape is around the phalanges of subject FC compared to subjects FP and BR, and even though moderately, with the atlas. The ankles of the atlas and subject FC are also quite different in shape compared to subject FP and BR, which are more prominent. In terms of internal morphology, BR has very narrow metatarsal bones compared to the other subjects. The talus of the atlas and subject FC are wider than the ones of subjects FP and BR. The navicular of the atlas and subject FC are located more medially than the ones of subjects FP and BR. The calcaneus of subject FC and BR are curved more medially than the ones of subject FP and the atlas.

The maximal Von Mises strains and cluster volumes for the atlas and patients are summarized in Table 1 (for a threshold of VM strains over 20%, representing long term lesion) and Table 2 (for a threshold of VM strains over 50%, representing short term lesion). For each patient (atlas, FC, FP and BR) the table is split vertically in two columns: maximal cluster volume for the considered VM strain (20% in table 1, 50% in table 2) and maximum VM strains inside that maximum cluster volume. Both tables are divided horizontally providing detailed strain information for each of the eight anatomical regions. The mean and standard deviation are also provided for the eight regions.

For the atlas, it can be seen that the highest cluster volume is located in the heel region (50.8 cm 3 with VM strains above 20 %). As for the VM strains, the highest values are located in region 6 comprising the second, third, fourth and fifth metatarsi (161 %), the top of the foot region (145 %), and the medial foot region (107 %). These values are highlighted in Figure 5, which illustrates, in a color code, the VM strains and the cluster volumes region by region. For each foot, the maximum value is coded as red and the minimum is blue. In Figure 5, the Achilles and top of the foot regions are not shown as the regions mostly at risk under plantar compression are below the foot. Figure 6 shows the differences in locations of the maximal VM strain and the largest cluster volume, compared to the calcaneus bone.

For patient FC, the maximum VM strain (348 %) is located in region 5 with the first metatarsal bone, while the maximum cluster volume (74.4 cm 3 for VM strains above 20 %) is in the heel region. For patient FP, these regions become the top of the foot region with values of 399 % (which is a numerical singularity due to a node compressed between two bones: the calcaneus and the cuboid) for the maximum VM strain and the medial foot region with a maximum cluster volume of 62.2 cm 3 (for VM strains above 20 %). Note that for this patient, a second region is highly at risk: the heel region, as it reaches a cluster volume of 57.5 cm 3 (for VM strains above 20 %). For patient BR, the maximum VM strains point to the medial foot region with 205 % for a high risk of pressure ulcer while the cluster volumes again highlight the heel region (63.0 cm 3 for VM strains above 20 %) as most at risk. Note that the cluster volumes for the medial region is only 14.1 cm 3 , therefore corroborating that this region is less at risk than the heel region.

A similar analysis can be made for all datasets using Table 2 and Figure 5 based on the VM strain threshold of 50 %.

Discussion

For the same pressure pattern applied below the foot of each patient, the variability of the results among patients presented in Tables 1, 2 and in Figure 5 clearly point out the influence of the patient's morphology. Whereas the simulations report a huge range of variations for maximum VM strains (between 161 to 399% from one patient to the other, with almost all foot regions affected, including the upper part of the foot which seems not relevant), the values reported for the maximum cluster volumes seem more coherent. Indeed, except a specific case (patient FP who has a medial foot maximal cluster volume barely higher than for the heel), the maximal cluster volume is located in the heel region for all patients, with values ranging between 50.8 and 74.4 cm 3 (for VM strains above 20 %) and between 0 and 21.6 cm 3 (for VM strains above 50%).

These results seem therefore to show that only monitoring the maximum VM strain is probably not the best option to evaluate the risk of pressure ulcer and its location. This observation seems similar to the one stemmed from another study on buttocks' ulcer analysis [START_REF] Luboz | Biomechanical Modeling to Prevent Ischial Pressure Ulcers[END_REF]. The non-realistic strain value of 399 % reported for patient FP is a good example of numerical "outliers" that can be generated by any FE model submitted to high pressures. This high value strain singularity is located at only one node that appears to be squeezed between the cuboid and calcaneus bones; the element associated to that node is therefore deformed a lot because of the mesh configuration and not because of the pressure applied below the foot. Furthermore, the high VM strains reported in the top of the foot (regions 1 and 2) for the atlas, patient FP and patient BR also show the limitations of using this criterion as high strain values should not be located in the top of the foot during a plantar load. Such singularities can be ignored by analyzing the maximum cluster volumes. Indeed, this analysis checks if the maximum strain is located at a single node (or few nodes), by computing the volume associated to contiguous elements exhibiting a strain over a given threshold. Outliers such as the one observed in our simulations will thus automatically be ignored by the cluster analysis since they only affect a much reduced volume of tissues. This is demonstrated especially for the atlas and patient FP where the heel region is not described as at risk by the maximal VM strain analysis while it is clearly at risk for the cluster volume one.

Given the variety of bone shapes observed in our small sample, see Figure 6 for an example on the calcaneus shape, we think that an accurate representation of the internal structures is necessary to capture accurately the behavior of soft tissues under compression for an articulated foot. A non-linear registration step is therefore compulsory in order to transfer the anatomical knowledge from the atlas to the patient. Using only steps 1 and 2 (rigid and affine) in the MMRep procedure described in paragraph 2.3 would result in an approximate deformation that is likely to miss morphological specificities that can possibly result in an injury, see Figure 6 for the location of the strains above the threshold of 20 % below the calcaneus of each patient. This is particularly true in the case of strongly pathological feet exhibiting large differences with the atlas mesh.

Another important conclusion is that, consistently with other pressure ulcer prevention studies [START_REF] Gefen | Plantar soft tissue loading under the medial metatarsals in the standing diabetic foot[END_REF] [START_REF] Gefen | The biomechanics of heel ulcers[END_REF] [START_REF] Luboz | Biomechanical Modeling to Prevent Ischial Pressure Ulcers[END_REF] [START_REF] Luboz | Influence of the Calcaneus Shape on the Risk of Posterior Heel Ulcer Using 3D Patient-Specific Biomechanical Modeling[END_REF], it appears impossible to build a pressure ulcer risk assessment scale by relying solely on interface pressures. The compression of soft tissues, with patient dependent thicknesses and parameters, under personalized bony prominences is key to an efficient personalized prevention strategy. This is for example the case for the medial foot region of patient FP which is subject to the same pressure as the other patients but turns to be at risk for patient FP and not for the others.

A large variability of results has been observed on a small sample (N=3) of patients and these trends would most likely be confirmed on a larger sample.

One limitation of our models is probably the fact that the same generic tissue constitutive law was proposed for each patient. The Neo Hookean law was chosen to simulate the quasiincompressibility of the soft tissues and to account for large deformations. The Young moduli of each tissue (muscle, skin, fat and fat pad) are based on the literature, as introduced in section 2.1. We assumed the same constitutive parameters to avoid hindering the influence of the anatomy in this study. However, the variations of the soft tissues properties would need to be specified for each patient, for example by using indentation or elastography techniques.

Another limitation identified during this study is the difficulty to compute a proper non-linear registration between the atlas and patient toes. Indeed, the differences between the toes' posture combined with the complexity of the shape drive the elastic registration algorithm towards a local minimum. Registration accuracy in this region is thus compromised as can be seen in Figure 5 where the toes in patients exhibit an unnatural distortion (see e.g. the fifth toe of patient FP and the first toe of patient BR). However, this phenomenon has only a regional effect and does not influence the results observed in other regions thanks to the realistic stress smoothing introduced by the foot tendons and ligaments included in our model and the different layers of material properties modelling the skin, fat and muscles. The discussion carried out on the heel and medial foot thus retains its relevance. An accurate registration in the toe region will require the individualization of each toe in order to rewrite the objective function used in MMRep to compute the volumetric deformation.

Finally, it must be noted that our comments are based on a static analysis (unipodal stance) while pressure ulcers can appear while walking. Gait can be decomposed into a number of foot positions on the ground and their corresponding pressure maps. Applying these pressure maps under the foot model would make it possible to simulate the deformation of the foot at each gait phase. It would therefore be possible to perform a dynamic analysis and estimate the risk of pressure ulcer at these stages exactly as it is done in this paper for unipodal standing alone.

In order to implement an efficient pressure ulcer prevention strategy, daily monitoring of plantar pressures underneath the foot (but also above the toes and wherever lesions may develop) is necessary, and as the present study points out, it also needs to be coupled with a predictive and personalized biomechanical model of the foot. Although plantar pressure monitoring is feasible in a laboratory or clinical setting using a heavy and expensive pressure sensor (such as the Zebris platform used in this article), it is impossible to implement for a large number of patients on a daily basis. To achieve this goal, a lighter and less expensive pressure sensor (such as the ones proposed by Novel, Tekscan, Vista Medical or Texisense) able to monitor the patient's foot pressures in his/her daily tasks must be used. For example, using a technology similar to that recently employed for the conception of the TexiCare device dedicated to the prevention of seated buttock pressure ulcers for people with spinal cord injury [START_REF] Chenu | An innovative embedded device for pressure ulcer prevention. Preliminary results with a paraplegic volunteer[END_REF], a "Smart Sock" (Bucki et al., 2011) (Perrier et al. 2014b), has been developed. It is made of a 100 % textile pressure sensing fabric wirelessly connected to a controller which can record and monitor the pressures all around the foot (not only under the sole). It can be used continuously during everyday activities. Once coupled to personalized biomechanical models such as the ones presented here, this device could be used to estimate the internal strains and to raise an alert, should the risk factor exceed a predefined personalized threshold. The main limitation at the moment lies in the implementation on a mobile platform of the complex biomechanical models such as the ones described here, where large non-linearities (mechanical, geometrical, contacts) need to be taken into account. At the moment the simulations require approximately three hours to converge to a steady state for a single pressure pattern applied to the foot sole. Our team is working towards drastically reducing the computational complexity of the models -while retaining most of their accuracy -in order to bring this technology into the clinical practice and benefit to the largest number of users: personalized prevention algorithms seamlessly embedded in wearable devices.

Conclusion

A workflow for generating a patient-specific biomechanical model of the foot has been presented and evaluated in this article in the context of pressure ulcer prevention. The technique implements an atlas based approach where anatomical knowledge is automatically transferred to the patient's modeling space using a non-linear registration algorithm. A new paradigm for the assessment of the level of tissue suffering in the context of pressure ulcer prevention has been proposed. It is based on the most recent consensus which relies on the measurement of Von Mises equivalent strains. The paradigm suggests looking at the volumes of "clusters" of elements undergoing a deformation greater than a predefined threshold (20% and 50% in the literature). This approach eludes the erratic results yielded by the monitoring of maximal VM strain values and opens the way to comparative assessment of a risk score that can be used to drive medical device design or clinical studies on a given population. The new paradigm however raises a new question which is "what is the minimal volume of tissue undergoing damage that is likely to lead to a lesion?" The answer to this question is beyond the scope of this study and will most likely not be a single figure but rather a threshold to be investigated based on the nature of the tissues (muscle or fat), the clinical condition and history of the patient, and other extrinsic factors that still have to be identified.

The approach was assessed on three patients and demonstrates the feasibility of patientspecific model generation. The evaluation was carried out by simulating the deformation of the personalized biomechanical models under the influence of a static common pressure pattern applied below the foot and by measuring the resulting internal strains. The results were further regionalized by dividing the foot into eight functionally meaningful regions. The results indicate that, for the chosen pressure pattern, the main risk of pressure ulceration is located below the heel for all four datasets (three patients and the atlas). The analysis shows that cluster analysis is an interesting alternative to the peak VM strain alone (as this value is strongly affected by numerical uncertainties inherent to numerical methods) and could be used to predict the risk of pressure ulcer and its localization within the foot regionalized representation. The study also confirms the influence of the patient's morphology on the range of the VM strains and associated cluster volumes: for the same pressure pattern, various values are obtained for both criteria, on all four datasets.

Before implementing this pressure ulcer prevention technique in a clinical workflow, some aspects of the approach still require improvement: the personalization of the patient's material properties for the various soft tissues layers, the precision of the registration on the toes, and the measurement of the pressure below the foot using a flexible textile sensor in real time to allow the patient to use this prevention tool on a daily basis. Figure 5 -Maximum strain and maximum cluster volume repartition for each anatomical region and for each patient, according to the two strain thresholds of 20 and 50 %. For the strains above a threshold of 20 % and above a threshold of 50 %, red colors mean strains above 200 %. For cluster volume above a threshold of 20 %, red means volume above 74 cm 3 . And for a cluster volume above a threshold of 50 %, red means volume above 6 cm 3 . Figure 6 -Location of the strains above the threshold of 20 % below the calcaneus of each patient (a blue color represents strains close to 20 % while red is the maximum, around 110 %). The morphological variation from one patient to another can also be observed.

Table 1 -Maximal Von Mises (VM) strains (in %) in the cluster of maximum volume and cluster volumes (in cm 3 ) above 20% for the atlas and the three considered patients (FC, FP and BR). The eight rows in the table correspond to the eight anatomical regions. The mean and standard deviation (in percentage points, ppt) for these eight regions are shown for both maximal VM and maximal cluster volumes.

Table 2 -Maximal Von Mises (VM) strains (in %) in the cluster of maximum volume and cluster volumes (in cm 3 ) above 50% for the atlas and the three considered patients (FC, FP and BR). The eight rows in the table correspond to the eight anatomical regions. The mean and standard deviation (in percentage points, ppt) for these eight regions are shown for both maximal VM and maximal cluster volumes.

Figure captions : Figure 1 -

 :1 Figure captions: Figure 1 -Cross section of the FE mesh representing the foot soft tissues: plantar skin layer, muscle layer, and the fat in-between. The white sections represent the locations of the rigid bodies modeling the bones.

Figure 2 -

 2 Figure 2 -Distribution of pressures applied under all foot models used in the study. The highest recorded pressures (14.5 N.cm -2 ) appear in red (below the heel). Lower values are shown in blue.

Figure 3 -

 3 Figure 3 -The eight anatomical regions defined to partition our foot model: (1) the Achilles tendon, (2) the top of the foot, (3) the heel, (4) the medial foot, (5) the first metatarsus,[START_REF] Chenu | An innovative embedded device for pressure ulcer prevention. Preliminary results with a paraplegic volunteer[END_REF] the four other metatarsi,[START_REF] Gefen | Plantar soft tissue loading under the medial metatarsals in the standing diabetic foot[END_REF] the hallux,[START_REF] Gefen | The biomechanics of heel ulcers[END_REF] and the four other toes.

Figure 4 -

 4 Figure 4 -The three personalized biomechanical models are derived from the atlas model using the MMRep algorithm which computes a non-linear correspondence function between atlas and patient anatomical landmarks.

  Figure 1

Table 2


Paula Méndez-Fernandez
email: paula.mendez_fernandez@outlook.com

Benoit Simon-Bouhet

Paco Bustamante

Tiphaine Chouvelon

Marisa Ferreira

Alfredo L López

Collin F Moffat

Graham J Pierce

Marie L Russell

Maria B Santos

Colin F Moffat

Jérôme Spitz

José V Vingada

Lynda Webster

Fiona L Read

Angel F González

Florence Caurant

Inter-species differences in polychlorinated biphenyls patterns from five sympatric species of odontocetes: Can PCBs be used as tracers of feeding ecology?

Keywords: Persistent organic pollutants, Biogeochemical tracers, Multivariate analysis, Cetaceans, Northwest Iberian Peninsula

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

Introduction

Marine mammals, occupying the upper trophic levels of marine food webs and possessing large lipid reserves, are susceptible to considerable bioaccumulation of lipophilic (fat-soluble) pollutants [START_REF] Tanabe | Capacity and mode of PCB metabolism in small cetaceans[END_REF][START_REF] Tanabe | Global contamination by persistent organochlorines and their ecotoxicological impact on marine mammals[END_REF][START_REF] Boon | A model for the interpretation of chlorinated biphenyl concentrations in marine mammals[END_REF][START_REF] Niimi | PCBs in Aquatic organisms[END_REF], such as the polychlorinated biphenyls (PCBs). Commercial PCB mixtures have been produced for a wide range of industrial applications because of their properties, which include resistance to breakdown by other chemicals [START_REF] Who | Evaluation of certain food additives[END_REF][START_REF] Who | Plychlorinated biphenyls and terphenyls[END_REF]. These properties have greatly contributed to the ubiquitous distribution of PCBs in the atmospheric, terrestrial and aquatic environments [START_REF] Niimi | PCBs in Aquatic organisms[END_REF]. In the aquatic environment, overall PCB profiles change from those of the industrial products as a result of distinct congener behaviour and relative rate of degradation through physico-chemical processes in the environment, and/or the metabolic action of organisms which ingest the contaminants from the food chain or acquire them directly from the water column and marine sediments [START_REF] Danis | Delineation of PCB uptake pathways in a benthic sea star using a radiolabelled congener[END_REF](Danis et al., , 2005a(Danis et al., , 2005b)).

Due to the variable number and position of the chlorine atoms on the biphenyl nucleus, individual PCB congeners follow different metabolic pathways. This results in the formation of diverse metabolites [START_REF] Letcher | Methyl sulfone and hydroxylated metabolites of polychlorinated biphenyls[END_REF] and different accumulation patterns [START_REF] Boon | A model for the interpretation of chlorinated biphenyl concentrations in marine mammals[END_REF]. Thus, the PCB concentration patterns observed in marine mammals differ not only from the patterns seen in the technical formulations originally released to the environment (e.g.

Arochlor, Clophen, Kanechlor, Pyralene) but also from those in their prey [START_REF] Muir | Organochlorine contaminants in Arctic food chains: Accumulation of specific polychlorinated biphenyls and chlordane-related compounds[END_REF].

Moreover, PCB patterns also differ among areas as a function of the type and of their distance from the source. This is a result of heavier congeners (i.e. those with higher degrees of chlorination) from regional sources adhering to organic particles and remaining closer to the source. The long-range signature includes a higher proportion of lighter congeners as a consequence of these being the more volatile PCBs and thus subject to atmospheric transport taking them some distance from their source (e.g. [START_REF] Staudinger | A critical review of Henry's law constants for environmental applications[END_REF][START_REF] Wania | Global fractionation and cold condensation of low volatility organochlorine compounds in polar regions[END_REF]. Therefore, PCBs can be used as a regional signature with coastal species being enriched in heavier congeners relative to oceanic species, which should present a greater proportion of the lighter compounds (e.g. [START_REF] Ross | Harbour seals (Phoca vitulina) in British Columbia, Canada, and Washington State, USA, reveal a combination of local and global polychlorinated biphenyl, dioxin, and furan signals[END_REF].

For seabirds and marine mammals, exposure to persistent organic pollutants through diet is the only relevant exposure pathway [START_REF] Borgå | Biological and chemical factors of importance in the bioaccumulation and trophic transfer of persistent organochlorine contaminants in arctic marine food webs[END_REF]. Thus, variation in prey preferences and/or in location of feeding grounds will result in varying tissue concentrations and patterns among different species, and among different individuals of the same species. In the same way, individual and species differences in the capacity to metabolize the different PCB congeners will result in intra-and inter-specific differences in PCB profiles in marine mammals [START_REF] Boon | The kinetics of individual polychlorinated biphenyl congeners in female harbor seals (Phoca vitulina), with evidence for structure related metabolism[END_REF][START_REF] Tanabe | Capacity and mode of PCB metabolism in small cetaceans[END_REF][START_REF] Wells | Techniques for pattern recognition of organochlorine residues in sea mammals from Scottish coastal waters[END_REF]. Previous studies, which investigated the relative metabolic degradation of PCBs in marine mammals, demonstrated that pinnipeds had a higher capacity to metabolize PCBs, especially the highly chlorinated congeners, than cetaceans [START_REF] Boon | The toxicokinetics of PCBs in marine mammals with special reference to possible interactions of individual congeners with the cytochrome P450-dependent monooxygenase system -An overview[END_REF][START_REF] Boon | Concentration dependent changes of PCB patterns in fish-eating mammals: structural evidence for induction of cytochrome P450[END_REF][START_REF] Weijs | Inter-species differences for polychlorinated biphenyls and polybrominateddiphenyl ethers in marine top predators from the Southern North Sea: Part 1. Accumulation patterns in harbour seals and harbour porpoises[END_REF]. Such differences seem to be the consequence of differences in cytochrome P450-mediated mono-oxygenase activities. Indeed, differences were not only found between pinnipeds and cetaceans, but also among different species of the same taxonomic group [START_REF] Goksøyr | Cytochrome P450 in seals: monooxygenase activities, immunochemical cross-reactions and response to phenobarbital treatment[END_REF][START_REF] Wells | Determination of individual chlorobiphenyls (CBs), including non-ortho, and mono-orthochloro substituted CBs in marine mammals from Scottish waters[END_REF].

There are sex-related differences in metabolism and transfer of PCBs during pregnancy and lactation that may also affect contaminant patterns in marine mammals. It is generally accepted that maternal transfer of lipophilic contaminants to the foetus and calves (during pregnancy and lactation) reduces the concentrations of these compounds in adult female marine mammals relative to males (e.g. [START_REF] Aguilar | Reproductive transfer and variation of body load of organochlorine pollutants with age in fin whales (Balaenoptera physalus)[END_REF]. Indeed, several studies suggested that lesser chlorinated PCB congeners were preferentially eliminated from females to calves through lactation (e.g. [START_REF] Subramanian | Bioaccumulation of Organochlorines (PCBs and p,p'-DDE) in Antarctic Adelie Penguins Pygoscelis adeliae Collected During a Breeding Season[END_REF][START_REF] Desforges | Transplacental transfer of polychlorinated biphenyls and Polybrominateddiphenyl ethers in arctic beluga whales (Delphinapterus leucas)[END_REF], resulting in a different PCB profile in females compared to males. These ecological, biological and physiological factors lead to different PCB signatures in individual marine mammal species [START_REF] Aguilar | Using organochlorines pollutants to discriminate marine mammal populations: a review and critique of methods[END_REF][START_REF] Aguilar | Reproductive transfer and variation of body load of organochlorine pollutants with age in fin whales (Balaenoptera physalus)[END_REF]. This potentially allows a set of persistent organic pollutants to be used as a tool for the identification of different ecological features of marine mammal species or for identifying populations within the species [START_REF] Aguilar | Using organochlorines pollutants to discriminate marine mammal populations: a review and critique of methods[END_REF][START_REF] Litz | Fine-scale spatial variation of persistent organic pollutants in Bottlenose dolphins (Tursiops truncatus) in Biscayne Bay, Florida[END_REF]Yunker et al., 2011).

In this context, the first objective of the present study was to assess whether the variation in feeding preferences and sources of contamination among five sympatric odontocete species is reflected in their PCB profiles. The second objective was to assess whether or not there is a minimum number of PCB congeners that together efficiently discriminated between species and can thus potentially be used as ecological tracers. To meet these objectives, we used previously published data on PCB concentrations, lipid content, stable isotope ratios (δ 13 C and δ 15 N), cadmium (Cd) concentrations, age and reproductive status [START_REF] Méndez-Fernandez | Foraging ecology of five toothed whale species in the Northwest Iberian Peninsula, inferred using δ 13 C and δ 15 N isotopic signatures[END_REF], 2013, 2014) from 120 individuals of five odontocete species (i.e. common dolphin Delphinus delphis, harbour porpoise Phocoena phocoena, bottlenose dolphin Tursiops truncatus, striped dolphin Stenella coeruleoalba and long-finned pilot whale Globicephala melas) from the Northwest Iberian Peninsula. Stable isotope and Cd values were used as proxies of habitat and diet. δ 13 C values are widely used in the literature as a proxy of feeding habitat, and more specifically neritic/benthic vs oceanic/pelagic (e.g. [START_REF] Deniro | Influence of diet on the distribution of carbon isotopes in animals[END_REF][START_REF] Fry | Stable isotope ecology[END_REF]. δ 15 N values indicate trophic level (e.g. [START_REF] Deniro | Influence of diet on the distribution of carbon isotopes in animals[END_REF][START_REF] Fry | Stable isotope ecology[END_REF] and also vary with feeding habitat (e.g. [START_REF] Chouvelon | Revisiting the use of δ 15 N in meso-scale studies of marine food webs by considering spatio-temporal variations in stable isotopic signaturesthe case of an open ecosystem: the Bay of Biscay (North-East Atlantic)[END_REF][START_REF] Ruiz-Cooley | Contrasting C and N isotope ratios from sperm whale skin and squid between the Gulf of Mexico and Gulf of California: effect of habitat[END_REF].

Cadmium concentrations can be used as a tracer of the type of prey eaten, specifically fish vs cephalopod consumption, since cephalopods concentrate more cadmium [START_REF] Honda | Heavy metal concentrations in muscle, liver and kidney tissue of striped dolphin, Stenella coeruleoalba, and their variations with body length, weight, age and sex[END_REF][START_REF] Bustamante | Cephalopods as a vector for the transfer of cadmium to top marine predators in the north-east Atlantic Ocean[END_REF][START_REF] Lahaye | Long-term dietary segregation of common dolphins Delphinus delphis in the Bay of Biscay, determined using cadmium as an ecological tracer[END_REF].

Material and methods

Study area and sampling

Fieldwork was carried out in the Northwest Iberian Peninsula (NWIP), from the Northern limit of the Galician coast in Spain (43° 31' N, 7° 2' W) to Nazaré on the Portuguese coast (39° 36' N, 9° 3' W; Fig. 1). Experienced members of the Spanish (Coordinadora para o Estudo dos Mamiferos Mariños, CEMMA) and Portuguese (Sociedade Portuguesa de Vida Salvagem, SPVS) stranding networks have been collecting stranded and by-caught cetaceans for over twenty-five years and over fifteen years, respectively. Animals were identified to species, measured, sexed and, if the decomposition state of the carcass allowed, full necropsies were performed and samples collected whenever possible. However, only animals recovered in a "fresh" state (a score of between 1 and 3 from the European Cetacean Society protocol: stranded alive, freshly dead or mildly decomposed; Kuiken and García Hartmann, 1991) were selected for analyses.

The age data are from [START_REF] Méndez-Fernandez | Ecological niche segregation among five toothed whale species off the NW Iberian Peninsula using ecological tracers as multi-approach[END_REF], 2014). The procedure for age determination consisted of counting Growth Layer Groups (GLGs) from tooth sections, assuming that one GLG equals 1 year (as described by [START_REF] Lockyer | A report on patterns of deposition of dentine and cement in teeth of pilot whales, genus Globicephala[END_REF][START_REF] Hohn | Protocol for obtaining age estimates from harbour porpoise teeth. Appendix 3, Report of the harbour porpoise age determination workshop[END_REF].

Lipid determination and polychlorinated biphenyl (PCB) analysis

Blubber samples were previously analysed for polychlorinated biphenyls and for lipid content [START_REF] Méndez-Fernandez | An assessment of contaminant concentrations in toothed whale species of the NW Iberian Peninsula: Part I. Persistent organic pollutants[END_REF]. Briefly, approximately 200 mg of blubber was cut (i.e. a vertical section of the full thickness), homogenised and mixed with pre-cleaned sodium sulphate (~ 20 g) and then extracted by Pressurised Liquid Extraction (PLE). Following PLE, lipid content was determined by gravimetry and the rest of the extract was concentrated by Syncore (fitted with a flushback module) to ~ 0.5 mL and passed through silica columns, before transferring, with washings, to amber glass gas chromatography (GC) vials. Then, the concentrations of 32 PCB congeners (International Union of Pure and Applied Chemistry PCB Nos 28,31,52,49,44,74,70,101,99,97,110,123,118,105,114,149,153,132,137,138,158,128,156,167,157,187,183,180,170,189,194 and 209) were determined by Gas Chromatography Electron Impact Mass Spectrometry (GC-EIMS) using an HP6890 series gas chromatograph interfaced with an HP5975 MSD.

The methods employed were validated by the replicate analysis of standards and samples, regular blank controls, and through spiking experiments or analysis of certified and laboratory reference materials (LRM; cod liver oil), following the same procedure as in [START_REF] Méndez-Fernandez | An assessment of contaminant concentrations in toothed whale species of the NW Iberian Peninsula: Part I. Persistent organic pollutants[END_REF].

Data treatment

The similarities/differences in PCB patterns for different species and/or different groups of individuals are commonly difficult to discern in graphical output arising directly from the chemical analysis. Data reduction (e.g. dimensionality reduction) is a way of identifying patterns more readily and to provide a descriptive overview. Multivariate analyses were used to infer patterns and differences in PCB profiles. These techniques (e.g. principal component, redundancy and discriminant function analysis) have been used in previous studies for detecting differences in organochlorine compound patterns between species of marine mammals, birds and other vertebrates as well as in marine invertebrates (e.g. [START_REF] Schwartz | Chemometric comparison of polychlorinated biphenyl residues and toxicological active polychlorinated biphenyl congeners in the eggs of Forster's Terns (Sterna fosteri)[END_REF][START_REF] Storr-Hansen | Patterns of chlorinated biphenyl congeners in harbour seals (Phoca vitulina) and their food. Statistical Analysis[END_REF][START_REF] Boon | Concentration dependent changes of PCB patterns in fish-eating mammals: structural evidence for induction of cytochrome P450[END_REF][START_REF] Borrell | Differences in DDT and PCB residues between common and striped dolphins from the Southwestern Mediterranean[END_REF].

All the PCB congener concentrations were normalized to the lipid-based content of the blubber and thus expressed in µg/g lipid weight (lw). For statistical analyses only those congeners for which the concentrations contributed more than 10% of the total sum of the 32 congeners were selected. Thus, a total of 22 PCBs were included: 31,52,49,44,70,101,99,110,118,105,149,132,138,156,187,183,180,170,189 and 194.

Differences in PCB patterns among species

Discriminant analysis (DA) was used to examine differences among the PCB patterns of the species using data from the 22 PCB congeners. However, in order to avoid problems due to collinearity of variables (i.e. concentrations of different congeners) we carried out a DA on principal component analysis (PCA) scores rather than on raw data. In addition, to better show the differences in the accumulation pattern of each PCB congener and to remove differences in absolute values of concentration between samples, PCBs were normalized to the concentration of 2',4,4',5, was used to calculate these ratios (i.e. CB ratios) since it is the dominant congener in the great majority of aquatic mammals, and because it has a molecular structure that makes it highly resistant to biotransformation [START_REF] Duinker | Individual chlorinated biphenyls and pesticides in tissues of some cetacean species from the North Sea and the Atlantic Ocean; Tissue distribution and biotransformation[END_REF][START_REF] Wells | Determination of individual chlorobiphenyls (CBs), including non-ortho, and mono-orthochloro substituted CBs in marine mammals from Scottish waters[END_REF][START_REF] Boon | The toxicokinetics of PCBs in marine mammals with special reference to possible interactions of individual congeners with the cytochrome P450-dependent monooxygenase system -An overview[END_REF].

Thus, the CB ratio was calculated as follows (values vary from 0 to 1):

CB ratiox = [CBx] / [CB-153]
Where [CBx] is the concentration of an individual CB congener x and is the concentration of CB-153.

Can we identify the most suitable PCBs to be used as tracers?

In order to identify the CB ratios that together efficiently discriminated between species and may thus be used as tracers; a re-sampling procedure was implemented based on the DA results.

With this aim, first the most relevant principal components from the PCA that were used as input data in the DA were selected using the "elbow method" (i.e. the principal components on the left of the visual break or elbow were selected). Then, a threshold of 10% of absolute contribution was used to select the CB ratios that contributed most to each principal component. Thus, our initial DA on principal components indicated that 9 CBs were most influential. Subsequently every possible DA using a combination of 9 CB ratios out of the 22 available (i.e. 497,420 combinations) was computed, and for each DA the same computation of assignment errors was made. Afterwards, an assignment errors histogram was plotted to compare the original assignment error with the distribution obtained.

Finally, in order to identify if there was a minimum number of congeners needed to efficiently discriminate pairs of species based on DA results, a resampling procedure was implemented.

For this, a discriminant analysis was produced with all CB ratios, and from 1 to 20 congeners were sequentially removed. Subsequently, for each random removal (i.e. from 1 to 20 congeners) the DA was repeated 1,000 times. For each DA, the mean Euclidean distance was calculated over the first 4 axes of the DA between the centroids of all pairs of species. Finally, the mean Euclidean distance of the pairs of species were plotted with their corresponding 95% confidence intervals as a function of the number of CB ratios used.

Factors influencing PCB patterns

To examine the relationship between CB ratios and the set of potential ecological and biological explanatory factors, redundancy analysis (RDA) was used. Redundancy analysis is a principal component analysis in which the axes are restricted to be linear combinations of explanatory variables [START_REF] Zuur | Principal component analysis and redundancy analysis[END_REF]. A RDA was performed using the combination of CB ratios leading to smallest numbers of assignment errors in the DA as response variables.

Moreover, RDA requires the number of samples to exceed the number of explanatory variables, preferably by a factor of two. In total, for the RDA we analysed 118 samples and we had eleven explanatory variables, which is an acceptable ratio. Those variables were: species, sex, age, individual reproductive status, lipid content, habitat (i.e. neritic/benthic vs oceanic/pelagic; represented by δ 13 C values), trophic level (i.e. high vs low trophic level; represented by δ 15 N values) and type of prey eaten (i.e. fish vs cephalopods; represented by renal cadmium concentrations). Sex and reproductive status were considered as categorical variables and species identity was coded as a series of dummy variables. Data for the ecological and biological variables used to perform the RDA were obtained from published studies conducted on the same individuals [START_REF] Méndez-Fernandez | Foraging ecology of five toothed whale species in the Northwest Iberian Peninsula, inferred using δ 13 C and δ 15 N isotopic signatures[END_REF], 2013, 2014). Significance testing in RDA is based on a permutation test so no assumption of normality is required and collinearity between explanatory variables is not an important issue [START_REF] Zuur | Principal component analysis and redundancy analysis[END_REF].

All the statistical analyses were performed using R version 3. 1.3 (R Development Core Team, 2014). The map (Fig. 1) was created using the MARMAP package [START_REF] Pante | Marmap: a package for importing, plotting and analysing bathymetric and topographic data in R[END_REF] and figures were created using the packages MASS and ggplot2 [START_REF] Venables | Modern Applied Statistics with S. Fourth Edition[END_REF][START_REF] Wickham | ggplot2: elegant graphics for data analysis[END_REF].

Results

Differences in PCB patterns among species

The average lipid content of blubber in the five species analysed was higher in females than in males (Table 1) although the difference was not statistically significant (p > 0.05). Total PCB concentrations (i.e. ΣPCB32) were higher in males than in females with the exception of harbour porpoises. Moreover, the species that showed the highest ΣPCB32 concentrations were bottlenose dolphin, pilot whale and harbour porpoise (Table 1).

The first two principal components of the PCA accounted for 84.0% of the CB ratio variability, while the first two-discriminant factors (i.e. LD1 and LD2) of the DA together explained 82.9% of the variance of the CB ratio data (Fig. 2). The DA revealed the existence of consistent differences among the five species, with harbour porpoises (shown in light blue) and striped dolphins (in orange) being clearly separated in this bi-plot. The three pilot whales were also distinctly placed relative to the striped dolphins, common dolphins and harbour porpoises, although there was overlap with one of the seven bottlenose dolphins. Between 66.7 and 100% of the individuals of each species were classified in the correct species group (Table 2). Individuals of common dolphin, bottlenose dolphin and harbour porpoise were all well classified to the correct species. However, one striped dolphin out of 15 (6.7%) was assigned to the common dolphins and one pilot whale out of three was assigned as a bottlenose dolphin.

For the first discriminant axis, which clearly separated striped dolphin and, to a lesser extent harbour porpoise, from the other species (Fig. 2), the most influential principal components of the PCA were, in decreasing order: 22 and 21 (Fig. S1a). For the second discriminant axis, which clearly separated pilot whales and harbour porpoise from the bottlenose dolphin and most of the common dolphins, the most influential principal components were, in decreasing order: 22, 21, 19, 20, 18 and 17 (Fig. S1b). Finally, CB-183, 189 and 70 contributed most to principal components 22 and 21,149,70,189,110,52,118 and 99 contributed most to principal components 19, 20, 18 and 17 (Table S1). All these 9 congeners had between 4 and 7 Cl atoms per biphenyl nucleus.

Can we identify the most suitable PCBs to be used as tracers?

The assignment errors histogram showed that the modal error rate using 9 congeners is around 12% and almost all combinations give error rates less than 20% (Fig. 3). Moreover, error rate of less than 6% was almost never achieved. However, some combinations of CB ratios led to very small numbers of assignment errors, including three combinations that led to only 4 assignment errors (3.3%). The PCBs present in these best combinations and that lead to low assignment errors were CB-105 (3 times), 110 (3), 170 (3), 180 (3), 187 (3), 194 (3), 49 (3), 189 (2), 31 (2), 138 [START_REF]32 congener concentrations (ΣPCB32 in µg/g lipid weight (lw) with range in parenthesis) measured in the blubber of males and females of five odontocete species from the Northwest Iberian Peninsula. n = sample sizes for each species and sex (total of 67 males and 51 females)[END_REF] and 183 [START_REF]32 congener concentrations (ΣPCB32 in µg/g lipid weight (lw) with range in parenthesis) measured in the blubber of males and females of five odontocete species from the Northwest Iberian Peninsula. n = sample sizes for each species and sex (total of 67 males and 51 females)[END_REF].

For all pairs of species, the mean Euclidean distance decreased with a decreasing number of congeners used in the DA (Fig. 4). Considering a specific number of congeners (e.g. 10), the mean Euclidean distances were different according to the pair of species compared, showing that some of them were easier to discriminate than others. As an example, the pairing of harbour porpoise with striped dolphin (Pp-Sc), and pilot whale with harbour porpoise (Gm-Pp), both exhibited high mean Euclidean distances. In comparison, the lowest mean Euclidean distances were shown when bottlenose dolphins (Tt) were paired with common dolphins, pilot whales or striped dolphins (Dd-Tt, Gm-Tt and Sc-Tt, respectively) and also when common dolphins were compared with striped dolphins (Dd-Sc), showing that these species were more difficult to discriminate. Approximately half of the comparisons presented in Figure 4 showed a clear plateau or at least a point of inflection, indicating that the analysis of more congeners would not result in a better discrimination. Where this was the case (e.g. Dd-Gm, Gm-Tr), the optimal number of congeners (showing the highest mean Euclidean distance) was between 10 and 15, which was consistent with the decreasing number of assignment errors related to the number of congeners used in the DA.

Factors influencing PCB patterns

The CB ratios that gave the best separation of the species and the smallest number of misassignment errors (i.e. 49,105,110,138,170,180,183,187,189 and 194) were used as response variables to compute a RDA together with the following ecological and biological parameters as explanatory variables: species, sex, age, reproductive status, lipid content, habitat, trophic level and prey type (see material and methods for details of explanatory variables). The first two axes explained 9.8% of the CB ratio variability (Fig. 5).

All explanatory variables tested had a significant influence (p < 0.05) with the exception of trophic level represented by δ 15 N values (Table 3). Additionally, either bottlenose dolphin or striped dolphin had no significant influence (p > 0.05) (Table 3).

The positions of individual animals in the bi-plot, represented by coloured points (Fig. 5), showed that harbour porpoises and bottlenose dolphins were positively correlated with the second axis as well as with lipid content and habitat. Most of the striped dolphins were associated with high renal cadmium concentrations and negatively correlated with habitat.

They were also associated with tetra-to hexa-chlorinated congeners (i.e. 99,110,118 and 138).

Examination of the bi-plot also highlighted a cluster of individuals in the lower part of the biplot corresponded to old common and striped dolphin females. These individuals were also associated in this part of the bi-plot with higher concentrations of the hepta-and octachlorinated congeners (i.e. 187 and 194), as well as with high cadmium concentrations and high age values. Finally, common and striped dolphin individuals showed the lowest segregation within the bi-plot (Fig. 5).

Discussion

Differences in PCB patterns and identification of most suitable PCBs to be used as tracers

The PCB profiles of the five odontocete species differed substantially and could be successfully distinguished from each other (Fig. 2), with all the common dolphins, bottlenose dolphins and harbour porpoises being correctly assigned to species (Table 2). Pilot whales exhibited the lowest percentage of successful classification. However, the sample size was very small for this latter species. As such, one out of three incorrectly classified animals resulted in only a 66.7% of well classification (Table 2). Harbour porpoise and striped dolphin were the species with the greatest separation on the DA plot and thus presented the most significant difference with respect to PCB profiles (Fig. 2). The habitat where these animals feed is likely to be a factor that strongly influences the PCB profiles. Higher pollutant burdens are generally found in species inhabiting coastal regions due to the close proximity of these animals to possible emissions, discharges and losses of persistent organic pollutants in temperate areas [START_REF] Storr-Hansen | Coplanar polychlorinated biphenyl congener levels and patterns and the identification of separate populations of Harbor Seals (Phoca vitulina) in Denmark[END_REF]. Moreover, the pattern of the PCBs may differ according to the distance from the source; the lighter congeners are more volatile, and thus are capable of being transported over longer distances [START_REF] Aguilar | DDT and PCB reduction in the western Mediterranean from 1987 to 2002, as shown by levels in striped dolphins (Stenella coeruleoalba)[END_REF]. Therefore, the proportion of the highly chlorinated congeners decreases with distance from the source. Additionally, the relative abundance of heavy to light carbon isotopes has also been used to discriminate between habitats, pelagic/offshore habitats where phytoplankton is the only source of organic carbon and vegetated inshore/benthic habitats where macrophytes may be an additional source of organic carbon [START_REF] Deniro | Influence of diet on the distribution of carbon isotopes in animals[END_REF][START_REF] Fry | Stable isotope ecology[END_REF]. Likewise,  15 N values also vary strongly with habitat (between inshore and offshore systems, with latitudes and between oceanic basins) and are used as an indicator of feeding habitat (e.g. [START_REF] Chouvelon | Revisiting the use of δ 15 N in meso-scale studies of marine food webs by considering spatio-temporal variations in stable isotopic signaturesthe case of an open ecosystem: the Bay of Biscay (North-East Atlantic)[END_REF][START_REF] Ruiz-Cooley | Contrasting C and N isotope ratios from sperm whale skin and squid between the Gulf of Mexico and Gulf of California: effect of habitat[END_REF]. Based on the stable isotope composition and direct observations at sea of the five species studied, we can affirm that in the waters off the NWIP the harbour porpoise is a coastal species, while the striped dolphin mainly inhabits offshore waters (Table 4). The discrimination between both species on the basis of their PCB profiles was consistent with the information available on features of their ecology. The other three species can be observed in both coastal and offshore waters, and are therefore likely to feed in both, as indeed may be inferred from their diets in this region (see Table 4 and references therein) and, based on this research, their PCB profiles.

Among the three combinations of PCBs that generated only four assignment errors, eleven of the selected PCBs (i.e. 110,170,180,187,194,49,189,31,138 and 183) presumably reflect interspecific differences in habitat but also in feeding preferences. In fact, a subsequent DA performed with this combination of PCBs confirm this statement giving a better separation of the species than if we consider the PCBs selected from the first DA (Tables S2 andS3).

Moreover, whatever the pair of species is considered, the mean Euclidean distance separating them decreases with the number of congeners randomly removed (Fig. 4). However, different pairs of species had different mean Euclidean distances between them as might be expected based on the eco-biological differences presented (see Table 4) and PCB profiles. The highest mean Euclidean distances between species correspond to harbour porpoise in comparison with striped dolphin and pilot whale, the species that are the most different from each other in terms of eco-biological traits. The lowest mean Euclidean distances correspond to the ecologically closest species, the common dolphin when compared with striped and bottlenose dolphins.

The assignment errors together with the mean Euclidean distances show that around 12 congeners is enough to well separate the species. This is a good result since PCB analysis requires considerable effort in terms of data reduction, quality assurance, and processing.

Potential quality control and comparability issues when applying this approach widely include variability in co-elution patterns, different sets of congeners analysed by different laboratories (although the ICES 7 52,101,118,153,138,180 is an internationally recognised set of PCBs), and normal inter-laboratory variation. Hence, the minimum number of congeners needed to allow discrimination of the species is an important issue since the cost of analyses can be greatly reduced by lowering the number of congeners that are analysed.

Although it was not possible to identify specific congeners as representing the most useful tracers, it is clear that among the set of 12 congeners selected, those that were more frequently repeated in the discriminant analyses performed (i.e. 138,183,189) should be included.

Factors influencing specific PCB patterns

The variability of PCB patterns in marine mammals basically depends on differences in their pollutant intakes, which is a direct consequence of the specific feeding preferences and associated habitats (e.g. [START_REF] Tanabe | Capacity and mode of PCB metabolism in small cetaceans[END_REF][START_REF] Borrell | Differences in DDT and PCB residues between common and striped dolphins from the Southwestern Mediterranean[END_REF], as well as capacities to metabolize and/or eliminate (e.g. maternal transfer) some congeners. However, identifying whether the separation is due to differences in PCB metabolism or differences in ecological and/or biological traits of the species is a difficult issue. This research revealed a good separation among species based on PCB profiles, that is in agreement with recent studies, using stable isotopes and other ecological tracers, showing that these species occupy different foraging niches within the NWIP [START_REF] Fernández | Interspecific and seasonal comparison of the niches occupied by small cetaceans off North-West Iberia[END_REF][START_REF] Méndez-Fernandez | Ecological niche segregation among five toothed whale species off the NW Iberian Peninsula using ecological tracers as multi-approach[END_REF]. Among the set of explanatory variables used to elucidate this issue (Fig. 5), only trophic level, represented by δ 15 N values, had no significant influence on CB ratios. More specifically, in the NWIP food web, these species have the same range of trophic level with no significant differences among them (Table 4). However, occupying the same trophic level does not necessarily mean that these species exploit the same food resource. This is reflected in the RDA since the prey type variable, which is represented by the renal cadmium concentrations of the specimens analysed, has a significant effect on the CB ratio data and is, together with habitat, the most important ecological explanatory variable (see Table 3).

Cadmium is a trace element used as an indicator of the type of prey eaten by the various odontocetes. More specifically, it helps distinguish fish vs cephalopod consumption, cephalopods being animals that concentrate more cadmium [START_REF] Honda | Heavy metal concentrations in muscle, liver and kidney tissue of striped dolphin, Stenella coeruleoalba, and their variations with body length, weight, age and sex[END_REF][START_REF] Bustamante | Cephalopods as a vector for the transfer of cadmium to top marine predators in the north-east Atlantic Ocean[END_REF][START_REF] Lahaye | Long-term dietary segregation of common dolphins Delphinus delphis in the Bay of Biscay, determined using cadmium as an ecological tracer[END_REF]. This would explain the positive correlation found between higher cadmium concentrations and the majority of the striped dolphins, some common dolphins and one of the three pilot whales (Fig. 5), which feed on cephalopods in the NWIP (Table 4). Moreover, common dolphins showed the highest dispersion, which is consistent with their mixed diet (Table 4), and the fact that the proportion of the different prey species in their diet can vary according to the oceanic or neritic origin of the individuals (Pusineri et al., 2007).

In the same way, [START_REF] Barone | PCBs and PCDD/PCDFs in fishery products: Occurrence, congener profile and compliance with European Union legislation[END_REF] compared PCB profiles of fishery products (i.e. fish, cephalopods and crustaceans) from Southern Italy, and highlighted a species-specific bioaccumulation of contaminants and differences in PCB profiles among the three different groups of seafood. Specifically, invertebrates (i.e. cephalopods and crustaceans) had a high percentage contribution of lower-chlorinated PCBs such as CB-28 and 52 with a slight predominance of CB-138, while in fish samples high-chlorinated congeners such as CB-180 were more prominent. Similarly, in the RDA, low chlorinated congeners such as CB-105, 49 31, 110 and 138 seems to be more associated with most of the striped and common dolphins (Fig. 5), and therefore with the mainly cephalopod feeders in the area. On the contrary, the high-chlorinated congeners CB-170 and 180 showed a positive relationship with most of the harbour porpoises and bottlenose dolphins, which are largely piscivorous in the NWIP (Fig. 5). [START_REF] Boon | Concentration dependent changes of PCB patterns in fish-eating mammals: structural evidence for induction of cytochrome P450[END_REF] reviewed the types of metabolic behaviour of several PCB congeners in five species of mammals (seals, otters and cetaceans) and revealed that CB170 and 180 are highly resistant to biotransformation and consequently difficult to metabolise. These congeners are also among those exhibiting a high octanol-water partition coefficient (log Kow ≈ 7) indicating their high liposolubility [START_REF] Jäntschi | Molecular descriptors family on structure activity relationships. Octanol-water partition coefficient of polychlorinated biphenyls[END_REF][START_REF] Walters | Trophic magnification of PCBs and its relationship to the octanol-water partition coefficient[END_REF]. The high lipid content exhibited by harbour porpoise and bottlenose dolphin (67.2 and 66.3% for males and 83.2 and 71.8% for females, respectively, Table 1) compared to the other species may thus partly explain the high ratio of these both congeners and the positive relationship with lipid content in these individuals (Fig. 5).

With respect to the various biological factors investigated, sex had an important effect on PCB accumulation and patterns, primarily due to the maternal transfer during gestation and lactation.

The physicochemical characteristics of the congeners govern this transfer, with the less lipophilic and lower molecular weight congeners being those primarily transferred to the foetus (e.g. [START_REF] Greig | Transplacental transfer of organochlorines in California sea lions (Zalophus californianus)[END_REF][START_REF] Desforges | Transplacental transfer of polychlorinated biphenyls and Polybrominateddiphenyl ethers in arctic beluga whales (Delphinapterus leucas)[END_REF]. Thus, the RDA bi-plot revealed a cluster of several individuals positively correlated with age and with kidney cadmium concentration (linked to diet), as well as with hepta-and octa-chlorinated CBs (CB-183, 187 and 194) (Fig. 5). Of the individuals belonging to this cluster, there are individuals of common and striped dolphin species, all of them are females and 60% were found to be mature. Thus, 187 and 194, which are high-chlorinated congeners, probably accumulates with age in females and may be not well transferred to offspring. Consequently, the older females, which have probably been pregnant and hence eliminated some less-chlorinated PCB congeners, may have high concentrations of these congeners as was previously reported in several marine mammal species (e.g. [START_REF] Desforges | Transplacental transfer of polychlorinated biphenyls and Polybrominateddiphenyl ethers in arctic beluga whales (Delphinapterus leucas)[END_REF][START_REF] Peterson | Effects of age, adipose percent, and reproduction on PCB concentrations and profiles in an extreme fasting North Pacific marine mammal[END_REF].

Despite the evident effects of eco-biological factors, the effect of metabolism on PCB patterns cannot be overlooked. Differences in metabolic capacities have been demonstrated between pinnipeds and cetaceans [START_REF] Tanabe | Capacity and mode of PCB metabolism in small cetaceans[END_REF][START_REF] Goksøyr | Cytochrome P450 in seals: monooxygenase activities, immunochemical cross-reactions and response to phenobarbital treatment[END_REF][START_REF] Wells | Chlorobiphenyl patterns in marine mammals from Northern European waters[END_REF], in which the iso-enzymes of the cytochrome P450 1A and 2B subfamilies (CYP1A and CYP2B) seem to play a key role in the biotransformation of PCBs. However, there are few comparable studies on odontocete species, and taxonomic proximity should not be taken as an indication of identical metabolic capacities. Indeed, substantial interspecific variation in the capacity to degrade persistent organic pollutants has been observed between common and striped dolphins [START_REF] Marsili | Organochlorine levels and mixed function oxidase activity in skin biopsy specimens from Mediterranean cetaceans[END_REF][START_REF] Borrell | Differences in DDT and PCB residues between common and striped dolphins from the Southwestern Mediterranean[END_REF], and between other taxonomically close vertebrates [START_REF] Fossi | Interspecies differences in mixed function oxidase activity in birds: Relationships between feeding habits, detoxification activities and organochlorine accumulation[END_REF]. The set of congeners used in the RDA, which resulted in a high segregation among species and low intra-species dispersion, could be involved in the species' metabolic capacities. Nevertheless, with the information currently available, it is impossible to distinguish the contribution of metabolism to the observed PCB patterns.

Conclusion

The results presented in this paper show clear differences in PCB patterns between different sympatric species of odontocetes. These differences are a presumably consequence of the effect of a mixture of metabolism and eco-biological parameters, with sex, age, habitat and the type of prey eaten probably being the most important among all that were tested. No single congener has been identified as a tracer of feeding ecology. However, 4 congeners (i.e. 138,183 and 189) from the 22 analysed seemed to be the prime ones and around 12 congeners appear to be enough to efficiently discriminate the species. Further studies are necessary to evaluate the potential generalisation of these results. Thus the same data treatment will have to be applied to other data sets of different species and/or different geographical areas in order to confirm the relevancy and the universality of these results. Moreover, the effectiveness of this data treatment for addressing issues associated with the conservation and management of wildlife could be significant. The importance of defining ecological diversity and populations has been the subject of much debate (ICES 2014), and genetic together with ecological segregation is the fundamental basis for identifying true populations. This is all the more important considering that the impact of localized anthropogenic threats, such as contamination, differs among populations. The profile of persistent organic pollutants, such as PCBs, was previously proposed as an additional tool to identify segregation of marine mammal populations (as in [START_REF] Borrell | Organochlorine compounds and stable isotopes indicate bottlenose dolphin subpopulation structure around the Iberian Peninsula[END_REF][START_REF] Herman | Feeding ecology of eastern North Pacific killer whales Orcinus orca from fatty acid, stable isotope, and organochlorine analyses of blubber biopsies[END_REF][START_REF] Krahn | Persistent organic pollutants and stable isotopes in biopsy samples (2004/2006) from Southern Resident killer whales[END_REF][START_REF] Pierce | Bioaccumulation of persistent organic pollutants in female common dolphins (Delphinus delphis) and harbour porpoises (Phocoena phocoena) from western European seas: Geographical trends, causal factors and effects on reproduction and mortality[END_REF]. In light of these results, the proposed data treatment can be used as a complementary tool in the identification of marine mammal populations. Adapted from [START_REF] Santos | Variability in the diet of bottlenose dolphin, Tursiops truncatus, in Galician waters, north-western Spain, 1990-2005[END_REF][START_REF] Santos | Long-term variation in common dolphin diet in relation to prey abundance[END_REF][START_REF] Santos | Quantifying the predation on sardine and hake by cetaceans in the Atlantic waters of the Iberian peninsula[END_REF][START_REF] Spitz | Diet overlap between harbor porpoise and bottlenose dolphin: an argument in favour of interference competition for food?[END_REF][START_REF] Spitz | Prey preferences among the community of deep-diving odontocetes from the Bay of Biscay, Northeast Atlantic[END_REF][START_REF] Pierce | Trends in cetacean sightings along the Galician coast, north-west Spain, 2003-2007, and inferences about cetacean habitat preferences[END_REF]Fernández et al., 2011a;[START_REF] Spyrakos | Spatio temporal patterns of marine mammal distribution in coastal waters of Galicia, NW Spain[END_REF][START_REF] Méndez-Fernandez | Foraging ecology of five toothed whale species in the Northwest Iberian Peninsula, inferred using δ 13 C and δ 15 N isotopic signatures[END_REF], 2013. 49,105,110,138,170,180,183,187,194 and 189)

Fig. 1 .Fig. 3 .

 13 Fig. 1. Map of the study area, Northwest Iberian Peninsula. The sampling area is delimited by Nazaré and Ribadeo. Depth contours (-100 m and -200 m) are shown in grey scale.

Fig. 4 .

 4 Fig. 4. Mean Euclidean distance and 95% confidence interval calculated between the centroids of all pairs of species when congeners are randomly removed. For each number of congeners removed we performed 1000 bootstraps. Common dolphin Delphinus delphis (Dde), harbour porpoise Phocoena phocoena (Pph), bottlenose dolphin Tursiops truncatus (Ttr), striped dolphin Stenella coeruleoalba (Sco) and long-finned pilot whale Globicephala melas (Gme).

Fig. 5 .

 5 Fig. 5. Results of redundancy analysis (RDA) on the best CB ratios combination (i.e. CB-31,

  in blubber of common dolphin Delphinus delphis, harbour porpoise Phocoena phocoena, bottlenose dolphin Tursiops truncatus, striped dolphin Stenella coeruleoalba and long-finned pilot whale Globicephala melas from the Northwest Iberian Peninsula. Bi-plot for axes 1-2 of significant explanatory and response variables. CdK (cadmium concentrations in kidney), C (δ 13 C values), N (δ 15 N values), AF (adult female), AM (adult male), JF (juvenile female), JM (juvenile male) and CBR (CB Ratios). The categorical variables (i.e. sex and reproductive status) are represented by empty black diamonds and CBR by red crosses.

Table 4 .

 4 Background on the feeding ecology of the odontocetes studied in the Northwest Iberian Peninsula, obtained from previous studies on: stomach contents analysis, land-based and sea surveys and chemical analyses (stable isotopes of carbon and nitrogen and renal cadmium concentrations). Trophic levels were calculated following VanderZanden and Rasmussen (2001) in[START_REF] Méndez-Fernandez | Foraging ecology of five toothed whale species in the Northwest Iberian Peninsula, inferred using δ 13 C and δ 15 N isotopic signatures[END_REF].

	Species	Habitat	Trophic level	Food source
	Common dolphin Delphinus delphis	Oceanic/neritic	4.7	Mixed feeder (blue whiting, sardine, hake and sand smelt, but also sepiolids, common and curled octopus)
	Harbour
	porpoise Phocoena	Neritic	5.3	Piscivorous (pouting, scad and blue whiting)
	phocoena
	Bottlenose dolphin Tursiops truncatus	Neritic (and offshore ecotype)	5.1	Mainly piscivorous (blue whiting and hake)
	Striped dolphin			Mainly teuthophagous (cephalopods and
	Stenella	Oceanic	4.3	crustaceans, but also blue whiting, sand
	coeruleoalba			smelt and scad)
	Long-finned
	pilot whale			Teuthophagous (common and curled
		Oceanic/neritic	4.9
	Globicephala			octopus)
	melas

Acknowledgements

The authors gratefully acknowledge the assistance of volunteers from the Galician (CEMMA) and Portuguese (SPVS) stranding networks. The authors would like to thank R. Gallois and C. Trichet for their participation on total lipid content analysis. P. Méndez-Fernandez was supported during the research period through a PhD grant from the Portuguese Foundation for Science and Technology (FCT) of the government of Portugal (SFRH/BD/36766/2007) and through a Science Without Borders (CSF) young talent postdoctoral grant of the Brazilian government. G. J. Pierce acknowledges support from the EU ANIMATE project (MEXC-CT-2006-042337), University of Aveiro and Caixa Geral de Depósitos (Portugal).

Species

N Age Lipid ΣPCB32

Male

Common dolphin 50 6.0 ± 4.9 59.0 ± 15.9 20.4 ± 16. 49, 105, 110, 138, 170, 180, 183, 187, 194 and 189) Table S3. Results of the discriminant analysis (DA) of the best CB ratios combination: 49,105,110,138,170,180,183,187,194 and