Linking morphological and molecular sources to disentangle the case of Xylodon australis

The use of different sources of evidence has been recommended in order to conduct species delimitation analyses to solve taxonomic issues. In this study, we use a maximum likelihood framework to combine morphological and molecular traits to study the case of Xylodon australis (Hymenochaetales, Basidiomycota) using the locate.yeti function from the phytools R package. Xylodon australis has been considered a single species distributed across Australia, New Zealand and Patagonia. Multi-locus phylogenetic analyses were conducted to unmask the actual diversity under X. australis as well as the kinship relations respect their relatives. To assess the taxonomic position of each clade, locate.yeti function was used to locate in a molecular phylogeny the X. australis type material for which no molecular data was available using morphological continuous traits. Two different species were distinguished under the X. australis name, one from Australia–New Zealand and other from Patagonia. In addition, a close relationship with Xylodon lenis, a species from the South East of Asia, was confirmed for the Patagonian clade. We discuss the implications of our results for the biogeographical history of this genus and we evaluate the potential of this method to be used with historical collections for which molecular data is not available.


Results
Phylogenetic analyses. A total of 66 new sequences were generated in this study: 35 sequences for ITS nrDNA region and 31 for nrLSU. Final alignments including sequences from EMBL/GenBank/DDBJ databases contained 110 ITS nrDNA sequences for a dataset length of 724 characters and 87 nrLSU sequences with 987 characters. All new sequences have been deposited in the EMBL/GenBank/DDBJ database and their accession numbers are presented in Table S1.
The identity of all samples named under X. australis was confirmed by molecular data (see Supplementary Figs. S1-S3 on line) and their phylogenetic position among other Xylodon species was according to Riebesehl et al. 17 . The number of constant, variable, uninformative and informative characters under maximum parsimony analyses are indicated in the legend of the Supplementary Figs. S1-S3.
Statistical tests of morphological characters. ANOVA on basidia and spore morphology was conducted on eight Australian and 12 Patagonian specimens (Table 1).
Significant differences were found for all measures: basidia length and width, and spore length, width and length/width relation ( Table 2; Fig. 1). Specimens from Australia showed longer and wider basidia than samples from Patagonia. In the same way, spores of specimens from Australia were longer and wider than those of Patagonian samples. Spore length/width ratio for the Australian lineage was lower than for the Patagonian clade, thus the former has spores narrowly ellipsoid or subcylindric compared to the latter.
Inferring the position of the nomenclatural type of Xylodon australis in the molecular phylogenetic tree. The ultrametric phylogenetic trees of 47 Xylodon specimens gave the same topology for the ITS nrDNA (not shown), nrLSU (not shown), and ITS nrDNA + nrLSU datasets (Fig. 2). Effective sample sizes for all parameters were higher than 200. Bayesian inference analyses showed that specimens under the Xylodon australis name were distributed in two non-directly-related and highly supported clades. All Australian collections were grouped in one clade, while the other clade included all the Chilean and Argentinean specimens. Three sequences of species Xylodon lenis (including from the type specimen) are the sister clade of the Chilean and Argentinean specimens; this relationship has strong support in all datasets (ITS nrDNA PP = 1.0, nrLSU PP = 0.99, and ITS nrDNA + nrLSU PP = 1). Xylodon lenis, X. australis from Australia and Xylodon gr. australis from Patagonia formed the crown clade for all X. australis specimens (ITS nrDNA + nrLSU PP = 0.97).
The topology of the combined dataset for this clade was used as subtree to infer the phylogenetic position of Xylodon australis type specimen using continuous morphological traits (Fig. 3a), as well as to the New Zealand samples from which no sequences were obtained (see Supplementary Fig. S4).
The type specimen of Xylodon australis was located in the Australian molecular lineage, according to the five continuous morphological traits used, under a maximum likelihood framework using the locate.yeti (Fig. 3a). www.nature.com/scientificreports/ Results of the randomization test to assess the accuracy of the method for our case study are shown in Fig. 3b,c. In general, trees reconstructed according to continuous mophological traits, that is, using locate.yeti, were more similar to the actual molecular tree than those reconstructed by random tip location (lower values of branch score and quadratic paths distances, Fig. 3b,c). The likelihood ratio test was conducted by constraining the type position to the Patagonian clade (the alternative position to our results). The hypothesis that the type collection of Xylodon australis belongs in the Patagonian clade was rejected by our analyses (P-value < 0.01; Table 3).

