This article appeared in a journal published by Elsevier. The attached copy is furnished to the author for internal non-commercial research and education use, including for instruction at the authors institution and sharing with colleagues. Other uses, including reproduction and distribution, or selling or licensing copies, or posting to personal, institutional or third party websites are prohibited. In most cases authors are permitted to post their version of the article (e.g. in Word or Tex form) to their personal website or institutional repository. Authors requiring further information regarding Elsevier’s archiving and manuscript policies are encouraged to visit: http://www.elsevier.com/copyright Author's personal copy Molecular Phylogenetics and Evolution 59 (2011) 587–602 Contents lists available at ScienceDirect Molecular Phylogenetics and Evolution journal homepage: www.elsevier.com/locate/ympev Complex patterns of speciation in cosmopolitan ‘‘rock posy’’ lichens – Discovering and delimiting cryptic fungal species in the lichen-forming Rhizoplaca melanophthalma species-complex (Lecanoraceae, Ascomycota) Steven D. Leavitt a,⇑, Johnathon D. Fankhauser b, Dean H. Leavitt c, Lyndon D. Porter d, Leigh A. Johnson a, Larry L. St. Clair a a Department of Biology and M.L. Bean Life Science Museum, Brigham Young University, Provo, UT 84602, USA Department of Plant Biology, University of Minnesota, 1445 Gortner Ave., St. Paul, MN 55108, USA c Department of Biology, San Diego State University, 5500 Campanile Drive, San Diego, CA 92182-4614, USA d USDA_ARS Vegetable and Forage Crops Research Unit, Prosser, WA 99350, USA b a r t i c l e i n f o Article history: Received 6 August 2010 Revised 16 December 2010 Accepted 17 March 2011 Available online 1 April 2011 Keywords: Coalescent theory Species delimitation Lichen species concepts Rhizoplaca Secondary metabolites Species delimitation Species tree phylogeny estimation Vagrant lichens a b s t r a c t A growing body of evidence indicates that in some cases morphology-based species circumscription of lichenized fungi misrepresents the number of existing species. The cosmopolitan ‘‘rock posy’’ lichen (Rhizoplaca melanophthalma) species-complex includes a number of morphologically distinct species that are both geographically and ecologically widespread, providing a model system to evaluate speciation in lichen-forming ascomycetes. In this study, we assembled multiple lines of evidence from nuclear DNA sequence data, morphology, and biochemistry for species delimitation in the R. melanophthalma species-complex. We identify a total of ten candidate species in this study, four of which were previously recognized as distinct taxa and six previously unrecognized lineages found within what has been thus far considered a single species. Candidate species are supported using inferences from multiple empirical operational criteria. Multiple instances of sympatry support the view that these lineages merit recognition as distinct taxa. Generally, we found little corroboration between morphological and chemical characters, and previously unidentified lineages were morphologically polymorphic. However, secondary metabolite data supported one cryptic saxicolous lineage, characterized by orsellinic-derived gyrophoric and lecanoric acids, which we consider to be taxonomically significant. Our study of the R. melanophthalma species-complex indicates that the genus Rhizoplaca, as presently circumscribed, is more diverse in western North American than originally perceived, and we present our analyses as a working example of species delimitation in morphologically cryptic and recently diverged lichenized fungi. Ó 2011 Elsevier Inc. All rights reserved. 1. Introduction Lichens are obligate symbiotic systems consisting of a filamentous fungus, a photosynthetic partner (eukaryotic alga and/or cyanobacterium), and, at least in some cases, non-photosynthetic bacteria (Cardinale et al., 2008; Grube et al., 2009; Hodkinson and Lutzoni, 2009; Selbmann et al., 2010). The lichenized condition has been extremely successful for many fungal lineages, with an estimated 40% of all ascomycetes forming lichens (Lutzoni et al., 2001). Traditionally, morphology and the expression of signature secondary metabolites have been used to infer taxonomic boundaries for lichen-forming fungi (Culberson, 1972; Hale, 1990; Huneck, 1999; Huneck and Yoshimura, 1996). However, these characters are often widely variable, and their homology has ⇑ Corresponding author. Present address: Department of Botany, Field Museum of Natural History, 1400 S. Lake Shore Drive, Chicago, IL 60605-2496, USA. E-mail address: sleavitt@fieldmuseum.org (S.D. Leavitt). 1055-7903/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.ympev.2011.03.020 proven difficult to assess between and within taxonomic groups (Blanco et al., 2004b; Crespo et al., 2007; LaGreca and Lumbsch, 2001; Lumbsch and Schmitt, 2001; Ott et al., 2004). A growing body of evidence suggests that in many cases lichen species diversity has been underestimated, especially within morphologically similar species with cosmopolitan distributions (Buschbom and Mueller, 2006; Crespo and Pérez-Ortega, 2009; Kroken and Taylor, 2001; O’Brien et al., 2009; Printzen, 2009; Wedin et al., 2009; Wirtz et al., 2008). Species represent fundamental units of analysis in various subdisciplines of biology, so accurate species diagnoses are critical, and reassessing current species delimitation is particularly relevant in lichenized fungi. Finding and applying the appropriate character sets and analytical tools is one of the greatest challenges with empirical species delimitation in lichen-forming fungi (Crespo and Pérez-Ortega, 2009; Wirtz et al., 2008). In spite of the complicated issues associated with attempts to empirically circumscribe species, contemporary species concepts share the Author's personal copy 588 S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 common view that species are segments of separately evolving metapopulation lineages (de Queiroz, 1998, 1999, 2007; Mayden, 1999). This concept allows researchers to delimit species using different empirical properties and facilitates the development of new methods to test hypotheses of lineage separation (de Queiroz, 2007). A rapidly growing interest in species delimitations has resulted in novel approaches to investigate species boundaries (Carstens and Dewey, 2010; Knowles and Carstens, 2007; O’Brien et al., 2009; O’Meara, 2010; Sites and Marshall, 2004; Vieites et al., 2009; Weisrock et al., 2010; Yang and Rannala, 2010), and more properties supporting putative lineages are associated with a higher degree of corroboration (de Queiroz, 2007). Although different data sets and operational criteria may give conflicting or ambiguous results due to multiple evolutionary processes occurring within and between populations, the use of several independent suites of characters, such as morphology, geographic range, host preference, chemistry, and cross-validation using inferences from multiple empirical operational criteria have been shown to establish robust species boundaries (Dayrat, 2005; Duminil et al., 2006; Hey et al., 2003; O’Brien et al., 2009; Roe and Sperling, 2007; Ruiz-Sanchez and Sosa, 2010). As traditional characters used to delimit lichen species may misrepresent mycobiont diversity, we feel it is important to address lichen species boundaries using multiple independent datasets and operational criteria to effectively identify and delimit lichen species. We selected the Rhizoplaca melanophthalma speciescomplex (Ascomycota, Lecanorales, Lecanoraceae) as a model system to assess species diversity for this study because of its broad ecological and geographical distribution, morphological, chemical and genetic diversity, and importance as a sensitive indicator of environmental health (Arup and Grube, 2000; Aslan et al., 2004; Dillman, 1996; Leuckert et al., 1977; Ugur et al., 2004; Zhou et al., 2006). This group has previously been identified as a wellsupported monophyletic lineage and includes the placodiod crustose taxon, Lecanora novomexicana H. Magn., the umblicate taxon R. melanophthalma (DC.) Leuckert & Poelt, and at least 4 vagrant, or obligatory unattached, species (Arup and Grube, 2000). R. melanophthalma sensu lato (s. l.) has a worldwide distribution, and in North America it ranges from the northern boreal zone to Mexico along the Rocky Mountain corridor. It is commonly found in the Intermountain Western United States growing in large populations on rocky substrates. Specimens are generally umbilicate (fixed to the substrate by a single point of attachment), but often appear squamulose or pulvinate (polyphyllous), and considerable chemical variation is found within the species (McCune, 1987; Ryan, 2001). However, taxonomic status of lichens with distinct morphologies and chemotypes within R. melanophthalma s. l. and closely related taxa within this group (including the placodiod L. novomexicana and vagrant Rhizoplaca) remains uncertain (Arup and Grube, 2000; McCune, 1987). The vagrant taxa are endemic to the high plains and mountains of the central and northern Rocky Mountains in western North America and are particularly susceptible to habitat fragmentation, altered fire dynamics, and agricultural conversion (McCune and Rosentreter, 2007; Rosentreter, 1993). Speciation in lichenized fungi is, in general, understudied, and we present our analyses of the R. melanophthalma species-complex to represent the larger focus of this study, which is robust species delimitation in morphologically cryptic lichenized fungi. In this study we followed the general lineage concept (GLC; de Queiroz, 1998, 1999) as our non-operational species definition. We analyzed molecular data within a phylogenetic framework to identify candidate species by examining monophyletic groups recovered in the topology, and assessed the putative lineages across individual gene trees to identify lineages that exhibited genealogical exclusivity, an expected pattern for divergent lineages (Avise and Ball, 1990; Baum and Shaw, 1995; Hudson and Coyne, 2002). Candidate species were also evaluated within a population-level framework to assess gene flow and genetic differentiation (O’Brien et al., 2009), and we used multi-locus sequence data to identify genetic clusters without a priori assignment of individuals (Groeneveld et al., 2009; Weisrock et al., 2010). The phylogenetic relationships among candidate species circumscribed in this study were also inferred within a coalescence framework (Heled and Drummond, 2010) and support of speciation events was assessed using a Bayesian species delimitation test (Yang and Rannala, 2010). Finally, we investigated patterns in morphological and chemical variation and geographical and ecological distributions for each candidate species. 2. Materials and methods 2.1. Taxon sampling We analyzed sequence data from 170 individual posy rock lichens. The focal group was represented by the four know species from the R. melanophthalma species-complex: L. novomexicana (6 from 4 localities); R. melanophthalma, including R. melanophthalma subsp. melanophthalma (DC.) Leuckert & Poelt (127 vouchers from 37 localities), R. melanophthalma subsp. cerebriformis Rosentreter & B. D. Ryan (1), R. melanophthalma subsp. crispa Rosentreter & B. D. Ryan (1); Rhizoplaca haydenii, including R. haydenii ssp. arbuscula Rosentreter (2 from a single locality) and R. haydenii spp. haydenii (Tuck.) W. A. Weber (6 from 4 localities); Rhizoplaca idahoensis Rosentreter & McCune (4 from 2 localities); and two species not formally described (see McCune and Rosentreter, 2007), Rhizoplaca cylindrica (1) and Rhizoplaca subidahoensis (1). Fig. 1 depicts the high degree of morphological variation within the sampled R. melanophthalma species-complex in western North America. The present study emphasized umbilicate saxicolous forms; therefore sampling of the lobate taxon L. novomexicana and vagrant taxa was relatively limited. Collections of R. melanophthalma s. l. were initially made in 1997 at ten, 9  15 m plots along an altitudinal gradient (2200–3400 m) at Thousand Lakes Mountain (TLM), Wayne County Utah, USA (Porter, 1998), and three additional 9  15 m plots (2200 m, 2800 m, and 3300 m) were collected on the neighboring Boulder Mountain Plateau (BM), Wayne and Garfield Counties, Utah, in 2008. Seven individual thalli were randomly chosen from collections made at each plot to assess ecological trends in distributions and reproductive isolation between candidate species identified in this study (see Section 3.3). We also sampled 39 additional specimens from the R. melanophthalma species-complex, collected from 24 populations throughout the Intermountain West, USA. Available internal transcribed spacer sequences obtained from GenBank, representing 20 individuals, were included to assess relationships within a broader taxonomic and phylogeographic context. Rhizoplaca subdiscrepans (Nyl.) R. Sant. and Rhizoplaca chrysoleuca (Sm.) Zopf were selected as outgroups (Arup and Grube, 2000; Cansaran et al., 2006; Zhou et al., 2006). Collection information for all included specimens is summarized in Supplementary Table S1, and new voucher material generated for this study is housed at the Brigham Young University Herbarium of Nonvascular Cryptogams (BRY), Provo, Utah, USA. 2.2. Molecular data and sequence alignment Total genomic DNA was isolated using either the E.Z.N.A. Plant DNA Kit (Omega Bio-Tek, Norcross, GA), following manufacturer’s instructions, or the Prepease DNA Isolation Kit (USB, Cleveland, OH), following the plant leaf extraction protocol. We generated new sequence data via polymerase chain reaction (PCR) for three Author's personal copy 589 S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 Fig. 1. Variation in morphology and habit within the Rhizoplaca melanophthalma species-complex (Lecanoraceae) in western North America: (A) the lobate, placodioid taxon Lecanora novomexicana; (B) Rhizoplaca melanophthalma sensu lato (s.l.), with distinct light colored, pruinose apothecia discs; (C) Rhizoplaca melanophthalma sensu lato (s.l.), umblicate form with distinct lobes and dark apothecia; (D) R. melanophthalma s.l., umblicate form lacking lobes with pruinose apothecia (E) R. melanophthalma s.l., erratic form completely lacking umbilicus growing free on soil from western Idaho, with apothecia. Images F–I vagrant taxa endemic to the high plains and mountains of the northern Rocky Mountains (F) R. melanophthalma ssp. crispa; (G) R. idahoensis; (H) R. haydenii; (I) R. haydenii ssp. arbuscula. Images A–C, F–I were provided by S. Sharnoff. Table 1 Primers used for PCR amplification and sequencing of the nuclear ribosomal IGS, ITS, and group I intron markers and nuclear markers b-tubulin and MCM7. Marker Primer name Forward primer sequence Annealing temperature (°C) Reference IGS IGS12 NS1R ITS1F ITS4 Bt3-LM Bt10-LM Bt_rhizo_F Bt_rhizo_R Mcm7–709for Mcm7–1348rev LecMCM7f LecMCM7r 50 -AGTCTGTGGATTAGTGGCCG-30 50 -GAGACAAGCATATGACTAC-30 50 -CTTGGTCATTTAGAGGAAGTAA-30 50 - TCCTCCGCTTATTGATATGC-30 50 -GAACGTCTACTTCAACGAG-30 50 -TCGGAAGCAGCCATCATGTTCTT-30 50 -GCAACAAGTATGTTCCTCGTGC-30 50 -GTAAGAGGTGCGAAGCCAACC-30 50 -ACIMGIGTITCVGAYGTHAARCC-30 50 -GAYTTDGCIACICCIGGRTCWCCC AT-30 50 -TACCANTGTGATCGATGYGG-30 50 -GTCTCCRCGTATTCGCATNCC-30 66–56 (touchdown) Carbone and Kohn (1999) Carbone and Kohn (1999) Gardes and Bruns (1993) White et al. (1990) Myllys et al. (2001) Myllys et al. (2001) Leavitt et al. (2011) Leavitt et al. (2011) Schmitt et al. 2009 Schmitt et al. 2009 This study This study ITS/group I intron b-tubulin MCM7 nuclear ribosomal loci and two low-copy protein-coding loci (Table 1). The nuRNA gene tandem repeat exists in large copy numbers facilitating amplification of ribosomal loci from older specimens (Thousand Lake Mountain collections, 1997). Although low levels of intragenomic variation in fungal rDNA repeats suggest convergent evolution in which homogenization is very rapid and effectively maintains highly similar repeat arrays through recombination (Ganley and Kobayashi, 2007), previous studies have confirmed the utility of the sampled ribosomal loci for speciesand population-level studies in lichenized ascomycetes (Blanco et al., 2004a,b; Brunauer et al., 2007; Buschbom and Mueller, 2006; Gutiérrez et al., 2007; Kroken and Taylor, 2001; Lindblom and Ekman, 2006; O’Brien et al., 2009; Thell, 1999; Wedin et al., 2009; Wirtz et al., 2008). A gene duplication of b-tubulin has 55–60 55–60 66–56 (touchdown) 56 66–56 (touchdown) occurred within Ascomycota, however the paralogs are easily distinguishable within the analyzed group, and the marker has been successfully employed to investigate species-level relationships in other lichenized ascomycetes (Buschbom and Mueller, 2006; O’Brien et al., 2009; Wedin et al., 2009). PCR cycling parameters used for amplifying the ITS, group I intron, and b-tubulin loci followed the methods of Blanco et al. (2004b) and Wedin et al. (2009); cycling parameters for amplifying the IGS followed a 66–56° touchdown reaction (Lindblom and Ekman, 2006); and PCR cycling parameters amplifying the MCM7 fragment followed Schmitt et al. (2009). PCR products were quantified on 1% agarose gel and stained with ethidium bromide. In cases where no PCR product was visualized for the b-tubulin and MCM7 loci, internally nested PCR reactions were performed Author's personal copy 590 S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 using 0.3 ll of PCR product from the original reaction and newly developed internal primers ‘BT-RhizoF’ and ‘BT-RhizoR’ for the tubulin fragment, and ‘LecMCM7f’ and ‘LecMCM7r’ for the MCM7 fragment. Nested PCR reactions followed the touchdown PCR cycling parameters described above used to amplify the IGS fragment. PCR fragments were cleaned using the PrepEase PCR Purification Kit (USB, Cleveland, OH), and complementary strands were sequenced using the same primers used for amplification. Sequencing reactions were performed using the Big Dye3 (Applied Biosystems, Foster City, CA), and products were run on an AB 3730xl automated sequencer at the DNA Sequencing Center, Brigham Young University Provo, Utah, USA. Sequences were assembled and edited using Sequencher version 3.1.1 (Gene Codes Corporation, Ann Arbor, MI), and sequence identity was confirmed with ‘megaBLAST’ search in GenBank (Wheeler et al., 2006). Sequences were aligned in Muscle version 3.6 (Edgar, 2004), using default settings. Because recombination within nuclear genes can lead to errors in estimated relationships, we tested for intragenic recombination events in all sampled loci using GENECONV (Padidam et al., 1999), Bootscan/Rescan (Martin et al., 2005a), MaxChi (Smith, 1992), Chimaera (Posada and Crandall, 2001a), SiScan (Gibbs et al., 2000), LARD (Holmes et al., 1999), Topal/DSS (McGuire and Wright, 2000), and 3Seq (Boni et al., 2007) methods implemented in Recombination Detection Program v3.44 (Heath et al., 2006; Martin et al., 2005b). Additionally, we searched alignments for evidence of recombination breakpoints using GARD (Kosakovsky Pond et al., 2006), implemented in www.datamonkey.org web server. 2.3. Nucleotide polymorphism analyses and population differentiation We used DnaSP 5.10 (Librado and Rozas, 2009) to calculate basic nucleotide polymorphism statistics, including numbers of haplotypes (H), total number of polymorphic sites (Npoly), average pair-wise diversity per site, (p; Nei, 1987) for each candidate species (see Section 3.3). In addition, genetic differentiation between candidate species was assessed by calculating FST values using DnaSP and counting the number of fixed nucleotides for all pairwise comparisons (O’Brien et al., 2009). F-statistic calculations were estimated from specimens with complete ITS, IGS, ß-tubulin, and MCM7 dataset (the group I intron was missing in all specimens assigned to a single candidate species, and this marker was therefore excluded from FST calculations). Aligned sequences were scanned for fixed characters between each candidate species and the remaining data matrix in DnaSP, and the total number of fixed nucleotide positions was tabulated for each candidate species. 2.4. Phylogenetic analyses Heterogeneity in phylogenetic signal among the sampled markers was assessed before combining the datasets. We performed maximum likelihood (ML) analyses of the concatenated ribosomal dataset (ITS, IGS, and group I intron), b-tubulin, and MCM7 markers separately in RAxML version 7.0.4 (Stamatakis, 2006; Stamatakis et al., 2008), using the ‘rapid bootstrapping’ option as implemented on the CIPRES Web Portal. We used the GTRGAMMA model, which includes a parameter (C) for rate heterogeneity among sites, and chose not to include a parameter for estimating the proportion of invariable sites because C mathematically account for this source of rate heterogeneity by using 25 rate categories (Stamatakis, 2006). Support values for ribosomal, b-tubulin, and MCM7 phylogenies were examined for well-supported (P70%) conflicts between data sets (Lutzoni et al., 2004). Phylogenetic relationships were estimated from the combined data set using mixed-model Bayesian inference (BI) in Mr.Bayes version 3.1.2 (Huelsenbeck and Ronquist, 2001), using the default priors, namely a flat Dirichlet prior for the relative nucleotide frequencies and rate parameters, a discrete uniform prior for topologies, and an exponential distribution (mean 1.0) for the gammashape parameter and branch lengths. We used MrModeltest version 2.3 (Nylander et al., 2004) to identify the best-fitting model of evolution for each marker using the Akaike Information Criterion (AIC; Posada and Crandall, 2001b), and we treated each marker as a separate partition. Four independent replicate searches were executed with eight chains; each run started with randomly generated trees and consisted of sampling every 1000 generations for 20,000,000 generations. To evaluate stationarity and convergence between runs, log-likelihood scores were plotted using TRACER version 1.5 (Rambaut and Drummond, 2007), ESS statistics, and the average standard deviation in split frequencies were assessed (Hall, 2007). Trees generated prior to stationarity were discarded as burn-in (Huelsenbeck et al., 2001). The results were summarized with a majority-rule consensus tree from the remaining trees from the four independent runs. Bayesian posterior probabilities (PP) were assessed at all nodes, and clades with PP P 0.95 were considered strongly supported (Huelsenbeck and Rannala, 2004). We conducted a ML analysis using RAxML 7.0.4, allowing a separate GTR model with unique parameter values for each locus (Stamatakis, 2006; Stamatakis et al., 2008). We used the GTRGAMMA model and did not include a parameter for the proportion of invariable sites (Stamatakis, 2006). A search combining 200 separate maximum likelihood searches (to find the optimal tree) and 1000 ‘‘fastbootstrap’’ replicates to evaluate nodal support was conducted on the complete dataset. Exploratory phylogenetic reconstructions including GenBank specimens represented solely by ITS sequences in the concatenated dataset resulted in reduced nodal support across the topology and important ambiguous relationships. Therefore we chose not to include GenBank accessions in the multilocus data matrix in order to minimize the effect of missing data (Baurain et al., 2007). To assess relationships within a broader geographic context we reconstructed the ITS gene tree using both BI and ML inference from all available ingroup ITS sequences, including 20 sequences retrieved from the GenBank database, with Rhizoplaca chrysolueca selected as the outgroup (Arup and Grube, 2000; Zhou et al., 2006). BI and ML reconstructions were performed for the complete ITS dataset as described above. Establishing a preliminary set of species boundaries using concatenation provides a reasonable starting point for identifying cryptic species and species tree inference (Le Gac et al., 2007; Leache, 2009; Šlapeta et al., 2006). We used the topology from the concatenated dataset to guide the identification of candidate species for this study. We chose to circumscribe a total of 10 candidate species to represent four currently accepted taxa and six phylogenetic lineages identified within the topology representing R. melanophthalma s. l (see Results, Section 3.3). Following recommendations of Sites and Marshall (2004) and de Queiroz (2007), we implemented multiple analytical approaches to assess species boundaries for independent corroboration of the candidate species identified in the current study. 2.5. Haplotype network reconstructions and genealogical concordance Relationships of candidate species were evaluated between individual haplotype networks to identify lineages that exhibited genealogical exclusivity across multiple loci (Avise and Ball, 1990; Hudson and Coyne, 2002). The presence of the same clades in the majority of single-locus genealogies is taken as evidence that the clades represent reproductively isolated lineages (Dettman et al., 2003; Pringle et al., 2005). We used statistical parsimony to assess the genealogical relationship of every individual and Author's personal copy S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 compare relationships of candidate species between genes. Networks were constructed under a 95% parsimony probability criterion (Templeton et al., 1992) from concatenated ribosomal sequences (ITS, IGS, intron), the b-tubulin, and the MCM7 fragments using the program TCS v1.21 (Clement et al., 2000). Gaps were treated as missing data for the ribosomal network reconstruction to include voucher specimens missing one of the three ribosomal loci. All protein-coding sequences were trimmed to the length of the fragment resulting from nested PCR reactions and a single sequence missing approximately half the fragment was removed from the -tubulin network analysis. All network uncertainties (i.e. closed loops) were treated following Templeton and Sing (1993). 2.6. Bayesian population structure analysis Individual-based approaches provide an alternative for identifying population structure and barriers to gene flow (Saisho and Purugganan, 2007), as analyses based on predefined groups may obscure patterns of differentiation (Latch et al., 2006; Rowe and Beebee, 2007). We used a Bayesian population assignment test implemented in STRUCTURE version 2.32 (Falush et al., 2003; Pritchard et al., 2000) to infer population structure based on a combined genotypic matrix from SNPs data from all loci (ITS, IGS, group I intron, b-tubulin, and MCM7), without using known geographic location or putative species classification of the individual as priors. An admixture model was used with correlated allele frequencies. We implemented 15 replicate runs for each number of assumed populations (K), with a range of K from 1 to12. Based on preliminary runs, all analyses used 30,000 MCMC generations to estimate the posterior distribution following a burn-in period of 15,000 generations. In some cases, independent runs for K values 3–12 appeared to converge on different parameter space, and longer burn-in or MCMC runs did not improve convergence. Therefore, we calculated the median log (ln) likelihood of each K value from the four best-scoring runs. We calculated the modal value (DK) based on the second order rate of change of the likelihood function between successive K values (Evanno et al., 2005). In some scenarios the DK method may fail detect the true number of clusters (K), and we examined subgroups created by the best individual assignments produced by STRUCTURE (Evanno et al., 2005; Groeneveld et al., 2009; Saisho and Purugganan, 2007; Weisrock et al., 2010). 2.7. Species tree inference and speciation probabilities We used the hierarchical Bayesian model implemented in BEAST v. 1.6.1 (Heled and Drummond, 2010) to estimate the species trees for the ten candidates species. ⁄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 (Heled and Drummond, 2010). We selected 2–8 individuals from each candidate species, representing the sampled genetic diversity. The group I intron was missing from all individuals recovered in clade IVd (Fig. 2) and was therefore excluded from the analysis. We specified the best-fitting models of molecular evolution for each locus using the same models used in the MrBayes v. 3.1.2 analysis under an uncorrelated relaxed lognormal clock (Drummond et al., 2006). We selected the lineage birth/death model and constant population size for species tree priors. Two independent MCMC analyses were run for a total of 50 million generations (sampling every 1000 steps and excluding the 10 million generations of each run as burnin). We assessed convergence by examining the likelihood plots through time using Tracer v. 1.5 (Rambaut and Drummond, 2007), and the effective sample sizes (ESS) of parameters of ⁄ 591 interest were all above 200. Posterior probabilities of nodes were computed from sampled trees after burn-in. We estimated the marginal posterior probability of speciation using the program BPP (Rannala and Yang, 2003; Yang and Rannala, 2010). This method accommodates the species phylogeny as well as lineage sorting due to ancestral polymorphism. We used a gamma prior h  G(2, 1000) on the population size parameters with algorithm 0. The age of the root in the species tree was assigned the gamma prior s0  G(2, 1000), while the other divergence time parameters were assigned the Dirichlet prior (Yang and Rannala, 2010). Because the prior distribution of h and s0 can result in strong support for models containing more species (Leaché and Fujita, 2010), we also explored a more conservative combination of priors that should favor fewer species by assuming large ancestral population sizes h  G(1, 10) and relatively shallow divergences among species s0  G(2, 2000). The species tree estimated in the ⁄BEAST analysis, representing the ten candidate species, was used as the fully resolved guide tree. Running the rjMCMC analysis for 1000,000 generations with a burn-in of 100,000 produced consistent results across independent analyses initiated with different starting seeds and species trees. Each analysis was run at least twice to confirm consistency between runs. 2.8. Morphological and biochemical comparisons Considering current studies (Arup and Grube, 2000; Cansaran et al., 2006; Ryan, 2001; Zheng et al., 2007; Zhou et al., 2006), a total of 14 morphological characters were quantified from all specimens to identify potentially diagnostic characters for candidate species identified in this study, including: point of attachment (distinctly umbilicate/squamulose), thallus form (polyphyllous/monophyllous), lobe morphology (distinct/intermediate/indistinct), appearance of upper surface (dull/shiny), upper surface texture (smooth/cracked), upper surface color (light to moderately greenish yellow/olive), lower surface texture (smooth/rough), lower surface edges (black near edges/not blackened edges), lower surface color (tan/brown), apothecia morphology (sessile/basally constricted), apothecia pruinosity (heavily pruinose/moderately pruinose/not pruinose), thallus margin (entire/crenate), spore shapes (ellipsoid/subglobose), spore size (continuous character). Lichen compounds were extracted from 0.02 g liquid nitrogenground specimens overnight in acetone at 4 °C. The supernatant was removed, dried, reconstituted in methanol, and analyzed using high-performance liquid chromatography (HPLC). Retention index values (RI) were calculated from benzoic acid and solorinic acid controls (Feige et al., 1993; Lumbsch, 2002). For HPLC, we used an Agilent Technologies 1200 series integrated system with a Zorbax Eclipse XDB8-CB column (4.6  150 mm, 5 lm) regulated at 30 °C, spectrometric detectors operating at 210, 254, 280, 310 nm, and a flow rate of 0.7 ml/min. Following established protocols (Feige et al., 1993; Lumbsch, 2002), two mobile phases, A and B, were used: 1% aqueous orthophosphoric acid (A) and methanol (B). The run started with 30% B for 1 min and was raised to 70% B within 15 min of the start time, then to 100% B during an additional 15 min, followed by isocratic elution in 100% B for the final 20 min. Mobile phase B was decreased to 30% within 1 min and the column was flushed with 30% B for 15 min following each run. UV spectra of each peak were recorded and computermatched against a library of ultraviolet spectra from authentic metabolites derived under identical conditions using Agilent Chemstation software. The correlation of UV spectra with the standards in the library was greater than 99.9% for each substance identified. When multiple library entries matched with this level of identity, calculated R/I values were used to discriminate between compounds. Author's personal copy 592 S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 Fig. 2. Relationships among sampled specimens collected from the Rhizoplaca melanophthalma group inferred from a maximum likelihood analysis of ribosomal and nuclear DNA sequence data (2600 bp, ITS, IGS, group I intron, b-tubulin, and MCM7). Values at each node indicate non-parametric-bootstrap support/posterior probability. Only support indices PBS 50/PP 0.50 are indicated. Clade numbers plotted to the right of the tree indicate candidate species. GenBank accessions represented solely by ITS sequences were not included. 3. Results For this study, 627 new sequences were generated, including 152 ITS, 129 IGS, 75 group 1 intron, 137 b-tubulin, and 134 MCM7 sequences. The combined data matrix consisted of 2639 aligned nucleotide position characters (Table 2; TreeBase Accession 11094). Missing data were generally limited to the outgroup taxa R. chrysoleuca and R. subdiscrepans. However, we were unable to generate group I intron sequences from all specimens recovered in clade IVd and single individuals from R. haydenii subsp. arbuscula (092f), R. idahoensis (093f) and clade II (693f) (defined below). All representative haplotypes of the five gene fragments have been deposited in GenBank under Accession Nos. HM576889– HM577515 (Supplementary Table S2). No evidence for recombination within any of the nuclear markers was detected. 3.1. Polymorphism statistics and estimates of genetic differentiation Polymorphism statistics are reported in Table 3. The greatest nucleotide diversity for candidate species was generally recovered Author's personal copy 593 S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 IV). Group I intron sequences were missing for all individuals assigned to clade IVd. for ribosomal loci. High levels of genetic differentiation between all pairs of candidate species were calculated from the combined data set, as measured by FST (Table 4). Fixed differences between candidate species circumscribed in this study were identified from ribosomal markers for all pair-wise comparisons, and fixed differences were identified in at least one of the protein-coding fragments for 40 of 45 pairwise comparisons (Table 4). The ribosomal data matrix showed the greatest number of fixed character differences between each candidate species compared to all remaining lineages; while the protein-coding matrixes generally did not reveal fixed character differences (Table 4). However, the b-tubulin fragment revealed nine fixed nucleotide positions in clade I and 1 fixed locus in clade IVb, and the MCM7 fragment revealed two fixed nucleotide positions in clade I and five fixed characters in R. idahoensis (clade 3.2. Phylogenetic reconstructions We detected limited discordance between the ribosomal, b-tubulin, and MCM7 topologies restricted to clades with relatively shallow evolutionary histories. Conflicting terminals are shown in individual gene trees (Supplementary Figs. 1–3). Given the overall congruence, the ribosomal, b-tubulin, and MCM7 gene regions were combined to maximize the total number of characters for phylogenetic analyses and branch length estimation (Rokas et al., 2003; Wiens, 1998). Major relationships and branch support values from a restricted dataset excluding the individuals showing incongruent positions when comparing the single gene phylogenies were highly similar to the phylogeny resulting from the complete dataset. The partitioned Bayesian analyses, summed from four independent runs, yielded a consensus tree with a negative ln harmonic mean of 11,092.49. All parameters converged within the first 25% of sampled generations, leaving a posterior distribution estimated from 15,000 trees per run (60,000 total post-burn-in sampled trees). The partitioned ML analysis yielded a single best scoring tree lnL = 10,755.758. As the recovered trees were similar across methods and the topologies did not show any strongly supported conflict; we present the results of the ML analysis with ML bootstrap (BS) and Bayesian posterior probability (PP) values in Fig. 2. The R. melanophthalma group is strongly supported as monophyletic and several other well-supported groups can be identified in the tree. Table 2 Genetic variability of sampled markers used in this study, including alignment length (number of basepairs); variable and parsimony-informative (PI) sites for each sampled locus; and locus-specific model of evolution identified using the Akaike information criterion in MrModeltest. Numbers in parentheses indicate the number of variable and parsimony-informative sites for the Rhizoplaca melanophthalma speciescomplex only. Locus Length ITS IGS Group I intron b-tubulin MCM7 Total # Variable sites # PI sites 561 374 269 819 616 163 (91) 138 (84) 98 (44) 165 (90) 158 (123) 127 (57) 103 (54) 84(30) 132(55) 123 (42) 2639 722 (432) 569 (238) Model selected GTR + G GTR + I SYM + G HKY + I + G GTR + G – Table 3 Polymorphism statistics for candidate species within the R. melanophthalma species-complex. N, number of individuals sampled, Npoly, number of polymorphic sites; h, number of unique haplotypes; p, estimate of 4 Nl per base pair using the average pairwise differences. ITS Clade Clade Clade Clade Clade Clade Clade Clade Clade Clade I (L. novomexicana) II III IV (R. haydenii) IV (R. h. spp. arbuscula) IV (R. idahoensis) IVa IVb IVc IVd Total IGS Intron MCM7 b-tubulin N/Npoly/h p N/Npoly/h P N/Npoly/h P N/Npoly/h p N/Npoly/h p 3/0/1 24/35/17 13/5/5 5/6/4 2/1/2 3/3/2 3/3/3 14/9/7 5/1/2 55/11/10 0 0.00930 0.00188 0.00475 0.00182 0.00367 0.00427 0.00235 0.00088 0.00162 4/0/1 21/37/18 13/1/2 4/4/4 2/1/2 3/1/2 3/2/3 13/3/4 5/3/3 55/19/18 0 0.01776 0.0014 0.00318 0.00272 0.00272 0.00363 0.00265 0.00381 0.01191 3/0/1 23/19/17 13/0/1 5/4/3 1/0/1 2/0/1 3/0/1 14/3/4 5/0/1 0/na/na 0 0.1089 0 0.00723 na 0 0 0.00327 0 na 4/2/3 24/34/17 13/3/2 5/2/2 1/0/1 3/4/2 3/0/1 13/9/9 5/5/3 55/32/8 0.00146 0.01430 0.00067 0.00117 na 0.0039 0 0.00285 0.00439 0.00266 2/11/2 23/10/8 13/4/4 5/7/2 2/0/1 37316 3/0/1 13/19/6 5/2/2 55/5/6 0.02041 0.00278 0.00157 0.00779 0 0.00124 0 0.01308 0.00148 0.00040 127/91/52 0.02221 122/84/54 0.02494 69/43/27 0.03521 127/71/40 0.01309 126/112/33 0.01486 Table 4 Fixed differences and fixation indices (FST) for all pairwise comparisons of candidate species identified within R. melanophthalma species-complex. Numbers across the top row correspond to candidate species numbers in the first column. Numbers of fixed differences (ribosomal/b-tubulin/MCM7 characters) are represented for all comparisons below the diagonal and FST values are represented above the diagonal. The last column indicates total number of fixed nucleotides identified between each candidate species and the remaining data matrix. Numbers within parentheses represent fixed ribosomal characters/fixed protein-coding characters. Accessions representing R. haydenii subspecies arbuscula were not included in FST calculations because of the small sample sizes and pairwise comparisons are not represented. Candidate species 1 2 3 4 5 6 7 8 9 10 Fixed characters 1. Clade I (L. novomexicana) 2. Clade II 3. Clade III 4. Clade IV (R. haydenii) 5. Clade IV (R. h. spp. arbuscula) 6. Clade IV (R. idahoensis) 7. Clade IVa 8. Clade IVb 9. Clade IVc 10. Clade IVd – 49 (31/13/5) 77(55/18/4) 77(51/20/6) 82(54/19/9) 71(53/8/10) 65(38/20/7) 76(54/19/3) 76(51/18/7) 61(36/18/7) 0.77102 – 32(28/0/4) 32(26/0/6) 36(28/1/7) 33(28/0/5) 36(29/0/7) 31(29/2/0) 35(28/0/7) 22(16/0/6) 0.89534 0.75792 – 55(36/11/8) 56(39/9/8) 55(38/8/9) 54(36/10/8) 48(39/7/2) 55(39/8/8) 45(29/8/8) 0.86359 0.732 0.90524 – 7(2/4/1) 12(1/0/11) 27(21/5/1) 13(4/8/1) 18(14/3/1) 10(6/3/1) na na na na – 15(2/0/13) 27(24/3/0) 11(5/6/0) 6(6/0/0) 8(8/0/0) 0.85763 0.69564 0.89291 0.58915 na – 38(23/2/13) 15(5/5/5) 18(5/0/13) 14(6/0/12) 0.88863 0.76148 0.9382 0.82339 na 0.84298 – 30(23/7/0) 24(24/0/0) 14(14/0/0) 0.85574 0.72461 0.88716 0.67851 na 0.6808 0.82136 – 9(4/5/0) 13(8/5/0) 0.88172 0.75139 0.9273 0.66667 na 0.71146 0.83333 0.67031 – 7(7/0/0) 0.88085 0.7426 0.92874 0.71894 na 0.75427 0.80228 0.72953 0.66841 – 32(21/11) 3(3/0) 15(15/0) 1(0/1) 0 (0/0) 7(1/6) 7(7/0) 3(2/1) 1(1/0) 1(1/0) Author's personal copy 594 S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 The ITS topology (Fig. 3) recovered most lineages identified in the combined analyses. GenBank accessions representing individuals collected in Austria (AF159935), China (AY509791, EF095286, and EF095297), and the United States (AF159929-Arizona and AF159935-Arizona) were recovered in a well-supported clade (91/1.0) corresponding to clade II identified in the combined analyses. Six accessions collected in China (EF095278, EF095280, EF095283, EF095285, EF095287, and EF095290) were recovered within a well-supported clade (81/0.98) corresponding to clade IVb from the combined analyses, and two accession representing R. cerebriformis (AF159942, Idaho, USA) and R. subidahoensis (AF159944, Idaho, USA) were recovered within a well-supported clade (90/1.0) corresponding to clade IVa from the combined analyses. A single accession representing R. cylindrical (AF159941, Idaho, USA) was recovered in a clade with high ML bootstrap support (82) and weak PP support (0.79) corresponding to clade IVd in the combined analyses. Two vagrant accessions representing R. idahoensis (AF159943-Idaho, USA) and R. haydenii subsp. haydenii (AF159937-Idaho, USA) were recovered in a well-supported clade (85/1.0) containing individuals all assigned to clades IVb, IVc, R. haydenii, R. haydenii ssp. arbuscula, and R. idahoensis in the combined analyses. L. novomexicana was recovered as polyphyletic in two well-supported lineages; one containing specimens collected in northeastern Utah, and the second (clade V, Fig. 3) in two Fig. 3. The maximum likelihood ITS topology obtained from all sampled specimens and available GenBank accessions collected from the Rhizoplaca melanophthalma speciescomplex. Values at each node indicate non-parametric-bootstrap support/posterior probability. Only support indices PBS 50/PP 0.50 are indicated. Clade numbers plotted to the right of the tree indicate lineages corresponding to candidate species shown in Fig. 2. Author's personal copy S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 GenBank accessions, one from Arizona (AF159923) and the other from New Mexico (AF159923). However, the relationship between the L. novomexicana lineages lacked strong statistical support. 3.3. Candidate species We circumscribed 10 candidate species based on the results from our phylogenetic reconstructions and current taxonomic boundaries for additional empirical testing of species boundaries. However, individuals represented by GenBank accessions were not included because they were represented by a single marker. L. novomexicana (clade I, Fig. 2) was recovered as a well-supported lineage (BS = 100/PP = 1.0), sister to the remaining R. melanophthalma taxa with weak nodal support. Clade II was recovered with high nodal support (95/1.0), and corresponds to a genetically and morphologically diverse assemblage of umbilicate saxicolous specimens collected throughout the intermountain western United States, all containing usnic and psoromic acids. However, the relationship of clade II with other well-supported sister lineages lacks strong nodal support (43/0.89). Clade III was also recovered with strong support (100/1.0), and is represented by umbilicate saxicolous individuals with little morphological or genetic variation collected from two plots (BM-3 and TLM-9) on the Aquarius Plateau in south central Utah, U.S.A. Clade III was recovered with strong nodal support (94/0.98) as sister to a fourth well-supported clade (99/ 1.0) containing a chemically diverse assemblage of umbilicate 595 and vagrant specimens (clade IV). Seven additional candidate species were identified within clade IV to accommodate currently described vagrant taxa and an exhaustive subdivision of the remaining specimens. All sampled vagrant taxa were recovered within a single monophyletic clade with weak nodal support (BS and PP < 50/ 0.50). R. idahoensis, R. haydenii subsp. haydenii, and R. haydenii subsp. arbuscula were treated as independent lineages based on current taxonomic circumspection. Both R. idahoensis and R. haydenii subsp. arbuscula were recovered as well-supported monophyletic lineages (94/1.0 and 81/1.0, respectively), while R. haydenii subsp. haydenii was found in two well-supported clades. A single saxicolous specimen with unique lobe morphology (715f) was recovered within the R. haydenii subsp. haydenii clade. In addition to the currently described vagrant taxa, four candidate species are circumscribed to accommodate exhaustive subdivision within the larger clade. Clade IVa (Fig. 2) was recovered with strong nodal support (100/1.0) and contains three morphologically and geographically diverse individuals. All specimens containing lecanoric or orsellinic acids were recovered within clade IVb with moderate to strong nodal support (BS = 83; PP = 0.93). Clade IVc (Fig. 2) was also recovered with strong support (82/ 1.0), and included five individuals; and clade IVd included the remaining 55 individuals. Although this lineage was recovered as monophyletic, it lacked strong support in the combined phylogenetic reconstructions. Fig. 4. Geographical distributions of candidate Rhizoplaca species in the Intermountain western USA. Colors refer to different lineages, indicated in key. Insert shows distributions of putative lineages along two altitudinal gradients in southern Utah, USA. A total of 7 individual were included from each plot and the proportion of candidate species recovered at each plot is represented. Author's personal copy 596 S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 Fig. 5. (A) Unrooted statistical parsimony haplotype networks at 95% probability of the ribosomal, MCM7, and b-tubulin loci representing relationship within the R. melanophthalma species-complex. Each candidate species is designated by a different color. Size of circles is proportional to the number of individuals of a given haplotype, and black dots represent inferred haplotypes not sampled. (B) Correspondence between candidate species identified from the combined maximum likelihood analysis and the population clusters identified using STRUCTURE. Numbers at nodes represent maximum likelihood bootstrap values and posterior probabilities, and relationships within candidate species are collapsed for ease of presentation (see Fig. 2 for detailed relationships). Candidate species are mapped to corresponding clusters in the STRUCTURE plot. Each population cluster is represented by a different color, and vertical bars within each cluster represent individuals and the proportion of a bar assigned to a single color represents the posterior probability that an individual is assigned to that cluster. The colors in the topology and STRUCTURE plot correspond to candidate species colors shown in Fig. 5A. phylogenetic hypothesis of relationships in the Rhizoplaca melanophthalma species-complex in western North America. Author's personal copy 597 S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 Geographic distributions of candidate species and the distribution of these species along altitudinal transects on Thousand Lakes Mountain and Boulder Mountain, Utah are summarized in Fig. 4. 3.4. Haplotype networks We recovered a total of five independent haplotype networks for the combined ribosomal data set, and two networks for both the -tubulin and MCM7 datasets (Fig. 5A). The ribosomal network haplotypes separated by up to 15 mutational steps had greater than 95% probability of being parsimoniously connected. In the b-tubulin and MCM7 distinct networks were connected by up to 11 or 10 mutational steps, respectively. For all markers, clade I (L. novomexicana) formed an independent network. In addition, clades II, III, and IVa formed independent networks constructed from the ribosomal dataset, while clades IVc, IVb, IVd, R. haydenii 3.5. Bayesian population structure The median ML values of the Bayesian clustering analysis using STRUCTURE with estimates of K = 1–12 are shown in Fig. 6A. These analyses reveal a general pattern of a plateau with a decrease in median maximum likelihood values above a K = 6 level. In contrast, the DK method indicates that a K = 2 model best fits the data (Fig. 6B; DK = 137.170 for K = 2; DK = < 25 for all other K values). However, the plateau in likelihood values around K = 6 suggest a -2000 B 140 -3000 120 -4000 ∆K value Likelihood (lnL) A subsp. arbuscula (clade IV), R. haydenii subsp. haydenii (clade IV), and R. idahoensis (clade IV), were found on a single network. In both the -tubulin and MCM7 datasets clades II, III, IVa, IVb, IVc, IVd, R. haydenii subsp. arbuscula (clade IV), R. haydenii subsp. haydenii (clade IV), and R. idahoensis (clade IV) were found on a single network. -5000 -6000 -7000 100 80 60 -8000 40 -9000 20 0 -10000 1 2 3 4 5 6 7 8 9 10 11 12 2 3 4 5 6 7 8 9 10 11 12 Inferred clusters (K) Inferred clusters (K) Fig. 6. Plots of calculations for K values 1–12 in STRUCTURE analysis of the combined dataset. (A) The mean log probability of the data for K = 1–12, calculated from the four best scoring runs for each K value. (B) DK values for K = 2–12. R. chrysoleuca clade I 1.0,1.0,1.0 clade III clade II 0.99,1.0,1.0 clade IVb 0.57,1.0,1.0 clade IVa 1.0,1.0,1.0 clade IVc 0.63,1.0,0.98 0.54,1.0,0.97 clade IVd 0.75,1.0,0.98 R. h. subsp. arbuscula 0.95,1.0,0.96 R. h. subsp. haydenii 0.54,0.95,0.62 R. idahoensis 0.02 Fig. 7. The coalescent-based species tree and Bayesian species delimitation results for the Rhizoplaca melanophthalma species-complex. Posterior probabilities at nodes indicate support from the BEAST analysis, speciation probabilities under prior means = 0.001, and prior means = 0.1. Author's personal copy 598 S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 Table 5 Chemotypic variation by candidate species in the R. melanophthalma species-complex based on HPLC analysis. Superscript number following acid nominal indicate acid occurrence: 1major or minor; 2major or not present; 3minor or not present; 4minor or trace; and 5trace or not present. Acid Clade I Clade II Clade III R. haydenii (clade IV) R. h. ssp. arbuscula (clade IV) R. idahonesis (clade IV) Clade IVa Clade IVb Clade IVc Clade IVd Usnic1 Psoromic2 Lecanoric2 Orsellinic3 Gyrophoric5 Constipatic3 Dehydroconstipatic3 Dehydroprotoconstipatic3 Subpsoromic acid3 20 -O-demethylsubpsoromic4 20 -O-demethylpsoromic3 1 1 0 0 0 0 0.25 0.25 0.25 0.75 0.75 1 0.91 0 0 0 0.91 0.91 0.7 0.43 0.52 0.39 1 1 0 0 0 0.64 0.36 0.36 1 1 0.82 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0.66 0 0 0 1 1 0.33 0 1 1 1 1 0.57 0.64 0.43 0.93 0.93 0.86 0.57 0.29 0.5 1 0.40 0 0 0 1 1 1 1 1 0.5 1 0.95 0 0 0 0.91 0.95 0.55 0.78 0.87 0.73 higher number of population clusters (Fig. 6A). A plot of individual membership coefficients for K = 6 reveals a high number of population clusters with average individual membership coefficients (i.e. posterior probabilities) greater than 0.9 (Fig. 5B). Population clusters inferred for K > 6 did not yield additional clusters with high membership coefficients. Therefore, we place our focus on K = 6 as an uppermost level of population structure. The K = 6 model is generally consistent with the identified candidate species. However, all vagrant species (R. haydenii subsp. haydenii, R. haydenii subsp. arbuscular, and R. idahoensis) were recovered within a single population cluster, along with all individuals assigned to clade IVc in the combined phylogenetic analysis. A total of three saxicolous accessions (554f, 556f, and 715F) and three erratic, or facultatively unattached, accessions (668f, 669f, 670f) were assigned to the cluster with vagrant taxa. Clades IVa and IVd were also recovered as a single population cluster; however, membership coefficients for individuals with posterior probabilities were <0.71 for clade IVa and P0.87 for clade IVd. 3.6. Species tree inference and speciation probabilities The species tree resulting from the analysis of the ten candidate species estimated with ⁄BEAST is shown in Fig. 7. Within the R. melanophthalma species-complex, clade I, representing L. novomexicana, was recovered with strong support (posterior probability >0.99) as the earliest diverging lineage. All candidate species recovered within clade IV in the concatenated analyses formed a wellsupported monophyletic clade and clade II was sister to this lineage, although without significant posterior support. Relationships among candidate species recovered within clade IV lack strong statistical support, although the three vagrant species formed a wellsupported clade (PP > 0.95), sister to clades IVa, IVc, and IVd. The Bayesian species delimitation results for the R. melanophthalma species-complex are shown in Fig. 7. Bayesian species delimitation supports the ten candidate species guide tree with speciation probabilities > 0.95 on all nodes. Exploratory analyses using different prior distributions for h and s0 did not affect speciation probabilities. 3.7. Morphology and chemistry We adopted the approach of Wiens and Penkrot (2002), suggesting that in order for characters to diagnose a lineage they must be invariant for alternative character states or show no overlap in trait values. Both vegetative morphology and reproductive characters, spore size and shape, were highly variable within some candidate species, and overall we were unable to identify morphological or reproductive characters corroborating candidate species following Wiens and Penkrot (2002). Occurrence of the 11 most common compounds identified in HPLC analyses within each lineage is summarized in Table 5. The majority of specimens belonged to the usnic/psoromic acids chemotype (119 specimens, including all specimens of L. n ovomexicana), having a broad geographical and ecological distribution; 9 specimens contained usnic, psoromic, and lecanoric acid; and 5 specimens contained usnic, psoromic, and orsellinic acid. All sampled vagrant specimens expressed usnic acid only. In addition to the previously reported psoromic acid, we found 2’-O-demethylsubpsoromic acid, 2’-O-demethylpsoromic acid, and the recently described b-orcinol depsidone, subpsoromic acid (Elix 2000). The dibenzofuran-derivative, usnic acid, was present in all samples, and some combination of the aliphatic acids – dehydroprotoconstipatic acid and constipatic acid – were present in all individuals, except the sampled vagrant taxa. We found gyrophoric (triorsellininc) acid and also the monocyclic-depside precursor, orsellinic acid, restricted to specimens assigned to clade IVb (defined in 3.3) in the combined molecular analyses, in addition to previous reports for lecanoric (diorsellinic) acid (Arup and Grube, 2000; McCune, 1987). 4. Discussion In this study, we assembled multiple lines of evidence to identify and delimit candidate species within the R. melanophthalma species-complex. Based on all of the available evidence, we circumscribe ten candidate species within this complex, and most fall within a nominal taxon currently recognized as a single cosmopolitan species, R. melanophthalma. Genetic patterns, generated by population-level processes operating within divergent lineages, provide an informative perspective about the process of speciation in the R. melanophthalma species-complex. Results of the empirical tests delimiting species are summarized in Table 6. Our results provide a compelling case of diversification within the R. melanophthalma species-complex using molecular data and multiple analytical tools, although most candidate species were not supported unambiguously by independent datasets. Besides the placodiod crustose taxon, L. novomexicana, we found that the greatest morphological and chemical variation was restricted to closely related lineages (sampled vagrant taxa and clades IVb and IVc), while morphological and chemical characters supporting more divergent groups were not identified. Ecological interactions are expected to drive phenotypic divergence during the early stages of lineage diversification when species richness is low and available niches are ‘‘open’’ (Schluter, 2000). The ecological transition from a saxicolous attached form to morphologically distinct vagrant forms appears to follow the ecological theory of adaptation (Funk et al., 2006). The STRUCTURE analysis assigned all vagrant forms to a single population cluster, suggesting a recent divergence Author's personal copy S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 599 Table 6 Summary of data supporting candidate species within the R. melanophthalma species-complex. Fixed characters, the total number of fixed nucleotide characters relative to the remaining data matrix (combined ribosomal loci-b-tubulin-MCM7); genealogical, candidate species recovered as an exclusive lineage in gene haplotype networks, ‘⁄0 indicate support from individual ribosomal, b-tubulin, and MCM7 network reconstructions; STRUCTURE, indicates if the candidate species was recovered as a unique population cluster in the Bayesian clustering analysis; BPP, posterior support for speciation events calculated using a Bayesian modeling approach; independent characters support, support from independent morphological or chemical data. Candidate species Fixed characters Genealogical exclusivity STRUCTURE BPP (PP > 0.95) Independent character support Clade Clade Clade Clade Clade Clade Clade Clade Clade Clade Yes Yes Yes Yes No Yes Yes Yes Yes Yes Yes⁄⁄⁄ Yes⁄-⁄ Yes⁄⁄⁄ No No No Yes⁄– Yes⁄⁄Yes⁄– Yes⁄-⁄ Yes Yes Yes =vagrant taxa & clade =vagrant taxa & clade =vagrant taxa & clade =clade IVa & IVd Yes =vagrant taxa & clade =clade IVa & IVd Yes Yes Yes Yes Yes Yes Yes Yes Yes Yes Lobate, placodioid thallus morphology Not identified Not identified Vagrant thallus morphology and usnic acid only Vagrant thallus morphology and usnic acid only Vagrant thallus morphology and usnic acid only Not identified Lecanoric/orsellinic acid are exclusive to this lineage Not identified Not identified I (L. novomexicana) II III IV (R. haydenii) IV (R. h. ssp. arbuscula) IV (R. idahonesis) IVa IVb IVc IVd (21–9–2) (3–0–0) (15–0–0) (0–0–1) (1–0–5) (7–0–0) (2–1–0) (1–0–0) (1–0–0) of morphologically diverse vagrant taxa. However, the inclusion of saxicolous attached taxa within this cluster suggests a recent divergence from saxicolous attached forms or an underlying genetic predisposition to vagrancy in at least some saxicolous lineages. Previous research (Leavitt et al., submitted for publication) indentified multiple independent origins of vagrancy within the lichen genus Xanthoparmelia (Parmeliaceae), but our data suggest that that vagrancy in the R. melanophthalma species-complex is limited to a single closely related lineage, even among morphologically distinct vagrant forms. However, a broader sample of vagrant individuals is essential to adequately addressing this question, particularly R. haydenii subsp. haydenii recently described in China (Zheng et al., 2007). Previous studies have used thin-layer-chromatography (TLC) to characterize lichen secondary metabolic products within Rhizoplaca. In this study HPLC provided a more sensitive approach to determine secondary metabolite diversity within the R. melanophthalma group, as many newly reported compounds here would be masked by other compounds, or likely found at levels undetectable by TLC. While data have supported the taxonomic use of some secondary metabolic characters for delimiting lichen taxa (Schmitt and Lumbsch, 2004; Tehler and Källersjö, 2001), other studies found no correlation between chemotypes and lineages identified using molecular phylogenetic reconstructions (Articus et al., 2002; Buschbom and Mueller, 2006; Nelsen and Gargas, 2009; Velmala et al., 2009). We have identified chemical characters corroborating some lineages identified within the R. melanophthalma group, including: clade IVb containing a combination of orsellinic, lecanoric, and gyrophoric acids; and R. haydenii, R. haydenii, ssp. arbuscula, and R. idahoensis all lack aliphatic acids related to constipatic acid. However, we were unable to identify secondary metabolic characters supporting most identified putative lineages, including the most genetically divergent clades. McCune (1987) suggested three hypotheses to explain chemical diversity in the genus Rhizoplaca: (1) chemotypes are sibling species that cannot or seldom hybridize assuming there are no reproductive barriers, (2) factors favoring polymorphism in chemistry do not differ markedly between regions, or (3) the polymorphism is neutral to natural selection. Although the present study was not designed to explicitly test these hypotheses, our results indicate within the usnic/psoromic acid race multiple lineages co-occur. The usnic/psoromic/lecanoric acid race appears to be a distinct lineage also containing specimens lacking lecanoric acid but expressing the lecanoric acid precursor, orsellinic acid. Additional studies will be needed to fully elucidate the relationship between R. melanophthalma s.l. containing lecanoric or orsellinic acids. Our sampling of the usnic acid chemical race in the R. melanophthalma species-complex was limited to a single IVc IVc IVc IVc saxicolous attached individual (715f) and all vagrant taxa. The saxicolous R. melanophthalma chemical race containing placodiolic acid was not sampled and its relationship to sampled taxa remains in question. Porter (1998) reported a correlation between some secondary metabolites and elevation in R. melanophthalma populations along an altitudinal gradient on Thousand Lakes Mountain, Utah. Besides the strict correlation of lecanoric and orsellinic acid with clade IVb, the present study did not identify any specific correlations between lineages identified from molecular data and expressed secondary metabolites on Thousand Lake Mountain, suggesting that the production of most minor compounds may be environmentally induced. A combination of species diversity in lichen-forming symbionts (alga and fungus) and ecological factors may explain secondary metabolite variation among the Thousand Lake Mountain populations (Brunauer et al., 2007). The current study was generally limited to the Intermountain region of western North America, and robust data from a broader geographic sampling will be essential to understand the general geographic distribution of the candidate species identified in this study. We anticipate that with improved sampling, additional lineages may be identified within the R. melanophthalma species-complex, particularly within L. novomexicana s. l. However, with the exception of L. novomexicana, GenBank accessions were recovered within the candidate species identified from our combined dataset set from samples in western North America in the ITS topology, suggesting our candidate species may have wide-spread distributions. The occurrence of cohesive cosmopolitan lineages found sympatrically with closely related divergent populations poses challenging questions about the processes that yield and maintain cohesive lineages within widespread lichenized ascomycetes. Clade-specific ecological or microhabitat differences considered alone do not appear to offer a plausible explanation of how sympatric diversification may occur in the candidate species. Some lineages exhibit extensive microsympatry (i.e., divergent lineages occurring within a single sampled plot), as well as the production of abundant perennial apothecia (sexual fruiting bodies) without detectable gene flow or hybridization between microsympatric individuals. This pattern suggests that candidate species may have achieved a significant level of reproductive isolation. However, the role of spatio-temporal isolation in lichenized fungal reproduction is relatively unexplored. It has been proposed that competition for symbiotic partners may be a major driver of diversity in mutualistic relationships (Bruns, 1995; O’Brien et al., 2009) and investigating competition for symbionts may provide insights into mechanisms that possibly drive sympatric speciation. Within lichenized fungi, gene trees have often been used to infer species boundaries, and the over-reliance on a single locus has Author's personal copy 600 S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 been problematic in delimiting species because gene duplication, horizontal gene transfer, and deep coalescence may create conflict between the sampled gene tree and the true species tree (de Queiroz and Donoghue, 1990; Maddison, 1997). Furthermore, simulations have shown that concatenation of independent genetic loci can mislead phylogeny estimation (Edwards et al., 2007). Interestingly, major topological differences between our concatenated and coalescent-based phylogeny were recovered (Figs. 2 and 7). Although our study emphasized species delimitation, phylogenetic relationships have crucial significance for understanding the role of biogeography, character evolution, and accurately dating divergences. The coalescent-based species tree assembled from the gene trees using ⁄BEAST provides less support for many of the relationships recovered, and these posteriors could more accurately reflect the level of uncertainty in these nodes (Leaché and Rannala, 2010). Our data suggest that a coalescent approach, explicitly incorporating the independent sorting of gene trees while reconstructing the species trees, is likely crucial in order to recover the best estimate of the R. melanophthalma complex species phylogeny. Improved sampling of the R. melanophthalma species-complex with additional genetic loci will be necessary to assemble an unambiguous species tree, stabilize the taxonomy of this group, and infer phylogeographical patterns. In some cases, rapidly evolving molecular characters may reach fixation in ephemerally isolated demes, with the potential to reticulate with other conspecific lineages at some point in the future (O’Hara, 1993). Additionally, phylogenetic structure can extend below the level of the species, particularly within asexual and haploid genomes (Birky et al., 1989, 1983; Davis, 1996; de Queiroz and Donoghue, 1990) making species limits based on molecular data within lichenized fungi particularly susceptible to excessive subdivision. Here we use a Bayesian modeling approach, incorporating the coalescent process and gene tree uncertainties, to generate the posterior probabilities of speciation (Yang and Rannala, 2010). Although the guide tree has been shown to play a critical role in the outcome of the species delimitation model (Leaché and Fujita, 2010), in this study significant posterior probabilities of speciation were recovered between even the closely related vagrant Rhizoplaca species. Given the fairly strict criteria the BPP assumes to designate species (Yang and Rannala, 2010), we feel that recognizing ten species provides a reasonable perspective on species diversity within the sampled R. melanophthalma species-complex. Our results indicate that in some cases morphological and chemical data provide additional lines of evidence of supporting distinct lineages. However, other candidate species are apparently morphologically and chemically cryptic. In spite of the limitations in delimiting taxa using molecular data, most of the candidate species identified in this study, were not supported by diagnostic morphological or chemical characters, and the effective use of genetic data appears to be an essential approach to appropriately identify natural groups in many fungal lineages (Crespo and Pérez-Ortega, 2009). The authors plan a detailed taxonomic revision for the R. melanophthalma species-complex in the near future, including additional taxonomic and morphological sampling to more fully characterize boundaries between candidate species. Results from this study suggest that robust taxon and molecular data sampling, using appropriate empirical operational criteria to delimit species, may provide an improved perspective on the diversification of lichenized fungi, compared to traditional morphological and chemical characters. However, we are not advocating the use of genetic data to the exclusion of other evidence for delimiting species; because corroboration of species boundaries via independent lines of evidence is important to the establishment of robust hypotheses of species diversity. Analysis of the R. melanophthalma species-complex comprises the larger focus of this study, which is using robust species delimitation in morphologically cryptic and recently diverged lichenized fungi. Rhizoplaca – as traditionally circumscribed – is a small, morphologically diverse lichen genus represented by 9 species (Arup and Grube, 2000; Zhou et al., 2006). This study indicates overall diversity within umbilicate Rhizoplaca species may be vastly underestimated, as multiple previously unrecognized lineages were identified within the R. melanophthalma group. Previous studies have identified well-supported lineages within R. chrysoleuca corresponding to two phenotypic groups (Zhou et al., 2006), and divergent, well-supported clades within the outgroup taxon R. chrysoleuca were also recovered in this study, which suggests an additional unrecognized diversity. Extending the present sampling of the R. melanophthalma species-complex to include a broader geographic context and robust sampling of underrepresented lineages will be critical to improve the understanding of the mechanisms driving speciation in lichenized fungi. Furthermore, an extension of the present sampling to other closely related cosmopolitan Rhizoplaca and Lecanora species-complexes will provide a potential opportunity for developing a comprehensive classification system for other closely related taxa. Additionally, continued investigation of independent characters supporting candidate lineages will be essential for generating robust hypotheses of species boundaries. Acknowledgments We thank Byron Adams (Provo), Eric Green (Provo), Roger Rosentreter (Boise), Imke Schmitt (Minnesota), and Jack Sites (Provo) for valuable discussion and comments on early versions of this manuscript; Christopher Jones and Peter Ririe for laboratory assistance; and LauraDawn Leavitt (Provo) and Gajendra Shrestha (Provo) for invaluable help in preparing figures. We would also like to thank Jack Elix (Canberra) for providing a digital HPLC library and Thorsten Lumbsch (Chicago) for a collection of authentic substances. This study was supported, in part, by funds from the University of Minnesota to Imke Schmitt (St. Paul), Brigham Young University graduate mentoring and graduate research fellowship awards to SDL, and a Walmart Foundation Internship Grant to JDF. The funding sources had no role in study design, data collection and analysis, preparation or decision to publish this manuscript. Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.ympev.2011.03.020. References Articus, K., Mattsson, J.E., Tibell, L., Grube, M., Wedin, M., 2002. Ribosomal DNA and beta-tubulin data do not support the separation of the lichens Usnea florida and U. Subfloridana as distinct species. Mycol. Res. 106, 412–418. Arup, U., Grube, M., 2000. Is Rhizoplaca (Lecanorales, lichenized Ascomycota) a monophyletic genus? Can. J. Botany 78, 318–327. Aslan, A., Budak, G., Karabulut, A., 2004. The amounts Fe, Ba, Sr, K, Ca and Ti in some lichens growing in Erzurum province (Turkey). J. Quant. Spectrosc. Radiat. Transf. 88, 423–431. Avise, J.C., Ball, A.S., 1990. Principles of genealogical concordance in species concepts and biological taxonomy. In: Futuyma, D., Antonovics (Eds.), Oxford Surveys in Evolutionary Biology. Oxford University Press, Oxford. Baum, D.A., Shaw, K.L., 1995. Genealogical perspectives on the species problem. In: Hoch, P.C., Stephenson, A.G. (Eds.), Experimental and Molecular Approaches to Plant Biosystematics. Missouri Botanical Garden, St. Louis, pp. 289–303. Baurain, D., Brinkmann, H., Philippe, H., 2007. Lack of resolution in the animal phylogeny: closely spaced cladogeneses or undetected systematic errors? Mol. Biol. Evol. 24, 6–9. Birky Jr., C.W., Fuerst, P., Maruyama, T., 1989. Organelle gene diversity under migration, mutation, and drift: equilibrium expectations, approach to equilibrium, effects of heteroplasmic cells, and comparison to nuclear genes. Genetics 121, 613–627. Author's personal copy S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 Birky Jr., C.W., Maruyama, T., Fuerst, P., 1983. An approach to population and evolutionary genetic theory for genes in mitochondria and chloroplasts, and some results. Genetics 103, 513–527. Blanco, O., Crespo, A., Divakar, P.K., Esslinger, T.L., Hawksworth, D.L., Thorsten Lumbsch, H., 2004a. Melanelixia and Melanohalea, two new genera segregated from Melanelia (Parmeliaceae) based on molecular and morphological data. Mycol. Res. 108, 873–884. Blanco, O., Crespo, A., Elix, J.A., Hawksworth, D.L., Lumbsch, H.T., 2004b. A molecular phylogeny and a new classification of parmelioid lichens containing Xanthoparmelia-type lichenan (Ascomycota: Lecanorales). Taxon 53, 959–975. Boni, M.F., Posada, D., Feldman, M.W., 2007. An exact nonparametric method for inferring mosaic structure in sequence triplets. Genetics 176, 1035–1047. Brunauer, G., Hager, A., Grube, M., Türk, R., Stocker-Wörgötter, E., 2007. Alterations in secondary metabolism of aposymbiotically grown mycobionts of Xanthoria elegans and cultured resynthesis stages. Plant Physiol. Biochem. 45, 146–151. Bruns, T.D., 1995. Thoughts on the processes that maintain local species diversity of ectomycorrhizal fungi. Plant Soil 170, 63–73. Buschbom, J., Mueller, G.M., 2006. Testing ‘‘species pair’’ hypotheses: evolutionary processes in the lichen-forming species complex Porpidia flavocoerulescens and Porpidia melinodes. Mol. Biol. Evol. 23, 574–586. Cansaran, D., Aras, S., Kandemir, I., Halici, M.G., 2006. Phylogenetic relations of Rhizoplaca Zopf. from anatolia inferred from ITS sequence data. Z. Naturforsch. Sect. C: J. Biosci. 61, 405–412. Carbone, I., Kohn, L.M., 1999. A method for designing primer sets for speciation studies in filamentous Ascomycetes. Mycologia 91, 553–556. Cardinale Jr., M., J.V.d.C., Müller, H., Berg, G., Grube, M., 2008. In situ analysis of the bacterial community associated with the reindeer lichen Cladonia arbuscula reveals predominance of Alphaproteobacteria. FEMS Microb. Ecol. 66, 63–71. Carstens, B.C., Dewey, T.A., 2010. Species delimitation using a combined coalescent and information-theoretic approach: an example from North American Myotis bats. Syst. Biol. 54, 400–414. Clement, M., Posada, D., Crandall, K.A., 2000. TCS: a computer program to estimate gene genealogies. Mol. Ecol. 9, 1657–1659. Crespo, A., Lumbsch, H.T., Mattsson, J.-E., Blanco, O., Divakar, P.K., Articus, K., Wiklund, E., Bawingan, P.A., Wedin, M., 2007. Testing morphology-based hypotheses of phylogenetic relationships in Parmeliaceae (Ascomycota) using three ribosomal markers and the nuclear RPB1 gene. Mol. Phylogenet. Evol. 44, 812–824. Crespo, A., Pérez-Ortega, S., 2009. Cryptic species and species pairs in lichens: a discussion on the relationship between molecular phylogenies and morphological characters. Anales del Jardin Botanico de Madrid 66 (S1), 71–81. Culberson, C.F., 1972. Improved conditions and new data for identification of lichen products by standardized thin-layer chromatographic method. J. Chromatogr. A 72, 113–125. Davis, J.I., 1996. Phylogenetics, molecular variation, and species concepts. Bioscience 46, 502–511. Dayrat, B., 2005. Towards integrative taxonomy. Biol. J. Linn. Soc. 85, 407–415. de Queiroz, K., 1998. The General Lineage Concept of Species, Species Criteria, and the Process of Speciation: A Conceptual Unification and Terminological Recommendations. Oxford University Press, Oxford. de Queiroz, K., 1999. The general lineage concept of species and the defining properties of the species category. Species, New Interdisciplinary Essays, 49–89. de Queiroz, K., 2007. Species concepts and species delimitation. Syst. Biol. 56, 879– 886. de Queiroz, K., Donoghue, M.J., 1990. Phylogenetic systematics and species revisited. Cladistics 6, 83–90. Dettman, J.R., Jacobson, D.J., Taylor, J.W., 2003. A multilocus genealogical approach to phylogenetic species recognition in the model eukaryote Neurospora. Evolution 57, 2703–2720. Dillman, K.L., 1996. Use of the lichen Rhizoplaca melanophthalma as a biomonitor in relation to phosphate refineries near Pocatello. Idaho Environ. Pollut. 92, 91–96. Drummond, A., Ho, S., Phillips, M., Rambaut, A., 2006. Relaxed phylogenetics and dating with confidence. PLoS Biol. 4, e88. Duminil, J., Caron, H., Scotti, I., Cazal, S.-O., Petit, R., 2006. Blind population genetics survey of tropical rainforest trees. Mol. Ecol. 15, 3505–3513. Edgar, R.C., 2004. MUSCLE: a multiple sequence alignment method with reduced time and space complexity. BMC Bioinformatics 5, 1–19. Edwards, S.V., Liu, L., Pearl, D.K., 2007. High-resolution species trees without concatenation. Proc. Natl. Acad. Sci. USA 104, 5936–5941. Evanno, G., Regnaut, S., Goudet, J., 2005. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol. Ecol. 14, 2611–2620. Falush, D., Stephens, M., Pritchard, J.K., 2003. Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics 164, 1567–1587. Feige, G.B., Lumbsch, H.T., Huneck, S., Elix, J.A., 1993. Identification of lichen substances by a standardized high-performance liquid chromatographic method. J. Chromatogr. 646, 417–427. Funk, D.J., Nosil, P., Etges, W.J., 2006. Ecological divergence exhibits consistently positive associations with reproductive isolation across disparate taxa. Proc. Natl. Acad. Sci. USA 103, 3209–3213. Ganley, A.R.D., Kobayashi, T., 2007. Highly efficient concerted evolution in the ribosomal DNA repeats: total rDNA repeat variation revealed by whole-genome shotgun sequence data. Genome Res. 17, 184–191. Gardes, M., Bruns, T.D., 1993. ITS primers with enhanced specificity for basidiomycetes – application to the identification of mycorrhizae and rusts. Mol. Ecol. Notes 2, 113–118. 601 Gibbs, M.J., Armstrong, J.S., Gibbs, A.J., 2000. Sister-scanning: a Monte Carlo procedure for assessing signals in recombinant sequences. Bioinformatics 16, 573–582. Groeneveld, L., Weisrock, D., Rasoloarison, R., Yoder, A., Kappeler, P., 2009. Species delimitation in lemurs: multiple genetic loci reveal low levels of species diversity in the genus Cheirogaleus. BMC Evol. Biol. 9, 30. Grube, M., Cardinale, M., de Castro Jr., J.V., Muller, H., Berg, G., 2009. Species-specific structural and functional diversity of bacterial communities in lichen symbioses. ISME J. 3, 1105–1115. Gutiérrez, G., Blanco, O., Divakar, P., Lumbsch, H., Crespo, A., 2007. Patterns of Group I intron presence in nuclear SSU rDNA of the Lichen Family Parmeliaceae. J. Mol. Evol. 64, 181–195. Hale, M.E., 1990. A Synopsis of the Lichen Genus Xanthoparmelia (Vainio) Hale (Ascomycotina, Parmeliaceae). Smithsonian Institution Press, Washington, DC. Hall, B.G., 2007. Phylogenetic Trees Made Easy: A How-To Manual. Sinauer Associates, Sunderland, MA. Heath, L., van der Walt, E., Varsani, A., Martin, D.P., 2006. Recombination patterns in Aphthoviruses mirror those found in other picornaviruses. J. Virol. 80, 11827– 11832. Heled, J., Drummond, A.J., 2010. Bayesian inference of species trees from multilocus data. Mol. Biol. Evol. 27, 570–580. Hey, J., Waples, R.S., Arnold, M.L., Butlin, R.K., Harrison, R.G., 2003. Understanding and confronting species uncertainty in biology and conservation. Trends Ecol. Evol. 18, 597–603. Hodkinson, B., Lutzoni, F., 2009. A microbiotic survey of lichen-associated bacteria reveals a new lineage from the Rhizobiales. Symbiosis 49, 163–180. Holmes, E.C., Worobey, M., Rambaut, A., 1999. Phylogenetic evidence for recombination in dengue virus. Mol. Biol. Evol. 16, 405–409. Hudson, R.R., Coyne, J.A., 2002. Mathematical consequences of the genealogical species concept. Evolution 56, 1557–1565. Huelsenbeck, J.P., Rannala, B., 2004. Frequentist properties of Bayesian posterior probabilities of phylogenetic trees under simple and complex substitution models. Syst. Biol. 53, 904–913. Huelsenbeck, J.P., Ronquist, F., 2001. MrBayes: bayesian inference of phylogenetic trees. Bioinformatics 17, 754–755. Huelsenbeck, J.P., Ronquist, F., Nielsen, R., Bollback, J.P., 2001. Bayesian inference of phylogeny and its impact on evolutionary biology. Science 294, 2310–2314. Huneck, S., 1999. The significance of lichens and their metabolites. Naturwissenschaften 86, 559–570. Huneck, S., Yoshimura, I., 1996. Identification of Lichen Substances. Springer, New York. Knowles, L.L., Carstens, B.C., 2007. Delimiting species without monophyletic gene trees. Syst. Biol. 56, 887–895. Kosakovsky Pond, S.L., Posada, D., Gravenor, M.B., Woelk, C.H., Frost, S.D.W., 2006. GARD: a genetic algorithm for recombination detection. Bioinformatics 22, 3096–3098. Kroken, S., Taylor, J.W., 2001. A gene genealogical approach to recognize phylogenetic species boundaries in the lichenized fungus Letharia. Mycologia 93, 38–53. LaGreca, S., Lumbsch, H.T., 2001. Three species of Lecanora new to North America, with notes on other poorly known Lecanoroid lichens. Bryologist 104, 204–211. Latch, E., Dharmarajan, G., Glaubitz, J., Rhodes, O., 2006. Relative performance of bayesian clustering software for inferring population substructure and individual assignment at low levels of population differentiation. Conserv. Genet. 7, 295–302. Le Gac, M., Hood, M.E., Fournier, E., Giraud, T., 2007. Phylogenetic evidence of hostspecific cryptic species in the anther smut fungus. Evolution 61, 15–26. Leache, A.D., 2009. Species tree discordance traces to phylogeographic clade boundaries in North American fence lizards (Sceloporus). Syst. Biol. 58, 547–559. Leaché, A.D., Fujita, M.K., 2010. Bayesian species delimitation in West African forest geckos (Hemidactylus fasciatus). Proc. Roy. Soc. B: Biol. Sci. 277, 3071–3077. Leaché, A.D., Rannala, B., 2010. The accuracy of species tree estimation under simulation: a comparison of methods. Syst. Biol. doi:10.1093/sysbio/syq073. Leavitt, S.D., Johnson, L.A., Goward, T., St. Clair, L.L., submitted for publication. Species delimitation in taxonomically difficult lichen-forming fungi: an example from morphologically and chemically diverse Xanthoparmelia (Parmeliaceae) in North America. Leavitt, S.D., Johnson, L.A., St. Clair, L.L., 2011. Species delimitation and evolution in morphologically and chemically diverse communities of the lichen-forming genus Xanthoparmelia (Parmeliaceae, Ascomycota) in western North America. Am. J. Bot. 98, 1–14. Leuckert, C., Poelt, J., Hahnel, G., 1977. Zur Chemotaxonomie der eurasischen Arten der Flechtengattung rhizoplaca. Nova Hedwigia 28, 71–129. Librado, P., Rozas, J., 2009. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452. Lindblom, L., Ekman, S., 2006. Genetic variation and population differentiation in the lichen-forming ascomycete Xanthoria parietina on the island Storfosna, central Norway. Mol. Ecol. 15, 1545–1559. Lumbsch, H.T., 2002. Analysis of phenolic products in lichens. In: Kranner, I., Beckett, R.P., Varma, A. (Eds.), Protocols in Lichenology. Springer, Berlin, pp. 281–295. Lumbsch, H.T., Schmitt, I., 2001. Molecular data suggest that the Lichen Genus Pertusaria is not monophyletic. Lichenologist 33, 161–170. Lutzoni, F., Kauff, F., Cox, C.J., McLaughlin, D., Celio, G., Dentinger, B., Padamsee, M., Hibbett, D., James, T.Y., Baloch, E., Grube, M., Reeb, V., Hofstetter, V., Schoch, C., Arnold, A.E., Miadlikowska, J., Spatafora, J., Johnson, D., Hambleton, S., Crockett, Author's personal copy 602 S.D. Leavitt et al. / Molecular Phylogenetics and Evolution 59 (2011) 587–602 M., Shoemaker, R., Sung, G.-H., Lucking, R., Lumbsch, T., O’Donnell, K., Binder, M., Diederich, P., Ertz, D., Gueidan, C., Hansen, K., Harris, R.C., Hosaka, K., Lim, Y.-W., Matheny, B., Nishida, H., Pfister, D., Rogers, J., Rossman, A., Schmitt, I., Sipman, H., Stone, J., Sugiyama, J., Yahr, R., Vilgalys, R., . Assembling the fungal tree of life: progress, classification, and evolution of subcellular traits. Am. J. Bot. 91, 1446–1480. Lutzoni, F., Pagel, M., Reeb, V., 2001. Major fungal lineages are derived from lichen symbiotic ancestors. Nature 411, 937–940. Maddison, W.P., 1997. Gene trees in species trees. Syst. Biol. 46, 523–536. Martin, D.P., Posada, D., Crandall, K.A., Williamson, C., 2005a. A modified bootscan algorithm for automated identification of recombinant sequences and recombination breakpoints. Anglais 21, 98–102. Martin, D.P., Williamson, C., Posada, D., 2005b. RDP2: recombination detection and analysis from sequence alignments. Bioinformatics 21, 260–262. Mayden, R.L., 1999. Consilience and a hierarchy of species concepts: advances towards closure on the species puzzle. J. Nematol. 31, 95–116. McCune, B., 1987. Distribution of chemotypes of Rhizoplaca in North America. Bryologist 90, 6–14. McCune, B., Rosentreter, R., 2007. Biotic Soil Crust Lichens of the Columbian Basin. Monographs in North American Lichenology, Corvalis. McGuire, G., Wright, F., 2000. Improved detection of mosaic sequences within multiple sequence alignments. Bioinformatics 16, 130–134. Myllys, L., Lohtander, K., Tehler, A., 2001. b-Tubulin, ITS and group I intron sequences challenge the species pair concept in Physcia aipolia and P. caesia. Mycologia 93, 335–343. Nei, M., 1987. Molecular Evolutionary Genetics. Columbia University Press, New York. Nelsen, M.P., Gargas, A., 2009. Assessing clonality and chemotype monophyly in Thamnolia (Icmadophilaceae). Bryologist 112, 42–53. Nylander, J.A.A., Ronquist, F., Huelsenbeck, J.P., Nieves-Aldrey, J., 2004. Bayesian phylogenetic analysis of combined data. Syst. Biol. 53, 47–67. O’Brien, H.E., Miadlikowska, J., Lutzoni, F., 2009. Assessing reproductive isolation in highly diverse communities of the lichen-forming fugnal genus Peltigera. Evolution 63, 2076–2086. O’Hara, R.J., 1993. Systematic generalization, historical fate, and the species problem. Syst. Biol. 42, 231–246. O’Meara, B.C., 2010. New heuristic methods for joint species delimitation and species tree inference. Syst. Biol. 59, 59–73. Ott, S., Brinkmann, M., Wirtz, N., Lumbsch, H.T., 2004. Mitochondrial and nuclear ribosomal DNA data do not support the separation of the Antarctic lichens Umbilicaria kappenii and Umbilicaria antarctica as distinct species. Lichenologist 36, 227–234. Padidam, M., Sawyer, S., Fauquet, C.M., 1999. Possible emergence of new geminiviruses by frequent recombination. Virology 265, 218–225. Porter, L., 1998. Rhiz. Brigham Young University, Provo, Utah. Posada, D., Crandall, K.A., 2001a. Evaluation of methods for detecting recombination from DNA sequences: computer simulations. Proc. Natl. Acad. Sci. USA 98, 13757–13762. Posada, D., Crandall, K.A., 2001b. Selecting the best-fit model of nucleotide substitution. Syst. Biol. 50, 580–601. Pringle, A., Baker, D.M., Platt, J.L., Wares, J.P., Latgé, J.P., Taylor, J.W., 2005. Cryptic speciation in the cosmopolitan and clonal human pathogenic fungus Aspergillus fumigatus. Evolution 59, 1886–1899. Printzen, C., 2009. Lichen systematics: the role of morphological and molecular data to reconstruct phylogenetic relationships. Progress in Botany 71. Springer Berlin Heidelberg, pp. 233–275. Pritchard, J.K., Stephens, M., Donnelly, P., 2000. Inference of population structure using multilocus genotype data. Genetics 155, 945–959. Rambaut, A., Drummond, A.J., 2007. Tracer v1.5. <http://beast.bio.ed.ac.uk/Tracer>. Rannala, B., Yang, Z., 2003. Bayes estimation of species divergence times and ancestral population sizes using DNA sequences from multiple loci. Genetics 164, 1645–1656. Roe, A.D., Sperling, F.A.H., 2007. Population structure and species boundary delimitation of cryptic Dioryctria moths: an integrative approach. Mol. Ecol. 16, 3617–3633. Rokas, A., Williams, B.L., King, N., Carroll, S.B., 2003. Genome-scale approaches to resolving incongruence in molecular phylogenies. Nature 425, 798–804. Rosentreter, R., 1993. Vagrant lichens in North America. Bryologist 96, 333–338. Rowe, G., Beebee, T.J.C., 2007. Defining population boundaries: use of three Bayesian approaches with microsatellite data from British natterjack toads (Bufo calamita). Mol. Ecol. 16, 785–796. Ruiz-Sanchez, E., Sosa, V., 2010. Delimiting species boundaries within the neotropical bamboo Otatea (Poaceae: Bambusoideae) using molecular, morphological and ecological data. Mol. Phylogenet. Evol. 54, 344–356. Ryan, B.D., 2001. Rhizoplaca. In: Nash, T.H., III., Ryan, B.D., Gries, C., Bungartz, F. (Eds.), Lichen Flora of the Greater Sonoran Desert Region. Lichens Unlimited, Tempe. Saisho, D., Purugganan, M.D., 2007. Molecular phylogeography of domesticated barley traces expansion of agriculture in the old world. Genetics 177, 1765–1776. Schluter, D., 2000. The Ecology of Adaptive Radiations. Oxford University Press, Oxford. Schmitt, I., Crespo, A., Divakar, P.K., Fankhauser, J.D., Herman-Sackett, E., Kalb, K., Nelsen, M.P., Rivas-Plata, E., Shimp, A.D., Widhelm, T., Lumbsch, H.T., 2009. New primers for promising single-copy genes in fungal phylogenies and systematics. Persoonia 23, 35–40. Schmitt, I., Lumbsch, H.T., 2004. Molecular phylogeny of the Pertusariaceae supports secondary chemistry as an important systematic character set in lichen-forming ascomycetes. Mol. Phylogenet. Evol. 33, 43–55. Selbmann, L., Zucconi, L., Ruisi, S., Grube, M., Cardinale, M., Onofri, S., 2010. Culturable bacteria associated with Antarctic lichens: affiliation and psychrotolerance. Polar Biol. 33, 71–83. Sites, J.W., Marshall, J.C., 2004. Operational criteria for delimiting species. Ann. Rev. Ecol. Evol. Syst. 35, 199–227. Šlapeta, J., López-García, P., Moreira, D., 2006. Global dispersal and ancient cryptic species in the smallest Marine eukaryotes. Mol. Biol. Evol. 23, 23–29. Smith, J.M., 1992. Analyzing the mosaic structure of genes. J. Mol. Evol. 34, 126–129. Stamatakis, A., 2006. RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics 22, 2688– 2690. Stamatakis, A., Hoover, P., Rougemont, J., 2008. A rapid bootstrap algorithm for the RAxML web servers. Syst. Biol. 57, 758–771. Tehler, A., Källersjö, M., 2001. Parmeliopsis ambigua and P. hyperopta (Parmeliaceae): species or chemotypes? Lichenologist 33, 403–408. Templeton, A.R., Crandall, K.A., Sing, C.F., 1992. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics 132, 619–633. Templeton, A.R., Sing, C.F., 1993. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping. IV. Nested Analyses with cladogram uncertainty and recombination. Genetics 134, 659–669. Thell, A., 1999. Group I intron versus ITS sequences in phylogeny of cetrarioid lichens. Lichenologist 31, 441–449. Ugur, A., Özden, B., Saç, M., Yener, G., Altinbasß, Ü., Kurucu, Y., Bolca, M., 2004. Lichens and mosses for correlation between trace elements and 210Po in the areas near coal-fired power plant at Yatağan, Turkey. J. Radioanal. Nucl. Chem. 259, 87–92. Velmala, S., Myllys, L., Halonen, P., Goward, T., Ahti, T., 2009. Molecular data show that Bryoria fremontii and B. Tortuosa (Parmeliaceae) are conspecific. Lichenologist 41, 231–242. Vieites, D.R., Wollenberg, K.C., Andreone, F., Köhler, J., Glaw, F., Vences, M., 2009. Vast underestimation of Madagascar’s biodiversity evidenced by an integrative amphibian inventory. Proc. Natl. Acad. Sci. 106, 8267–8272. Wedin, M., Westberg, M., Crewe, A.T., Tehler, A., Purvis, O.W., 2009. Species delimitation and evolution of metal bioaccumulation in the lichenized Acarospora smaragdula (Ascomycota, Fungi) complex. Cladistics 25, 161–172. Weisrock, D.W., Rasoloarison, R.M., Fiorentino, I., Ralison, J.M., Goodman, S.M., Kappeler, P.M., Yoder, A.D., 2010. Delimiting species without nuclear monophyly in Madagascar’s mouse lemurs. PLoS ONE 5, e9883. Wheeler, D.L., Barrett, T., Benson, D.A., Bryant, S.H., Canese, K., Chetvernin, V., Church, D.M., DiCuccio, M., Edgar, R., Federhen, S., Geer, L.Y., Kapustin, Y., Khovayko, O., Landsman, D., Lipman, D.J., Madden, T.L., Maglott, D.R., Ostell, J., Miller, V., Pruitt, K.D., Schuler, G.D., Sequeira, E., Sherry, S.T., Sirotkin, K., Souvorov, A., Starchenko, G., Tatusov, R.L., Tatusova, T.A., Wagner, L., Yaschenko, E., 2006. Database resources of the National Center for Biotechnology Information. Nucl. Acids Res., gkl1031. White, T.J., Bruns, T.D., Lee, S., Taylor, J. (Eds.), 1990. Amplification and Direct Sequencing of Fungal Ribosomal RNA Genes for Phylogenetics. Academic Press, San Diego. Wiens, J., Penkrot, T., 2002. Delimiting species using DNA and morphological variation and discordant species limits in spiny lizards (Sceloporus). Syst. Biol. 51, 69–91. Wiens, J.J., 1998. Combining data sets with different phylogenetic histories. Syst. Biol. 47, 568–581. Wirtz, N., Printzen, C., Lumbsch, H.T., 2008. The delimitation of Antarctic and bipolar species of neuropogonoid Usnea (Ascomycota, Lecanorales): a cohesion approach of species recognition for the Usnea perpusilla complex. Mycol. Res. 112, 472–484. Yang, Z., Rannala, B., 2010. Bayesian species delimitation using multilocus sequence data. Proc. Natl. Acad. Sci. 107, 9264–9269. Zheng, X.-L., Sheng, H.-M., An, L.-Z., 2007. Phylogenetic analysis of lichen-forming fungi Rhizoplaca Zopf from China based on ITS data and morphology. Z. Naturforsch. Sect. C: J. Biosci. 62c, 757–764. Zhou, Q.-M., Guo, S.-Y., Huang, M.-R., Wei, J.-C., 2006. A study of the genetic variability of Rhizoplaca chrysoleuca using DNA sequences and secondary metabolic substances. Mycologia 98, 57–67.