Emmanuel Corse
email: emmanuel.corse@gmail.com

Emese Megl Ecz

| Ga€ It Archambaud

| Morgane Ardisson

Jean-Franc ßois Martin

| Christelle Tougard

| R Emi Chappaz

Vincent Dubut
email: vincent.dubut@imbe.fr

A from-benchtop-to-desktop workflow for validating HTS data and for taxonomic identification in diet metabarcoding studies

Keywords: cytochrome c oxidase I, diet studies, HTS data filtering, metabarcoding, taxonomic assignment

The main objective of this work was to develop and validate a robust and reliable "from-benchtop-to-desktop" metabarcoding workflow to investigate the diet of invertebrate-eaters. We applied our workflow to faecal DNA samples of an invertebrate-eating fish species. A fragment of the cytochrome c oxidase I (COI) gene was amplified by combining two minibarcoding primer sets to maximize the taxonomic coverage. Amplicons were sequenced by an Illumina MiSeq platform. We developed a filtering approach based on a series of nonarbitrary thresholds established from control samples and from molecular replicates to address the elimination of crosscontamination, PCR/sequencing errors and mistagging artefacts. This resulted in a conservative and informative metabarcoding data set. We developed a taxonomic assignment procedure that combines different approaches and that allowed the identification of ~75% of invertebrate COI variants to the species level. Moreover, based on the diversity of the variants, we introduced a semiquantitative statistic in our diet study, the minimum number of individuals, which is based on the number of distinct variants in each sample. The metabarcoding approach described in this article may guide future diet studies that aim to produce robust data sets associated with a fine and accurate identification of prey items.

| INTRODUCTION