Discussion
The exhaustive study carried out by Gresbelin et al. 3 was not enough to consider Xylodon australis and Xylodon magallanesii as two different species due to the lack of evidence in addition to morphology. However, these authors pointed toward a speciation process due to the differences found in spore size and shape of samples from each area. Our phylogenetic analyses not only confirmed the identity of two species under the X. australis name, but also revealed the relation of Xylodon lenis as the sister species of X. magallanesii (Fig. 2).
Our microscopic studies agree with Gresbelin et al. 3 , and showed statistically significant differences in basidia and spore size and shape between Xylodon australis and Xylodon magallanesii (Fig. 1). Spores of X. magallanesii were in general smaller, as were their basidia. This correlation between spore and basidia size has been traditionally reported 18,19 and, therefore, our results could be expected. However, spore shape can be related to more complex responses, including dispersal abilities, bioclimatic fitness, or life history characteristics 11,20 , so a specific study should be carried out to explain differences in basidiospore shape between X. australis and X. magallanesii. No molecular data were obtainable from New Zealand samples identified under the X. australis name, probably due to problems in material conservation (M. Padamsee, PDD Fungarium Curator, pers. comm.); for this reason, these specimens were not included in our statistical analysis in order to compare only those clades confirmed by molecular data. Taking into account our results from X. magallanesii, and despite the fact that our ML analyses using the locate.yeti function included all New Zealand samples in the Australian clade, future analyses should be carried out in order to confirm the identity of the X. australis samples from New Zealand. In contrast to the type specimen of X. australis, samples from PDD addressed in this study are available for destructive sampling, and therefore new approaches for DNA extraction could be successful.
The inferred position by locate.yeti function for the type of Xylodon australis, in the molecular tree arranged with Australian samples, supported the designation of Xylodon magallanesii samples from Chile and Argentina as the new species (Fig. 3a). Our analysis resulted in a high performance in the validation test, and was able to correctly locate the pruned molecular taxon better than randomness for our study group (Fig. 3b,c). These results may be due to choosing morphological traits to infer the phylogenetic position, which are known to be taxonomically informative in differentiating closely related species in Xylodon 21 . Since the performance of the www.nature.com/scientificreports/ method is known to increase with the amount of data 11 , the quality of morphological traits could be a key factor when the number of characters is small or the number of taxa in the molecular tree is less than 20, as in our case. In addition, the likelihood ratio test rejected the alternative position proposed for the Xylodon australis type in the X. magallanesii clade. Therefore, its position with the Australian samples is strongly supported. These results are also in accord with geographic evidence, since the type specimen of X. australis was collected in Tasmania and therefore their connection with specimens from Australia was expected. The close phylogenetic relation between Xylodon australis, X. lenis and X. magallanesii is reported in our study for the first time. Xylodon lenis was described as Hyphodontia mollis by Wu, in 1990 from Taiwan 22 . Though the X. australis characteristic color change with the application of KOH was not described for X. lenis, Wu 22 highlighted the presence of granular material over the hyphal system that dissolves in KOH, a character also shown by X. australis and X. magallanesii. In an additional macromorphological inspection conducted on an X. lenis isotype  for this study, we could observe color change from orange toward violet after the application of KOH. Therefore, as in other fungal groups 23,24 , this character emerges as a useful trait for taxonomic classification when closely related species are compared. Another morphological character that points to a relation among these three species is the cracked hymenial surface, described in several studies 3,22 , and also shown by X. magallanesii.
From a biogeographical point of view, these phylogenetic relations are a challenge due to the distribution pattern of the species. The spatial structure inside each species remains congruent and a geographic isolation between them is shown in our results (Fig. 2). The sister relation between Xylodon magallanesii from Patagonia and Xylodon lenis from continental China and Taiwan could be explained as an example of long distance www.nature.com/scientificreports/ dispersion 25 . This ability has been confirmed for many fungi, even for the same South America/Asia/Australia pattern in some cases, such as for the Ganoderma applanatum-australe species complex 26 . Other possible explanations for the disjunct distribution of sister species may be due to incomplete taxon sampling, or to the extinction of lineages that had linked these species in the past. Epitypification and neotypification have been proposed as possible solutions to address taxonomic confusion in those cases where: type specimens are damaged, characters used for species identification are not manifest, or the material is not available 27 . However, these solutions should be applied cautiously, since they have many other associated risks 27 . Some other authors 10 have argued that, although the best possible solution would be to examine the type material, one option could be to describe new and well-documented taxa and to ignore old species names, but other problems such as taxonomic inflation, could arise with this practice 28 .
Our study shows how new methodological frameworks can help to solve old taxonomic problems that have become more evident during the DNA era. The possibility to place a type collection into a molecular tree, using phenotypic traits, increases the value of herbaria and museum collections. This is especially important in groups such as Xylodon, in which new species and combinations are being proposed every year 19 , and taxonomy is www.nature.com/scientificreports/ quickly changing 29 . The study of type materials is essential to avoid bad taxonomy that could lead to important ecological and economic losses 30 .

