Coalescent-Based Species Delimitation Approach Uncovers High Cryptic Diversity in the Cosmopolitan Lichen-Forming Fungal Genus Protoparmelia (Lecanorales, Ascomycota)

Species recognition in lichen-forming fungi has been a challenge because of unsettled species concepts, few taxonomically relevant traits, and limitations of traditionally used morphological and chemical characters for identifying closely related species. Here we analyze species diversity in the cosmopolitan genus Protoparmelia s.l. The ~25 described species in this group occur across diverse habitats from the boreal -arctic/alpine to the tropics, but their relationship to each other remains unexplored. In this study, we inferred the phylogeny of 18 species currently assigned to this genus based on 160 specimens and six markers: mtSSU, nuLSU, ITS, RPB1, MCM7, and TSR1. We assessed the circumscription of species-level lineages in Protoparmelia s. str. using two coalescent-based species delimitation methods – BP&P and spedeSTEM. Our results suggest the presence of a tropical and an extra-tropical lineage, and eleven previously unrecognized distinct species-level lineages in Protoparmelia s. str. Several cryptic lineages were discovered as compared to phenotype-based species delimitation. Many of the putative species are supported by geographic evidence.


Introduction
Lichens are symbiotic organisms consisting of a fungal partner (mycobiont), one or more photosynthetic partners (photobionts; [1]), and diverse bacterial communities [2]. Lichens contribute to ecosystem functioning by nutrient recycling [3], weathering rocks, preventing soil erosion, and acting as pioneer species in barren areas. They inhabit diverse ecosystems from the arctic to the tropics and commonly form an integral part of terrestrial biodiversity [4]. Lichens are preferred model systems for ecological, evolutionary, phylogeographic and population genetic studies of symbiotic associations on the account of their wide, often cosmopolitan, distribution, intriguing eco-physiological interdependence and co-evolutionary and adaptive strategies [5]. Almost one fifth of all known fungi and half of all ascomycetes are lichenized, consisting of approximately 28,000 species worldwide [6,7]. However, studies suggest that the estimate of existing lichen diversity might represent only 50-60% of the real diversity [8,9], as current species recognition in lichen-forming fungi appears to vastly underestimate the true number of species. According to Galloway [8], the number of known taxa in different genera has increased from 20% (Parmelia sensu stricto, [10]) to 86% in the New World Oropogon [11]. Recent molecular studies have demonstrated the presence of many distinct lineages subsumed under a single species name (e.g., [12][13][14][15]). In the basidiolichen fungus Dictyonema glabratum a single taxon was found to be composed of at least 126 species [9], thus showing a tremendous amount of unexplored diversity in lichen-forming fungi.
Species recognition in lichen-forming fungi has been a challenge because of i) the few taxonomically relevant characters (reviewed by [16,17]), ii) unsettled species concepts [18,19], iii) and unexplored regions containing high levels of diversity, especially in the tropics [20,21]. Morphological and chemical characters that have commonly been used to circumscribe species may not be useful for identifying closely-related species and often fail to accurately characterize species-level diversity [2,19,22,23]. Accurate species delimitation may be obscured by cryptic speciation [24,25], incongruence between morphology and molecular data [26,27], or incongruence between gene trees and species trees [28]. Moreover, morphological and chemical variations may constitute morpho-or chemotypes of the same species with no molecular differentiation, thus blurring our understanding of species boundaries [29,30]. The implementation of molecular techniques and availability of markers for amplifying phylogenetically informative loci have provided great insights into otherwise unrecognized species complexes. Improved species recognition has important implications for understanding diversity, ecological and biogeographical patterns, factors promoting diversification, and for devising better conservation policies [31].
Different studies have utilized varied combinations of the available techniques for unraveling hidden diversity. For example, Harrington and Near [32] used STEM [33] to explore the independent evolutionary lineages within snubnose darters (Etheostoma simoterum species complex). Giarla et al. [15] used two coalescent-based approaches (BP&P and spedeSTEM) for delimiting species in Andean mouse opossums (Thylamys spp) using three nuclear loci and found three additional lineages than previously recognized. Leavitt et al. [34,35] used Bayesian population clustering, genealogical concordance, Bayesian species delimitation, and a DNA barcode approach to support the presence of five previously unrecognized species in the lichen-forming fungus Rhizoplaca melanophthalma species-complex (Lecanoraceae). Parnmen et al. [36] used a 4-locus phylogenetic approach, combined with GMYC [37,38] and STEM [33] and found at least 12 species in the Cladia aggregata complex. Mounting evidence continues to support the perspective that traditional phenotype-based species boundaries fail to adequately characterized species-level diversity in many lichen-forming fungi (reviewed in [22]).
We implemented a molecular approach for species recognition in the cosmopolitan lichenforming genus Protoparmelia s. str., combining phylogenetic trees and coalescent-based species delimitation methods. The phylogenetic relationships of the heterogeneous genus Protoparmelia have been a matter of debate. Morphological and anatomical characters of this genus show similarity to both Lecanoraceae and Parmeliaceae. Protoparmelia was initially placed in Lecanoraceae because it includes crustose lichens, with one-celled hyaline ascospores and Lecanora-type ascus [39,40]. Later, secondary metabolite profiles showing the presence of lobaric acid brought into question its placement in Lecanoraceae [41]. Studies on the ascoma ontogeny [42,43] further showed the presence of a typical character of Parmeliaceae in Protoparmelia, i.e. cupular exciple, a cup-shaped structure below the hymenium [44]. DNA sequence-based studies suggested Protoparmelia to be the sister-group to Parmeliaceae [45][46][47]. Tropical species of Protoparmelia with multispored asci were previously placed in the genus Maronina [48]. The authors indicated a close relationship of Protoparmelia and Maronina on the basis of similar ascus types, and suggested the former to be the multi-spore derivative of Protoparmelia. Subsequently, Papong et al. [49] proposed the inclusion of Maronina in Protoparmelia based on molecular data. However, the tropical clade differs from other species in Protoparmelia in being predominantly corticolous, having alectoronic acid as a major compound, and containing many isidiate or sorediate species, whereas most species in the traditional circumscription of Protoparmelia are saxicolous and occur in boreal-arctic/alpine and temperate regions.
Protoparmelia s.l. offers an interesting study system for a variety of reasons. This genus is morphologically and chemically heterogeneous [43,50], and in a previous study [47], we showed that Protoparmelia s.l. is polyphyletic. In addition, the relationships of most taxa to each other remain largely unexplored. Members of this genus inhabit ecologically diverse habitats, such as boreal-arctic/alpine, temperate, Mediterranean, subtropical, and tropical regions and also vary greatly in their distribution range with some species being cosmopolitan (e.g., P. badia, P. memnonia), whereas other, mainly tropical species being locally restricted (e.g., P. orientalis, P. multifera). Furthermore, congeners occur on various substrates, with some species growing on bark or decorticated wood, and others on rocks. Protoparmelia species exhibit varied life styles. For example, some species are lichenicolous and parasitize other lichenforming fungi during early parts of their life cycle [50]. Sexual reproduction is common in some species (P. badia and P. orientalis), whereas others propagate mainly via asexual propagules (P. isidiata, P. corallifera and P. capitata) with or without any sexual reproduction.
The heterogeneity of characters makes Protoparmelia s.l. [51] an interesting candidate for testing species delimitation scenarios using multi-locus DNA sequence data. Protoparmelia s. str. [47] although being a small genus, is sister to the largest family of lichen-forming fungi, i.e., Parmeliaceae [45][46][47], consisting of approximately 2,800 species distributed in 80 genera [52,53]. Resolving relationships of Protoparmelia s. str. may contribute to understanding character evolution in an important clade of lichen-forming fungi. The aims of the current study are two-fold: 1) exploring the phylogenetic relationships of Protoparmelia s.l. species by constructing a multi-locus phylogeny, and 2) assessing the circumscription of lineages in Protoparmelia s. str. based on multi-locus species-tree inference and coalescent approaches.

