Hidden Genetic Diversity in an Asexually Reproducing Lichen Forming Fungal Group

Asexual species with vegetative propagation of both symbiont partners (soredia) in lichens may harbor lower species diversity because they may indeed represent evolutionary dead ends or clones. In this study we aim to critically examine species boundaries in the sorediate lichen forming fungi Parmotrema reticulatum–Parmotrema pseudoreticulatum complex applying coalescent-based approaches and other recently developed DNA-based methods. To this end, we gathered 180 samples from Africa, Asia, Australasia, Europe, North and South America and generated sequences of internal transcribed spacer of nuclear ribosomal DNA (ITS) and DNA replication licensing factor MCM7 (MCM7). The dataset was analysed using different approaches such as traditional phylogeny–maximum likelihood and Bayesian–genetic distances, automatic barcode gap discovery and coalescent-based methods–PTP, GMYC, spedeSTEM and *Beast–in order to test congruence among results. Additionally, the divergence times were also estimated to elucidate diversification events. Delimitations inferred from the different analyses are comparable with only minor differences, and following a conservative approach we propose that the sampled specimens of the P. reticulatum–P. pseudoreticulatum complex belong to at least eight distinct species-level lineages. Seven are currently classified under P. reticulatum and one as P. pseudoreticulatum. In this work we discuss one of only few examples of cryptic species that have so far been found in sorediate reproducing lichen forming fungi. Additionally our estimates suggest a recent origin of the species complex–during the Miocene. Consequently, the wide distribution of several of the cryptic species has to be explained by intercontinental long-distance dispersal events.


Introduction
Lichenized fungi form mutualistic associations either with green algae or cyanobacteria and reproduce sexually by ascospores or asexually by diaspores most commonly being isidia or soredia, which disperse both partners. Traditional criteria to delimit species are based on phenotypical characters, but in groups with limited morphological features, i.e. lichen-forming fungi, this approach has shown to underestimate species diversity in specific cases. Molecular studies have advanced our knowledge of cryptic diversity in lichens (reviewed in [1][2][3]). Identification leads to appearance of genetically distinct clusters, which are separated by long internal branches. A number of studies have employed this coalescent approach to successfully delineate species and detect cryptic species in various groups of organisms [33,[36][37][38].
SpedeSTEM uses maximum likelihood and information theory to evaluate phylogenetic models of lineage composition, enables biologists to identify distinct evolutionary lineages shortly after their formation. This approach requires operational taxonomic units (OTUs) to be designated a priori.
Additionally, putative species were evaluated based on the threshold of 0Á016 substitutions per site (s/s) that separate intra-and interspecific distances in parmelioid lichens [39].
Moreover, estimates of divergence dates of species is a major step to elucidate potential ecological, biogeographic and climatic patterns that drive diversification and explain the current distribution of species. In lichen-forming fungi, limited data are available on species turnover, diversification events and subsequently on factors driving diversification [8,9,[40][41][42][43][44]. This is mainly due to the poor fossil record for lichen-forming fungi and due to uncertainties in the interpretation of the few known fossil records [45][46][47]. It has been demonstrated that dating speciation events using priors on substitutions rates in a full probabilistic coalescent species tree framework-the approach implemented in Ã Beast-is appropriate for studies on taxonomic groups with poor fossil records [48,49].
We here aim to critically examine species diversity in the locally frequent P. reticulatum-P. pseudoreticulatum complex. Specifically, we aim to: i) delimit species applying coalescentbased species delimitation methods, and ii) elucidate divergence times. To this end, we gathered 180 samples from Africa, Asia, Australasia, Europe, North and South America and generated sequences of internal transcribed spacer of nuclear ribosomal DNA (ITS) and DNA replication licensing factor MCM7 (MCM7).

