Several presentations addressed the interplay between metabolic and inflammatory risk factors in the development of obesity, type 2 diabetes (T2D) and atherosclerosis. The immune system impacts on these processes at several levels and inflammation in adipose tissue contributes to the development of insulin resistance. Gwen Randolph (Mount Sinai School of Medicine, New York, USA) showed that antigen-presenting dendritic cells (DCs) and macrophages within the muscular wall of collecting lymphatic vessels extend pseudopods into the lumen of adipose lymphatic vessels. Through these, they import lymph-derived antigens into adipose tissues in a CCL21/CCR7-dependent manner. In the absence of CCR7 expression, lymph-derived antigens were not found in adipose DCs or macrophages and inflammation was abated, suggesting that the collecting lymphatic ducts communicate with the surrounding adipose tissue and macrophages to control the immune response. Jaap Neels (University of Nice, France) discussed the roles of the IL8 (CXCL8)-CXCR2 chemokine-chemokine receptor system and paxillin/4-integrin interactions in monocyte/macrophage rolling and arrest on adipose tissue endothelial vessel walls. He demonstrated their important role in the subsequent macrophage infiltration of adipose tissue and development of insulin resistance. He also highlighted the role of CD11c + cells in the induction of high fat dietinduced insulin resistance and the expression of inflammatory markers. Deficiency or inhibition of 5-lipoxygenase (LO)-mediated leukotriene production altered M1/M2 content in adipose tissue and protected against insulin resistance. Diane Mathis (Harvard Medical School, USA) identified a specific subset of regulatory T cells (Tregs) in abdominal fat that express higher levels of PPAR1 than other Tregs. Abdominal fat Tregs produce high amounts of IL-10 and their deletion leads to increased inflammation, obesity and insulin resistance upon high-fat diet feeding. Luca Mazzarella (European Institute of Oncology, Milan, Italy) reported that insulin or IGF at high levels inhibit the induction of human Foxp3 + Tregs. This may lead to an increased inflammatory state.

Several other mechanisms were reported to be involved in inflammation in adipose tissue.

The plasma concentration of CD40L is elevated in T2D and metabolic syndrome patients.

Marjorie Poggi (Cardiovascular Research Institute, Maastricht, The Netherlands) reported that antibody-mediated depletion of CD40L improved metabolic control in a mouse model of a high-fat diet, Jean-François Tanti (INSERM U895, Nice, France) showed that macrophagederived IL-1 and TNF increase adipocyte expression of TPL2 (MAP3K8). This controls cytokine-induced lipolysis and down-regulates insulin signaling. Conversely, Tpl2 inhibition in adipocyte-macrophages co-cultures reduces production of inflammatory cytokines.

Anna Krook (Karolinska Institutet, Sweden) discussed cytokine control of skeletal muscle glucose metabolism. She showed that IL-6 and IL-13 are induced during exercise in muscle, inducing myotube formation, glycogen synthesis and glucose uptake in skeletal muscle cells.

Intriguingly, the skeletal muscle of T2D patients displays IL-6-resistance and secretes less IL-13, due to induction of the inhibiting microRNA Let7. David Patsouris (U. Lyon Sud, France) found that the skeletal muscle of obese mice displays increased levels of inflammatory markers (TNF, MCP1) along with the presence of macrophages. David Simar (U. New South Wales, Australia) compared insulin signaling in the monocytes of lean and obese individuals, the latter having decreased cell surface GLUT4 expression related to increased JNK, IKK and IRS1-Ser phosphorylation.

Several presentations addressed the regulation of high-density lipoprotein (HDL) production, which is protective against cardiovascular disease. Laurent Yvan-Charvet (Columbia U., USA) suggested an intriguing and novel regulatory role for ABCA1 and ABCG1. These are well known controllers of cholesterol efflux from macrophages. Both ABCA1 and ABCG1 genes are highly expressed on hematopoietic stem cells (HSC). The genetic deletion of these transporters in mice leads to proliferation of HSC, resulting in increased levels of circulating and infiltrating neutrophils, monocytes and eosinophils, with pro-atherogenic consequences. The pleiomorphism of the monocyte-macrophage lineage and a precise definition of their activation status were also major topics of discussion. David Hume (U. Edinburgh, UK) performed gene expression profiling in various monocyte populations to identify specific transcription networks. Surprisingly, many surface markers thought to be subpopulationspecific, such as CD11c (a DC marker), are rather more specific for the inflammatory response than the cell subset. Extensive expression profiling of mouse and human CSF-1treated macrophages showed that 10% of genes are discordantly expressed and 25% showed significant differences in expression levels. Many gene pathways were found to be differentially expressed in human and mouse macrophages, including the arginase and tryptophan metabolic pathways. Graham Thomas (U. Edinburgh, UK) studied modifications of gene expression profiles in alternative (M2) macrophages induced by nematode infection.

Parasite mainly affected pathways involved in energy metabolism, in particular arachidonic acid metabolism, with a large increase in 5-and 15-LO expression. Sophie Lotersztajn (Henri Mondor Hospital, Créteil, France) discussed the protective role of the endocannabinoid showed that activation of S6K, a kinase downstream of mTOR, is delayed in 'exhausted' CD8 + T cells responding to chronic viral infection, but is partially restored by ligation of the inhibitory receptor PD-1. This meeting highlighted the increasingly complex bi-directional interactions between immune cells and inflammation with metabolic alterations and cancer. Oxidative metabolism and nuclear receptors mechanistically link processes previously considered to be unconnected. Email: Bart.Staels@pasteur-lille.fr; david.dombrowicz@pasteur-lille.fr; ugrohmann@tin.it