Taxon sampling
This study includes a total of 160 samples of Protoparmelia s.l. from 18 currently described species. About 70% of the total described species were included in this study. Additionally, three unidentified species, most likely new to science, were also included in the study. We selected 73 taxa from reportedly close relatives of Protoparmelia s.l. [45,47], namely Parmeliaceae (40 taxa), Lecanoraceae (4 taxa), Gypsoplacaceae ([54]; 2 taxa), Miriquidica group (12 taxa), and Ramboldia (10 taxa) to infer the relationship of Protoparmelia s.l. with other taxa within related groups within Lecanorales. Cladoniaceae (5 taxa) were selected as outgroup. Details of the study material and GenBank accession numbers are given in S1 Table. DNA amplification and sequencing Genomic DNA was extracted from lichen thalli using the CTAB method [55]. PCR amplification was performed using general, previously published primers for RPB1, TSR1, MCM7, nuLSU, mtSSU and ITS (Table 1). For some species of Protoparmelia s.l. and Miriquidica group specific primers were designed (Table 1). PCR reactions were carried out in a volume of 25 μl. Each reaction mix contained 2.5 μl buffer, 0.13 μl (0.65 U) Ex Taq polymerase, 1.0 μl dNTP mix (2.5 mM each), 1.0 μl each (10 mM) of the primer set, ca. 20 ng of template, and 16 μl H 2 O. Reactions were performed with the following cycling conditions: initial denaturation at 95°C for 4 min, followed by 35 cycles of 95°C for 30 s, 50°C for 40 s, 72°C for 1 min, and final elongation at 72°C for 5 min. PCR products were checked for amplification on 1% agarose gels. Bands of expected size were extracted using the peqGOLD Gel Extraction Kit. All PCR products were labeled with the Big Dye Terminator v3. For each locus, consensus sequences were assembled separately and aligned using MAFFT [56] as implemented in Geneious v5.4 [57], followed by manual editing. Gaps were treated as missing data and ambiguously aligned nucleotides were excluded.

