Morphological analyses

Fig 1 illustrates the morphological diversity of Ramaria species. S1 Text describes the morphometric analyses.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1.

Morphological diversity of Ramaria species: (A) R. abetonensis AH 48006; (B) R. arcosuensis MA-Fungi 48049; (C) R. bataillei MA-Fungi 49453; (D) R. broomei MA-Fungi 79903; (E) R. cedretorum MA-Fungi 48074; (F) R. cokeri MA-Fungi 79893; (G) R. fennica var. fumigata AH 47767; (H) R. flavescens AH 48029; (I) R. ignicolor MA-Fungi 47978; (J) R. mediterranea MA-Fungi 39877; (K) R. praecox AH 47804; (L) R. rubrievanescens AH 47481; (M) R. stricta MA-Fungi 48068; (N) R. subbotrytis MA-Fungi 48088; (O) R. subdecurrens AH 48370; (P) Ampulliform septum of R. aff. capucina AH 48381; (Q) Cystidium on rhizomorphs from R. quercus-ilicis MA-Fungi 47984; (R) Bipyramidal crystal on tomentum from R. mediterranea MA-Fungi 39877; (S) Rosette crystals on rhizomorph hyphae from R. comitis MA-Fungi 47970; (T) Star-shaped crystals on rhizomorph hyphae from R. flaccida MA-Fungi 48020; (U) Spores of R. comitis MA-Fungi 47970; (V) Spores of R. cokeri MA-Fungi 79893; (W) Spores of R. praecox AH 47732; (X) Spores of R. botrytis MA-Fungi 48056. Scale bar = 5 μm.

https://doi.org/10.1371/journal.pone.0237507.g001

Sequence analyses

Our Total dataset contained 647 ITS sequences, with 231 sequences as “Ramaria sp.” S1 Table lists the PCI dataset used in the barcode analysis, which contained the remaining 416 sequences, with 382 unique sequences and 416–382 = 34 replicates. Our Long dataset was the subset of the PCI dataset whose sequences had length at least 500 bp. It contained 381 sequences, with 347 unique sequences and 34 replicates. Our Short dataset was the complement of the Long dataset within the PCI dataset. Our Taxonomic dataset was the subset of the PCI dataset whose sequences had expert taxonomic species identification. It consisted of 153 samples (107 samples that we identified, plus 46 GenBank sequences annotated as species types), with 151 unique sequences and 2 replicates. S2 Text contains “Sequence analyses” with summary statistics of the sequence length and nucleotide compositions of subsets of the PCI dataset. S1 Table lists the PCI dataset highlighting its Taxonomic and Short subsets

Barcode analysis

Call a sample or species a “singleton”, if the relevant dataset contains only a single DNA sequence from the species. The PCI dataset retained singletons, because they provide decoys that complicate the species identification of other samples, mimicking practical difficulties in barcoding [16]. To quantify barcoding accuracy, the present article defines four variants of the probability of correct identification (PCI), all of which are averages of species PCIs over non-singleton species. The Methods section describes the four PCIs, and the Theory section defines them mathematically.

The PCI dataset had 416 samples, distributed among 75 singletons and 71 non-singleton species. Of these 75 singletons, two were from the genus Gloeocantharellus; one, from the genus Gomphus; and one, from the genus Schildia. None of these singletons shared its sequence with another species. The PCI dataset also contained four samples from Gomphus clavatus; two, from Gomphus ludovicianus; and three, from Turbinellus floccosus, so these non-singleton species contributed to PCIs. Ramaria samples contributed all other Species PCIs for the PCI dataset.

Presented as mean ± one standard error mean, the PCIs averaged over the n = 71 non-singleton species in the PCI dataset were: the Barcode Gap PCI p_B = 0.324±0.056; the Unanimous PCI p_U = 0.451±0.059; the Average Unanimous PCI p_A = 0.691±0.040; and the Fractional PCI p_F = 0.731±0.038. On any fixed dataset like the PCI dataset, the PCI inequality 0≤p_B≤p_U≤p_A≤p_F≤1 always holds (as demonstrated by the mathematical “Proof of the PCI Inequality” in the Theory section). As the PCIs increase, they convey increasingly optimistic assessments of the same set of species identifications. Thus, the PCI inequality shows that the PCIs reflect progressively weaker criteria for correct species identification.

