{'fraction_non_alphanumeric': 0.030845265769373523, 'fraction_numerical': 0.025512297083471705, 'mean_word_length': 5.035369016741335, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Morphology and metabarcoding: a test with stream diatoms from Mexico highlights the complementarity of identification methods

Diatoms are among the most commonly used bioindicators. Correct taxonomic identifications are critical to their use as bioindicators because closely related diatom species can respond differently to water physicochemical characteristics and pollutants. However, diatom identification based on morphology can be time consuming, and requires highly specialized taxonomic skills. To optimize diatom identification, DNA metabarcoding is increasingly used because it is generally less time consuming and may be more accurate than morphological identification. To date, however, neither DNA metabarcoding nor DNA barcoding diatom studies have been conducted in Mexico. Thus, we studied epilithic diatoms from streams in Central Mexico with a combination of morphological and metabarcoding techniques, and compared the diatoms identified and quantified by each method. We also assembled a barcode reference library based on clonal culturing. This library is composed of 190 strains that belong to 72 species in 24 genera. The morphological analysis of environmental samples resulted in the identification of 204 infrageneric taxa in 42 genera, but clonal culturing from the same samples retrieved 12 additional infrageneric taxa and 1 additional genus, thereby revealing concealed diversity. The metabarcoding approach resulted in the identification of 266 infrageneric taxa that belonged to 35 genera. Together, these methods detected 49 genera. Of these genera, 14 were identified only by morphology, 29 were identified by both methods, and 6 were only identified by metabarcoding. Of the 266 taxa we retrieved with metabarcoding, we confidently assigned 94 infrageneric taxa because a direct morphological or barcode sequence correlation was possible. Thirty-four of these 94 taxa were only detected with the metabarcoding method. One-fourth (23) of the assignments were only possible because of the barcode reference library we developed during this study, because there were no existing barcode sequences that matched these barcodes in the International Nucleotide Sequence Database Collaboration databases. Large disparities existed between relative abundances based on valve counts and sequence reads of the most abundant taxa, even after we corrected for cell biovolume. Overall, we conclude that the combination of morphological and molecular methods increases the detection and identification of diatoms.

quality (Kelly 1998, Prygiel 2002, Potapova and Charles 2007. Taxonomic identification is a crucial step in recording accurate and meaningful biomonitoring results because closely related and cryptic species can have different environmental optima and tolerances (Morales et al. 2001, Poulíčková et al. 2008. However, taxonomic identification of diatoms based on morphology can be time consuming, needs in-depth knowledge of the diatom diversity in the study area, and requires awareness of the morphological plasticity species may display under varying environmental conditions and throughout their life-cycle (Zimmermann et al. 2015). DNA barcoding (Hebert et al. 2003) and DNA metabarcoding (Taberlet et al. 2012a) based on DNA barcoding have been proposed as alternatives to overcome the impediments to morphology-based identification of organisms.
DNA barcoding is a method used to identify species based on DNA sequences that are linked to morphologically identified specimens. Ideally, barcoding sequences are short enough to be sequenced in 1 Sanger sequencing run with universal primers and can unambiguously identify a species independent of its life-cycle stage (Hebert et al. 2003, Moritz and Cicero 2004, Zimmermann et al. 2011. DNA metabarcoding, or metabarcoding, is a method that can identify multiple species from bulk samples by relying on reference barcode sequences for species identification. DNA for metabarcoding can be extracted from bulk samples of soil, water, or air, containing a 'soup of biodiversity' (Taberlet et al. 2012b, Yu et al. 2012.
The standard metabarcoding approach consists of several steps that involve processing environmental samples (biofilm, water, soil, sediment) to obtain DNA sequences of the organisms present in those samples via High-Throughput Sequencing (HTS). This is followed by bioinformatics treatments that result in lists of species or Molecular Operational Taxonomic Units (MOTUs). These lists can then be used for applications such as biomonitoring. In most cases, however, the diversity measured by metabarcoding exceeds the diversity measured by morphological analyses (Zimmermann et al. 2015, Groendahl et al. 2017. In some cases, this reveals diversity that is indiscernible with morphological methods.
The choice of a barcoding marker that gives the desired taxonomic resolution is of critical importance in diatom metabarcoding, because most diatom studies require specieslevel resolution. In diatoms, the preferred DNA metabarcoding markers are rbcL and the 18S V4 region (Kermarrec et al. 2014, Zimmermann et al. 2015. The V4 region of the 18S rRNA gene is considered the universal barcode for protists including diatoms (Pawlowski et al. 2012), although the quest for the ideal marker for diatoms continues. The main constraint that results from DNA-based identification approaches based on barcodes is the natural intraspecific and intragenomic variability and interspecific divergence of the barcoding marker. This issue is particularly problem-atic when a single, traditionally recognized species or bioindicator taxon has a variety of genotypes at the barcoding region. This can lead to sequences that correspond to different genotypes within the same taxon clustering into different MOTUs, thereby artificially inflating the taxonomic richness relative to morphological techniques (Brown et al. 2015, Bálint et al. 2016. In diatoms, the main sources of the barcodes used to identify sequences in metabarcoding studies are clonal cultures obtained from single cell isolations. Clonal culturing has the advantage of allowing barcode sequences to be correlated with morphologically analyzed valves , Stachura-Suchoples et al. 2016. Unfortunately, however, diatom culturing is time-consuming. Single cell PCR amplifications may be an alternate way to obtain barcodes for taxa that are difficult to culture (Lang and Kaczmarska 2011, Chen et al. 2013, Hamilton et al. 2015. However, corroborating taxon identity with singlecell amplification is difficult because the valves of the isolated cells are normally destroyed in the DNA extraction process. According to Skibbe et al. (2018), however, it is sometimes possible for several genetically identical cells to be isolated, for example by isolating various cells originating from the same stalks (e.g., Gomphoneis tegelensis R. Jahn et N. Abarca), which can provide enough material for both DNA extractions and morphological examinations. A recent study proposed the use of HTS data as an additional way to obtain barcodes, and suggested criteria to ensure the barcodes truly correspond to the species observed by microscopy .
To date, neither DNA metabarcoding nor DNA barcoding studies have been conducted for diatoms in Mexico, despite the potential use of these methods in freshwater diversity studies and in biomonitoring. To compare the performance of morphology and metabarcoding in the identification and quantification of diatom abundances, our objectives were to: 1) compare the number of taxa retrieved by morphological analysis and metabarcoding of environmental samples, 2) create a regional vouchered barcode reference library to aid taxonomic assignments from sequences derived from the metabarcoding approach, 3) compare taxon abundances derived from morphology and metabarcoding approaches, and 4) test the suitability of HTS data as a source of barcode sequences.

General approach
The goal of this study was to compare the performance of morphology and metabarcoding methods in the identification and quantification of diatom species. We did this by analyzing epilithic diatoms from streams collected from Central Mexico. For the morphological approach, we used the data from light microscopy (LM) and scanning electron microscopy (SEM) observations obtained by Mora  (Proft et al. 2017) in combination with a phylogenetic-based coalescent model. We also generated the first barcode reference library of diatoms from Mexico based on clonal culturing to increase the number of reference sequences for identifying taxa to species level with the metabarcoding approach. We compared the abundances of the most common taxa obtained by morphology (LM counts) and metabarcoding (sequence read numbers). Finally, we examined the potential of retrieving barcode sequences from our HTS dataset. For the comparison of morphology and metabarcoding, we analyzed a subset of 18 samples from the 42 samples collected by Mora et al. (2017). For clonal culturing, we used all 42 samples.

Study area
The Lerma-Chapala River Basin is located in Central Mexico and covers an area of 53,590 km 2 . It lies within 2 biodiversity hotspots-Mesoamerica and the Madrean Pine-Oak Woodlands (Myers et al. 2000, Cotler et al. 2006. It is geologically and climatically heterogeneous and has well defined rainy (June-October) and dry seasons (November-May). This basin is one of the most important areas in the country for agriculture and industry and has a population of >15 million inhabitants, but it is also one of the most environmentally degraded basins in the country (Aparicio 2001, Wester et al. 2005. The studied streams (Fig. 1, Table 1) can be divided into 3 groups according to the physical and chemical composition of their waters : streams in the 1 st group have acidic waters and low specific conductivity; the 2 nd group has circumneutral waters, low specific conductivity, and the highest P concentration on average among the 3 groups; the 3 rd group is characterized by well mineralized waters with circumneutral waters and the lowest average N concentrations of the 3 groups.

Sampling
Composite epilithon samples were collected by Mora et al. (2017) by brushing cobbles across a transversal section of the streams with disposable toothbrushes, suspending the brushed material in a total volume of 60 ml. Samples were then homogenized and divided into 3 subsamples of 20 ml each, which were used for 3 different purposes: I) frozen (-247C) for HTS, II) no prior treatment for the establishment of clone cultures to build the barcode library, and III) fixed in 70% alcohol for morphological analyses.

Morphological analysis of environmental samples
Subsamples (subsample group III) were cleaned by heating with 35% hydrogen peroxide at 807C, followed by a se-ries of rinses with distilled water. Permanent slides were mounted with the high refraction index mounting medium Naphrax ® . Observations were performed by LM and SEM.
For SEM observations, cleaned material was mounted on stubs and sputter-coated with gold-palladium. Further taxa were identified to the lowest taxonomic level possible with the same identification references cited in Mora et al. (2017). Quantification of taxon abundance was done by counting a minimum of 500 valves in LM. For a detailed description of these methods see Mora et al. (2017).
Barcode reference library of the Lerma-Chapala River Basin Isolation, culturing, and harvesting of clonal cultures We isolated single-cells from aliquots of each environmental subsample (subsample group II) with micro-capillary glass pipettes under either a LM or a stereo LM. We used both high (LM) and low (stereo LM) magnification to increase the size diversity of the taxa isolated. Prior to isolating cells, we diluted the samples to decrease organism density and make it easier to isolate individual cells. We also did a series of isolation and re-isolation on microscope slides to ensure single-cell isolations. The isolated cells were then transferred to 5-cm diameter Petri dishes that contained culture medium (Alga-Gro; Carolina Biological Supply Company, Burlington, North Carolina). Most cells were grown on culture medium at the concentration recommended by the manufacturer, but we also used culture media that was either ½ or ¼ of the manufacturer-recommended concentration to increase the number of taxa that grew in our cultures (Ferris et al. 1996, Connon andGiovannoni 2002). The cultures were grown in Memmert ® Growth Chambers (Memmert, Schwabach, Germany) at 17 to 207C with a 12h day: night photoperiod. After a successful clonal culture had been established, i.e., growing at exponential phase and noncontaminated by other diatom or protist cells, the culture was divided into 3 subsamples to be used for A) DNA extraction, B) a reserve, and C) morphological analysis. Each of the subsamples was then transferred to a new Petri dish, maintained in the growth chamber for 1 to 3 wk, and subsequently harvested.

Molecular analysis of clonal cultures
We transferred the cultured material from subsample group A to 15-mL plastic centrifuge tubes for molecular analysis. We then centrifuged the tubes at 2000 rpm for 10 min, removed the supernatant, and transferred the pellets to 1.5-mL tubes. We extracted the DNA with a NucleoSpin ® Plant II Mini Kit (Macherey and Nagel, Düren, Germany) following manufacturer instructions. We checked the DNA fragment size and concentrations with gel electrophoresis (1.5% agarose gel) and NanoDrop ® (Peqlab Biotechnology LLC; Erlangen, Germany) and stored the DNA samples at 2207C for future use. We amplified the V4 region of the 18S rRNA gene, which has ∼390 to 410 base pairs (bp), with the primers and PCR regime from Zimmermann et al. (2011). We visualized PCR products in 1.5% agarose gel and cleaned them with MSB Spin PCRapace ® (Invitek LLC, Berlin, Germany) following manufacturer instructions. We measured DNA concentrations with NanoDrop ® (Peqlab Biotechnology) and normalized samples to a total DNA content >100 ng/lL for sequencing. M13 tails were used as sequencing primers following Zimmermann et al. (2014) and Ivanova et al. (2007). Sanger sequencing was conducted by Starseq ® (GENterprise LLC; Mainz, Germany). We edited the sequences in PhyDE version 0.9971 (Müller et al. 2005). The final DNA sequences were submitted to the ENA (European Nucleotide Table 1. Streams sampled in the Lerma-Chapala Basin, Central Mexico, including name, geographical coordinates, elevation (m asl) and the numbers assigned to samples collected in each stream. Samples 1 to 9 were taken during the dry season (February 2014), and samples 10 to 18 during the rainy season (September 2014).
We centrifuged subsample group B and transferred them to 1.5-mL tubes. We stored the tubes in the freezer (2247C) as reserve material in case further DNA extractions were necessary.
Morphological analysis of clonal cultures We separately transferred the cultivated material from the C subsamples to 15-mL plastic centrifuge tubes and filled them with 35% hydrogen peroxide to oxidize the organic material. After 2 d, we removed the peroxide residue by rinsing the samples with distilled water 4 Â-an initial rinse followed by 3 additional rinses every other day. We used the cleaned samples to make 1 permanent slide per subsample with Naphrax ® . We examined and photographed the diatoms on the permanent slides under a Zeiss Axio Imager M2 LM that was connected to a camera (AxioCam HRc; Zeiss, Oberkochen, Germany).
Aliquots of cleaned, cultivated material were air dried, mounted on stubs, and examined under a Hitachi FE 8010 (Hitachi, Tokyo, Japan) scanning electron microscope (SEM) operated at 1.0 kV. Taxon identification was conducted with the same identification references cited in Mora et al. (2017).

DNA metabarcoding by HTS
Subsample group I from the field were defrosted, transferred to 15-mL plastic centrifuge tubes, and centrifuged at 5000 rpm for 5 min. We then removed the supernatant and transferred the pellets to 1.5-mL tubes. We extracted the DNA with a NucleoSpin ® Plant II Mini Kit (Macherey and Nagel) following manufacturer instructions. We quantified DNA concentrations with a Qubit 2.0 fluorometer (Invitrogen, Carlsbad, California) and adjusted the volume to a concentration of 20 ng/lL. We used PCR to amplify the V4 region (18S) with the Nextera primers DIV4for: 5 0 -GCGGTAATTCCAGCTCCAATAG-3 0 and DIV4rev3: 5 0 -CTCTGACAATGGAATACGAATA-3 0 following Zimmermann et al. (2011) with a modification for 300-bp paired-end sequencing for Illumina MiSeq (Visco et al. 2015). PCR amplifications were done in duplicate for each sample, and each amplification had a total volume of 25 lL that was comprised of 0.5 lL dNTP mix (25 mM each dNTP), 0.25 lL BSA (10 mg/mL), 0.25 lL DMSO, 1 lL of each forward and reverse primers (10 pm/lL), 0.4 lL of Herculase II Fusion DNA Polymerase (Agilent Technologies Inc., Santa Clara, California), 5 lL Herculase II reaction buffer, 1 lL of template DNA (20 ng/lL), and 15.6 lL of HPLC grade water. The PCR regime included an initial denaturation at 947C (2 min), 35 cycles of denaturation at 947C (45 seconds), annealing at 527C (45 s), elongation at 727C (1 min), and a final elongation at 727C (10 min). We visualized the PCR products with electrophoresis on 1% agarose gels and pooled duplicate PCR products to make a final volume of 50 lL. Aliquots of 25 lL of the amplicons were purified with HighPrep PCR paramagnetic beads (Magbio Genomics, Gaithersburg, Maryland).
We then did a 2 nd PCR (indexing PCR) on the purified samples to ligate a unique combination of tags to the 5 0 end of the primer. The indexing PCR reactions of 25 lL were comprised of 0.25 lL dNTP mix, 1 lL DMSO, 0.625 lL of each primer, 0.25 lL of Herculase, 5 lL Herculase II reaction buffer, 10 lL of template DNA, and 7.25 lL of HPLC grade water. We started the indexing PCR regime with denaturation at 947C (2 min), 8 cycles of denaturation at 957C (20 s), annealing at 527C (30 s), elongation at 727C (30 s), and a final elongation at 727C (3 min). We purified the PCR products with HighPrep PCR paramagnetic beads and quantified them with Quant-iT PicoGreen dsDNA Assay Kits (Invitrogen). We prepared the library with MiSeq Reagent Kit V3 (Illumina, San Diego, California) following manufacturer instructions, such that samples were normalized to equal nM DNA concentrations. We then pooled the samples, denatured them to 4 nM, diluted them to 20 pM, mixed them with 5% denatured and diluted PhiX (30 lL of PhiX and 570 lL of library), and loaded them onto the MiSeq cartridge.
We used MetBaN 1.01 (Proft et al. 2017) for bioinformatic analyses. MetBaN is a bioinformatics pipeline that implements a modular and flexible phylogenetic based species delimitation approach by streamlining metabarcoding and phylogenetic software packages. Within MetBaN, we only used the 1 st modules, mostly from the OBITOOLS package (Boyer et al. 2016). We first merged the samples that consisted of paired-end reads with the illuminapairedend module, and retained only the merged reads that were >150 bp and had complete primer sequences. We then removed the primers with the ngsfilter module function and merged identical sequences to prevent redundant classifications. Singleton, chimeric, and low-quality sequences were then filtered out. Finally, we pooled all filtered sequences and clustered them into MOTUs at a 6 bp difference identity threshold. Subsequent taxonomic assignments were done by matching the sequences to the EMBL nucleotide sequence database (Kanz et al. 2005). MOTUs without hits were retained as unclassified. Demultiplexed fastq files were added to the Zenodo repository (https://doi.org/10.5281 /zenodo.1318593).
We then divided our dataset of EMBL-assigned sequences into smaller datasets. These datasets were created to refine those taxonomic assignments with a phylogenetic-based coalescent model approach (PCMA) after Zimmermann et al. (2015). This method identifies taxonomic boundaries from the variation in branching rates of a phylogenetic tree (Monaghan et al. 2009). We constructed 12 datasets: Achnanthidiaceae, Bacillariaceae, centrics, Cocconeidaceae, Cymbellales, Eunotiaceae, Fragilariophycidae, Mastogloiales, Sellaphoraceae-Pinnulariaceae-Caloneis-Mayamaea-Stauroneidaceae, Naviculales (excluding Sellaphoraceae-Pinnulariaceae-Caloneis-Mayamaea-Stauroneidaceae), Surirellales-Rhopalodiales, and Bacillariophyta for genusunassigned sequences. Each dataset contained the environmental sequences generated by HTS, reference sequences from the barcode reference library we made via clonal culturing, reference sequences from the BGBM Diatom Sequence Reference Database (unpublished), and annotated diatom sequences from the NCBI nucleotide database. We aligned the datasets in the software MEGA version 6.06 (Tamura et al. 2013), which uses the MUSCLE (Edgar 2004) alignment algorithm following Zimmermann et al. (2015). We then visualized and manually improved the alignments in PhyDE version 0.9971 (Müller et al. 2005).
We conducted phylogenetic analyses on those datasets with Maximum Likelihood (ML), as implemented by RAxML version 8 (Stamatakis 2006, Stamatakis et al. 2008 in the CIPRES platform (Miller et al. 2010). We used a model of sequence evolution that was general time reversible, had a gamma distribution (G), and included an estimate of the proportion of invariable sites I (Tavaré 1986). We did 1000 replicates of this model for bootstrap analysis. We then used FigTree version 1.4.2 to draw the phylogenetic trees (Rambaut 2014). To taxonomically assign MOTUs unambiguously at species level, we only considered well-supported clades (≥60% bootstrap support) that were also correlated to morphological data (from environmental samples or clonal cultures). To assign less supported clades (<60% bootstrap support) we followed Zimmermann et al. (2015): sequences in clades clustered together with reference sequences from the International Nucleotide Sequence Database Collaboration (INSDC) databases (i.e., DDBJ, EMBL-EBI and NCBI) but with no morphological match were given the abbreviation cf. (confer) before the epithet, e.g., Anomoeoneis cf. sphaerophora; and sequences in distinct clades were considered as unspecified members of a genus (e.g., Tryblionella sp.). We used a 95% identity threshold to assign distinct genera after our experience setting this threshold at different values because the 95% threshold balanced the actual intrageneric variation of the amplified region and the artefactual variation created during sequencing (Brown et al. 2015).
We visualized how well morphology (LM and SEM of environmental samples and cultures) and DNA metabarcoding were able to identify taxa with Venn diagrams.

Comparison of abundance data
We evaluated whether the morphological and metabarcoding methods gave similar estimates of relative diatom abundance by comparing the abundances of the 5 most common taxa found with each method. To make the results from the 2 methods comparable, we first transformed the valve counts and sequence reads into relative abundances. After observing large disparities between the abundances found by each method, we assessed whether applying Correction Factors (CFs) for biovolume to the metabarcoding data improved abundance estimates compared to abundances obtained by morphology. We used the CFs calculated by V. Vasselon et al. (unpublished), even though there were taxonomic differences, because the diatoms they studied were similar in size and therefore in biovolume. We used the CF for Achnanthidium minutissimum for Achnanthidium sp. 115 (we merged A. sp. 1 and sp. 5 in our analyses because the barcoding marker did not differentiate these 2 distinct morphodemes) and for Achnanthidium cf. tropicocatenatum because these taxa are similar in size. We used the CF for Cocconeis placentula for Cocconeis sp. 2 and the CF for Navicula cryptotenella for Navicula notha, also because of size similarities.

Retrieval of barcode sequences from HTS data
We explored the potential to obtain barcode sequences for individual taxa from our HTS data following some of the guidelines from Rimet et al. (2018) for the rbcL gene. Sequences that are good candidates for barcodes include those that: 1) are among the most abundant in a sample, 2) are phylogenetic neighbors of the same neighbor taxa expected from morphological observations, and 3) have neither indels nor stop codons. The 18S rRNA gene is a non-coding region, unlike rbcL, so we relaxed the criteria requiring no indels to a maximum of 1 indel after aligning the putative barcode sequences with sequences from closely related taxa. We followed the same protocols for sequence alignment, phylogenetic tree construction, and tree drawing that we followed for the PCMA.

Barcode reference library of the Lerma-Chapala River Basin
We isolated 111 cultures (strains) from the 18 samples we used to compare morphology and metabarcoding in this study. These 111 cultures corresponded to 46 infrageneric taxa in 21 genera. In combination with the strains cultured from all 42 samples taken by Mora et al. (2017), we cultivated a total of 190 strains that belong to 72 infrageneric taxa in 24 genera (Table S1). Of the 190 strains we produced, 100 strains have sequences that are novel for the 18S V4 rRNA gene, because no identical sequences were available in the International Nucleotide Sequence Database Collaboration (INSDC) databases (DDBJ, EMBL-EBI and NCBI). Of the 24 total genera we identified, 5 genera were underrepresented in the INSDC databases because they had <10 entries for the 18S locus. Of those 5 underrepresented genera, the 9 sequences generated from the 9 strains of

Morphological analysis of environmental samples
We found a total of 204 infrageneric taxa (species and varieties) in 42 genera from the data analyzed by Mora et al. (2017) with further observations made in the present study that improved taxon detection and identification (Table 2, Fig. 2A). From that total, we found 148 infrageneric taxa in 38 genera in the LM count data to determine relative abundances. The additional 56 infrageneric taxa were found by scanning the whole slides under LM to look for rare taxa, and also by SEM examinations (Table 2, Fig. 2A). Among the additional 56 infrageneric taxa, we detected 4 genera not recorded in the 38 genera from the LM counts. Those 4 genera were Achnanthes Bory and Neidium Pfitzer, which were only observed under LM, and Cymbella Agardh and Iconella Jurilj, which were found with both LM and SEM (Table 2). Of the counted taxa, the 5 most abundant across all samples were, in decreasing order, Achnanthidium sp. 5, Achnanthidium sp. 1, Gomphonema parvulum (Kützing) Kützing, Cocconeis sp. 2, and Achnanthidium cf. tropicocatenatum Marquardt, C. E. Wetzel et Ector (Table S2).    Morphological diversity detected by clonal culturing We identified 46 infrageneric taxa in the clonal cultures, but only 34 of them were also found in the microscopy examinations of environmental samples, which means that 12 taxa were found only after culturing. This result increased the total richness to 216 infrageneric taxa in 43 genera ( Fig. 2A, Table 2) from the 204 infrageneric taxa in 42 genera identified only from the environmental samples. Simonsenia strains were neither observed by LM nor by SEM from environmental samples, adding 1 genus to the 42 found in environmental samples.

DNA metabarcoding by HTS
The Illumina MiSeq sequencing run generated 2,738,628 reads from our full data set. After we deleted singleton and chimeric reads, we retained 1,156,360 quality reads. Diatoms comprised 3.8% (43,703) of the quality reads, according to the BLASTn (Altschul et al. 1990) of the EMBL nucleotide database. We obtained a total of 2181 MOTUs from 43,703 diatom reads.
The phylogenetic-based coalescence model approach (PCMA) resulted in 350 taxonomic units (hereafter referred to as infrageneric taxa). To further remove potential sequencing noise, we removed infrageneric taxa that were comprised of only 1 doubleton or 1 tripleton, and had no association with morphology (no valves observed) or a reference sequence (no reference sequences in the corresponding tree branch). These removals reduced the number of taxa to 331 in 35 genera. We additionally removed 65 of these infrageneric taxa because we could not assign them to a described genus based on the 95% identity threshold, which resulted in a richness of 266 infrageneric taxa in 35 genera. From the remaining 266 infrageneric taxa we were able to confidently assign 94 of them unambiguously because a morphological correlation was possible, or a correlation to a sequence from our own barcode reference database or to sequences from NCBI in well-supported clades (≥ 60 bootstrap support).

Comparison of diatom composition from morphology and metabarcoding
The morphological analysis recovered 216 taxa, in contrast to the 266 we recovered with metabarcoding. We found 43 genera based on morphological identification and 35 based on metabarcoding. In total, we found 14 genera only by morphological identification, 29 with both morphological identification and metabarcoding, and 6 only by metabarcoding (Fig. 2B, Table 2). The combination of the total morphological richness of 216 taxa with the 94 infrageneric taxa unambiguously assigned to species level by metabarcoding resulted in a total of 250 infrageneric taxa (Fig. 2C). We detected 60 of these taxa with both methods and 34 only by metabarcoding. The barcode reference library we present here allowed the assignment of 23 infrageneric taxa by infrageneric taxa identified after additional scanning of the slides after the 500-valve counts were completed under LM, as well as taxa observed during SEM examinations. R 5 infrageneric taxa identified from clonal cultures isolated from the 18 samples that were the focus of this study. B.-Genera identified by morphology (M) and metabarcoding (HTS). C.-Infrageneric taxa identified by morphology and metabarcoding (only assigned taxa to species level from metabarcoding shown). M 5 taxa identified by morphology, HTS 5 taxa identified by metabarcoding, and R 5 taxa retrieved by metabarcoding whose taxonomic assignment was only possible with the reference library of barcode sequences established in this study. For a color version of this figure see Fig. S1.

Volume 38
September 2019 | 455 metabarcoding. Without our reference library, those 23 infrageneric taxa would have been left unassigned because no matching reference sequences were available in the INSDC databases before our study (Fig. 2C).

HTS data as a source of barcodes
After in-depth examination of HTS data obtained in our study, we proposed barcode sequences of the V4 region (18S) for 2 taxa, Iconella delicatissima (F. W. Lewis) Ruck et Nakov (Fig. 5A, B) and Navicula notha Wallace (Fig. 5C-G).
Iconella delicatissima (F. W. Lewis) Ruck et Nakov (Fig. 5 A, B) ; Stenopterobia delicatissima (Lewis) Van Heurck The genus Iconella has been recently resurrected to accommodate Stenopterobia and the 'robustoid' members of Surirella and Campylodiscus (Ruck et al. 2016a, b, Jahn et al. 2017a. We found Iconella sequences only in sample 18. The morphological examination of this sample confirmed the presence of Iconella delicatissima. Surirella angusta was the only other member of the Surirellales found in this sample, but this species belongs to Surirella sensu stricto (Ruck et al. 2016b). A total of 52 sequence reads were obtained for Iconella delicatissima from the 6873 total reads obtained from sample 18 (Laguna de Servín 2), making it the 13 th most abundant taxon in this sample. After BLASTing this sequence on NCBI, the 3 most similar sequences belonged to the genus Stenopterobia, with 98% similarity to the 2 most similar sequences (Stenopterobia pumila and S. curvula) and 95% similarity to the 3 rd most similar sequence (Stenopterobia sp. 50). Other sequences with 95% similarity include 1 sequence of Surirella sp. and 2 sequences of Campylodiscus levanderi. After we aligned our sequence with 17 sequences from the resurrected genus Iconella (Ruck et al. 2016b), including those of the aforementioned taxa, our sequence had only 1 indel in comparison with the 3 available sequences of Stenopterobia and no indels compared with other included surirelloid species. In the unrooted ML phylogenetic tree, Iconella delicatissima clusters with Stenopterobia pumila, S. curvula, and S. sp. 50 in a well-supported clade (bootstrap value 5 96%) (Fig. S2). This new barcode sequence has been submitted to the ENA under accession number LS990839.

C-G)
We found sequences of N. notha in 17 out of the 18 samples we analyzed with HTS. The read abundance of this taxon across all samples was the 2 nd highest after G. parvulum s.l. We also found N. notha in our morphological examinations of these samples. We used samples 4, 5, and 14 as the barcode sources, because N. notha was the only Navicula representative found from valve counts in these samples. Further, N. notha reached high relative valve abundances in these samples (2.9% in sample 4; 4.5% in sample 5; 5.4% in sample 14). The relative abundance of sequence reads in those samples was as high as 16% (sample 4), 34% (sample 5), and 22% (sample 14). The sequences from these 3 sites were identical, so we only submitted 1 sequence to the ENA (from site 5, accession number LS990785). The closest sequence that corresponded to N. notha on the NCBI database was Navicula cryptotenelloides, which had a 4 bp difference (99% similarity). Other sequences with high similarity (98%) to our sequence included Navicula cryptotenella and N. reinhardtii. However, we did not observe N. cryptotenelloides or N. reinhardtii in our samples, so we are confident that our sequence does not correspond to these taxa. Further, we are confident that our sequence does not correspond to N. cryptotenella even though we observed this taxon by microscopy because it was only present at 1 site, it was not recorded during the 500-valve counts to calculate relative abundances, and we only detected it after additional scans of the slides were performed. We aligned our N. notha sequence to 29 Navicula sequences present in the NCBI nucleotide database. This alignment had no indels. Navicula notha clustered with N. cryptotenelloides, N. cryptocephala, N. cryptotenella, and N. reinhardtii (Fig. S3) in a well-supported clade (bootstrap value 5 84%) in the unrooted ML phylogenetic tree.

DISCUSSION
Barcode reference library of the Lerma-Chapala River Basin The regional vouchered barcode reference library presented here is the 1 st of its kind for stream diatoms from Mexico. Only 13 strains of epicontinental diatoms from Mexico had previously been cultured, sequenced, and published (Zimmermann et al. 2011, Jahn et al. 2017b). There are other entries in the INSDC databases that refer to diatoms from epicontinental locations in Mexico, but all of them correspond to uncultured organisms with no vouchered material. These facts highlight the importance of having vouchered material to allow traceabil- ity of the data and ensure the availability of reserve material for both morphological and molecular studies for other diatom studies (e.g., phylogenetic and biogeographic). From the 190 strains established here, 87 strains corresponded to already described species, 49 strains were closely related (cf. confer) to already described species, 7 were similar (aff. affinis) to other species, and 47 strains were only named at the genus level even after a thorough morphological and bibliographical examination. Some of these unnamed, similar, or closely related strains should probably be described as new species. Our finding of such a large fraction (∼50%) of unidentified taxa is not surprising because our samples were collected in Central Mexico, a region not extensively studied for diatoms and for which no monographs have been published (Mora Hernández 2018). A similar proportion of unidentified taxa was found for a small sample (n 5 26 strains) of polar diatoms (58%) (Stachura-Suchoples et al. 2016). Even in thoroughly investigated regions like Berlin, 10% of the diatom species identified by Zimmermann et al. (2014) were newly described.

Diatom composition detected by microscopy and metabarcoding
The 500-valve counts based on LM from Mora et al. (2017) led to the identification of 148 taxa. Fifty-six further taxa were found by these authors and us after additional scanning slides under LM and during SEM examinations, raising the richness to 204 infrageneric taxa in 42 genera. This increase in species with increasing sampling effort (in this case, increasing the number of slides observed under microscopy) is a well-known property (Von Falkenhayn 2008, Gotelli and Colwell 2011) of samples such as ours that are characterized by high diversity and a high proportion of rare taxa.
We found a higher taxonomic richness after the PCMA of our metabarcoding data than we did from morphological analyses: 266 from metabarcoding versus 216 taxa identified morphologically (including 12 taxa only identified from cultures). When we took only the unambiguously assigned infrageneric taxa into account (94), ⅔ corresponded to taxa found by microscopy in this study. The assignment of ¼ of those taxa (23) was only possible with our barcode reference sequences. This number represents a Figure 5. Taxa for which barcodes were retrieved from HTS data. A and B 5 Iconella delicatissima from Laguna de Servín 2 (sample 18) observed under light microscopy. C-G 5 Navicula notha. C-E were observed under light microscopy and F and G were observed with scanning electron microscopy. All specimens were collected from La Laborcilla 1 (sample 5). The black scale bar 5 10 lm. fraction of the diatom diversity of the region, but the barcode reference library we present here represents a milestone in documenting diatom taxonomic composition in streams of Mexico. In our study alone, the reference library we created allowed us to confidently assign 23 distinct taxa that would have otherwise been left unassigned at the infrageneric level.
Metabarcoding could lead to overestimation of diversity because it normally retrieves more species or MOTUs than morphologically identified taxa (see below; Zimmermann et al. 2015, Vasselon et al. 2017b). However, metabarcoding is a good way to screen for biodiversity because it can show gaps in the described groups of organisms and in barcode reference libraries. Thus, metabarcoding can lead to the refinement of diversity assessments. The main reason cited in the literature for incongruence of taxa lists obtained from morphology and metabarcoding analyses is the incompleteness and lack of accuracy of reference databases (e.g., Zimmermann et al. 2014, 2015, Apothéloz-Perret-Gentil et al. 2017, Vasselon et al. 2017b, which hinders correct taxonomic assignment of environmental sequences. Missing taxa in reference databases would not normally be identified in environmental sequences, whereas sequences with the wrong taxonomy in databases will generate inaccurate identifications (Kermarrec et al. 2014, Lejzerowicz et al. 2015. The proportion of sequence reads of diatoms in this study (3.8%) was low compared with the total number of high-quality reads we obtained (43,703 out of 1,156,360). This could be the result of an actual low abundance of diatoms in the sampled sites or a low proportion of live diatom cells. Alternatively, the relatively low proportion of diatom sequences could be because the primers we used in this study had a higher affinity for other protists during PCR, such as Chrysophytes, which comprised 70% of the total reads. This result is unsurprising because Chrysophytes are an abundant and diverse group in microbial freshwater food webs (del Campo andMassana 2011, Grossmann et al. 2016).

Concealed diversity revealed by clonal culturing
We found 12 infrageneric taxa in 7 genera in the cultures that we did not find by microscopy in our environmental samples. For example, we successfully cultivated clones from the genus Simonsenia, but we did not detect this genus in the environmental samples analyzed by LM and SEM. These results give rise to the question: how much of the diversity in a sample remains concealed even after exhaustive microscopy examination? In our study, 6% of the diversity (12 out of 216 taxa) was only detected through laboratory culturing. This result could have occurred if the culture media and culturing conditions (i.e., light, day/ night cycle and temperature) allowed taxa to grow that were otherwise too rare to be detected through even thorough microscopy examinations. In this case, the suitable cultur-ing conditions would have enabled these taxa to reach abundances high enough to be observed and picked up during the cell isolations. Previous studies of cyanobacteria (Ferris et al. 1996) and marine bacteria (Connon and Giovannoni 2002) have reported that culture media at very low concentrations (e.g., 3 orders of magnitude lower than commonly used) can lead to the identification of taxa undetected by morphological examinations and by standard culturing techniques. Here, some strains were cultured in media at concentrations ½ and ¼ of the manufacturer recommendations, in addition to the recommended concentration. These lower concentrations could partially explain why this concealed diversity was uncovered.

Richness overestimation
The taxa list retrieved by metabarcoding (266 genusassigned plus 65 unassigned to genus) was larger than the 216 morphology-based list even after a thorough morphological examination of samples. Generally, the number of species or MOTUs that are generated by DNA metabarcoding deviates from the number of taxa observed morphologically (Cowart et al. 2015, Groendahl et al. 2017), normally recovering more taxa than morphology-based approaches. Several biological, environmental, and technical factors contribute to this.
The most important biological factor that could cause DNA-based approaches to overestimate richness is the natural intraspecific and intragenomic variability of the barcoding marker. This variability is particularly problematic when a single traditionally recognized species or bioindicator taxon has multiple genotypes at the barcoding region (Brown et al. 2015, Bálint et al. 2016. When a single taxon has multiple genotypes at the barcoding region, members of that taxon may cluster into different MOTUs, artificially inflating taxonomic richness. High intraspecific genetic variation is common in nearly all bioindicator groups, such as aquatic insects (Alp et al. 2012, Elbrecht et al. 2014) and diatoms (Rynearson andArmbrust 2000, Trobajo et al. 2009). Moreover, some taxa show high intragenomic polymorphism, such as nematodes (Bik et al. 2013), foraminifera (Weber and Pawlowski 2014), and prokaryotes (Sun et al. 2013). Intragenomic polymorphism has not been widely assessed in diatoms (Alverson and Kolnick 2005) but it could also contribute to MOTU inflation.
MOTU richness can also be artificially inflated through technical errors at different steps of sample processing, especially during amplification (Fonseca et al. 2012, Kermarrec et al. 2013b, Bálint et al. 2016, Elbrecht et al. 2017b) and sequencing (Meacham et al. 2011, Schirmer et al. 2015. In contrast to amplification and sequencing, DNA extraction methods do not affect MOTUs richness significantly in diatoms (Vasselon et al. 2017a).
The MOTU delimitation approach is another factor that influences richness estimation and interpretation. However, MOTUs do not necessarily correspond to species, and can fail to identify meaningful ecological or phylogenetic units in a straightforward manner (Ryberg 2015, Bálint et al. 2016. To ameliorate this issue, MOTUs have been further analyzed by phylogenetic-based approaches that assign MOTUs to specific taxa (Monaghan et al. 2009, Zimmermann et al. 2015. Morphology-based assessments can also lead to richness overestimation because diatoms are oxidized before microscopy slides are prepared. Thus, taxonomic identifications, counts to calculate abundance, and diversity estimations rely on the valves of dead cells that could be transported from locations other than the target assemblage (Sawai 2001, Potapova andCharles 2005). The proportion of live diatoms found in lotic environments varies greatly, ranging from 2 to 98% (Gillett et al. 2011). Despite this large variation, the oxidation method for slide preparation gives taxonomic confidence because it allows the visualization of the fine structure of the cell walls, which is needed for identification (Gillet et al. 2009).

Differences in abundance data
Discrepancies in abundance estimates based on valve counts and sequence reads are widely debated by the metabarcoding community, particularly in terms of the use of the number of sequence reads to infer taxa abundances and its application to biomonitoring Leese 2015, Vasselon et al. 2018). The barcoding marker and its ability to discriminate among closely related species, as well as primer specificity, are very important and can hinder the otherwise straightforward use of sequence reads to determine species abundances (Elbrecht et al. 2017a). These issues were evident in our results, especially for Gomphonema and Nitzschia, because the 18S V4 barcode region does not have the discriminatory power to differentiate among some closely related species within these genera. The existence of semi-cryptic diversity within G. parvulum has been documented (Kermarrec et al. 2013a), and only a multi-marker phylogeny coupled with detailed micromorphology examinations are able to disentangle some of the species within this species complex that was once thought to have a cosmopolitan distribution . For N. palea, both morphology and metabarcoding found 3 taxa, but we treated this species as a complex because it was not possible to dissociate the nominate variety from the varieties N. palea var. debilis and N. palea var. tenuirostris we obtained from the morphological analysis from the 3 taxa we found through metabarcoding. Morphological, genetic, and mating studies of N. palea concluded that this taxon is comprised of 3 or more species but molecular and mating experiments do not separate these taxa into the varieties traditionally recognized morphologically (Trobajo et al. 2009(Trobajo et al. , 2010. The disparities observed in abundances with each method (Fig. 3) for N. palea and G. parvulum might have conse-quences for their use as bioindicators because it was impossible to match taxa assigned by metabarcoding with species recognized morphologically. This limitation might be overcome with taxonomy-free approaches to bioindication that avoid any taxonomic assignments to morphologically recognized taxa (Apothéloz-Perret-Gentil et al. 2017, Tapolczai et al. 2019); this is done by a direct calibration and calculation of ecological values of the MOTUs, without prior species assignment.
Cell size may be another factor in the abundance disparity we observed, as shown in our results for Achnanthidium and Ulnariataxa that differ markedly in cell size. On one hand, the Achnanthidium species compared here have a size of 7 to 23 lm and represented 28.5% of the total diatom abundance in the morphological analyses, but comprised only 1.4% of the total sequence read counts. On the other hand, the cells of Ulnaria cf. ulna analyzed here ranged in size from 70 to 200 lm and represented only 1.1% of the total valve abundance, but consisted of 6.6% of the total read counts. There may be an association between cell size, biovolume, and gene copies of the SSU rDNA , Godhe et al. 2008) that could partially explain the disparities in abundances between these methods. Thus, the contrast in cell size, and therefore in biovolume, could explain why Achnanthidium spp. were underrepresented in the read counts relative to valve abundance, compared with the relatively large Ulnaria cf. ulna, which was overrepresented in read abundances relative to valve counts.
Applying CFs based on the cell biovolume of diatoms can reduce the disparity between morphological and molecular abundance data by as much as 45% in both mock communities and environmental samples . However, when we applied CFs for the 18S gene (V. Vasselon et al. unpublished) we saw no overall improvement in the disparity among the abundances recorded by morphology and metabarcoding. One of the reasons for this lack of overall improvement may be the 7 CFs we applied, because only 2 corresponded to the actual species for which they were proposed. The other 5 CFs we used are from taxa that belong to the same genus and are similar in size. Our results point out that species-specific CFs may be needed. The technical errors generated during amplification and sequencing previously described may also affect the disparities observed from both methods.

HTS data as a source of barcodes
Several challenges and limitations make it difficult to establish complete barcode libraries for diatoms, which currently rely mostly on clonal cultures. A primary challenge is the time-consuming process of single cell isolation and culture maintenance. Furthermore, culturing can often be unsuccessful because of recalcitrance of species to culturing conditions Chepurnov 2004, Rimet et al. 2018). HTS data is an alternative source of barcodes ). As we demonstrate here with I. delicatissima and N. notha, HTS data can be used as a source of barcodes if data are analyzed carefully.
Retrieving I. delicatissima was relatively straightforward because it was the only representative of the genus Iconella we found. In this sample, another representative of the Surirellales, S. angusta, was observed which might have hindered our findings, but we had reference sequences for S. angusta obtained via culturing and Sanger sequencing, so we were able to easily distinguish these species.
Retrieving N. notha was more difficult because of the high abundance of Navicula taxa in our samples, as evidenced not only by sequences derived from HTS but also from microscopy observations. Thus, we obtained sequences of this taxon from 3 samples in which N. notha was the only representative of the genus Navicula as detected by microscopy observations.

Conclusion
Our study demonstrates that the combination of morphological (LM and SEM) and molecular (metabarcoding via HTS) methods applied to environmental samples, in combination with a regional barcode reference library based on clonal culturing, increases the detection and identification of diatom species. Thus, our work highlights the complementary aspects of classical taxonomy and DNA metabarcoding (i.e., the importance of their reciprocal illumination). Even with major advances in the development and standardization of molecular tools for diversity assessments and monitoring , the role of morphology in species detection and identification remains central in the ecogenomic era.