1

  EMBO Workshop "Immunology and metabolism" Bart Staels, David Dombrowicz & Ursula Grohmann Recent findings have identified roles for the immune system in what had been considered prototypic metabolic diseases, while metabolic control has emerged as an important determinant of immune function. Giovanna Chimini, Lee Leserman, Diane Mathis and Philippe Naquet organized the EMBO Workshop on Immunology and Metabolism, which took place in January 2011 at the Centre d'Immunologie de Marseille-Luminy in Marseille, France. The meeting brought together around 100 scientists to discuss interactions between metabolism and inflammation, and was sponsored by Sanofi-Aventis and the CIML.

  Eric Moses (Texas Biomedical Research Institute, USA) showed that variants of the gene coding for the pantetheinase Vanin 1, might determine plasma HDL-C concentrations and therefore atherogenesis. Laure Perrin-Cocon (INSERM U851, Lyon, France) showed that HDL phospholipids inhibit DC maturation and LPS-induced NF-B activity, identifying a novel immune-modulatory and potentially anti-atherosclerotic mechanism of HDL.Ziad Mallat (U. Cambridge, UK) addressed the contribution of B lymphocytes to atherogenesis, and reported the beneficial effect of anti-CD20 mediated B cell depletion. This led to decreased atherosclerosis, reduced IFN- levels and IgG anti-oxidized LDL antibody production, but enhanced IL-17. Interestingly, CD20-depletion did not markedly influence anti-oxLDL IgM-producing B cells.Kenneth Rock (U. Massachusetts, USA) discussed the induction of inflammation by sterile stimuli, such as urate and cholesterol crystals, via activation of the NRLP3 inflammasome that matures IL-1. Phagocytosis of these crystals by macrophages triggers the inflammasome via lysosomal rupture and cathepsin B activation, a mechanism that may contribute to the induction of atherosclerosis. Host sterol metabolism is essential for controlling viral replication. Peter Ghazal (U. Edinburgh, UK) presented data on the impact on sterol metabolism of virus-induced activation of pattern recognition receptors (PRR), which suggests the existence of an innate immune mediated cross-regulation between infection and sterol metabolism.

  receptor CBR2 in alcoholic liver disease. The receptor prevents alcohol-induced M1 polarization of Kupffer cells and favors transition into an M2 phenotype. Réjane Paumelle (INSERM U1011, Lille, France) showed that deficiency of the anti-oncogene and cell cycle regulator p16INK4A polarizes macrophages from M1 to M2 in a manner independent of its cell cycle controlling activities. Liver Kupffer cells also belong to the monocyte/macrophage lineage. Ilse Scroyen (INSERM U626, Marseille, France) presented data linking fatty liver disease and LPS-induced plasminogen activator inhibitor-1 (PAI-1) with its hypoglycemic effect. Kupffer cells aggravate this process, possibly through the NFκB pathway, although they are not primary actors.Christopher Glass (U. California San Diego, USA) reviewed the crucial positioning of nuclear receptors at the crossroads between metabolic and inflammatory pathways in his Keynote address. He summarized the mechanisms of anti-inflammatory activities via trans-repression of PPAR, LXR and GR. Trans-repression by LXRs and PPAR occurs via their SUMOylated forms, which inhibit the removal of the co-repressor NCoR from inflammatory response gene promoters. Recent findings have identified Coronin2A (Cor2A) as an NCoR complex protein. Cor2A interacts with SUMO-LXR and is thus a SUMO-LXR docking site. Cor2A knockdown also resulted in decreased inflammatory gene expression and Cor2A was found to be necessary for NCoR clearance from the promoter via the recruitment of actin upon LPS stimulation. Bart Staels (INSERM U1011, Lille, France) discussed the counterintuitive presence of M2 macrophages in human atherosclerotic lesions. These macrophages display decreased cholesterol handling due to low LXR expression. Consequently, PPAR activation does not enhance ABCA1 expression, but stimulates phagocytic activity via regulating thrombospondin-1 expression. Antonio Castrillo (IIBM CSIC-UAM, Madrid, Spain) showed that LXR -/-mice lack marginal zone macrophages, potentially due to differences in bone marrow progenitors in LXR, but not LXR-null mice. As shown by Andres Hidalgo (CNIC, Madrid, Spain), the control of senescent neutrophil efferocytosis by macrophages enables LXR to have an unexpected role in the regulation of circadian oscillations in haematopoietic stem and progenitor cell (HSPC) transport, which is negatively correlated with bone marrow CXCL12 levels. Christine Delprat (CNRS UMR 5239, Lyon, France) showed that IL-17A promotes foam cell formation by regulating lipid metabolism genes as well as LXR and PPAR in Langerhans cell histiocytosis. Percy Knolle (U. Bonn, Germany) identified a role for PPAR in experimental autoimmune encephalomyelitis (EAE). PPAR activation delays disease onset via its effects on CD4 + T lymphocytes. PPAR -/-CD4 + T cells are prone to develop into Th17 cells due to increased expression of the nuclear receptor RORt. Signals upstream of the glucocorticoid receptor also play a role in metabolic diseases linked to inflammation. Katia Karalis (Children's Hospital, Harvard Medical School, USA) showed that LPS treatment of corticotropin-releasing hormone-deficient mice induces anorexia. Another theme of the meeting concerned the role of oxidative stress and the generation of reactive oxygen species (ROS) in the modulation of inflammatory responses. Anna Rubartelli (IST, Genova, Italy) showed that ROS production and subsequent activation of the antioxidant system are required for Toll-like receptor-induced IL-1 secretion by myeloid cells. An improper redox response results in altered inflammasome activation and IL-1 secretion. Shyam Biswal (Johns Hopkins U., USA) reported roles for Nuclear factor erythroid-2 related factor (Nrf2) and its inhibition by Kelch-like ECH-associated protein 1 in the expression of antioxidant and detoxification enzymes. Nrf2 inactivation ameliorated atherosclerosis due to decreased activation of metabolic enzymes, decreased glucose uptake and oxidation, and decreased fatty acid synthesis. In contrast, constitutive Nfr2 activation in cancer cells enhanced metabolism and inflammation. Alice Carrier (INSERM U624, Marseille, France) showed that ROS levels are controlled by Tumor Protein 53-Induced control of weight and triglycerides, whereas TORC1 drives organ-specific autoimmune/inflammatory diseases such as EAE. Andy Johnson (U. Pennsylvania, USA)

  The IDO-mTOR axis may represent a novel link between metabolic and immune/inflammatory pathways. Analyses of transcripts from particular immune subpopulations and cell-specific inactivation or overexpression of proteins controlling these pathways should yield exciting discoveries in this emerging field.Bart Staels and David Dombrowicz are at the Inserm U1011, University of Lille Nord de France and Institut Pasteur de Lille, France. Ursula Grohmann is at the University of Perugia, Italy.

Nuclear Protein 1, a p53-target which is repressed by miR155 in cancer cells, further establishing a link between redox state, metabolism and cancer. Leonore Herzenberg (Stanford U., USA) discussed the effect of N-acetylcysteine, an acetaminophen detoxifier, which increases GSH levels in individual cases of cystic fibrosis and HIV. She emphasized the beneficial effect of culturing cells under physiologic conditions (5% 0 2 ). Samuel Rommelaere (CIML, Marseille. France) studied the complex action of epithelial vanin 1, a pantetheinase that provides tissue cysteamine. Vanin 1-deficiency aggravates T1D but protects animals from acute colitis and modulates resistance to infectious disease. Amino acid catabolism is an ancestral survival strategy that controls immune responses in mammals. Indoleamine 2,3-dioxygenase (IDO) catalyzes the first step of tryptophan (Trp) catabolism. Ursula Grohmann (U. Perugia, Italy) showed that IDO may have evolved in mammals to acquire signaling functions for the induction and maintenance of long-term immune tolerance. In plasmacytoid DCs, IDO signaling requires TGF--induced tyrosine phosphorylation of two immune tyrosine-based inhibitory motifs (ITIM) contained in the noncatalytic domain of the enzyme. As shown by Andrea Nino Castro (U. Bonn, Germany), IDO is a key feature of regulatory DCs that exert microbicidal activity in a model of chronic inflammation sustained by Listeria cells. The effect does not depend on Trp starvation, but is mediated by tryptophan metabolites. Interestingly, human obesity is often associated with reduced plasma Trp and increased adipose and hepatic IDO expression (Odile Poulain-Godefroy; CNRS UMR 8199, Lille, France). A shift in Tregs/Th17 balance was associated with IDO increase in the subcutaneous adipose compartment, but was not observed in visceral adipose tissue. In T cells, IDO inhibits the mammalian target of rapamycin (mTOR), a central component of nutrient and hormone sensitive signaling pathways. mTOR signaling proceeds via two downstream complexes, TORC1 and TORC2. Using TORC1 -/-and TORC2 -/-mice, Jonathan Powell (Johns Hopkins U., USA) demonstrated that TORC2 is involved in the


Najla Dsouli
email: dsouli2005@yahoo.fr

Frédéric Delsuc

Johan R Michaux

Eric De Stordeur

Arnaud Couloux

Michel Veuille

Gérard Duvallet

Najla Dsouli-Aymes

Inornatus Grünbreg

S Transvittatus Villeneuve

S Omega Newstead

S Pallidus

Phylogenetic analyses of mitochondrial and nuclear data in haematophagous flies support the paraphyly of the genus Stomoxys

Introduction

Stomoxys flies are principally haematophagous, and are associated with livestock and wildlife throughout the world. The most studied species, Stomoxys calcitrans, is an economically important pest of cattle and several studies have attempted to estimate its impact on cattle production [START_REF] Miller | Effect of stable flies on feed intake and milk production of dairy cows[END_REF][START_REF] Campbell | Effects of stable flies (Diptera-Muscidae) on weight gains and feed-efficiency of calves on growing and finishing rations[END_REF][START_REF] Campbell | Effects of stable flies (Diptera: Muscidae) on weight gains of grazing yearling cattle[END_REF]. These flies represent a serious nuisance not only because of their painful bites and blood predation, but also because they are involved in the mechanical transmission of several pathogens, such as the Capripoxvirus causing a cattle disease (Lumpy skin disease), Anaplasma marginale the causative agent of severe bovine anaplosmosis, and Dermatophilus congolensis the causative agent of dermotophilosis [START_REF] Zumpt | The Stomoxyinae biting flies of the world. Taxonomy, biology, economic importance and control measures[END_REF][START_REF] D'amico | Are stable fly (Diptera: Stomoxyinae) vectors of Trypanosoma vivax in the Central African Republic[END_REF][START_REF] Foil | Mechanical transmission of disease agents by arthropods[END_REF].

The genus Stomoxys belongs to the tribe Stomoxyini in the subfamily Muscinae [START_REF] De Carvalho | Classificação de Muscidae (Diptera): una proposta através da análise cladística[END_REF][START_REF] Couri | Systematic relations among Philornis Meinert, Passeromyia Rodhain & Villeneuve and allied genera (Diptera, Muscidae)[END_REF]. This subfamily is part of the large family Muscidae comprising about 4,500 described species classified in 180 genera [START_REF] De Carvalho | A catalogue of the Muscidae (Diptera) of the Neotropical region[END_REF]. Within Muscidae, phylogenetic analyses have been conducted at several taxonomic levels to assess the relationships among constitutive species [START_REF] Couri | Cladistic analysis of Coenosiini (Diptera: Muscidae: Coenosiinae)[END_REF][START_REF] Couri | Systematic relations among Philornis Meinert, Passeromyia Rodhain & Villeneuve and allied genera (Diptera, Muscidae)[END_REF][START_REF] Schuhli | Revision and cladistics of the Neotropical genus Pseudoptilolepis Snyder (Diptera, Muscidae)[END_REF][START_REF] De Carvalho | Taxonomy, cladistics and biogeography of the South American Brachygasterina Macquart (Diptera: Muscidae)[END_REF][START_REF] Schuhli | Molecular phylogenetics of the Muscidae (Diptera:Calyptratae): new ideas in a congruence context[END_REF]). Recently, [START_REF] Nihei | Phylogeny and classification of Muscini (Diptera, Muscidae)[END_REF] carried out a cladistic analysis to assess the monophyly of the Muscini tribe. However, since the monograph of [START_REF] Zumpt | The Stomoxyinae biting flies of the world. Taxonomy, biology, economic importance and control measures[END_REF] who proposed the monophyly of Stomoxyini on the basis of morphological characters, and despite the medical and economic importance of these parasitic flies, no phylogenetic analysis has been carried out to assess the relationships within Stomoxyini. The flies of this tribe are easily recognized by their typical piercing/sucking mouthparts which differentiate them from the common housefly (Musca domestica) and relatives from the Muscini tribe. The Stomoxyini tribe consists of 10 genera and about 39 species reviewed by [START_REF] Zumpt | The Stomoxyinae biting flies of the world. Taxonomy, biology, economic importance and control measures[END_REF], the most medically and economically important species are members of the genera Haematobosca, Haematobia and Stomoxys.

The genus Stomoxys, which originated from the Old World, includes 18 species [START_REF] Zumpt | The Stomoxyinae biting flies of the world. Taxonomy, biology, economic importance and control measures[END_REF]. Among these species, only S. calcitrans has a worldwide distribution and is a synanthropic fly. All other species are exclusively tropical, twelve of which are located on the African continent, four on the Asian continent, and one species, S. sitiens Rondani, has been reported in both Africa and Asia [START_REF] Zumpt | The Stomoxyinae biting flies of the world. Taxonomy, biology, economic importance and control measures[END_REF] (Table 1).

In this study, we address the phylogenetic relationships of the Stomoxys genus using mitochondrial (COI, CytB) and nuclear ribosomal (ITS2) nucleotide sequences. Mitochondrial DNA (mtDNA) has been widely used in systematic and many universal PCR primers are available for genes like COI and CytB [START_REF] Avise | Molecular Markers, Natural History and Evolution[END_REF]. The ITS2 region is easy to amplify and have been used in previous phylogenetic studies of Diptera [START_REF] Hwang | Revisited ITS2 phylogeny of Anopheles (Anopheles) Hyrcanus group mosquitoes: re-examination of unidentified and misidentified ITS sequences[END_REF][START_REF] Thanwisai | Molecular phylogeny of black flies (Diptera: Simuliidae) from Thailand, using ITS2rDNA[END_REF]. This work explores congruence and information content within the different molecular datasets analysed using probabilistic methods of phylogenetic reconstruction. Phylogenetic results are discussed and taken into consideration to propose a taxonomic revision of the group. Also, we propose estimations of divergence times for major clades within Stomoxyini based on a relaxed molecular clock approach. This allowed discussing the evolutionary history of Stomoxys species in its biogeographical context.

Materials and methods

Taxon sampling