To assess annotations misidentifying species within the PCI dataset, we partitioned the PCI dataset into two non-overlapping sets: a Taxonomic dataset (153 samples, distributed among 74 singletons and 26 non-singleton species) and a GenBank dataset (263 samples). The GenBank dataset consisted of every GenBank sequence not annotated as having derived from a species typus. We then modified the log-rank test in survival analysis for taxonomic purposes (see details under “A Taxonomic Unanimity Test” in S2 Text). The data for the unanimity test are a series of 2x2 tables. In the context of survival analysis, each table represents a different time of observation, and the table for each time counts the patients in each of the control and treatment groups that survived or died since the previous observation. In contrast, the taxonomic unanimity test does not impose any temporal relationships among the 2x2 tables of the log-rank test. It retains, however, the formal mathematical manipulations of the log-rank test to produce a z-score with an approximately Gaussian distribution. In our taxonomic application, each table represents a different species, and it counts the samples in each of the Taxonomic and GenBank datasets that were correctly or incorrectly identified (i.e., the identifications agreed or disagreed with annotation). Applied to the 71 non-singleton species, the taxonomic unanimity test gave a continuity-corrected z-score Z = 1.055–0.260 = 0.795 (written as uncorrected z-score minus its continuity correction) with a one-sided p-value p = 0.213. If the p-value fell below a significance threshold, we could declare at that threshold that the Taxonomic dataset contained significantly better species identifications than the GenBank dataset. Within our PCI dataset, therefore, GenBank identifications were not significantly worse than expert taxonomy. The implicit species structure of the statistic suggests the explicit construction of S3 and S4 Tables in S2 Text, where barcoding results highlight possibly misannotated samples for further taxonomic scrutiny.

In comparing barcoding pipelines for species identification, the taxonomic unanimity test can also indicate when a pipeline significantly improves species identification. Short ITS sequences might have degraded species identification, so we considered removing them from our analysis pipeline. To assess the degradation, we partitioned the PCI dataset into a Long dataset (381 samples, distributed among 73 singletons and 67 non-singleton species) and a Short dataset (35 samples). The Short dataset contained every non-singleton sample sequence with length less than 500 bp. The resulting continuity-corrected z-score Z = 0.509–0.301 = 0.208 (again, written as uncorrected z-score minus its continuity correction) from the taxonomic unanimity test yielded a one-sided p-value p = 0.418. If the p-value fell below a significance threshold, we could declare at that threshold that the Long dataset contained significantly better species identifications than the Short dataset. The actual p-value shows that the estimated barcoding accuracy of species identification (using the fixed sampling density of the PCI dataset) did not differ significantly between the Short and Long datasets. Like S3 Table in S2 Text for the GenBank dataset, S4 Table in S2 Text lists the samples in the Short dataset, where barcoding results single out possibly misannotated samples for further taxonomic scrutiny.

Implicitly, the taxonomic unanimity test on the Short dataset above used the following bioinformatics pipeline A&E (“align and extract”): (1) perform a multiple sequence alignment of the PCI dataset; and then (2) extract the Long dataset while maintaining the original multiple sequence alignment. We had concerns, however, that pipeline A&E might be inferior and degrade species identification relative to another pipeline, E&A, which: (1) extracts the Long dataset first; and only then (2) aligns the extracted sequences. Accordingly, we performed a Wilcoxon matched-pair signed-rank test [22, p. 75–83] on the pairs of Species PCIs for each species from the two pipelines, A&E and E&A.

Table 1 shows the (Overall) PCIs under the headings Pipeline A&E and for Pipeline E&A (an Overall PCI is the average of Species PCIs over all species). The column headed by n counts the number of species where each pipeline produced a different Species PCI from the other (e.g., for the Fractional PCI p_F, the pair of Species PCIs p_F,s was different for 18 species). None of the PCIs showed any statistical significance, so reusing sequence alignments from the PCI dataset did not significantly reduce the accuracy of species identification within the Long dataset. In fact, PCIs decrease from pipeline A&E to pipeline E&A, suggesting that reusing sequence alignments from the PCI dataset may be a better bioinformatics strategy than realigning subsets.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Comparison of different pipelines on the LONG dataset.

https://doi.org/10.1371/journal.pone.0237507.t001

Just as with the Long dataset, pipelines A&E and E&A may both be applied to the Taxonomic dataset. S5 Table in S2 Text gives results for the Taxonomic dataset in the same format as Table 1. In fact, the PCIs there increase from pipeline A&E to pipeline E&A, so the reuse of sequence alignments from the PCI dataset is not a consistently superior strategy.

To quantify the effect of database size on our four PCIs, we resampled subsets of different sizes from the PCI dataset under uniform sampling without replacement (for details, see “Resampling the PCI dataset” in S2 Text). In Fig 2, the squares correspond to average PCIs in the resampled datasets of 50, 70, 90, …, 390 samples. Error bars display the sample standard error of each resampled PCI. They decrease to 0 as the resampled size approaches the size of the PCI dataset (416) because sampling without replacement becomes less random (e.g., at a size of 416, resampling without replacement simply selects the PCI dataset every time). Fig 2 complies with the PCI inequality 0≤p_B≤p_U≤p_A≤p_F≤1 and displays a partition of the PCIs into two groups: (1) the Fractional PCI p_F and the Average Unanimous PCI p_A; and (2) the Unanimous PCI p_U and the Barcode Gap PCI p_B. Compared to p_U and p_B, p_F and p_A remain relatively constant as the resampled dataset increases in size; moreover, their error bars exhibit progressively less random variation.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Resampling PCIs from the PCI dataset.

The X-axis shows the number of samples in the dataset after resampling; the Y-axis, the probability of correct identification (PCI), with the points indicating mean and sample standard deviation of the resampled distribution. The PCIs are the Barcode Gap PCI p_B (bottom, purple), the Unanimous PCI p_U (second from bottom, blue), the Average Unanimous PCI p_A (second from top, orange), and the Fractional PCI p_F (top, red).

