Raphaël Morard
email: rmorard@marum.de.

Kate F Darling

Frédéric Mahé

Stéphane Audic

Yurika Ujiié

Agnes K F Weiner

Aurore André

Heidi A Seears

Chris M Wade

Frédéric Quillévéré

Agnes K M Weiner

Christophe J Douady

Gilles Escarguel

Thibault De Garidel-Thoron

Michael Siccha

Michal Kucera

Colomban De Vargas

PFR 2 : a curated database of planktonic Foraminifera18S ribosomal DNA as a resource for studies of plankton ecology, biogeography, and evolution

Keywords: Planktonic foraminifera, 18S ribosomal DNA, molecular ecology, genetic diversity, molecular taxonomy, sequence database

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

Introduction

Despite their ubiquity and the critical role they play in global biogeochemical cycles, unicellular eukaryotes (protists) remain the most poorly known domain of life (e. g., [START_REF] Pawlowski | CBOL protist working group: barcoding eukaryotic richness beyond the animal, plant, and fungal kingdoms[END_REF]. Because of their extreme morphological and behavioral diversity, the study of even relatively narrow lineages requires a high degree of taxonomic expertise (e. g., Guillou et al., 2012, Pawlowski and[START_REF] Pawlowski | A plea for DNA barcoding of foraminifera[END_REF]. As a result, the knowledge of protistan ecology and evolution is limited by the small number of taxonomists, resulting in scarcity of taxonomically well-resolved ecological data. As an alternative approach, numerous studies have demonstrated the potential of identification of protists by means of short DNA sequences or barcodes (e. g., [START_REF] Saunders | Applying DNA barcoding to red macroalgae: a preliminary appraisal holds promise for future applications[END_REF][START_REF] Sherwood | Universal primers amplify a 23S rDNA plastid marker in eukaryotic algae and cyanobacteria[END_REF][START_REF] Hollingsworth | A DNA barcode for land plants[END_REF]Nossonova et al., 2010;[START_REF] Pawlowski | Short rDNA barcodes for species identification in foraminifera[END_REF][START_REF] Hamsher | Barcoding diatoms: exploring alternatives to COI-5P[END_REF][START_REF] Stern | Environmental barcoding reveals massive dinoflagellate diversity in marine environments[END_REF][START_REF] Schoch | Nuclear ribosomal internal transcribed spacer (ITS) region as a universal DNA barcode marker for Fungi[END_REF], both at the single-cell and metacommunity levels (e. g., [START_REF] Sogin | Microbial diversity in the deep sea and the underexplored "rare biosphere[END_REF][START_REF] Logares | Patterns of rare and abundant marine microbial eukaryotes[END_REF]. Such barcoding/metabarcoding approaches critically rely on the fidelity of the marker gene with respect to specificity (avoiding ambiguity in identification), comprehensiveness (assuring all taxa in the studied group are represented in the reference barcode database) and accuracy (assuring that barcode assignments are consistent with a coherent, phenotypic taxonomic framework; e. g., [START_REF] Zimmermann | Taxonomic reference libraries for environmental barcoding: a best practice example from diatom research[END_REF]). These three pre-requisites are rarely found in protists, where classical morphological taxonomy is often challenging, DNA extraction and sequencing from a single cell is prone to contamination, and a large portion of the diversity in many groups remains unknown (e. g., [START_REF] Mora | How many species are there on Earth and in the ocean?[END_REF]. In this respect, planktonic foraminifera represent a rare exception.

Planktonic foraminifera are ubiquitous pelagic marine protists with reticulated pseudopods, clustering within the Rhizaria [START_REF] Nikolaev | The twilight of Heliozoa and rise of Rhizaria, an emerging supergroup of amoeboid eukaryotes[END_REF]. The group is marked by a rather low number of extant morphospecies (47; [START_REF] Hemleben | Modern Planktonic Foraminifera[END_REF], which can be distinguished using structural characteristics of their calcite shells. Their global geographic distribution, seasonal dynamics, vertical habitats and trophic behavior have been thoroughly documented by analyses of plankton hauls (e.g., [START_REF] Bé | Ecology of planktonic foraminifera and biogeographic patterns of life and fossil assemblages in the Indian Ocean[END_REF], sediment trap series (e.g., Zaric et al., 2005) and thousands of surface sediment samples across the world oceans (e.g., [START_REF] Kucera | Reconstruction of sea-surface temperatures from assemblages of planktonic foraminifera: multi-technique approach based on geographically constrained calibration data sets and its application to glacial Atlantic and Pacific Oceans[END_REF]. Their outstanding preservation in marine sediments resulted in arguably the most complete fossil record, allowing comprehensive reconstruction of the evolutionary history of the group [START_REF] Aze | A phylogeny of Cenozoic macroperforate planktonic foraminifera from fossil data[END_REF]. Over the last two decades, the morpho-taxonomy and phylogeny of the group have been largely confirmed by molecular genetic analyses (e.g., Aurahs et al., 2009a) based on the highly informative, ~1,000 bp fragment at the 3'end of the 18S rDNA gene.

These analyses confirmed that the morphological characters used to differentiate planktonic foraminifera taxa are phylogenetically valid both at the level of morphological species and at the level of higher taxa. The studied gene fragment contains six hypervariable expansion segments, some unique to foraminifera, providing excellent taxonomic resolution [START_REF] Pawlowski | Short rDNA barcodes for species identification in foraminifera[END_REF]. Analyses of this fragment revealed the existence of genetically distinct lineages within most of the morphospecies, which likely represent reproductively isolated units [START_REF] Darling | Molecular Phylogeny of the planktic foraminifera[END_REF][START_REF] Darling | Planktic foraminiferal molecular evolution and their polyphyletic origins from benthic taxa[END_REF][START_REF] Darling | The Diversity and Distribution of Modern Planktic Foraminiferal Small Subunit Ribosomal RNA Genotypes and their Potential as Tracers of Present and Past Ocean Circulations[END_REF][START_REF] Darling | Molecular evidence for genetic mixing of Arctic and Antarctic subpolar populations of planktonic foraminifers[END_REF][START_REF] Darling | Seasonal distribution of genetic types of planktonic foraminifer morphospecies in the Santa Barbara Channel and its paleoceanographic implications[END_REF][START_REF] Darling | Molecular evidence links cryptic diversification in polar planktonic protists to Quaternary climate dynamics[END_REF][START_REF] Darling | A resolution for the coiling direction paradox in Neogloboquadrina pachyderma[END_REF][START_REF] Darling | Global molecular phylogeography reveals persistent Arctic circumpolar isolation in a marine planktonic protist[END_REF][START_REF] Darling | Surviving mass extinction by bridging the benthic/planktic divide[END_REF][START_REF] Darling | The genetic diversity of planktic foraminifera and the global distribution of ribosomal RNA genotypes[END_REF][START_REF] Wade | Early Evolutionary Origin of the Planktic Foraminifera Inferred from Small Subunit rDNA Sequence Comparisons[END_REF][START_REF] De Vargas | Phylogeny and rates of molecular evolution of planktonic foraminifera: SSU rDNA sequences compared to the fossil record[END_REF][START_REF] Poseidon | Poseidon[END_REF], 2001, 2002[START_REF] De Vargas | Molecular versus taxonomic rates of evolution in planktonic foraminifera[END_REF][START_REF] Stewart | Genotypic variability in subarctic Atlantic planktic foraminifera[END_REF]Aurahs et al., 2009b[START_REF] Aurahs | A revised taxonomic and phylogenetic concept for the planktonic foraminifer species Globigerinoides ruber based on molecular and morphometric evidence[END_REF][START_REF] Ujiié | Cryptic diversity in planktonic foraminifera in the northwest Pacific Ocean[END_REF][START_REF] Ujiié | Molecular evidence for an independent origin of modern triserial planktonic foraminifera from benthic ancestors[END_REF][START_REF] Ujiié | Longitudinal differentiation among pelagic populations in a planktic foraminifer[END_REF][START_REF] Morard | Morphological recognition of cryptic species in the planktonic foraminifer Orbulina universa[END_REF][START_REF] Morard | Worldwide genotyping in the planktonic foraminifer Globoconella inflata: implications for life history and paleoceanography[END_REF][START_REF] Morard | Ecological modeling of the temperature dependence of cryptic species of planktonic foraminifera in the Southern Hemisphere[END_REF][START_REF] Seears | Ecological partitioning and diversity in tropical planktonic foraminifera[END_REF][START_REF] Quillévéré | Global scale same-specimen morpho-genetic analysis of Truncorotalia truncatulinoides: A perspective on the morphological species concept in planktonic foraminifera[END_REF][START_REF] Weiner | Vertical niche partitioning between cryptic sibling species of a cosmopolitan marine planktonic protist[END_REF][START_REF] Weiner | Phylogeography of the tropical planktonic foraminifera lineage Globigerinella reveals isolation inconsistent with passive dispersal by ocean currents[END_REF][START_REF] André | SSU rDNA Divergence in Planktonic Foraminifera: Molecular Taxonomy and Biogeographic Implications (V Ketmaier[END_REF]. In order to assess the ecology and biogeography of such cryptic species, large numbers of rDNA sequences from single-cell extractions collected across the world oceans have been generated for most morphospecies (Figure 1). Due to this extensive single-cell rDNA sequencing, the genetic and morphological diversity of planktonic foraminifera have been linked together to a degree that now allows for transfer of taxonomic expertise. The knowledge of the genetic and morphological taxonomy of the group allows the establishment of an exceptionally comprehensive reference genetic database that can be further used to interpret complex data from plankton metagenomic studies with a high level of taxonomic resolution. Because planktonic foraminifera are subject to the same ecological forcing as other microplankton, including the dominance of passive transport in a relatively unstructured environment, huge population sizes, and basin-scale distribution of species, they can potentially serve as a model for the study of global ecological patterns in other groups of pelagic protists, whose diversity remains largely undiscovered [START_REF] Mora | How many species are there on Earth and in the ocean?[END_REF].

By early 2014, 1,787 partial 18S rDNA sequences from single-cell extractions of planktonic foraminifera were available in public databases. However, their NCBI taxonomy is often inconsistent, lacking standardization. It includes (and retains) obvious identification errors, as discussed by Aurahs et al. (2009a) and [START_REF] André | SSU rDNA Divergence in Planktonic Foraminifera: Molecular Taxonomy and Biogeographic Implications (V Ketmaier[END_REF], and their annotation lacks critical metadata. In addition, an equivalent number of rDNA sequences not deposited in public databases have been generated by the co-authors of the present study. Collectively, the existing rDNA sequences from single cells collected throughout the world oceans cover the entire geographic and taxonomic range of planktonic foraminifera. This collection unites the current morphological, genetic, ecological, and biogeographic knowledge of the group and may serve as a Rosetta Stone/Philae Obelisk for interpreting metabarcoding data (Pawlowski et al., 2014). To pave the way for future exploitation of this resource, we combined all published and unpublished planktonic foraminifera rDNA sequence data and curated the resulting database with a semiautomated bioinformatics pipeline. The resulting Planktonic Foraminifera Ribosomal Reference database (PFR 2 ) is a highly resolved, fully annotated and internally entirely consistent collection of 18S rDNA sequences of planktonic foraminifera, aligned and evaluated in a way that facilitates, among others, direct assessment of barcoding markers.

Material and Methods

Primary database assembly

A total of 1,787 18S rDNA sequences of planktonic foraminifera were downloaded from the GenBank query portal (http://www.ncbi.nlm.nih.gov/; release 201) on the 14 th of May 2014. The taxonomic path and metadata for these sequences were extracted from NCBI and supplemented by information in original papers when available. The metadata associated to each sequence consisted of: (i) their organismal origin (specimen voucher, taxonomic path, infra specific genetic type assignment), (ii) their methodological origin (direct sequencing or cloning), and (iii) their spatio-temporal origin (geographic coordinates, depth, and time of collection). Metadata were described using standard vocabularies and data formats. For 47 sequences, the coordinates of the collection site could not be recovered, in which case the locality was described in words (Supplementary Material 1).

We next compiled all unpublished 18S rDNA sequences generated by the co-authors of this paper and linked them with the same suite of metadata. These sequences originate from singlecell extractions of planktonic foraminifera collected by stratified or non-stratified plankton net hauls, in-situ water pumping, as well as SCUBA diving. After collection, the specimens were individually picked under a stereomicroscope, cleaned, taxonomically identified and transferred into DNA extraction buffer or air-dried on cardboard slides and stored at -20°C or -80°C. DNA extractions were performed following the DOC [START_REF] Holzmann | Preservation of foraminifera for DNA extraction and PCR amplification[END_REF], the GITC* [START_REF] Morard | Morphological recognition of cryptic species in the planktonic foraminifer Orbulina universa[END_REF], or the Urea [START_REF] Weiner | Phylogeography of the tropical planktonic foraminifera lineage Globigerinella reveals isolation inconsistent with passive dispersal by ocean currents[END_REF] protocols. Sequences located at the 3' end of the 18S rDNA were obtained following the methodology described in de [START_REF] Darling | Molecular Phylogeny of the planktic foraminifera[END_REF], 1997[START_REF] De Vargas | Phylogeny and rates of molecular evolution of planktonic foraminifera: SSU rDNA sequences compared to the fossil record[END_REF], Aurahs et al. (2009b), [START_REF] Morard | Worldwide genotyping in the planktonic foraminifer Globoconella inflata: implications for life history and paleoceanography[END_REF] and [START_REF] Weiner | Phylogeography of the tropical planktonic foraminifera lineage Globigerinella reveals isolation inconsistent with passive dispersal by ocean currents[END_REF]. A total of 820 new planktonic foraminiferal sequences were analyzed and annotated for this study. In addition, 925 unpublished sequences analyzed in [START_REF] Darling | Molecular evidence for genetic mixing of Arctic and Antarctic subpolar populations of planktonic foraminifers[END_REF][START_REF] Darling | Seasonal distribution of genetic types of planktonic foraminifer morphospecies in the Santa Barbara Channel and its paleoceanographic implications[END_REF][START_REF] Darling | Molecular evidence links cryptic diversification in polar planktonic protists to Quaternary climate dynamics[END_REF][START_REF] Darling | A resolution for the coiling direction paradox in Neogloboquadrina pachyderma[END_REF][START_REF] Darling | Global molecular phylogeography reveals persistent Arctic circumpolar isolation in a marine planktonic protist[END_REF], [START_REF] Darling | The genetic diversity of planktic foraminifera and the global distribution of ribosomal RNA genotypes[END_REF], [START_REF] Seears | Ecological partitioning and diversity in tropical planktonic foraminifera[END_REF][START_REF] Weiner | Phylogeography of the tropical planktonic foraminifera lineage Globigerinella reveals isolation inconsistent with passive dispersal by ocean currents[END_REF] were also included. All unpublished sequences, except 177 sequences shorter than 200 bp, were deposited in GenBank under the accession numbers KM19301 to KM194582. Overall, PFR 2 contains data from 460 sites sampled during 54 oceanographic cruises and 15 near shore collection campaigns between 1993 and 2013. It covers all oceanic basins, all seasons, and water depths ranging between the surface and 700 meters (Figure 1; Supplementary Material 1).

Taxonomy

Morphological taxonomy

As the first step in the curation process, the primary taxonomic annotations of all 3,532 18S rDNA sequences gathered from NCBI and our internal databases were harmonized. The identification of planktonic foraminifera is challenging especially for juvenile individuals, which often lack diagnostic characters [START_REF] Brummer | Planktonic foraminiferal ontogeny and new perspectives for micropalaeontology[END_REF]. Thus, many of the published and unpublished 18S rDNA sequences were mislabeled or left in open nomenclature. In some cases the same taxon has been recorded under different names, reflecting inconsistent use of generic names, synonyms and misspelling. To harmonize the taxonomy, we first carried out a manual curation of the original annotations to remove the most obvious taxonomic conflicts in the primary database. To this end, the sequence annotations were aligned with a catalog of 47 species names based on the taxonomy used in [START_REF] Hemleben | Modern Planktonic Foraminifera[END_REF], but adding Globigerinoides elongatus following [START_REF] Aurahs | A revised taxonomic and phylogenetic concept for the planktonic foraminifer species Globigerinoides ruber based on molecular and morphometric evidence[END_REF] and treating Neogloboquadrina incompta following [START_REF] Darling | A resolution for the coiling direction paradox in Neogloboquadrina pachyderma[END_REF]. Thus, the 109 sequences labelled as Globigerinoides ruber (pink) and the 63 labelled as Globigerinoides ruber (white) were renamed as Globigerinoides ruber. The 113 sequences of Globigerinoides ruber and Globigerinoides ruber (white) attributed to the genotype II were renamed Globigerinoides elongatus following [START_REF] Aurahs | A revised taxonomic and phylogenetic concept for the planktonic foraminifer species Globigerinoides ruber based on molecular and morphometric evidence[END_REF]. The 12 sequences labelled Globigerinella aequilateralis were renamed Globigerinella siphonifera following [START_REF] Hemleben | Modern Planktonic Foraminifera[END_REF]. The 7 sequences corresponding to the right-coiled morphotype of Neogloboquadrina pachyderma were renamed Neogloboquadrina incompta following [START_REF] Darling | A resolution for the coiling direction paradox in Neogloboquadrina pachyderma[END_REF]. All taxonomic reassignments were checked by sequence similarity analyses to the members of the new group. Next, we attempted to resolve the attribution of sequences with unresolved taxonomy and searched manually for obviously misattributed sequences. This refers to sequences that are highly divergent from other members of their group but identical to sequences of other well-resolved taxa. Overall, these first steps of manual curation led to the taxonomic reassignment of 124 sequences. All corrections and their justification are documented in the Supplementary Material 1.

Annotation of genetic types

In order to preserve the information on the attribution of 18S rDNA sequences to genetic types (potential cryptic species), we harmonized the existing attributions at this level for species where extensive surveys have been carried out and published. A total of 1,356 sequences downloaded from NCBI were associated with a genetic type label, which was always retained. In addition, 19 sequences labelled as Globigerinoides ruber, 15 as Globigerinoides sacculifer, 36 as Globigerinita glutinata, 6 as Globigerinita uvula, 9 as Globorotalia inflata, 10 as Neogloboquadrina incompta, 6 as Neogloboquadrina pachyderma, 5 as Orbulina universa, 5 as Pulleniatina obliquiloculata, 30 as Hastigerina pelagica, and 32 as Globigerinella siphonifera have been analyzed after their first release in the public domain by Aurahs et al. (2009), [START_REF] Ujiié | Longitudinal differentiation among pelagic populations in a planktic foraminifer[END_REF], [START_REF] Weiner | Vertical niche partitioning between cryptic sibling species of a cosmopolitan marine planktonic protist[END_REF][START_REF] Weiner | Phylogeography of the tropical planktonic foraminifera lineage Globigerinella reveals isolation inconsistent with passive dispersal by ocean currents[END_REF], and [START_REF] André | The cryptic and the apparent reversed : lack of genetic differentiation within the morphologically diverse plexus of the planktonic foraminifer Globigerinoides sacculifer[END_REF][START_REF] André | SSU rDNA Divergence in Planktonic Foraminifera: Molecular Taxonomy and Biogeographic Implications (V Ketmaier[END_REF], and were attributed to a genetic type by these authors. These attributions differ from those in the NCBI label, but were retained in the PFR² database. In case of multiple attributions of the same sequence to different genetic types by several authors, we retained the molecular taxonomy that was based on the study presenting the most resolved and comprehensive attribution. In addition, 877 unpublished sequences belonging to Orbulina universa, Globigerina bulloides, Neogloboquadrina incompta, Neogoboquadrina dutertrei, Neogloboquadrina pachyderma, and Turborotalita quinqueloba received a genotypic attribution following [START_REF] De Vargas | Molecular evidence of cryptic speciation in planktonic foraminifers and their relation to oceanic provinces[END_REF] and [START_REF] Darling | Molecular evidence links cryptic diversification in polar planktonic protists to Quaternary climate dynamics[END_REF][START_REF] Darling | A resolution for the coiling direction paradox in Neogloboquadrina pachyderma[END_REF][START_REF] Darling | Global molecular phylogeography reveals persistent Arctic circumpolar isolation in a marine planktonic protist[END_REF][START_REF] Darling | The genetic diversity of planktic foraminifera and the global distribution of ribosomal RNA genotypes[END_REF]. Most of these sequences have been produced and identified within earlier studies, but were not originally deposited on NCBI. Their PFR² genotypic assignment is therefore entirely consistent with the attribution of the representative sequences of the same genetic type that were deposited on NCBI.

PFR 2 final taxonomic framework

As a result of the first manual curation and annotation to the genetic type level, the original 3,532 18S rDNA sequences were re-assigned to 33 species names and 2,276 sequences were annotated to the level of genetic types (Supplementary Material 1). For all sequences, we established a ranked taxonomy with six levels: 1-Morphogroup, 2-Genus, 3-Species, 4-Genetic type level 1, 5-Genetic type level 2, 6-Genetic type 3. For the "Morphogroup" rank we used the taxonomical framework of [START_REF] Hemleben | Modern Planktonic Foraminifera[END_REF], dividing the extant planktonic foraminifera species into five clades based on the ultrastructure of the calcareous shell: Spinose, Non-spinose, Microperforate, Monolamellar and Non-spiral. The "Genus" and "Species" ranks follow the primary annotation as described above. For the "Genetic type level 1", "Genetic type level 2"

and "Genetic type level 3" ranks, we used the hierarchical levels presented in the labels of the genetic types of Globigerinoides ruber, Globigerinoides elongatus, Globigerinella siphonifera, Globigerinella calida, Hastigerina pelagica, Globigerina bulloides, Neogloboquadrina dutertrei, Pulleniatina obliquiloculata, and Turborotalita quinqueloba. Genetic type attributions lacking hierarchical structure were reported in the rank "Genetic type level 1". After this step, the Primary Reference Database (Figure 2) of 3,532 sequences contained 113 different taxonomic paths (Supplementary Material 1).

Sequences partitioning into conserved and variable regions

Because PFR² is a resource not only for taxonomic assignment but also for ecological and biogeographical studies, all planktonic foraminiferal 18S rDNA sequences were included irrespective of length, as long as they contained taxonomically relevant information. As a result, the length of the sequences included in the annotated primary database ranges between 33 and 3,412 bp. To evaluate their coverage and information content, all sequences were manually aligned using Seaview 4 [START_REF] Gouy | SeaView version 4: A multiplatform graphical user interface for sequence alignment and phylogenetic tree building[END_REF] to the borders of each variable region of the 18S rDNA fragment. The positions of the borders were determined according to the SSU rDNA secondary structure of the monothalamous foraminifera Micrometula hyalostera presented by [START_REF] Pawlowski | Short rDNA barcodes for species identification in foraminifera[END_REF], except for the region 37/f where a strict homology was difficult to establish for all sequences. Instead, we defined the end of this region by the occurrence of a pattern homologous to the series of nucleotides "CUUUCACAUGA" located at the 3' end of Helix 37. We also noticed that the short conserved fragment located between the variable regions 45/e and 47/f was difficult to identify across all sequences. We thus merged the regions 45/e, 46 and 47/f into a single region that we named 45E-47F (Table1). As a result, the position and length of six conserved (32-37, 37-41, 39-43, 44-45, 47-49, 50) and five variable (37F, 41F, 43E, 45E-47F, 49E) regions were identified for all sequences (Figure 2). The remaining part of the 18S rDNA sequence, only present in sequences EU199447, EU199448 and EU199449 and located before the motive "AAGGGCACCACAAGA" has not been analyzed in this way. All regions fully covered in a sequence and containing sequence motives observed at least twice in the whole dataset were labelled as "complete". Regions fully covered but containing a sequence motive that was observed only once in the whole dataset were labelled as "poor". This is because we consider sequencing/PCR errors as the most likely cause for the occurrence of such unique sequence motives. We realize that using this procedure, even genuine unique sequences may be discarded from the analysis, but this would be the case only if such sequences deviated in all regions. In all other cases, the regions were labelled as "partial" when only a part of the region was present or "not available" if they did not contain any fragment of the sequence. As a result we obtain the Partitioned Primary Reference Database (Figure 2). The coverage of each individual region in the Partitioned Primary Reference Database is given in Supplementary Material 1, and all sequence partitions are given in Supplementary Material 2.

Semi-automated iterative curation pipeline for optimal taxonomic assignment

The consistency of taxonomic assignments within the annotated database of partitioned sequences was assessed using a semi-automated process (Figures 2 and3). All "complete" regions of sequences with the same taxonomic assignment at the morphospecies level were automatically aligned using global pairwise alignment [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF], as implemented in the software needle from the Emboss suite of bioinformatics tools [START_REF] Rice | EMBOSS: The European Molecular Biology Open Software Suite[END_REF]. To detect annotation inconsistencies, mean pairwise similarities were computed for each "complete" region of each sequence against all other sequences with the same taxonomic assignment from the finest annotation level "Genetic type level 3" up to the "Species level" rank.

Results are provided in Supplementary Material 1 and were visualized using R (R Development Core Team, 2014) and the ggplot2 library [START_REF] Wickham | ggplot2: elegant graphics for data analysis[END_REF]. The resulting plots are given in Supplementary Material 3. If all annotations are consistent and there is no variation within taxa, each sequence within the analyzed taxon should only find an exact match and the mean pairwise similarity for that taxon should be 1. However, beyond sequencing/PCR errors introducing spurious sequence differences, there are several reasons why the mean pairwise similarity within a taxon may be lower. First, if a sequence has been assigned the wrong name, its similarity to all other sequences labelled with that name will be low, thus decreasing the resulting mean pairwise similarity. Second, if a sequence has been assigned to the correct taxon, but the taxon comprises multiple sequence motives, that sequence will find a perfect match within the taxon but the mean pairwise similarity will also be lower than 1.

In order to deconvolve the different sources of sequence variability within taxa, we followed a three-step iterative approach, which was repeated for each of the 11 "complete" regions of the analyzed SSU rDNA fragment. First, we considered the distribution of mean pairwise similarities for all sequences within each region assigned to one taxon at the finest rank of "Genetic type level 3". Assuming that misidentifications are rare and result in large pairwise distances, we manually searched for sequences whose mean pairwise similarity deviates substantially from the rest of the sequences within the taxon. Such sequences were initially "invalidated", whereas all other sequences analyzed at this level were "validated". We then repeated the same procedure for the higher ranks of "Genetic type level 2", "Genetic type level 1" and finally "Species level", always starting with the full database (Figures 2 and3A). Thus, at each level, we expected a misidentified sequence to have a pairwise similarity markedly lower than the mean of pairwise similarities between correctly assigned sequences (Figure 3B). This procedure had to be repeated for every rank, because not all sequences in the database are assigned to all ranks. Nevertheless, once "validated", a sequence cannot be "invalidated" during analyses of higher rank taxa, because it represents an accepted variability within that taxon. In taxa where all sequences within a region show low mean pairwise similarities all attributions are initially invalidated (this would be typically the case for a "wastebasket taxa"; Figure 3C).

In the second step, all sequences invalidated during step 1 were reconsidered based on their pairwise similarities with 'validated' sequences from the same region. The main goal of the curated taxonomy being to achieve correct taxonomic assignment at the species level, the pairwise comparison was carried out at this rank. If the best match is a "validated" sequence with the same initial species attribution as the invalidated sequence, this sequence is "validated" at the species level and its assignment at the "genetic type" level is then deleted. Such a situation can only occur when the sequence was initially assigned to the wrong genetic type within the correct species. If the pairwise comparisons of all regions analyzed match sequences with different (but consistent) species attributions than the invalidated sequence, the sequence is reattributed to that species. If the pairwise comparisons indicate that the analyzed sequence has no close relative in the validated part of the database, the initial attribution is retained, provided that the initial attribution is not yet in the validated dataset. This case occurs when all sequences of one species have been initially invalidated because the same species name was associated with highly divergent sequences. When the sequence has no close relative but its initial attribution is represented in the validated part of the dataset, the initial attribution is discarded and the sequence receives an artificial attribution derived from the nearest higher rank that matches the pairwise comparisons. In all cases, the erroneous attributions are replaced by the corrected ones in the database (Figure 2, Supplementary Material 1).

In the third step, sequences that received new attributions were reanalyzed as described in step 1.

If inconsistencies in the distribution of mean pairwise similarities remain, steps 2 and 3 are repeated until no inconsistency is observed.

As a final diagnosis we performed leave-one-out analyses to evaluate the robustness and potential limitations of the curated taxonomy, as well as a monophyly validation by Neighbor-Joining using only sequences that are covering the 6 conserved and 5 variable regions of the 5' end fragment. First, each individual sequence included in the first version of PFR² was blasted against the remaining part of the database including n-1 sequences using SWIPE [START_REF] Rognes | Faster Smith-Waterman database searches with inter-sequence SIMD parallelisation[END_REF].

The sequences among the "n-1 PFR² database" returning the highest score were retrieved and their taxonomic attribution compared to the one of the blasted sequence (Supplementary Material 1). Second, we retrieved all sequences covering the 5 variable and 6 conserved regions and divided them according to their assignment to higher taxa (here simplified by the morphogroups Monolamellar, Non-Spinose, Spinose, and Microperforates + Benthic). Each subset was automatically aligned using MAFFT v.7 [START_REF] Katoh | MAFFT multiple sequence alignment software version 7: improvements in performance and usability[END_REF] and the subsequent alignments were trimmed off on the edges to conserve only homologous position, finally leading to 41, 583, 271, and 100 analyzed sequences for the Monolamellar, Non-Spinose, Spinose, and Microperforates + Non-spiral morphogroups, respectively. For each alignment, a tree was inferred using a Neighbor-Joining approach with Juke and Cantor distance while taking into account gap sites as implemented in SEAVIEW 4 (Supplementary Material 4) with 100 pseudoreplicates. The scripts used to perform the different curation steps are available as Supplementary Material 5.

Results

Of the 3,532 planktonic foraminiferal 18S rDNA partial sequences analyzed, 3,347 (94.8%) contained at least one "complete" gene region making possible the curation process. The remaining 185 sequences included 33 singletons (rare motives or poor quality sequences) and Amongst the 3,347 curated sequences, the taxonomic assignment of 84 was initially invalidated.

Of these, 3 represent cases where the morphospecies attribution was correct, but the attribution to a genetic type was erroneous. In 46 cases, the invalidated sequences found a perfect match with a different taxon and thus their taxonomic assignment was changed. In all of these cases, the novel taxonomic assignment corresponded to a morphologically similar morphospecies, explaining the original misidentification of the sequenced specimen. In 14 cases, the original assignment was retained because the sequences did not find any match and their original attribution did not appear in the validated part of the dataset. All of these sequences were labelled as Hastigerinella digitata. This species name had been entirely invalidated in the first step because of inconsistent use of the homonymous species name Beella digitata. Finally, 17 sequences received an unresolved artificial assignment. These represent six different sequence motives diverging substantially from all sequences in the validated part of the database and also between each other. Because the original attribution upon collection was obviously wrong, we could not reassign these sequences to the species level. In two cases, we could identify the most likely generic attribution, but four sequences are left with an entirely unresolved path. Finally, our procedure captured one sequence with a spelling error in its path and three sequences that appear to have been attributed correctly but represent small variants within species. After resolution of the 84 conflicts described above, the re-annotated dataset was subjected to a second round of the curation process for verification. All sequences were validated.

Based on this internally consistent taxonomic annotation for all 3,347 18S rDNA sequences from individual planktonic foraminifera, we generated the Planktonic Foraminiferal Ribosomal Reference or PFR 2 database. Of the 3,347 sequences, 25 were shorter than 200 bp, and could not be deposited in NCBI (see Supplementary Material 1). The PFR 2 1.0 database thus includes 3,322 reference sequences assigned to 32 morphospecies and 6 taxa with unresolved taxonomy (Figure 2), and contains 119 unique taxonomic paths when including all three levels of genetic types.

The leave-one-out BLAST evaluation applied on the first version of PFR² to assess its robustness returned an identical taxonomic path for 2,509 sequences. For 614 sequences, the BLASTdetermined taxonomic paths were identical between the "morphogroup" and "species" rank but displayed a different resolution between the ranks "genetic type level 1" and "genetic type level 3". This reflects a situation where some sequences belonging to one species are annotated to the level of a genetic type, whereas others are not. Finally, 19 sequences were assigned to the correct species but to a different genetic type. This illustrates the case of genetic types represented by only one sequence in the database, which were logically assigned to the closest genetic type within the same species by the leave-one-out procedure. Thus, 94.5 % of the sequences in the PFR 2 database find a nearest neighbor with a correct taxonomic assignment at the species target level. For the remaining 180 sequences, the returned taxonomic path was inconsistent at the species level. In two cases, the sequences were assigned to a morphologically and phylogenetically close sister species (Globorotalia ungulata and Globorotalia tumida), reflecting insufficient coverage in the database for these species. Two cases involved singleton sequences with unresolved taxonomy, which find no obvious nearest neighbor. Finally, 176 cases of inconsistent identification refer to sequences of Globigerinella calida and Globigerinella siphonifera, whose species names have been used interchangeably in the literature [START_REF] Weiner | Phylogeography of the tropical planktonic foraminifera lineage Globigerinella reveals isolation inconsistent with passive dispersal by ocean currents[END_REF] and the clade has been shown to be in need of a taxonomic revision [START_REF] Weiner | Genetic and morphometric evidence for parallel evolution of the Globigerinella calida morphotype[END_REF].

The leave-one-out evaluation thus reveals excellent coverage of PFR 2 and confirms that the curated taxonomy is internally entirely consistent.

To further confirm the validity of morphospecies level taxonomy, we constructed NJ trees for the five clades including only the long sequences (Supplementary Material 4). This analysis confirmed the monophyly of all morphospecies, except the Globigerinella calida/Globigerinella siphonifera plexus. All clades were strongly supported except for the sister species Globorotalia tumida and Globorotalia ungulata and the monolamellar species Hastigerina pelagica and Hastigerinella digitata. In the first case, the poor support reflects the lack of differentiation between these two species in the conserved region of the gene, thus decreasing the bootstrap score; in the second case the extreme divergence of two genetic lineages of Hastigerina pelagica renders the phylogenetic reconstruction difficult [START_REF] Weiner | Vertical niche partitioning between cryptic sibling species of a cosmopolitan marine planktonic protist[END_REF]).

An analysis of the taxonomic annotations retained in PFR 2 reveals that the database covers at least 70-80% of the traditionally recognized planktonic foraminiferal species in each clade. The species represented in PFR 2 constitute the dominant part of planktonic foraminifera assemblages in the world oceans. Compared with a global database of census counts from surface sediments (MARGO database, [START_REF] Kucera | Reconstruction of sea-surface temperatures from assemblages of planktonic foraminifera: multi-technique approach based on geographically constrained calibration data sets and its application to glacial Atlantic and Pacific Oceans[END_REF], the species covered by PFR 2 account for >90% of tests larger than 150 μm found in surface sediments (Figure 4). In cold and temperate provinces, PFR 2 species account for almost the entire assemblages, while in warmer subtropical and tropical waters, only up to 4% of the sedimentary assemblages are not represented in PFR 2 . Evidently, PFR 2 reference sequences cover most of the ecologically relevant portion of the morphological diversity and the taxa that are not yet represented in PFR 2 are small, rare or taxonomically obscure. It is possible that some of these taxa may correspond to the six sequences with still unresolved taxonomy. If so, PFR 2 may be considered to cover up to 38 of the 47 recognized species.

Finally, for each species present in PFR 2 , we evaluated the ecological coverage of the global sampling effort (Figure 4). Morphospecies of planktonic foraminifera are known to be distributed zonally across the world oceans, reflecting the latitudinal distribution of sea surface temperature (e. g., [START_REF] Bé | Distribution and ecology of living planktonic foraminifera in surface waters of the Atlantic and Indian Oceans[END_REF]. A comparison between the temperature range of each species as indicated by their relative abundance in surface sediment samples [START_REF] Kucera | Reconstruction of sea-surface temperatures from assemblages of planktonic foraminifera: multi-technique approach based on geographically constrained calibration data sets and its application to glacial Atlantic and Pacific Oceans[END_REF] and the temperatures measured at sampling localities shows that a large portion of the ecological range of the species is covered by the reference sequences in PFR 2 (Figure 4).

The PFR 2 web interface

To facilitate data download and comparative sequence analyses, PFR² has been implemented into a dedicated web interface, available at http://pfr2.sb-roscoff.fr. The website provides:

(1) a search/browse module, which allows the user to download parts of the database either by taxonomic rank (morphogroup name, genus name, species name), geographic region (e. g., North Atlantic, Mediterranean Sea, Indian Ocean) or collection (cruise name) ;

(2) a classical BLAST/Similarity module that facilitates identification of unknown sequences;

(3) a map module displaying the localities for all sequences present in PFR² and facilitating download of all data from each single locality;

(4) a download section with direct access to all data included in PFR². All sequences and sequence partitions are available in FASTA format and the metadata are available in a tabulated file.

Discussion

Comprehensive databases of ribosomal RNA sequences with curated taxonomy are available for Protists (Protist ribosomal reference database, PR²; [START_REF] Guillou | The Protist Ribosomal Reference database (PR2): a catalog of unicellular eukaryote small sub-unit rRNA sequences with curated taxonomy[END_REF] and for the major domains of life (SILVA; [START_REF] Yilmaz | The SILVA and "All-species Living Tree Project (LTP)" taxonomic frameworks[END_REF]. These databases include sequences of planktonic foraminifera. However, they are used mainly as benchmarks to annotate complex environmental datasets (e.g., [START_REF] Logares | Patterns of rare and abundant marine microbial eukaryotes[END_REF] at the morphological species level. In contrast, PFR² has been designed and implemented in a way that facilitates other applications.

First, because of structural limitations PR² contains "only" 402 sequences of planktonic foraminifera (based on Released 203 of GenBank, October 2014), compared to PFR², which contains for now 3,322 SSU rDNA sequences. Second, 2276 of the sequences present in PFR² have an assignation to the genetic type level and as far as possible, the sequences are associated with metadata related to the origin of each specimen and the conditions where it was collected, thus forming a basis for ecological modelling. Third, most importantly, using planktonic foraminifera as a case study, we propose and implement an annotation scheme with unmatched accuracy and full tracking of changes. This is only possible because of the narrower focus of PFR² combined with high-level expert knowledge of their taxonomy. The fidelity of the annotations will facilitate a qualitatively entirely different level of analysis of eDNA libraries.

For example, the design of PFR² allows to incorporate advances in classical and molecular taxonomy, particularly at the level of genetic types (e.g., [START_REF] André | SSU rDNA Divergence in Planktonic Foraminifera: Molecular Taxonomy and Biogeographic Implications (V Ketmaier[END_REF], which can be reevaluated depending of the criteria used to delineate molecular OTUs. Further, by retaining information on clone attribution to specimens (vouchers), PFR² allows to evaluate intra-genomic polymorphism, which offers excellent opportunity to identify the taxonomically relevant level of variability [START_REF] Weber | Wide occurrence of SSU rDNA intragenomic polymorphism in foraminifera and its implications for molecular species identification[END_REF]. Finally, the modular structure of PFR 2 (i.e., its partitioning into variable and conserved regions) is particularly suitable for the evaluation of existing barcodes or the design of new barcoding systems needed to capture total or partial planktonic foraminiferal diversity within complex plankton assemblages. Indeed, an examination of the length polymorphism in the 11 regions of the 18S rDNA fragment that have been aligned for all PFR 2 sequences reveals that next to the variable 37/f region identified as a barcode for benthic foraminifera [START_REF] Pawlowski | Short rDNA barcodes for species identification in foraminifera[END_REF], several other regions may be suitable as targets for barcoding of planktonic foraminifera (Figure 5).

The main difference between PFR² and classical databases is in the association of sequence data with environmental and collection data. Such level of annotation is not feasible in large databases, which have to rely on the completeness and level of metadata details provided in GenBank. The association of metadata to PFR² sequences facilitates an assessment of biogeography and ecology of genetic types (potential cryptic species). This is significant for studies of evolutionary processes in the open ocean such as speciation and gene flow at basin scale, but also for paleoceanography, which exploits ecological preferences of planktonic Foraminferal species to reconstruct climate history of the Earth (e.g., [START_REF] Kucera | Reconstruction of sea-surface temperatures from assemblages of planktonic foraminifera: multi-technique approach based on geographically constrained calibration data sets and its application to glacial Atlantic and Pacific Oceans[END_REF].

Modeling studies showed that the integration of cryptic diversity into paleoceanographic studies will improve their accuracy [START_REF] Kucera | Cryptic species of planktonic foraminifera: their effect on palaeoceanographic reconstructions[END_REF][START_REF] Morard | Ecological modeling of the temperature dependence of cryptic species of planktonic foraminifera in the Southern Hemisphere[END_REF]. Together with the MARGO database [START_REF] Kucera | Reconstruction of sea-surface temperatures from assemblages of planktonic foraminifera: multi-technique approach based on geographically constrained calibration data sets and its application to glacial Atlantic and Pacific Oceans[END_REF], which records the occurrence of morphospecies of planktonic foraminifera in surface sediments and the CHRONOS/NEPTUNE database [START_REF] Spencer-Cervato | How synchronous are neogene marine plankton events?[END_REF]; http://www.chronos.org/), which records their occurrence through geological time, PFR² represents the cornerstone to connect genetic diversity to the fossil record in an entire group of pelagic protists.

Conclusion and perspectives

The PFR 2 database represents the first geographically and taxonomically comprehensive reference barcoding system for an entire group of pelagic protists. It constitutes a pivotal tool to investigate the diversity, ecology, biogeography, and evolution in planktonic foraminifera as a model system for pelagic protists. In addition, the database constitutes an important resource allowing reinterpretation and refinement of the use of foraminifera as markers for stratigraphy and paleoceanography. In particular, PFR 2 can be used to: (i) annotate and classify newly generated 18S rDNA sequences from single individuals; (ii) study the biogeography of cryptic genetic types; (iii) design rank-specific primers and probes to target any group of planktonic foraminifera in natural communities; and (iv) assign accurate taxonomy to environmental sequences from metabarcoding or metagenomic datasets. This last point is particularly worth noting. Indeed, future global metabarcoding of planktonic foraminifera covering comprehensive spatio-temporal scales will likely reveal the full extent and complexity of species diversity and ecology in this group, serving as a model system for studies of the evolutionary dynamics of the plankton and its interaction with the Earth system. Figure 2 Workflow to constitute PFR². In step I the sequences, metadata and taxonomic information are retrieved from public databases and literature or from the internal databases of the co-authors to constitute the Primary Reference Database. In step II, the coverage of each sequence is evaluated by alignment with structural regions of the 18S RNA secondary structure derived for the species Micrometula hyalostera [START_REF] Pawlowski | Short rDNA barcodes for species identification in foraminifera[END_REF]. In step III, the consistency of the annotation is checked from the most exclusive level of annotation "genetic type 3" up to the species level (Phase 1) to detect annotation inconsistencies (See Figure 3). Sequences with wrong annotation are invalidated, compared to the validated part of the dataset (Phase 2) and reannotated depending on the best hit out of the valid dataset. The consistency of all annotations is then checked again following the same procedure as in Phase 1 (Phase 3), to ensure that no taxonomic inconsistency remains. In step IV, all sequences which have been subjected to the curation process are integrated in the Planktonic Foraminifera Ribosomal Reference database (PFR²). The results of all steps are given in Supplementary Material 1.

Figures

Figure 3

Annotation inconsistency detection. The procedure followed to identify annotation inconsistencies is exemplified by three cases. Each graph represents variability in pairwise similarities observed across each region of all sequences sharing the same annotation level. The names of the taxon and annotation level are given above the plot with the number of sequences in parenthesis. Each vertical line represents one region with the variability represented as box plot, the number of "complete" regions is given at the bottom of the line. The case "A" describes the annotation validation process starting from the most exclusive rank of "genetic type level 3" to the "species" rank. After the validation at one rank level, the sequences with valid annotation are merged into a taxonomic unit of a higher rank, this now including multiple sequence motifs which decreases the average similarity level of each region, thus leading to higher variability in higher ranks. Case "B" represents the occurrence of obvious outliers at the species level, which are invalidated. Case "C" represents the co-occurrence of divergent sequences under the same taxonomic attribution, which are consequently all invalidated. Box plots for all ranks can be found in Supplementary Material 3 and the pairwise similarities calculated for each taxonomic level are given in Supplementary Material 1. [START_REF] Kucera | Reconstruction of sea-surface temperatures from assemblages of planktonic foraminifera: multi-technique approach based on geographically constrained calibration data sets and its application to glacial Atlantic and Pacific Oceans[END_REF] and the mean annual SST (MODIS Aqua, NASA, Greenbelt, MD, USA). The grey dots highlight the mean annual SST at the location where the living planktonic foraminifera yielding sequences were sampled.

The number of sequences available for each species as well as the number of taxonomic paths above the species level is shown next to the graphs. Also shown is the cumulative mean relative abundance in the sediments of all species included in PFR² plotted against the mean annual SST in discrete 1°C intervals. Vertical bars represent 95% confidence intervals for each 1°C bin.

Scientific cruise

Near-Shore Collection In step I the sequences, metadata and taxonomic information are retrieved from public databases and literature or from the internal databases of the co-authors to constitute the Primary Reference Database. In step II, the coverage of each sequence is evaluated by alignment with structural regions of the 18S RNA secondary structure derived for the species Micrometula hyalostera [START_REF] Pawlowski | Short rDNA barcodes for species identification in foraminifera[END_REF]. In step III, the consistency of the annotation is checked from the most exclusive level of annotation "genetic type 3" up to the species level (Phase 1) to detect annotation inconsistencies (See Figure 3). Sequences with wrong annotation are invalidated, compared to the validated part of the dataset (Phase 2) and re-annotated depending on the best hit out of the valid dataset. The consistency of all annotations is then checked again following the same procedure as in Phase 1 (Phase 3), to ensure that no taxonomic inconsistency remains. In step IV, all sequences which have been subjected to the curation process are integrated in the Planktonic Foraminifera Ribosomal Reference database (PFR²).

The results of all steps are given in Supplementary Material 1. . Annotation inconsistency detection. The procedure followed to identify annotation inconsistencies is exemplified by three cases. Each graph represents variability in pairwise similarities observed across each region of all sequences sharing the same annotation level. The names of the taxon and annotation level are given above the plot with the number of sequences in parenthesis. Each vertical line represents one region with the variability represented as box plot, the number of "complete" regions is given at the bottom of the line. The case "A" describes the annotation validation process starting from the most exclusive rank of "genetic type level 3" to the "species" rank. After the validation at one rank level, the sequences with valid annotation are merged into a taxonomic unit of a higher rank, this now including multiple sequence motifs which decreases the average similarity level of each region, thus leading to higher variability in higher ranks. Case "B" represents the occurrence of obvious outliers at the species level, which are invalidated. Case "C" represents the co-occurrence of divergent sequences under the same taxonomic attribution, which are consequently all invalidated. Box plots for all ranks can be found in Supplementary Material 3 and the pairwise similarities calculated for each taxonomic level are given in Supplementary Material 1. [START_REF] Kucera | Reconstruction of sea-surface temperatures from assemblages of planktonic foraminifera: multi-technique approach based on geographically constrained calibration data sets and its application to glacial Atlantic and Pacific Oceans[END_REF] and the mean annual SST from the World Ocean Atlas (Locarnini, 2005). The grey dots highlight the mean annual SST at the location where the living planktonic Foraminifera yielding sequences were sampled. The number of sequences available for each species as well as the number of taxonomic paths above the species level is shown next to the graphs. Also shown is the cumulative mean relative abundance in the sediments of all species included in PFR² plotted against the mean annual SST in discrete 1°C intervals. Vertical bars represent 95% confidence intervals for each 1°C bin.

Figure 1

 1 Figure 1 Sampling Map. Location of the 460 oceanic stations sampled over 20 years for single-cell genetic studies of planktonic foraminifera. Each symbol corresponds to a scientific cruise or near shore collection site. Cruise names and dates of the collection expeditions are indicated in the legend. Grey shading shows ocean bathymetry.

Figure 4

 4 Figure 4Taxonomic and ecological coverage of PFR². For each morphogroup (Spinose, Non-Spinose, Microperforates, Monolamellar and Non-Spiral) the number of species included in PFR² is given in the filled bar while the number of species not present is indicated in the adjacent open bar. The relative abundance in the sediments of each species included in PFR² is given in a log-scale value against mean Sea Surface Temperature (SST) at the sampling station. Relative abundances in sediments are derived from the MARGO database[START_REF] Kucera | Reconstruction of sea-surface temperatures from assemblages of planktonic foraminifera: multi-technique approach based on geographically constrained calibration data sets and its application to glacial Atlantic and Pacific Oceans[END_REF] and the mean annual

Figure 5

 5 Figure 5 Length polymorphism. Each rectangle represents the length polymorphism within each region of the analyzed 18S rDNA fragment across all resolved taxonomic units in PFR². The regions are based on the rRNA secondary structure and are named following Pawlowski and Lecroq (2010).

Figure 1 .Figure 2 .

 12 Figure 1. Sampling Map. Location of the 460 oceanic stations sampled over 20 years for single-cell genetic studies of planktonic Foraminifera. Each symbol corresponds to a scientific cruise or near shore collection site. Cruise names and dates of the collection expeditions are indicated in the legend. Grey shading shows ocean bathymetry.

Figure 4 .

 4 Figure 4. Taxonomic and ecological coverage of PFR².For each morphogroup (Spinose, Non-Spinose, Microperforates, Monolamellar and Non-Spiral) the number of species included in PFR² is given in the filled bar while the number of species not present is indicated in the adjacent open bar. The relative abundance in the sediments of each species included in PFR² is given in a log-scale value against mean Sea Surface Temperature (SST) at the sampling station. Relative abundances in sediments are derived from the MARGO database[START_REF] Kucera | Reconstruction of sea-surface temperatures from assemblages of planktonic foraminifera: multi-technique approach based on geographically constrained calibration data sets and its application to glacial Atlantic and Pacific Oceans[END_REF] and the mean annual SST from the World Ocean Atlas(Locarnini, 2005). The grey dots highlight the mean annual SST at the location where the living planktonic Foraminifera yielding sequences were sampled. The number of sequences available for each species as well as the number of taxonomic paths above the species level is shown next to the graphs. Also shown is the cumulative mean relative abundance in the sediments of all species included in PFR² plotted against the mean annual SST in discrete 1°C intervals. Vertical bars represent 95% confidence intervals for each 1°C bin.

Table 1 .

 1 Flanking conserved sequences of the 5 variable regions in planktonic foraminifera. The minimum and maximum lenght of each region are given as well as their coverage in the database (See details in the text).

	Region	Specificity	Begining	End	Min lenght Max lenght Not available Partial Poor Complete
	32-37	Eukaryotes	-	-	-	-	949	2583	0	0
	37F	Foraminifera 5'-GGAUUGACA CUUUCACAUGA-3'	38	132	800	272	249
	37-41	Eukaryotes	-	-	68	72	547	403	138
	41F	Foraminifera	5'-AAUUGCG	GCAACGAA-3'	58	322	349	346	282
	39-43	Eukaryotes	-	-	27	29	460	34	57
	43E	Eukaryotes	5'-CUUGUU	AACUAGAGGG-3'	33	195	401	263	265
	44-45	Eukaryotes	-	-	113	123	487	1288 136
	45E-47F Euk -Forams	5'-CAGUGAG	GGUGGGG-3'	179	312	1660	187	386
	47-49	Eukaryotes	-	-	140	148	1827	425	152
	49E	Eukaryotes	5'-GUGAG	CGAACAG-3'	27	127	2251	130	125
	50	Eukaryotes	-	-	-	-	2389	1143	0	0

Acknowledgments:

We thank all crew members and scientist for their help in the collection of planktonic foraminifera that were used to generate the database. We thank Erica de Leau for her help in compiling the data and Dominique Boeuf and ABIMS for their help in designing and hosting the PFR² website. This work was supported by grants from ANR-11-BTBR-0008 OCEANOMICS, ANR-09-BLAN-0348 POSEIDON, ANR-JCJC06-0142-PALEO-CTD, from Natural Environment Research Council of the United Kingdom (NER/J/S2000/00860 and NE/D009707/1), the Leverhulme Trust and the Carnegie Trust for the Universities of Scotland, from DFG-Research Center/Cluster of Excellence "The Ocean in the Earth System" and from the Deutsche Forschungsgemeinschaft KU2259/19 and DU1319/1-1. This study is a contribution to the effort of the SCOR/IGBP Working Group 138 "Modern planktonic Foraminifera and ocean changes".

Data Accessibility

Sequences, NCBI accession numbers and metadata are available in Supplementary Material 1 and 2 and on the PFR² website at http://pfr2.sb-roscoff.fr. The custom scripts used to perform the curation procedure are available in Supplementary Material 5; the results of the curation process are given in Supplementary Material 1 and 2.

Author contribution

KFD, CdV, YU, RM, TdG, AKMW, HAS, MK, AA, MS participated in sample collection, CdV, MK, KFD, CMW, CJD, FQ, GE, TdG provided laboratory infrastructure, KFD, YU, RM, AKMW, AA, HAS participated in laboratory work. FM and RM conceived and designed the bioinformatics pipeline, FM performed the computational work, SA built the website. RM wrote the manuscript with help from MK and CdV. All authors read, edited and approved the final manuscript.

Supplementary Material.

Supplementary Material 1

Information on all consecutive steps followed to constitute the PFR². All fields are explained in the file.

Supplementary Material 2

FASTA files of sequences used to build the PFR². FASTA files are provided for the full sequences and individual partitions.

Supplementary Material 3

Box plots showing pairwise similarities for each taxonomic level. See Figure 3 for explanations of the content of the plots.

Supplementary Material 4

Neighbor-joining trees showing the monophyly of each morphospecies present in PFR².

Supplementary Material 5

Custom scripts used to perform the different curation steps.


Karam Naser

Vincent Ricordel

Patrick Le Callet
email: patrick.le-callet@univ-nantes.fr

PERFORMANCE ANALYSIS OF TEXTURE SIMILARITY METRICS IN HEVC INTRA PREDICTION

The visual signal is highly occupied by regions of homogeneous and repetitive patterns known as Textures. Textures have a common property that their similarity highly deviates from point by point comparison, i.e, two textures can look very similar even if they have some shift, rotation and difference in their distribution.

In the context of compression, All of the MPEG reference encoders (including HEVC), aim at minimizing the bitrate at a certain distortion level measured in terms of pixel comparison. For textures, this kind of distortion measure does not usually reflect the amount of perceived distortion. For this reason, we investigate the use of state of the art perceptual similarity metrics as a replacement for this measure. In other words, we aim at optimization the bitrate such that we minimize the perceptual distortion rather that the pixels difference. We used two metrics (Local Radius Index and Structure Texture Similarity Metric) in selecting the best intra prediction mode and block partitioning. Experimental results showed that these metrics try always to retain some structural properties of the textures. These metrics also showed a better rate-distortion performance when the distortion is measured via a distance metric based on texture features.

INTRODUCTION

Visual signal is a rich source of information. It can be decomposed into different classes and components. In this work, we consider it as two components, structure and texture. Structure contains the semantic meaning of the scene (edges, lines, corners, etc.) and textures fill the gap between structures. According to this, Textures range from very simple ones, like a DC block, to more complex ones. They can be also classified as regular and stochastic. Each texture type possesses different spectral, statistical and perceptual properties. For this reason, encoding them without considering their properties does not end up with the best rate-quality performance, which is the optimal goal of video compression.

In terms of visual similarity, humans are less sensitive to variations in textures, i.e, two textures can look very similar even if they have some deviations in scale, orientation and repetition. In contrast, they can look very different even if they have small average pixel difference. Therefore, assessing the texture similarity by the means of pixel comparison is avoided. Instead, the comparison can be done between structural information (encoded in frequency subband channels) and/or statistical information. this problem is an active research topic in both engineering and psychology. Details of different texture similarity metrics can be found in [START_REF] Pappas | Image analysis: Focus on texture similarity[END_REF].

Video coding standards, such as HEVC, aim at minimizing the bitrate within a certain distortion level. The measured distortion is typically the mean squared error. This type of distortion, which is based on comparing pixel values, does not proportionally reflect the amount of perceived distortion (especially for texture components). For this reason, many approaches consider different weight for the distortion computed at each block according to its texture complexity (ex. [START_REF] Yu | A perceptual bit allocation scheme for h. 264[END_REF]). Others consider replacing this measure by more appropriate ones (ex. [START_REF] Ou | Ssim-based perceptual rate control for video coding[END_REF]). In contrast, some approaches try to replace some textures with a synthesized ones which look visually very similar (ex. [START_REF] Zhang | A parametric framework for video compression using region-based texture models[END_REF] [START_REF] Balle | Models for static and dynamic texture synthesis in image and video compression[END_REF]).

In this paper, we investigate the possibility of encoding static textures in such a way that the encoder try to minimize the bitrate while maximizing visual similarity between the reconstructed and the original signal. We borrow two similarity metrics from the texture retrieval problem. The two metrics, known as Local Radius Index (LRI) and Structure Texture Similarity Index (STSIM), are the most recent and successful ones in texture retrieval. We used them to select the best intra prediction mode and block partitioning. By doing so, we do not violate the HEVC standard, the HEVC decoder can thus be directly used to decode the resulting bitstream.

The rest of the paper is organized as follows: Section 2 gives an overview of the texture similarity metrics used in this work. Section 3 presents the procedure carried out to evaluate each metric. In Section 4.1, the experimental results are provided and discussed with a conclusion given in Section 5.

OVERVIEW OF TEXTURE SIMILARITY METRICS

Texture similarity metrics exist in various forms. Some of them compare the statistics of textures in the spatial domain and others in the subband frequency domain (details can be found in [START_REF] Pappas | Image analysis: Focus on texture similarity[END_REF]. In this paper, we consider LRI and STSIM as being recent and successful texture similarity metrics.

STSIM

STSIM was presented in [6] and further improved in [START_REF] Zujovic | Structural texture similarity metrics for image analysis and retrieval[END_REF]. It is based on comparing set of statistics in the subband decomposition. These statistics consist of mean, standard deviation, and horizontal and vertical auto-correlation of each subband. Beside that, it computes also the cross correlation between subbands with the same scale or subbands with the same orientation. This set of statistics provides a solid description of a given texture and thus can well characterize the similarity between two textures.

LRI

LRI is a successive to STSIM. LRI is much less computationally expensive as compared to STSIM and performs better in the context of similar texture retrieval [START_REF] Zhai | Local radius index -a new texture similarity feature[END_REF]. It computes local index for each pixel in the spatial domain, beside that, it also computes the local binary pattern, standard deviation of each subband in the subband frequency domain and an intensity penalization term. Thus it is a combination between the analysis in frequency decomposition and spatial domain.

Adaptation of Similarity Metrics in HEVC

STSIM

One interesting property of STSIM is that it is bounded between one and zero, where one is the maximum similarity index. Embedding this metric in HEVC is straightforward.

As we want to use it as a distortion metric, we consider the distortion function as:

D ST SIM = 1 -I ST SIM (1)
where I ST SIM is the STSIM similarity index.

LRI

LRI, in contrast, is not bounded. It computes information divergence between the distribution of the local indexes of two textures. This divergence, as well as the one used inside the local binary pattern, can easily have infinite value when used as block based distortion measure (specially for small blocks comparison) as both of them compute the logarithm of the probability distribution. To overcome this problem, we modified both of them in the following manner:

• We use a normalized Kullback Leiber Divergence (KLD) in LRI. Since the Gibbs inequality assures that KLD value cannot accede the entropy of the first element of two compared distributions, we normalize KLD by dividing it by this entropy. The only exception is when any compared distributions have a probability value of zero when the other does not. In this case, we assume that the two distributions are different and we return the maximum value which is one.

• We normalize also the the local binary pattern term which computes the log-likelihood function between two distributions. The function is divided it by the entropy of the second distribution. The similarity function is then obtained by subtracting the normalized LBP from one. Similarly, we return the maximum value (one) when the distributions are assumed to be different.

We also eliminate the intensity penalization term in our work. This term was initially added to penalize the differences in the intensity values. Since we normalized LRI, the IP term has a greater impact on the overall metric and causes the metric to be more pixel dependent, which is far away from our goal.

PERFORMANCE EVALUATION

To evaluate the performance of the metrics, we implemented both STSIM and LRI and integrated them in HEVC encoder. These metrics were used inside the cost function for selecting the best intra prediction mode and block partitioning. To elaborate more, the details of HEVC intra prediction mode selection is given below:

Intra Mode Selection in HEVC reference Encoder

HEVC defines 33 possible directional prediction modes. Beside this, it also defines DC and Planar prediction modes.

For each mode, the residual block can be directly encoded (transform, qunatization and entropy coding) or further partitioned into quadtree. Each time, the cost function is computed. The particular mode and partition which minimizes the cost is then selected.

In the reference encoder, HEVC considers only 3 most probable modes for the full rate-distortion optimization. These 3 modes are the ones among the 35 modes that have a minimum Sum of Absolute Transformed Difference (SATD) between the original and prediction blocks. For selecting the best mode, HEVC computes the Sum of Squared Difference (SSD) as a distortion metric in the cost function.

Replacing HEVC Distortion Metric

In this work, we replaced the SATD and SSD of HEVC reference encoder by distortion metrics evaluated with perceptual similarity metrics. We experimented both STSIM and LRI and adapted in HEVC as described in 2.3. To keep the range of SATD and SSD, the metrics were multiplied by:

• 255 × N × N for SATD • 255 2 × N × N for SSD.
where N is the dimension of the prediction block.

EXPERIMENTS AND RESULTS

We have experimented the use the two metrics in HEVC for coding static textures. We used Brodazt textures downloaded from USC-SIPI dataset [START_REF]USC-SIPI Dataset[END_REF]. This contains 13 different gray scale textures (see Fig. 1) which are extensively used in textures analysis for engineering and psychophysical experiments. We used HM 9.0 [START_REF]/WG[END_REF] as a host encoder. In the following subsections, we provide the details of each experiment.

Quality of the Decoded Textures

To study the impact of using perceptual similarity metrics in HEVC, we encoded all textures for different QP values.

First we compared the visual quality when using the default metric or the perceptual ones. For low QP values (low compression), no change has been observed. For high QP, however, significant changes can be seen. Examples of this are shown in the figures below.

The first example is shown in Fig. 2. We can see that encoding a highly structured texture with a large QP value results in loosing most of its semantics. This is because many blocks are replaced by DC values. Using either one of the metrics can retain the overall structure of the texture. One can also notice that there exists many wrong directions, but the overall quality is much more better.

Another example is shown in Fig. 3. In this figure, the effect of wrong prediction direction is more clear when LRI is used. On the other hand, the right part of the texture is completely eliminated when the default metrics are used.

Fig. 4 gives an example of a extensive compression. In this case, the details of the texture are lost in all three coding configuration. But one thing to notice that there is different kind of compression artifacts when the perceptual metrics are used. We noticed that whatever hard the compression is, the perceptual metrics try to provide a structure of a textures rather than simple DC blocks. This can be particularly interesting if we consider the no reference quality of the decoded texture. This means that the original texture is somehow replaced by another texture which looks more natural than when the default metrics are used.

Encoder Partitioning Behavior

To understand more the effects of each metric on the prediction mechanism, we measured the number of splitting depths that the encoder uses to encode each texture. We ran a simulation on the 13 textures from the used dataset with a quantization parameter (QP) taken range of [22,27,32,37,43,47,51]. The corresponding histograms are shown in Fig. 5. The splitting depth of zero means that the prediction block has its maximal size (64x64). Increasing the splitting depth by one corresponds to partition the block into four sub-blocks. In this figure, the histograms were scaled by the number of smallest blocks (4x4) that the corresponding splitting depth contains. This was done to have a fair comparison between splitting depths as each splitting occupies different areas of the frame. We observe from these histograms that when the default metrics are used, the encoder uses small prediction blocks for low compression (low QP) to better approximate the input signal. For high compression, it tries to approximate large prediction blocks (mostly with DC values) to have better compression. The behavior is totally different when LRI or STSIM is used. The encoder behavior does not change much as the compression changes. It uses always large block sizes to approximate the input signal and small block sizes (less than 16x16) are rarely chosen. The is because these metrics compare statistics of different distributions. For small block sizes, there is always a lack of enough statistics and usually these metrics return a high value of distortion in such a condition.

Rate Distortion Analysis

Until now, we considered only the quality of the decoded textures. In video coding, however, the comparison between two different schemes is a rate-distortion based. For this work, the rate-distortion comparison is not straightforward. This is because, up to our knowledge, no reliable quality metric has been designed designed for textures. The usual approach of comparing PSNR value does not provide a useful information as it compare pixel values, which is far away from the purpose of this work. Beside that, the type of distortion when the perceptual metrics is used is totally different and cannot be assessed by PSNR (ex. see Fig 2).

Again, we follow the same approach of borrowing a distance metric from texture retrieval problem. We avoid using LRI or STSIM as this may result in biased assessment. We used another metric [START_REF] Manjunath | Texture features for browsing and retrieval of image data[END_REF] which is based on comparing features of textures in the frequency domain. This metric is compared to both LRI and STSIM ( [START_REF] Zhai | Local radius index -a new texture similarity feature[END_REF]) and provides close by results in retrieval rate. This metric compares the energy and mean of frequency subbands (computed using Gabor filters). The metric was downloaded from the authors website and used as a distortion metrics in our work.

By calculating the distance measure by this metric to What we observe from these is that in most cases, LRI and STSIM provides better score than the default metrics in the low bitrate region. For high bitrate region, no gain when using LRI or STSIM. This can prove that even for high compression, both the similarity metrics try to keep the structural information rather than replacing the texture blocks with DC blocks.

CONCLUSION

In this paper, we studied the effect of replacing some of the default distortion metrics in HEVC by perceptual ones. These metrics were used to select the best prediction mode and block partitioning. We used recent perceptual similarity metrics, namely LRI and STSIM which have been devel-oped in the context of texture retrieval. The direct benefit of this approach is that it is its compatibility with HEVC standard, which means no modification to the decoder is needed. When these metrics are used, the structural information of the textures is tried to be retained in contrast to the default metrics which tries to smooth the contents and potentially replace them by DC values. For severe compression, the decoded textures can have a noisy structure as HEVC intra prediction cannot provide anything better than parallel lines of the directional prediction. Using both metrics, wrong prediction directions might be selected. This is because these metrics are less sensitive to pixel by pixel comparison. LRI, as compared to STSIM, is much less computationally expensive. But it results in more wrong prediction directions than STSIM. The reason behind this is the approximation used when the metric is adapted to this work.

The encoder behavior is changed when these metrics are used. It always tends to use large prediction block sizes for all range of compression. This is mainly because in small blocks, there is a lack of enough statistics to compare and the metrics will return high dis-similarity values. Experimental results showed that using these metrics can improve the overall perceived quality of the decoded textures. The small details of textures are better preserved and the decoded textures look more pleasant.

In terms of texture similarity, the rate-distortion curves show that both metrics perform better than the default metric. The distortion metric that was used compares the energy and mean of subband frequency channels obtained by Gabor filter.

As a conclusion, the use of texture similarity metrics can generally give a better no-reference quality than the default one. A possible further improvement will be rate-distortion optimization and perceptual post processing to reduce the effects of wrong prediction direction. Beside that, since the metrics usually work for large block size, a hybrid mechanism can be used to reduce the complexity of their use. That is, for blocks larger than 16x16, the perceptual metrics can be used where the default metrics can be used for the rest.

Fig. 1 .

 1 Fig. 1. Texture dataset used in this paper (from USC-SIPI dataset [9])

Fig. 2 .

 2 Fig. 2. Example of using LRI and STSIM inside HEVC (QP=51). Top left: original Image, top right: encoded with default HEVC metrics, down left: encoded with LRI, down right: encoded with STSIM.

Fig. 3 .

 3 Fig. 3. Example of using LRI and STSIM inside HEVC (QP=43). Top left: original Image, top right: encoded with default HEVC metrics, down left: encoded with LRI, down right: encoded with STSIM.

Fig. 4 .

 4 Fig. 4. Example of using LRI and STSIM inside HEVC (QP=51). Top left: original Image, top right: encoded with default HEVC metrics, down left: encoded with LRI, down right: encoded with STSIM.

Fig. 5 .

 5 Fig. 5. Histograms of splitting depths vs QP. Each depth is scaled by the number of 4x4 block that it has

  his work was supported by the Marie Sktodowska-Curie under the PROVISION (PeRceptually Optimized VIdeo Com-presSION) project bearing Grant Number 608231 and Call Identifier: FP7-PEOPLE-2013-ITN.