Phylogenetic analyses
Model selection. Model selection was performed to find the best-fitting model for each data set. We used the Corrected Akaike Information Criterion (AICc) [58] as implemented in jModelTest v2.1.1 [59].
Congruence among loci. To test the level of congruence among loci, we used the Congruence Among Distance Matrices test (CADM, [60]), as implemented in the package ape in R. The null hypothesis assumes that all tested phylogenetic trees are completely incongruent. Incongruence here refers to phylogenetic trees with different topologies among loci, which suggests completely distinct evolutionary histories. The level of congruence ranges from 0 to 1. In addition, maximum likelihood (ML) analyses were performed individually on each locus with RAxML-HPC BlackBox v8.1.11 [61] on the Cipres Science gateway [62] using the default GTR + G model with 1,000 bootstrap (BS) replicates. Conflicts were considered significant if individuals grouped in a clade with ! 70% BS support in one data set, but in a different clade with high support in another locus.
Phylogeny of Protoparmelia s.l. Since no supported conflicts were observed in single locus trees and CADM analysis rejected the hypothesis of incongruence among loci, data sets were concatenated (see Results). The maximum likelihood search was performed on the concatenated 6-locus data set including all the relatives of Protoparmelia s.l. with RAxML-HPC BlackBox v8.1.11 [61] on the Cipres Scientific gateway [62]. Only those taxa for which the sequence information was available for at least three loci were included in the concatenated data set. The default GTR + G model was used as the substitution model and data was partitioned according to the different genes. RPB1, TSR1 and MCM7 data were also partitioned by codon position. Bayesian inference was performed using the best fitting model as inferred by jModelTest, for the single as well as concatenated data sets as implemented in MrBayes v3.2.1 [63,64] on the Cipres Scientific gateway [62]. Two parallel MCMCMC runs were performed each using four chains and 5,000,000 generations, sampling trees every 100th generation. A 50% majority rule consensus tree was generated from the combined sampled trees of both runs after discarding the first 25% as burn-in (12,500 trees, likelihoods below stationary level).
Phylogeny of Protoparmelia s. str. Maximum likelihood analysis was performed individually on each locus of Protoparmelia s. str. (excluding Lecanoraceae, Parmeliaceae, Miriquidica group and Ramboldia clades), with RAxML-HPC BlackBox v8.1.11 [61] on the Cipres Science gateway [62], using the default GTR + G model, with 1,000 BS replicates. Gypsoplacaceae was used as outgroup. Only taxa for which sequence information was available for at least three loci were included in the concatenated data set. The default GTR + G model was used as the substitution model and the data was partitioned according to the different genes. For RPB1, TSR1 and MCM7 data were also partitioned by codon position. Since no supported conflicts were observed in single locus trees and CADM analysis rejected the hypothesis of incongruence among loci, data sets were concatenated. Maximum likelihood search was then performed on the concatenated 6-locus data set using RAxML-HPC BlackBox v8.1.11 [61] on the Cipres Scientific gateway v3.3 [62].
We performed jModelTest for each locus on the reduced data set to select the best locusspecific models of evolution.
Bayesian inference was performed using the best fitting model as suggested by jModelTest, for the single and concatenated data sets separately as implemented in MrBayes v3.2.1 [63,64] on the Cipres Scientific gateway [62]. Two parallel MCMCMC runs were performed each using four chains and 5,000,000 generations, sampling trees every 100th generation. A 50% majority rule consensus tree was generated from the combined sampled trees of both runs after discarding the first 25% as burn-in (12,500 trees). Ã BEAST as implemented in BEAST v2.1 [65] was used to estimate the species tree for BP&P [66]. We used a Birth-Death process and gamma-distributed population sizes for the species tree prior and a pairwise linear population size model with a constant root. Ã BEAST incorporates the coalescent process and the uncertainty associated with gene trees and nucleotide substitution model parameters and estimates the species tree directly from the sequence data. For each locus, the closest model to the best-suggested model from jModelTest under the AICc criterion was selected as the best substitution model for Ã BEAST. Two independent Markov chain Monte Carlo (MCMC) analyses were performed for a total of 100,000,000 generations, sampling every 5,000 steps. Default values were used for the remaining priors. Convergence of the runs to the same posterior distribution and the adequacy of sampling (using the Effective Sample Size [ESS] diagnostic) were assessed with Tracer v1.4 [67]. After removing the first 20% of the samples as burn-in, all runs were combined to generate posterior probabilities of nodes from the sampled trees using TreeAnnotator v1.7.4 [68]. The species tree produced by Ã BEAST was subsequently used for inferring speciation probabilities by BP&P [66].
Species delimitation in Protoparmelia s. str.
For testing the species boundaries in Protoparmelia s. str. [47], currently accepted taxa were taken as putative species (12 described species). In addition well-supported (BS ! 70%, PP ! 0.94) monophyletic clades from ML and Bayesian phylogenies were taken as putative species, resulting in a 25-species scenario (Figs 1 and 2). The marginal posterior probability of 25-species scenario suggested by molecular data was estimated using the program BP&P v3 [66]. BP&P utilizes reversible-jump Bayesian Markov chain Monte Carlo (MCMC) algorithms for analyzing phylogenetic data from multiple loci to generate speciation probabilities of assigned species. It takes into account uncertainties due to unknown gene trees and ancestral coalescent processes. This method accommodates the species phylogeny as well as incomplete lineage sorting due to ancestral polymorphism. Species tree from Ã BEAST was used to infer the speciation probabilities by BP&P. BP&P v3 incorporates nearest-neighbor interchange (NNI) algorithm allowing changes in the species tree topology, eliminating the need for a fixed user-specified guide tree [66]. BP&P gives the posterior probability of each delimited species and the posterior probability for the number of delimited species. A gamma prior G (1, 10), with mean 1/10 = 0.1 (one difference per 10 bp) was used on the population size parameters (s). The age of the root in the species tree (τ 0 ) was assigned the gamma prior G (2, 2000) which means 0.1% of sequence divergence, while the other divergence time parameters were assigned the Dirichlet prior [66]. Each analysis was run twice to confirm consistency between runs.
We also used spedeSTEM for calculating probabilities of the species scenario. SpedeSTEM [69] is based on the multilocus species-tree method STEM [33]. It assumes all putative species as separate lineages and estimates gene trees in PAUP Ã [70]. It then calculates the likelihood for alternative species trees in various permutations and combinations of subpopulations by collapsing two or more species into a single lineage using previously estimated gene trees. Species boundaries are then compared using Akaike information criteria and gives probabilities of different species scenarios. We used θ = 0.05 and each analysis included 500 replicates. We tested all 25 possible permutations for clustering within taxonomic species.