https://doi.org/10.1371/journal.pone.0237507.g002

During analysis, the PCI dataset underwent a small revision (with S1 and S2 Tables giving the final datasets). After completion of the barcode analysis (with results given verbatim above), e.g., the species Ramaria rubribrunnescens provided a concrete example of the effects of misannotation on barcode analysis. All its six samples in the PCI dataset were from GenBank. We discovered, however, that another GenBank sample (KY354750) was in fact from Ramaria rubribrunnescens and had been misannotated as Ramaria rubribrunescens. Under barcode computations, therefore, it was a singleton sample from a non-existent species. KY354750 shared its sequence with JX310406 (also Ramaria rubribrunnescens), and together, their identical sequences provided the two nearest neighbors of EU652352 (also Ramaria rubribrunnescens). In fact, barcode species identification within Ramaria rubribrunnescens was perfect, i.e., the species PCIs for s = Ramaria rubribrunnescens were p_B,s = p_U,s = p_A,s = p_F,s = 1 (see “Mathematical Definitions of the Four PCIs” in the Theory section). The misannotated sample, however, had spuriously reduced the species PCIs to p_B,s = p_U,s = 0, p_A,s = 4/6≈0.667, and p_F,s = 4.5/6 = 0.750, with the species PCIs having the most stringent criteria for correct species identification, the Barcode PCI p_B,s and the Unanimous PCI p_U,s, suffering the greatest distortion. The PCI dataset had 71 non-singleton species, however so the misannotation perturbed overall PCIs by at most (1–0)/71≈0.014, well within statistical errors (see the PCIs at the start of “Barcode analysis” in the Results).

Supplementary parsimony and Bayesian analyses