The material used in this study belongs to the Laboratoire de Zoogéographie of the Université Paul-Valéry in Montpellier (France), and to the Department of Entomology of the Natural History Museum in London (United-Kingdom). For Stomoxys we sampled 11 representative species or subspecies. As outgroups, we incorporated Prostomoxys saegerae, the unique representative species of Prostomoxys, and Haematobosca croceicornis, a newly described species from Gabon [START_REF] Pont | A new species of Haematobosca Bezzi from Gabon (Diptera: Muscidae)[END_REF]. For this study, we also included sequences of Haematobia irritans and Musca domestica (Muscidae), as well as species belonging to the more distantly related Drosophilidae, Syrphidae and Dolichopodidae families in order to provide calibration points for molecular dating analyses [START_REF] Wiegmann | Time Flies, a New Molecular Time-Scale for Brachyceran Fly Evolution without a Clock[END_REF]. Geographical origins of studied specimens and sequences accession numbers are listed in Table 2. Specimens are available upon request from the corresponding author and are stored in the Laboratoire de Zoogéographie of the Université Paul-Valéry in Montpellier (France).

DNA extraction, amplification and sequencing

Genomic DNA was extracted using the DNAeasy tissue Kit (QIAGEN) to a final volume of 180 µl. Amplifications by PCR, using the specific primer pairs described in Table 3, led to amplicons of different lengths for the COI (753 bp), CytB (587 bp) and ITS2 (~ 430 bp) regions. All PCR amplifications were performed in a 30µl reaction volume containing, at final concentrations, 200 µM dNTPs (diNucleotide Tri Phosphate), 10X buffer, 25µM of each primers and 0.5µl of Taq polymerase (Eurogentec Red GoldStar ® ), and 3 µl of purified DNA.

Thermal cycling conditions for PCR were as followed: initial denaturation at 94°C for 4 min., 35 cycles of denaturation at 94°C for 40 sec., annealing at 48-50°C for COI, 57-58°C for CytB, and 60-62°C for ITS2, and extension at 72°C for 1 min. A final elongation step at 72°C for 10 min completed the DNA amplification process. Ten specimens for each species were sequenced using Sanger sequencing on an ABI 3730 automatic sequencer at the Centre National de Séquençage (Génoscope) in Evry (France).

Phylogenetic analyses

The nucleotide sequences from COI, CytB, and ITS2 were automatically aligned using the multiple alignment program ClustalW 1.4 [START_REF] Thompson | CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position specific gap penalties and weight matrix choice[END_REF] using default parameters. Multiple sequence alignments were then adjusted by visual inspection, taking the sequences of S. calcitrans available in GenBank as a reference for each gene portion.

Alignments were cleaned from problematic alignment blocks using Gblocks 0.91 [START_REF] Castresana | Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis[END_REF] Probabilistic analyses were carried out on each individual datasets (COI, CytB, and ITS2). Maximum likelihood (ML) reconstruction was conducted using PAUP* 4.0b10 [START_REF] Swofford | PAUP*: Phylogenetic Analysis Using Parsimony (*and Other Methods)[END_REF]. The best-fitting models of sequence evolution for different partitions were determined based on the Akaike Information Criterion (AIC) as implemented in jModelTest [START_REF] Posada | jModelTest: phylogenetic model averaging[END_REF] using PHYML [START_REF] Guindon | A simple, fast and accurate algorithm to estimate large phylogenies by maximum likelihood[END_REF] for calculating likelihood scores.

ML heuristic searches were conducted with PAUP* using Tree Bisection Reconnection (TBR) branch-swapping on a Neighbor-Joining (NJ) starting tree using the best-fitting model and associated parameters selected by jModelTest. ML bootstrap proportions were obtained by repeating the same ML heuristic search on 100 pseudo-replicated datasets in order to evaluate the confidence for each node of the tree topology.

Crossed statistical SH tests [START_REF] Shimodaira | Multiple comparisons of log-likelihoods with applications to phylogenetic inference[END_REF] of congruence between the three genes were performed in PAUP* by testing the best ML topology obtained from each gene against the topology inferred from their concatenation. These tests were run using each individual gene dataset and the concatenated dataset (see Table 4 for details).

Bayesian phylogenetic inference was conducted using MrBayes 3.1.2 [START_REF] Ronquist | MrBayes 3: Bayesian phylogenetic inference under mixed models[END_REF]. Three a priori partition schemes were tested for analyzing the concatenation of the three genes (COI, CytB, and ITS2): one single partition, a 3 gene partition, and a 7 partition scheme distinguishing 6 partitions for each codon position of the two coding genes (COI and CytB) plus a single partition for ITS2. For each partition, we used the best-fitting model selected by jModelTest. Bayesian inference under each partition scheme was conducted with two independent runs of four incrementally-heated Metropolis Coupled Markov Chain Monte Carlo (MCMCMC) starting from a random tree. MCMCMC were run for 3,000,000 generations with trees and associated model parameters being sampled every 300 generations. The initial 2,000 trees in each run were discarded as burn-in samples and the harmonic mean of the likelihood was calculated by combining the two independent runs. The harmonic means were then used to compute the Bayes factor for the three possible partition comparisons (no partition versus 3 partitions by gene, no partition versus 7 partitions by codon and gene, and 3 partitions by gene versus 7 partitions by codon and gene). The best partition strategy was then determined by the 2ln Bayes factor criterion as recommended by [START_REF] Brandley | Partitioned Bayesian analyses, partition choice, and the phylogenetic relationships of scincid lizards[END_REF]. The 50% majority-rule consensus tree was then computed from the 16,000 trees sampled in the two independent runs under the best model.

Posterior probabilities greater than or equal to 95 % are generally regarded as strong support for a clade [START_REF] Wilcox | Phylogenetic relationships of the dwarf boas and a comparison of Bayesian and bootstrap measures of phylogenetic support[END_REF], but the correspondence with bootstrap support values is not exact [START_REF] Douady | Comparison of Bayesian and maximum likelihood bootstrap measures of phylogenetic reliability[END_REF].

Molecular datings

Bayesian estimation of divergence times under relaxed-clock models was conducted using the Multidivtime package [START_REF] Thorne | Divergence time and evolutionary rate estimation with multilocus data[END_REF]. We used the previously estimated Bayesian topology as the best hypothesis for Stomoxys phylogeny onto which divergence dates were estimated. The inclusion of species belonging to Drosophilidae (Drosophila virilis, D. melanogaster and D. yakuba), Syrphidae (Simosyrphus grandicornis, Cheilosia longula and C. naruska) and Dolichopodidae (Dolichopus longula and D. nubilus), which represent other brachyceran groups, allowed calibrating the tree. Indeed, based on a previously established timescale for brachyceran flies from Bayesian analyses of 28S rRNA data [START_REF] Wiegmann | Time Flies, a New Molecular Time-Scale for Brachyceran Fly Evolution without a Clock[END_REF], our tree was calibrated by using (1) the occurrence of Schizophora 87 million years ago (Mya) as the a priori expected number of time units between tip and root, and (2) the estimated date of the Drosophila/Musca split between 48 and 51 Mya as the calibration constraint. This previous molecular estimate was used in the absence of any relevant fossil calibration for Stomoxyini.

First, the program Baseml of the PAML package version 4.2b [START_REF] Yang | PAML 4: a program package for phylogenetic analysis by maximum likelihood[END_REF] was used to estimate the nucleotide frequencies, transition/transversion ratio, and rate heterogeneity among sites for the concatenated dataset. These values parameterize the F84+G nucleotide substitution model. Second, we calculated the branch lengths of the constrained topology and the associated variance-covariance matrix under this model using with the program Estbranches. Finally, the program Multidivtime was used to run a Markov chain Monte Carlo (MCMC) for estimating mean posterior divergence times on nodes with associated standard deviations and 95% credibility intervals (95% CredI) from the variance-covariance matrix produced by Estbranches. The MCMC was sampled 10,000 times every 100 cycles after a burn-in stage of 100,000 cycles. The prior for the expected number of time units between tip and root was set at 87 Mya (SD = 43 Mya). The estimated branch lengths obtained by Estbranches were used to estimate the median amount of evolution between the root and all the tips of the ingroup. Other priors for gamma distribution of the rate at root node (rtrate and rtratesd) and the Brownian motion constant describing the rate variation (brownmean and brownsd) were derived from the median branch length of the phylogram. The highest possible number of time units between tip and root (Bigtime) was set to 119 Mya which corresponds to the upper bound of the 95% CredI for Schizophora as estimated by [START_REF] Wiegmann | Time Flies, a New Molecular Time-Scale for Brachyceran Fly Evolution without a Clock[END_REF]. The age of the node Musca/Drosophila was constrained using an upper bound of 51 Mya and a lower bound of 48 Mya [START_REF] Wiegmann | Time Flies, a New Molecular Time-Scale for Brachyceran Fly Evolution without a Clock[END_REF].

Results

Sequence characteristics

Alignment statistics and models selected for phylogenetic analyses for each data partitions are summarized in Table 5. Alignment of the mitochondrial genes was straightforward, as no indels were introduced. The alignment of COI sequences resulted in

Phylogenetic relationships

We evaluated topological congruence among the individual genes by computing crossed SH tests in which the highest-likelihood topologies obtained with each individual datasets and their concatenation were compared against each other (Table 4). While there is apparent incongruence among the individual markers, none of the three individual datasets in fact significantly rejects the ML topology supported by their concatenation (see last line of Table 4). This indicates that combining the three individual gene datasets leads to a phylogenetic estimate that is compatible with the signal contributed by each individual gene.

We therefore decided to concatenate the three datasets in order to maximize the number of characters analysed and to enhance the phylogenetic signal as advocated by [START_REF] Delsuc | Molecular phylogeny of living xenarthrans and the impact of character and taxon sampling on the placental tree rooting[END_REF].

The total number of nucleotide positions in the concatenation was 1635. Bayes factor comparisons showed that the best-fitting model for the whole concatenation was the use of a single GTR+G model. Indeed, the 2ln Bayes factor criterion [START_REF] Brandley | Partitioned Bayesian analyses, partition choice, and the phylogenetic relationships of scincid lizards[END_REF] was always in favour of H0 in the three partition scheme comparisons (2lnBF < 0.03) indicating that it was not worth partitioning in our case. The phylogenies obtained from ML and Bayesian methods under a single GTR+G were identical, but the Bayesian analysis provided higher support values for most of the nodes (Fig. 1). Our results support the monophyly of bengalensis being sister-groups (PP = 1.0; BP ML = 91) (Fig. 1).

Estimation of Divergence times

The divergence time estimates are presented as a chronogram where branching nodes correspond to the mean age estimated from the posterior distribution and its associated 95% CredI (Fig. 2). The divergence time between the Stomoxys genus and its Haematobia + Haematobosca sister-clade is estimated around 30.8 Mya (95% CredI: 40.3-

22.3). The age estimate for the early emergence of S. indicus (clade A) within the genus

