domain: sdv sde sdv.bv

text:
Charlie Pauvert

Marc Buée

Valerie Laval

Véronique Edel-Hermann

Laure Fauchery

Angelique Gautier

Isabelle Lesur

Jessica Vallance

Corinne Vacher
email: corinne.vacher@inra.fr

Valérie Laval

Véronique Edel -Hermann

Angélique Gautier

Bioinformatics matters: The accuracy of plant and soil fungal community data is highly dependent on the metabarcoding pipeline

Keywords: Bioinformatics, Environmental DNA, Fungi, Illumina MiSeq, Internal transcribed spacer (ITS), Metabarcoding, DADA2, USEARCH, VSEARCH, LULU

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

Introduction

Fungal communities associated with soil and plant tissues have a significant impact on plant fitness and ecosystem function [START_REF] Dighton | The Fungal community: its organization and role in the ecosystem[END_REF]Buée et al. 2009a;[START_REF] Rodriguez | Fungal endophytes: diversity and functional roles[END_REF][START_REF] Hacquard | Towards a holistic understanding of the beneficial interactions across the Populus microbiome[END_REF][START_REF] Vandenkoornhuyse | The importance of the microbiome of the plant holobiont[END_REF]Vacher et al. 2016a;[START_REF] Baldrian | Forest microbiome: diversity, complexity and dynamics[END_REF]. Identification of the fungal species present is a prerequisite for understanding these complex communities, but this task is challenging, due to the cryptic nature, microscopic characters and morphological variability of many fungal species [START_REF] Hibbett | Fungal systematics: is a new age of enlightenment at hand?[END_REF][START_REF] Yahr | Scaling up discovery of hidden diversity in fungi: impacts of barcoding approaches[END_REF]. Sequence-based taxonomic identification of fungal community members, or metabarcoding, has been the standard technique for the last 10 years (Buée et al. 2009b;[START_REF] Hibbett | Fungal ecology catches fire[END_REF][START_REF] Jumpponen | Massively parallel 454 sequencing indicates hyperdiverse fungal communities in temperate Quercus macrocarpa phyllosphere[END_REF][START_REF] Öpik | Large-scale parallel 454 sequencing reveals host ecological group specificity of arbuscular mycorrhizal fungi in a boreonemoral forest[END_REF][START_REF] Cordier | Spatial variability of phyllosphere fungal assemblages: genetic distance predominates over geographic distance in a European beech stand (Fagus sylvatica)[END_REF][START_REF] Hibbett | Fungal systematics: is a new age of enlightenment at hand?[END_REF][START_REF] Schmidt | Illumina metabarcoding of a soil fungal community[END_REF][START_REF] Hibbett | Sequence-based classification and identification of fungi[END_REF]). The internal transcribed spacer (ITS) region is now recognized as the universal barcode for fungi (Schoch et al. 2012), and is conventionally used to sequence fungal communities [START_REF] Lindahl | Fungal community analysis by high-throughput sequencing of amplified markers -a user's guide[END_REF][START_REF] Bálint | An Illumina metabarcoding pipeline for fungi[END_REF], in combination with large taxonomic reference databases, such as UNITE [START_REF] Abarenkov | The UNITE database for molecular identification of fungi -recent updates and future perspectives[END_REF][START_REF] Kõljalg | Towards a unified paradigm for sequence-based identification of fungi[END_REF].

Despite their widespread use, metabarcoding approaches suffer from various biases due to the sampling process, molecular biology steps and bioinformatic analyses used [START_REF] Lindahl | Fungal community analysis by high-throughput sequencing of amplified markers -a user's guide[END_REF][START_REF] Schmidt | Illumina metabarcoding of a soil fungal community[END_REF][START_REF] Bálint | Millions of reads, thousands of taxa: microbial community structure and associations analyzed via marker genes[END_REF][START_REF] Sommeria-Klein | Inferring neutral biodiversity parameters using environmental DNA data sets[END_REF][START_REF] Palmer | Non-biological synthetic spike-in controls and the AMPtk software pipeline improve mycobiome data[END_REF]). These biases can prevent accurate recovery of the fungal community. For instance, the fungi identified may differ according to the barcode region chosen and the primers used for its amplification [START_REF] Tedersoo | Shotgun metagenomes and multiple primer pair-barcode combinations of amplicons reveal biases in metabarcoding analyses of fungi[END_REF], the sequencing platform [START_REF] Motooka | Fungal ITS1 deep-sequencing strategies to reconstruct the composition of a 26-species community and evaluation of the gut mycobiota of healthy Japanese individuals[END_REF], the method used to assemble reads [START_REF] Nguyen | Parsing ecological signal from noise in next generation amplicon sequencing[END_REF], the sequence clustering method [START_REF] Cline | Moving beyond de novo clustering in fungal community ecology[END_REF][START_REF] Halwachs | Critical issues in mycobiota analysis[END_REF]) and the filters subsequently applied to the operational taxonomic unit (OTU) table [START_REF] Bokulich | Quality-filtering vastly improves diversity estimates from Illumina amplicon sequencing[END_REF][START_REF] Brown | Scraping the bottom of the barrel: are rare high throughput sequences artifacts?[END_REF]. Fungal ecologists thus face difficult decisions at every stage in metabarcoding studies [START_REF] Alberdi | Scrutinizing key steps for reliable metabarcoding of environmental samples[END_REF]).

