Transoceanic Dispersal and Subsequent Diversification on Separate Continents Shaped Diversity of the Xanthoparmelia pulla Group (Ascomycota)

In traditional morphology-based concepts many species of lichenized fungi have world-wide distributions. Molecular data have revolutionized the species delimitation in lichens and have demonstrated that we underestimated the diversity of these organisms. The aim of this study is to explore the phylogeography and the evolutionary patterns of the Xanthoparmelia pulla group, a widespread group of one of largest genera of macrolichens. We used a dated phylogeny based on nuITS and nuLSU rDNA sequences and performed an ancestral range reconstruction to understand the processes and explain their current distribution, dating the divergence of the major lineages in the group. An inferred age of radiation of parmelioid lichens and the age of a Parmelia fossil were used as the calibration points for the phylogeny. The results show that many species of the X. pulla group as currently delimited are polyphyletic and five major lineages correlate with their geographical distribution and the biosynthetic pathways of secondary metabolites. South Africa is the area where the X. pulla group radiated during the Miocene times, and currently is the region with the highest genetic, morphological and chemical diversity. From this center of radiation the different lineages migrated by long-distance dispersal to others areas, where secondary radiations developed. The ancestral range reconstruction also detected that a secondary lineage migrated from Australia to South America via long-distance dispersal and subsequent continental radiation.