Stomoxys is estimated at about 27 ). The divergence between clade B and clade C occurred around 20.8 ). According to this inferred timescale, the major lineages within the genus Stomoxys were present by the late Oligocene, whilst the greatest amount of cladogenesis occurred during the Early Miocene. The two subspecies of S. niger (S. niger niger and S. niger bilineatus) separated around 16.3 Mya (24.2-10.3), a divergence time almost exactly similar to the one inferred for the separation between the species S. inornatus and S. transvittatus estimated at 16.4 Mya (24.6-10.0). Finally, the divergence between P. saegerae and S. varipes is estimated to have occurred 14.2 Mya (21.6-8.5) concomitantly with the separation of S. calcitrans from its S. bengalensis/S. sitiens sister-clade at 14.1 Mya (21.5-8.5).

Discussion

Phylogenetic relationships and taxonomy of the Stomoxys group

The phylogenetic analyses performed in the present study allowed reconstructing a phylogenetic framework for the major constitutive species of the genus Stomoxys and closely related genera. Our analyses strongly support the monophyly of Stomoxyini including a monophyletic group consisting of Haematobia and Haematobosca as a sister-clade to all remaining Stomoxyini. However, according to our results the genus Stomoxys sensu stricto is found to be paraphyletic due to the unexpected position of Prostomoxys saegerae which appears to be well nested within the Stomoxys group as a sister-group to S. varipes (Fig. 1). [START_REF] Zumpt | The Stomoxyinae biting flies of the world. Taxonomy, biology, economic importance and control measures[END_REF] created the genus Prostomoxys for the sole species P. saegerae. According to Zumpt's identification key for Stomoxyini, Prostomoxys is characterized by maxillary palps that are as long as the proboscis, whereas Stomoxys is characterized by palps that are shorter than half the length of the proboscis. This key includes no other diagnostic character for the distinction of this genus from Stomoxys. [START_REF] Zumpt | The Stomoxyinae biting flies of the world. Taxonomy, biology, economic importance and control measures[END_REF] mentioned that the plesiomorphic form of the Stomoxyini mouthparts can be accepted as the maxillary palps being about as long as the proboscis. Based on our phylogenetic study there appears to be no reason to incorporate this species in a separate genus as P. saegerae is in fact closely related to Stomoxys species. Prostomoxys should thus be synonymised with Stomoxys, and the species Prostomoxys saegerae must be renamed into Stomoxys saegerae. The long palp characters would be better considered as mere specific traits, since they are likely to be plesiomorphic.

Within Stomoxys lato sensu (including P. saegerae), our phylogenetic analyses strongly support the placement of S. indicus as the sister-group to the remaining Stomoxys species. Taxonomically, S. indicus is considered to be highly variable, with a number of synonyms used in the entomological literature. As mentioned by [START_REF] Zumpt | The Stomoxyinae biting flies of the world. Taxonomy, biology, economic importance and control measures[END_REF]: "it is possible that a numerical taxonomic study based on great numbers of specimens from various populations may reveal that several subspecies may be retained in the future". This species is described as the most common Stomoxys species in the Oriental region after S. calcitrans and has been recorded from India to the Pacific islands of Fiji and Samoa, including China and Japan. Its geographic range also reaches the Palaearctic region on the occidental side.

Moreover, it has also been recorded in a number of South-East Asian countries such as Thailand [START_REF] Masmeatathip | First survey of seasonal abundance and daily activity of Stomoxys spp. (Diptera: Muscidae) in Kamphaengsaen campus, Nakornpathom Province Thailand[END_REF], Malaysia, Viet Nam, Taiwan, Borneo, Sumatra, and Java [START_REF] Zumpt | The Stomoxyinae biting flies of the world. Taxonomy, biology, economic importance and control measures[END_REF]. Since both of our samples are from Thailand origin, they provide no information on the genetic diversity among populations of the species, but they reveal its evolutionary distinctiveness from other members of the genus Stomoxys.

The remaining Stomoxys species included in our study can be divided in two main clades. Clade B contains only African species and clade C groups both African and more geographically widespread species (Fig. 1). Clade B regroups S. inornatus, S. transvitattus, S. omega and S. pallidus sampled from different parts of the Ethiopian region. These species are ubiquitous and are restricted mainly to forest ecotones where they feed especially on wildlife fauna [START_REF] Mavoungou | Ecology of Stomoxys flies (Diptera: Muscidae) in Gabon. First survey in different ecological areas[END_REF]. This large clade appears to be divided in two main subclades associating S. inornatus and S. transvittatus on one side, and S. omega and S. pallidus, on the other side. S. inornatus and S. transvittatus have a diurnal activity whereas S. omega and S. pallidus have a crepuscular activity (Duvallet, unpublished data). These observations would fit well with the proposed phylogeny. However, the information content of ecological data for phylogeny is probably limited since activity levels in Stomoxys have been shown to be highly dependent upon geographical origin (primary forest, secondary forest and human-modified area) [START_REF] Mavoungou | Ecology of Stomoxys flies (Diptera: Muscidae) in Gabon. First survey in different ecological areas[END_REF], climate (temperature, humidity, and solar radiation level) [START_REF] Kangwagye | The seasonal incidence of biting flies (Diptera) in Rwenzori National Park and Kigezi Game Reserve, Uganda[END_REF][START_REF] Charlwood | The age-structure and biting behaviour of Stomoxys calcitrans (L.) (Diptera: Muscidae) from Manaus, Brazil[END_REF], trapping method, and physiological state of individuals (Simmond, 1944;[START_REF] Breque | Graphic models as intermediate steps to computerized simulations of stable fly populations[END_REF].

The third clade C can be divided into three subclades, two of them being exclusively African, and one being more geographically widespread. The first subclade groups S. niger niger and S. niger bilineatus, considered as two subspecies of S. niger according to [START_REF] Zumpt | The Stomoxyinae biting flies of the world. Taxonomy, biology, economic importance and control measures[END_REF], whereas the second subclade associates P. saegerae and S. varipes. Solving species boundaries between closely related species is notoriously difficult in these parasitic flies. Based on the tree topology, S. niger niger and S. niger bilineatus sequences form distinct phylogenetic clusters. The pairwise distance between S. niger niger and S. niger bilineatus based on COI sequences is 8.2%. The use of DNA sequences for species delimitation has been widely criticized [START_REF] Tautz | DNA points the way ahead in taxonomy[END_REF][START_REF] Desalle | The unholy trinity: taxonomy, species delimitation and DNA barcoding[END_REF][START_REF] Meier | DNA Barcoding and taxonomy in Diptera: a tale of high intraspecific variability and low identification success[END_REF]. The wide overlap between intraspecific and interspecific variability observed in Diptera COI sequences (0% to 15.5%) [START_REF] Meier | DNA Barcoding and taxonomy in Diptera: a tale of high intraspecific variability and low identification success[END_REF]) is especially problematic here, since the pairwise distance for these subspecies falls into this overlapping area. However, the genetic distance between S. niger niger and S. niger bilineatus based on the nuclear ITS2 is also quite large with 14.3 %. Moreover, our molecular dating analysis estimates the separation between the two subspecies around 16.3 Mya (Fig. 2). This is almost as old as the separation between S. inornatus and S. transvittatus (c.a. 16.4 Mya) and is comparable with the divergence observed between D. melanogaster and D. yakuba (c.a 16.2 Mya). It is also more ancient than the divergence between S. omega and S. pallidus (c.a. 12.4 Mya), and even more ancient than the split between S. varipes and P. saegerae (c.a. 14.2 Mya) (Fig. 2).

Morphologically, S. niger niger and S. niger bilineatus are distinguished only by the colour of their tibiae and tarsi. [START_REF] Mavoungou | Ecology of Stomoxys flies (Diptera: Muscidae) in Gabon. First survey in different ecological areas[END_REF] described S. niger bilineatus and S. niger niger as sympatric species, the former being abundant in savannas with abundant wild fauna, while the latter is more associated with anthropized area. [START_REF] Mihok | Phenology of Stomoxinae in a Kenya forest[END_REF] also pointed out differences in habitat affinities, sex ratio and activity patterns between these species. Given morphological observations, ecological isolation and molecular data, it appears very likely that S. niger niger and S. niger bilineatus belong to different species.

Consequently, we suggest raising them to full specific status and to use S. niger and S. bilineatus as species names.

Finally, the third subclade within clade C associates S. calcitrans, S. sitiens and S. bengalensis Picard. S. bengalensis occurs only in Asia, S. sitiens is found in both the African and Asian continents, while the livestock pest S. calcitrans has a worldwide distribution. S. calcitrans and S. sitiens appeared to be strictly human commensals since they have invariably been found in association with human activity, including inside buildings.

Molecular timescale and biogeography of the Stomoxys group

Our molecular estimate of the separation between Stomoxyini and Muscini falls into the Late Eocene epoch (c.a. 34.9 Mya) (Fig. 2). This divergence time is compatible with the record of the oldest Muscidae fossils in the Eocene and the Lower Miocene, 50-20 Mya [START_REF] Evenhuis | Catalogue of the fossil flies of the world (Insecta: Diptera)[END_REF]. The oldest Muscini fossil (about 20-15 Mya) was described from Dominican amber [START_REF] Pont | Three species of Muscidae (Diptera) from Dominican amber[END_REF], whereas no Stomoxyini fossil has been reported yet. The results presented here show that the radiation of the major Stomoxyini groups largely overlaps with the hypothetical time span (mid-Oligocene to mid-Miocene) during which temperate forests developed in the Northern Hemisphere [START_REF] Guo | Onset of Asian desertification by 22 Myr ago inferred from loess deposits in China[END_REF]. C). The collision between Eurasia and the Arabo-African plate was initiated by the Oligocene, but the Early Miocene period seems to have been key for the establishment of intercontinental pathways permitting faunal exchanges between the Afro-Arabian and Asia plates [START_REF] Bernor | A consideration of some major topics concerning Old World Miocene mammalian chronology, migrations and paleogeography[END_REF]. This biogeographical event fits well with our estimates of the split between clades B and C occurring within this period. Continental exchanges through migration events continued until the Late Miocene and could explained the occurrence of S. calcitrans and S. sitiens in both the African and Asian continents.