Many pipelines have been developed that can be used for processing fungal ITS sequence data. These pipelines include MOTHUR [START_REF] Schloss | Introducing mothur: open-source, platform-independent, communitysupported software for describing and comparing microbial communities[END_REF], QIIME [START_REF] Caporaso | QIIME allows analysis of high-throughput community sequencing data[END_REF], SCATA [START_REF] Durling | SCATA -An efficient bioinformatic pipeline for species identification and quantification after high-throughput sequencing of tagged amplicons[END_REF], CLOVR-ITS [START_REF] White | CloVR-ITS: Automated internal transcribed spacer amplicon sequence analysis pipeline for the characterization of fungal microbiota[END_REF], VSEARCH [START_REF] Rognes | VSEARCH: a versatile open source tool for metagenomics[END_REF], FROGS [START_REF] Escudié | FROGS: Find, Rapidly, OTUs with Galaxy Solution[END_REF], PIPITS [START_REF] Gweon | PIPITS: an automated pipeline for analyses of fungal internal transcribed spacer sequences from the Illumina sequencing platform[END_REF] and DADA2 [START_REF] Callahan | DADA2: High-resolution sample inference from Illumina amplicon data[END_REF]. However, the variable length of fungal ITS sequences between taxa and high levels of sequence variability render analysis and interpretation of the data particularly difficult [START_REF] Tedersoo | Shotgun metagenomes and multiple primer pair-barcode combinations of amplicons reveal biases in metabarcoding analyses of fungi[END_REF][START_REF] Halwachs | Critical issues in mycobiota analysis[END_REF][START_REF] Palmer | Non-biological synthetic spike-in controls and the AMPtk software pipeline improve mycobiome data[END_REF]. Fortunately, comprehensive guidelines have been developed, to help fungal community ecologists to make the most appropriate choices [START_REF] Lindahl | Fungal community analysis by high-throughput sequencing of amplified markers -a user's guide[END_REF][START_REF] Bálint | An Illumina metabarcoding pipeline for fungi[END_REF][START_REF] Bálint | Millions of reads, thousands of taxa: microbial community structure and associations analyzed via marker genes[END_REF]). These guidelines suggest, for example, that sequence clustering yields the best results with ITS extraction tools such as ITSx [START_REF] Bengtsson-Palme | Improved software detection and extraction of ITS1 and ITS2 from ribosomal ITS sequences of fungi and other eukaryotes for analysis of environmental sequencing data[END_REF][START_REF] Lindahl | Fungal community analysis by high-throughput sequencing of amplified markers -a user's guide[END_REF][START_REF] Bálint | An Illumina metabarcoding pipeline for fungi[END_REF]. The removal of rare OTUs, which may be artifacts, is also generally recommended [START_REF] Bálint | Millions of reads, thousands of taxa: microbial community structure and associations analyzed via marker genes[END_REF]), but there is no consensus concerning the threshold number of sequences below which an OTU can be considered rare. The proposed thresholds range from 1 to 10 sequences [START_REF] Brown | Scraping the bottom of the barrel: are rare high throughput sequences artifacts?[END_REF] or depend on the relative abundance of OTUs [START_REF] Bokulich | Quality-filtering vastly improves diversity estimates from Illumina amplicon sequencing[END_REF]. Fungal mock communities have also recently been used for the development of guidelines [START_REF] Nguyen | Parsing ecological signal from noise in next generation amplicon sequencing[END_REF][START_REF] Cline | Moving beyond de novo clustering in fungal community ecology[END_REF][START_REF] Bakker | A fungal mock community control for amplicon sequencing experiments[END_REF]. [START_REF] Nguyen | Parsing ecological signal from noise in next generation amplicon sequencing[END_REF] showed, for example, that single forward reads could be used to recover all of the 25 well-amplified species of their mock community, whereas only 23 of these species were recovered with assembled pairedend reads. The clustering algorithm of USEARCH [START_REF] Edgar | Search and clustering orders of magnitude faster than BLAST[END_REF]) has also been recommended, based on the demonstration that it recovered the expected number of mock species [START_REF] Cline | Moving beyond de novo clustering in fungal community ecology[END_REF].

There is currently no clear consensus in the scientific community concerning the most appropriate bioinformatic approach for analysis of the fungal ITS regions sequenced on Illumina MiSeq platforms. We aimed to fill this gap, by creating and sequencing a mock community of 189 Dikarya strains commonly found in agricultural and forest soils and in plant tissues. As advised by [START_REF] Nguyen | Parsing ecological signal from noise in next generation amplicon sequencing[END_REF], this mock community had a large taxonomic breadth and some genera were represented by several closely related strains (Fig. 1). We compared the ability of 360 combinations of bioinformatic softwares and parameters to recover the fungal strains present in this mock community, in the expected proportions. In particular, we investigated whether clustering-free software packages that identify exact sequence variants of amplicons (ASVs) rather than clustering similar sequences into OTUs [START_REF] Callahan | DADA2: High-resolution sample inference from Illumina amplicon data[END_REF]) outperformed conventional clustering approaches, by fully exploiting molecular barcode resolution [START_REF] Callahan | Exact sequence variants should replace operational taxonomic units in marker-gene data analysis[END_REF]. We also tested novel post-clustering curation tools (Frøslev et al. 2017). We provide new guidelines, based on our results, for researchers using metabarcoding approaches for the analysis of fungal community richness and composition.

Materials and methods

Fungal mock community

The mock community consisted of an equimolar mixture of DNA extracted from 189 pure fungal strains isolated from soils, sporocarps or plant tissues. All the strains belonged to the superkingdom Holomycota [START_REF] Tedersoo | High-level classification of the Fungi and a tool for evolutionary ecological analyses[END_REF]): 87 Ascomycota strains, 99 Basidiomycota strains and 3 Mucoromycota strains, corresponding to 181 different species, 97 genera, 67 families, 30 orders and 11 classes. Altogether, 30 genera were represented by several species and 4 species were represented by several strains (Fig. 1 and Table S1).