Taxon sampling
Sequence data of the ITS and the protein-coding MCM7 genes were analyzed in 28/21 specimens of P. pseudoreticulatum and 152/121 of P. reticulatum, collected from distant geographic regions throughout the species distributions. Specimens were collected in locations outside of national parks or protected areas therefore no specific permissions were required. Also this study does not involve endangered or protected species.
Parmotrema cetratum was selected as outgroup since it has been shown to be closely related [50]. For ITS, data from 94 individuals were generated for this study and 88 sequences downloaded from GenBank. For MCM7, data from 144 sequences was generated for this study. Detailed collection information and GenBank accession numbers are listed in S1 Table. final extension step at 72°C for 8 min was added, after which the samples were kept at 4°C. Amplification products were visualized on 1% agarose gels stained with SYBR 1 Safe DNA (Life Technologies Corporations, USA) gel stain (10 000× concentrated in DMSO) and subsequently purified using the enzyme exoSAP-IT (GE Healthcare, UK) according to the manufacturer's instructions. Both complementary strands were sequenced using Big Dye Terminator reaction kit (ABI PRISM, Applied Biosystems). Cycle sequencing reactions were performed with the same sets of primers used in the amplification step. Sequencing reactions were electrophoresed on a 3730 DNA analyzer (Applied Biosystems) at the Unidad de Genómica (Parque Científico de Madrid).

Sequence alignment
Sequence fragments generated for this study were assembled and edited using the program SeqMan v.7 (Lasergene R, DNASTAR, Madison, Wisconsin, USA). Sequence identity was confirmed using the mega-BLAST search function in GenBank. We used the program MAFFT v.6 [54] with the parameters set to default values to align the DNA sequences for each data set separately. Ambiguously aligned positions were identified and removed using G-Blocks implementing the options for a less stringent selection [55].