Migration until the Late Miocene may have been facilitated by a forest connection that existed between the African and Oriental regions at this epoch [START_REF] Moreau | Vicissitudes of the African biomes in the late Pliocene[END_REF].

Our study is a first step into the molecular phylogenetic analysis of the Stomoxys genus. Yet, our results point to the need for a taxonomic revision of the Stomoxyini tribe.

Since very few data on the biology and ecology of this group are currently available it would be necessary to advance the knowledge in these basic research areas. This is a prerequisite for better understanding how these species originated and diversified in a biogeographical context. Epochs are indicated following the geologic timescale 2004 of the Geological Society of America [START_REF] Gradstein | A Geologic Timescale 2004[END_REF]. Plio.: Pliocene; Pleist.: Pleistocene. Log-likelihood values of ML topologies inferred from each individual gene (CO1, CYTB and ITS2) and from their concatenation were computed using each of the four data matrices and then compared with the corresponding highest log-likelihood value (in bold). The difference in log-likelihood values derived from these crossed comparisons are indicated. An asterisk (*) signifies that the tested topology is significantly worse at the 5% level than the best ML topology inferred from the corresponding dataset. ns means that this difference is not statistically significant.

Figure legends

  using the following parameters: Minimum number of sequences for a conserved position = 38; Minimum number of sequences for a flanking position = 38; Maximum number of contiguous nonconserved positions = 8; Minimum length of a block = 5; Allowed gap positions = with half.

  753 nucleotide sites of which 228 (30.3%) are variable, and 198 (26.3 %) are parsimony informative. jModelTest identified the GTR+G+I model as the best-fitting model for the COI gene (lnL= -3610.24) based on the AIC. The gamma distribution shape parameter (alpha) was estimated to 0.2 which reveals strong among-site rate heterogeneity in this barcoding gene. The CytB alignment contained 587 sites of which 204 (34.8%) are variable and 161 (27.4 %) are informative. jModelTest also identified the GTR+G+I model as the best-fitting nucleotide substitution model for the CytB gene (lnL= -3121.46) with an estimated alpha of 0.62. Finally, the initial alignment of ITS2 sequences totalized 426 sites including indels. 131 sites corresponding to ambiguously aligned hypervariable regions were excluded from subsequent analyses by applying Gblocks. The remaining 295 sites contain 168 (57 %) variable sites of which 147 (50 %) are parsimony informative which makes this marker the most variable of our study. jModelTest identified the HKY+G model as the best-fitting model for the ITS2 marker (lnL= -2454.70). The gamma distribution shape parameter alpha was estimated to 0.43 in this ribosomal internal transcribed spacer.

  Haematobia and Haematobosca (PP = 1.0; BP ML = 71), and this clade appears as the sistergroup to the remaining Stomoxyini. The genus Stomoxys is rendered paraphyletic with strong support (PP = 1.0; BP ML = 90) due to the inclusion of Prostomoxys saegerae within it, as a sister-group to S. varipes (Fig.1). Strong support is obtained (PP = 1.00; BP ML = 97) for S.indicus as representing the sister-group to all remaining species within the Stomoxys group (Fig.1). Three major clades can be distinguished according to the Bayesian consensus tree within Stomoxyini. The first clade (clade A) is represented by S. indicus alone which is an Asian species. The second clade (clade B) is well supported (PP = 0.99 and BP ML = 63) and groups species from African origin (S. pallidus, S. omega, S. transvittatus, and S. inornatus).Finally, the third clade (clade C), only supported by the Bayesian analysis (PP = 0.92), regroups all the remaining species containing cosmopolitan, African, and Asian species. This clade can be further divided into three subclades. The first subclade associates the two subspecies S. niger niger and S. niger bilineatus (PP = 0.98), the second subclade associates P. saegerae and S. varipes (PP = 0.89), and the third subclade composed of S. calcitrans, S. sitiens, and S. bengalensis (PP = 1.0; BP ML = 56) with S. sitiens and S.

  The genus Stomoxys appears to diverge from other Stomoxyini genera in the Oligocene (c.a. 30 Mya). An Oriental origin of the Stomoxys genus lato sensu is tentatively suggested, owing to the strong support obtained for the basal branching of S. indicus (Clade A) and its absence from Africa. Under such a scenario, the common ancestor of the remaining Stomoxys members would have then split into two groups in the Early Miocene (c.a. 21 Mya), with one group being distributed in the Ethiopian region (Clade B), while the other colonized Asia (Clade

Fig. 1 .

 1 Fig. 1. Phylogenetic relationships among 33 dipteran species inferred from the concatenation of COI, CytB and ITS2 (1635 nucleotide sites). This phylogram is the 50% majority rule consensus tree obtained with Bayesian inference under the GTR+G model. Numbers at nodes indicate Posterior Probabilities (PP) / Maximum Likelihood Bootstrap Proportions (BP ML ). Black circles indicate nodes receiving maximum values of PP and BP and dashes mark nodes above the 50% level. The three main lineages identified within the Stomoxys group are labeled A, B, and C.

Fig. 2 .

 2 Fig. 2. Chronogram resulting from the relaxed molecular clock Bayesian analysis of the concatenation of COI, CytB and ITS2 (1635 nucleotide sites). Mean divergence times are indicated at nodes and node bars represent 95% Credibility Intervals. Scale is in Millions of years before present. Species names belonging to the genus Stomoxys have been abbreviated (S.). The Drosophila/Musca calibration constraint is shown in black. Tertiary

Figure 2

 2 Figure 2

Table 2 .

 2 Species sampling, source localities, and GenBank accession numbers.

	Species	Origin	COI	CytB	ITS2
	Stomoxys calcitrans	India	EU836073-	EU851301	NA
			EU836074
		Thailand	EU836075-	EU851303-	EU851200-
			EU836082	EU851308	EU851208
		Cameroun EU836070	EU851302	EU851199
	Stomoxys indicus	Thailand	EU836083-	EU851309-	EU851209-
			EU836092	EU851318	EU851218
	Stomoxys sitiens	Thailand	EU836131-	EU851357-	EU851251-
			EU836138	EU851364	EU851257
		Burkina	FJ386382	FJ386381	NA
	Stomoxys bengalensis	Faso Thailand	EU836060-	EU851281-	EU851191-
			EU836068	EU851288	EU851198
	Stomoxys niger niger	Gabon	EU836110-	EU851327-	EU851228-
			EU836115	EU85133	EU851237
	Stomoxys niger bilineatus	Gabon	EU836101-	EU851289-	EU851219-
			EU836108	EU851295	EU851227
	Stomoxys varipes	Ethiopia	EU836147-	EU851372-	EU851189-
			EU836148	EU851374	EU851190
	Stomoxys transvittatus	Gabon	EU836139-	EU851365-	EU851258-
			EU836145	EU851370	EU851265
	Stomoxys inornatus	Gabon	EU836093-	EU851319-	EU851179-
			EU836100	EU851326	EU851188
	Stomoxys omega	Gabon	EU836116-	EU851339-	EU851238-
			EU836121	EU851346	EU851244
	Stomoxys pallidus	Gabon	EU836122-	EU851347-	EU851245-
			EU836130	EU851356	EU851250
	Prostomoxys saegerae	Gabon	EU836055-	EU851275-	EU851266-
			EU836059	EU851280	EU851270
	Haematobosca croceicornis	Gabon	EU836053-	EU851273-	EU851271-
			EU836054	EU851274	EU851272
	Haematobia irritans	-	AY526195*	DQ029097*	DQ437515*
	Musca domestica	-	AF104622*	DQ657064*	Z28417*
	Drosophila yakuba	-	NC_001322* NC_001322* Z28416*
	Drosophila melanogaster	-	NC_001709* NC_001709* EU306667*
	Drosophila virilis	-	DQ471577*	AY646771*	Z28415*
	Cheilosia naruska	-	DQ417498*	NA	FJ028661*
	Cheilosia longula	-	FJ158631*	NA	FJ158631*
	Simosyrphus grandicornis	-	NC_008754* NC_008754* NA
	Dolichopus nubilus	-	AY958244*	AY958244*	NA
	Dolichopus excisus	-	AY958245*	AY958245*	NA

Table 3 .

 3 COI, CytB and ITS2 Primers used for amplifications and sequencing.

	Primers	Sequences (5'->3')	Reference
	C1-J-2813 (direct)	CAACATTTATTTTGATTTTTTGG	Simon et al.,1994;
	TL2-N-3014 (reverse)	TCCATTGCACTAATCTGCCATATTA	2006
	CB-J10933 (direct) CB-N11526 (reverse)	GTTTTACCTTGAGGACAAATATC TTCAACTGGTCGAGCTCCAATTCA	Simon et al., 1994
	ITS2A (direct) ITS2B (reverse)	TGTGAACTGCAGGACACAT TATGCTTAAATTCAGGGGGT	Sharpe et al., 2000

Table 4 .

 4 Statistical tests of topological congruence among the three genes using crossed SH tests[START_REF] Shimodaira | Multiple comparisons of log-likelihoods with applications to phylogenetic inference[END_REF].

		Datasets
	Topologies	CO1	CYTB	ITS2	Concatenation
	CO1	5207.93	42.75*	72.91*	70.32 ns
	CYTB	145.87*	4174.34	144.34*	297.63*
	ITS2	162.32*	200.42*	2452.65	323.39*
	Concatenation 12.89 ns	23.69 ns	3.86 ns	12265.83

Table 5 .

 5 Best fitting models and associated maximum likelihood parameters obtained for the gene partitions used in this study.

	Dataset	Sites	A	T	C	G	Model Best	Alpha Pinv
	COI	753	0.31	0.41	0.14	0.14	GTR+I+G 0.2	0.01
	First Pos.	251	0.46	0.49	0.04	0.01	TrN+G	0.45	-
	Second Pos.	251	0.29	0.32	0.12	0.27	GTR+G	0.16	-
	Third Pos.	251	0.21	0.41	0.24	0.14	HKY+I	-	0.80
	CytB	587	0.32	0.42	0.14	0.12	GTR+I+G 0.62	0.35
	First Pos.	196	0.30	0.32	0.16	0.21	TIM+G	0.24	-
	Second Pos.	196	0.20	0.46	0.29	0.11	HKY+I	-	0.68
	Third Pos.	195	0.45	0.48	0.06	0.01	HKY+G	0.77	-
	ITS2	295	0.36	0.37	0.10	0.17	HKY+G	0.43	-
	Concatenation	1635	0.33	0.11	0.12	0.44	GTR+I+G 0.95	0.44

The sequences marked with an asterisk (*) were obtained from GenBank, (NA) no available data.

Acknowledgements

We thank Jacques-François Mavoungou (Gabon), Roundthip Masmeatathip (Thailand), Alekaw Sinshaw Tegegne (Ethiopia) and the British Museum for providing some Stomoxys samples. We extend our thanks to Gwenaëlle Mondor (CBGP) for technical assistance. We also thank two anonymous reviewers for constructive comments. This study was partially supported by Agence Nationale de la Recherche through the ANR-Biodiv grant IFORA (Iles Forestières Africaines) and by the CNS-Genoscope "Barcoding insects for identification" sequencing grant to Michel Veuille. During the course of this research, N. D.-A.

was partially supported by a grant from Société Entomologique de France. This is contribution ISEM 2011-XXX of the Institut des Sciences de l'Evolution (UMR 5554-CNRS). Table 1. Stomoxys species distributions as mentioned by Zumpt (1973). Species Distribution S. bengalensis Picard India to Java S. boueti Roubaud Benin (ex-Dahomey), Congo S. calcitrans (Linnaeus) Cosmopolitan S. indicus Picard Oriental region and neighbouring Paleartic territories S. inornatus Grünberg Tropical Africa S. luteolus Villeneuve Central and East Africa S. niger Macquart S. niger niger Ethiopian and Madagascan regions S. niger bilineatus Ethiopian and Madagascan regions S. ochrosoma Speiser Central and East Africa S. omega Newstead Ethiopian Region S. pallidus Roubaud Tropical Africa S. pullus Austen India S. sitiens Rondani Ethiopian and Oriental regions, in Egypt it reaches the Palearctic region S. stigma Van Emden Uganda -Congo S. taeniatus Bigot Ethiopian region S. transvittatus Villeneuve Southern and Central Africa S. uruma Shinonaga & Kano Oriental region S. varipes Bezzi East and Central Africa southward to Rhodesia S. xanthomelas Roubaud Congo -Tanzania -Uganda


Julien Fleureau

Amar Kachenoura
email: amar.kachenoura@univ-rennes1.fr

Laurent Albera

Jean-Claude Nunes

Lotfi Senhadji

Multivariate empirical mode decomposition and application to multichannel filtering

Keywords: Mono-and multivariate empirical mode decomposition, filter bank structure, electroencephalography data analysis

Empirical Mode Decomposition (EMD) is an emerging topic in signal processing research, applied in various practical fields due in particular to its data-driven filter bank properties. In this paper, a novel EMD approach called X-EMD (eXtended-EMD) is proposed, which allows for a straightforward decomposition of mono-and multivariate signals without any change in the core of the algorithm. Qualitative results illustrate the good behavior of the proposed algorithm whatever the signal dimension is. Moreover, a comparative study of X-EMD with classical mono-and multivariate methods is presented and shows its competitiveness. Besides, we show that X-EMD extends the filter bank properties enjoyed by monovariate EMD to the case of multivariate EMD. Finally, a practical application on multi-channel sleep recording is presented.

Introduction

The Empirical Mode Decomposition (EMD) is raising great interest since its first appearance in the nineteens [1]. Indeed, it received wide attention in various fields such as biomedical engineering [2, 3], space research [START_REF] Coughlin | 11-Year solar cycle in the stratosphere extracted by the empirical mode decomposition method[END_REF], hydrology [START_REF] Huang | Analysis of daily river flow fluctuations using Empirical Mode Decomposition and arbitrary order Hilbert spectral analysis[END_REF], synthetic aperture radar [START_REF] Zhou | Narrow-Band Interference Suppression for SAR Based on Complex Empirical Mode Decomposition[END_REF], speech enhancement [START_REF] Chatlani | EMD-based Noise Estimation and Tracking (ENET) with application to speech enhancement[END_REF], watermarking [START_REF] Bi | Robust Image Watermarking Based on Multiband Wavelets and Empirical Mode Decomposition[END_REF], etc. EMD aims at decomposing sequentially a given signal into the sum of Amplitude and Frequency Modulated (AM/FM) zero-mean oscillatory signals, called Intrinsic Mode Functions (IMFs), plus a non-zero mean low-degree polynomial remainder. Each IMF is computed by using an iterative procedure, named sifting process, applied to the residual multicomponent signal. One sifting iteration consists: i) in computing the local mean of the residual signal, and ii) in subtracting it from the residue. This process is repeated until convergence to a designated IMF. Generally, the local mean is calculated as the half sum of the upper and the lower envelopes, obtained by interpolating between local minima and maxima, respectively, using cubic splines [1].