Introduction
Methods for delimiting species, the fundamental taxonomic unit, have always fascinated evolutionary biologists [1][2][3]. Understanding the circumscription of species is important for biological and ecological studies and for conservation issues. However, the main challenge is to recognize species in organisms with relatively simple morphologies. In lichenized fungi traditional species circumscriptions are based on phenotypic characters, such as thallus and ascomatal morphology or chemical characters. However, there is a growing body of evidence from molecular studies that the traditional morphology-based species circumscriptions are insufficient to represent the diversity in lichenized ascomycetes . A number of DNA sequence-based phylogenetic studies revealed the presence of distinct lineages within currently delimited species. Subsequent, detailed studies often revealed previously overlooked morphological subtleties or chemical differences among those clades and authors often refer to these as ''semi-cryptic'' species [8].
On a par with the phenotypic-based species circumscription, researchers often accepted wide distribution ranges for species occurring on different continents. This was at least partially due to a common belief by lichenologists in the ''everything is everywhere'' hypothesis [27,28] applied to fungi, as discussed elsewhere [8,29]. In several cases molecular data assisted in a better understanding of the biogeography of lichen-forming fungi where taxa were shown to represent different species on different continents, e.g. the Leptogium furfuraceum aggr. on different continents [18], Melanelixia glabra s. lat. in Europe and North America [14], Parmelina quercina s. lat. on different continents [4], Physcia aipolia aggr. in Europe and Australia [9], Xanthoparmelia spp. in North America and Australia [15,16]. In the Leptogium furfuraceum aggr. complex sister-group relationships were found between populations from the same hemispheres which were incongruent with previous classifications based on morphological differences [18], and the dated phylogeny indicated that the species had migrated via transoceanic dispersal to different continents.
Here we report another case of a group of lichenized fungi where transoceanic dispersal to different continents is correlated with the phylogenetic lineages. The group studied here is the Xanthoparmelia pulla group which belongs to the family Parmeliaceae. This family represents one of the largest families of lichenized fungi [30,31]. The main clade of the family is the parmelioid clade with almost 2000 species [32] currently classified in 27 genera, with Xanthoparmelia being the largest with over 800 accepted species [33]. The species in this genus characteristically occur on siliceous rocks or soil, predominantly in arid to subarid regions, with a center of distribution in the southern hemisphere. The genus is characterized by having cell wall polysaccharides of the Xanthoparmelia-type, small ascospores with an arachiform vacuolar body [34], and the presence of a pored epicortex [33,35]. It has been hypothesized that the genus diversified in a rapid radiation following a shift towards drier habitats at the base of the Xanthoparmelia clade [36] leading to the high current diversity.
Previously, the Xanthoparmelia pulla group has been classified within the separate genus Neofuscelia based on the different cortical chemistry (having melanoid pigments and lacking usnic acid or atranorin, characteristic of the majority of Xanthoparmelia species) [37,38]. A subsequent molecular study showed that the genus Neofuscelia was polyphyletic, with its clades scattered within Xanthoparmelia [33]. Consequently, the genus Neofuscelia was reduced to synonymy with Xanthoparmelia, as have other genera previously distinguished by cortical chemistry or growth form [33,[39][40][41][42][43]. The Xanthoparmelia pulla group is a monophyletic clade within the complete Xanthoparmelia clade, that includes the former Esslinger's Xanthoparmelia pulla species and other related species [44].
Although the clades were largely incongruent with the current species circumscription, we found a correlation of the main clades of X. pulla group with their geographical distribution and chemical profile ( Fig. 1, 2). Clade 1 includes specimens from California, Macaronesia and areas around the Mediterranean basin, all of which contain depsides and depsidones derived from the orcinol pathway or with aliphatic acids; clade 2 includes specimens from Australia (subclade 2.1) and South America (subclade 2.2) with depsides and depsidones derived from the orcinol pathway; clade 3 derives from South African specimens containing olivetoric acid; clade 4 specimens with hypostictic acid from California (subclade 4.1) and South America (subclade 4.2); and clade 5 specimens with physodic acid (orcinol depsidones) from South Africa.
The species delimitation within the Xanthoparmelia pulla group is currently based on a combination of morphological and chemical characters ( Table 1). The morphological characters include the color of the lower surface, shape of the lobes, attachment to the substrate, and presence of vegetative propagules while the chemical differences pertain to upper cortical and medullary secondary metabolites. A number of the currently accepted species have a wide distribution spanning several continents. To address the species delimitation in this group and to test the hypothesis of widely distributed species we have generated a data set using two loci (nuITS rDNA, nuLSU rDNA) from specimens collected on different continents. The molecular data were used to perform phylogenetic reconstructions in a maximum likelihood (ML) and Bayesian (B/MCMC) framework. We have also estimated the timing of the diversification events leading to the main clades found in our study to discriminate between vicariance and longdistance dispersal as possible explanations for the current distribution patterns. A Bayesian-based approach of ancestral range reconstruction was used to identify potential areas in which the group and major clades within the group originated.

Phylogenetic analyses
One hundred sixty-eight DNA sequences of ITS and nuLSU rDNA of 88 representative specimens of Xanthoparmelia were assembled. One hundred forty of these sequences were newly generated in this study. The specimens included 25 currently accepted species in the Xanthoparmelia pulla group, four unassigned specimens, and six samples of four species as outgroup. A data matrix of 1283 unambiguously aligned characters, with 454 characters in the ITS and 829 characters in the nuLSU rDNA was used for phylogenetic analyses. The data set included 1081 constant characters. The general time-reversible model with a gamma distribution and invariant model of rate heterogeneity (GTR+I+G) was employed for analyses of the single-loci and concatenated data sets. Since no strongly supported conflicts between the two single-locus ML phylogenetic trees were detected, a combined data set was analyzed. In the B/MCMC analysis of the combined data set, the likelihood parameters in the sample had the following averaged values for the partitioned data set The phylogenetic estimates of the ML and B/MCMC analyses were congruent, hence only the ML tree ( Fig. 1) is shown here. Specimens of the Xanthoparmelia pulla group form a strongly supported monophyletic group with five main, mostly wellsupported, clades (Fig. 1). The clades do not agree with the current species circumscription, with 11 species being polyphyletic, five of them with specimens from different continents entering different clades. For example, all the Northern Hemisphere specimens identified as X. luteonotata, X. pulla, X. delisei or X. glabrans belong to clade 1 while all the Australian specimens of the same species belong to clade 2.1. Similarly, specimens of X. imitatrix from South America belong to clade 2.2 and are not directly related to the South African specimen.

Estimates of divergence times and ancestral range reconstructions
A Bayesian phylogenetic tree was dated to estimate the age of the X. pulla group and its main clades. The results of the divergence time analysis are summarized in Fig. 2, and the whole parmelioid tree is shown as supp. mat. (Fig. S1). The Xanthoparmelia pulla group started to diversify around 11.61 Ma (7.61 -16.50 Ma), the age of the crown node of clade 1 was estimated at 5.31 Ma (3.01 -8.38 Ma), the ancestor of clade 2 around 8.10 Ma (5.13 -11.74 Ma), and the crown of clade 4 around 3.44 Ma (1.56 -6.07 Ma).
The results of the ancestral range reconstruction analyses are summarized in Figs. 1 and 2. This established that South Africa was the most likely origin of the X. pulla group, with a marginal probability of 0.745, indicating localized uncertainty. The four other areas explored (South America, Australia, California and the Mediterranean basin) were rejected with probabilities below 0.05. For the base of clade 1, the Mediterranean basin was reconstructed as the most likely ancestral range with a marginal probability of 0.555, but California could not be rejected (probability of 0.104). For clade 2, Australia was recovered as the most likely ancestral area with marginal probability of 0.672; while South America, the other area from which specimens of this clade occur, was rejected as potential ancestral area (p,0.05); similarly, California, the Mediterranean basin and South Africa were also rejected as ancestral areas. South America was found to be the most likely origin for clade 4 (which also includes specimens from California and South America) with a marginal probability of 0.48, although neither California nor South Africa were rejected (probabilities of 0.106 and 0.083, respectively). Australia and the Mediterranean basin were rejected as ancestral ranges for clade 4. Table 1. Main differences of the species of Xanthoparmelia pulla group studied in this paper [38,44,47,59,60].

Discussion
Understanding the diversity and delimiting species in lichenized fungi has been a long standing challenge and current studies using molecular data have dramatically changed our ability to distinguish species in this group [8,23,45]. The Xanthoparmelia pulla group is a good example for illustrating the difficulties in in distinguishing species by morphology due to the remarkable  [44,46,47]. Following the current classification using a combination of vegetative morphology and secondary chemistry, a number of species have broad geographical distributions spanning several continents.
Here we have used molecular data to investigate the current classification within the group and attempt to explain their distribution. We used likelihood-based and Bayesian approaches to investigate the evolutionary origin of the group and timing of speciation events. Hopefully such data will reveal evolutionary patterns so we may develop a framework for their taxonomic classification which better reflects the phylogenetic relationships in the X. pulla group. Our results clearly indicate that the species as currently delimited are polyphyletic (Fig. 1). This is consistent with results from other studies of Xanthoparmelia species believed to occur on different continents which were subsequently found to represent distinct lineages [15,16]. Further, similar patterns have been found in other groups of lichenized fungi [4,14,18,48].
The ancestral range reconstruction points to South Africa as the most likely origin of the X. pulla group. Although there is a certain degree of uncertainty in the reconstruction (marginal probability of 0.745), the analysis rejected other areas as potential ancestral areas for the group. Interestingly, South Africa has the highest morphological and chemical diversity within the group and the specimens studied here belong to different, unrelated lineages ( Fig. 1). South African specimens containing olivetoric acid cluster in clade 3 and those with physodic acid in clade 5. The phylogenetic relationships of other specimens from South Africa with hypostictic acid, physodic acid or other orcinol depsides and depsidones are still unresolved. The South African specimens show remarkable morphological variability, including subcrustose and foliose species. Further, many Xanthoparmelia species occur in arid climates and the diversification of X. pulla group occurred around 11.61 Ma (7.61 -16.50). At this time the Cape Region underwent a major aridification [49], which may be responsible for the rapid radiation and current richness of the Cape flora. Thus, it is likely that the X. pulla group originated in South Africa around the same time. Unlike most Cape region elements in flowering plants, the species of the X. pulla group subsequently extended their distribution by transoceanic dispersal.
Within the X. pulla group, the five lineages identified are characterized by the presence of different substance classes and in some cases they diverged secondarily in different geographical areas. The correlation between chemical pathways and the lineages found in molecular studies has also been found in Pertusariaceae among lichenized fungi [50][51][52].
Clade 1 includes specimens containing orcinol depsides and depsidones that occur in California and the Mediterranean basin. Neither internally supported subclades nor a geographical pattern was found within this clade, and specimens with different phenotypical characters from different geographical areas are  intermingled. In fact, shifting between orcinol depsides and depsidones can occur by one-step transformations [46]. By contrast, in other genera (e.g. Melanelixia, Parmelina, Leptogium) with similar disjunct distributions (North America and the Mediterranean basin), the geographical distribution correlates with the clades found in the molecular study. In the X. pulla group this pattern was not found, possibly due to insufficient sampling or absence of a phylogenetic signal in the markers used. This might be due to the slower evolutionary rates of lichenized fungi from arid and subarid regions compared to oceanic parmelioid lichens [36]. The diversification age of clade 1 was estimated at 5.31 Ma (3.01 -8.38 Ma), at the end of the Miocene, which was a geological period when numerous groups radiated in arid conditions. The Mediterranean region was suggested as the ancestral area of this clade by the ancestral range analysis, but the result was poorly supported. Clade 2 also contains species with orcinol depsides and depsidones and comprises two disjunct lineages, one occurring in Australia (clade 2.1) and the other in South America (2.2). The estimated age for clade 2 (8.10 Ma; 5.13 -11.74 Ma) rules out the possibility of vicariance, since the breakup of Australia, Antarctica and South America occurred between 35-52 Ma ago [53]. The ancestral range reconstruction points to Australia as the ancestral area of clade 2. This would be consistent with long distance dispersal from Australia to South America, a phenomenon frequently found in many plants groups [54]. Within the Australian clade several strongly supported lineages are not consistent with the current species delimitation of the group, indicating that the phenotypical characters used to distinguish species in the group have limited phylogenetic validity. Similar disparities between phylogenetic relationships and current species delimitations were found within the yellow Xanthoparmelia species from western North American [55,56].
Clade 4 includes specimens containing hypostictic acid from California and South America. Here again all the South American specimens form a monophyletic group. Specimens of X. subhosseana occurring in different continents are not closely related. The most likely ancestral origin of clade 4 is South America (marginal probability of 0.48), although neither California nor South Africa could be rejected.
Our study indicates that the X. pulla group started to radiate during the Miocene in South Africa, where the highest diversity of this group is found. From this region, different lineages with distinct secondary metabolites belonging to different chemical pathways were dispersed to other regions, where they experienced rapid and more recent radiations. In some cases our results showed that the sympatric species of the X. pulla group in an area belong to distantly related groups. For example, the Californian X. pulla flora includes species from clade 1 and clade 4.1, the latter most probably having migrated from South America. Indeed our study indicated that the current taxonomic circumscription of species in the group does not agree with the evolutionary hypotheses inferred by molecular markers. The incongruence of phenotype-based classification and molecular phylogeny is a challenge for the classification of these fungi. Additional studies will be needed to determine whether the lineages found here represent cryptic species or whether new phenotypical characters can be found to distinguish these distinct lineages (as has been found in some other Parmeliaceae [10,57]). Future research should address how such parallel evolution of phenotypical characters in lichenized fungi could be explained in order to provide a better framework to test the adaptive value of these characters [58]. Our results here have important implications for conservation and ecological issues, since species were found to have much more restricted distribution than previously thought.

Taxon sampling
Eighty two specimens of the Xanthoparmelia pulla group from California, the Mediterranean basin, Macaronesia, South America, Australia and South Africa were used for the phylogenetic study. The specimens were identified following the current species delimitations [38,44,47,59,60]. Chemical constituents were identified using thin layer chromatography (TLC) [61][62][63][64], and gradient-elution high performance liquid chromatography (HPLC) [65]. The major medullary compounds were classified into four major groups based on their chemical structure: 1) Orcinol depsides: olivetoric, divaricatic, stenosporic and glomelliferic acids. 2) Orcinol depsidones, physodic, alectoronic and glomelliferonic acids. 3) b-Orcinol depsidones: stictic acid (only present in the outgroup) and hypostictic acid. 4) aliphatic acids: constipatic acid. All necessary permits were obtained for the described field studies. Collecting permits in Australia were all obtained by J.A. Elix (ca. 50 permit numbers for each states and over several years) and in Chile by W. Quilhot. For European locations specific permission was not required, since the locations were neither in privatelyowned or protected areas. The field studies did not involve endangered or protected species.

