Abstract
Our understanding of fungal diversity is far from complete. Species descriptions generally focus on morphological features, but this approach may underestimate true diversity. Using the morphological species concept, Hesperomyces virescens (Ascomycota, Laboulbeniales) is a single species with global distribution and wide host range. Since its description 120 years ago, this fungal parasite has been reported from 30 species of ladybird hosts on all continents except Antarctica. These host usage patterns suggest that H. virescens could be made up of many different species, each adapted to individual host species. Using sequence data from three gene regions, we found evidence for distinct clades within Hesperomyces virescens, each clade corresponding to isolates from a single host species. We propose that these lineages represent separate species, driven by adaptation to different ladybird hosts. Our combined morphometric, molecular phylogenetic and ecological data provide support for a unified species concept and an integrative taxonomy approach.
Similar content being viewed by others
Introduction
What is a species? This is a perennial question in evolutionary biology. The answer is complex and has been intensely argued for decades. Different species concepts corresponding to multiple biological properties provide a means to recognize, delineate and describe species. These properties include differences in morphological traits, nucleotide divergence and monophyly, reproductive isolation, ecological niches or adaptive zones, mate recognition or mating systems, geographic range, exclusive coalescence of alleles, etc. However, biologists from various research fields have advocated different and sometimes incompatible species concepts, leading to varying conclusions regarding delimitation of species and their numbers. Rather than disagreeing on the conceptual agreement of what is a species (a separately evolving metapopulation lineage1), de Queiroz2,3 argues that each species concept emphasizes different properties. In evolutionary biology, “species” are hypotheses for which evidence can be sought by the study of multiple properties3. The absence of a certain property does not provide evidence contradicting any given species hypothesis. This is the unifying species concept3.
Fungi have essential functions in ecosystems, they are virtually everywhere, even in the most extreme habitats, and associate with diverse organisms (algae, plants, invertebrates and other fungi)4. Currently, about 135,000 species of fungi have been described5; still many localities, habitats and taxonomic groups remain poorly sampled. In the pre-molecular era, Hawksworth6 estimated the number of fungal species to be 1.5 million, based on the ratio of vascular plants and fungi on the British Isles, which he accepted as 1:6. An ITS-based evaluation of soil fungal diversity of two temperate plots and the vascular plant richness in those plots led to an extrapolation of global (soil) fungal species richness estimates from 3.5 to 5.1 million7. Taylor and colleagues8, using a large fungal dataset from a boreal ecosystem with well-established plant diversity, suggested up to 6 million species of fungi as a global estimate. Understanding how these millions of fungal species have come to existence has stimulated widespread interest. The challenges of diversity studies are posed especially for fungi, which produce propagules that are microscopic in size, have sometimes worldwide distributions and use a multitude of host species. Pringle and colleagues9 postulated that morphology is a poor means to distinguish species of this magnitude given these dispersal potentials and patterns of host usage. As a result, species hypotheses about microscopic organisms with global distributions or multiple host ranges should be treated with care.
Many fungal species have been described based on morphological traits; representatives of any given species share a set of morphological characteristics. However, this morphological species concept is a poor means of species delimitation when phenotypic plasticity allows for overlapping morphologies in distinct species or when morphological traits have not yet arisen in the process of speciation or incipient speciation. For example, the genus Protoparmelia s. str. (Ascomycota, Lecanorales) consists of 12 species based on morphological and chemical features but a phylogenetic-coalescent approach recognizes 23 species10. Another widely cited example is that of Dictyonema glabratum (Basidiomycota, Agaricales), a single morphological species that constitutes 126 species using a Generalized Mixed Yule Coalescent (GMYC) analysis of a large dataset of the internal transcribed spacer (ITS) DNA region, and even more than 400 species based on predictive modeling11.
Many species of fungi form associations with other organisms and these associations may be critical in species recognition. As a result, fungal species may be circumscribed based on the property of host associations. Host specificity represents an ecological condition; it entails resource availability and niche specialization. The concept of “ecological species” generally refers to reproductive isolation evolved through adaptation to different environments. The micro-evolutionary process of natural selection among diverging populations or subsets of a single population acts in contrasting directions between environments and leads to the fixation of alleles, which may be advantageous in one environment but not in others12,13,14. The ecological species concept dates from the 1940s, when Theodosius Dobzhansky15 wrote that “[s]peciation in Drosophila proceeds mainly through evolving physiological complexes which are successful each in its environment.” An interesting case is Leccinum (Basidiomycota, Boletales), a genus of ectomycorrhizal fungi forming associations with many plant hosts. A recent study16 found high host specificity in all species included, except for the generalist L. aurantiacum. The authors raised the point that ecological information on its own (“the ability to grow on a new host”) does not a priori provide evidence for a species hypothesis. More recently, three species were described within the ant-parasitic Ophiocordyceps unilateralis species complex (Ascomycota, Hypocreales) based on the combination of molecular, micro-morphological and ecological (host specificity) data17. All this is in line with de Queiroz’s (2007) view3 that multiple properties provide evidence for lineage separation, that is, divergence of populations, and, thus, speciation.
In this paper, we explore species limits in an enigmatic group of microscopic fungi, the Laboulbeniales. These are obligate, microscopic ectoparasites of arthropods. Around 2200 species in 141 genera are known to infect various groups in three arthropod subphyla (Chelicerata, Hexapoda, Myriapoda) and they are known from all continents except Antarctica18,19,20,21. Laboulbeniales never form mycelia; the ascospores do not form germ tubes but rather divide mitotically after attachment to the host to form thalli of up to thousands of cells by determinate growth. Hesperomyces virescens (Figs 1, 2) is a species of Laboulbeniales that has been reported to parasitize around 30 species of ladybirds (Coleoptera, Coccinellidae) in all continents but Antarctica. It grows exclusively on adult ladybird hosts in 21 genera in 5 subfamilies21. Since its discovery on the invasive ladybird Harmonia axyridis, biologists have discussed H. virescens as a candidate model for studying host-parasite co-evolution and biological control programs22. Contrary to what is the consensus for Laboulbeniales, H. virescens has been reported to have negative effects on its hosts21,22. Based on intra- and interspecific transmission experiments, Cottrell and Riddick23 suggested that different lineages of H. virescens exist and that each of these lineages may have a high degree of host specificity. The question whether H. virescens truly is a single species or an assemblage of morphologically similar but genetically distinct species19 has provided the starting point for the present research study. To identify H. virescens, mycologists have used morphological characters that can be compared across a range of different host species. Here, we combine morphological, molecular and ecological data as independent lines of evidence to infer the number of species within H. virescens, following the unified species concept proposed by de Queiroz2,3.
Hesperomyces virescens parasitizing Harmonia axyridis. A group of thalli of Hesperomyces virescens attached to the elytral tips of a specimen of Harmonia axyridis (from Watermaal-Bosvoorde, Belgium). Image provided by Gilles San Martin. Inserted is a drawing of a single adult thallus (D. Haelew. 601c, FH 00313615, Byron, Georgia, USA), with the following structures annotated: haustorium; cells I, II, and III of the receptacle; the appendage (app.) with antheridia (anth.); and the perithecium (perith.) with its terminal projections (proj.). Drawing by André De Kesel.
Results
Morphometric approach
Details of measurements and ratios for all 181 thalli included in the analysis are given in Supplementary Table S1. For a majority of variables, the best model to explain differences in measurements contained host species (Mod1 in Table 1). Inclusion of host species as an explanatory variable considerably improved model performance. Of the 35 studied variables, 10 did not differ significantly between host species: W cell I, L cell II, W cell II, L cell VI, lngst. proj., L/W cell II, total L rec./total L, tier II/L perith., tier III/L perith., and lngst. proj./L perith.
We only considered ratios for PCA to focus on shape rather than natural variation in absolute size. Significant differences were observed for the following ratios: L/W cell I, L/W cell III, total L app./total L, L/W cell VI, tier IV/L perith., lobes/L perith., L/W perith., and L perith./total L. Statistical processing of these ratios revealed two principal components (PCs) that together accounted for 81.54% of the observed variation in thallus morphology of H. virescens between C. propinqua, H. axyridis and O. v-nigrum. PC1, 48.39% variation explained, represents L/W cell I, L/W perith., and L/W cell VI (Fig. 3a,b). PC2, 33.15% variation explained, represents L/W cell VI and L/W cell I (Fig. 3a,c). In the morphospace formed by the two first PCs, clouds of individuals from the 3 different host species overlap partly, but they also occupy a considerable part of the morphospace without overlap (Fig. 4).
Results of Principal Component Analysis. (a) Morphospace formed by the first two PCs of the PCA showing the importance of ratios. (b,c) Contributions of included ratios to PC1 (b) and PC2 (c) separately. The dashed line is a reference corresponding to the expected value if the contributions were uniform. Contributions above the reference line are considered as important.
PCA using morphometric variables showing variation in thallus shape. Each symbol represents an individual thallus in the two-dimensional morphospace formed by the first two PCs. Thalli are colored by host species (yellow circles Cheilomenes propinqua, green triangles Harmonia axyridis, purple squares Olla v-nigrum).
Nucleotide alignment datasets
We generated 93 Hesperomyces sequences during this study, of the SSU (31), ITS (37) and LSU (25) regions. Our ITS dataset comprised 1068 characters, of which 769 were constant and 229 were parsimony-informative. A total of 51 ITS sequences were included, of which 35 have been newly generated during the course of this study, complemented by 16 sequences that we retrieved from GenBank (http://www.ncbi.nlm.nih.gov/genbank/): 43 sequences of H. virescens, 2 of H. coleomegillae, 2 of H. palustris and 4 of Herpomyces spp. (outgroup). Isolates of H. virescens originated from 9 host species: Adalia bipunctata (5 isolates), A. decempunctata (2), Azya orbigera (1), Cheilomenes propinqua (5), Cycloneda sanguinea (2), Halyzia sedecimguttata (1), Harmonia axyridis (16), Olla v-nigrum (8) and Psyllobora vigintimaculata (3). Our LSU dataset consisted of 22 newly generated sequences of H. virescens complemented by 14 sequences downloaded from GenBank (9 of H. virescens, 4 of Herpomyces spp. and 1 of Pyxidiophora cf. microspora as outgroup), and 1051 characters, of which 769 were constant and 222 were parsimony-informative. Isolates of H. virescens originated from 7 host species: Adalia bipunctata (5 isolates), A. decempunctata (1), Azya orbigera (1), Cheilomenes propinqua (2), Harmonia axyridis (13), Olla v-nigrum (6) and Psyllobora vigintimaculata (3).
Our concatenated SSU+ITS+LSU dataset included 2767 characters and 52 isolates representing 6 species (GenBank accession numbers in Table 2). Of all characters, 2183 were constant and 446 were parsimony-informative. Taxonomic sampling covered 3 genera in the Laboulbeniomycetes: Herpomyces, Hesperomyces and Pyxidiophora. In addition to 43 isolates of Hesperomyces virescens, we included Herpomyces chaetophilus, H. periplanetae and Pyxidiophora cf. microspora (outgroup). Isolates of H. virescens originated from 9 different host species: Adalia bipunctata (5 isolates), A. decempunctata (2), Azya orbigera (1), Cheilomenes propinqua (5), Cycloneda sanguinea (2), Halyzia sedecimguttata (1), Harmonia axyridis (16), Olla v-nigrum (8) and Psyllobora vigintimaculata (3).
Phylogenetic inferences
The Hesperomyces clade has maximum support in the ITS and LSU datasets. In both datasets, each monophyletic clade within the H. virescens complex consists of isolates from thalli removed from a single host species. There is one exception, the clade consisting of isolates from two host species in the same genus, Adalia bipunctata and A. decempunctata. Henceforward, in the text we will refer to these distinct clades by the first letters of the host genus and species names: Ab for Adalia bipunctata, Ad for Adalia decempunctata, Ao for Azya orbigera, Cp for Cheilomenes propinqua, Cs for Cycloneda sanguinea, Ha for Harmonia axyridis, Hs for Halyzia sedecimguttata, Ov for Olla v-nigrum and Pv for Psyllobora vigintimaculata. In addition, figures are consistently color-coded by host species where appropriate.
In the ITS dataset, 8 clades are recognized, in addition to H. coleomegillae and H. palustris, which are positioned basally compared to the other clades (Fig. 5a). Of the 8 clades, 6 are strongly supported (BS ≥ 84): Ab+Ad, Cp, Cs, Ha, Ov and Pv. Only the singleton clades (Ao and Hs) lack support. In the LSU dataset, 6 clades are recognized (Fig. 5b). All of these clades have strong support (BS ≥ 81). In the three-gene dataset, again, the Hesperomyces clade has strong support and consists of 10 clades (Fig. 6). Of these, nine have high support (ML BS ≥ 93, pp ≥ 0.94). Only the Hs clade is unsupported in the ML analysis (BS = 65) whereas the Bayesian analysis provides moderate support (pp = 0.82). This clade consists of only a single isolate for which we were able to generate a 784-bp ITS sequence only. ML inference agrees fully with the Bayesian analysis. In both reconstructions, H. coleomegillae and H. palustris are placed basally compared to the other clades within the genus Hesperomyces.
Phylogenetic reconstruction of Hesperomyces virescens sensu lato. (a) Phylogram reconstructed from the ITS dataset. (b) Phylogram reconstructed from the LSU dataset. Topologies are the result of maximum likelihood inference. For each node, the ML bootstraps (if >70) are presented above the branch leading to that node. The arrowhead in both reconstructions points at the node leading to the Hesperomyces clade. Monophyletic clades are color-coded by host species. All host species included in the analysis are pictured left and right of the phylogenies. Photo credits: Adalia bipunctata, A. decempunctata, Halyzia sedecimguttata, Harmonia axyridis, Gilles San Martin (Flickr); Azya orbigera, Pavel Kirillov (Wikimedia Commons); Cheilomenes propinqua, Sally Adam (iNaturalist); Cycloneda sanguinea, Damon Tighe (iNaturalist); Olla v-nigrum, Roberto Güller (Flickr); Psyllobora vigintimaculata, Katja Schulz (iNaturalist).
Three-gene phylogeny showing results of sequence-based species delimitation methods. Maximum clade credibility tree of the Hesperomyces virescens complex, reconstructed from the concatenated SSU+ITS+LSU dataset. The tree is the result of a Bayesian analysis performed in BEAST. For each node, the ML bootstraps (if >70) and posterior probabilities (if >0.70) are presented above/below the branch leading to that node. Monophyletic clades are color-coded by host species. To the left of the phylogeny, all host species included in the analysis are pictured. Photo credits: Adalia bipunctata, A. decempunctata, Halyzia sedecimguttata, Harmonia axyridis, Gilles San Martin (Flickr); Azya orbigera, Pavel Kirillov (Wikimedia Commons); Cheilomenes propinqua, Sally Adam (iNaturalist); Cycloneda sanguinea, Damon Tighe (iNaturalist); Olla v-nigrum, Roberto Güller (Flickr); Psyllobora vigintimaculata, Katja Schulz (iNaturalist). To the right of the phylogeny, the results of species delimitation analyses are summarized, from left to right: ABGD analysis of the SSU+ITS+LSU alignment; bPTP analysis of the ITS topology; bPTP analysis of the LSU topology; and GMYC analysis of the ITS, LSU and SSU+ITS+LSU ultrametric trees (without outgroups) generated in BEAST. Hatching implies lack of support, whereas no coloration means that clade was absent in that analysis.
Species delimitation analyses
Results of the sequence-based methods for species delimitation are summarized in Fig. 6 and Table 3. The ABGD analysis resulted in 8 distinct groups, irrespective of distance metrics or gap width values. Eight species were identified within the Hesperomyces virescens clade from the bPTP analysis of the ITS topology. The bPTP analysis of the LSU topology resulted in 6 species (no LSU sequences were generated for isolates from Cycloneda sanguinea and Halyzia sedecimguttata). Support was lacking for the Cp clade in this analysis. The GMYC analyses of the ITS and concatenated trees resulted in 8 species, all with high support except for clades Ao and Hs (which comprised a single isolate only). The GMYC analysis of the LSU tree resulted in 6 recognized species, without support for Ao, Hs (each comprising a single isolate), and Cp and Ha.
Discussion
Our results illustrate that H. virescens encompasses several unique genetic lineages. Each of these lineages occurs on a single host species (or two host species in the case of Adalia), regardless of geographic origin of the collection. Some of the clades in our phylogenies have no or only moderate support but these are the clades for which only a single isolate is available (Ao and Hs). Some ladybird species were only recently discovered as hosts and others are not frequently found. For example, during fieldwork in Panama in 2015 we found Hesperomyces thalli on Azya orbigera, which had previously not been reported as a host. Out of 151 ladybirds, only 10 infected ones were found, each individual carrying a single thallus (1 individual carried 3 thalli). We tried two extraction protocols, each with 1 thallus. The extraction using the Extract-N-Amp PCR Plant Kit failed; the one with the REPLI-g Single Cell Kit was successful and we were able to generate SSU, ITS and LSU sequences (isolate D. Haelew. 928 g). In another case, a single Hesperomyces-infected individual of Halyzia sedecimguttata was found in the Netherlands in 2015. It bore 13 adult and 4 juvenile thalli. We chose to use 10 adult thalli for DNA isolation with the Extract-N-Amp Plant PCR Kit (isolate D. Haelew. 955b) but for this isolate we were only able to generate an ITS sequence. Since this report, no further infected specimens of H. sedecimguttata have been collected. Finally, we only had 3 infected individuals of Psyllobora vigintimaculata available for study, but these specimens carried sufficient thalli for both morphological study and molecular work.
Some of our host species were collected only from a single population. This is the case for Cheilomenes propinqua, Olla v-nigrum and P. vigintimaculata. However, specimens of O. v-nigrum, although from a single locality, were collected on multiple occasions in 2014, 2015 and 2016 (also from a laboratory colony for many generations). Adalia bipunctata and H. axyridis are the host species with the widest geographical range included in this study. Infected specimens of A. bipunctata were collected in Denmark, Italy and Sweden; specimens of H. axyridis were collected on different continents. Even so, both clades Ab+Ad and Ha each form two monophyletic lineages, in all datasets (ITS, LSU, SSU+ITS+LSU). In other words, there is no geographic signal. We conclude that phylogenetic structure is primarily determined by host specialization. Based on intra- and interspecific transmission experiments, Cottrell and Riddick23 proposed that “isolates/strains of H. virescens may exist under field conditions and only infect closely related Coccinellidae or even a single species.” Based on the results of our species delimitations analyses, we propose that these lineages (or clades, as we refer to them) represent distinct species.
The clade Ab+Ad is the only clade consisting of H. virescens isolates taken from more than a single host species. In this case, the isolates originated from either Adalia bipunctata or A. decempunctata. Ladybirds in the genus Adalia are popular for studies in ecology and population genetics because they vary distinctively in their elytral and pronotum color patterns. Each of these color forms, or phenotypes, corresponds to a unique genotype24. Adalia bipunctata and A. decempunctata are very similar in ecology and habitat use; both coexist throughout temperate regions in Europe. They are also closely related genetically, and even some A. bipunctata × A. decempunctata hybrids have been reported25,26. This Ab+Ad clade makes for an interesting case because one could pose the question whether the specificity of H. virescens relates to host species-level or genus-level. Based on the currently available data we cannot solve this question, but we will continue to collect ladybirds and thus increase our dataset of ladybird host species.
Comparison of species delimitation methods
Whereas molecular data provide a valuable tool to validate morphology-based species descriptions, the application of species delimitation methods can increase confidence if several methods offer congruent estimates of species diversity within a given dataset. Incongruences in results imply that multiple methods differ in their delimitation power. Alternatively, it is also possible that users make incorrect assumptions when employing a given species delimitation approach27. In the event of incongruent results, it is better to be conservative, rather than to designate entities that do not actually represent evolutionary metapopulation lineages as species. In any case, the multiple species delimitation analyses that we used in our study identified congruent species boundaries. The combination of BS, pp and species delimitation support provides strong evidence for H. virescens being a complex of multiple species.
In the bPTP analysis of the LSU topology, the Cp clade was not supported. PTP models speciation in terms of number of nucleotide substitutions28. Upon manual inspection of the multiple alignment, for all clades it is the case that the number of nucleotide substitutions within the clade is zero. The only exception is the Cp clade, with 2 substitutions between the isolates of this clade. The number of substitutions between Cp and other clades ranges from 5 to 13. The PTP model did probably not interpret the Cp clade as a distinct species because of these within-clade substitutions for the Cp clade. Our PTP analyses based on single-gene trees are consistent with the results obtained by ABGD. We also performed a bPTP analysis on the concatenated SSU+ITS+LSU dataset (only Bayesian support values shown, Table 3). Although the number of species is the same as in the PTP analysis of single gene topology and the ABGD and GMYC approaches, the Bayesian support dropped significantly; none of the delimited clades have support higher than 0.52. We repeat earlier findings28,29 that PTP is most accurate with single gene trees.
Kekkonen and Hebert30 put forward that GMYC usually delimits more species compared to other methods. In our analysis, the results from GMYC are congruent with ABGD and bPTP of the ITS topology (and bPTP of the LSU topology, noting that two clades were missing in this analysis). Two clades that lack support – Ao and Hs – consist of single isolates only. GMYC looks at intraspecific branching versus interspecific branching, and thus it is no surprise that these singleton clades have no support from this approach. The low support for the Ha clade in the GMYC analyses of the LSU and three-gene topologies may be explained by the fact that for a number of isolates coming from Harmonia axyridis, sequences are incomplete (missing or only partial SSU, ITS or LSU). Because this clade holds many isolates, missing sequence data for a number of these isolates may influence generating an ultrametric tree, which is a computationally intensive and error-prone process. Since GMYC is dependent on the accuracy of this input tree, any alterations will strongly influence species delimitation analyses.
Upon the introduction of PTP, it was noted that delimited groups represent “putative” species28. PTP uses phylogenetic reconstructions inferred from single gene datasets, which are gene trees rather than species trees. Also GMYC is based on a single-gene tree. As a consequence, more data should be collected to confirm and validate the species boundaries set by these delimitation approaches, in an integrative taxonomy framework across disciplines31,32,33,34,35. Note that this framework is in line with the unified species concept3.
Comparison of ITS and LSU as barcode markers
Molecular identification of fungi relies on the availability of good DNA barcode markers. Currently, DNA-based identifications focus on genes that code for ribosomal RNA (rDNA), because these regions have many copies in the genome and thus are well-suited target regions for PCR amplification. The internal transcribed spacer (ITS) region has been proposed as universal barcode for Fungi36. This means that for a majority of fungi, the interspecific variation at this marker should exceed the intraspecific variation, and for over 70% of fungi the ITS is indeed effective in recognizing species. A number of considerations have been made since the acceptance of this barcode marker36,37,38,39: (1) RPB1 is actually better in discriminating species but its amplification success is much lower; (2) whereas ITS performs best overall across the fungal kingdom, its identification power is equal to LSU for subphyla Pezizomycotina and Saccharomycotina (Ascomycota); (3) ITS does not contain enough variation to discriminate between species for some groups of fungi, such as Aspergillus and economically important plant pathogens in the genera Alternaria, Diaporthe, Fusarium and others; (4) arbuscular mycorrhizal-forming species in Glomeromycota are multinucleate and extremely intraspecific divergent in their ribosomal DNA. These challenges have driven mycologists to developing other, lineage- or genus-specific barcodes. These secondary barcodes are often more difficult to amplify (because of the lack of universal primers) but have a better delimiting power than the ITS40. Examples of secondary barcodes include calmodulin (CaM) for Aspergillus; the translation elongation factor (TEF1), topoisomerase I (TOP1) and phosphoglycerate kinase (PGK) for Fusarium; the LSU rDNA region for Amanita; and the Apn2-Mat1–2 intergenic spacer and partial mating type (Mat1–2) gene (ApMat) and glutamine synthetase (GS) combined for Colletotrichum38,41,42.
We experience low amplification success for the ITS region with the Laboulbeniales using general fungal-specific primers. There are many possible reasons for failed ITS amplification, ranging from simple primer mismatch as is the case in the Archaeorhizomycetes43 to significant intragenomic44,45 or intraspecific variability46. Although the 5.8S region is highly conserved, both spacer regions (especially ITS1) appear to be rapidly evolving47. Previously generated sequences of Laboulbeniales suggest that the ITS differs significantly among genera and we have currently no idea of the extent of this variability. Still, ITS may be useful and important as a marker to study intrageneric relationships. During our studies of Hesperomyces we designed and are routinely using primers that specifically target conserved regions of the ITS (ITShespL, 5′-CTCCTGTAGAACCTACACATC-3′ and ITShespR, 5′-CAAATTTAAGCTTTTGCCGC-3′).
During the course of this study, we constructed phylogenies based on single genes (SSU, ITS, LSU) and on a concatenated SSU+ITS+LSU dataset. The SSU gene is very conservative and has no discriminative power at the species level. But, both the ITS and LSU datasets result in high support for the individual clades of the H. virescens complex (Fig. 6). In addition to its discrimination power, amplification of the LSU region poses virtually no problem within the Laboulbeniomycetes so far investigated. The commonly used fungal primers for the LSU region, such as LR0R/LR5 and LIC24R/LR348,49, generally work well for most species of Laboulbeniomycetes. Based on these preliminary results, the LSU region should be further investigated as barcode for species delimitation in Laboulbeniomycetes.
Hesperomyces virescens, a complex of cryptic species?
Recent molecular (phylogenetic) studies point at a dazzling diversity of the kingdom Fungi. However, it is not always possible to infer this diversity from morphological features. Species that “have been classified as a single nominal species because they are at least superficially morphologically indistinguishable” are referred to as cryptic species50. Many (or almost any51) species studied using molecular, incompatibility, secondary metabolites have been shown to “mask” several biological species. Examples are found in diverse groups of fungi – among Ascomycota: Eurotiales (Aspergillus9), Helotiales (Phialocephala52), Lecanorales (Protoparmelia10); among Basidiomycota: Agaricales (Cortinarius53; Tricholomopsis54), Polyporales (Ganoderma55), Russulales (Lactifluus56, Russula57), Ustilaginales (Tranzscheliella58). In this study, we have observed at a superficial level cryptic species. Consequently, we employed landmark-based geometric morphometry59 followed by principal component analysis of shape variation, aimed at finding characters, if any, to circumscribe species with.
To date, morphometric methods in Laboulbeniales have only been applied in studies dealing with the genus Laboulbenia. Morphological plasticity of L. flagellata from different carabid hosts (Coleoptera, Carabidae) was studied by De Kesel and Van Den Neucker60. The general habitus of thalli was stable, but size was related to host species, habitat of the host and position of the fungus on the host integument. Subsequently, De Kesel and Haelewaters61 tested differences in thallus shape and dimension between two morphologically similar species of Laboulbenia. Most variables were significantly different between both species, particularly the shape of the receptacle was different regardless of size or growth position. In this study, we generated the largest morphometric dataset to date for any species of Laboulbeniales, including measurements and ratios of 181 thalli from 3 host species. Our PCA suggests that the shapes of cell I, cell VI and the perithecium contribute most to the observed variation within the dataset. If we were to formally describe clades Cp, Ha and Ov as separate species, we expect to find most descriptive features in cell I of the receptacle and in the perithecium and its basal cell (VI).
In both the ITS and three-gene phylogenies, we retrieved H. coleomegillae and H. palustris62 (on Coleomegilla maculata) basal to the H. virescens clade. The basal-most clade within H. virescens is Ao in the ITS and three-gene analysis. In terms of morphology, the thalli on Azya orbigera are structurally quite different compared to thalli from other host species. For example, the appendage is only 3-celled with the third cell carrying two antheridia. This structure is completely different than thalli from any of the other hosts in our dataset and also than H. coleomegillae and H. palustris21,62,63,64. The appendage of these thalli consists of a single row of at least 4 cells, and every cell starting from the second carries an antheridium (the distal-most cell carries 2 antheridia). In the LSU phylogeny, the Ov clade is basal-most positioned. However, this tree lacks internal support, except for the sister relationship between clades Ao and Pv. Also the ITS phylogeny lacks internal support and all internal clades collapse to a basal polytomy. On the other hand, incongruence between trees is possible as a result of incomplete lineage sorting. Gene trees are not equal to species trees; a multi-locus approach is preferred for investigating species divergence and gauging relationships between species65. In the three-gene phylogeny, the backbone relationships within the H. virescens clade are well resolved, as indicated by high pp values. The sister relationship between clades Cs and Ov is without support, as is that between clades Ab+Ad and Cp. This may indicate that there is a lack of taxon sampling, which is no surprise given the broad host range of H. virescens21. Given the available data, we accept the clades Ab+Ad, Ao, Cp, Cs, Ha, Hs, Ov and Pv to be part of the H. virescens complex, or H. virescens sensu lato, and to represent independent evolutionary lineages. Following the strict host specificity detected by this study, we propose to restrict H. virescens sensu stricto to those thalli found on Chilocorus stigma, the host species on which the fungus was originally described64.
Conclusions
Through DNA isolation, PCR amplification, sequencing and analysis methods, thousands of characters became available for minute fungi that do not have many morphological features and do not grow in culture. These remarkable improvements in the collection of character data will help us answer questions about the validity of “worldwide” and “cosmopolitan” geographic distributions ascribed to many morphological forms of the Laboulbeniomycetes. Here, we provided answers in the case of Hesperomyces virescens, which we have shown to be a complex of different species, each with its own host (genus). We are only starting to unravel patterns of speciation in this group of fungi. The findings of this paper are not only promising for future studies, but they also emphasize the necessity for an integrative approach in taxonomic research. We hope with this contribution to include the Laboulbeniales ectoparasitic fungi in contemporary discussions considering molecular evolution and speciation patterns, rather than treating them as obscure fungi for specialists only.
Methods
Collection of hosts
Our main field site for the collection of ladybirds was the 480-ha land of the USDA Southeastern Fruit and Tree Nut Research Laboratory in Georgia, USA66, where we collected specimens of Harmonia axyridis and Olla v-nigrum in 2014–2015. In addition, ladybirds were collected by the first author or collaborators at different sites in four continents: Africa, Asia, Europe and North America (details in Supplementary Table S2). Sampling of ladybirds was done using a variety of standard entomological methods: Tedders pyramidal traps, light traps, hand collecting and sweeping in stands of weedy vegetation along the banks of swamps and small lakes and at the sides of roadways. Long-time preservation of ladybird specimens was in 95% ethanol at −20 °C. In addition to field-collected material, pinned ladybirds in dried insect collections were screened for the presence of Laboulbeniales. The Coccinellidae collection of the Boston Harbor Islands All Taxa Biodiversity Inventory project at the Harvard Museum of Comparative Zoology (Cambridge, MA), the Division of Invertebrate Zoology at American Museum of Natural History (New York City, NY) and the Florida State Collection of Arthropods (Gainesville, FL) were primary sources for infected ladybirds.
Collection of Laboulbeniales
Preserved insects were examined for the presence of Laboulbeniales under a dissecting microscope at 10–50× magnification. Hesperomyces thalli were removed from their hosts using Minuten Pins (BioQuip, #1208SA, Rancho Dominguez, CA) inserted onto wooden rods. Following standard procedure by Benjamin67, we removed thalli or groups of thalli and mounted them in Amann’s medium, a liquid solution. Before applying Amann’s medium and to facilitate microscopic observations, thalli first had to be arranged and fixed onto the microscope slide. To make thalli a bit sticky, they were first placed in a droplet of Hoyer’s medium (30 g arabic gum, 200 g chloral hydrate, 16 mL glycerol, 50 mL ddH20). Next, thalli were individually picked up and arranged in one or two rows. After a brief period of drying, the slide was closed using a cover slip with a drop of Amann’s medium (drop facing downward) and subsequently sealed with nail polish or B-72 in acetone (Gaylord, #AB72, Syracus, NY). We viewed mounted specimens at 400–1000× magnification using an Olympus BX40 microscope equipped with an XC50 camera (Olympus, Waltham, MA). Identification was done using descriptions from Thaxter64,68, Santamaria69 and De Kesel63. Slides are deposited at Farlow Herbarium (FH; Harvard University, Cambridge, MA).
Morphological studies
To assess morphological variation in thalli we took measurements of 22 parameters per thallus (Fig. 1): total length of the thallus including haustorium (total L w foot, point a—point x in Fig. 1), total length of the thallus (total L, b—x), length of cell I (L cell I, b—d), width of cell I (W cell I, c—o), length of cell II (L cell II, m—o), width of cell II (W cell II, l—n), length of cell III (L cell III, d—f), width of cell III (W cell III, e—l), total length of receptacle (total L rec., b—f), length of basal cell of the appendage (L bas. app., f—g), total length of appendage (total L app., f—k), length of longest antheridium (L lngst. anth., h—j), length of longest antheridial neck (L anth. neck, i—j), length of cell VI (L cell VI, m—z), width of cell VI (W cell VI, p—y), perithecium length (L perith., w—z), perithecium width (W perith., r—x), length of second tier of perithecial wall cells (tier II, q—r), third tier (tier III, r—s), fourth tier (tier IV, s—t), length of lobes (lobes, t—w) and length of longest projection (lngst. proj., u—v). To correct for natural variation in length and width, these ratios were calculated: L/W cell I, L/W cell II, L/W cell III, total L rec./total L, total L app./total L, L/W cell VI, tier II/L perith., tier III/L perith., tier IV/L perith., lobes/L perith., L/W perith., L perith./total L and lngst. proj./L perith.
Measurements were made at 400–1000× magnification with cellSens Standard 1.14 software (Olympus) using the Polyline measuring tool. We measured at least 30 adult thalli from each host population. Maturity was judged by the presence of ascospores within the perithecium. To exclude potential position-induced morphological variation62, only thalli from the elytra were measured and used in this study.
We analyzed variation in morphology of thalli from different host species and populations using generalized linear mixed models (GLMM), implemented in the R package lme470. Random effects for insect specimen were included, because we measured several thalli from the same host individuals. Hypothesis testing was done using likelihood ratio tests, with P-values calculated based on chi-squared distributions, declaring an effect significant when P ≤ 0.05. Two models were compared for each variable, the null model (Mod0) and the model with host species as explaining variable (Mod1). Model selection happened using the Akaike Information Criterion71. For a selection of variables with significant differences between host species in the GLMMs, principal component analysis (PCA) followed by exploratory biplots were made. PCA was only done for ratios to visualize variation in shape and structure independent of size. PCA and biplots were obtained using the R package factoextra72.
DNA extraction methods
We extracted DNA from 1–18 Hesperomyces thalli either using the QIAamp DNA Micro Kit (Qiagen, Stanford, CA), a modified Extract-N-Amp Plant PCR kit (Sigma-Aldrich) procedure73 or a modified REPLI-g Single Cell Kit (Qiagen, St. Louis, MO) protocol74. The QIAamp DNA Micro Kit protocol was followed as per the manufacturer’s instructions. One major change we implemented was the increase of the incubation time at 56 °C for complete lysis, to several days. With the Extract-N-Amp Plant PCR kit, thalli were removed at the foot with a tiny drop of Hoyer’s medium or glycerin at the tip of a Minuten Pin and placed in a 0.5 µL PCR tube filled with 20 µL of Extraction Solution. The tube was incubated at room temperature for 10–30 min and then at 95 °C for 20 min. The extract was then diluted with 60 µL of Dilution Solution (3% Bovine Serum Albumin). The REPLI-g Single Cell Kit is different from the previous protocols because it adds a whole-genome amplification (WGA) step to the DNA isolation, thus providing a considerable benefit when material is scarce. A Minuten Pin was submerged in glycerin to remove a single thallus from its host and place it in a droplet of glycerin on a microscope slide. The thallus was carefully placed in a 0.2 mL PCR tube with 2 µL of phosphate-buffered saline (PBS). These steps were done at 40× magnification under a stereomicroscope. After adding 1.5 µL of prepared D2 buffer, the tube was incubated at 65 °C for 20 min. Subsequent steps followed the manufacturer’s instructions. All steps of this procedure were performed under a laminar flow hood to ensure sterile conditions.
For a majority of our isolates, we applied pre-treatments to increase the likelihood of successful isolation and subsequent PCR amplification. These pre-treatments included subsequent cycles of freezing on liquid nitrogen and heating to 95 °C, prolonged incubation at 56 °C in 180 µL ATL buffer + 20 µL proteinase K or in 20 µL Extraction Solution using a Shake ‘N Bake Hybridization Oven (Boekel Scientific, model no. 136400–2, Feasterville, PA) and homogenization in a FastPrep FP120 Cell Disrupter at 5.0 m/sec for 15 s (Thermo Fisher Scientific, Waltham, MA). For both the QIAamp DNA Micro Kit and the Extract-N-Amp Plant PCR kit, we often manually crushed thalli in 1.5 mL tubes using a 1.5 mL pellet pestle (Kimble, #749521-1500, Vineland, NJ). In the REPLI-g Single Cell Kit, the single thallus was often cut in half through the perithecium using a sterile no. 10 surgical blade on a disposable Bard-Parker handle (Aspen Surgical, Caledonia, MI) before placing it in the 0.2 mL PCR tube.
PCR amplification and DNA sequencing
We amplified the nuclear small and large ribosomal subunits (SSU and LSU) and the internal transcribed spacer region of the ribosomal DNA (ITS). Primer combinations used were NS1/NS2, NS1/NS4, SL122/SR4, NSL1/NSL2, SL122/NSL2 and SL344/NS6 for SSU73,75,76; ITS1f/ITS4, ITS1f/ITS4A and ITShespL/ITShespR for ITS74,75,77,78; and LIC24R/LR3 and LR0R/LR5 for LSU48,49. PCR reactions consisted of 2.5 µL of each 10 µM primer, 13.3 µL of RedExtract Taq polymerase (Sigma-Aldrich), 5.7 µL of ddH2O and 1 µL of DNA extract. For all amplifications, an Eppendorf Mastercycler ep gradient thermocycler was used with initial denaturation at 94 °C for 3:00 min; followed by 35 cycles of denaturing at 94 °C for 1:00 min, annealing at 50 °C for 0:45 min, extension at 72 °C for 1:30 min; and a final extension step of 72 °C for 10:00 min. When PCR reactions were unsuccessful, we optimized PCR conditions to include multiple annealing temperatures79: initial denaturation at 95 °C for 10 min; followed by 30 cycles at 95 °C for 1 min, 62 °C for 1 min (decreasing 1 °C every 3 cycles) and 72 °C for 1:30 min; then 30 cycles at 95 °C for 30 s, 55 °C for 30 s and 72 °C for 1 min; and a final extension step of 72 °C for 7 min. PCR products that showed clear bands on agarose gel were purified using the QIAquick PCR Purification Kit (Qiagen) and subsequently sequenced. We prepared 10 μL reactions with the same primers and 3–5.5 μL of purified PCR product. The sequencing reactions were performed using the Big Dye® Terminator v3.1 Cycle Sequencing Kit (Life Technologies, Carlsbad, CA). Generated sequences were assembled, trimmed and edited in Sequencher 4.10.1 (Gene Codes Corporation, Ann Arbor, MI).
Sequence alignment and phylogenetic analyses
We constructed 3 datasets, ITS, LSU and a concatenated SSU+ITS+LSU dataset, to investigate phylogenetic structure within H. virescens. We aligned sequences of each region separately using Muscle v3.780, implemented on the Cipres Science Gateway81. For SSU and LSU, ambiguously aligned regions and uninformative positions were detected and removed using trimAl v1.382 with 60% gap threshold and minimal coverage of 50%. In the ITS dataset, we manually removed the ITS1 (positions 1–525) and ITS2 (687–1067) spacer regions for those sequences other than Hesperomyces because they were too variable to align. The data for each region were concatenated in MEGA783 to create a matrix of 2767 bp with phylogenetic data for 52 isolates.
Maximum likelihood (ML) analyses were run using PAUP on XSEDE84 in Cipres81. Nucleotide substitution models were selected statistically with the help of jModelTest85 by considering the Akaike Information Criterion (AIC). For the ITS dataset, the TVM+G model was selected (lowest -lnL = 3566.8229). For the LSU dataset, the TIM1+G model gave the best scoring tree (-lnL = 3151.3620). ML was inferred for each dataset under the appropriate model; rapid bootstrap (BS) analysis was implemented with 100 replicates.
We performed ML and Bayesian analyses for the SSU+ITS+LSU dataset. For ML, the dataset was divided into three partitions. The best substitution model was selected using jModelTest85 by considering the Akaike Information Criterion (AIC). The following models were selected: TPM2uf (SSU), TVM+G (ITS) and GTR+G (LSU). [We did not use all LSU isolates for the three-gene concatenated dataset, hence the selection of a different substitution model.] ML was inferred under partitioned models using IQ-tree86,87; ultrafast bootstrap (BS) analysis was implemented with 1000 replicates88. Also for Bayesian inference, the three-gene dataset was divided into partitions. Analyses were done with a Markov chain Monte Carlo (MCMC) coalescent approach implemented in BEAST89, with an uncorrelated lognormal relaxed molecular clock allowing for rate variation across the tree. We selected the Birth-Death Incomplete Sampling speciation model90 as tree prior with JC (for SSU), TPM2uf+G (for ITS) and TIM1+G (for LSU) nucleotide substitution models (considering the Bayesian Information Criterion from jModelTest) and a lognormal ucld.mean (mean = 5.0, stdev = 1.0). Five independent runs were performed from a random starting tree for 10 million generations, with a sampling frequency of 1000. Prior settings were entered in BEAUti89 to generate an XML file, which was run using the BEAST on XSEDE tool in Cipres81 (three runs) and locally from the command line (two runs). The resulting log files of the five runs were entered in Tracer v1.6.091 to check trace plots for convergence and effective sample size (ESS). Burn-in was adjusted to achieve ESS values of ≥200 for the majority of sampled parameters. While removing a portion of each run as burn-in, log files and trees files were combined in LogCombiner. TreeAnnotator was used to generate consensus trees with 0% burn-in and to infer the maximum clade credibility tree, with the highest product of individual clade posterior probabilities. Final trees with bootstrap values (BS) and posterior probabilities (pp) were visualized in FigTree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/).
Species delimitation analyses
Morphology-based identification of H. virescens thalli63,64,68,69 may mask multiple species within a geographical context or with strict host specificity. Therefore, we used 3 species delimitation methods to validate species hypotheses28,92,93: Automatic Barcode Gap Discovery (ABGD) and General Mixed Yule Coalescent methods (GMYC), and a Poisson tree processes (PTP) model approach.
ABGD is based on the detection of a “barcode gap,” which is observed when nucleotide divergence among isolates of the same species is smaller than divergence among isolates of different species in a given multiple alignment. Gaps are identified and used to partition (or: split) the data into the maximum number of groups, which represent species hypotheses92. We used the web version of ABGD (at http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html) to identify barcode gaps in the SSU+ITS+LSU dataset. Genetic distances were calculated using both available distance metrics (JC69, K80) with the following parameters: Pmin = 0.001, Pmax = 0.01, steps = 10 and Nb bins = 20. To assess consistency of the species recognized by ABGD, we evaluated results for four gap width values (X): 0.1, 0.5, 1.0 and 1.5.
In the PTP model approach, the number of nucleotide substitutions is directly used to model speciation rate. The underlying assumption is that the number of substitutions between species is significantly higher than the number of substitutions within species28. Compared with GMYC, Zhang and colleagues28 found that PTP performs best especially when the evolutionary distances between species are small. PTP is intended for the delimitation of species in single-gene trees. As a result, we applied this method to both the ITS and LSU phylogenetic reconstructions separately. As input for the PTP model approach, we used phylogenetic trees generated by Bayesian analyses for the two datasets. The MCMC analyses were done under a strict molecular clock, with the Yule speciation tree prior and the appropriate nucleotide substitution model, as selected by the Bayesian Information Criterion from jModelTest 2.1. For the ITS dataset, the TVM2uf+G model was selected (lowest -lnL = 3573.5434). For the LSU dataset, the K80+I model gave the best scoring tree (-lnL = 1672.9035). Two independent runs were performed from a random starting tree for 10 million generations, with a sampling frequency of 1000. The two resulting log files were combined in LogCombiner with 1% burn-in. Consensus trees with 0% burn-in were generated and the maximum clade credibility tree was constructed in TreeAnnotator. We used the bPTP web server (http://species.h-its.org). The “b” in bPTP stands for the Bayesian support values that are added to delimited species. The different parameters were set as default (number of MCMC generations, thinning, burn-in, seed). For both analyses the outgroups were removed from the dataset prior to constructing the phylogenetic tree.
GMYC uses a fully-resolved ultrametric tree93 inferred from a single marker to model processes at the population level (coalescence) and processes at the species level (speciation). As input we used the ITS and LSU maximum clade credibility trees generated for PTP. In addition, we reconstructed a maximum clade credibility tree in BEAST89 using the concatenated SSU+ITS+LSU dataset. For this analysis we removed the outgroups from the dataset, because the inclusion of distantly related species makes it more difficult for GMYC to detect closely related species. The MCMC analysis was done under an uncorrelated lognormal relaxed molecular clock, with the Birth-Death Incomplete Sampling speciation model89 tree prior, the GTR+I+G nucleotide substitution model (as selected under the Bayesian Information Criterion) and a lognormal ucld.mean (mean = 5.0, stdev = 1.0). Two independent runs were performed from a random starting tree for 80 million generations, with a sampling frequency of 8000. The two resulting log files were combined in LogCombiner with 10% burn-in. The maximum clade credibility tree was again constructed in TreeAnnotator. Species were delimited based on this generated ultrametric tree with the GMYC method using R packages rncl94 and splits95.
Data Availability
Voucher specimens of infected ladybirds will be deposited at the Harvard Museum of Comparative Zoology, Cambridge, MA and the Brabant Museum of Nature in Tilburg, The Netherlands. Voucher slides of Hesperomyces virescens are deposited at Farlow Herbarium at Harvard University (FH) (barcodes available in Supplementary Table S1). All generated sequences have been uploaded to GenBank (accession numbers MG745336–MG745358, MG757798–MG757831, MG760581–760611). The following data are available from figshare https://doi.org/10.6084/m9.figshare.c.3944749.v1: sequence alignments generated during this study (in NEXUS format), input XML and output log and trees files from Bayesian analyses, R code used for the GMYC analysis.
References
Simpson, G. G. Principles of animal taxonomy (Columbia University Press, 1961).
de Queiroz, K. The general lineage concept of species, species criteria, and the process of speciation: A conceptual unification and terminological recommendations in Endless forms. Species and speciation (eds Howard, D. J. & Berlocher, S. H.) 57–75 (Oxford University Press, 1998).
de Queiroz, K. Species concepts and species delimitation. Syst. Biol. 56, 879–886 (2007).
Blackwell, M. The Fungi: 1, 2, 3 … 5.1 million species? Am. J. Bot. 98, 426–438 (2011).
Hibbett, D. et al. Sequence-based classification and identification of Fungi. Mycologia 108, 1049–1068 (2016).
Hawksworth, D. L. The fungal dimension of biodiversity: magnitude, significance, and conservation. Mycol. Res. 95, 641–655 (1991).
O’Brien, H. E., Parrent, J. L., Jackson, J. A., Moncalvo, J. M. & Vilgalys, R. Fungal community analysis by large-scale sequencing of environmental samples. Appl. Environ. Microbiol. 71, 5544–5550 (2005).
Taylor, D. L. et al. A first comprehensive census of fungi in soil reveals both hyperdiversity and fine‐scale niche partitioning. Ecol. Monogr. 84, 3–20 (2014).
Pringle, A. et al. Cryptic speciation in the cosmopolitan and clonal human pathogenic fungus Aspergillus fumigatus. Evolution 59, 1886–1899 (2005).
Singh, G. et al. Coalescent-based species delimitation approach uncovers high cryptic diversity in the cosmopolitan lichen-forming fungal genus Protoparmelia (Lecanorales, Ascomycota). Plos One 10, e0124625, https://doi.org/10.1371/journal.pone.0124625 (2015).
Lücking, R. et al. A single microlichen constitutes hundreds of unrecognized species. Proc. Natl. Acad. Sci. 111, 11091–11096 (2014).
Schluter, D. Ecology and the origin of species. Trends Ecol. Evol. 16, 372–380 (2001).
Schluter, D. The ecology of adaptive radiation (Oxford University Press, 2015).
Rundle, H. D. & Nosil, P. Ecological speciation. Ecol. Lett. 8, 336–352 (2005).
Dobzhansky, T. Complete reproductive isolation between two morphologically similar species of Drosophila. Ecology 27, 205–211 (1946).
den Bakker, H. C., Zuccarello, G. C., Kuyper, T. H. & Noordeloos, M. E. Evolution and host specificity in the ectomycorrhizal genus Leccinum. New Phytol. 163, 201–215 (2004).
Araújo, J. P. M., Evans, H. C., Geiser, D. M., Mackay, W. P. & Hughes, D. P. Unravelling the diversity behind the Ophiocordyceps unilateralis (Ophiocordycipitaceae) complex: Three new species of zombie-ant fungi from the Brazilian Amazon. Phytotaxa 220, 224–238 (2015).
Weir, A. & Hammond, P. M. Laboulbeniales on beetles: host utilization patterns and species richness of the parasites. Biodivers. Conserv. 6, 701–719 (1997).
Weir, A. & Blackwell, M. Fungal biotrophic parasites of insects and other arthropods in Insect-fungal associations: Ecology and evolution (eds Vega, F. E. & Blackwell, M.) 119–145 (Oxford University Press, 2005).
Santamaria, S. et al. First Laboulbeniales from harvestmen: the new genus. Opilionomyces. Phytotaxa 305, 285–292 (2017).
Haelewaters, D. & De Kesel, A. De schimmel Hesperomyces virescens, een natuurlijke vijand van lieveheersbeestjes. Entomol. Ber. 77, 106–118 (2017).
Haelewaters, D. et al. Parasites of Harmonia axyridis: current research and perspectives. BioControl 62, 355–371 (2017).
Cottrell, T. E. & Riddick, E. W. Limited transmission of the ectoparasitic fungus Hesperomyces virescens between ladybirds. Psyche 2012, 814378, https://doi.org/10.1155/2012/814378 (2012).
Majerus, M. E. N. Ladybirds. (Harper Collins, 1994).
Ireland, H., Kearns, P. W. E. & Majerus, M. E. N. Interspecific hybridisation in the coccinellids: some observations on an old controversy. Entomol. Rec. J. Var. 98, 181–185 (1986).
Majerus, M. E. N. Interspecific hybridisation in the Coccinellidae. Entomol. Rec. J. Var. 109, 11–23 (1997).
Carstens, B. C., Pelletier, T. A., Reid, N. M. & Satler, J. D. How to fail at species delimitation. Mol. Ecol. 22, 4369–4383 (2013).
Zhang, J., Kapli, P., Pavlidis, P. & Stamatakis, A. A general species delimitation method with applications to phylogenetic placements. Bioinformatics 29, 2869–2876 (2013).
Leavitt, S. D. et al. Who’s getting around? Assessing species diversity and phylogeography in the widely distributed lichen-forming fungal genus Montanelia (Parmeliaceae, Ascomycota). Mol. Phylogenet. Evol. 90, 85–96 (2015).
Kekkonen, M. & Hebert, P. D. DNA barcode‐based delineation of putative species: efficient start for taxonomic workflows. Mol. Ecol. Resour. 14, 706–715 (2014).
Dayrat, B. Towards integrative taxonomy. Biol. J. Linn. Soc. 85, 407–415 (2005).
Zamora, J. C., Calonge, F. D. & Martín, M. P. Integrative taxonomy reveals an unexpected diversity in Geastrum section Geastrum (Geastrales, Basidiomycota). Persoonia 34, 130–165 (2015).
Padial, J. M., Miralles, A., De la Riva, I. & Vences, M. The integrative future of taxonomy. Frontiers Zool. 7, 16, https://doi.org/10.1186/1742-9994-7-16 (2010).
Leavitt, S. D., Moreau, S. D. & Lumbsch, H. T. The dynamic discipline of species delimitation: progress toward effectively recognizing species boundaries in natural populations in Recent advances in lichenology (eds Upreti, D. K., Divakar, P. K., Shukla, V. & Bajapi, R.) 11–44 (Springer, 2015b).
Sousa, J. O. et al. More than one fungus in the pepper pot: Integrative taxonomy unmasks hidden species within Myriostoma coliforme (Geastraceae, Basidiomycota). Plos One 12, e0177873, https://doi.org/10.1371/journal.pone.0177873 (2017).
Schoch, C. L. et al. Nuclear ribosomal internal transcribed spacer (ITS) region as a universal DNA barcode marker for Fungi. Proc. Natl. Acad. Sci. 109, 6241–6246 (2012).
Krüger, M., Krüger, C., Walker, C., Stockinger, H. & Schüβler, A. Phylogenetic reference data for systematics and phylotaxonomy of arbuscular mycorrhizal fungi from phylum to species level. New Phytol. 193, 970–984 (2012).
Samson, R. A. et al. Phylogeny, identification and nomenclature of the genus. Aspergillus. Stud. Mycol. 78, 141–173 (2014).
Crous, P. W., Hawksworth, D. L. & Wingfield, M. J. Identifying and naming plant-pathogenic fungi: past, present, and future. Ann. Rev. Phytopathol. 53, 247–267 (2015).
Yahr, R., Schoch, C. L. & Dentinger, B. T. Scaling up discovery of hidden diversity in fungi: impacts of barcoding approaches. Phil. Trans. R. Soc. B 371, 20150336, https://doi.org/10.1098/rstb.2015.0336 (2016).
Liu, F. et al. Unravelling Colletotrichum species associated with Camellia: employing ApMat and GS loci to resolve species in the C. gloeosporioides complex. Persoonia 35, 63–86 (2015).
Al-Hatmi, A. M. S. et al. Evaluation of two novel barcodes for species recognition of opportunistic pathogens in Fusarium. Fungal Biol. 120, 231–245 (2016).
Rosling, A. et al. Archaeorhizomycetes: unearthing an ancient class of ubiquitous soil fungi. Science 333, 876–879 (2011).
Kovács, G. M., Jankovics, T. & Kiss, L. Variation in the nrDNA ITS sequences of some powdery mildew species: do routine molecular identification procedures hide valuable information? Eur. J. Plant Pathol. 131, 135, https://doi.org/10.1007/s10658-011-9793-3 (2011).
Lindner, D. L. et al. Employing 454 amplicon pyrosequencing to reveal intragenomic divergence in the internal transcribed spacer rDNA region in fungi. Ecol. Evol. 3, 1751–1764 (2013).
den Bakker, H. C., Gravendeel, B. & Kuyper, T. W. An ITS phylogeny of Leccinum and an analysis of the evolution of minisatellite-like sequences within ITS1. Mycologia 96, 102–118 (2004).
Nilsson, R. H., Kristiansson, E., Ryberg, M., Hallenberg, N. & Larsson, K. H. Intraspecific ITS variability in the kingdom Fungi as expressed in the international sequence databases and its implications for molecular species identification. Evol. Bioinform. 4, 193–201 (2008).
Vilgalys, R. & Hester, M. Rapid genetic identification and mapping of enzymatically amplified ribosomal DNA from several Cryptococcus species. J. Bacterial. 172, 4238–4246 (1990).
Miadlikowska, J. & Lutzoni, F. Phylogenetic revision of the genus Peltigera (lichen-forming Ascomycota) based on morphological, chemical, and large subunit nuclear ribosomal DNA data. Int. J. Plant Sci. 161, 925–958 (2000).
Bickford, D. et al. Cryptic species as a window on diversity and conservation. Trends Ecol. Evol. 22, 148–155 (2007).
Hawksworth, D. L. Fungal diversity and its implications for genetic resource collections. Stud. Mycol. 50, 9–18 (2004).
Grünig, C. R., Duo, A., Sieber, T. N. & Holdenrieder, O. Assignment of species rank to six reproductively isolated cryptic species of the Phialocephala fortinii s.l.-Acephala applanata species complex. Mycologia 100, 47–67 (2008).
Stefani, F. O. P., Jones, R. H. & May, T. W. Concordance of seven gene genealogies compared to phenotypic data reveals multiple cryptic species in Australian dermocyboid Cortinarius (Agaricales). Mol. Phylogenet. Evol. 71, 249–260 (2014).
Olariaga, I., Laskibar, X. & Holec, J. Molecular data reveal cryptic speciation within Tricholomopsis rutilans: description of T. pteridicola sp. nov. associated with Pteridium aquilinum. Mycol. Progress 14, 21, https://doi.org/10.1007/s11557-015-1040-4 (2015).
Hong, S. G. & Jung, H. S. Phylogenetic analysis of Ganoderma based on nearly complete mitochondrial small-subunit ribosomal DNA sequences. Mycologia 96, 742–755 (2004).
Van de Putte, K., Nuytinck, J., De Crop, E. & Verbeken, A. Lactifluus volemus inEurope: three species in one–revealed by a multilocus genealogical approach, Bayesian species delimitation and morphology. Fungal Biol. 120, 1–25 (2016).
Adamcik, S. et al. A molecular analysis reveals hidden species diversity within the current concept of Russula maculata (Russulaceae, Basidiomycota). Phytotaxa 270, 71–88 (2016).
Li, Y.-M., Shivas, R. G. & Cai, L. Cryptic diversity in Tranzscheliella spp.(Ustilaginales) is driven by host switches. Sci. Rep. 7, 43549, https://doi.org/10.1038/srep43549 (2017).
Zelditsch, M. L., Swiderski, D. L. & Sheets H. D. Geometric morphometrics for biologists: a primer (Academic Press, 2012).
De Kesel, A. & den Neucker, V. T. Morphological variation in Laboulbenia flagellata (Ascomycetes, Laboulbeniales). Belg. J. Bot. 138, 165–172 (2006).
De Kesel, A. & Haelewaters, D. Laboulbenia slackensis and L. littoralis sp. nov. (Ascomycota, Laboulbeniales), two sibling species as a result of ecological speciation. Mycologia 106, 407–414 (2014).
Goldmann, L., Weir, A. & Rossi, W. Molecular analysis reveals two new dimorphic species of Hesperomyces (Ascomycota, Laboulbeniomycetes) parasitic on the ladybird Coleomegilla maculata (Coleoptera, Coccinellidae). Fungal Biol. 117, 807–813 (2013).
De Kesel, A. Hesperomyces (Laboulbeniales) and coccinellid hosts. Sterbeeckia 30, 32–37 (2011).
Thaxter, R. Supplementary note on North American Laboulbeniaceae. Proc. Am. Acad. Arts Sci. 25, 261–270 (1891).
Maddison, W. P. Gene trees in species trees. Syst. Biol. 46, 523–536 (1997).
Riddick, E. W. & Cottrell, T. E. Is the prevalence and intensity of the ectoparasitic fungus Hesperomyces virescens related to the abundance of entomophagous coccinellids. Bull. Insectol. 63, 71–78 (2010).
Benjamin, R. K. Introduction and supplement to Roland Thaxter’s contribution towards a monograph of the Laboulbeniaceae. Bibl. Mycol. 80, 1–155 (1971).
Thaxter, R. Contribution towards a monograph of the Laboulbeniaceae. Mem. Am. Acad. Arts Sci. 12, 187–429 (1896).
Santamaría, S. L. II Acompsomyces–Ilyomyces. Flora Mycol. Iberica 5, 1–344 (2003).
Bates, D., Maechler., M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).
Akaike, H. New look at statistical-model identification. IEEE Trans. Autom. Control 19, 716–723 (1974).
Kassambara, A. factoextra: Visualization of the outputs of a multivariate analysis. R package version 1.0.1 https://cran.r-project.org/web/packages/factoextra/index.html (2015).
Haelewaters, D. et al. Bringing Laboulbeniales to the 21st century: enhanced techniques for extraction and PCR amplification of DNA from minute ectoparasitic fungi. IMA Fungus 6, 363–372 (2015).
Haelewaters, D. Studies of the Laboulbeniomycetes: Diversity, evolution, and speciation. Ph.D. dissertation thesis, Harvard University (2018).
White, T. J., Bruns, T. D., Lee, S. B. & Taylor, J. W. Analysis of phylogenetic relationships by amplification and direct sequencing of ribosomal RNA genes in PCR Protocols: a guide to methods and applications (eds Innis, M. A., Gelfand, D. H., Sninsky, J. J., White, T. J.) 315–322 (Academic Press, 1990).
Landvik, S., Egger, K. N. & Schumacher, T. Towards a subordinal classification of the Pezizales. Nord. J. Bot. 17, 403–418 (1997).
Gardes, M. & Bruns, T. D. ITS Primers with enhanced specificity for Basidiomycetes – application to the identification of mycorrhizae and rusts. Mol. Ecol. 2, 113–118 (1993).
Larena, I., Salazar, O., González, V., Julián, M. C. & Rubio, V. Design of a primer for ribosomal DNA internal transcribed spacer with enhanced specificity for ascomycetes. J. Biotechnol. 75, 187–194 (1999).
Don, R. H., Cox, P. T., Wainwright, B. J., Baker, K. & Mattick, J. S. ‘Touchdown’ PCR to circumvent spurious priming during gene amplification. Nucleic Acids Res. 19, 4008, https://doi.org/10.1093/nar/19.14.4008 (1991).
Edgar, R. C. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792–1797 (2004).
Miller, M. A., Pfeiffer, W. & Schwartz, T. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. Proc. Gateway Comp. Environ. Workshop 14 Nov. 2010, 1–8 (2010).
Capella-Gutiérrez, S., Silla-Martínez, J. M. & Gabaldón, T. TrimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics 25, 1972–1973 (2009).
Kumar, S., Stecher, G. & Tamura, K. MEGA7: Molecular Evolutionary Genetics Analysis version 7.0 for bigger datasets. Mol. Biol. Evol. 33, 1870–1874 (2016).
Swofford, D. L. PAUP: Phylogenetic Analysis Using Parsimony, Version 3.1 (1991).
Darriba, D., Taboada, G. L., Doallo, R. & Posada, D. jModelTest 2: more models, new heuristics and parallel computing. Nature Methods 9, 772, https://doi.org/10.1038/nmeth.2109 (2012).
Nguyen, L.-T., Schmidt, H. A., von Haeseler, A. & Minh, B. Q. IQ-TREE: A fast and effective stochastic algorithm for estimating maximum likelihood phylogenies. Mol. Biol. Evol. 32, 268–274 (2015).
Chernomor, O., von Haeseler, A. & Minh, B. Q. Terrace aware data structure for phylogenomic inference from supermatrices. Syst. Biol. 65, 997–1008 (2016).
Hoang, D. T., Chernomor, O., von Haeseler, A., Minh, B. Q. & Vinh, L. S. UFBoot2: Improving the ultrafast bootstrap approximation. Mol. Biol. Evol. 35, 518–522 (2017).
Drummond, A. J., Suchard, M. A., Xie, D. & Rambaut, A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol. Biol. Evol. 29, 1969–1973 (2012).
Stadler, T. On incomplete sampling under birth–death models and connections to the sampling-based coalescent. J. Theor. Biol. 261, 58–66 (2009).
Rambaut, A., Suchard, M. A., Xie, D. & Drummond, A. J. Tracer v1.6. http://tree.bio.ed.ac.uk/software/tracer/ (2014).
Puillandre, N., Lambert, A., Brouillet, S. & Achaz, G. ABGD, Automatic Barcode Gap Discovery for primary species delimitation. Mol. Ecol. 21, 1864–1877 (2012).
Pons, J. et al. Sequence-based species delimitation for the DNA taxonomy of undescribed insects. Syst. Biol. 55, 595–609 (2006).
Michonneau, F., Bolker, B., Holder, M., Lewis, P. & O’Meara, B. rncl: an interface to the nexus class library. R package version 0.6.0 http://CRAN.R-project.org/package = rncl (2015).
Ezard, T., Fujisawa, T. & Barraclough, T. G. splits: SPecies’ LImits by Threshold Statistics. R package version 1.0-14/r31 http://R-Forge.R-project.org/projects/splits/ (2009).
Acknowledgements
This work was done with the help and support of many researchers, collaborators and citizen scientists. For reading and commenting on previous versions of the manuscript: Ph.D. committee members Meredith Blackwell, Charles C. Davis, Brian D. Farrell, Annemieke Verbeken; OEB 210′s Andrew D. Richardson, Ambika Kamath and fellow students M. Holly Elmore, Morgan E. Furze, Patrick S. Gorring, Sofia R. Prado-Irwin, Anna Waldeck, Qian Wan. For providing infected insects: Susanna Clusella-Trullas, Ted E. Cottrell, Emily A. Kemp, Ingrid A. Minnaar, Katie M. Murray, Gabriele Rondoni, Tove Steenberg, Damon Tighe, Simon Tragust, Maura Tyrrell, Paul S. van Wielink. For assistance in the molecular lab: Michał Gorczak, Ina Kodra, Hamidah Mahmud, Julie Park. For assistance with fieldwork in Panama: Juan A. Bernal Vega, Maria Fuentes, Leila González, Lauren A. Meckler, Tomás A. Ríos Gonzáles, Kirk A. Silas, Sarah J.C. Verhaeghen, Rosa V. Villarreal-Saucedo. For providing permission to use their photographs of ladybirds: Sally Adam, Roberto Güller, Pavel Kirillov, Gilles San Martin, Katja Schulz, Damon Tighe. Funding was provided by the following sources (to D.H.): David Rockefeller Center for Latin American Studies (2015 and 2017 Summer Research Travel Grants, 2016 Term-Time Research Fellowship), Smithsonian Tropical Research Institute (2017 Short-Term Research Fellowship), Torrey Botanical Society (2015 Graduate Student Research Fellowship). Panama’s Secretaria Nacional de Ciencia, Tecnología e Innovacion (SENACYT) provided funding for the international workshop “Hongos Asociados a Insectos: Un Aporte al Conocimiento de la Biodiversidad” in July 2016 (No. APY-GC-2015-57, to Rosa V. Villarreal-Saucedo, Universidad Autónoma de Chiriquí), which allowed D.H. to travel to Panama to conduct additional fieldwork.
Author information
Authors and Affiliations
Contributions
D.H. initiated the project, generated molecular sequence data, performed molecular phylogenetic analyses, prepared all figures and wrote the manuscript with input from co-authors A.D.K. and D.H.P. D.H.P. provided expertise at all stages of research.
Corresponding author
Ethics declarations
Competing Interests
The authors declare no competing interests.
Additional information
Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Haelewaters, D., De Kesel, A. & Pfister, D.H. Integrative taxonomy reveals hidden species within a common fungal parasite of ladybirds. Sci Rep 8, 15966 (2018). https://doi.org/10.1038/s41598-018-34319-5
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41598-018-34319-5
Keywords
This article is cited by
-
Singleton-based species names and fungal rarity: Does the number really matter?
IMA Fungus (2024)
-
New insights into the DNA extraction and PCR amplification of minute ascomycetes in the genus Laboulbenia (Pezizomycotina, Laboulbeniales)
IMA Fungus (2024)
-
Beetlehangers.org: harmonizing host–parasite records of Harmonia axyridis and Hesperomyces harmoniae
Arthropod-Plant Interactions (2024)
-
Parasitism of ladybirds by Hesperomyces in the Canary Islands
Arthropod-Plant Interactions (2024)
-
Integrated evidence sheds light on the taxonomy of the widespread Tantilla melanocephala species complex (Serpentes: Colubridae) and indicates the existence of a new species from southern South America
Organisms Diversity & Evolution (2024)
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.