Since the original work on EMD [1], several studies have been presented. Some ones aim at giving a theoretical framework to monovariate EMD [START_REF] Delechelle | Empirical mode decomposition: an analytical approach for sifting process[END_REF][START_REF] Kopsinis | Investigation and Performance Enhancement of the Empirical Mode Decomposition Method Based on a Heuristic Search Optimization Approach[END_REF][START_REF] Hawley | Some Properties of an Empirical Mode Type Signal Decomposition Algorithm[END_REF] in order to make easier its analysis. Other ones underline the behavior of original EMD [1] on simulated data such as its data-driven filter bank structure in the case of broadband noise [START_REF] Wu | A study of the characteristics of white noise using the empirical mode decomposition method[END_REF][START_REF] Flandrin | Empirical Mode Decomposition as a Filter Bank[END_REF]. Eventually, EMD procedure devoted to bivariate [START_REF] Tanaka | Complex empirical mode decomposition[END_REF][START_REF] Bin Altaf | Rotation invariant complex empirical mode decomposition[END_REF][START_REF] Rilling | Bivariate empirical mode decomposition[END_REF], trivariate [START_REF] Ur Rehman | Empirical Mode Decomposition for Trivariate Signals[END_REF] and, recently, multivariate signals [START_REF] Ur Rehman | Multivariate empirical mode decomposition[END_REF][START_REF] Fleureau | 3A-EMD: a generalized approach for monovariate and multivariate EMD[END_REF][START_REF] Fleureau | Turning Tangent Empirical Mode Decomposition: a unified framework for mono-and multivariate signals[END_REF] have been proposed. The first bivariate approach in [START_REF] Tanaka | Complex empirical mode decomposition[END_REF] uses the real parts of analytic and anti-analytic components of the bivariate signal which leads to two sets of complex-valued IMFs corresponding to the positive and negative frequency components of the signal, respectively. The second bivariate technique presented in [START_REF] Bin Altaf | Rotation invariant complex empirical mode decomposition[END_REF] defines the extrema of a complex signal as points where the angle of the first derivative of the signal vanishes.

In [START_REF] Rilling | Bivariate empirical mode decomposition[END_REF], two bivariate algorithms based on more geometrical considerations are developed. These two last approaches compute a local mean by averaging the monodimensional EMD's of projections of the original bivariate signal onto different angular planes. Rehman et al. proposed a trivariate approach [START_REF] Ur Rehman | Empirical Mode Decomposition for Trivariate Signals[END_REF] which is based on the computation of the extrema by projecting the input signal in multiple directions in 3D spaces via a quaternion rotation framework. The multivariate extension of EMD, developed in [START_REF] Ur Rehman | Multivariate empirical mode decomposition[END_REF], can be seen as a generalization of the concept employed in the Rilling's bivariate [START_REF] Rilling | Bivariate empirical mode decomposition[END_REF] and the Rehman's trivariate [START_REF] Ur Rehman | Empirical Mode Decomposition for Trivariate Signals[END_REF] In recent works, [START_REF] Fleureau | 3A-EMD: a generalized approach for monovariate and multivariate EMD[END_REF][START_REF] Fleureau | Turning Tangent Empirical Mode Decomposition: a unified framework for mono-and multivariate signals[END_REF], we proposed an EMD algorithm, called 2T-EMD (Turning Tangent EMD), ables to process mono-and multivariate signals whatever the signal dimension D is, without any change in the algorithm.

In addition, 2T-EMD algorithm seems to require a smaller amount of calculations than classical EMD algorithms. To do so, the signal mean trend is redefined as the signal obtained by interpolating the barycenters of particular oscillations, called elementary oscillations (see section II [START_REF] Fleureau | 3A-EMD: a generalized approach for monovariate and multivariate EMD[END_REF][START_REF] Fleureau | Turning Tangent Empirical Mode Decomposition: a unified framework for mono-and multivariate signals[END_REF]). More pre-cisely, 2T-EMD algorithm is based on two main steps: i) identification of all elementary oscillations and ii) computation of the barycenters of each associated elementary oscillations and interpolation between all these barycenters to obtain the signal mean trend. Recall that, in [START_REF] Fleureau | 3A-EMD: a generalized approach for monovariate and multivariate EMD[END_REF][START_REF] Fleureau | Turning Tangent Empirical Mode Decomposition: a unified framework for mono-and multivariate signals[END_REF], an elementary oscillation of a given function s with values in Ê D (D ≥ 1) is considered as a piece of s defined between two consecutive oscillation extrema of s.

In this paper, we propose to compute the local mean trend without calculation of the oscillation barycenters (as proposed in 2T-EMD), giving rise to a new multivariate EMD algorithm, named X-EMD (eXtended EMD). Indeed, in X-EMD algorithm (as described in section 2), the local mean is directly computed by averaging two envelopes generated directly from oscillation extrema. In other words, the exploitation of oscillation extrema is clearly different between X-EMD and 2T-EMD algorithms and makes X-EMD algorithm easier both to implement and to use. The numerical complexity of the proposed approach is then analyzed in section 3 and compared with the complexity of four existing EMD methods. A performance study is proposed in section 4 on synthetic mono-, bi-, and trivariate signals showing: i) the competitive behavior of X-EMD versus existing methods, and ii) its ability to process any output signal dimension. Finally, data-driven spectrum analyzer properties satisfied by X-EMD are highlighted and a biomedical application is considered in section 5. Indeed, empirical filter bank structure enjoyed by monovariate EMD [START_REF] Flandrin | Empirical Mode Decomposition as a Filter Bank[END_REF] is shown to be preserved by X-EMD in multivariate context. Such results is used to analyze a quadrivariate ElectroEncephalo-Graphic (EEG) sleep recording.

An eXtended EMD approach: X-EMD