Molecular study
Total DNA was extracted from frozen lobes of thalli crushed with sterile glass pestles, using the DNeasy Plant Mini Kit (Qiagen) following the manufacturer's instructions and modifications of Crespo et al. [66]. The following primers were used: ITS1-LM [67] and ITS2-KL [68] for nuITS rDNA, and LR0R and LR5 [69] for nuLSU rDNA.
For each amplification we used a reaction mixture of 25 mL, containing: 2.5 mL of 10x DNA polymerase buffer (including MgCl 2 2mM) (Biotools), 1.25 mL of each primer, 0.75 mL of DNA polymerase (1U/mL), 0.5 mL of dNTPs containing 10 mM of each base (Biotools), 5 mL of DNA (third elution of DNA extraction) and 13.5 mL dH 2 O. Amplifications were carried out in an automatic thermocycler (Techne Progene 3000) with the following steps: an initial denaturation at 94uC for 5 min; 35 cycles of 94uC for 1 min, 58uC (nuITS rDNA) or 56uC (nuLSU rDNA) for 1 min, and 72uC for 1.5 min; a final extension at 72uC for 5 min. PCR products were cleaned with DNA Purification Kit (Flavorgen) and sequenced with the same primers using the ABI Prism Dye Terminator Cycle Sequencing Ready reaction kit (Applied Biosystems) with the following program: initial denaturation at 94uC for 3 min, 25 cycles at 96uC for 10s, 50uC for 5s and 60uC for 4 min. Sequencing reactions were electrophoresed on a 3730 DNA analyzer (Applied Biosystems). The sequence fragments were assembled with Bioedit v. 7.0 [70] and manually adjusted.