Model test
Bayesian analysis on the complete data set and the reduced Protoparmelia s. str. data set was performed using the best fitting model for each locus in the concatenated sequence as shown in Table 2.  For Ã BEAST, the first available best fitting model for each locus in the concatenated data set, from the models suggested by jModelTest v2.1.1 were the following: RPB1: GTR, TSR1: HKY, MCM7: HKY, nuLSU: GTR, mtSSU: HKY, and ITS: GTR.

Congruence among loci
CADM results showed no significant incongruence among loci, thus allowing concatenation. The null hypothesis of complete incongruence among loci was rejected for both complete (W = 0.75; p<0.0001) and reduced (W = 0.84; p<0.0001) data sets.

Phylogeny of Protoparmelia
Protoparmelia s.l. Nuclear and mitochondrial gene partitions supported the same overall topology. The concatenated six-locus data set contained 233 specimens. Gene partitions had the following lengths: 696 bp for RPB1, 756 for TSR1, 655 bp for MCM7, 1064 bp for nuLSU rDNA, 834 bp for mtSSU and 807 bp for ITS rDNA. The total length of the concatenated alignment was 4812 bp (dryad doi:10.5061/dryad.0q515). The ML tree for the concatenated data set is presented in Fig 1. Nodes with BS ! 70% and Bayesian posterior probability (PP) !0.94 were considered as supported.
Protoparmelia s. str. The concatenated six-locus data set contained 138 specimens, including two taxa from outgroup Gypsoplacaceae. Gene partitions had the following lengths: 696 bp for RPB1, 754 for TSR1, 672 bp for MCM7, 972 bp for nuLSU rDNA, 839 bp for mtSSU and 787 bp for ITS rDNA. The total length of the concatenated alignment was 4720 bp. Most species as currently circumscribed were monophyletic, except P. isidiata, which formed three independent lineages within the tropical clade (P. isidiata A-C, D and E), and the cosmopolitan species P. badia, which contained multiple supported lineages and formed a species complex with P. memnonia (Fig 2). We found evidence for cryptic species-level diversity in the nominal taxa P. badia, P. montagnei, and P. isidiata (clade P. isidiata A-E). Cryptic diversity corresponded to biogeographic patterns in P. isidiata (clades A-C representing North America, South America and Asia, respectively). Within P. badia, the largest lineage (clade P. badia A) was cosmopolitan, whereas the other supported lineages had a Mediterranean, or Iberian distribution (Fig 2).