The proposed X-EMD method aims at providing a unique tool able to process mono-and multivariate signals without any modification. The notion of oscillation extrema, originally presented in [START_REF] Fleureau | 3A-EMD: a generalized approach for monovariate and multivariate EMD[END_REF][START_REF] Fleureau | Turning Tangent Empirical Mode Decomposition: a unified framework for mono-and multivariate signals[END_REF], is exploited to redefine the local mean operator. However, before going further, let's briefly recall the definition of oscillation extrema.

Oscillation extremum. In [1], the local mean is given by taking average of the upper and the lower envelopes, which are computed by interpolating between the local maxima and minima, respectively. However, such points cannot be used for higher output signal dimensions, since they are not defined directly in multivariate contexts. To deal with this problem, we propose [START_REF] Fleureau | 3A-EMD: a generalized approach for monovariate and multivariate EMD[END_REF][START_REF] Fleureau | Turning Tangent Empirical Mode Decomposition: a unified framework for mono-and multivariate signals[END_REF] a new and general way to compute extrema in multidimensional output space Ê D .

Let s be a function defined from Ê to Ê D and differentiable with a continuous first derivative, say in class C 1 . Let define the function tangent vector to s by:

T s : t -→ [1, ds dt (t)] T (1)
Based on the following definition:

α s : t -→ lim h→0 T s (t -h), T s (t + h) (2)
where • , • denotes the Euclidean inner product of Ê D+1 , the function α s allows us to identify direction changes in function s. Due to the continuity of the inner product, it is noteworthy that function α s can be rewritten as:

∀t ∈ Ê, α s (t) = lim h→0 T s (t -h), lim h→0 T s (t + h) (3)
Now, since the considered function s is in C 1 , α s takes the following form:

∀t ∈ Ê, α s (t) = T s (t) 2 = 1+ ds dt (t) 2 (4)
where . abusively denotes the Euclidean norm of both Ê D+1 and Ê D . So, oscillation extrema may be defined as the local minima of the following β s function:

β s : t -→ β s (t) = ds dt (t) 2 (5)

Numerical complexity analysis

This section aims at relating the X-EMD algorithm and some existing EMD methods from a computational complexity point of view. Numerical complexity is calculated in terms of number of floating point operations (flops). A flop corresponds to a multiplication followed by an addition. But in practice, only multiplications are counted since, in general, there are about as many (and slightly more) multiplications as additions. The first considered EMD method is Huang's solution [1] and will be named Huang in the sequel.

The second studied standard EMD approach is the second geometrical bivariate solution proposed in [START_REF] Rilling | Bivariate empirical mode decomposition[END_REF] and called Rilling hereafter. As far as the cases of multivariate signals are concerned, two methods, namely 2T-EMD algorithm [START_REF] Fleureau | 3A-EMD: a generalized approach for monovariate and multivariate EMD[END_REF][START_REF] Fleureau | Turning Tangent Empirical Mode Decomposition: a unified framework for mono-and multivariate signals[END_REF] and Rehman's method [START_REF] Ur Rehman | Multivariate empirical mode decomposition[END_REF], are considered. )

Huang [1] 1 18L + 15M H (d n,k ) Rilling [16] 2 L(11P + 2) + 15 P/2 p=1 M R (d n,k , p) Rehman [18] AE * LP (2D + 18) + 15 P p=1 M RM (d n,k , p) 2T-EMD [19, 20] AE * D(19L + 16M 2T (d n,k )) + M 2T (d n,k ) X-EMD AE * D(19L + 15M X (d n,k ))
Table 1: Computational complexity for one sifting iteration of X-EMD and four existing EMD methods. L and D represent the data length and the data dimension, respectively.

Where EMD. For given values of L and D, the computational complexity, F (d n,k+1 ), of one sifting iteration necessary to obtain d n,k+1 from d n,k is given in table 1 for Huang, Rilling, Rehman, 2T-EMD and X-EMD methods. Note that these results were obtained for a standard tridiagonal implementation of the spline interpolation and a signal from Ê to Ê D . Finally, the total computational cost of the five methods is obtained by summing the elementary complexities given in table 1 over the number of iterations, K n , and the number of IMFs, N, as shown below:

M H (d n,k ), M R (d n,k , p), M RM (d n,k , p), M 2T (d n,k
C(s) = N n=1 Kn k=1 F (d n,k ) (6)

Performance evaluation of X-EMD

The objective of this section is to compare the performance of X-EMD algorithm with classical approaches and to illustrate its ability to process multivariate signals whatever their output dimension. Qualitative multivariate decompositions. This section is dedicated to qualitatively illustrate the behavior of X-EMD in various multivariate contexts. Thus, the decompositions of the monovariate signal s 12 , the bivariate signal s 21 , and the trivariate signal s 31 computed on the time interval [0; 1] are respectively depicted in figures 1, 2 and 3. Figure 1 presents the exact (dashed line) and estimated (plain line) IMFs. The correct linear shape of the residue illustrates the good behavior of the X-EMD algorithm in a monovariate case. Regarding the bivariate signal s 21 depicted in figure 2, the two rotated planar IMFs are well extracted, with a low residue. As far as the trivariate signal is concerned, the figure 3 Comparative study. This section adopts a more quantitative point of view and is dedicated to the comparison of X-EMD performance with those obtained by Huang's algorithm [1] in monovariate case, Rilling's method [START_REF] Rilling | Bivariate empirical mode decomposition[END_REF] in a bivariate context, Rehman's approach [START_REF] Ur Rehman | Multivariate empirical mode decomposition[END_REF] for trivariate signals, and 2T-EMD [START_REF] Fleureau | 3A-EMD: a generalized approach for monovariate and multivariate EMD[END_REF][START_REF] Fleureau | Turning Tangent Empirical Mode Decomposition: a unified framework for mono-and multivariate signals[END_REF] methods for all tested signals. More precisely, the second geometric approach proposed in [START_REF] Rilling | Bivariate empirical mode decomposition[END_REF] is used as reference for bivariate decompositions with P = 8 projection planes. For Rehman's method we take P = 2D projection directions on the associated hypersphere. Note that, the sifting process termination criterion and the border effect management are identical for all compared methods. This comparison is achieved in terms of performance using criteria e [0,1] (s) and e [-1,0]∪[1,2] (s), but also in terms of numerical complexity using criterion C(s). Results are provided in table 2.

s 3i = n d (
It appears that for no border effects (I = [0, 1]), the monovariate signal s 11 with saddle points is better decomposed by our approach in comparison to the Huang's method. In the bivariate case, Rilling's algorithm seems to be slightly more efficient on the rotating signal s 22 but X-EMD seems to offer a better performance on the signal s 21 made of rotated planar components. Regarding the trivariate signal s 31 , X-EMD method leads to better result than Rehman's algorithm. As far as the border effects are concerned, X-EMD seems to have quasi-identical performance than Huang's method, and provides a more accurate management of border effects in comparison to the Rilling's and Rehman's methods. Concerning the computational complexity, table 2 clearly shows that for D = 1 the computational complexities of X-EMD and Huang's method are quasi-equivalent. On the contrary, for D = 2 and D = 3, X-EMD requires less sifting iterations and less computational operations than Rilling's and Rehman's methods. When focusing on the comparison between X-EMD method and 2T-EMD method (say the two methods which exploit, in different manner, the oscillation extrema), table 2 exhibits that the behavior of the two methods, using criteria e [0,1] (s) and e [-1,0]∪[1,2] (s), are globally equivalent. More finely, X-EMD seems to provide a more accurate management of border effect for AM-FM signal s 12 , where as 2T-EMD gives better border effects management in case of signal s 11 with saddle points. However, in accordance with the formulas presented in table 1, the proposed X-EMD generally requires less computational operations than 2T-EMD method.

Filter bank structure and application to multichannel sleep recording

One important property of the classical monovariate EMD [1] especially highlighted in [START_REF] Wu | A study of the characteristics of white noise using the empirical mode decomposition method[END_REF] is its empirical filter bank structure observed when EMD is applied on broadband noises. Filter bank structure. The goal of this section is to show how this property is preserved by X-EMD for higher output signal dimensions. Nevertheless, the notion of frequential content for multivariate signals may be ambiguous.

As suggested in [START_REF] Flandrin | Empirical Mode Decomposition as a Filter Bank[END_REF], one alternative way to deal with IMFs frequency in a monovariate context is to consider their number of zero-crossings or alternatively their number of extrema, N ext (as those quantities differ at most by one in the definition of an IMF). In the same sense, we propose here to evaluate the frequential content of a mono-or multivariate IMF by considering its number of oscillation extrema. It is worth to point out that all the results presented hereafter do not take into account the first mode, d n,k (k=1), because as shown in [START_REF] Wu | A study of the characteristics of white noise using the empirical mode decomposition method[END_REF][START_REF] Flandrin | Empirical Mode Decomposition as a Filter Bank[END_REF] the filter associated to the first mode is essentially highpass. In a first simulation, 500 uncorrelated mono-, bi-and trivariate white noises (1024 samples, with normal distribution) are generated and decomposed using X-EMD. For each IMF d n,k (k>1) and for each signal, N ext is computed and a linear regression model is then fitted on the log 2 averaged number of oscillation extrema (log 2 (N ext )) as depicted in figure 4. The linear regression quality shows how N ext exponentially decreases with the number of IMFs, and suggests a filter bank structure whatever the considered dimension is. This bank appears close to a dyadic filter bank for a monovariate decomposition (-0.65 slope in log 2 coordinates) and the exponential decrease seems to be lower with the dimensionality increase (-0.17 slope for the hexavariate D = 6 decomposition). This result strongly suggests that a higher number of IMFs could be expected for decompositions with higher dimensionality.

In a second simulation, the influence of spatial correlation between signal channels on the filter bank structure is studied. To do so, the previous simulation is repeated in a bivariate context for different values of the spatial covariance factor ρ. More precisely, for 11 values of ρ regularly varying from 0 (uncorrelated channels) to 1 (the two channels are identical), 500 bivariate white noises (1024 samples, bidimensional-normal distribution with unitary variance and ρ covariance) are generated. For each signal, the averaged number of oscillation extrema, N ext , is computed and fitted in log 2 coordinates. Fitting errors and associated slope values are depicted in figure 5. The low fitting error suggests that the spatial covariance doesn't influence the filter bank structure. Regarding the slope, it clearly seems to vary with ρ values. For ρ = 0, an uncorrelated bivariate white noise is decomposed and the slope value (around -0.48) is close to the one found in the first bivariate simulation. For ρ = 1, the two signal channels are identical, and the decomposition is "equivalent" to a monovariate uncorrelated decomposition.

The slope value (around -0.62) is thus close to the one found in the first monovariate simulation.

Application to multichannel sleep recording. The filter bank structure property enjoyed by X-EMD strongly suggests the use of the proposed algorithm in order to analyze the frequential content of real multivariate signals.

One example of application dealing with multi-channel sleep recording is proposed here. Therefore, an ambulatory sleep staging system [START_REF] Porée | Blind Source Separation for Ambulatory Sleep Recording[END_REF] takes advantage of the mutual information contained in the four channels, preserves the frequential contents on each channel. Such results may be very interesting in practical contexts but they require a priori regarding the processed signal in order to select the suited IMFs. More particularly, an analysis of the broadband noise present in the signal to process should be performed with X-EMD beforehand. In the biomedical context considered here the broadband noise is essentially due to the neuronal electrical activity.

It is noteworthy that if no broadband noise is present in the considered application, it could be artificially added to the data.

Conclusion and perspectives

A new algorithm called X-EMD based on a redefinition of the local mean operator is proposed. The obtained results show that, under certain assumptions on the processed signal, this alternative definition enables to decompose both mono-and multivariate signals without any modification in the core of the algorithm. This last point is the main difference regarding the existing approaches of the literature. The comparative study also suggests that X-EMD seems to offer competitive performance whereas the computational complexity analysis suggests that no loss in terms of numerical performance should be observed in comparison with the existing EMD approaches. Moreover filter bank properties of the classical EMD seem to be preserved and applicable to multivariate signal denoising. A restricted scope of application mainly due to the use of the first derivative in the algorithm remains the most important limitation of our approach. However, as shown in the application of X-EMD in quadrivariate sleep recording, those limitations should be not too much restrictive in practical fields. More simulated and real data decompositions will be performed in future works.

Acknowledgment

We which acknowledge all the reviewers and the associate editor for their valuable comments. This work is supported by the French Government under an ANR Contract, namely mv-EMD (BLAN07-0314-02).

[1] N. Huang, Z. Shen, S. Long, M. Wu, H. Shih, Q. Zheng, N.-C. Yen, C. Tung, and H. Liu, "The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis," Proc. R.

Soc. London, Ser. A (1998).

[2] L. Hualou, L. Qiu-Hua, and J. Chen, "Application of the empirical mode decomposition to the analysis of esophageal manometric data in gastroe-

  EMD. The mean envelope is obtained by averaging multiple envelopes, generated by taking projections of the multivariate signal along multiple directions on hyperspheres. Nevertheless, none of these techniques is a straightforward extension of the original monovariate EMD concept [1] to multivariate signals, say signals defined from Ê to Ê D with D > 1. It is probably due to the difficulty of interpreting the notion of extremum when multivariate signals are considered.

For a given EMDD

  algorithm, let N, K n , and L be the number of extracted IMFs, the number of sifting iterations performed to extract the n-th IMF, and the data length, respectively. Now, denote by d n,k the n-th IMF computed at the k-th iteration of the sifting process, by: i) M H (d n,k ) the number of extrema detected in d n,k by Huang [1], ii) M R (d n,k , p) the number of extrema detected in the p-th projection of d n,k using Rilling [16] when P projection planes are used, iii) M RM (d n,k , p) the number of extrema detected in the p-th projection of d n,k by Rehman [18] when P projection directions are used, iv) M 2T (d n,k ) the number of oscillation barycenters detected in d n,k by 2T-EMD [19, 20], and V) M X (d n,k ) the number of oscillation extrema detected in d n,k by X-Method Numerical complexity F (d n,k+1

  ) and M X (d n,k ) denote the number of extrema detected in d n,k by Huang, the number of extrema detected in the p-th projection of d n,k by Rilling (when P projection planes are used), the number of extrema detected in the p-th projection of d n,k by Rehman (when P projection directions are used), the number of barycenters detected in d n,k by 2T-EMD, and the number of oscillation extrema detected in d n,k by X-EMD, respectively.

  3i) n , defined on the time interval T = [-1; 2] and sampled at f s = 1kHz, are tested. The mono-and bivariate signals are given by (succes-sive IMFs are in curly brackets): s 11 (t) = d