Fungal DNA was obtained from the inner flesh of sporocarps, or from aerial mycelium scraped aseptically from the surface of pure cultures grown on PDA (Potato Dextrose Agar), MA (Malt Agar) or Pachlewski's medium [START_REF] Martin | Nitrogen assimilation in mycorrhizas. I. Purification and properties of the nicotinamide adenine dinucleotide phosphate-specific glutamate dehydrogenase of the ectomycorrhizal fungus Cenococcum graniforme[END_REF]. The mycelium was lyophilized for 24 h in an Edwards Modulyo 4K lyophilizer (Edwards, United Kingdom) and 100 mg of lyophilized mycelium was then placed in a Fast-Prep tube (2 mL) containing 130 mg glass beads (4.5 mm in diameter; Dutscher, France) and ground with a FastPrep® machine (MP Biomedicals, France) for 30 s at maximum shaking frequency. DNA was extracted with a DNeasy Plant Minikit (Qiagen, France), in accordance with the manufacturer's instructions, except that the incubation time was extended to 1 h at 65°C, and the volumes of buffers AP1 and P3 were doubled. DNA from all strains was quantified with a Qubit® 2.0 Fluorometer (Life Technologies, USA) and pooled in an equimolar mixture. Pooling was performed three times (replicates A, B and C). The fungal ITS1 region was amplified from each replicate with the ITS1F (5'-CTTGGTCATTTAGAGGAAGTAA-3', [START_REF] Gardes | ITS primers with enhanced specificity for basidiomycetes -application to the identification of mycorrhizae and rusts[END_REF] and ITS2 (5'-GCTGCGTTCTTCATCGATGC-3', [START_REF] White | Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics[END_REF] primers. These primers are considered as universal fungal primers and are commonly used in fungal community analyses (Buée et al. 2009b;[START_REF] Cordier | Spatial variability of phyllosphere fungal assemblages: genetic distance predominates over geographic distance in a European beech stand (Fagus sylvatica)[END_REF][START_REF] Nguyen | Parsing ecological signal from noise in next generation amplicon sequencing[END_REF][START_REF] Palmer | Non-biological synthetic spike-in controls and the AMPtk software pipeline improve mycobiome data[END_REF]. They are known to amplify Ascomycota, Basidiomycota and Mucoromycota [START_REF] Bellemain | ITS as an environmental DNA barcode for fungi: an in silico approach reveals potential PCR biases[END_REF]Schoch et al. 2012). PCR was performed with a GeneAmp PCR System 2700 (Applied Biosystems, USA). The reaction mixture (20 µL final volume) consisted of 1x of PCR buffer, 0.56 mg mL -1 of bovine serum albumin (A2153-10G, Sigma, USA), 0.2 mM of each dNTP, 0.2 µM of each primer, 0.05 U µL -1 Taq DNA polymerase (D1806, Sigma-Aldrich) and 5 ng of DNA template. The following cycling parameters were then used for amplification: enzyme activation at 94°C for 3 min; 35 cycles of denaturation at 94°C for 30 s, annealing at 53°C for 30 s, extension at 72°C for 45 s, and a final extension at 72°C for 10 min. The quality of the PCR products was checked by electrophoresis on 2% agarose gels. PCR products were purified (CleanPCR, MokaScience), multiplex identifiers and sequencing adapters were added, and library sequencing on an Illumina MiSeq platform (v3 chemistry, 2x250 bp) and sequence demultiplexing (with exact index search) were performed at the Get-PlaGe sequencing facility (Toulouse, France).

Full-length ITS sequences were also obtained by Sanger sequencing for the 189 fungal strains.

PCR was performed with the ITS1F and ITS4

(5'-TCCTCCGCTTATTGATATGC-3', [START_REF] White | Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics[END_REF]) primers, with the same PCR mixture as described above. The PCR program consisted of an initial denaturation at 95°C for 3 min, followed by 30 cycles of denaturation at 95°C for 30 s, annealing at 55°C for 30 s, and elongation at 72°C for 45 s. Sequencing reactions were performed by Genewiz (Takeley Essex, UK) and sequences from both strands were assembled with MultAlin [START_REF] Corpet | Multiple sequence alignment with hierarchical clustering[END_REF]) and manually curated.

Bioinformatic approaches

We analyzed the MiSeq sequences with 360 combinations of bioinformatic softwares and parameters (referred to hereafter as bioinformatic approaches). These approaches differed in (1) the paired-end read assembly algorithm used (Assembly),

(2) the fungal ITS1 extraction method (Extraction), (3) the method of sequence variation analysis (Variation), (4) the treatment of chimeric sequences (Chimeras) and

(5) the final filtering of the community (Filtering). The various steps of the bioinformatic analyses are described below and in Figure 2.

(1) Paired-end read assembly algorithm (Assembly)

Two paired-end read assembly algorithms were compared (Fig. 2). Paired-end sequences were joined with the FASTQ-JOIN function of QIIME v1.8.0 (Assembly=FASTQ-JOIN_OL) or with PEAR v0.9.10 (Assembly=PEAR_OV) [START_REF] Caporaso | QIIME allows analysis of high-throughput community sequencing data[END_REF][START_REF] Zhang | PEAR: a fast and accurate Illumina Paired-End reAd mergeR[END_REF]. For each algorithm, three minimum overlapping lengths (OLs) between forward and reverse sequences were tested (50 bp, 100 bp and 150 bp). For the FASTQ-JOIN algorithm, no mismatch was allowed in the overlap region. We also considered the use of single forward (R1) sequences (Assembly=QUALITY_R1) (Fig. 2).

(2) Fungal ITS1 extraction (Extraction)

The ITS1 region was either extracted (Extraction=YES) from the high-quality sequences with ITSx v1.0.10 (Bengtsson-Palme et al. 2013), or not extracted (Extraction=NO) (Fig. 2). ITSx uses an alignment of conserved ribosomal genes to identify and delineate highly variable regions, such as the ITS1 region, accurately. The minimum length of the region between the binding sites for the ITS1F-ITS2 primers is about 100 bp [START_REF] Motooka | Fungal ITS1 deep-sequencing strategies to reconstruct the composition of a 26-species community and evaluation of the gut mycobiota of healthy Japanese individuals[END_REF][START_REF] Palmer | Non-biological synthetic spike-in controls and the AMPtk software pipeline improve mycobiome data[END_REF]), but the ITS1 region sensu stricto (as defined by ITSx) is shorter, as it does not include portions of the 18S and 5.8S flanking regions. We thus discarded sequences of less than 100 bp in length in cases in which the ITS1 region was not extracted, or 50 bp in cases in which it was extracted.

(3) Sequence variation analysis (Variation)

Two clustering algorithms were compared (Fig. 2). Fungal ITS1 sequences displaying more than 97% similarity were clustered into OTUs with the popular USEARCH v7.0 program (Variation=USEARCH) [START_REF] Edgar | Search and clustering orders of magnitude faster than BLAST[END_REF] or with the open-source alternative VSEARCH v2.5.2 (Variation=VSEARCH) [START_REF] Rognes | VSEARCH: a versatile open source tool for metagenomics[END_REF]. A similarity threshold of 97% was chosen as this threshold is commonly used in fungal metabarcoding studies (e.g. [START_REF] Bakker | A fungal mock community control for amplicon sequencing experiments[END_REF][START_REF] Durand | Environmental metabarcoding reveals contrasting belowground and aboveground fungal communities from poplar at a Hg phytomanagement site[END_REF]) and has been shown to perform well on a mock community [START_REF] Tedersoo | Shotgun metagenomes and multiple primer pair-barcode combinations of amplicons reveal biases in metabarcoding analyses of fungi[END_REF], despite its tendency to aggregate closely related species [START_REF] Ryberg | Molecular operational taxonomic units as approximations of species in the light of evolutionary models and empirical data from Fungi[END_REF][START_REF] Bálint | Millions of reads, thousands of taxa: microbial community structure and associations analyzed via marker genes[END_REF]. All other settings were left to default.

Sequences with a Phred score greater than 30 over 75% of the read length were included in the clustering process. Quality filtering was performed with the QIIME script split_libraries_fastq.py (Fig. 2).

We also used the R package DADA2 [START_REF] Callahan | DADA2: High-resolution sample inference from Illumina amplicon data[END_REF] to correct sequencing errors and to infer exact amplicon sequence variants (Variation=DADA2) (Fig. 2). We retained only reads with less than one expected error (given the quality scores; Edgar and Flyvbjerg 2015). Quality filtering was performed with the fastqFilter function.