In ecology and conservation, reliable diet data sets are critical to the understanding of prey/habitat relationships and feeding habitats. In this perspective, the use of DNA-based approaches in trophic ecology has grown during the last years (e.g., [START_REF] Razgour | High-throughput sequencing offers insight into mechanisms of resource partitioning in cryptic bat species[END_REF][START_REF] Soininen | Highly overlapping winter diet in two sympatric lemming species revealed by DNA metabarcoding[END_REF], especially due to the advent of high-throughput sequencing (HTS), leading to the development of metabarcoding [START_REF] Taberlet | Towards next-generation biodiversity assessment using DNA metabarcoding[END_REF]. However, diet metabarcoding approaches face four main challenges: (i) amplification bias related to the degradation of DNA [START_REF] Sint | Optimizing methods for PCR-based analysis of predation[END_REF]; (ii) taxonomic coverage of primers [START_REF] Gibson | Simultaneous assessment of the macrobiome and microbiome in a bulk sample of tropical arthropods through DNA metasystematics[END_REF]; (iii) taxonomic identification and resolution of DNA barcode sequences [START_REF] Richardson | Evaluating and optimizing the performance of software commonly used for the taxonomic classification of DNA metabarcoding sequence data[END_REF]; and (iv) filtering of HTS data to eliminate artefacts (e.g., PCR/ sequencing errors, mistagging, contamination) that produce false positives and constitute low-frequency noise (LFN; sensu De [START_REF] De Barba | DNA metabarcoding multiplexing and validation of data accuracy for diet assessment: Application to omnivorous diet[END_REF]. In diet studies, metabarcoding primers should therefore target short regions (i.e., <300 bp) of multicopy DNA to tackle the degradation of DNA [START_REF] Pompanon | Who is eating what: Diet assessment using next generation sequencing[END_REF]. Moreover, binding sites of primers should be sufficiently conserved to minimize biases in taxonomic coverage [START_REF] Clarke | Environmental metabarcodes for insects: In silico PCR reveals potential for taxonomic bias[END_REF][START_REF] Deagle | DNA metabarcoding and the cytochrome c oxidase subunit I marker: Not a perfect match[END_REF]. Alternatively, taxonomic coverage can be improved by amplifying several loci [START_REF] De Barba | DNA metabarcoding multiplexing and validation of data accuracy for diet assessment: Application to omnivorous diet[END_REF] or using several sets of primers that target the same locus [START_REF] Gibson | Simultaneous assessment of the macrobiome and microbiome in a bulk sample of tropical arthropods through DNA metasystematics[END_REF]. In both cases, the choice of the PCR primers is critical for optimizing detection of all prey in faeces.

As for the taxonomic identification, a considerable trade-off between accuracy and sensitivity should be considered when selecting an assignment procedure [START_REF] Richardson | Evaluating and optimizing the performance of software commonly used for the taxonomic classification of DNA metabarcoding sequence data[END_REF]. Moreover, both the confidence and the resolution of taxonomic classifiers are highly dependent on the richness of reference sequence databases of the targeted loci [START_REF] Gibson | Simultaneous assessment of the macrobiome and microbiome in a bulk sample of tropical arthropods through DNA metasystematics[END_REF][START_REF] Porter | Rapid and accurate taxonomic classification of insect (class Insecta) cytochrome c oxidase subunit 1 (COI) DNA barcode sequences using a na€ ıve Bayesian classifier[END_REF]. Additionally, producing robust data sets is critical for conducting reliable ecological studies. However, most of the current clustering-based methods for filtering HTS need specific (and partly arbitrary) parameterization and often overestimate the real number of taxa in samples [START_REF] Brown | Divergence thresholds and divergent biodiversity estimates: Can metabarcoding reliably describe zooplankton communities?[END_REF][START_REF] Clare | The effects of parameter choice on defining molecular operational taxonomic units and resulting ecological analyses of metabarcoding data[END_REF][START_REF] Flynn | Toward accurate molecular identification of species in complex environmental samples: Testing the performance of sequence filtering and clustering methods[END_REF].

In fact, even if clustering reads into operational taxonomic units can partially address the overestimation of taxa, clustering-based methods do not account for some LFNs, such as cross-sample contamination and mistagging.

In this context, and following the recommendations by [START_REF] Murray | From benchtop to desktop: Important considerations when designing amplicon sequencing workflows[END_REF], we developed a "from-benchtop-to-desktop" metabarcoding workflow to investigate the diet of invertebrateeaters. We particularly focused on two main objectives. Our first objective was data reliability, which involves both minimizing false-negatives and false-positives. This was achieved by ensuring the efficiency of the combined use of two primer sets, by performing several PCR replicates and by developing a clustering-free filtering method based on control samples and nonarbitrary thresholds. Second, we developed an approach that should maximize the taxonomic resolution of molecular identification of prey by combining four assignment procedures and three reference databases. Our workflow was applied to the biodiversity assessment of prey ingested by Zingel asper (Linnaeus, 1758), a critically endangered benthic freshwater fish. Using this model, we demonstrated the application of our workflow on faecal samples and illustrated its interest for a better characterization of feeding habitats.

| MATERIAL AND METHODS

| Faecal sample collection

Thirty-five Z. asper specimens from which faeces could be collected were sampled on 5 September 2014 in the Durance River (France: 44°20 0 14″N, 5°54 0 46″E). Fishes were caught by electrofishing and their abdomen was squeezed by hand in order to drain their faeces.

Collected faeces were stored in a 1.5-ml vial containing 96% ethanol.

After sampling, individuals were immediately released within the fishing area. Faeces were stored at À20°C until DNA extraction.

Additionally, five faecal samples from a fish species living in brackish-water habitats, Pomatoschistus microps (Krøyer, 1838), were also analysed to assess the versatility of the workflow and to control for mistagging in our HTS data set (see below). These individuals were sampled in the Vaccar es Lagoon (French Mediterranean coast: 44°20 0 14″N, 5°54 0 46″E).

| Faecal DNA extraction and controls

All faecal DNA extraction steps were conducted in a room dedicated to the handling of degraded DNA ("Plateforme ADN D egrad e" of the Institut des Sciences de l'Evolution de Montpellier, France) and following the specific safety measures described by [START_REF] Monti | Being cosmopolitan: Evolutionary history and phylogeography of a specialized raptor, the Osprey Pandion haliaetus[END_REF].

Before DNA extraction, faeces were dried using the Eppendorf Concentrator Plus (Eppendorf, Germany). One volume of dried faeces, one volume of zirconium oxide beads (0.5 mm) and ½ volume of sterile water were mixed to crush samples using a Bullet Blender (Next Advance, USA). The DNeasy â mericon Food Kit (QIAGEN, Germany) was used to extract DNA from faecal samples to minimize the level of co-extracted products and improve PCR success (Zarzoso-Lacoste, [START_REF] Zarzoso-Lacoste | Improving PCR detection of prey in molecular diet studies: Importance of group-specific primer set selection and extraction protocol performances[END_REF]. Each extraction series included (i) 23 faecal samples, (ii) a negative control for extraction (T ext ) that consisted of 50 ll of DNA-free water subjected to DNA extraction protocol and (iii) a negative control for DNA aerosols (T pai ) that consisted of a 1.5-ml vial containing 50 ll of DNA-free water that remained open but otherwise untouched during the extraction protocol. DNA concentrations were quantified using a Qubit Fluorometer (Invitrogen, Darmstadt, Germany) and standardized to max. 20 ng/ll.

| Local DNA library construction

A local (noncomprehensive) DNA library (lcDNA samples) was constructed using invertebrates sampled in the Durance River and invertebrate samples from our laboratory collections. The sample composition and laboratory protocols are detailed in Table S1 and Appendix S1, respectively. We successfully sequenced 301 samples, representing 209 distinct species.

| Minibarcoding protocol and taxonomic coverage

Two DNA primer pairs were initially selected. They both amplify a short fragment from the 5 0 end of the mitochondrial cytochrome c oxidase I (COI): ZBJ-ArtF1c and ZBJ-ArtR2c (Arthropods "universal"; [START_REF] Zeale | Taxon-specific PCR for DNA barcoding arthropod prey in bat faeces[END_REF], hereafter abbreviated as ZF and ZR; and Uni-MinibarF1 and Uni-MinibarR1 (Eukaryotes "universal"; [START_REF] Meusnier | A universal DNA mini-barcode for biodiversity analysis[END_REF], hereafter abbreviated as MF and MR.

All four pairwise combinations of the primers were tested for PCR success on lcDNA samples, and only MFZR and ZFZR were selected as a result (see Appendix S2). These two primer pairs produce overlapping amplicons (from ~210 to ~230 bp including primers) with MFZR amplicons being 18 bp longer in the 5 0 region. The COI region of the predator (Z. asper) was not amplified using either ZFZR or MFZR; however, the MFZR primer pair amplified an unspecific fragment (~850 bp). Therefore, a blocking primer (BlNupRan) was developed to inhibit the amplification of the nontargeted locus (see Appendix S2). The taxonomic coverage of both primer sets was assessed by in vitro tests using the lcDNA samples and in silico using EcoPCR, EcoTaxStat and EcoTaxSpecificity [START_REF] Ficetola | An In silico approach for the evaluation of DNA barcodes[END_REF] (see Appendix S2).

| PCR-based enrichment of prey DNA and high- throughput sequencing

To track amplicons back to the samples and to avoid flashes of light during HTS, 12-14 nucleotide-long sequence tags were added onto the 5 0 end of each primer (eight distinct tags for the ZF and MF primers, and 12 tags for ZR), creating 96 forward and reverse potential tag combinations (Table S2). Three PCR replicates were conducted in a volume of 25 ll for each minibarcode primer pair, resulting in six independent PCR enrichments per sample. Template DNA consisted of 2.5 ll of faecal standardized DNA extracts. For the primer pair MFZR, blocking primer BlNupRan was added in the PCR mix at 400 nM. Additionally, several control samples were subjected to the PCR enrichment step: (i) T ext and T pai (see above) are controls for cross-sample or exogenous contaminations during the DNA extractions; (ii) T PCR indicates the level of cross-contamination during the preparation of the PCR mix and plates (tagged primers but no DNA template in the PCR vial); (iii) one negative control (T tag ) was included to assess the level of mistagging due to the recombination of sequences from different samples (see [START_REF] Schnell | Tag jumps illuminated-reducing sequence-to-sample misidentifications in metabarcoding studies[END_REF]. The T tag consisted of an empty vial on the PCR plate that did not contain any tagged primers or DNA. This creates an extra tag combination that is not used for any samples or controls (as in: [START_REF] Esling | Accurate multiplexing and filtering for high-throughput amplicon-sequencing[END_REF]; (iv) finally, two mock faecal samples were also amplified (T pos1 and T pos2 ) and they consisted of identical artificial mixes of the DNA obtained from six potential prey (Caenis pusilla, Baetis rhodani, Orthocladiinae sp., Chironomidae sp., Hydropsyche pellucidula, Phoxinus cf. phoxinus) and from Z. asper. The DNA concentration of T pos1 and T pos2 was 0.2 ng/ll for each potential prey, and 0.8 ng/ll for Z. asper. T pos1 and T pos2 were used to gauge sequencing or PCR artefacts and evaluate reliability of our metabarcoding analyses (see De [START_REF] De Barba | DNA metabarcoding multiplexing and validation of data accuracy for diet assessment: Application to omnivorous diet[END_REF]. The DNA from the specimens used for our mock samples was extracted in a room free of DNA handling. Moreover, the COI reference sequences of each specimen were originally obtained by amplification using CK4 primer set (see Appendix S1) and Sanger sequencing (except the Chironomidae sp. sample that failed to be amplified with CK4). Amplicons were checked by gel electrophoresis and were then pooled by replicate series, that is, 48 samples for each of the six replicate series (two primer pairs, three replicates each). Electrophoretic migrations were carried out using 20 ll of each pool on an agarose gel at 1.25%. The amplicons with expected sizes were isolated using a sterile scalpel, and the DNA was purified using the PureLink â Quick Gel Extraction Kit (Life Technologies, Germany). About 20 ng of purified DNA from each replicate series was used to generate six Illumina sequencing libraries using the TruSeq â Nano DNA Sample Preparation Kit (Illumina, San Diego, CA, USA). The library preparation included DNA end-repairing, A-tailing, Illumina adapter ligation, and limited amplification (12 PCR cycles). Illumina libraries were distinguished by the ligation of distinct adapters. The libraries were controlled for size and quality using the Agilent Bioanalyzer DNA 1000 Kit (Agilent Technologies, Palo Alto, CA, USA) and for DNA concentrations with the Kapa Library Quantification Kit for Illumina â platforms (KapaBiosystems, Wilmington, MA, USA). The six libraries were pooled at equimolar concentration (4 nM) and sequenced on an ILLUMINA MISEQ v3 platform as paired-end 250nucleotide reads.

| Sequence processing and filtering

We developed a pipeline (see Figure 1) to filter the obtained MiSeq data (see Appendix S3 for a detailed version of the filtering pipeline).

Unless specified otherwise, bioinformatics processing of reads was performed using custom Perl scripts (Dryad; https://doi.org/10. 5061/dryad.f40v5) and statistical analyses were conducted using R software (R Development Core Team 2014). Reads from the different primer pairs and replicate series were sorted according to the Illumina adapter sequences. PEAR v0.9.5 [START_REF] Zhang | PEAR: A fast and accurate Illumina Paired-End reAd mergeR[END_REF] was used to merge read pairs and discard lowquality reads (Figure 1, Step 1). Merged reads were then assigned to samples and replicates using a BLAST-based approach, and tags and primers were trimmed from the reads. Only reads that had a perfect match to tags were accepted. This constituted an additional step to eliminate low-quality reads. Strictly identical trimmed reads within each of the six replicate series were pooled into variants (i.e. dereplicated reads) and singletons (i.e. only one read in a replicate series) were discarded (Figure 1, Step 2). Variants of the three replicate series obtained from the same primer set (MFZR or ZFZR) were pooled.

The read number associated with each variant-replicate combination, however, was kept.

The following filtering steps were done separately for MFZR and ZFZR. Variants that did not comply with our BLAST conditions (Evalue threshold: 1e-10; minimum query coverage: 80%) against our custom COI database (COI-filtering-DB; Dryad, https://doi.org/10. 5061/dryad.f40v5) were considered as non-COI sequences and discarded (Figure 1, Step 3).

Then, LFN filters were used to discard variants with low read counts, likely originating from contamination, mistagging or sequencing/PCR errors bias. All variant-replicate combinations were considered independently. For each replicate of a sample, all the variants associated with a read count or a relative frequency that were under one of the LFN thresholds (i.e., not distinguishable from the noise inherent to the HTS data) were removed. Three different thresholds were considered. In fact, a variant can be rare in a given replicate (i) compared to the total number of reads in the replicate (N repl ), (ii) compared to its total number of reads in the run (N var ) or (iii) if it has few reads in a replicate in absolute terms (N var-repl ). The first threshold (LFN pos = 0.3%) was based on the relative frequency (N var-repl /N repl ) of the least frequent expected variant in all mock community replicates. This threshold helps to eliminate most low-frequency variants due to the bias cited above (Figure 1, Step 4a).

The second threshold (LFN tag = 0.25%) was determined based on (i) freshwater prey variants present in samples of brackish-water fish samples (P. microps), (ii) brackish-water prey variants present in samples of freshwater fish samples (Z. asper) and (iii) unexpected variants present in the mock samples. The LFN tag allowed discarding variants that appear by mistagging or cross-sample contamination due to their high frequency in the run (Figure 1 When occurring early in PCR cycles, errors can generate variants with quite high frequency compared to the error-free sequences.

We therefore used the Obiclean program (Obitools package; [START_REF] Boyer | Obitools: A unix-inspired software package for DNA metabarcoding[END_REF] to further filter out variants resulting from PCR and/or sequencing errors. Based on the known variants of mock samples, a 20% threshold was determined, and all variants classified as "internal" (see: [START_REF] Shehzad | Prey preference of snow leopard (Panthera uncia) in South Gobi, Mongolia[END_REF][START_REF] De Barba | DNA metabarcoding multiplexing and validation of data accuracy for diet assessment: Application to omnivorous diet[END_REF]Appendix S3) were discarded (Figure 1 Step 4b

LFN tag

Filter out mistagging from highly amplified variants, where LFN pos does not work

Step4

All LFN

Step8

Variants present in at least 2 replicates per sample

Ensure variant consistency between replicates

Step9

Concensus over replicates

Step 4c

LFN neg

Eliminate low frequency contamination from replicates with low read count (e.g., in T neg )

Step10

Chimera filter (UCHIME2)

Step11

Pseudogene filter

Step1

Merge read pairs with PEAR

Includes quality filtering

Step2

Assign reads to replicates Dereplicate Delete singletons Accept only reads with perfect match to tags and delete singletons

Step3

Filter out non-COI reads In parallel, the data set that passed filtering steps 1-3 was also analysed by UNOISE2 [START_REF] Edgar | UNOISE2: Improved error-correction for Illumina 16S and ITS amplicon sequencing[END_REF] for denoising reads. UNOISE2

F I G U R E
intends to eliminate all noise coming from sequencing and PCR errors. However, UNOISE2 does not deal with mistagging and contamination. We therefore further filtered the data by UNCROSS [START_REF] Edgar | UNCROSS: Filtering of high-frequency cross-talk in 16S amplicon reads[END_REF]. The UNOISE2/UNCROSS filtering procedure is roughly comparable to our filtering steps 4 and 6. To make the results of this second filtering approach comparable to our pipeline, we further applied our filtering steps 5, 7 and 8 to ensure repeatability within samples, as well as filtering steps 10 and 11 to eliminate chimeras and pseudogenes variants, respectively. Second, we used an automatic procedure that assigned each variant/contig to the lowest taxonomic group (LTG) of BLAST hits against a custom-built local database (Taxassign-DB; Appendix S4). The principle of our LTG approach is similar to the lowest common ancestor developed by [START_REF] Huson | Integrative analysis of environmental sequences using MEGAN4[END_REF]. All variants/contigs were BLASTed against the Taxassign-DB (BLASTN, Evalue 1e-10, minimum coverage of the query sequence: 90%, subject is annotated to a family or lower level). The LTG was determined based on different similarity cut-offs (S) from 100% to 70%. At each S, the LTG was defined as the lowest taxonomic group that contained at least 90% of the selected hits and contains at least three taxa for S < 97%.

| Taxonomic assignation of variants

The LTG corresponding to the highest S was selected.

Third, we used the Identification Request of BOLD Systems Tools (www.boldsystems.org; Ratnasingham & Hebert, 2007) using the COI "All Barcode Records" database (performed in September 2015). This approach involved all the sequences of the BOLD Systems, even those not published yet.

Fourth, when a cross-comparison between the assignment analyses revealed inconsistencies or when the taxonomic assignments had low resolution, we performed complementary phylogenetic analyses of the variants and contigs. For each of these ambiguous variants/contigs, we selected ten sequences from the Taxassign-DB by retaining: (i) the best BLAST hits (BLASTN, E-value 1e-10, minimum coverage of the query sequence: 90%); (ii) a maximum of three sequences per taxon;

and (iii) only sequences annotated to at least the Family level. Neighbor-joining trees were constructed based on K2P distances with 1,000 bootstraps using MEGA 6.06 [START_REF] Tamura | MEGA6: Molecular Evolutionary Genetics Analysis version 6.0[END_REF]. The topology of the phylogenetic trees, especially the apices of branches, was used to resolve the taxonomy of variants and contigs. When the tree topologies did not resolve the taxonomic ambi-

| Additional analyses

The variants and contigs validated in the Z. asper and P. microps samples were used to further evaluate the coverage of MFZR and ZFZR primer pairs. To standardize this evaluation, a complete-linkage clustering [START_REF] Sneath | Numerical taxonomy: The principles and practice of numerical classification[END_REF] was used to delineate molecular operational taxonomic units (MOTUs; [START_REF] Blaxter | Defining operational taxonomic units using DNA barcode data[END_REF] among validated variants and contigs. A maximum of 3% divergence between variants/contigs was allowed within a MOTU. This maximal divergence was previously used as a proxy for delineating invertebrates' species (e.g., [START_REF] Clare | The diet of Myotis lucifugus across Canada: Assessing foraging quality and diet variability[END_REF][START_REF] Vesterinen | What you need is what you eat? Prey selection by the bat Myotis daubentonii[END_REF]. For taxa that present a higher within-species sequence divergence, however, this 3% threshold will allow taking into account possible differential intraspecific coverage.

In diet studies, the quantification of individuals per prey would add considerable ecological information. However, due to PCR biases, read counts are generally considered as inappropriate to assess the relative abundance of prey [START_REF] Elbrecht | Can DNA-based ecosystem assessments quantify species abundance? Testing primer bias and biomass-Sequence relationships with an innovative metabarcoding protocol[END_REF].

Alternatively, we used a variant-centred approach: we assumed that Two categories of prey were taken into consideration during our analysis: invertebrates (excluding Rotifera) and microorganisms (including diatoms, red and green algae, and Rotifera). In this work, we focused on invertebrates as we collected faeces from invertebrateeating species. In contrast, microorganisms are less relevant since some uncertainty remains about their ingestion (secondary predation and/or nontargeted ingestion when predating invertebrates).

| RESULTS

| HTS data filtering

High for the preparation of the mock samples (Figure 4a, Table S3), and one unexpected variant (MFZR_082136) assigned to an undetermined Eukaryota (Fig. S4). This variant is present in both mock community samples but absent from all other samples. Therefore, it is likely that it comes from an organism ingested by or attached to one of the prey individuals used for the construction of the mock sample.

When using the pipeline based on UNOISE2 and UNCROSS for filtering our data set, a total of 18 variants were retained for each mock community sample (17 present in both). While this pipeline allowed the validation of the six variants expected in the mock samples, it retained 12 additional and unexpected variants (including MFZR_082136).

| Taxonomic identification and resolution

The 99 variants and contigs identified for Z. asper and P. microps faecal samples were analysed using the SAP, BOLD, LTG and FTA analyses (see Table S3). Regarding the variants and contigs assigned to invertebrates in Z. asper samples (n = 42), the taxonomic resolution of each analysis was quantified using the IR index (Figure 2a).

The IR of FTA (IR FTA = 5.8 AE 0.3) was significantly higher than those calculated for SAP (IR SAP = 4.1 AE 1.1), BOLD (IR BOLD = 5.2 AE 1.2) and LTG (IR LTG = 5.7 AE 0.3). The FTA approach resulted in the assignment of most variants and contigs to the species level (Figure 2b). It should be noted that a very close performance (although significantly lower) was obtained by the fully automatic LTG approach.

During the cross-validation step, the taxonomic assignments of SAP, LTG and BOLD were compared: (i) in 6% of the cases, all three analyses suggested the same taxon; (ii) in 87% of cases, the analyses were consistent but assigned the variants/contigs to different taxonomic levels and the lowest taxonomic level was retained;

(iii) in 7% of cases, the three approaches produced conflicting assignments. Variants or contigs with unsatisfactorily taxonomic levels or conflicting assignments were subjected to a phylogenetic analysis. All variants and contigs assigned to Simulium sp. were also included in the phylogenetic analyses, since at least one of them showed a conflicting assignment. Phylogenetic analyses were conducted for a total of 33 variants and contigs (sequence alignment deposited in Dryad; https://doi.org/10.5061/dryad.f40v5), which resulted in (i) the resolution of the assignment of nine variants and

(ii) the refinement of the identification of a further 14 variants (Figs S3 andS4). In one case, biogeographical data were also considered for the FTA decision (Table S3; contig_0260 identified as Antocha vitripennis since it is recorded in sampling area, whereas Austrophorocera Janzen04 sequences originate from Costa-Rica).

T

| Taxonomic coverage of the minibarcode primer set

In vitro tests showed that 87% and 85% of the taxa were amplified by ZFZR and MFZR, respectively. The combination of both primer pairs leads to a mean amplification success of 94% (Table 2): 100%

for Arthropoda, 80% for Mollusca, 90% for Annelida and 65% for Teleostei. The PCR success was generally similar among samples of the same species, with a few exceptions (Table S1).

In silico tests showed that the taxonomic coverage at the species level for Arthropoda was 63% and 66% for MFZR and ZFZR, respectively (Table S4). When combining the coverage of the two primer pairs, the taxonomic coverage increased to 77% for Arthropoda, varying from 62% (Hymenoptera) to 100% (Lepidoptera).

The discriminating power of the minibarcoding amplicons was also estimated in silico: 79%-100% of the analysed species could be differentiated by at least one mutation for the analysed taxa (Table S4). To further evaluate the complementarity of the two primer pairs (ZFZR and MFZR), we looked at their performance in detecting the 81 MOTUs identified in the Z. asper and P. microps faecal samples (Table S5): ~30% MOTUs were detected with both primer pairs, 23% were detected with ZFZR only, and 47% were detected with MFZR only (Figure 3a). When focusing on invertebrates, 53% of MOTUs were detected with one of the two primer sets only (Figure 3b). The ZFZR pair was more efficient in detecting Diptera (both for Chironomidae and Simuliidae), whereas MFZR detected more Baetidae (Ephemeroptera) and Crustacea.

SAP

| Sample composition

The number of variants and contigs in the Z. asper faecal samples varied from 1 to 24 (6.6 AE 5.5). At least one variant of invertebrates was present in all Z. asper samples, with a mean value of 3.4 (AE 1.7) (Table S3; Figure 4a). These samples also varied in terms of prey composition: whereas Baetis was detected in 79% of the samples, the frequency of other prey varied widely. Ephemeroptera, Diptera and Trichoptera appeared to be the most abundant prey at the population scale (Figure 4b). Four Ephemeroptera families were identified with a predominance of Baetidae (five species), which also displayed the highest MNI (22% of the total MNI in the Z. asper samples; Figure 4b). The main groups of Diptera that were preyed upon were Simuliidae and Orthocladiinae (subfamily of Chironomidae). Noninvertebrate taxa were mainly composed of metazoan microorganisms: Bacillariophyta, Rotifera, Rhodophyta and Oomycetes. Furthermore, none of the dietary MOTU accumulation curves approached an asymptote (Figure 4c).

Regarding the P. microps samples, prey DNA was detected in samples P3 and P5 only (one and four variants/contigs, respectively), and the DNA of P. microps was detected in all five faecal samples (Table S3).

| DISCUSSION

| Producing robust HTS data sets

The choice of the primer sets is crucial in diet metabarcoding as it may lead to false negatives due to insufficient taxonomic coverage [START_REF] Pompanon | Who is eating what: Diet assessment using next generation sequencing[END_REF]. In this study, we used two primer sets to amplify a short but informative COI region (MFZR and ZFZR). In vitro tests on our lcDNA samples showed that the primer sets were complementary and enabled coverage of a large taxonomic spectrum, especially for invertebrates (100% of items), whereas in silico tests suggested lower coverage. Our Z. asper faecal samples displayed a taxonomic diversity comparable to a previous morphological-based diet study [START_REF] Cavalli | Diet and growth of the endangered Zingel asper in the Durance River[END_REF]; summary in Fig. S5). Interestingly, the complementarity of primer sets was higher in the Z. asper faecal samples (53% of invertebrates MOTUs detected by one of the two primer sets only) compared to in vitro tests (11%). Increased complementarity in multiplexed samples is likely to be the result of preferential primer binding (e.g. [START_REF] Thomas | Quantitative DNA metabarcoding: Improved estimates of species proportional biomass using correction factors derived from control material[END_REF], which justifies the use of more than one primer pair. In our case, the selection of two primer pairs led to the detection of all prey in the positive samples and to a high prey diversity within and among Z. asper samples. Furthermore, the successful detection of prey in P. microps faecal samples suggested that these primer sets are applicable for species that live in nonfreshwater environments. However, in vitro tests on lcDNA samples revealed that the combined use of ZFZR and MFZR displayed insufficient taxonomic coverage in some groups (e.g. 80% of Mollusca). Nevertheless, our workflow could easily be adapted to other predators and environments, thanks to the availability of several "universal" minibarcodes located in the COI region that we have targeted (e.g., [START_REF] Brandon-Mong | DNA metabarcoding of insects and allies: An evaluation of primers and pipelines[END_REF][START_REF] Hajibabaei | Environmental barcoding: A next-generation sequencing approach for biomonitoring applications using river benthos[END_REF][START_REF] Leray | A new versatile primer set targeting a short fragment of the mitochondrial COI region for metabarcoding metazoan diversity: Application for characterizing coral reef fish gut contents[END_REF][START_REF] Shokralla | Massively parallel multiplex DNA sequencing for specimen identification using an Illumina MiSeq platform[END_REF].

The most commonly used filtering pipelines when dealing with HTS metabarcoding data, such as mothur [START_REF] Schloss | Introducing mothur: Open-source, platform-independent, community-supported software for describing and comparing microbial communities[END_REF] and QIIME [START_REF] Caporaso | QIIME allows analysis of highthroughput community sequencing data[END_REF], rely on the clustering of reads into MOTUs. In these clustering-based approaches however, almost each step and each parameter in the bioinformatics pipeline influence the outcome, and the number of taxa is often overestimated (e.g. [START_REF] Brown | Divergence thresholds and divergent biodiversity estimates: Can metabarcoding reliably describe zooplankton communities?[END_REF][START_REF] Clare | The effects of parameter choice on defining molecular operational taxonomic units and resulting ecological analyses of metabarcoding data[END_REF][START_REF] Majaneva | Bioinformatic amplicon read processing strategies strongly affect eukaryotic diversity and the taxonomic composition of communities[END_REF]. Therefore, each data set will need a specific analysis method and parameters should be tailored to the purpose of the study [START_REF] Flynn | Toward accurate molecular identification of species in complex environmental samples: Testing the performance of sequence filtering and clustering methods[END_REF]. We therefore favoured and developed a clustering-free pipeline that relies on a series of stringent filtering , 2014). For this purpose, the thresholds used in our filtering approach were inferred from mock community samples (Step 4a) and negative controls (Step 4c). It has to be noted that this conservative approach will discard all variants for which frequencies are below the LFN thresholds. In the context of complex community samples, such as samples collected for biodiversity assessment and monitoring (e.g., Elbrecht & Leese, 2017;[START_REF] Lanz En | High-throughput metabarcoding of eukaryotic diversity for environmental monitoring of offshore oil-drilling activities[END_REF][START_REF] Leray | DNA barcoding and metabarcoding of standardized samples reveal patterns of marine benthic diversity[END_REF], some low-abundance taxa may not be validated. In fact, their read count may fall below one or more LFN thresholds, making the corresponding variant undistinguishable from noise. In this case, the analysis of several biological replicates (i.e., distinct DNA extractions corresponding to different fractions of a given sample) should help with the detection and the validation of low-abundance taxa (Lanz en, Lekang, Jonassen, [START_REF] Lanz En | DNA extraction replicates improve diversity and compositional dissimilarity in metabarcoding of eukaryotes in marine sediments[END_REF][START_REF] Zhan | Reproducibility of pyrosequencing data for biodiversity assessment in complex communities[END_REF]. In addition, we recommend that the mock samples approximate the complexity of the communities sampled, by approximating their expected taxonomic composition and even the differential abundance of taxa.

Ignoring the mistagging bias [START_REF] Schnell | Tag jumps illuminated-reducing sequence-to-sample misidentifications in metabarcoding studies[END_REF] may lead (i) to overestimate the number of taxa in a sample and (ii) to blur the differentiation between samples with respect to their taxonomic composition. To assess and overcome mistagging bias, two strategies were recently suggested: (i) increasing the number of unused tag combinations using Latin square design [START_REF] Esling | Accurate multiplexing and filtering for high-throughput amplicon-sequencing[END_REF]; and (ii) using fusion primers [START_REF] Herbold | A flexible and economical barcoding approach for highly multiplexed amplicon sequencing of diverse target genes[END_REF] to avoid the creation of intersample chimeras when pooling samples during the sequencing library preparation. In the last case however, PCR efficiency is reduced substantially [START_REF] Schnell | Tag jumps illuminated-reducing sequence-to-sample misidentifications in metabarcoding studies[END_REF] dependent approach (LFN tag ; Step 4b) to control mistagging. In our case, the threshold of the LFN tag filter is based on mock community samples and on the co-analysis of samples from different habitats (brackish-water and freshwater) in the same HTS run.

In our workflow, LFN thresholds appear as the most critical parameters for the filtering and the validation of HTS data: most variants that are expected to be false positives were discarded at Step 4 (see Table 1). Nevertheless, subsequent filtering steps (steps 5-11) are also important for the validation of the data. Expectedly, omitting the filtering steps 5-11 will decrease the robustness of the data and inflate the taxonomic and/or genetic diversity within samples, biasing biodiversity estimations. In our case, for example, omitting filtering steps 5-11 is expected to inflate the MNI estimator (see below). As a matter of fact, further false positives were controlled, such as variants with PCR/sequencing errors (Step 6) and chimeras (Step 10). Moreover, our workflow included three PCR replicates for each primer pair to avoid falsepositive detections (see [START_REF] Ficetola | How to limit false positives in environmental DNA and metabarcoding[END_REF]. These PCR replicates were sequenced separately and used at filtering steps 5, 7

and 8 for validating reproducible variants only. Several authors, however, used PCR replicates differently in the case of complex community samples (e.g., [START_REF] Lanz En | High-throughput metabarcoding of eukaryotic diversity for environmental monitoring of offshore oil-drilling activities[END_REF][START_REF] Leray | Random sampling causes the low reproducibility of rare eukaryotic OTUs in Illumina COI metabarcoding[END_REF]: PCR replicates were pooled and then sequenced jointly to limit the impact of random sampling amplification biases and hence minimizing false negatives due to low-abundance taxa. Nevertheless, this approach does not allow minimizing the false-positive detections. In our case, the size of faecal samples allowed one single DNA extraction reaction. However, in the case of larger or more complex samples, several biological replicates (i.e., distinct fractions of the same sample) may be analysed for controlling false negatives (see [START_REF] Lanz En | DNA extraction replicates improve diversity and compositional dissimilarity in metabarcoding of eukaryotes in marine sediments[END_REF][START_REF] Zhan | Reproducibility of pyrosequencing data for biodiversity assessment in complex communities[END_REF].

The combined use of biological replicates and of PCR replicates (sequenced separately) should ensure both avoiding false-positive detections and maximizing the coverage of taxa diversity.

Alternatively, for complex community samples, more relaxed filtering thresholds (at steps 4a and 4c) and more relaxed reproducibility criteria (at steps 5, 7 and 8) might be considered if one wants to maximize the detection low-abundance taxa. In this case however, the presence of these taxa in the final data set will remain uncertain since confounded with the low-frequency noise and the false positives. These low-frequency taxa will therefore require to be validated by some complementary field data (for estimating the probability of false negatives, as suggested by: [START_REF] Leray | Random sampling causes the low reproducibility of rare eukaryotic OTUs in Illumina COI metabarcoding[END_REF] [START_REF] Gonz Alez-Tortuero | The quantification of representative sequences pipeline for amplicon sequencing: Case study on within-population ITS1 sequence variation in a microparasite infecting Daphnia[END_REF][START_REF] Shokralla | Massively parallel multiplex DNA sequencing for specimen identification using an Illumina MiSeq platform[END_REF] or prey biomass [START_REF] Jo | Discovering hidden biodiversity: The use of complementary monitoring of fish diet based on DNA barcoding in freshwater ecosystems[END_REF]. In this study, thanks to rigorous filtering procedures we produced a robust data set of variants and contigs that is reliable and relatively free from false positives and artefacts. Moreover, the variants corresponding to pseudogenes were filtered out.

Therefore, we are confident that the number of COI variants/contigs can be used to estimate MNIs. However, in some cases, that are prey taxa prone to heteroplasmy or that present tissue-specific mitochondrial variants, an overestimation bias should be considered when using

MNIs. The MNI statistic may therefore deserve a further evaluation by using appropriate mock samples as controls.

In the case of Z. asper prey, the biological analyses suggested that the final variants can be used as a reliable estimation of their DNA diversity in faeces and therefore that the MNI statistic adequately reflects relative prey abundance. In fact, congruence between the MNI and previously observed prey diversity based on morphological gut-content analysis in Durance River [START_REF] Cavalli | Diet and growth of the endangered Zingel asper in the Durance River[END_REF] was observed. More specifically, the family Baetidae exhibited the highest MNIs in Z. asper faecal samples, which was also the most abundant prey found in morphological gut-content analysis (Fig. S5).

| Taxonomic assignment of prey

Taxonomic assignment procedures are critical in metabarcoding studies, and several innovative approaches were recently developed (e.g., [START_REF] Porter | Rapid and accurate taxonomic classification of insect (class Insecta) cytochrome c oxidase subunit 1 (COI) DNA barcode sequences using a na€ ıve Bayesian classifier[END_REF][START_REF] Somervuo | Unbiased probabilistic taxonomic classification for DNA barcoding[END_REF][START_REF] Somervuo | Quantifying uncertainty of taxonomic placement in DNA barcoding and metabarcoding[END_REF]. However, most metabarcoding studies are based on one single taxonomic assignment approach, and yet using different procedures as well as different reference databases should improve the reliability and the accuracy of MOTU identification. For this study, we combined different assignment procedures to combine their respective advantages. Although our strategy can be time-consuming-some of the steps cannot be fully automated-the IR statistic showed that it did significantly improve the precision of the assignments. This resulted in a very low proportion of high-level taxonomic assignments (e.g., Eukaryota), and a high proportion of low-level taxonomic assignments for variants of invertebrates (species level for most). In contrast, most invertebrates identified by morphological-based gut-content analysis were identified at the Family level in [START_REF] Cavalli | Diet and growth of the endangered Zingel asper in the Durance River[END_REF]. Moreover, our detailed assignment of Z. asper prey supports the use of COI as the favoured target gene for invertebrates (see also [START_REF] Brandon-Mong | DNA metabarcoding of insects and allies: An evaluation of primers and pipelines[END_REF].

The COI appears as the one of the most appropriate gene for invertebrate metabarcoding because (i) it reveals species-level variation [START_REF] Elbrecht | Testing the potential of a ribosomal 16S marker for DNA metabarcoding of insects[END_REF], and (ii) a huge amount of annotated sequences is available in public databases [START_REF] Deagle | DNA metabarcoding and the cytochrome c oxidase subunit I marker: Not a perfect match[END_REF].

The taxonomic assignment of invertebrates indicated that the Baetidae found in Z. asper faeces belonged mostly to the genus Baetis although other genera were also present in the sampling area (e.g., Alainites, Acentrella, Centroptilum and Procloeon). According to [START_REF] Tachet | Invert ebr es d'eau douce: Syst ematique[END_REF], Baetis has a higher affinity towards coarse substrates (e.g., stones, pebbles) with higher water velocity than other Baetidae, which in turn prefer epiphyte lifestyle on macrophyte or algae substrate. Additionally, the genera Hydropsyche (Hydropsychiidae) and Simulium (Simuliidae), all epibenthic and rheophilic taxa, occurred at non-negligible frequencies in Z. asper faecal samples.

Consequently, the cumulated frequency based on MNI of Baetis, Hydropsyche and Simuliidae was ~70% (see Figure 4). This suggests that Z. asper actively selected rheophilic and epibenthic macroinvertebrates, which constitutes accurate information related to habitat use of this critically endangered species.

| CON CLUSION AND PERSPECTIVE S

Diet analyses are critical to gain a better understanding of prey/habitat relationships and feeding habitats (e.g., [START_REF] Corse | A PCR-based method for diet analysis in freshwater organisms using 18S rDNA barcoding on faeces[END_REF][START_REF] Sanchez-Hernandez | Age-related differences in prey-handling efficiency and feeding habitat utilization of Squalius carolitertii (Cyprinidae) according to prey trait analysis[END_REF], especially when a fine and accurate taxonomic identification of prey can be achieved (e.g., Adamczuk & Mieczan, 2015). In this study, stringent wet-laboratory conditions, carefully selected primer sets, PCR replicates and nonarbitrary filtering thresholds based on control samples led to the validation of a robust data set dedicated to the study of the diet of a critically endangered fish species. The robustness of our data set allowed us to take into account the prey genetic variability and to obtain a semiquantitative estimate of diet through the use of MNI. Furthermore, a species-level identification for most of the invertebrate prey was obtained through a complementary taxonomic assignment approach. On the whole, our approach produced a robust data set for ecological analyses and opens perspective for more precision on feeding habitats of invertebrate-eaters. This new information would in turn improve conservation strategies such as habitat restoration or the selection of optimal re-introduction sites.

Finally, our results reinforced previous findings that suggested diet metabarcoding can be a powerful tool in trophic ecology as it allows the determination of large scale and highly resolved networks [START_REF] Evans | Merging DNA metabarcoding and ecological network analysis to understand and build resilient terrestrial ecosystems[END_REF]. It may also produce a level of precision that reveals unexpected food web structures [START_REF] Roslin | The use of DNA barcodes in food web construction-terrestrial and aquatic ecologists unite![END_REF]. This article proposes a from-benchtop-to-desktop workflow that provides an efficient tool for the study of invertebrateeater diets and will, we hope, stimulate and inspire future trophic works using metabarcoding approaches.

ACKNOWLEDG EMENTS

We thank three anonymous reviewers for their constructive comments, which improved our manuscript significantly. We warmly thank

  , Step 4b). Finally, for replicates that had a low number of reads, contamination could not necessarily be discarded with LFN pos and LFN tag filters. Therefore, a last threshold (LFN neg ) based on the maximum read count (N var-repl ) among all variants of all negative controls replicates was used (LFN neg : 31 and 53 for MFZR and ZFZR, respectively; Figure 1, Step 4c). The three LFN filters were run in parallel on the merged and dereplicated COI variants and only variant-replicate combinations that passed all three LFN filters were retained (Figure 1, Step 4), except those that were present in only one replicate within a sample (Figure 1, Step 5).

  , Step 6). By comparing replicates within samples, we were able to assess the repeatability of the experimental procedure in two different ways: (i) only variants present in at least two replicates of a sample were retained (Figure 1, Step 5), and (ii) distance between replicates was taken into account (Figure 1, Step 7). To this latter end, we followed the strategy developed by De Barba et al. (2014) and used the Renkonen distance (RD) to compare distances between replicates of the same sample. The threshold was set at the 10% upper tail of the distribution of the RDs, which corresponds in our data set to RD = 0.157 for MFZR, and RD = 0.088 for ZFZR (Fig. S1). All those replicates separated from the other replicates by RD above the defined threshold were discarded. Furthermore, samples with only one replicate were excluded (Figure 1, steps 7-8). After the filtering, consensus samples were created by averaging read counts over replicates (Figure 1, Step 9). The pool of all remaining variants and contigs was filtered for chimeras using UCHIME2 (Edgar, 2016a; Figure 1, Step 10) and for potential pseudogenes

1

  High-throughput sequencing filtering pipeline(Figure 1, Step 11). Variants obtained from ZFZR and MFZR primer pairs that perfectly matched on their overlapping regions were merged into contigs. At this stage, the read count values associated with variants were disregarded, and only the presence/absence of the variants/contigs was considered in a given sample.

  Four different approaches were used for the taxonomic assignment of variants and contigs. More detailed descriptions of the methods are given in Appendix S4. First, the phylogenetic approach implemented in the Statistical Assignment Package (SAP;[START_REF] Munch | Statistical assignment of DNA sequences using Bayesian phylogenetics[END_REF] was used to build phylogenetic trees for each of the variants/contigs and their homologues in GenBank using ≥95%, ≥85% and ≥70% sequence identity thresholds in three consecutive runs. The posterior probability (pp) for the query sequence to belong to a clade was estimated at different taxonomic levels.

  guities, biogeographical data were then considered[START_REF] Tachet | Invert ebr es d'eau douce: Syst ematique[END_REF]; OPIE-Benthos database: www. opie-benthos.fr). The combination of the cross-comparison between SAP, LTG and BOLD assignment analyses and the phylogenetic and biogeographical approaches led to a final taxonomic assignment (FTA).To assess the efficiency of the different assignment methods, we focused on the invertebrates' variants and contigs identified in Z. asper faecal samples. An Identification Resolution index (IR; Zarzoso-Lacoste et al., 2016) was calculated for each Z. asper sample and for each assignment analysis. A score was attributed to each variant/contig according to its taxonomic rank assignment, where the maximal value is given to the species level (i.e., Species = 6, Genus = 5, Family = 4, Order = 3; Class = 2, Phylum = 1, Kingdom or NA = 0) and the IR represented the mean score among the variants of a given sample. Identification Resolution indices of the different assignment methods were compared by pairwise nonparametric Wilcoxon rank sum, controlled for multiple comparisons with the Benjamini & Hochberg (1995) procedure.

  the number of distinct variants/contigs represented the Minimal Number of Individuals (MNI; White, 1953) ingested by the predator, providing a semiquantitative estimation of the diet composition. The MNI was used to summarize the proportion of prey items or MOTUs in each faecal sample as well as at the population level. Furthermore, to assess the relation between sample size and MOTU diversity, we constructed two sample-based rarefaction curves (for invertebrates and for all prey) using the observed-richness function implemented in EstimateS v9.1.0 (http://viceroy.eeb.uconn.edu/estimates/).

  Taxonomic assignment analysis. (a) The distribution of index resolution (IR; using invertebrates only) among Zingel asper samples was represented by box plots for each taxonomic assignment procedure. Significant differences between taxonomic analyses were indicated at the top; n.s., nonsignificant; *p < .05; **p < .005; ***p < .0005). (b) The taxonomic levels of the Z. asper invertebrate prey identified using the FTA procedure T A B L E 2 Taxonomic coverage of the metabarcoding primers (in vitro tests)

  U R E 3 Complementarity of minibarcoding primer sets based on Zingel asper and Pomatoschistus microps samples steps based on read counts of each variant in each replicate (for a similar approach, see De Barba et al.

  and template-specific amplification bias may be inflated (O'Donnell, Kelly, Lowell, & Port, 2016). Following Esling et al. (2015), we chose a variant frequency-The diet composition of Zingel asper. (a) Invertebrate composition by sample based on absolute minimal number of individuals (MNIs). (b) Invertebrate composition at the population level, and proportions of microorganisms and invertebrates based on cumulative MNIs. (c) Dietary molecular operational taxonomic unit accumulation curves for Z. asper samples (AE95% confidence intervals)

  Caroline Costedoat, Andr e Gilles, Georges Olivari and the technicians of the French Office National de l'Eau et des Milieux Aquatiques (ONEMA) for their assistance in the field, especially Guillaume Verdier. This work is part of the French Plan National d'Action en faveur de l'apron du Rhône (2012-2016), supervised by the Direction Régionale pour l'Environnement, l'Am enagement et le Logement de Rhône-Alpes and coordinated by the Conservatoire d'Espaces Naturels Rhône-Alpes. E.C. was supported by a postdoctoral grant from Electricit e de France (EDF) and ONEMA. This study was funded by the Syndicat Mixte d'Am enagement du Val Durance (SMAVD), the Agence de l'Eau Rhône-M editerran ee-Corse (AERMC) and the Conseil R egional de Provence-Alpes-Côte d'Azur. This work was conducted in accordance with permits from the French Direction D epartementale des Territoires des Hautes-Alpes (DDT 05). Data used in this work were partly produced through the molecular facilities of LabEx CeMEB (Montpellier) and CIRAD (Montpellier).

  -throughput sequencing of the six amplicon libraries (288 PCR products: 40 faecal samples, 2 T pos , 2 T ext , 2 T pai , 1 T PCR , 1 T tag ; 2 primer pairs; 3 replicates/sample/primer pair) generated a total of about 8.6 million (M) paired-end reads (per base read quality plots available in Fig.S2). The number of reads, variants, replicates and samples that were validated by each filtering step is reported in Table1, for the whole data set and for T pos1 only (as a sample example). After the initial quality filtering and assignment steps (steps 1 and 2), the number of reads per replicate varied between 137 and 64,718 (median: 26,943) for faecal samples, 21,613 and 32,490 (me-mock samples nor in the faecal samples. At the end of the HTS filtering process, all negative controls were eliminated, and 38 of the 40 initial faecal samples were represented by at least one variant or contig (samples 14Ben05 and 14Ben09 were not validated). Seven variants/contigs were identified in the mock samples (T pos1 and T pos2 ): six that correspond to the six organisms used

	cate per sample (Step 5), one sample for MFZR (14Ben09) and two
	samples for ZFZR (14Ben09, 14Deo04) were discarded. The Obiclean
	step (Step 6) reduced further the number of variants (by 32% for
	MFZR and 35% for ZFZR). When applying the Renkonen filter (Step
	7) and eliminating variants that occur in only one replicate per sam-
	ple (Step 8), one MFZR sample (14Ben05) and two ZFZR samples
	(14Ben02, 14Ben05) were further discarded. UCHIME2 detected four
	MFZR and 17 ZFZR variants as chimeras (Step 10) and four MFZR
	and two ZFZR variants were potential pseudogenes (Step 11). The
	final 81 MFZR and 61 ZFZR variants represented <0.3% of the num-
	ber of variants validated at Step 2, but they corresponded to over
	70% of the reads initially assigned to samples.

dian: 26,943) for positive controls, and 16 and 230 (median: 38) for negative controls in the MFZR data set. For ZFZR, read counts per replicate varied between

83 and 76,279 (median: 26,603)

for faecal

samples, 19,144 and 42,281 (median: 36,645)

for positive controls, and 15 and 167 (median: 48) for negative controls. The total number of variants excluding singletons was 27,921 and 23,619 for MFZR and ZFZR, respectively (Table

1

). During the filtering stage, the three LFN filters drastically reduced the number of variants: 99.2% and 99.5% of the variants were discarded. Nevertheless, the remaining variants still represented 75% (MFZR) and 78% (ZFZR) of the reads present before the LFN filtering (Step 4). All negative control replicates for both markers were eliminated at this step. Additionally, three P. microps samples (P1, P2 and P4) were retained for MFZR data set only. After eliminating the variants present in only one repli-After the contigation of the MFZR and ZFZR variants, 38 distinct contigs and 66 distinct variants (43 amplified by MFZR only and 23 amplified by ZFZR only) were retained. A total of 93 variants and contigs were found in the Z. asper faecal samples and six in the P. microps samples. The DNA of Z. asper was neither detected in the

  restituting the composition of our mock samples. In fact, the bioinformatics pipeline we developed retained all the expected variants and filtered out all but one nonexpected variants in mock controls (this latter variant was likely introduced during construction of the mock samples; see "Results"), as well as all variants in negative controls. The final data set contained only 0.3% of the original variants, but these variants represented over 70% of the reads, reinforcing that the eliminated variants could be considered to be noise.

	4.2 | Towards a quantitative approach: the MNI
	Due to PCR biases, the number of reads is a very poor estimator of
	abundance in metabarcoding studies (Elbrecht & Leese, 2015), and
	only presence/absence of MOTUs can be obtained with clustering-
	based approaches, since different alleles coming from the same taxon
	are confounded in the same MOTU. Alternatively, several studies
	highlighted the pertinence of determining the number of distinct
	sequences belonging to the same taxon when assessing genetic diver-
	sity (Gonz alez-
	and/or by
	the use of site occupancy-detection models (that can account for
	the presence of false positives; e.g., Lahoz-Monfort, Guillera-Arroita,
	& Tingley, 2016).
	Recently, new clustering-based methods have been developed
	for denoising HTS reads (e.g., UNOISE2: Edgar, 2016b; Swarm: Mah e,
	Rognes, Quince, De Vargas, & Dunthorn, 2015). These methods

avoid clustering the reads based on fixed similarity level and their outcome depends on a very few parameters. We therefore used a modified version of our pipeline, where our denoising steps (LFN and Obiclean, steps 4 and 6) were replaced by UNOISE2 and UNCROSS. This filtering approach retained 12 unexpected variants in the mock community samples, showing that our pipeline is more accurate for

DATA ACCESSIBILI TY

The lcDNA sequences were deposited in GenBank (GenBank ID: MF458551-MF458851). Supplementary data deposited in Dryad (https://doi.org/10.5061/dryad.f40v5) included: (i) custom COI database used during HTS filtering named COI-filtering-DB; (ii) Taxassign-DB; (iii) Perl script with the .csv files indicating the samples/tag combination correspondence; (iv) unfiltered HTS data; and (v) two alignments of COI sequences used during phylogenetic analysis.


Abdeldjalil Aïssa-El-Bey

Karim Abed-Meraim

Yves Grenier
email: grenier@tsi.enst.fr

UNDERDETERMINED AUDIO SOURCE SEPARATION USING FAST PARAMETRIC DECOMPOSITION

In this paper, we consider the problem of underdetermined blind source separation using modal decomposition. Indeed, audio signals and, in particular, musical signals can be well approximated by a sum of damped sinusoidal (modal) components. Based on this representation, we propose a two steps approach consisting of a signal analysis (extraction of the modal components) followed by a signal synthesis (pairing of the components belonging to the same source) using vector clustering. Our contributions in this paper are a new separation method with relaxed assumption and reduced computational cost compared to other existing algorithms. Simulation results are given to assess the performance of the proposed algorithm.

INTRODUCTION

The objective of blind source separation (BSS) is to extract the original source signals from their mixtures using only the information within the observed mixtures with no, or very limited knowledge about the source signals and the mixing matrix. BSS problem arises in many fields, such as noise reduction, radar and sonar processing, speech enhancement, separation of rotating machine noises, biomedical signal processing and even in optical tracking system [START_REF] Nandi | Blind estimation using higher-order statistics[END_REF]. This problem has been intensively studied in the literature and many effective solutions have been proposed so far [START_REF] Nandi | Blind estimation using higher-order statistics[END_REF]. In the particular case where the number of sources is larger than the number of observed mixtures (underdetermined BSS case (UBSS)), the separation can be achieved only if side information about the sources is available (sparseness, W-disjointness, finite alphabet sources, etc). In the case of non-stationary signals (including the audio signals), certain solutions using time-frequency (TF) analysis of the observations and the sources TF-orthogonality exist for the underdetermined case [START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF][START_REF] Yilmaz | Blind separation of speech mixtures via time-frequency masking[END_REF]. In this paper, we propose an alternative approach using modal decomposition (MD) of the received signals [START_REF] Aïssa-El-Bey | Blind separation of audio sources using modal decomposition[END_REF]. More precisely we propose to decompose the signal into its various modes. The audio signals and more particularly the musical signals can be modeled by a sum of damped sinusoids [START_REF] Boyer | Audio modeling based on delayed sinusoids[END_REF] and hence are well suited for our separation approach. We propose here to exploit this last property for the separation of audio sources by means of modal decomposition. To start, we review the MD-UBSS approach presented in [START_REF] Aïssa-El-Bey | Blind separation of audio sources using modal decomposition[END_REF], then we propose an improved algorithm that reduces the computational cost and relax some of the working assumptions. In this paper, we use bold upper and lower case letters for matrices and vectors, respectively. The remaining notational conventions and major symbols are listed as follows:

(•) * Complex conjugation. (•) T Transpose. (•) H Transpose conjugate.
• Frobenius norm. I Identity matrix.

DATA MODEL

The blind source separation model assumes the existence of N independent signals s1(t), . . . , sN (t) and M observations x1(t), . . . , xM (t) which represent the mixtures. These mixtures are supposed linear and instantaneous, i.e.

xi(t) = N X j=1 aij sj(t) i = 1, . . . , M . (1)
This can be represented compactly by the mixing equation

x(t) = As(t) (2)
where s(t)

. , aMi]

T contains the mixture coefficients. We assume that for any pair (i, j) with i = j, the vectors ai and aj are linearly independent. It is known that BSS is only possible up to some scaling and permutation [START_REF] Cardoso | Blind signal separation: statistcal principles[END_REF]. We take advantage of these indeterminacies to assume without loss of generality that the column vectors of A are of unit norm i.e. ai = 1 for 1 ≤ i ≤ N . The considered source signals are supposed to be decomposable in a sum of modal components c j i (t), i.e:

si(t) = l i X j=1 c j i (t) t = 0, . . . , T -1 . (3)
The usual source independence assumption is replaced here by a quasi-orthogonality assumption of the modal components, i.e.

c j i |c j ′ i ′ c j i c j ′ i ′ ≈ 0 for (i, j) = (i ′ , j ′ ) (4)
where

c j i |c j ′ i ′ def = T -1 X t=0 c j i (t)c j ′ i ′ (t) * (5)
and

c j i 2 = c j i |c j i . (6)

MD-UBSS ALGORITHM

Based on the previous model, we propose an approach in two steps consisting of:

• An analysis step: in this step, one applies an algorithm of modal decomposition to the sensor outputs in order to extract all the harmonic components from them. • A synthesis step: this is to group together the modal components corresponding to the same source in order to reconstitute the original signal. This is done by observing that all modal components of a given source signal 'live' in the same spatial direction. Therefore, the proposed clustering method is based on the component's direction evaluated by correlation of the extracted (component) signal with the observed antenna signal.

Parametric signal analysis

The source signal and hence the observations are modeled as sum of damped sinusoids:

x k (t) = ℜe ( L X l=1 α l,k z t l ) (7)
where α l,k represents the complex amplitude and z l = e d l +ω l is the l th pole where d l is the negative damping factor and ω l is the angular-frequency.

× (T -D) data Hankel matrix [H(x k )]n 1 n 2 def = x k (n1 + n2), D being a window parameter chosen in the range T /3 ≤ D ≤ 2T /3.
We use Kung's algorithm given in [START_REF] Kung | Spacetime and singular-value decomposition based approximation methods for the harmonic retrieval problem[END_REF] that can be summarized in the following steps:

1. Form the data Hankel matrix H(x k ).

Estimate the

↓ Ψ = U (L) ↑ ⇔ Ψ = U (L)# ↓ U (L) ↑ (8)
where Ψ = Φ∆Φ -1 , Φ being a non-singular 2L × 2L matrix and ∆ = diag(z1, z * 1 , . . . , zL, z * L ). () # denotes the pseudo-inversion operation and arrows ↓ and ↑ denote respectively the last and the first row-deleting operator. 4. Estimate the poles as the eigenvalues of matrix Ψ. 5. Estimate the complex amplitudes by solving the least squares fitting criterion

min α x k -Zα 2 ⇔ α = Z # x k (9)
where

x k = [x k (0) . . . x k (T -1)
] T is the observation vector, Z is a Vandermonde matrix constructed from the estimated poles and α is the vector of complex amplitudes.

Signal synthesis using vector clustering

For the synthesis of the source signals one observes that thanks to the quasi-orthogonality assumption, one has:

x|c j i c j i 2 def = 1 c j i 2 2 6 4 x1|c j i . . . xM |c j i 3 7 5 ≈ ai
where ai represents the i th column vector of A. We can then associate each component b c k j to a space direction (vector column of A) that is estimated by

b a k j = x|b c k j b c k j 2 .
Two components of a same source signal are associated to the same column vector of A. Therefore, we propose to gather these components by clustering the vectors b a k j into N classes 1 . One will be able to rebuild the initial sources up to a constant by adding the various components within a same class.

Existing MD-UBSS algorithm [4]

This algorithm applies the previous analysis and synthesis steps to each signal output x k (t). By doing so, one obtains M estimates of each source signal with an estimation quality varying significantly from one sensor to another. Indeed, this latter depends strongly on the mixing matrix coefficients and, in particular, on the signal to interference ratio (SIR) of the desired source. Consequently, a blind selection method to choose a 'good' estimate among the M available is proposed in [START_REF] Aïssa-El-Bey | Blind separation of audio sources using modal decomposition[END_REF]. First, the source estimates are paired together by associating each source signal extracted from the first sensor to the (M -1) signals extracted from the (M -1) other sensors that are maximally correlated with it. The correlation factor of two signals s1 and s2 is evaluated by s 1 |s2 s 1 s 2 . Once, the source pairing achieved, the source estimate of maximal energy is selected, i.e.

b si(t) = arg max b s j i (t) {E j i = T -1 X t=0 |b s j i (t)| 2 , j = 1, • • • , M } (10)
where E j i represents the energy of the i th source extracted from the j th sensor b s j i (t).

Proposed MD-UBSS algorithm

We propose here to improve the previous algorithm w.r.t the computational cost and the estimation accuracy when Assumption 4 (Equation ( 4)) is poorly satisfied 2 . First, in order to avoid repeated estimation of modal components for each sensor output, we use all the observed data to estimate (only once) the poles of the source signals. Hence, we apply the ESPRIT-like technique on the averaged data covariance matrix H(x) define by:

H(x) = M X i=1 H(xi)H(xi) H (11)
and we apply steps 1 to 4 of Kung's algorithm described in Section 3.1 to obtain all the poles zi, i = 1, . . . , L. This way, we reduce significantly the computational cost and avoid the problem of 'best source estimate' selection of the previous algorithm. Now, to relax Assumption 4, we can re-write the data model as:

Γz(t) = x(t) (12)
where

Γ def = [γ 1 , γ * 1 , • • • , γ L , γ * L ], γ i = βie φ i bi,
where bi is a unit norm vector representing the spatial direction of the i th component (i.e. bi = a k / a k if the i th component belongs to the k th source signal) and z(t

) def = [z t 1 , (z * 1 ) t , • • • , z t L , (z * L ) t ] T .
The estimation of Γ using the least-squares fitting criterion leads to : min

Γ X -ΓZ 2 ⇔ Γ = XZ # (13) where X = [x(0), • • • , x(T -1)] and Z = [z(0), • • • , z(T - 1)
]. After estimating Γ, we estimate the phase of each pole as:

φi = arg(γ H 2i γ 2i-1 ) 2 (14)
The spatial direction of each modal component is estimated by:

b vi = γ 2i-1 e -φ i + γ 2i e φ i = 2βibi . ( 15
)
Finally, we group together these components by clustering the vectors b vi into N classes. After clustering, we obtain N classes with N centroids b a1, . . . , b aN corresponding to the estimates of the column vectors of the mixing matrix A. If the pole zi belongs to the j th class, then according to (15), its amplitude can be estimated by:

βi = b a H j b vi 2 . ( 16
)
One will be able to rebuild the initial sources up to a constant by adding the various modal components within a same class Ci as follow:

b si(t) = ℜe

NON-DISJOINT SOURCES CASE

We consider here the case where a given component c k j (t) can be shared by several sources. This is the case, for example, for certain musical signals such as those treated in [START_REF] Rosier | Unsupervised classification techniques for multipitch estimation[END_REF]. To simplify, we suppose that a component belongs to at most two sources. Thus, let us suppose that the component c k j (t) is present in the sources sj 1 (t) and sj 2 (t) with the amplitudes αj 1 and αj 2 , respectively. It follows that the spatial direction associated with this component as estimated by (15), is given by:

b vi ≈ αj 1 aj 1 + αj 2 aj 2 . (18)
It is now a question of finding the indices j1 and j2 of the two sources associated with this component, as well as the amplitudes αj 1 and αj 2 . With this intention, one proposes an approach based on subspace projection. Let assume that M > 2 and matrix A is known and satisfies the condition that any triplet of its column vectors are linearly independent. Consequently, we have:

P ⊥ e A vi = 0, (19)
if and only if e A = [aj 1 aj 2 ], e A being a matrix formed by a pair of column vectors of A and P ⊥ e A represents the matrix of orthogonal projection on the orthogonal range space of e A, i.e.

P ⊥ e A = I -e A " e A H e A " -1 e A H . (20)
In practice, by taking into account the noise, one detects the columns j1 and j2 by minimizing:

(j1, j2) = arg min (l,m) n P ⊥ e A b vi | e A = [a l am] o . (21)
Once e A found, one estimates the weightings αj 1 and αj 2 by:

» αj 1 αj 2 - = e A # b vi. (22)
In the simulation, the optimization problem of (21) is solved using exhaustive search. This is computationally tractable for small vector array sizes but would be prohibitive if M is very large. In this paper, we treated all the components as being associated to two source signals. If ever a component is present only in one source, one of the two coefficients estimated in (22) should be zero or close to zero. Also, in what precedes, the mixing matrix A is supposed to be known. This means, it has to be estimated before applying a subspace projection. This is performed here by clustering all the spatial direction vectors in (15) as for the preview MD-UBSS algorithm. Then, the i th column vector of A is estimated as the centroid of Ci assuming implicitly that most modal components belong mainly to one source signal. This is confirmed by our simulation experiment shown in Figure 2.

SIMULATION

We present here some simulation results to illustrate the performance of our blind separation algorithms. For that, we consider a uniform linear array with M = 3 sensors receiving the signals from N = 4 audio sources. The angles of arrival of the sources are chosen randomly. The sample size is set to T = 10000 samples (the signals are sampled at a rate of 22 kHz). The observed signals are corrupted by an additive white noise of covariance σ 2 I (σ 2 being the noise power). The separation quality is measured by the normalized mean squares estimation errors (NMSE) of the sources evaluated over 200 Monte-Carlo runs and defined as:

N M SEi def = 1 Nr Nr X r=1 min α " αb si,r -si 2 si 2 « (23) N M SEi = 1 Nr Nr X r=1 1 - " b si,rs H i b si,r si « 2 (24)
N M SE = 1 N N X i=1 N M SEi . (25)
where si def = [si(0), . . . , si(T -1)], b si,r ( defined similarly) represents the r th estimate of source si and α is a scalar factor that compensate for the scale indeterminacy of the BSS problem. In Figure 1, we compare the separation performance obtained by existing MD-UBSS algorithm and the new MD-UBSS algo-rithm. We observe a performance gain in favor of the new MD-UBSS due mainly to the fact that it does not rely on the quasiorthogonality assumption. This plot also highlights the problem of 'best source estimate' selection related to the MD-UBSS as we observe a performance loss between the results given by the proposed energy-based selection procedure and the optimal3 one using the exact source signals. performance of the mixing matrix A using proposed clustering method. The observed good estimation performance translates the fact that most modal components belong 'effectively' to one single source signal. In Figure 3, we compare the performance of new MD-UBSS algorithm and the same algorithm with subspace projection. One can observe that using the subspace projection leads to a performance gain at moderate and high SNRs. At low SNRs, the performance is slightly degraded due to the noise effect. Indeed, when a given component belongs 'effectively' to

  def = [s1(t), . . . , sN (t)] T is a N × 1 column vector collecting the source signals, vector x(t) similarly collects the M observed signals, and the M × N mixing matrix A def = [a1, . . . , aN ] with ai = [a1i, . .

  ℜe(•) represents the real part of a complex entity. For the extraction of the modal components, we propose to use the ESPRIT-like (Estimation of Signal Parameters via Rotation Invariance Technique) technique that estimates the poles of the signals by exploiting the row-shifting invariance property of the D

  2L-dimensional signal subspace U (L) = [u1 . . . u2L] of H(x k ) by means of the SVD (u1 . . . u2L are the principal left singular vectors of H(x k )). 3. Solve (in the least squares sense) the shift invariance equation U (L)

Fig. 1 .

 1 Fig. 1. NMSE versus SNR for 4 audio sources and 3 sensors: comparison of the performance of MD-UBSS algorithms with and without quasi-orthogonality assumption.

Figure 2 Fig. 2 .

 22 Fig. 2. Mixing matrix estimation: NMSE versus SNR for 4 speech sources and 3 sensors.

Fig. 3 .

 3 Fig. 3. NMSE versus SNR for 4 audio sources and 3 sensors: comparison of the performance of modified MD-UBSS algorithm and the same algorithm with subspace projection.

In the simulation, we have used the k-means algorithm in[START_REF] Frank | The data analysis handbook[END_REF] for vector clustering.

This is the case when the modal components are closely spaced or for modal components with strong damped factors.

Clearly, the optimal selection procedure is introduced here just for performance comparison and not as an alternative selection method since it relies on the exact source signals that are unavailable in our context. only one source signal, equation (22) would provide a non zero amplitude coefficient for the second source due to noise effect which explains the observed degradation.

CONCLUSION

This paper introduces a new MD-UBSS algorithm of audio sources. The main advantages over the proposed MD-UBSS algorithm are, a reduced the computational cost and relaxed a quasi-orthogonality assumption. Moreover, this algorithm is extended to the non-disjoint sources case using an approximate subspace projection technique. Simulation results illustrate the effectiveness of our algorithm compared to the one [4].


N Linh-Trung
email: linh-trung.nguyen@cnes.fr

A Aïssa-El-Bey

K Abed-Meraim

A Belouchrani
email: adel.belouchrani@enp.edu.dz

UNDERDETERMINED BLIND SOURCE SEPARATION OF NON-DISJOINT NONSTATIONARY SOURCES IN THE TIME-FREQUENCY DOMAIN

This paper considers the blind separation of nonstationary sources in the underdetermined case, in which we have more sources than sensors. A recently proposed algorithm applied time-frequency distributions (TFDs) to this problem and gave good separation performances in the case where sources were disjoint in the time-frequency (TF) plane. However, in the non-disjoint case, the method simply relied on the interpolation at the intersection TF points implicitly performed by a TF synthesis algorithm, instead of directly treating these points. In this paper, we propose a new algorithm that combines the abovementioned method with subspace projection in order to explicitly treat non-disjoint sources. Another contribution of this paper is the estimation of the mixing matrix in the underdetermined case.

INTRODUCTION

Blind source separation (BSS) considers the recovery of unobserved original sources from several mixtures which are observed at the output of a set of sensors. Each mixture contains a combination of the sources that is resulted from the mixing medium between the sources and the sensors. The term "blind" indicates that no a priori knowledge of both the sources and the medium structure is available. To compensate for this lack of information, the sources are usually assumed to be statistically independent. BSS has application in different areas, such as communication, speech and image processing and biomedical engineering [START_REF] Nandi | Blind Estimation Using Higher-Order Statistics[END_REF].

A challenging problem of BSS arises when there are more sources than sensors; this is now referred to as underdetermined blind sources separation (UBSS). A TFbased UBSS algorithm has been recently proposed in [START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF] to successfully separate nonstationary sources using TFDs. This algorithm gave good performances when the sources were disjoint in the TF plane. It also provided the separation of TF quasi-disjoint sources, that is the sources were allowed to have a small degree of overlapping in the TF plane.

However, the intersection TF points were not directly treated. In particular, a point at the intersection of two sources is clustered "by chance" to belong to one of the sources. As a result, the source that has picked up this point now has some information of the other source while the later lost some information of its own. The lost information can be fortunately recovered to some extent by the interpolation at the intersection point using TF synthesis. However, for the other source, there is an interference at this point, hence the separation performance may degrade if no treatment is provided for this. The more the number of intersection points not to be treated, the worse the interpolation result and the interference would become, hence the worse the final separating performance.

In this paper, we propose another algorithm that combines the above TF-UBSS algorithm with subspace projection, offering an explicit treatment of the intersection points. The main assumption used in this algorithm is that the number of sources simultaneously present at any intersection point must be smaller than the total number of sensors.

DATA MODEL AND ASSUMPTIONS

Let N -dimensional vector s(t) = [s 1 (t), . . . , s N (t)]

T represent N nonstationary source signals. The source signals are transmitted through a medium so that an array of M linear sensors picks up a set of mixed signals represented by an M -dimensional vector x(t) = [x 1 (t), . . . , x M (t)]

T . Consider the instantaneous mixing model, as given by:

x(t) = As(t) + w(t), (1)
where A = [a 1 , . . . , a N ] is the mixing matrix and w(t) = [w 1 (t), . . . , w M (t)] T is the observation noise vector. The goal of BSS is to recover s(t) from x(t). When M < N , the problem becomes UBSS. In this case, we assume that any M column vectors of A are linearly independent.

In [START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF], sources are assumed to be disjoint in the TF plane; that is, for any pair of two sources s 1 (t) and s 2 (t) whose TF signatures are denoted by Ω 1 and Ω 2 respectively, we have Ω 1 ∩ Ω 2 = ∅. The notion of TF disjoint is illustrated in Fig. 1-a.

In this paper, we relax the above assumption by allowing the sources to be generally non-disjoint, but limiting the degree of disjoint to the following two conditions. First, there are at most (M -1) sources present at any TF point. This allows us to apply the subspace projection approach for the source TFDs estimation, to be shown in Section 4. Second, for each source signal, there exists a 00000 00000 00000 00000 00000 region in the TF plane where the source exists alone; in other words, the energy of the other sources are negligible at the TF points within the considered region. This is needed for the estimation of the mixture matrix, also to be shown in Section 4. In the case the latter condition is not satisfied, one can always estimate the mixing matrix using algebraic [START_REF] Lathauwer | ICA techniques for more sources than sensors[END_REF] or other approaches. The TF non-disjoint is illustrated in Fig. 1-b.

SOURCE SEPARATION USING TFD

In this section, we briefly review the TF-UBSS method proposed in [START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF]. This method uses the spatial TFD (STFD) that is defined in [START_REF] Belouchrani | Blind source separation based on time-frequency signal representations[END_REF] by the following:

D xx (t, f ) = +∞ l=-∞ +∞ k=-∞ φ(k, l)• x(t + k + l)x H (t + k -l)e -j4πf l , (2)
where φ(k, l) is the TFD time-lag kernel and the superscript ( H ) denotes the complex conjugate transpose operator. The matrix D xx (t, f ) varies with respect to t and f . When evaluated at a particular TF point, its (i, j)-element has the value of the cross-TFD of x i (t) and x j (t) at this point.

Applying [START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF] to the linear data model in (1), while assuming that the additive noise is not present, leads to the following expression:

D xx (t, f ) = AD ss (t, f )A H , ( 3
)
where D xx (t, f ) is now the mixture STFD matrix and D ss (t, f ) is the source STFD matrix. After computing the above mixture STFD matrix (we did use here the Wigner-Ville distribution (WVD) or the Modified-WVD (MWVD) [START_REF] Boashash | Time Frequency Signal Analysis and Processing: Method and Applications[END_REF]) and using noise thresholding, one proceeds then to the selection of auto-source TF points; these are the points at each of which only one source is active. This selection is done following the procedure in [START_REF] Belouchrani | Blind separation of nonstationary sources[END_REF].

The structure of the mixture STFD is as such: at an auto-source TF point, there is only one value not equal to zero among its diagonal elements. Therefore, if all sources are disjoint in the TF plane then, for all TF points belong to the TF signature Ω i of a source s i (t), we have

D xx (t, f ) = D s i s i (t, f )a i a H i , (4)
where D s i s i (t, f ) is the TFD of s i (t). In [START_REF] Lathauwer | ICA techniques for more sources than sensors[END_REF], a i represents the principal eigenvector of D xx (t, f ) and

D s i s i (t, f ) is
(up to a constant) the corresponding eigenvalue, or equivalently the trace value, of D xx (t, f ). Therefore, all the auto-source points that belong to one particular source must have the same spatial direction. In other words, if we cluster these points into classes corresponding to different spatial directions, then these classes represent the individual source signals.

Given the class C i representing the source s i (t), the TFD estimate of s i (t) is then computed (up to a constant factor) by:

Dsisi (t, f ) = = trace{D wvd xx (t a , f a )}, ∀(t a , f a ) ∈ C i 0, otherwise . (5)
Having obtained the source TFD estimates Ds i s i , we then use an adequate source synthesis procedure to estimate the source waveform s i (t). The recovery of the source waveform from its TFD is made possible thanks to the following inversion property of the WVD [START_REF] Boashash | Time Frequency Signal Analysis and Processing: Method and Applications[END_REF] x

(t) = 1 x * (0) ∞ -∞ ρ wvd x ( t 2 , f ) e j2πf t df , ( 6
)
which implies that the signal can be reconstructed to within a complex exponential constant e jα = x * (0)/|x(0)| given that |x(0)| = 0. The method in [START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF] uses the synthesis algorithm that is proposed in [START_REF] Boudreaux-Bartels | Timevarying filtering and signal estimation using Wigner distributions[END_REF]. This algorithm recovers the source waveform from its WVD estimates.

In brief, the TF-UBSS algorithm in [START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF] can be summarized in the following four steps:

Tab.1: TF-UBSS algorithm S1: STFD computation and noise thresholding S2: Auto-source TF point selection S3: Clustering and source TFD estimation S4: Source signal synthesis

TF-UBSS USING SUBSPACE PROJECTION

As explained in the introduction, the TF-UBSS method reviewed in Section 3 did not treat the intersection points properly in the case where the sources are non-disjoined in the TF plane. We propose here to use an appropriate subspace projection to estimate the TFDs of the individual sources from the selected auto-source points, under the assumption that there are at most (M -1) sources simultaneously present at any given point.

In particular, consider an auto-term point (t a , f a ) at which sources s i 1 , . . . , s i I contribute (I < M ) and define s = [s i 1 , . . . , s i I ] T and A = [a i 1 , . . . , a i I ] T . We have then D xx (t a , f a ) = AD s s (t a , f a ) A H [START_REF] Belouchrani | Blind separation of nonstationary sources[END_REF] and consequently (assumed that D s s (t a , f a ) is of full rank)

Range(D xx (t a , f a )) = Range( A).
Let us now assume that A is known (or already estimated), and proceed with the estimation of the source TFDs. One can obtain information about contributing sources at this point by observing that:

P ⊥ A a i j = 0, for j = 1, . .

. , I P ⊥

A a i = 0, otherwise .

Above, P ⊥

A represents the orthogonal projection matrix onto the noise subspace of D xx (t a , f a ) and can be computed by

P ⊥ A = I -VV H ( 8
)
where V is the matrix of the I principal singular eigenvectors of D xx (t a , f a ). In practice, to take into account estimation noise, one detects the I column vectors of A as the vectors of A corresponding to the I smallest values of the set { P ⊥ A a i } where i = 1, . . . , N . Once A is obtained, we estimate the TFDs of I sources at point (t a , f a ) as the diagonal entries of:

A # D xx (t a , f a ) A # H ≈ D s s (t a , f a ). (9)
In the estimation of source TFDs above, we have assumed that A had been known a priori or already been estimated. Along with other existing methods, we propose here a method to estimate A by using the assumption that for each source signal s i there exists a TF region R i where

D xx (t, f ) = D sisi (t, f )a i a H i , ∀(t, f ) ∈ R i .
Based on that observation we estimate A as follows.

• First, detect the TF points belonging to the region R = N i=1 R i using the test criterion

(t, f ) ∈ R iff λ max {D xx (t, f )} trace{D xx (t, f )} -1 < ,
where is a small threshold (typically ≤ 0.1) and λ max {D xx (t, f )} denotes the maximum eigenvalue of D xx (t, f ).

• Then, for each point (t, f ) in R, estimate the vector a(t, f ) as the principal eigenvector of D xx (t, f ).

• Finally, cluster the set of vectors { a(t, f )}, where (t, f ) ∈ R, into N classes1 using any vector clustering technique existing in the literature [START_REF] Frank | The data analysis handbook[END_REF]. In this paper, we have used the algorithm k-means. The column vectors of A are estimated as the N centroids of the N clusters.

Table 2 presents a summary of the subspace projection based TF-UBSS algorithm.

Tab.2: subspace-based TF-UBSS algorithm S1: STFD computation and noise thresholding S2: Auto-source TF point selection S3: Selection of rank-1 STFD matrices and computation of their respective principal eigenvectors S4: Clustering of the previous set of vectors and estimation of column vectors of A as the centroid points S5: For all auto-source TF point, perform the subspace based TFD estimation S6: Source signal synthesis

SIMULATION RESULTS

In the simulation, we use a uniform linear array of M = 3 sensors having half wavelength spacing. It receives signals from N = 4 independent linear frequency-modulated sources, each has 256 samples, in the presence of additive Gaussian noise where the SNR=20 dB. The effect of cross-terms on the WVD representation of one the mixture (as shown in Fig. 2-(a)) was reduced by using MWVD (as shown in Fig. 2-(b)). Fig. 2-(d) displays the auto-source TF points using MWVD. Fig. 2-(c) re-displays these points but with the TFD values extracted from the WVD in Fig. 2-(a); these TFD values will be used for TFD estimation. We compare the algorithm proposed in [START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF] (now referred to as cluster-based TF-UBSS) with the algorithm proposed in this paper, now referred to as subspace based TF-UBSS. Note that Fig. e,h,k) show the estimated source TFDs using the cluster-based algorithm, whereas Fig. f,i,l) are those obtained by the subspace based algorithm.

From Fig. 3-(b,e), we can see that the intersection points between source s 1 (t) and source s 2 (t) were picked up by source s 2 (t) by the cluster-based algorithm. On the other hand, using the subspace-based algorithm, the intersection points have been redistributed to the two sources (Fig. 3-(c,f)).

In Fig. 4, we provide a statistical evaluation of the performance gain we observed in the previous experiment.

For that, we run N r = 500 Monte Carlo trials with the same sources but random noise realization at each trial. The quality of source extraction is measured by the fol-

lowing normalize MSE MSE = 1 N r N r r=1 s r -s 2 s 2 ,
where s r represents the estimate of the sources at the r-th trial. Note that for comparison we remove the scalar and permutation indeterminacies. Even though the considered sources are almost disjoint in the TF domain, we observe an improvement of the separation method thanks to subspace projection. We expect this improvement to be more significant if the sources have larger intersection region in TF plane.

CONCLUSION

This paper introduces a new approach for blind separation of non-disjoint and nonstationary sources using TFDs. The proposed method can separate more sources than sensors and provides, in the case of non-disjoint sources, a better separation quality than the method proposed in [START_REF] Linh-Trung | Separating more sources than sensors using time-frequency distributions[END_REF]. This method is based on a vector clustering procedure that estimates the mixing matrix A, and subspace projection to separate the sources at the intersection points in the TF plane.

Fig. 1 .

 1 Fig. 1. (a)-TF disjoint, (b) TF non-disjoint

Fig. 2 .

 2 Fig. 2. TFD choices and auto-source selection.

Fig. 3 .

 3 Fig. 3. TFD estimation.

Fig. 4 .

 4 Fig. 4. MSE versus SNR.

There exist techniques that perform both the clustering and the estimation of the number of classes. For simplicity, we assumed here the number of sources known.