=

  {0.5e 0.23(1+t) sin(20πt)} + {e 0.23(1+t) sin(10πt))}s 12 (t) = d

=

  {2 sin(50πt)} + {e 0.23(1+t) cos(π(1.41t + 9.90) 2 )} +{3 cos(π(t + 6) 2 )} + 3ts 21 (t) = d

  sin(80πt + 1.5)} + {e -iπ 4 sin(2π[2.5(1 + t)] 2 )} s 22 (t) = d

=

  {0.3 cos(56πt) + 0.8i sin(54πt)} + {1.7 cos(20πt) + 3.2i sin(22πt)} + {e 0.23t-0.46 [2.5 cos(π(0.82t + 4.49) 2 ) + 1.1i sin(π(0.82t + 4.49) 2 )]} and the trivariate signal is defined by s 31 = d ) = 2.5 sin(100πt), 2 cos(100πt), e 0.11t+0.70 sin(100πt) T d (31) 2 (t) = [5e 0.14(1+t) cos(π(t + 13) 2 ), 2e 0.14(1+t) sin(π(t + 13) 2 ), 0.5 sin(28πt + 1.2)] T d (31) 3 (t) = [4 sin(12πt), cos(12πt), 6 sin(14πt + 1.5)] T From a general point of view, the proposed synthetic signals are the sum of multiple generic AM-FM components. More precisely, the monovariate signal s 11 is the sum of two exponential AM components containing saddle points. The signal s 12 is a more general AM-FM signal with a linear residue. Regarding the bivariate signals, s 21 is composed of two orthogonal AM-FM planar components whereas s 22 is the sum of different AM-FM elliptic rotat-ing components. Finally the trivariate signal s 31 is the sum of three general AM-FM components.

  (a) represents the original signal (dark line) and the associated local mean (gray line) on the restricted [0.5; 0.6] time interval and the figure 3 (b) displays the exact (dashed line) and the estimated (plain line) IMFs of the signal projected on the three main axis, namely X, Y and Z. The local mean seems to nicely go through the original signal, and the comparison between the estimated and exact projections of the IMFs shows the ability of the X-EMD method to perform trivariate decomposition without any change in the core of the algorithm.

Figure 1 : 21 Figure 2

 1212 Figure 1: X-EMD monovariate decomposition of the signal s 12 . Dashed line: real IMF / plain line: estimated IMF.

  Estimated IMF's projected on the three main axis : Signal s 31

Figure 3 :

 3 Figure 3: X-EMD trivariate decomposition of the signal s 31 on [0; 1]: (a) original signal (dark) and its X-EMD local mean (gray), and (b) exact (dashed line) and estimated (plain line) IMFs projected on the three main axis X, Y and Z.

Table 2 :

 2 Comparative study of X-EMD versus Huang (1D), Rilling (2D), Rehman (2D) and 2T-EMD (1D, 2D and 3D) reference methods.

16 Figure 4 :

 164 Figure 4: X-EMD filter bank structure: i) top (first and second lines): number of oscillation extrema, N ext , with boxplots as a function of the index of IMFs for uncorrelated white noises of dimensions : D = 1, D = 2 and D = 3,D = 4, D = 5 and D = 6, ii) bottom (third and fourth lines): averaged number of oscillation extrema, N ext (cross) as a function of the index of IMFs in log 2 coordinates and the associated linear regression model (continuous line).

Figure 6 :

 6 Figure 6: Application of X-EMD to multichannel sleep recording.

  Before applying X-EMD on test signals, it is prudent to give some details about the set of signals that X-EMD can decompose successfully. As previously mentioned, the considered signal has to be in class C 1 or at least, in the case of irregular signals, an appropriate numerical estimation of the derivative has to be proposed. For instance, the first order derivative may be computed by means of where dn denotes the estimate of the exact n-th IMF d n of signal s. More precisely, if we consider the interval, I = [0, 1], where border effects should be low, equation[START_REF] Chatlani | EMD-based Noise Estimation and Tracking (ENET) with application to speech enhancement[END_REF] allows us to compute a performance independent of any border effects. On the contrary, if I = [-1, 0] ∪ [1, 2] and for the same border management of all algorithms, equation (7) allows us to evaluate the ability of a given algorithm to minimize border effects. Indeed, such effects are often critical, especially, in real life data and their minimization should facilitate the practical exploitation of the IMFs. The second criterion, C(s), evaluates the numerical complexity of each EMD algorithms, with respect to the formulas listed in table1.Signal selection. To evaluate the behavior of X-EMD algorithm, two monovariate (s 11 , s 12 ) signals of the form s 1i = n d

				(1i) n , two bivariate signals (s 21 ,
	s 22 ) of the form s 2i =	n d	(2i) n	and a trivariate (s 31 ) signal of the form
				N
				e I (d n )	(7)
				n=1

a centered finite difference when forward and backward finite differences have the same signs and set to 0 otherwise. In addition, if {d n } 1≤n≤N represents the theoretical set of IMFs composing the signal to analyze, it is obvious that any d n with a piecewise constant function β dn (5) is not visible by X-EMD (a piecewise constant function has no local minimum). It concerns, for instance, monovariate signals with piecewise linear IMFs, and bivariate signals with purely circular rotating IMFs. Performance criteria. Two criteria are used to evaluate the X-EMD algorithm and compare it to the existing EMD methods. The first criterion allows us to evaluate the efficiency of EMD algorithms to accurately extract one or all expected IMFs. Given I a subinterval of T = [-1, 2], let's define the following quadratic errors: e I (d n ) = I ||d n (t) -dn (t)|| 2 dt I ||d n (t)|| 2 dt and e I (s) =

could practically take advantage of the X-EMD algorithm. A quadrivariate application of X-EMD to the quadrivariate observation O is therefore performed on four IMFs in a supervised way. The sum of the first four obtained

T and the as-

The channels of the EMG signal estimated by using X-EMD (line 2) and the EMG signals obtained with a highpass filtering (line 4) are, visually, quasi-equivalent. Indeed, the correlation coefficients computed between each channel are comprised between 0.8095 and 0.8320. This high correlation is also observed on the associated low-frequency residual parts (line 3 and 5).

Those last results show that the direct quadrivariate X-EMD method, which