Species delimitation in Protoparmelia s. str.
We treated terminal clades supported by ! 70% BS and ! 0.94 PP (Figs 1 and 2) as putative species for species delimitation analyses. This resulted in a 25-species scenario for Protoparmelia s. str., in contrast to the current 12-species scenario for Protoparmelia s. str., based on morphological and chemical characters. The 25-species scenario in Protoparmelia s. str. was then investigated for species delimitation using BP&P and spedeSTEM. BP&P supported the presence of 23 species with highest probability (PP = 0.41127). Posterior probability of each delimited species is given in  Table 3). Sixteen putative species (P. memnonia, P. hypotremella, P. oleagina, P. montagnei C, P. orientalis, P. multifera, P. pulchra, P. capitata, P. corallifera, P. sp. KE, P. sp. ZA and the five cryptic isidiate lineages in P. isidiata) were supported as separate lineages by both BP&P and spedeSTEM (Table 4), therefore we suggest these clades to be evolutionary independent. We found conflicting speciation scenarios for P. ochrococca A & B, P. badia A, B1, B2, & C and P. montagnei A & B by the two species delimitation approaches (Fig 2).

Discussion
The genus Protoparmelia is more diverse than the traditional taxonomy suggests. This diversity comprises several previously undescribed species, and cryptic lineages within currently accepted species. Most species of Protoparmelia belong to Protoparmelia s. str., consisting of a tropical and an extra-tropical clade. The tropical clade includes several taxa having multispored asci, which were formerly classified in the genus Maronina [48,49,71,72]. All of its members, except the undescribed South African and Keniyan species, are corticolous. Most members of the tropical clade reproduce vegetatively, although a limited number of species propagate predominantly via sexual reproduction. All supported genetic species-level lineages in the tropical clade are congruent with biogeographic origin of the specimens. Evolutionary rates, i.e. rates of base substitutions in an evolutionary lineage over time, appeared to be accelerated in the tropical clade. This phenomenon has also been previously observed in tropical lichens, and attributed to shorter generation times, higher metabolic rates, continuous physiological activity of a poikilohydric organism in a moist environment, and lack of sexuality [73]. The extra-tropical clade   contains mostly saxicolous taxa, most of which reproduce sexually. Within this group, while some species show restricted distribution, some other have wide geographic distributions, such as the cosmopolitan P. 'badia A' and P. hypotremella which occurs in Europe and North America. Five previously described species and one species putatively new to science group with members of the genus Miriquidica. In contrast to members of Protoparmelia s. str., which produce lobaric or alectoronic acids, these taxa synthesize norstictic acid as major secondary metabolite. Many of these species parasitize other lichens during at least parts of their life cycle [50], a lifestyle not known from members of Protoparmelia s. str. Close affiliations between Miriquidica and Protoparmelia based on shared morphological characteristics have been suggested before [74,75], and a recent molecular study confirmed the close relationship of the P. atriseda-group and Miriquidica [47]. A revision of the genus Miriquidica based on molecular data is currently under way by our colleagues (Timdal, pers. comm.).