S1 and S2 Tables list the datasets contributing to the phylogenetic analysis (with the correct species annotation for KY354750). In our analyses (Fig 3 and S1 Fig), including two Gloeocantharellus as outgroup, 73 clades and 74 singletons were obtained, of which 97 samples could be assigned to known taxa (3 Gomphus, 92 Ramaria, 1 Turbinellus, and 1 Schildia), and 50 might be undescribed species, but also possibly known species but misidentified, or related to known taxa (e.g. R. conjunctipes-clade 2, R. gracilis 1, R. thiersii 1, R. neoformosa 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. The 50% majority rule Bayesian tree inferred from ITS nrDNA assuming HSK+ I + G model of Gomphales included in the PCI analyses.

Colors: Blue, unique species; Green, one species and related infraspecific taxa; Red, clades including two or more species; Orange: species that occur in two or more clades. See S1 Table and S1 Fig for clades and singleton names. In Fig 3, the left circle corresponds to clade identifications before DNA barcoding; the right, after DNA barcoding. Essentially, therefore, the left circle displays morphological species; the right circle displays the clades after DNA barcoding has increased the apparent number of monophyletic species.

https://doi.org/10.1371/journal.pone.0237507.g003

As shown in Fig 3 and S1 Fig (the same tree with sample identifiers at its leaves), 31 clades and 50 singletons (in blue) represent single species; 5 clades (in green) included related infraspecific taxa; 10 clades (in red) included sequences with different species; and 27 clades and 24 singletons (in orange) belong to 23 taxa that occur in two or more clades.

Some examples of clades with only one species (S1 Fig, in blue) are Ramaria suecica, R. stuntzii, R. maculatipes or R. amyloidea; and the singletons R. cyanocephala, R. boreimaxima or R. quercus-ilicis. In the case of green clades, such as R. flavescens, the sequence obtained from the type of R. flavescens var. suaveolens is included; the other four green clades also include a species and its varieties. The most significant example of red clades is in the top of S1 Fig, and involve seven names (Ramaria abietina, R. abietina 1, R. apiculata, R. apiculata 1, R. apiculata var. brunnea, R. rubella, and R. tsugina). In the case of orange clades, a well-known species in temperate regions, such as R. abietina, appears not only in the red clades already mentioned, but also in the red clades named by us as R. flaccida and two orange clades, one of them the R. abietina sensu stricto. The right of Fig 3 displays how information from DNA barcoding improved tentative initial morphological identifications shown on the left of Fig 3 to increase clade consistency, particularly in the final clades belonging to a single species (in blue).

Of the 35 short sequences included in the analysis, 24 matched with expected clades (singletons excluded); and of 12 with ambiguous nucleotides, only one matched with expected clades (R. argentea JQ408231 from GenBank).

Barcode analysis

Typically in DNA barcoding, probabilities of correct identification (PCIs) provide the figure of merit for comparing variant bioinformatics pipelines [17]. The barcoding literature contains several types of PCI (e.g., [18]), and unjustifiably, some studies have even compared different types of PCIs directly, as though their status as probabilities permitted direct comparison. The PCI Inequality 0≤p_B≤p_U≤p_A≤p_F≤1, proved mathematically shows that for a fixed dataset, some types of PCIs are always larger than others. The type of PCI used is therefore critical in evaluating the results of a study.

PCIs similar to species PCIs can be calculated at other taxonomic levels like genus, family, etc. As figures of merit for taxonomic identification, the corresponding overall PCIs have the property of “taxonomic normalization”, described here only for the species level. In contrast, several studies [19, 23–25] used a measure related to PCI, namely, the fraction of sequences with correctly identified species. To demonstrate by exaggeration the defects of the measure, if 75% of samples come from a single species, the figure is thoroughly skewed by the dominant species. On the other hand, an overall PCI that averages species PCIs automatically normalizes datasets against overrepresented species [26]. Sometimes, normalization may be undesirable, e.g., if species have differing importance to the purpose at hand, an explicitly weighted average of species PCIs may be pertinent. Regardless, the ability of the overall PCI to normalize datasets taxonomically and reduce skewing by over-represented species led to its broad acceptance among botanists [16] and mycologists [15].

Taxonomic normalization has other desirable consequences. Measures without it require detailed ancillary data to support their relevance, e.g., rather than just a species count or enumeration, those measures require a separate demonstration that the sequences are about evenly distributed across the relevant species. When comparing different datasets, a figure of merit requiring detailed ancillary data can seriously impede human evaluation (e.g., Section 3.2 in [27] criticizes the ROC_n in sequence retrieval because the ROC_n requires ancillary information to verify its validity). To avoid distracting humans with detailed ancillary data, the figure of merit then requires a rigid dependence on benchmarking with standard datasets to ensure scientific progress.

As indicated in the Introduction and Fig 2, a figure of merit should remain stable as a database grows. If adventitious database features like species overrepresentation can profoundly alter a figure of merit, the figure is not fulfilling its purpose.

For any specific purpose, the four types of PCI may have widely differing utilities. The Barcode Gap PCI (p_B) reflects the most stringent criterion for correct species identification, namely, the presence or absence of a barcode gap for a species. Because of its stringency, and because DNA barcodes might be used to construct phylogenetic trees as well as to identify species, botanists [16] and mycologists [15] used the Barcode Gap PCI to select their barcode markers. Table 1 indicates that compared to other types of PCIs, the Barcode Gap PCI is less broadly informative of species behavior in a bioinformatics pipeline, however. The column labeled n in Table 1 shows that only three species had different Species Barcode Gap PCIs for pipelines A&E and E&A. As criteria for correct species identification became less stringent (i.e., as PCIs became larger), n also increases, reaching a maximum value n = 18 for the Fractional PCI. S5 Table in S2 Text of the SI displays a similar trend.

The Barcode Gap PCI displays impressive improvements in Table 1 and S5 Table in S2 Text of the SI, but a single misannotated sample can spuriously destroy the barcode gap in a species s. The Results for Ramaria rubribrunnescens show that a single sample can change the PCI for species s violently, from p_B,s = 1 to p_B,s = 0, so the Barcode Gap PCI displays an undesirable sensitivity to errors in the data. The Unanimous PCI (p_U) suffers a similar defect, because all samples in a species s are either correctly identified or not (i.e., either p_U,s = 1 or p_U,s = 0).

Insensitivity to subtle changes across several species and vulnerability to misannotation or otherwise anomalous taxonomy probably make the Barcode Gap and Unanimous PCIs among the worst PCIs for assessing gradual improvements in a barcoding pipeline. After the proof of the PCI inequality, the Theory section contains a second mathematical proof: if the Unanimous PCI p_U,s can detect that procedural pipeline has improved identification within a species s, the Average Unanimous PCI p_A,s can also. The statement virtually disqualifies the Unanimous PCI as a desirable figure of merit when improving barcoding pipelines.

A violent discontinuity in a PCI has other undesirable consequences. Fig 2 shows results from resampling the PCI dataset to produce smaller datasets. Resampling (e.g., the bootstrap [28]) is a popular method to make inferences about complex situations, although its inferences about maxima and minima carry some caveats (see particularly, counter-example 2 in [29]). When querying a barcode database with a new sample, many practical computer algorithms for species identification depend on the new sample’s nearest neighbors, database sequences that minimize some distance to the sample sequence. This article measured genetic variation between two aligned DNA sequences with p-distance [30], the fraction of unambiguous nucleotide pairs that were mismatches. Caveats notwithstanding, Fig 2 indicates that the present values of p_B and p_U are less predictive than p_A and p_F of their future values as a barcoding database grows. As new samples enter a barcoding database, a single sample in a species s can reduce the Barcode Gap or Unanimous PCI for species s (p_B,s or p_U,s) violently, from 1 to 0, whereas its effects on the other PCIs are more gradual (e.g., see again the Results for Ramaria rubribrunnescens).

When evaluating a pipeline for taxonomic identification, a singleton species may cause difficulties if only one sequence dataset must serve both as a reference database and a source of queries. Following the pattern of previous studies [15, 16], our PCIs considered only non-singleton species and omitted each query in turn from the reference database. Each query then had a species representative in the reference database available as a possible nearest neighbor. Consequently, the evaluation of the pipeline conferred a special status on singletons, not our identification algorithm. At first, the distinction between evaluating a pipeline and the pipeline itself can be subtle, and the reader should reread this paragraph if its point is unclear. The early paragraph in the Introduction, starting “To avoid confusion…” is also relevant.

When the evaluation of a pipeline uses a single dataset as a reference database and a source of queries, PCIs should be accompanied by a separate count of species singletons having unique sequences. There is nothing objectionable in counting the singletons separately, but calculating a single combined PCI from non-singleton species and the singletons [25] has the undesirable consequences mentioned during the discussion of resampling: violent change and lack of predictability as a barcoding database grows. In particular, a single PCI combined with singletons typically views a singleton as correctly identified if its sequence is unique within the database. Another sample from the same species entering the reference database then lowers the species PCI from 1 to 0, and no third sample from its species modulates the change.

On one hand, although we examined the genus Ramaria as a case study, we expect our observations to be most pertinent to other taxa with relatively low PCIs. On the other hand, if taxonomic identification is straightforward, then methodology becomes less important.

The present article also introduced a taxonomic unanimity test. The tests and tables of samples derived from the test (e.g., S3 and S4 Tables in S2 Text of the SI) may help when cleaning taxonomy from and reporting errors to a database like GenBank or allaying fears about using sequences with ambiguous nucleotides. Similarly, the signed-rank test applied species by species showed that expedients like reusing sequence alignments did not significantly worsen results. With other datasets, it might have shown that improvements to a barcoding pipeline are statistically significant. We chose to use it, because it has greater statistical power than its competitor, the sign test, making it an orthodox choice. In a taxonomic context, however, the signed-rank test gains its power by using ranks to break the symmetry of the species under scrutiny. Studies rarely sample species evenly, however, so sampling biases can subtly and systematically influence the ranking of species, and therefore the ranked-sign test. Sampling biases do not affect the sign test, which maintains symmetry between species throughout the analysis. In the taxonomic context, therefore, the sign test has some features to recommend it, but its use would not have changed the scientific conclusions in this article.

In addition, although not performed here because of insufficiently deep sampling of two closely related species, a Fisher exact test for individual 2x2 tables (two species, two experts) could flag statistically significant expert disagreement on species boundaries, prompting closer examination of systematic expert disagreement on characters defining the species. The basic motivation for all statistical tests suggested here is that non-parametric tests [31] are readily applied to taxonomic groups of equal weight (here, each species PCI contributes equally to the p-value).

Statistical findings were uniformly negative within Ramaria. Within the PCI dataset, barcoding species identification was not significantly worsened by: (1) non-expert annotations in GenBank; or (2) short sequences (less than 500 bp). For both the Taxonomic and Long datasets, the pipeline reusing the sequence alignment from the PCI dataset did not have significantly worse species identification than the pipeline that extracted the datasets and then aligned them. The negative statistical findings therefore allayed our initial concerns that the PCI dataset required further filtering before barcoding analysis. In addition, they confirmed that the difficulties that Ramaria present to taxonomists extend to the nucleotide realm.

Parsimony and Bayesian analyses

About a third part of the accepted Ramaria species (96/300) were included in our analyses, and 50 clades could be related to cryptic, misidentified or unknown species. As in other groups, some of them are known species, but they do not appear with the correct name in GenBank: in 2003, about 20% of the named sequences were misidentified [7, 32, 33]; more recently, nearly 30% [14]. In the opposite case, some sequences could appear misidentified, because they belong to different ITS copies of the same species; this was one the criticisms of using ITS as a universal barcode for fungi [34]. If species annotations are correct, ITS sequences from Ramaria conjunctipes, R. botrytis, R. formosa or R. stricta display a large molecular variation; and their clades could be species complexes or contain cryptic species.

Ramaria and related genera exemplify the importance of sound taxonomic knowledge before starting molecular analyses in a group. As a case study, they also show how molecular analysis helps taxonomists to improve the identification of doubtful samples, and how “the taxonomic feedback loop” links both disciplines synergistically [35]. During the last years, phylogenetic analyses using ITS sequences have helped us to elucidate many groups of fungi (e.g. gasteroid and corticioid) and to describe new species (e.g. [36–38]). Unfortunately, however, even if molecular analysis detects a new species, some authors do not follow through by describing and naming it. Moreover, they also avoid the responsibility of curating their sequences in DNA repositories as well as in publications. These repositories contain a many unidentified sequences, even though many of them are connected to vouchers and are not environmental samples. In our study, 647 sequences from Gomphales, mainly Ramaria, were downloaded, but 231 had a non-specific epithet (e.g., “sp.”), and our PCI analysis had to exclude them. Although not shown in this paper, many of these Ramaria sp. could be assigned to a known species, however. On one hand, GenBank places sample KT824242 under Ramaria sp. (S2 Table), e.g., but its specimen voucher (KD-14-006) is the holotype of R. subalpina [39]. On the other hand, sample MH216040 (S2 Table, holotype) was under Ramaria sp. PA-2018 when we did the analyses, but was later renamed as R. parabotrytis [40].

Many species names appear to be added after a direct BLAST search, identified only by macromorphology, without microscopic examination, e.g., R. apiculata clade contains many sequences under R. abietina. Microscopical features differentiate the two species: R. apiculata has rosette oxalate crystals in the mycelium, ellipsoid and verrucose spores, and ornamented ampulliform septa, whereas R. abietina has stellate oxalate crystals in the mycelium, amygdaloid and nearly spinose spores, and smooth ampulliform septa (Fig 4). This is a good example showing that DNA repositories need to include high quality sequences from well identified specimens [41, 42].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Comparative of features of Ramaria abietina and R. apiculata.

Ramaria abietina: (A) Basidioma (AH 48373); (B) Spores (MA-Fungi 49419); (C) Crystals from mycelium (MA-Fungi 49419); (D) Ampulliform septum (MA-Fungi 49419). Ramaria apiculata: (E) Basidioma (AH 47751); (F) Spores (MA-Fungi 48462); (G) Crystal from mycelium (MA-Fungi 48485); (H) Ampulliform septum (MA-Fungi 47981). Scale bars = 10 μm.

https://doi.org/10.1371/journal.pone.0237507.g004

In addition, many mistakes in the DNA repositories are related with contaminations or mixed-up samples. Among the sequences obtained for this study, after the phylogenetic analyses we detected a contamination; our sequence AF442098 identified as R. praecox grouped in the R. flava clades. However, R, praecox is very different to R. flava. In R. praecox spores show conspicuous ornamentation, hyphae without clamps, and vernal phenology, whereas in R. flava, spores are slightly ornamented, hyphae clamped, and mainly autumnal phenology (Fig 5). Just as in metagenomic sequencing [43], contamination can arise from previous amplifications. Here, although we also include negative controls to each set of amplifications as a routine, contamination occurred.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5.

Comparative of features of Ramaria flava and R. praecox. Ramaria flava: (A) Basidioma (AH 48375); (B) Spores (MA-Fungi 48465); (C) Clamped hyphae from context (MA-Fungi 48484); (D) Hyphae from mycelium (MA-Fungi 48484). Ramaria praecox: (E) Basidioma (AH 47804); (F) Spores (AH 47732); (G) Unclamped hyphae from context (AH 47732); (H) Hyphae from mycelium with crystals (AH 47804) or mucilaginose drops (AH 47732). Scale bars = 10 μm.

https://doi.org/10.1371/journal.pone.0237507.g005

A preliminary analysis detected possible mix-ups in the R. aff. spinulosa sequence and one sequence of R. botrytis, so we amplified and sequenced them again. The sequences were then correct: AF442904 corresponded to R. aff. spinulosa; and AJ292294, to R. botrytis. Only taxonomic characters of the species can expose such errors, so the metadata in DNA repositories should include the voucher for sequences [44]. (R. botrytis has grooved spores from a fleshy and white basidioma with reddish apices whereas R. aff. spinulosa has verrucose spores from a delicate brownish basidioma with concolor apices).

Moreover, we must remember questions related to nomenclature, e.g., sequences in DNA repositories that belong to one species, but appearing under two or more different names. During our morphological study we have unified some of these names. Thus, for example, R. luteovernalis is synonymous with R. praecox; or R. fennica var. violacea synonymous with R. cedretorum (S1 Table). In other cases (S1 Fig and S1 Table), some nomenclature questions to be addressed after the phylogenetic analyses included in this study, e.g., R. claviramulata and R. celerivirescens published as new species in [45], where our analyses grouped them into the same clades. Exeter et al. [46] suggested that R. celerivirescens is the correct name to the species. We studied the unique collection from R. claviramulata (holotype), and the basidiomata show a morphology that could correspond with a phytoplasma infection, as we have shown in other Ramaria species [47].

Mathematical definitions of the four PCIs

Fix a species s, with n samples in s. For each sample σ = 1, 2,…,n from s, define the nearest neighbors of σ to have the minimum p-distance from σ (in the relevant dataset). Let the fraction of nearest neighbors of σ within s be f_σ(0 ≤ f_σ ≤ 1). Define the unanimity function U(f_σ) = 1 if f_σ = 1 (all nearest neighbors of σ are in s); and 0, otherwise. For species s, define the Species Fractional PCI p_F,s = (f₁ + f₂ +…+ f_n)/n; the Species Average Unanimous PCI p_A,s = (U(f₁)+U(f₂)+…+U(f_n))/n; and the Species Unanimous PCI p_U,s = min{U(f₁), U(f₂),…,U(f_n)}. The Species Barcode Gap PCI p_B,s = 1 if the maximum intraspecies p-distance is strictly less than the minimum interspecies distance (i.e., species s has a barcode gap); and 0, otherwise. Note for any species s: (1) p_U,s and p_B,s only take the values 0 and 1; and (2) all four of the Species PCIs p_s ∈ (p_F,s, p_A,s, p_U,s, p_B,s) satisfy 0 ≤ p_s ≤ 1. For any of the four Species PCIs p_s, define a corresponding (Overall) PCI p by averaging the Species PCI over all species. In a standard notation from physics, we denote the average over species with angle brackets, e.g., p_B = 〈p_B,s〉_s. S1 Table contains several examples of these calculations for samples in species in the PCI dataset.

Proof of the PCI inequality

(1)

By definition, 0≤p_B,s. To prove p_B,s ≤ p_U,s, there are two cases: p_B,s = 0 or p_B,s = 1. If p_B,s = 0, then automatically p_B,s ≤ p_U,s. Otherwise, p_B,s = 1, and the barcode gap for species s implies that for every sample σ in the species s, all nearest neighbors belong to s, i.e., f_σ = 1 and U(f_σ) = 1. Thus, p_U,s = min{U(f₁), U(f₂), …,U(f_n)} = 1, so again p_B,s ≤ p_U,s. To prove p_U,s ≤ p_A,s, the minimum p_U,s of the numbers U(f_σ) can never exceed their average p_A,s. To prove p_A,s ≤ p_F,s, note that U(f_σ) ≤ f_σ (on one hand, if f_σ = 1, then U(f_σ) = f_σ = 1; on the other hand, if f_σ < 1, then 0 = U(f_σ) ≤ f_σ). Average over all samples in the species s to derive (2)

Finally, p_F,s ≤ 1 to finish the proof of Eq (1). In Eq (1), average over all species s to derive (3)

The definitions p_B = 〈p_B,s〉_s, etc., yield (4) finishing the proof.

The Average Unanimous PCI p_A,s is more sensitive than the Unanimous PCI p_U,s for detecting improvements in a procedural pipeline

Consider a procedural pipeline (represented by unprimed symbols), later improved to another pipeline (represented by primed symbols). Then for any species s, if , then .

Proof: If , then (because p_U,s is either 0 or 1). After improvement on one hand, . Thus, and . Originally on the other hand, p_U,s = min{U(f₁), U(f₂), …,U(f_n)} = 0, so at least one i = 1, 2, …, n satisfied U(f_i) = 0. Consequently, p_A,s = (U(f₁) + U(f₂)+…+U(f_n))/n<1. Because , we have proved .

A taxonomic unanimity test

A taxonomic unanimity test can detect improvements in barcode species assignment, as follows.

Consider a specific barcoding sequence dataset, in which every sequence bears a species identification from one of two sources. It might be useful to have objective evidence that one source (e.g., an expert taxonomist) is systematically more accurate than the other. In a related situation, one might be concerned that species identification is superior in one of two disjoint sets of samples, e.g., that sequences longer than 500 bp identify species better than sequences shorter than 501 bp. To be specific, consider the case of expert identification.

Consider a sequence dataset of interest and any species s with at least two samples in the dataset. The nearest neighbors of a sample from s are the sequences in the dataset, excluding the sample sequence itself, with the minimum p-distance to the sample sequence. We say that a sample from species s has unanimous neighbors if its nearest neighbors are all from s.

Let a dataset of N (more than 2) sequences be assembled from two sources, e.g., with and without expert species taxonomy. Construct a 2x2 table corresponding to the species s (see “A Taxonomic Unanimity Test” in S2 Text), as follows. Of the sequences with expert taxonomy, let A count the sequences with unanimous neighbors; C, the sequences without. Similarly, of the sequences without expert taxonomy, let B count the sequences with unanimous neighbors; D, the sequences without. Thus, of the N = A+B+C+D samples, A+C had expert species taxonomy, whereas B+D had not. Similarly, each sample within species s either had unanimous neighbors or had not, and of the N samples, A+B were unanimous, whereas C+D were not. One might hope that expert taxonomy displays itself by systematically increasing the number of unanimous neighbors.

If the unanimity of neighbors is probabilistically independent of the source of taxonomic identification, then in species s, of the A+C samples with expert identification, the expected count of samples with unanimous neighbors is E_s = (A + C)(A + B)/N, because (A + B)/N is fraction of samples with unanimous neighbors. The observed count of samples with expert taxonomy and unanimous neighbors is in fact O_s = A. Conditioned on the marginals A+C, B+D, A+B, and C+D, if the source of taxonomic identification is independent of the samples of species s, A follows a hypergeometric distribution, so its variance is (5)

Sum the excess O_s—E_s over all species s and denote the sum by X. Sum similarly, and denote the result by σ². Under the reasonable approximation that species are probabilistically independent, X has an approximate normal distribution with variance σ². Thus, the upper tail of a standard normal cumulative distribution function beyond Z = X/σ gives a one-sided p-value for testing whether expert taxonomy is significantly more accurate than inexpert species identification.

Sequence datasets

The Total dataset consisted of 647 sequences including those under Ramaria sp. in the GenBank. Moreover, Gomphus, Schildia, and Turbinellus sequences were included due to the relationship with Ramaria mentioned in the introduction. Two samples of Gloeocantharellus were chosen as outgroup to the posterior phylogenetic analyses. Before proceeding to the phylogenetic analyses, searches with BLAST+ 2.8.1 checked dubious sequences and genera misidentifications [48, 49]. The PCI dataset consisted of 416 sequences (specifically excluding Ramaria sp. in the Total dataset). S2 Text contains summary statistics describing the lengths and nucleotide compositions of the sequence data. To assess the effects of short ITS sequences on species identification, we subdivided the PCI dataset into two datasets: (1) the Long dataset, sequences exceeding 500 bp; and (2) the Short dataset containing the remaining sequences. To assess the effect of faulty taxonomy in GenBank sequences on species identification, we subdivided the PCI dataset into two datasets: (1) the Taxonomic dataset, our samples plus species type-derived sequences in GenBank; and (2) the GenBank dataset containing the remaining GenBank sequences.

Multiple sequence alignment and p-distance

To prepare for multiple sequence alignment, we removed replicate sequences, so all aligned sequences were unique. Most multiple sequence alignment programs are progressive, so they add sequences sequentially and then align them. Progressive sequence alignment programs can therefore align replicate sequences differently, spuriously inferring a positive phylogenetic distance between identical sequences. In general, the removal of any replicate sequence before sequence alignment improves the stability of multiple sequence alignments [50], and surprisingly, the workflows in some multiple sequence alignment programs do not remove replicate sequences from their inputs before sequence alignment. Muscle 3.8.31 under its default parameters [51, 52] aligned our unique sequences in the Total, PCI, Long, and Taxonomic datasets.

The principle basic to DNA barcoding is that the genetic variation between species usually exceeds the variation within a species (e.g. [6, 53]). Many distances can measure genetic variation between sequences, but the choice of distance rarely influences species identifications (e.g., [13, 54–56]). For simplicity, the p-distance between two aligned DNA sequences [30], the fraction of unambiguous nucleotide pairs that were mismatches, measured genetic variation here.

Four types of the probability of correct identification (PCI)

Many variants of PCI appear in the taxonomic literature (e.g., [18, 57]). Given a dataset, we calculated the following four types of Species PCIs for all non-singleton species s. (The Theory section contains a mathematically formal “Mathematical Definitions of the Four PCIs”).

The Species Barcode Gap PCI (p_B,s) is 1 if species s has a barcode gap (i.e., its minimum interspecies p-distance exceeds its maximum intraspecies p-distance), and 0 otherwise.

Fix a dataset, and within species s, consider a sample’s nearest neighbors under p-distance. Call a sample unanimous (i.e., unanimously correctly identified) if its nearest neighbors under p-distance are all from species s. Each unanimous sample has a Sample Unanimous PCI of 1; all other samples have a Sample Unanimous PCI of 0. Define the Species Unanimous PCI (p_U,s) as 1 if every sample in species s is unanimous, and 0 otherwise.

Averaging provides another procedure for turning Sample Unanimous PCIs into Species PCIs. Define the Species Average Unanimous PCI (p_A,s) as the average of the Sample Unanimous PCIs over all samples in the species s.

Finally, let the Sample Fractional PCI to be the fraction of nearest neighbors belonging to species s. The Species Fractional PCI (p_F,s) is the average of the Sample Fractional PCIs over all samples in species s.

Each of the four Species PCIs (p_B,s, p_U,s, p_A,s, p_F,s) may be averaged over all (non-singleton) species s to produce the corresponding (Overall) PCI (p_B, p_U, p_A, p_F).

Consider any PCI obtained by averaging: (e.g., the (Overall) PCI is an average of Species PCIs; the Species Average Unanimous PCI is an average of Sample Unanimous PCIs; etc.). If the PCIs p₁, p₂,…, p_N are uncorrelated, the unbiased sample variance (6) and the standard error mean estimates the error in . Thus, the error in the (Overall) PCIs (p_B, p_U, p_A, p_F) can be estimated from the corresponding Species PCIs (p_B,s, p_U,s, p_A,s, p_F,s).

Multiple testing

All p-values are given without multiple-test correction. As only unpublished preliminary datasets had uncorrected p-values p ≤ 0.05, the planned (Bonferroni) multiple-test correction for our published datasets is irrelevant.

Supplementary parsimony and Bayesian analyses

To evaluate whether the species were recovered as monophyletic groups, the sequence alignment of the PCI dataset was analyzed under maximum parsimony (MP) using the heuristic search option in PAUPver4.0a [58], with a default setting for stopping the analysis. Phylogenetic trees were rooted with two sequences under Gloeocantharellus persicinus and Gloeocantharellus purpurascens. Gaps were treated as missing data [13]. Branch lengths equal to zero were collapsed to polytomies. Nonparametric bootstrap support [59] for each clades was assessed with the fast-step option, using 10,000 replicates, yielding a composite consistency index [60], retention index [61], and homoplasy index [61]. PAUP*Version 4.a selected as the best model for each partition ITS1, 5.8S, and ITS2 of sequences, as well as for their union ITS1+5.8+ITS2, so we had MrBayes 3.2 [62] do a Bayesian analysis [63, 64] with the model HSK+I+G. Two independent and simultaneous analyses starting from different random trees were run for two million generations with 12 parallel chains, and trees model scores saved every 100^th generation. The initial 1000 trees were discarded as burn-in before calculating the 50% majority-rule consensus tree and the posterior probability (PP) of the nodes, as described elsewhere [65]. A combination of bootstrap proportions and posterior probabilities was used to assess the level of confidence for a specific node [66]. The phylogenetic trees were viewed with FigTree v. 1.3.1 (http://tree.bio.ed.ac.uk/sotware/figtree/) and edited with Adobe Illustrator CS3 v.11.0.2 (Adobe Systems).

Perl Programs for calculating the PCIs are publicly available as open source software at https://go.usa.gov/xyuut.