Phylogenetic analyses
The separate alignments and the combined data set were analyzed using maximum likelihood (ML) and a Bayesian Markov chain Monte Carlo (B/MCMC) approach.
We examined nodes to identify well-supported (ML bootstrap values > 70%) conflicts among the individual ITS and MCM7 phylogenies before combining the alignments [56,57].
The ML analysis was performed using an online version of the program RAxML v.7.0.4 (http://phylobench.vital-it.ch/raxml-bb/) [58,59], assuming a GTRGAMMA model which includes a parameter (Γ) for rate heterogeneity among sites and chose not to include a parameter for estimating the proportion of invariable sites. Nodal support was evaluated using 1000 bootstrap pseudoreplicates. For the combined data set we use a locus-specific model partitions in RAxML. Both loci were treated as separate partitions.
Separate and concatenated data sets were also analyzed with Bayesian inference as implemented in MrBayes v3.1.2 [60]. Models of DNA sequence evolution for each locus were selected with the program jModeltest v0.1 [61], using the Akaike information criterion (AIC) [62]. The model TIM2ef+G was selected for ITS; and SYM+I for the MCM7 region. The combined data set was partitioned into the two parts (ITS, MCM7), and each partition was allowed to have its own parameters [63]. No molecular clock was assumed. Two parallel runs of 3 million generations were done, starting with a random tree and employing 12 simultaneous chains each. Every 1000 th tree was saved into a file. The first 300,000 generations (i.e., 3000 trees) were deleted as the "burn-in" of the chains. The outputs of MrBayes were examined with the program Tracer v1.5 [64] to check for convergence of different parameters, determine the approximate number of generation at which log likelihood values stabilized and identify the effective sample size (ESS) for each parameter. Topological convergence in the two independent MCMC runs was checked with the "compare" plots in the program AWTY [65]. Posterior probabilities (PPs) of clades were obtained from the 50% majority-rule consensus of sampled trees after excluding the initial 10% as burn-in. Only clades that received bootstrap support equal to or above 70% in ML analyses and posterior probabilities equal to or above 0.95 were considered as strongly supported.

Calculation of genetic distances
Pairwise ML distances (given as the number of nucleotide substitutions per site) among the ITS rDNA sequences of the Parmotrema reticulatum-P. pseudoreticulatum complex were calculated (S1 Table). Genetic distances were calculated with TREE-PUZZLE 5.2 [67] using the GTR model of nucleotide substitution, assuming a discrete gamma distribution with six rate categories. The program generates an output file which consists of a triangular matrix with all pairwise distances between all the samples included. This matrix was visualized with the Microsoft Office program Excel 2000 and genetic distances between different specimens of the P. reticulatum-P. pseudoreticulatum complex were manually identified. Candidate species were proposed based on the threshold of 0Á016 substitutions per site (s/s) that separate intra-and interspecific distances in parmelioid lichens [39]. The distance values in the matrix 0Á016 s/s have been considered the values between the samples of the single species. The filter provided by Microsoft Excel was applied to separate values 0Á016, obtaining for every specimen included in the analysis the group of specimens with which it shares the values that characterize the species range.

Automatic Barcode Gap Discovery (ABGD) analyses
This is an automatic procedure that sorts the sequences into hypothetical species based on the barcode gap. This method automatically finds the distance where the barcode gap is located [32].
The ABGD method was carried out for the ITS dataset using the Web interface at http:// wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html. Default parameters were chosen using Kimura 2-parameter (K2P) distances that correct for transition rate bias (relative to transversions) in the substitution process. The default for the minimum relative gap width was set to different values between 0.1 and 0.15.

Poisson tree processes (PTP)
PTP does not require an ultrametric tree, as the transition point between intra-and inter-specific branching rates is identified using directly the number of nucleotide substitution [31]. PTP incorporates the number of substitutions in the model of speciation and assumes that the probability that a substitution gives rise to a speciation event follows a Poisson distribution. The branch lengths of the input tree are supposed to be generated by two independent classes of Poisson events, one corresponding to speciation and the other to coalescence. The ML phylogeny obtained with RAxML of the ITS data set, was used as the input trees to run PTP species delimitation analysis in the PTP webserver (http://species.h-its.org/ptp/). We ran the PTP analysis for 100,000 MCMC generations, with a thinning value of 100 a burn-in of 25%. Outgroup taxa were removed for species delimitation.

General Mixed Yule Coalescent (GMYC) species delimitation
This method requires a fully resolved tree with branch lengths estimates. It is based on the differences in branching rates between interspecific branching events and intraspecific relationships in a chronogram.
For the analyses, the ML tree obtained from a RAxML search using the ITS was used to infer an ultrametric tree using the program BEAST v.1.8.0 [68]. For this analysis identical haplotypes were removed. We used a site-specific GTR substitution matrix and a gamma distributed model of among-site rate heterogeneity with four discrete rate categories. We implemented an uncorrelated relaxed lognormal clock [69], and selected a Yule tree prior. Default values were used for remaining priors. MCMC analysis was run for a total of 10 million generations, sampling every 1000steps and excluding the first 10 million generations of each run as burn-in. Convergence was assessed by examining the likelihood plots through time using TRACER v1.4 [64].
GMYC requires a fully dichotomous chronogram and thus we used multdivtime to convert our chronogram into a fully dichotomous chronogram with internal branches of length zero, where appropriate. The ITS chronogram was then analyzed using the GMYC package in SPLITS in R (version 2.10, www.cran.r-project.org), using single and multiple threshold approaches [33,34].
The two outgroup samples (P. cetratum) were excluded from the data set using the drop.tip command in ape [70].
After optimization, we plotted the lineage through time (LTT) plot [71] with the threshold indicated and a chronogram that had the putative species indicated. Finally, we used the summary command to summarize the output statistics, including the results of the Likelihood ratio test (LRT) and the indication of the numbers of clusters and entities.

Validation approach: SPEDEStem
We used spedeSTEM v1.0 [35] as a validation approach to evaluate support for the candidate species previously delimited by PTP, ABGD, genetic distances and GMYC. This method calculates the probability of different models of lineage composition using maximum likelihood and evaluates these models using information theory [72].
SpedeSTEM takes as input ultrametric gene trees from independent loci that were inferred from the ML trees obtained from a RAxML search and using the program BEAST, as described above. SpedeSTEM requires an estimate of θ in order to scale the branch lengths in the species trees it produces. Using DnaSPv5 [73] we computed the average θ across the two included loci (theta = 0.03), with a scaling (segregating sites/total sites) of 1:1.6 for ITS and MCM7 respectively.
Analyses were run online in the SpedeSTEM server at: https://spedestem.osu.edu/ runspedestem. A table of models ranked by model probability is returned.
A total of seventeen species delimitation scenarios were tested and validated through spe-deSTEM. These candidate species were related with those proposed in our previous study [21]. Twelve of these putative species were selected from those identified in the GMYC multiple approach from the combined data set, as the less conservative scenario (data not shown). The remaining five were selected from the supported groups nested in the clade A2 in the phylogenetic concatenated tree.

Species tree and Divergence time estimates
We used the coalescent-based hierarchical Bayesian model Ã BEAST implemented in BEAST 1.8.0, as described elsewhere [49], to estimate a species tree following our proposed species delimitation scenario and to infer divergence dates using a coalescent-based species tree approach. Ã BEAST estimates the species tree directly from the sequence data, and incorporates the coalescent process, uncertainty associated with gene trees, and nucleotide substitution model parameters [74]. Further, species tree methods incorporating the process of gene lineage coalescence likely provide a more biologically realistic framework for dating divergence events [48]. In the absence of relevant fossil evidence for the P. reticulatum-P. pseudoreticulatum complex, we used the molecular evolution rates for the ITS marker [2.43 x 10− 9 substitution/ site/year (s/s/y)] recently reported for the related lichen-forming genus Melanelixia [9] to estimate the time to the most recent common ancestor (MRCA) for the clades. Implementing an uncorrelated relaxed lognormal clock, we selected a Yule process and gamma-distributed population sizes for the species-tree prior and a piecewise linear population size model with a constant root. Default values were used for remaining priors.
Two independent Markov chain Monte Carlo (MCMC) analyses were run for a total of 50 million generations, sampling every 1000 steps and excluding the 12500 trees as burn-in. We assessed convergence by examining the likelihood plots through time using Tracer and compared summarized tree topologies from separate runs; the effective sample sizes (ESS) of parameters of interest were all above 200. The posterior probabilities of nodes were computed from the sampled trees (excluding burn-in samples) using TreeAnnotator 1.8.0 [68].
While most of the putative species were represented by multiple individuals, the other three had two to three individuals. In such cases, species delimitation approach could be biased and caution must be taken. "Bayesian Phylogenetics and Phylogeography" (BP&P) [75], implementing the nearest neighbor interchange algorithm, could be a reliable method to handle disproportionate sampling, however, this approach requires a multi-locus dataset. Since we had only two loci in our dataset, we refrained from using this approach.

Phylogenetic analyses
The ITS data matrix included 180 OTUs of the P. reticulatum-P. pseudoreticulatum complex and 462 unambiguously aligned nucleotide positions (Tree-BASE No. S19485). 94 sequences were newly generated (S1 Table). For the Bayesian analysis the LnL value was -1851.772 with a standard deviation of ±318.04, and for ML the LnL value was -1553.1812.
The MCM7 data matrix included 142 OTUs of the P. reticulatum-P. pseudoreticulatum complex and 610 unambiguously aligned nucleotide positions (Tree-BASE No. S19485). 144 sequences were newly generated (S1 Table). For the Bayesian analysis the LnL value was -1885.326 with a standard deviation of ±316.98, and for ML the LnL value was -1581.497.
The partitioned ML analysis of the concatenated data matrix yielded the optimal tree with Ln likelihood value = -3321.011. The mean LnL value of the two parallel run of the Bayesian analysis for the two combined loci was -3682.9 with a standard deviation of ±0.55.
Since the topologies of the trees estimated from ML and Bayesian methods did not present any well-supported conflict, only ML topologies are shown with bootstrap and posterior probability values indicated on these topologies (Fig 1, S1 Fig and S2 Fig).
The topologies of the trees show that the samples of P. reticulatum and P. pseudoreticulatum split into multiple statistically supported clades. Statistical support for the different clades recovered in the concatenated analyses and in the independent ITS and MCM7 topologies are summarized in S2 Table. In the concatenated topology, specimens of P. reticulatum were spread between two main monophyletic groups (clades A1 and A2) with strong support in both analyses. These two clades split in different monophyletic groups (Fig 1).
It is remarkable to note that clade A1.1 splits in three different supported clusters in MCM7 topology, without showing any conflict with the ITS and concatenated topologies. Samples from the clades A1.2 and A1.3, in the concatenated topology, nested in a monophyletic group supported in the Bayesian analysis. In the MCM7 topology, samples from both clades are distributed on intermixed mode in two independent supported clades.
In clade A2.2 from the MCM7 topology are nested some samples that, in ITS, are grouped with samples of clade A2.1a. Additionally the two samples from clade A2.1a that are These conflicts did not influence the compatibility of the concatenated data set.
In the three topologies, two samples from Australia and two from Chile formed a different well-supported clade (clade D). Sixteen samples from the Canary Islands clustered in two independent supported groups (clades C and E).
Clade A2 further split in two supported clades, A2.1b and A2.3, which are present in both single locus genealogies. Clades C, D and E are also present in both single-locus genealogies.
Specimens of P. pseudoreticulatum split into two well-supported monophyletic groups (clades B and F) in ITS and in the concatenated tree. Additionally in the combined topology, clades B and F are grouped in a supported clade. In MCM7 topology, all the samples of P. pseudoreticulatum are nested in a single supported monophyletic group, in which samples from clades B and F are intermixed.
Phylogenetic relationships among the clades were not resolved with statistical support.

Identifying candidate species
The RAxML tree obtained from the concatenated data set was used to illustrate the delimitation of putative species recognized by the different approaches conducted with the ITS and with the concatenated dataset (Fig 1). ABGD and PTP analyses applied to ITS dataset detected 10 candidate species, which correspond to the well supported clades A1.1, A1.2 A1.3, A2, B, C, E and F obtained in phylogenetic analyses. However, clade D splits in two different candidate species: the first includes samples from Chile, and the second from Australia. The GMYC approach employing the multiple threshold (GMYC multiple ) on ITS dataset suggests almost the same candidate species as ABGD and PTP analyses with the only difference that all samples nested in clade D are included in a single species.
Genetic distance analyses in which a barcode gap is applied (0.016 s/s), also suggested the same clusters as those obtained in ABGD and PTP analyses, with one exception. Clade A2 is split into two different putative species, one of them includes the samples from the well-supported clade A2.1a, and the other grouped the well-supported clades A2.1b, A2.2 and A2.3 (Fig  1). Sample number 133 from clade A2.1a is included in both candidate species.
In both versions of the GMYC method (single and multiple threshold), the likelihood of the GMYC model was significantly higher than the likelihood of the null model. However the likelihood values of the single and multiple threshold analyses did not differ significantly. Therefore we show the results from multiple threshold analyses, which represent a more comparative approach.
All the samples nested in clade E were collected in different localities from the Canary Islands. Clades marked with an asterisk are present in both single locus analyses.
probabilities ! 0.95) in any of two analyses RaxML, and B/MCMC are in boldface. The branch that received strong support only in the ML analysis is indicated by a blue boldface line, whereas branches that received strong support only in the B/MCMC analysis are indicated by a black boldface line. Branches that were strongly supported in both analyses are indicated by a red boldface line. Species delimitation scenarios obtained from different methods are indicated in columns to the right (details discussed in text). The 8-species scenario recovered in species tree estimation using *BEAST is also indicated in the right column. Finally, a species tree was estimated with Ã BEAST following our proposed species delimitation scenario by SPEDEStem. An 8-species scenario was recovered in species tree estimation, which corresponds to the species validated by SpedeSTEM with the exception of samples from clade A2.3 that are not validated as a separate species but they are included in the clade A2 species.

Molecular dating analysis
The estimated timing of diversification events in the P. reticulatum-P. pseudoreticulatum complex is shown in the rate calibrated species tree inferred using Ã BEAST (Fig 2). The mean node ages and divergence date ranges (95% highest posterior density intervals, HPD) of the clades are shown in the figure.
The crown of the complex was estimated at 9.6 Ma (95% HPD = 5.

Discussion
Previous molecular phylogenetic studies have provided evidence for the existence of several cryptic lineages hidden under the names of P. reticulatum and P. pseudoreticulatum [15,21]. The specific aim of this paper was to evaluate eight cryptic lineages in P. reticulatum and two in P. pseudoreticulatum that were identified and proposed as candidate species in our previous work [21]. From an extended taxon sampling and a two locus data set, we have applied recently developed approaches used to delimitate and validate species boundaries, and in addition, we studied the timing of diversification events in the group.
In the current study six major well-supported clades were found in the combined analyses of P. reticulatum-P. pseudoreticulatum complex. These clades correspond to those previously discussed, with the exception of clade C. Additionally, in our previous studies samples of P. pseudoreticulatum fell in two independent well supported clades (B and F), whereas in the current study they formed one statistically supported monophyletic clade. Relationships among clades remained largely unsupported as in our previous studies that employed ITS data only [15,21].
Phylogenetic analyses alone are insufficient to draw conclusions on species boundaries in this complex and additional studies have been applied. To achieve confidence in species delimitation we used an integrative taxonomic approach that uses complementary sources of data as ecology, chemistry, morphology, biogeography, and multiple molecular analytical approaches. However, in the P. reticulatum-P. pseudoreticulatum complex these data are limited as in other cryptic lineages of lichenized-fungi. Only in two clades (C and E) a geographical pattern has been found with all samples coming from the Canary Islands. For the other lineages we could not observe any correlation with geographical, morphological, chemical or ecological patterns. An increasing number of studies have revealed the presence of cryptic species in lichen-forming fungi without any recognizable phenotypical or biogeographical characterization of these clades [2,3]. In these cases, coalescent-based approaches and other recently developed DNA-based methods have been widely used as an objective measure to delimitate and validate species boundaries [38,[76][77][78].
Delimitations inferred from the different analyses are comparable with only minor differences. The analyses supported between 9 and 11 candidate species in the P. reticulatum-P. pseudoreticulatum complex. Differences among methodological approaches are common [76] due to different statistical power to detect independent lineages and could indicate that assumptions of one or more of the methods are violated. In this work, differences in the circumscriptions of clades A2 and D were detected.
In clade A2 there is a conflict between PTP / ABGD and our manually calculated genetic distances. Also, the lineage labelled as "b" and group A2.3 showed gene genealogical concordance, which is taken as evidence that these clades represent reproductively isolated lineages [2,3]. Therefore it seems that the speciation process has not been completed within this clade. When number of nucleotide substitutions are considered to assign DNA sequences to species relied on the existence of a barcode gap [76], one might expect to find overlap between inter and intra-specific divergences for clades that are speciating [76]. Genetic distance calculations proposed two species instead of the one proposed by PTP and ABGD, due to a higher statistical power to detect species lineages as we consider a specific barcode gap calculated for parmelioid lichens [31,79].
Another difference was found in clade D. In this case, GMYC regarded it as one species, whereas other methods identified two different candidate species. In our study GMYC suggested more conservative species delimitations, which differs from other studies in which oversplitting in comparison to ABGD or PTP was found [31,79].
Based on the proposed candidate species by the different methods tested in this work and in the proposed candidate species in our previous work, we have implemented a validation coalescentbased approach through spedeSTEM analysis, which pointed 9 species scenario as the most likely.
Finally, an 8-species scenario was recovered in species tree estimation by Ã BEAST. This largely corresponds to the species validated by SpedeSTEM with the exception of samples from clade A2.3 that are not validated as a separate species however they are included in the clade A2 species.
Ã BEAST is as a coalescent-based species tree method that can accurately depict relationships even in cases where incomplete lineage sorting and gene tree heterogeneity obscures phylogenetic relationships among species [29,30]. However spedeSTEM may be less accurate in cases of recent speciation events [76].
With the data at hand we are following a conservative approach advocated by Carsten and coworkers [76] who argued that it would be better to fail to delimit species to falsely delimit entities that do not represent independent lineages.
We propose that the sampled specimens of the P. reticulatum-P. pseudoreticulatum complex belong to at least eight distinct species-level lineages. Seven currently classified as P. reticulatum and one as P. pseudoreticulatum. Although most analyses proposed the presence of two species lineages within P. pseudoreticulatum, we tentatively follow a conservative approach as suggested by spedeSTEM and Ã BEAST. The formal description of the new species will be done in a companion paper.
The Parmotrema reticulatum complex has traditionally been considered a purely asexually reproducing fungal group. Further, asexually reproducing species in lichens-sorediate-and in filamentous fungi in general has largely been observed as an evolutionary dead end or clones doi:10.1371/journal.pone.0161031.g002 [80]. By contrast, we here demonstrate high cryptic species diversity within the P. reticulatum complex (Fig 1). Moreover, recently, there has been a growing body of evidence in many other fungal groups with a supposedly asexual lifestyle that they indeed have the ability to endure sexual reproduction; suggesting cryptic sexuality [81,82]. Numerous criteria, such as observation of initial morphological stages of sexual reproductive organs, finding of recombination and identification of mating-type genes involved in sexual reproduction in the genome have been used to provide evidence of cryptic sexuality in putative asexual species [83,84]. In the P. reticulatum complex, ascomata or ascomatal primordia are rarely found [22,23]. However, given the evidence of cryptic sexuality in ascomycetes we hypothesize that cryptic sexuality may contribute to the genetic diversity and speciation in the P. reticulatum complex. A test of our hypothesis can be the search for mating-type genes responsible for sexual reproduction in the genome of that species complex.
Our divergence time estimates suggest the origin of the species complex in the Miocene (Tortonian) around 9.6 Ma. This is consistent with the estimated age of the genus Parmotrema [40] at around 13 Ma. Our data suggest that the major diversification events within the complex occurred during the Pliocene with more recent radiations at the tips during the Pleistocene. Although rate-calibrated divergence estimates must be interpreted with caution, these results are comparable with those found in Melanohalea [8]. Miocene and Pliocene were characterized by cooling and a series of climate fluctuations including ice ages and uplifting of major mountain systems in Eurasia [76]. These events accelerated turnover of the terrestrial biota and several studies have shown that the Miocene has been an important time for diversification of species groups in Parmeliaceae [76]. Specifically the origin of the genus and the species complex corresponds to an event of increasing aridity in the Miocene . This climatic event, together with habitat changes due to previously indicated mountain uplift and associated shifts in vegetation patterns could have played a key role in the origin of the genus and the diversification of the complex. Radiations within the complex occurred from the end of the Pliocene when the climate became cooler, drier and seasonal.
The cryptic species identified in our study are currently widely distributed in different areas from the Northern and Southern Hemisphere, including Europe, North Africa and Canary Islands, Australia, South America and Asia. This wide distribution has to be explained by intercontinental long-distance dispersal events given the recent origin of the species complex.  Table. Collection information for all the specimens included in the present study. Gen-Bank accession numbers for the two sampled loci: nuclear ribosomal internal transcribed spacer region (ITS), and DNA replication licensing factor MCM7 (mcm7). Newly generated sequences for this study are indicated in boldface. (DOC) S2 Table. Statistical support for the clades recovered in the concatenated, ITS and mcm7 topologies. Clades do not recovered in the topology are indicated with an asterisk ( Ã ). (DOC)