Sequence alignment and selection of the substitution model
We used a dataset of 2 loci of 82 specimens representing 25 species of the Xanthoparmelia pulla group and 6 specimens as outgroup. The sequences were mainly generated in this study (140 sequences) and the others taken from our previous studies [33,41]. The outgroup selection was based on previous phylogenetic studies [41]. GenBank accession numbers and sources of the specimens are listed in Table 2.
The two loci were aligned separately with Muscle 3.6 [71] and the ambiguous positions were removed manually. The general time reversible model including estimation of invariant sites (GTR+I+G) was selected by jModelTest v 0.1.1 [72] as the most appropriate nucleotide substitution model for both loci.

Phylogenetic Analyses
Potential conflict between the two loci was assessed by comparison of the ML analyses obtained with Garli 0.96 [73] for each locus, using 100 pseudoreplicates for the bootstrap analyses. The phylogenetic analyses of the combined matrix were done using maximum likelihood (ML) and a Bayesian approach. ML analysis was performed using Garli 0.96 [73] with default settings and 100 replicates for the bootstrap analyses. The Bayesian analysis was performed using MrBayes 3.1.1 [74] using the GTR+I+G model, and the data set partitioned into nu ITS and nu LSU. Each partition was allowed to have its own parameter values [75]. Heating of chains was set to 0.2, with 5 million generations sampled every 500th tree. The first 1000 trees were discarded as burn in. We used AWTY [76] to compare splits frequencies in the different runs and to plot cumulative split frequencies to insure that stationarity was reached. Of the remaining 18000 trees (9000 from each of the parallel runs) a majority rule consensus tree with average branch lengths and posterior probabilities was calculated using the sumt option of MrBayes. Clades with bootstrap support equal or above 70 % under ML and/or posterior probabilities $0.95 in the Bayesian analysis were considered as strongly supported. Phylogenetic trees were visualized using the program FigTree [77].