Speciation analyses and cryptic diversity
We validated the 25-species scenario for Protoparmelia s. str., which was based on the previously defined species and a few new clades suggested by molecular data (phylogenetic species concept). Based on our sampling, this study largely supported traditionally circumscribed Protoparmelia s. str. species as distinct lineages. However, exceptions included P. isidiata, an asexual tropical species, and P. badia, a sexually reproducing, boreal-arctic/alpine cosmopolitan species. The former was found to be polyphyletic and separated into three distinct lineages, while the later was paraphyletic and formed a species complex with P. memnonia (Figs 1 and 2). The combined use of species-tree topology and coalescent methods revealed the presence of several cryptic lineages in Protoparmelia s. str. This is in concordance with other studies in which molecular markers in combination with statistical tools revealed many genetically distinct lineages hidden under a single taxon [9,36,[76][77][78]. Studies suggest that cosmopolitan species such as P. badia may reveal high cryptic diversity [79,80], which may or may not correlate to geography. In our study we found that the cosmopolitan P. badia as currently delimited consists of at least four independent evolutionary lineages. Among these newly recognized lineages only P. badia A turns out to be cosmopolitan, inhabiting boreal-arctic/alpine habitats in North America, Europe, New Zealand and Australia. The other lineages of P. badia (P. badia B1, B2 and C) have a more limited distribution, having been collected so far on siliceous substrates in Spain and Italy. Cryptic lineages within P. isidiata (clades A-C) also correspond to broad biogeographic patterns, while lineages identified within P. montagnei co-occur in the Mediterranean region (Fig 2). Thus, geographic evidence supports species delimitation suggested by coalescent-based speciation analyses in most cases. However, current sampling in many lineages is relatively sparse and does not allow conclusions about finer-scale biogeographic patterns, such as endemism. It remains to be seen whether sympatrically-occurring cryptic lineages identified in this study are supported by additional, previously overlooked morphological or chemical characteristics. We have preliminary evidence that the currently recognized P. montagnei chemotypes [81] correspond to the three molecular clades and may thus indeed represent closely related, but separate species.
Conflicts between different methodological approaches to species delimitation are common [13,15,78,82]. In general we follow the approach of adopting the speciation scenario that is supported by both the analyses, in our case 16 species [83]. For some clades, i.e. P. ochrococca A & B, P. badia A, B1, B2 & C, P. montagnei A & B, the most likely speciation scenario given by spedeS-TEM deviates from BP&P, and contradicts supported branching patterns in the phylogeny (Figs  2 and 3). For P. badia A, B1, B2 & C, P. montagnei A & B phylogenetic tree and BP&P supported these clades to be evolutionary independent, whereas spedeSTEM suggested them to be a single species. For P. ochrococca A & B phylogenetic tree and spedeSTEM supported these clades to be evolutionary independent, whereas BP&P suggested them to be a single species. Recent studies indicated that spedeSTEM may be less accurate than other species delimitation methods in cases of recent speciation events [84]. For the clades supported by BP&P and not spedeSTEM, we preferred BP&P results as BP&P has been shown to perform well even when putative species were modeled to have diverged from one another only very recently [84]. In addition, BP&P has been shown to outperform other coalescent-based species delimitation approaches especially when using multi-locus DNA sequence data and a modest number of individuals per species [69,83]. Previously the reliability of BP&P has been suggested to be dependent on the accuracy of the user-provided guide tree. However, in the latest version of BP&P the authors addressed this issue and applied the NNI algorithm, which allows flexibility in the species tree. Moreover BP&P is suggested to be conservative in delimiting species, with high probability to be a reliable indicator of evolutionary independence of the lineages [66]. Therefore in case of conflicts we considered BP&P to be more accurate and suggested the lineages supported by BP&P as distinct species.
Our analyses suggest that the sampled specimens of the tropical Protoparmelia s. str. group belong to five distinct species. Two sexually reproducing (apotheciate) species, P. multifera and P. orientalis, traditionally distinguished by having different minor secondary metabolites [49] were supported as different species and were not sister to each other. In fact, the sexually reproducing species P. pulchra was sister to P. orientalis. In addition, we found four distinct asexually reproducing (isidiate) species of Protoparmelia s. str. Two of these species (P. 'isidiata D' and P. 'isidiata E') occur sympatrically in Australia. Several studies have shown the occurrence of phylogenetically unrelated but morphologically similar lineages thus indicating the presence of high hidden diversity in lichen-forming fungi [25,27,34,85,86].

Conclusions
Our analyses support the presence of 23 distinct lineages in Protoparmelia s. str. in contrast to 12 currently delimited species, revealing much more diversity than currently suggested for this genus. Our study shows that the sister group of the largest family of lichen-forming fungi may harbor a considerable amount of cryptic lineages which can be identified using molecular data. These data highlight the presence of substantial phylogenetic diversity especially in the tropics, and the need for careful re-evaluation of morphological and chemical characters in the group.
Supporting Information S1 Table. Specimens used in this study including voucher information and GenBank accession numbers. (XLSX)