Material and methods
Taxon sampling and morphological studies. A total of 37 specimens of Xylodon australis from five different herbaria and private collections were analyzed in this study: Australian National Herbarium (CANB), New Zealand Fungal & Plant Disease Collection (PDD), Real Jardín Botánico from Madrid (MA-Fungi), and Alina Greslebin and Mario Rajchenberg private collections (Table 1). In addition, the type specimen of X. australis (under Grandinia australis Berk.) from Royal Botanic Gardens Kew K(M) was also studied morphologically. Basidioma colors were recorded according to Kelly and Judd 31 . Color changes were examined with 3% aqueous KOH. Microscopic measurements were made from sections mounted in aqueous solutions of 3% KOH and 1% aqueous solution of ammoniacal Congo red or 1% aqueous floxine. Sections were examined at magnifications up to 1250× using an Olympus BX51 microscope. Six basidia were measured from each sample. The width (W) and length (L) of 10 spores were also measured and length/width ratios (Q) were calculated. Average values of each character were calculated for each specimen. Additional morphological measurements were performed to provide a general description of each species. Drawings were made with the aid of a drawing tube.
DNA extraction, amplification and sequencing. Genomic DNA isolation was performed using DNeasy Plant Mini Kit (Qiagen, Valencia, California, USA) following the manufacturer's instructions, except in three steps: the incubation with the RNAase was done overnight at 65 °C, a second drying at 20,000×g was done for 2 min after cleaning with AW buffer, and elution buffer was preheated to 60 °C. When this extraction was not successful, FTA Indicating Micro Cards (Cat Nº WB120211, Whatman, Maidstone, England) were used following the protocol in Telleria et al. 32 . DNA amplifications, purifications and sequencing protocols are deposited in protocols.io (https ://doi.org/10.7504/proto cols.io.wpdfd i6). Polymerase chain reactions (PCR) were performed to amplify DNA from two loci using the following primer combinations: ITS5/ITS4 33 were used to obtain DNA amplifications of the nuclear ribosomal internal transcribed spacer regions ITS1 and ITS2, including 5.8S, ITS nrDNA barcode 6 and LR0R/LR7r for nrLSU region (1-1583) 34,35 . When these pairs of primers failed, both regions were amplified in two parts: in the case of ITS nrDNA, the region ITS1, including part of the 5.8S, with primers ITS5 and ITS2 33 , and part of 5.8S and the region ITS2 with primers ITS3 and ITS4 33 ; in the case of nrLSU, one region between the pair of primers LR0R and LR5 33 and another region between the primers LR3R and LR7r 35 were amplified. When neither direct nor amplification by parts gave good amplicons (above 20 ng/μL concentration), two semi-nested or nested PCR was used. For ITS nrDNA, a first amplification was done with ITS1F 36 and ITS4B 36 primers, amplifying part of the 18S and 28S nuclear ribosomal genes, and a second amplification was done with the pair of primers ITS5/ITS4 (nested-PCR) or only one of the inner primer and one external primer (semi-nested PCR). For nrLSU the first amplification was done with LR0R and LR7r primers, and the second amplifications were done with LR0R/LR5, and LR3R/LR7r primers. Individual reactions to a final volume of 25 μL were carried out using Illustra PureTaq Ready-To-Go PCR Beads (GE Healthcare, Buckinghamshire, UK) with a 10 pmol μL −1 primer concentration following the thermal cycling conditions used in Martín and Winka 37 . Negative controls lacking fungal DNA were run for each experiment to check for contamination.
The PCR products were subsequently purified using two different methods. When the quality of the DNA was low, due to the presence of multiple bands, QIAquick Gel Extraction kit (QiaGen) was used following manufacturer's instructions. When the quality of the DNA was high (a unique amplicon of above 20 ng/μL concentration), purifications were done using Exosap, Illustra ExoStar-1-Step (GE Healthcare, Buckinghamshire, UK) following the instructions of the manufacturers. Purified amplicons with a concentration of 20 ng/μL or more were sent to Macrogen (Korea) for Sanger sequencing with primers used in the amplification.
Phylogenetic analyses. Consensus sequences were obtained using Geneious version 9.0.2 http://www. genei ous.com 38 . Subsequently, they were subjected to a BLAST search with megablast option and compared against the sequences in the National Center for Biotechnology Information (NCBI) nucleotide databases 39 to check for contamination. Evaluation of EMBL/GenBank/DDBJ databases for ITS nrDNA and nrLSU sequences of a large set of Xylodon species was performed to provide a phylogenetic framework to X. australis and to maximize the molecular information available for these taxa (See Suplementary Table S1 on line). Three specimens of the sister genus Lyomyces P. Karst. were included as outgroup in the phylogenetic analyses 40 The maximum parsimony (MP), maximum likelihood (ML), and Bayesian inference analyses of specimens in the S1 Table are also deposited in protocols.io under the doi mentioned above. The ML and Bayesian analyses were done with the general time reversible model 41 , including estimation of invariant sites and a discrete gamma distribution with six categories (GTR + G), as selected by PAUP*Version 4.0b10.
Individual datasets of ITS nrDNA, nrLSU and a combined alignment of ITS nrDNA + nrLSU regions were used to compare specimens of Xylodon australis from South America and Australia with other Xylodon species. In the combined alignment, only samples with complete ITS nrDNA and nrLSU sequences were included to improve the resolution power of the obtained phylogeny.
Statistical tests for morphological traits. Basidia and spore morphology were analyzed in order to assess the morphological difference found in previous studies for Xylodon australis specimens from different locations 3 . One-way ANOVA tests were performed between X. australis lineages determined by molecular phylogenetic analyses.  43 with a burnin of 5000 trees for each run. Second, the subtree formed by the crown group for all Xylodon australis specimens analyzed was selected as a base tree for the next analyses. Five continuous morphological traits usually known as taxonomically informative for Hymenochaetales were selected: basidia length and width, spore length and width, and spore length-width ratio 45 . These traits were measured for the X. australis type specimen and for all samples in the subtree. To place the type specimen of X. australis into the molecular tree using the continuous morphological traits, the function locate.yeti from the R package phytools (v0.6.60) was used 16 . This function adapts the approach proposed by Felsenstein 14,16 to estimate phylogeny from continuous traits using a maximum likelihood framework. To include more than one continuous trait, a phylogenetic principal component analysis is performed first 46 . Then, these principal components are used to identify the optimal position of the type specimen in the phylogenetic tree applying the maximum likelihood criterion (see Equation 1 in Revell et al. 12 ). The method relies on the assumption that the characters have evolved along the tree and that the morphological differences between species are mostly due to inherited genetic differences.
In order to assess the performance of this approach, a randomization test with the original pure molecular subtree was conducted. We ran 100 replicates in which we pruned one tree-tip at random per replicate. Then, the tip was included again in the tree in two ways: randomly located or using the locate.yeti function. These two trees were compared with the original molecular tree by computing branch-score 47 and quadratic path distances 48 using the R package phangorn 49 . A lower distance means more similarity of the reconstructed tree with the actual molecular tree. In addition, to test the hypothesis that the inferred ML phylogenetic position of the type specimen through locate.yeti function is significantly better than alternative locations, a likelihood ratio test was conducted by comparing the likelihood score of tree in which type position was constrained, to trees with unconstrained type locations. Simulations of continuous traits were performed on the constrained tree, and then type position was inferred without constraint. These simulations were used to generate a null distribution to check for significance of the likelihood ratio test 12 .
The same methodology was applied to each specimen from New Zealand (PDD Herbarium) for which no molecular data were obtained.