Calibration of nodes and dating analysis
The ages of the X. pulla group and its major clades were estimated by a divergence time analysis based on a calibrated phylogeny of the parmelioid lichens [78]. We used a matrix of two loci (nu ITS and LSU) with a proportional number of samples of each parmelioid clade to have a representative tree and trend in speciation through time [79]. The matrix included 299 specimens of parmelioid lichens and 3 specimens of the genus Usnea (as outgroup); 62 new sequences of Xanthoparmelia species outside the X. pulla group were included. GenBank accession numbers with the specimens of the dating analysis are listed in Table S1. The sequence alignment, selection of the nucleotide substitution model, and phylogenetic analyses were done using the same procedures used for the Xanthoparmelia pulla dataset (see above).
The divergence time analyses were performed using BEAST v.1.6.1 [80]. We used a starting tree obtained from a ML analysis using Garli 0. 96 [73] of the concatenated dataset, calculated an ultrametric tree using nonparametric rate smoothing (NPRS) implemented in TreeEdit v.10a10 [81]. The age of the crown node of the parmelioid lichens was calibrated at 60 Ma, following Amo de Paz et al. [78]. The starting tree was topologically congruent with the parmelioid phylogeny presented in Crespo et al. [32].
For the divergence time analyses we used two points of calibration: the age of the crown node of the parmelioid lichens set at 60.28 Ma (49.81 -73.55 Ma) [78], and the age of the crown node of the genus Parmelia (dated from the fossil Parmelia ambra from the Dominican amber, 15-45 Ma, [82] as discussed previously) [78].
The BEAST analysis was performed using the GTR+I+G substitution model, a Birth-Death process tree prior, and a relaxed clock model (uncorrelated lognormal) for the concatenated dataset. Calibration points were defined as prior distribution, minimal ages and calibrated with a lognormal distributions: 1) the parmelioid crown node at uniform distribution between 49 -73 Ma; 2) the Parmelia crown node at log-normal mean = 2.77, offset = 14, lognormal standard deviation = 0.5. The analysis was run for 10 million generations, with parameter values sampled every 1000 generation. We checked the stationary plateau with Tracer v. 1.5 [83]. We discarded 10% of the initial trees as burn in and the consensus tree was calculated using Tree Annotator v 1.6.1 [80]. The results were visualized with FigTree v. 1.3.1 [84]. Ages of the X. pulla clades were estimated for nodes with more than 0.95 of posterior probability in the BEAST runs and in the previous Bayesian analysis.