Quality data were lost during the extraction step. We therefore applied DADA2 only to the fungal ITS1 sequences not extracted with ITSx. This analysis strategy is different from that recommended in the DADA2 tutorial (http://benjjneb.github.io/dada2/tutorial.html). We, therefore, also included an approach adhering to the strategy described in the tutorial (Assembly=CUTADAPT_MERGED) (Fig. 2). Primers were removed from both forward and reverse reads, with Cutadapt v1.13 [START_REF] Martin | Cutadapt removes adapter sequences from high-throughput sequencing reads[END_REF]. The forward and reverse reads were then truncated (at 200 bp and 180 bp, respectively) and we retained only reads with fewer than two expected errors (as in the default parameters of the filterAndTrim function). Reads were merged after the inference of sequence variation as described in the tutorial.

(4) Treatment of chimeric sequences (Chimeras)

The chimeric sequences identified were either removed (Chimeras=Removed) or retained in the dataset (Chimeras=Retained) (Fig. 2). Detection with QIIME script identify_chimeric_seqs.py was performed on the demultiplexed reads before USEARCH clustering (Fig. 2), as recommended in the QIIME tutorial (http://qiime.org/scripts/identify_chimeric_seqs.html). Following VSEARCH clustering, we combined de novo and reference-based strategies for chimera detection. The de novo strategy used the UCHIME [START_REF] Edgar | UCHIME improves sensitivity and speed of chimera detection[END_REF] algorithm implemented in VSEARCH. The reference-based strategy used the ITS1-only UNITE-UCHIME dataset v7.2 (as of 2017-10-10) as a reference [START_REF] Nilsson | A comprehensive, automatically updated fungal ITS sequence dataset for reference-based chimera control in environmental sequencing efforts[END_REF]. Following DADA2 sequence variation analysis, chimeric sequences were removed with the removeBimeraDeNovo function, using the consensus option.

(5) Final filtering of the community (Filtering)

Finally, we filtered the OTU and ASV tables (Fig. 2). Five filtering methods were compared: Filtering=1 involved removing OTUs (or ASVs) composed of a single sequence; Filtering=10 involved removing OTUs (or ASVs) for which less than 10 sequences were obtained when all three replicates were considered; Filtering=RA involved removing OTUs (or ASVs) with a relative abundance lower than 0.005% of the total number of sequences; Filtering=LULU used the LULU curation algorithm (Frøslev et al. 2017) to collapse erroneous OTUs (or ASVs) into their parent OTUs;

Filtering=All involved keeping all OTUs or ASVs, regardless of the number of sequences obtained. Representative sequences were assigned to taxa with the QIIME script assign_taxonomy.py, with BLAST v2.2.22 [START_REF] Altschul | Basic local alignment search tool[END_REF]) and default QIIME parameters (e-value < 0.001; identity ≥ 90%) against the local database of Sanger sequences for the fungal strains of the mock community. The LULU curation algorithm was applied with both default settings and a set of parameters adjusted to the features of the mock community. Three parameters can be tuned in LULU: the minimum sequence similarity between a 'potential daughter' and its 'potential parent' (default 84%), the minimum ratio of parent OTU abundance to daughter OTU abundance in all samples (default 1) and their minimum co-occurrence rate across samples (default 95%). Increasing the first parameter is only advised when the barcode region has little variation or when few PCR and sequencing errors are expected, and changing the second parameter is generally not recommended (Frøslev et al. 2017).

Therefore we tuned the third parameter. We lowered its value to 66.6% to account for the small number of samples in our study (3 replicate samples per bioinformatic approach).

Comparison criteria

We defined three criteria for comparisons of the ability of the bioinformatic approaches to recover the mock community: sensitivity, precision and compositional similarity.

Sensitivity and precision were defined as the true positive rate TP/(TP+FN) and the positive predictive value TP/(TP+FP), respectively, where TP is the number of truepositive OTUs (or ASVs), FN is the number of false-negative OTUs (or ASVs) and FP is the number of false-positive OTUs (or ASVs). True-positive OTUs corresponded to fungal strains present in the mock community and identified by the bioinformatic approach considered. False-negative OTUs corresponded to fungal strains present in the mock community but not detected by the bioinformatic approach considered. False-positive OTUs corresponded to all other OTUs. If several OTUs were assigned to the same fungal strain of the mock community (i.e. 'split' OTUs), only the most abundant was considered to be a true-positive OTU, the others being considered false-positive OTUs. Compositional similarity was defined as the Bray-Curtis similarity [START_REF] Odum | Bird populations of the Highlands (North Carolina) Plateau in relation to plant succession and avian invasion[END_REF] between the community recovered and the mock community. It was calculated as 1-BC, where BC is the Bray-Curtis dissimilarity obtained from the vegdist function of the R vegan package [START_REF] Oksanen | vegan: Community Ecology Package[END_REF], assuming a uniform distribution of sequences between the fungal strains in the mock community. The expected number of sequences per fungal strain in the mock community was calculated for each replicate and each bioinformatic approach as the total number of high-quality sequences (obtained after the Filtering step) divided by the total number of fungal strains in the mock community. Ribosomal RNA gene copy number information was not available for each strain and was not used to adjust the expected number of sequences.

All three criteria theoretically range between 0 and 1. They equal 1 when the algorithm successfully identifies all members of the mock community. However, maximum sensitivity may be below 1 if the sequences of some fungal strains are absent from the raw Illumina dataset. We, therefore, estimated the total number of strains present in the raw dataset, by aligning the forward and reverse MiSeq sequences with the ITS1 Sanger sequences, with a similarity threshold of 100% and an alignment length threshold of 90% of the length of the shorter sequence. Alignments were performed with VSEARCH (--usearch_global) [START_REF] Rognes | VSEARCH: a versatile open source tool for metagenomics[END_REF], using the following parameters: --id 1 -userout --userfields query+target+qcov+tcov+id --maxaccepts 20 -top_hits_only.

Results

Assessment of maximum sensitivity

The manually curated Sanger database contained the ITS1 sequences of the 189 fungal strains of the mock community (Fig. 1 Only 160 of these 175 unique ITS1 sequences were detected in the raw Illumina MiSeq dataset (Table S1), suggesting that the other 15 strains were either not amplified by the ITS1F-ITS2 primer pair, not sequenced or were sequenced with errors. Eleven of these strains were detected in the Illumina data when the similarity threshold between Sanger and Illumina sequences was lowered to 93.5% (Table S1), suggesting that their apparent absence was caused by mismatches between the Sanger sequence and the Illumina sequences. Four of these eleven strains had ambiguous bases in the Sanger sequence, preventing a perfect match with Illumina sequences. These ambiguous bases might be due to the within-strain polymorphism of the ITS region that exists for some fungi [START_REF] Fiers | Genetic diversity of Rhizoctonia solani associated with potato tubers in France[END_REF]. The other four ITS1 sequences absent from the raw Illumina dataset came from the following species: Lepiota clypeolaria, Mycena abramsii, M. galopus and Panellus stipticus. The first species was successfully amplified with the ITS1F-ITS4 primer pair before Sanger sequencing and possessed the exact sequence of the ITS2 primer, suggesting that its absence from the Illumina dataset was caused by DNA pooling biases rather than lack of amplification. In contrast, the sequence of the ITS2 primer was detected with some mismatches for the three last species. Their absence could be due to a lack of amplification by the ITS1F-ITS2 primer pair (Table S1). The maximum sensitivity attainable by any bioinformatic approach (that is, the maximum proportion of fungal strains that could actually be found) therefore ranged from 84.7% to 90.1%. The lower bound was obtained by considering that the raw Illumina dataset contained the ITS sequences of 160 strains (of 189), while the upper bound also took into account the nine strains sequenced with some errors.

Influence of Assembly, Extraction, Variation, Chimera and Filtering on sequencing data (1) Influence of Assembly

In total, we obtained 143873 paired-end Illumina reads of 250 bp each (Table S2). We obtained 43352, 56406 and 44115 sequences for the three replicates. The mean quality of the forward (R1) reads was slightly higher than that of the reverse (R2) reads (35.59 versus 34.31, respectively) (Table S2).

The choice of paired-end read assembly algorithm strongly influenced the number, length and quality of the consensus sequences (Table S2 and Figure S1). FASTQ-JOIN retained on average 53.1% of the raw reads (whatever the minimum overlapping length), whereas PEAR retained on average 96.3% of the raw reads (Table S2). Mean sequence quality was also higher for PEAR than for FASTQ-JOIN (Table S2 and Figure S1). More than 90% of assembled reads passed the quality filter, whatever the assembly algorithm used (Table S3). Thus, PEAR generated twice as many highquality assembled reads as FASTQ-JOIN (Table S3).

(

2) Influence of Extraction

The extraction of the ITS region with ITSx (Bengtsson-Palme et al. 2013) retained 97% to 99% of the assembled reads, but only 59% of the forward reads (Table S4). Reads were 118 nucleotides shorter, on average, after extraction. Most reads were between 200 and 300 bp long before extraction (Table S2), versus 100 to 200 bp after extraction (Table S4).

(

3) Influence of Variation and Chimera

The total number of OTUs (or ASVs) varied by several orders of magnitude, depending on the method used to analyze sequence variation. For example, USEARCH identified 878 non-chimeric OTUs and 71 chimeric OTUs on average with the following parameters, Assembly=QUALITY_R1 and Extraction=No. VSEARCH identified 577 non-chimeric OTUs and 315 chimeric OTUs with the same parameters. DADA2 identified 157 non-chimeric ASVs and 40 chimeric ASVs. These striking differences were found for all mock replicates (Figure S2).

(4) Influence of Filtering

The final filtering step also strongly influenced the number of OTUs. For instance, removing OTUs with less than 10 sequences (Filtering=10) reduced the number of non-chimeric OTUs identified by USEARCH from 878 to 329, and the number of nonchimeric OTUs identified by VSEARCH from 577 to 257, for Assembly=QUALITY_R1

and Extraction=No. The LULU curation algorithm reduced even more the number of non-chimeric OTUs but, in contrast to other filtering methods, it did not lose any sequence (Table S5). ASV tables were more robust than OTU tables to variations in the filtering methods (Table S5).

Comparison of the bioinformatic approaches on the basis of sensitivity, precision and compositional similarity criteria

Bioinformatic analyses generated 360 matrices containing the number of sequences per OTU (or ASV) for the three replicates (Fig. 2). The matrices differed considerably.

For example, the mean number of OTUs (or ASVs) per replicate ranged from 57 to 1562, depending on the bioinformatic approach used (Table S6). Sensitivity, precision and compositional similarity values were calculated at replicate level. The ranking of approaches according to these criteria differed between the three mock replicates (Table S6), but the approaches that performed very well according to a given criterion for one replicate generally also performed well for the other two replicates. We, therefore, used the mean value of the criteria over the three replicates to rank the bioinformatic approaches.

The sensitivity of the 360 bioinformatic approaches ranged from 22% to 87% (Table S6). The 10 most sensitive approaches are listed in Table 1. All used USEARCH or VSEARCH to cluster sequences into OTUs and did not extract the fungal ITS1 region with ITSx. All these approaches produced very large numbers of OTUs, up to seven times more than the actual number of fungal strains in the mock community (Table 1 and Fig. 3A). They, thus, recovered most of the fungal strains of the mock community but also generated many false-positive OTUs (Fig. 3B). The precision of these approaches was, therefore, very low (Fig. 4). However, they displayed a high degree of compositional similarity to the mock community despite the large number of falsepositive OTUs (Table 1).

Precision (the proportion of OTUs (or ASVs) corresponding to true strains), ranged from 9% to 98% (Table S6). The 10 most precise approaches are listed in Table 2. All used DADA2 to identify amplicon sequence variants (Fig. 4), did not extract the ITS region with ITSx and used the LULU curation algorithm. The 8 first approaches used assembled reads as input data (Table 2). Removing chimeras with DADA2 before LULU curation appeared to be unfavourable, as it slightly reduced the sensitivity of all these top-ranking approaches (Table 2). Adjusting the minimum co-occurrence threshold of the LULU algorithm did not influence the results (data not shown). Unlike the most sensitive approaches, the most precise approaches yielded fewer ASVs than there were fungal strains in the mock community (Table 2) and produced very few false-positive ASVs (Fig. 3B). However, they did not recover all mock strains. The most precise approach, P1, recovered only 36% of the fungal strains of the mock community (Table 2 and Fig. 3B).

Compositional similarity to the mock community ranged from 0.15 to 0.396 (Table S6).

The 10 best approaches according to this criterion are listed in Table 3. Like the most precise approaches, these 10 approaches used DADA2 to identify amplicon sequence variants (Fig. 2) and did not extract the ITS region with ITSx. However, unlike the most precise approaches, most used non-assembled reads as input data and they did not use the LULU curation algorithm. The approach with the best performance according to the similarity criterion (Si1) used R1 reads as input data, retained chimeras and applied no filters to the final ASV table. This approach recovered 77.4% of the fungal strains from the mock community but had a relatively low precision (Table 3). In contrast, the Si5 approach, which used the same options but with the removal of chimeras, had a precision increased by 17%. This is because chimera removal efficiently discarded false positive ASVs, lowering their number from 51 to 13 (Fig. 3B).

As a side effect, chimera removal triggered the loss of 3 true positive ASVs, slightly reducing the sensitivity of the Si5 approach (Fig. 3B). The removal of primers (as recommended in the DADA2 tutorial) did not improve the performance of these two top-ranking approaches. It slightly lowered the precision of the Si1 approach (Table S7). The compositional similarities of the Si1 and Si5 approaches were 0.396 and 0.393, respectively (Table 3). These values were among the highest obtained, but were far from the maximal value of 1 indicating an exact match between the observed and expected community. This difference resulted from the huge variability in the number of sequences per ASV, contrasting with the expected uniform distribution of reads between fungal strains (Fig. S3). The expected number of reads for each fungal strain was then multiplied by the number of fungal strains with an identical ITS1 sequence (Fig. S3), which increased compositional similarity values (Table S6) but did not change the ranking of the bioinformatic approaches (Spearman ρ = 0.99; p < 2.2e-16).

Finally, comparison of the bioinformatic approaches revealed that some steps that are commonly recommended, such as ITS extraction and chimera removal, can have positive effects but also negative ones. For instance, the extraction of the ITS1 region with ITSx before USEARCH and VSEARCH clustering increased significantly precision but it decreased sensitivity (Fig. S4), suggesting that ITS extraction discarded some false-positive OTUs but also some true-positive OTUs. Similarly, bioinformatic approaches that kept chimeras after DADA2 sequence correction (Si1-Si4 in Table 3)

were slightly more sensitive than approaches that remove chimeras (Si5-Si8 in Table 3), indicating that chimera removal discarded some true-positive OTUs. This negative effect of chimera removal also occurred in the USEARCH and VSEARCH pipelines, but to a lower extent (Fig. S2).

Discussion

Metabarcoding approaches have revolutionized fungal ecology over the last decade [START_REF] Hibbett | Fungal ecology catches fire[END_REF]) and have become the gold standard for describing the richness and composition of communities and the networks of associations between community members [START_REF] Bálint | Millions of reads, thousands of taxa: microbial community structure and associations analyzed via marker genes[END_REF]. They have been so successful that fungal ecologists are struggling to cope with the boom in sequencing platforms, bioinformatic pipelines, taxonomic databases and community analysis tools. Benchmark studies and methodological reviews are required to help them make the most appropriate choices (e.g. [START_REF] Lindahl | Fungal community analysis by high-throughput sequencing of amplified markers -a user's guide[END_REF][START_REF] Bálint | Millions of reads, thousands of taxa: microbial community structure and associations analyzed via marker genes[END_REF][START_REF] Weiss | Correlation detection strategies in microbial data sets vary widely in sensitivity and precision[END_REF][START_REF] Pollock | The madness of microbiome: attempting to find consensus "best practice" for 16s microbiome studies[END_REF]. In this study, we focused on one aspect of the metabarcoding approach, bioinformatic analysis, assessing its effect on the recovery of community richness and composition.

We compared the ability of 360 bioinformatic approaches to recover a mock community of fungal strains commonly found in soils and plants and including 97 genera from subkingdom Dikarya. This mock community was much larger than the fungal mock communities analyzed in previous studies [START_REF] Amend | Quantifying microbial communities with 454 pyrosequencing: does read abundance count?[END_REF][START_REF] Ihrmark | New primers to amplify the fungal ITS2 region -evaluation by 454-sequencing of artificial and natural communities[END_REF][START_REF] Nguyen | Parsing ecological signal from noise in next generation amplicon sequencing[END_REF][START_REF] Taylor | Accurate Estimation of Fungal Diversity and Abundance through Improved Lineage-Specific Primers Optimized for Illumina Amplicon Sequencing[END_REF][START_REF] Cline | Moving beyond de novo clustering in fungal community ecology[END_REF][START_REF] Bakker | A fungal mock community control for amplicon sequencing experiments[END_REF]) and covered both the Ascomycota and Basidiomycota clades (Fig. 1).

We selected three criteria for comparing bioinformatic approaches: sensitivity, precision and compositional similarity to the mock community. The first two criteria are related to the number of OTUs (or ASVs) recovered and are commonly used in benchmark studies (see [START_REF] Weiss | Correlation detection strategies in microbial data sets vary widely in sensitivity and precision[END_REF]. The third takes relative abundance into account and has been used by [START_REF] Bakker | A fungal mock community control for amplicon sequencing experiments[END_REF]. We believe that this third criterion is very important, particularly if the fungal metabarcoding data are to be used to reconstruct fungal association or interaction networks for biocontrol [START_REF] Poudel | Microbiome networks: a systems framework for identifying candidate microbial assemblages for disease management[END_REF]Vacher et al. 2016b;[START_REF] Hassani | Microbial interactions within the plant holobiont[END_REF] or biomonitoring applications [START_REF] Bohan | Next-generation global biomonitoring: large-scale, automated reconstruction of ecological networks[END_REF][START_REF] Karimi | Microbial diversity and ecological networks as indicators of environmental quality[END_REF][START_REF] Derocles | Biomonitoring for the 21st Century: integrating next-generation sequencing into ecological network analysis[END_REF]. Indeed, network inference requires the most accurate possible recovery of microbial species and their abundances [START_REF] Faust | Microbial interactions: from networks to models[END_REF][START_REF] Friedman | Inferring correlation networks from genomic survey data[END_REF][START_REF] Berry | Deciphering microbial interactions and detecting keystone species with co-occurrence networks[END_REF][START_REF] Weiss | Correlation detection strategies in microbial data sets vary widely in sensitivity and precision[END_REF].

Our comparison revealed huge discrepancies between bioinformatic approaches, thereby confirming the importance of carefully selecting the most appropriate method for the analysis of fungal metabarcoding data [START_REF] Nguyen | Parsing ecological signal from noise in next generation amplicon sequencing[END_REF][START_REF] Cline | Moving beyond de novo clustering in fungal community ecology[END_REF][START_REF] Anslan | Great differences in performance and outcome of high-throughput sequencing data analysis platforms for fungal metabarcoding[END_REF]). The number of operational taxonomic units (OTUs) or amplicon sequence variants (ASVs) identified by the bioinformatic approaches compared ranged from 57 to 1562, even though there were only 189 strains in the mock community.

These results confirm that fungal community analyses should not focus on absolute values of richness estimated from metabarcoding data, but rather on the relative changes in richness between samples [START_REF] Cline | Moving beyond de novo clustering in fungal community ecology[END_REF]. The percentage of fungal strains recovered by the bioinformatic approaches ranged from 22% to 87.5%. This second value may be considered a very good result, because we estimated the maximum sensitivity attainable by a bioinformatic approach, given our data, at 90.1%.

Indeed, not all the strains in the mock community could be distinguished on the basis of their ITS1 sequences, and several strains were either not amplified at all or not accurately amplified. Our analyses revealed that four fungal species (Lepiota clypeolaria, Mycena abramsii, M. galopus, and Panellus stipticus) were absent from the sequence dataset (Fig. 1), suggesting a lack of amplification by the so-called "universal" primers [START_REF] Bellemain | ITS as an environmental DNA barcode for fungi: an in silico approach reveals potential PCR biases[END_REF][START_REF] Tedersoo | Fungal identification biases in microbiome projects[END_REF] or a sequencing failure [START_REF] Nguyen | Parsing ecological signal from noise in next generation amplicon sequencing[END_REF][START_REF] Palmer | Non-biological synthetic spike-in controls and the AMPtk software pipeline improve mycobiome data[END_REF].

DADA2, a clustering-free software package [START_REF] Callahan | DADA2: High-resolution sample inference from Illumina amplicon data[END_REF], effectively recovered the composition of the mock community. The top ten bioinformatic approaches in terms of performance for the compositional similarity criterion all used DADA2 to identify amplicon sequence variants. The total number of ASVs generated by these 10 approaches ranged from 148 to 197, which was therefore of the same order of magnitude as the total number of strains in the mock fungal community (i.e.

189). We highlighted several options for increasing the efficiency of DADA2 for fungal metabarcoding datasets. Firstly, our results confirm that the use of single forward (R1) reads as input data is a good option [START_REF] Nguyen | Parsing ecological signal from noise in next generation amplicon sequencing[END_REF]. This made it possible to ensure that strains with longer ITS regions (such as those of the genus Cantharellus, for instance; [START_REF] Feibelman | Length variation in the internal transcribed spacer of ribosomal DNA in chanterelles[END_REF] were not excluded. Based on our results, we also recommend retaining the primers for fungal communities amplified with the ITS1F-ITS2 primer pair. Indeed, we found that primer removal did not improve the recovery of mock community composition. These findings may be accounted for by the absence of degenerate nucleotides in the ITS1F-ITS2 primer pair. Primer retention may be relevant in this case, because non-degenerate primers have no impact on the denoising step of DADA2. The merging of reads after sequence variation inference, as recommended in the DADA2 tutorial (http://benjjneb.github.io/dada2/tutorial.html), did not improve the recovery of the mock community either.

The Si5 approach represented one of the best trade-offs between the three selection criteria among the 360 bioinformatic approaches compared. The Si5 approach used single forward (R1) reads as input. Quality filtering, sequence variation analysis and chimera removal were performed with DADA2 [START_REF] Callahan | DADA2: High-resolution sample inference from Illumina amplicon data[END_REF]. ITS1 extraction [START_REF] Bengtsson-Palme | Improved software detection and extraction of ITS1 and ITS2 from ribosomal ITS sequences of fungi and other eukaryotes for analysis of environmental sequencing data[END_REF]) and downstream OTU table filtering were not required.

The Si5 approach recovered the ITS1 regions of 80 out of 87 Ascomycota strains, 83 out of 99 Basidiomycota strains and all 3 Mucoromycota strains (Fig. 1), suggesting that there was no detection bias against Ascomycota strains despite the intron insert downstream of the ITS1F primer binding site that might impair their amplification (see [START_REF] Taylor | Accurate Estimation of Fungal Diversity and Abundance through Improved Lineage-Specific Primers Optimized for Illumina Amplicon Sequencing[END_REF]. We recommend the use of this simple bioinformatic approach in ecological studies of fungal communities, for the following reasons: (i) it did not overestimate the number of fungal strains, (ii) it was among the ten best bioinformatic approaches in terms of recovery of the composition of the mock community and (iii) it performed very well according to the two other criteria used for comparison (precision and sensitivity). Based on these results, the Si5 approach appears to be an appropriate bioinformatic approach for studies involving whole-community profiling and network inference.

By contrast, the clustering algorithms of USEARCH [START_REF] Edgar | Search and clustering orders of magnitude faster than BLAST[END_REF] Overall, our study highlights the importance of carefully selecting the bioinformatic approach to be used according to the objective of the metabarcoding study. Indeed, the ability of bioinformatic approaches to recover fungal strains and the relative abundances of the strains recovered varied greatly. Some approaches detected almost all strains of the mock community but overestimated community richness, whereas others retrieved the actual richness and composition of the mock community more accurately. The former are more appropriate for the detection of target species, whereas the latter are more appropriate for community ecology studies. However, none of the bioinformatic approaches compared recovered the mock community perfectly. In particular, none of the approaches found the expected distribution of sequences between fungal strains. This may be due to differences in the number of ribosomal RNA gene repeats between fungal species [START_REF] Ganley | Highly efficient concerted evolution in the ribosomal DNA repeats: Total rDNA repeat variation revealed by whole-genome shotgun sequence data[END_REF], and imperfections in equimolar pooling of DNA samples, together with biased amplification for pooled species [START_REF] Palmer | Non-biological synthetic spike-in controls and the AMPtk software pipeline improve mycobiome data[END_REF]. Because of these biases, current fungal community analyses should not focus on the within-sample distribution of taxa abundance, but rather on the changes in taxa abundance between samples. Future methodological developments should focus on reducing biases caused by molecular biology steps [START_REF] Nichols | Minimizing polymerase biases in metabarcoding[END_REF][START_REF] Porter | Scaling up: A guide to high-throughput genomic approaches for biodiversity analysis[END_REF] and on improving the bioinformatics pipelines to better recover the abundances of fungal strains. Our comparison of bioinformatics approaches could be extended, since the 360 bioinformatic approaches compared here constitute only a small fraction of the approaches that could be used to analyze fungal metabarcoding data. Other approaches may give better results, and their ranking may vary with sequence data quality [START_REF] Nguyen | Parsing ecological signal from noise in next generation amplicon sequencing[END_REF]. Future bioinformatic approach comparisons should therefore be based on multiple mock communities sequenced independently. They could also include error-correction methods alternative to that of DADA2, such as UNOISE2 [START_REF] Edgar | UNOISE2: improved error-correction for Illumina 16S and ITS amplicon sequencing[END_REF], or recent clustering approaches, such as OptiClust [START_REF] Westcott | OptiClust, an improved method for assigning amplicon-based sequence data to operational taxonomic units[END_REF] or SeekDeep [START_REF] Hathaway | SeekDeep: single-base resolution de novo clustering for amplicon deep sequencing[END_REF], or consider reference-based clustering approaches [START_REF] Cline | Moving beyond de novo clustering in fungal community ecology[END_REF][START_REF] Halwachs | Critical issues in mycobiota analysis[END_REF]; but see [START_REF] Westcott | De novo clustering methods outperform referencebased methods for assigning 16S rRNA gene sequences to operational taxonomic units[END_REF]. All the data required for the extension of our methodological comparison are provided. Illumina dataset are indicated with a cross and those that were recovered by the recommended bioinformatic approach (Si5, Table 3) are indicated with a black circle. (or ASVs) retrieved by the most sensitive approach (Se1; Table 1), the most precise approach (P1; Table 2), the approach with the best performance in terms of compositional similarity to the mock community (Si1; Table 3) and the bioinformatic approach recommended in this study (Si5; Table 3). The black horizontal line indicates highlighted with different colors and symbols. Se1 (Table 1), P1 (Table 2) and Si1

(Table 3) correspond to the most sensitive approach, the most precise approach and the approach with the best performance in terms of compositional similarity to the mock community, respectively. The bioinformatic approach recommended in this study is Si5 (Table 3).

Table 1 -List of the 10 most sensitive approaches. Sensitivity, precision and similarity values were averaged over the three replicates for each bioinformatic approach.

Richness is defined as the mean number of OTUs identified by the bioinformatic approach.

Figure 1 -

 1 Figure 1 -Taxonomic composition of the artificial fungal community. The terminal nodes of the tree are the ITS1 sequences (n = 175) of the fungal strains (n = 189) that constitute the mock community. The ITS1 sequences that were not present in the raw

Figure 2 -

 2 Figure 2 -Overview of the 360 bioinformatic approaches compared in this study. The Illumina MiSeq sequences were (1) assembled with FASTQ-JOIN (Caporaso et al. 2010) or PEAR (Zhang et al. 2014) with three minimum overlapping lengths (50 bp, 100 bp or 150 bp), or not assembled. In this latter case, single forward (R1) reads were used. After quality filtering, (2) the ITS1 region was extracted from the reads with ITSx (Bengtsson-Palme et al. 2013), or not extracted. (3) Sequence variations were then analyzed with DADA2 (Callahan et al. 2016), USEARCH[START_REF] Edgar | Search and clustering orders of magnitude faster than BLAST[END_REF] or VSEARCH[START_REF] Rognes | VSEARCH: a versatile open source tool for metagenomics[END_REF]) and (4) chimeras were either retained or removed. (*: for USEARCH, chimera detection was performed before clustering). Finally, (5) the datasets were either filtered by removing rare or erroneous OTUs (or ASVs), or left unfiltered. Filtering thresholds (T) were based on the number of sequences per OTU, or on their relative abundance (RA), or OTU curation was performed using the LULU algorithm(Frøslev et al. 2017). §: When DADA2 was used, an alternative method of read processing (CUTADAPT_MERGED) was included.

Figure 3 -

 3 Figure 3 -Richness estimates for the top four approaches. (A) Total number of OTUs

  the expected richness. (B) Total number of OTUs (or ASVs) per bioinformatic approach depending on OTU (or ASV) category. TP = true positives, FN = false negatives and FP = false positives. Results were averaged over the three replicates and rounded for clarity.

Figure 4 -

 4 Figure 4 -Values of precision and (A) sensitivity or (B) compositional similarity to the mock fungal community, for all 360 bioinformatic approaches. Each dot corresponds to the mean value obtained for an approach over the three replicates. The methods used to analyze sequence variation (DADA2, USEARCH or VSEARCH) are

  and Table S1). Several strains had identical Sanger sequences for ITS1: two strains from the genus Alternaria, six from the genus Botrytis (B. calthae, B. pseudocinerea, B. ranunculi and three strains of B. cinerea), two strains from two different species of Colletotrichum (C. destructivum and C. higginsianum), two strains of Craterellus cornucopioides, four pairs of strains from the genus Fusarium (F. acuminatum and F. avenaceum, F. langsethiae and F.

sporotrichioides, F. oxysporum and F. commune, F. verticillioides and another species of the F. fujikuroi species complex), two strains from two different species of Lepista (L. irina and L. nuda) and two strains from the genus Zymoseptoria. The Sanger database, therefore, contained 175 unique ITS1 sequences.

  and VSEARCH[START_REF] Rognes | VSEARCH: a versatile open source tool for metagenomics[END_REF] should be favored in studies in which species detection is the main goal. These clustering algorithms generally overestimated the actual number of fungal strains, but were able to retrieve almost all detectable strains. Their sensitivity[START_REF] Agtmaal M Van, Straathof | Exploring the reservoir of potential fungal plant pathogens in agricultural soil[END_REF]). On the other hand, if the purpose of a study is to focus only on fungal species present with high certainty (i.e. on a precise but incomplete community), then DADA2 and LULU(Frøslev et al. 2017) should be combined and applied to assembled sequences. The Pi3 approach, that merges reads after sequence variation inference as recommended in the DADA2 tutorial, seems to be a good compromise in this case.

	was close to the maximum value. The most sensitive approach, Se1, used single
	forward (R1) reads as input and clustered them with the VSEARCH algorithm. ITS1
	extraction (Bengtsson-Palme et al. 2013), chimera removal and downstream OTU
	table filtering were not required. In general, our comparison revealed that the steps of
	ITS extraction and chimera removal can eliminate fungal strains that are actually
	present in the community and should not be systematically used. The second most
	sensitive approach, Se2, used the USEARCH clustering algorithm. These two highly
	sensitive bioinformatic approaches are potentially useful for the early detection of
	invasive species (Comtet et al. 2015), including fungal pathogens (Munck and Bonello
	2018), for the detection of emerging pathogens accounting for the decline or death of
	host populations (Ricciardi et al. 2017), and for exploring environmental reservoirs of
	pathogens (

Acknowledgments

We thank Matthieu Barret, Martial Briand, Lucas Auer, Gregory Gambetta, Guilherme Martins, Frédéric Barraquand and Tania Fort for useful discussions on the first draft of the manuscript. We also thank the editor, Peter Kennedy and one anonymous reviewer for helpful comments of the submitted version. We are grateful to the INRA MEM metaprogramme (Meta-Omics of Microbial Ecosystems) for financial and scientific support. The mock community sequencing was funded by the INRA MEM MetaBAR project (PI: MB) and the bioinformatic analyses were performed as part of the INRA MEM Learn-biocontrol project (PI: CV). Additional funding was received from the LABEX COTE (ANR-10-LABX-45) and the LABEX CEBA (ANR-10-LABX-25-01). CP's PhD grant was funded by the INRA and Bordeaux Sciences Agro (BSA). We thank the Genotoul sequencing facility (Get-PlaGe) for sequencing the mock community and the Genotoul bioinformatics facility (Bioinfo Genotoul) for providing computing and storage resources. We also thank Julie Sappa from Alex Edelman & Associates for English language revision.

Data Availability

The raw sequence data were deposited in Dataverse and are available in the FASTQ format at https://doi.org/10.15454/8CVWRR. The code is available as an archive at https://doi.org/10.15454/VKTWKR.

Author contributions

CP performed the bioinformatic work, analyzed the results in accordance with the recommendations of IL and CV, and wrote the first draft of the article in collaboration with JV and CV. MB coordinated the design and sequencing of the mock community. MB, VL, VEH and LF provided fungal DNA for the mock community and performed the molecular biology work. AG and VL provided the Sanger sequence database. CV conceived the study in collaboration with MB, supervised the work and made a major contribution to the writing of the manuscript. All authors revised the manuscript.

Supplementary Materials

Table S1 -List of the fungal strains of the mock community and their ITS Sanger sequences.

Table S2 -Sequence summary statistics for the raw data and the assembled reads, before quality filtering.

Table S3 -Sequence summary statistics for the raw data and the assembled reads, after quality filtering.