Ancestral range reconstructions
The biogeographical analysis to reconstruct the ancestral area was performed using an indirect Bayesian approach to character state reconstruction [85] implemented in SIMMAP v1.5 [86] following [87]. This analysis integrates the combination of the uncertainty in the tree, branch lengths and the substitution models using Markov chain Monte Carlo. We treated the biogeographic regions as discrete characters. The major areas in which X. pulla species are distributed were categorised broadly into five areas: California, Mediterranean basin (including Macaronesia), South America, South Africa, and Australia. Presence/absence was coded as binary states and each area was given equal probability. We performed the ancestral state reconstruction analysis on a subsample of 1000 trees derived from the MrBayes tree sampling of the Xanthoparmelia pulla group.
We also performed ancestral range reconstruction analysis using dispersal-extinction-cladogenesis (DEC) implemented in Lagrange program [88]. The results were inconclusive due to the lack of confidence in parts of the X. pulla phylogeny and hence the results are not included in this paper. Figure S1 Chronogram of parmelioid lichens focusing in Xanthoparmelia pulla group. Calibration points: A, inferred age of radiation of parmelioid lichens and B, age of the Parmelia fossil. The Xanthoparmelia pulla group is highlighted by a box and the dated clades are indicated by a branch in bold.

(TIF)
Table S1 GenBank accession numbers of parmelioid lichens (except X. pulla group, see Table 2) used for divergence time analysis. New sequences are in bold. (DOC)