C G Boluda1,2, V J Rico1, P K Divakar1, O Nadyeina2, L Myllys3, R T McMullin4, J C Zamora1,5, C Scheidegger2, D L Hawksworth6. 1. Departamento de Farmacología, Farmacognosia y Botánica (U.D. Botánica), Facultad de Farmacia, Universidad Complutense, Plaza de Ramón y Cajal s/n, Madrid 28040, Spain. 2. Biodiversity and Conservation Biology, Swiss Federal Research Institute WSL, Zürcherstrasse 111, Birmensdorf 8903, Switzerland. 3. Botanical Museum, Finnish Museum of Natural History, P.O. Box 7, 00014 University of Helsinki, Finland. 4. Research and Collections, Canadian Museum of Nature, Ottawa, ON K1P 6P4, Canada. 5. Museum of Evolution, Uppsala University, Norbyvägen 16, 75236 Uppsala, Sweden. 6. Department of Life Sciences, The Natural History Museum, Cromwell Road, London SW7 5BD, UK; and Comparative Plant and Fungal Biology, Royal Botanic Gardens, Kew, Surrey TW9 3DS, United Kingdom.
Abstract
In many lichen-forming fungi, molecular phylogenetic analyses lead to the discovery of cryptic species within traditional morphospecies. However, in some cases, molecular sequence data also questions the separation of phenotypically characterised species. Here we apply an integrative taxonomy approach - including morphological, chemical, molecular, and distributional characters - to re-assess species boundaries in a traditionally speciose group of hair lichens, Bryoria sect. Implexae. We sampled multilocus sequence and microsatellite data from 142 specimens from a broad intercontinental distribution. Molecular data included DNA sequences of the standard fungal markers ITS, IGS, GAPDH, two newly tested loci (FRBi15 and FRBi16), and SSR frequencies from 18 microsatellite markers. Datasets were analysed with Bayesian and maximum likelihood phylogenetic reconstruction, phenogram reconstruction, STRUCTURE Bayesian clustering, principal coordinate analysis, haplotype network, and several different species delimitation analyses (ABGD, PTP, GMYC, and DISSECT). Additionally, past population demography and divergence times are estimated. The different approaches to species recognition do not support the monophyly of the 11 currently accepted morphospecies, and rather suggest the reduction of these to four phylogenetic species. Moreover, three of these are relatively recent in origin and cryptic, including phenotypically and chemically variable specimens. Issues regarding the integration of an evolutionary perspective into taxonomic conclusions in species complexes, which have undergone recent diversification, are discussed. The four accepted species, all epitypified by sequenced material, are Bryoria fuscescens, B. glabra, B. kockiana, and B. pseudofuscescens. Ten species rank names are reduced to synonymy. In the absence of molecular data, they can be recorded as the B. fuscescens complex. Intraspecific phenotype plasticity and factors affecting the speciation of different morphospecies in this group of Bryoria are outlined.
In many lichen-forming fungi, molecular phylogenetic analyses lead to the discovery of cryptic species within traditional morphospecies. However, in some cases, molecular sequence data also questions the separation of phenotypically characterised species. Here we apply an integrative taxonomy approach - including morphological, chemical, molecular, and distributional characters - to re-assess species boundaries in a traditionally speciose group of hair lichens, Bryoria sect. Implexae. We sampled multilocus sequence and microsatellite data from 142 specimens from a broad intercontinental distribution. Molecular data included DNA sequences of the standard fungal markers ITS, IGS, GAPDH, two newly tested loci (FRBi15 and FRBi16), and SSR frequencies from 18 microsatellite markers. Datasets were analysed with Bayesian and maximum likelihood phylogenetic reconstruction, phenogram reconstruction, STRUCTURE Bayesian clustering, principal coordinate analysis, haplotype network, and several different species delimitation analyses (ABGD, PTP, GMYC, and DISSECT). Additionally, past population demography and divergence times are estimated. The different approaches to species recognition do not support the monophyly of the 11 currently accepted morphospecies, and rather suggest the reduction of these to four phylogenetic species. Moreover, three of these are relatively recent in origin and cryptic, including phenotypically and chemically variable specimens. Issues regarding the integration of an evolutionary perspective into taxonomic conclusions in species complexes, which have undergone recent diversification, are discussed. The four accepted species, all epitypified by sequenced material, are Bryoria fuscescens, B. glabra, B. kockiana, and B. pseudofuscescens. Ten species rank names are reduced to synonymy. In the absence of molecular data, they can be recorded as the B. fuscescens complex. Intraspecific phenotype plasticity and factors affecting the speciation of different morphospecies in this group of Bryoria are outlined.
Accurate identification and characterization of species is the basis of communication, conservation, resources management, and material used in biological research. However, in groups of relatively recent origin, species delimitation is often difficult (Jakob & Blattner 2006, Leavitt et al. 2011, Lumley & Sperling 2011). Organisms are always evolving, changing in response to either selective pressures or genetic drift, so that delimiting units to accord species names is not always clear (Naciri & Linder 2015). Several phenomena can hinder species delimitation: phylogenetic/phenotypic mismatches (Articus et al. 2002, Mark et al. 2016, Pino-Bodas et al. 2016), ‘intermediate’ specimens between generally accepted taxa (Seymour et al. 2007), hybridization (Konrad et al. 2002, Steinová et al. 2013), an absence of delimited clades (Jakob & Blattner 2006, Lumley & Sperling 2011), or incomplete lineage sorting (Saag et al. 2014, Leavitt et al. 2016). Long-term reproductive isolation may produce structured, non-overlapping lineages, whereas an intraspecific phylogeny, as well as a recent or contemporary speciation event, may produce reticulated lineages (Abbott et al. 2016).The family Parmeliaceae is one of the most studied amongst lichenised fungi. It contains many genera with species delimitation problems, such as Cetraria aculeata (Lutsak et al. 2017); Letharia (Altermann et al. 2014), the Parmotrema reticulatum complex (Del-Prado et al. 2016), and Pseudephebe (Boluda et al. 2016). In some cases, a lack of correlation between genotypes and phenotypes has led to the recognition of cryptic species within morphologically indistinguishable or scarcely indistinguishable morphospecies (Molina et al. 2011a, b, Leavitt et al. 2012a, b, Singh et al. 2015, Boluda et al. 2016, Del-Prado et al. 2016), and so far, more than 80 cryptic lineages have been detected in Parmeliaceae (Crespo & Lumbsch 2010, Divakar et al. 2010). However, in other cases there is a mismatch between lineages revealed by standard DNA-barcoding markers and long-accepted morphospecies (Articus et al. 2002, Seymour et al. 2007, Velmala et al. 2014, Mark et al. 2016, Kirika et al. 2016a, b, McMullin et al. 2016).In the morphologically similar ‘beard’ and ‘hair’ lichens of the Alectoria sarmentosa, Bryoria sect. Implexae, and Usnea barbata species complexes (Velmala et al. 2014, Mark et al. 2016, McMullin et al. 2016), DNA sequences from standard barcoding markers show that what were considered well delimited morphospecies are found admixed in a single lineage that may be interpreted as a single phylogenetic species. In such situations, many processes may be operative, including environmental plasticity (Boluda et al. 2016), hybridisation, ancestral polymorphisms, incomplete lineage sorting (Joly et al. 2009), limited value of neutral markers (Bekessy et al. 2003), or morphological variability mediated by low selective pressure, genetic drift, or huge population sizes (Hartl & Clark 2007). In these cases, the use of additional markers, especially highly variable ones (e.g., microsatellites), may contribute to an explanation of the underlying phenomena.Chemical characters, mainly the production of polyketides, were accorded major importance in species delimitation in lichen-forming fungi in the 1960s and 1970s (Hawksworth 1976, Lumbsch 1988). These compounds are formed by the fungal partner, and that expression can differ according to the position in a thallus or in pure culture. For almost 50 years, chemical products, generally linked to minor morphological differences, have been used to circumscribe species in Bryoria (Hawksworth 1972, Brodo & Hawksworth 1977, Myllys et al. 2011, Velmala et al. 2014). The advent of molecular phylogenetics has enabled such species concepts to be tested, and they have proved particularly wanting in one group of species, those placed in Bryoria sect. Implexae (Myllys et al. 2011, Velmala et al. 2014, Boluda et al. 2015). Velmala et al. (2014) provided DNA sequence data for 11 species in the section, and with the exception of B. glabra, all the other species were intermixed in clades with diverse, and not concordant, chemical and morphological features. Genetically indistinguishable taxa (with the markers used), maintain distinctive phenotypes even when growing in physical contact with one another (Velmala et al. 2014, Boluda et al. 2015), so the variation cannot be attributed solely to ecological factors.A study on the morphospecies B. fuscescens in central Spain (Boluda et al. 2015) revealed specimens with the same nuclear internal transcribed spacer sequence (nuITS) but different extrolites (compounds formed on the surface of or excreted from hyphae). Subsequent fieldwork across Europe has revealed further combinations of extrolites, and also specimens sharing characters of additional morphospecies. In order to understand the evolutionary processes involved in B. fuscescens and related species we have adopted an integrative approach including morphological, distributional, and chemical data together with DNA sequences from three standard loci (Schoch et al. 2012), two newly tested loci, and eighteen microsatellite (SSRs) markers (Nadyeina et al. 2014). We then analysed these datasets in a rigorous statistical framework to effectively integrate an evolutionary perspective into a revised and defensible taxonomic treatment. These studies are reported here, and we anticipate that the experience gained in this group of lichens will inform how other species complexes with similarly discordant datasets can be addressed.
MATERIALS AND METHODS
Sampling
We examined 142 specimens from 14 countries in Europe, the Mediterranean Basin, and North and South America, representing 11 named morphospecies in Bryoria sect. Implexae (Table 1). Our dataset included 91 of the 97 specimens used by Velmala et al. (2014) in their revision of B. sect. Implexae. Newly obtained sequences are shown in bold in Table 1. Bryoria furcellata was used as outgroup to root the tree (Velmala et al. 2014). Names used in the analyses follow the species concepts adopted in Velmala et al. (2014).
Specimen information and GenBank accession numbers of the Bryoria sect. Implexae samples used in this study. Newly obtained sequences are in bold.
The newly studied specimens (Table 1) were examined morphologically under a Nikon SMZ-1000 dissecting microscope, and hand-cut sections studied with a Nikon Eclipse-80i compound microscope equipped with bright field and differential interference contrast (DIC). Habit photographs were taken with a Nikon 105 mm f/2.8D AF Micro-Nikkor Lens coupled to a Nikon D90 camera with daylight. Spot tests (K, C, and PD) and TLC were carried out following Orange et al. (2010). Solvent system C (200 ml toluene / 30 mL acetic acid) was used for TLC, with concentrated acetone extracts at 50 °C spotted onto silica gel 60 F254 aluminium sheets (Merck, Darmstadt, Germany). Spotted sheets were dried for 10 min in an acetic acid atmosphere to maximize resolution. Segments from the same lichen branch were used for both TLC and DNA extraction to avoid the possible risk of taking samples from mixed collections. Morphological and thin layer chromatographic (TLC) analyses of the samples used in Velmala et al. (2014; Table 1) were taken from that study.
DNA dataset
The molecular dataset comprised DNA sequences and SSRs frequencies. DNA extraction was performed with the DNeasy Plant Mini Kit (Qiagen, Barcelona, Spain), following the manufacturer’s instructions.Eighteen fungal-specific microsatellites markers (Bi01, Bi02, Bi03, Bi04, Bi05, Bi06, Bi07, Bi08, Bi09, Bi10, Bi11, Bi12, Bi13, Bi14, Bi15, Bi16, Bi18 and Bi19) were amplified following Nadyeina et al. (2014) using fluorescently labelled primers. Fragment lengths were determined on an ABI PRISM® 3130 Genetic Analyser (Life Technologies, Carlsbad, CA, USA). Genotyping was performed using GeneScan-500 LIZ as the internal size standard and GeneMapper v. 3.7 (Applied Biosystems, Foster City, CA, USA).For DNA sequencing, five loci were selected (Table 2), three commonly used as standard markers in fungi (ITS, IGS, and GAPDH), which were also used in Velmala et al. (2014), and two microsatellite flanking regions tested here for the first time (FRBi15 and FRBi16). Microsatellite flanking regions are variable non-coding DNA fragments that can contain phylogenetical signal through a neutral molecular evolution (Zardoya et al. 1996, Chatrou et al. 2009). To explore this possibility, the flanking regions of the 18 microsatellite markers were checked upstream and downstream in the 454 pyrosequencing contigs used for microsatellite searching in Nadyeina et al. (2014). The variability of each region was assessed with the number of variable sites in contigs supported by 2–16 copies. From the 36 regions (two for each of the 18 microsatellites), the most variable flanking regions were in Bi15 and Bi16, and specific primers were designed for those loci (Table 2).
Table 2
Primer information used in Bryoria sect. Implexae.
Marker
Description
Primer forward (5’–3’)
Source
Primer reverse (5’–3’)
Source
ITS
Internal transcribed spacers of the nuclear rDNA including the 5.8S region
Flanking region of Bryoria sect.Implexae microsatellite marker 15
FRBi15f:
This paper
FRBi15r:
This paper
GTCATAAGGGTATCAATCC
TGAAAAGGTTTGGTGACTC
FRBI16
Flanking region of Bryoria sect. Implexae microsatellite marker 16
FRBi16f:
This paper
FRBi16r:
This paper
CGAGGTTTCAGGAAAGGGAA
AGGAAGTGATGTCGAGGT
New DNA sequences (Table 1) were obtained using polymerase chain reactions (PCRs) as follows: a reaction mixture of 25 μL, containing 12 μL sterile water, 9 μL JumpStartTM REDTaq ReadyMix PCR Reaction Mix (Sigma-Aldrich, St Louis, MI, USA), 1.25 μL of each primer (forward and reverse) at 10 μM, and a 1.5 μL DNA template. Cycling conditions for ITS, GAPDH, FRBi15, and FRBi16 were 2 min at 94 °C; 35 cycles of 30 s at 94 °C; 30 s at 56 °C; 2 min at 72 °C; and a final extension of 5 min at 72 °C. For IGS, the cycling process was: 2 min at 94 °C; 15 cycles of 30 s at 94 °C, 30 s at 55 °C (decreasing 1 °C each cycle down to 40 °C), 2 min at 72 °C, then 35 cycles of 30 s at 94 °C, 30 s at 55 °C; 90 s at 72 °C, and a final extension of 5 min at 72 °C. PCR products were checked and quantified on 1 % agarose gel stained with ethidium bromide and cleaned using Exonuclease I and FastAP Thermosensitive Alkaline Phosphatase (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s instructions. Sequencing was performed with labelling using BigDye Terminator v. 3.1 Kit (Applied Biosystems) as follows: 25 cycles of 20 s at 96 °C, 5 s at 50 °C, and 2 min at 60 °C. PCR products were cleaned-up with the BigDye XTerminator Purification Kit (Applied Biosystems) according to the manufacturer’s instructions. Sequences were obtained in an ABI PRISM 3130 Genetic Analyser (Life Technologies) and manually adjusted using DNA Workbench v. 6 (CLC bio, Aarhus, Denmark) and MEGA5 (Tamura et al. 2011). Newly generated sequences were deposited in GenBank (Table 1 in bold).
Clustering methodologies
Phenetic analyses
Two presence/absence (1/0) matrices were constructed, one for the extrolites detected by TLC, and another with morphology and geography data (Appendix 1). Morphological characters scored comprised those traditionally used to separate morphospecies in the group:pale/dark thallus colour;branching angles (acute/obtuse/mixed);soralia (absent/fissural/tuberculate/both); andpseudocyphellae (conspicuous/inconspicuous).For distributions, Old World vs New World was used. The R package cluster (Maechler et al. 2013) was used to obtain the dissimilarity matrix, and then the pvclust package (Suzuki & Shimodaira 2006) was run to obtain a phenogram (Zamora et al. 2013). Multiscale bootstrap resampling with 10 000 bootstrap (bp) replicates was used to obtain approximately unbiased (au) p-values for branch supports. Groups were considered as supported when bp values exceeded 70 or au values exceeded 95.
Phylogenetic tree
Alignments for each locus were performed using MAFFT v. 7 (http://mafft.cbrc.jp/alignment/server/; Katoh & Standley 2013) with the G-INS-i alignment algorithm, a ‘1PAM/K = 2’ scoring matrix, with an offset value of 0.1, and the remaining parameters set as default. Alignments were deposited in TreeBASE under accession nos TB2:S20007 (ITS, IGS, and GAPDH), TB2:S20005 (FRBi15), and TB2:S20004 (FRBi16). RDP v. 4 (Martin et al. 2010) was used to detect potential recombination events, through the methods RDP (Martin & Rybicki 2000), GENECONV (Padidam et al. 1999), Chimaera (Posada & Crandall 2001), Maxchi (Maynard-Smith 1992), Bootscan (Gibbs et al. 2000, Martin et al. 2005), SiScan (Weiller 1998, Gibbs et al. 2000), PhylPro (Weiller 1998), and 3Seq (Boni et al. 2007). Partitionfinder (Lanfear et al. 2012) was used to detect possible intra-locus substitution model variability, resulting in the splitting of the ITS region into ITS1, 5.8S, and ITS2, and coding each codon position separately in GAPDH. Models of DNA sequence evolution for each locus partition were selected with jModeltest v. 2.0 (Darriba et al. 2012), using the Akaike information criterion (AIC, Akaike 1974). The best-fit model of evolution obtained was: ITS1 = TIM2, 5.8S = K80, ITS2 = TIM2ef + G, IGS = TrN + I, GAPDH 1st position = TrN + I, GAPDH 2nd position = F81 + I, GAPDH 3th position = TPM3uf, FRBi15 = TPM3uf + I, FRBi16 = TPM3uf + G. To detect possible topological conflicts among loci, the CADM test (Legendre & Lapointe 2004, Campbell et al. 2011) was performed using the function ‘CADM.global’ implemented in the library ‘ape’ of R (Paradis et al. 2004). As loci FRBi15 and FRBi16 were not congruent among them and neither with the remaining loci, three alignments were used, resulting in three trees, one for each FRBi region and another for the concatenated dataset including loci ITS, IGS and GAPDH. For the concatenated matrix, specimens with more than one missing locus were excluded. Datasets were analysed using maximum likelihood (ML) and Bayesian (B/MCMCMC) approaches with gaps treated as missing data.For ML tree reconstruction, we used RAxML v. 8.2.10 (Stamatakis 2006) implemented in CIPRES Science Gateway (https://www.phylo.org/; Miller et al. 2010) with the GTRGAMMA model (Stamatakis 2006, 2014, Stamatakis et al. 2008). Support values were assessed using the ‘rapid bootstrapping’ option with 1 000 replicates. For the Bayesian reconstruction, MrBayes v. 3.2.1 (Ronquist & Huelsenbeck 2003) was used. Two simultaneous runs with 10 M generations each, starting with a random tree and employing 12 simultaneous chains, were executed. Every 500th tree was saved to a file. Preliminary analysis resulted in an overestimation of branch lengths and to correct this we used the uniform compound Dirichlet prior brlenspr = unconstrained : gammadir (1, 1, 1, 1; Zamora et al. 2015). We plotted the log-likelihood scores of sample points against generations using Tracer v. 1.5 (Rambaut et al. 2014) and determined that stationarity had been achieved when the log-likelihood values of the sample points reached an equilibrium and ESS values exceeded 200 (Huelsenbeck & Ronquist 2001). Posterior probabilities (PPs) were obtained from the 50 % majority rule consensus of sampled trees after excluding the initial 25 % as burn-in. The phylogenetic tree was drawn with FigTree v. 1.4 (Rambaut 2009).
STRUCTURE
STRUCTURE v. 2.3.4 (Pritchard et al. 2000, Falush et al. 2003) was run with the SSRs data matrix. Analysis was computed with 100 000 burn-in generations and 100 000 iterations using a K value from 1 to 12 (i.e., the putative number of species we may have) and 20 replicates for each K. To combine the 20 runs of each K in a single result, CLUMMP v. 1.1.2 (Jakobsson & Rosenberg 2007) was used and visualised replacing the CLUMMP output values in a STRUCTURE output of the same K, and then plotted using the STRUCTURE software. To show the probability of each K value, STRUCTURE HARVESTER (Earl & Von Holdt 2012), with the ΔK method (Evanno et al. 2005) was used, considering the most probable K the first one that appears close to 0 in the output graphic.
Principal coordinate analysis
Principal coordinate analysis (PCoA) was carried out with the SSRs length data in GenAlEx 6.5. The results of the three first axes were plotted in a three-axis graph using The Excel 3D Scatter Plot v. 2.1, in which the graphic can be moved in 3D to obtain a better understanding of how the plots are distributed in the space. Since the projection of this 3D graph on a paper is necessarily confusing, PCoA results were plotted on two different 2D graphs showing axes 1 and 2, and 1 and 3, respectively.
Haplotype network
Haplotype network reconstruction was performed using TCS v. 1.2.1 (Clement et al. 2000) with the concatenated sequences matrix, excluding the outgroup, using gaps as missing data, and a 95 % connection limit. Specimens differing only by missing or ambiguous characters were not counted as haplotypes.
Species delimitation analyses
In order to examine species delimitation, four computational approaches not requiring prior hypothesis of a putative number of species were used:Automatic Barcode Gap Discovery ABGD (Puillandre et al. 2011) based on barcode gaps using genetic distances;Poisson Tree Processes PTP (Zhang et al. 2013), based on gene trees;The Generalized Mixed Yule coalescent approach GMYC, which combines a coalescent model of intraspecific branching with a Yule model for interspecific branching (Pons et al. 2006, Monaghan et al. 2009); andDISSECT (Jones et al. 2014) based on the multispecies coalescent model for species delimitation.ABGD and PTP were carried out using the online servers http://wwwabi.snv.jussieu.fr/public/abgd/ and http://species.h-its.org, respectively. GMYC was analysed with the gmyc function in the SPLITS package in R (v. 2.10, www.cran.r-project.org), employing the single (GMYCs) and multiple (GMYCm) threshold methods. Because GMYC needs a strictly ultrametric and bifurcating tree with no zero branch lengths, identical sequences were deleted and an ultrametric tree was generated using BEAST v. 1.8.2 software (Drumond et al. 2012), with the evolutionary models explained in the Bayesian phylogenetic reconstruction. A run of 100 M iterations logging every 1 000th iteration was conducted. Consensus tree was generated with TreeAnotator v. 1.8.2 after discarding the initial 10 % trees as burn-in. ESS values above 200 were ensured using Tracer v. 1.6 (Rambaut et al. 2014).DISSECT analysis was implemented in STARBEAST (*BEAST, Drumond et al. 2012) using the concatenated DNA matrix after removing identical sequences and following the instructions of Jones et al. (2014). First, we used BEAUti (Drumond et al. 2012) to produce the xml file, with every individual encoded as if it was a separate species. Sites, clocks and trees were released as unlinked. Nucleotide substitution models and other parameters (as in the Bayesian analysis, see above), were encoded using BEAUti if possible, or manually entered. For the ITS locus, a substitution rate of 0.0033 substitutions per site per million years was introduced (Leavitt et al. 2012a, b), setting other loci as estimated with a lognormal relaxed clock. A birth-death-collapse prior that controlled the minimal split heights for the putative resulting species was manually added to the xml file. This prior contained the parameters CollapseHeight (ɛ) with a value of 0.0001 and CollapseWeight (ω), set as estimated using a Beta distribution with values 10 and 1.5. Selected parameters provide a highest probability density around 4–5 clusters, the most probable number of taxa meriting separation according to other analyses performed for this paper. However, this prior is diffuse and allows to obtain a different number of putative taxa if they adjust better to the data. The xml file was executed in BEAST with 250 M MCMC iterations, sampling every 10 000th iteration. Tracer v. 1.6 (Rambaut et al. 2014) was used to assess ESS values above 140. The resulting *BEAST species tree output was then treated with SpeciesDelimitationAnalyzer (Jones et al. 2014), with a burn-in of 5 000 trees (20 % of the total generated), a collapse height of 0.001 (one fraction lower than in the *BEAST analysis) and a simcutoff value of 1 to ignore this parameter, as according to sequence variability, we expected very similar putative species to emerge. The resulting similarity matrix was plotted with R v. 2.15.1 (R Core Team 2014) following the method of Jones et al. (2014).
Divergence time estimation
Two divergence time estimations were performed, one only with the ITS region and a defined substitution rate, and the other with the concatenated data matrix of ITS, IGS, and GAPDH loci. A rate of 3.30 × 10–9 s·s–1·yr–1 for the ITS region as a whole was used, with a GTR + G + I substitution model (Leavitt et al. 2012a, b). In the concatenated matrix analysis, as no previous literature on substitution rates for IGS and GAPDH in lichen-forming fungi is available, these were set as estimated in the Bayesian phylogenetic analysis. A *BEAST analysis was executed, using a relaxed clock model (uncorrelated lognormal), a birth-death model prior for the node heights and unlinked substitution models, clocks and trees for each partition. Clades G, Ko, NA, and WD were selected as potential species, forcing them to remain monophyletic (Fig. 6). No calibration points could be used, as no fossils or previous dating of this species complex are available. To avoid stochastic events, two independent analyses were run, each with 200 million generations, sampling each 5 000 trees, and discarding the first 10 000 trees (25 %) as burn-in. Tracer v. 1.6 (Rambaut et al. 2014) was used to ensure ESS parameter values above 115 in the concatenated matrix, and 185 for the ITS analysis. Different priors were tested but no higher ESS values could be obtained, which we suspect was due to the very similar sequences, and the uncertain topology of the backbone connecting the groups Ko, NA, and WD. The two runs performed for each input were merged with logcombiner v. 1.8.2, and the resulting trees merged in a consensus tree using TreeAnnotator v. 1.8.2 (Drumond et al. 2012). FigTree v. 1.4 (Rambaut 2009) was used to display the ITS and the concatenated dated species trees.
Fig. 6
Integrated assessment of results in Bryoria sect. Implexae. Tree topology depicts the result of the Bayesian Markov chain Monte Carlo (B/MCMC) analysis. Posterior probabilities and bootstrap analysis for the supported nodes (≥ 0.95 and ≥ 70 %) are indicated at the main nodes. For each specimen, the extrolites detected, and the putative number of species predicted by the different methodologies is shown. The left top corner tree shows the results of the molecular dating analysis. — a. Bryoria implexa morphotype (Spain, Asturias, 2013, Boluda, MAF-Lich. 20749); b. B. capillaris morphotype (Spain, Navarra, 2013, Boluda & Villagra, MAF-Lich. 20748); c. B. fuscescens morphotype (Morocco, Ifrane, 2013, Boluda, MAF-Lich. 20751). — Ale. = Alectorialic acid; Bar. = Barbatolic acid; Fum. = Fumarprotocetraric acid; G = Glabra clade; Gyr. = Gyrophoric acid; Hap. Net. = Haplotype Network; HPD = Highest Posterior Density; Ko = Kockiana clade; Mya = Million years ago; NA = North American clade; Nor. = Norstictic acid; PCoA = Principal Coordinates Analysis; Pso. = Psoromic acid; STRUCT. = Structure; WD = Wide Distributed clade.
Demography
Changes in population sizes through time were estimated using the Bayesian skyline analysis (Drumond et al. 2005) with BEAST. Only clades Ko, NA, and WD, isolated and merged, were studied, as they show a clock-like tree topology and adequate sampling sizes.Following the methods used for divergence time estimation analysis, the demography analyses were run using the ITS region without partitioning, with the GTR + G + I model of nucleotide substitution and a substitution rate of 3.30 × 10–9 s·s–1·yr–1, and with a strict molecular clock model (Leavitt et al. 2012a, b). Additionally, the same analysis was repeated with the concatenated data matrix using the ITS substitution rate, estimating the other loci rates with a relaxed clock model and using the nucleotide substitution models for IGS and GAPDH explained in the Bayesian phylogenetic reconstruction. Four independent runs for each input were processed with 50 M MCMC generations, sampling parameter values every 5 000th generation, using the Bayesian Skyline tree prior model, six discreet changes in population size and the linear growth option. ESS values were checked with Tracer v. 1.5 (Rambaut et al. 2014), and the two best of the four runs were combined, obtaining values usually above 200, with some exceptions with a lower limit of 100. Skyline plots were drawn with Tracer v. 1.5.To support the Bayesian skyline test, a neutrality test was performed to infer if populations are in mutation-drift equilibrium. Tajima’s D (Tajima 1989) and Fu’s Fs (Fu 1997) were calculated with DnaSP v. 5.10 (Librado & Rozas 2009). A significantly positive D is interpreted as a diversifying selection or a recent bottleneck, whereas a negative significant D shows purifying selection or a recent expansion. If D is not significantly different from 0, a mutation-drift equilibrium may be occurring. Fu’s Fs can be interpreted in the same way.
RESULTS
Morphological and chemical clustering
The wide geographical range of collections revealed a combination of characters not previously reported in the Bryoria fuscescens complex, especially those from the previously less-studied Mediterranean Basin. Specimens with intermediate morphologies amongst traditionally accepted species were recognized, and the application of species names according to the current taxonomy was ambiguous. Individuals connecting the phenotypes and chemotypes of the taxa currently recognized as Bryoria fuscescens, B. implexa, B. kuemmerleana, and B. vrangiana were not rare in the Mediterranean Basin. For example, chemotypes thought to be diagnostic for a particular taxon were detected in specimens morphologically belonging to other taxa, as well as specimens producing extrolites characteristic of different taxa in a single thallus. Thin-layer chromatography (TLC) revealed seven different extrolites: alectorialic, barbatolic, fumarprotocetraric, gyrophoric, norstictic, and psoromic acids, atranorin, and sometimes also related substances such as chloroatranorin, protocetraric, or connorstictic acids. Atranorin, a typical accessory substance in the genus Bryoria, was not used in the posterior analyses because it appears in trace amounts in many samples and is often difficult to unequivocally discern if it is present or absent by TLC alone.The chemical presence/absence matrix resulted in the phenogram shown in Appendix 2a. The matrix included specimens with as many as four extrolites, something not previously reported in the complex. Chemical characters were separated into two main groups:specimens that contain benzyldepsides (i.e., alectorialic and barbatolic acids), substances traditionally used to separate B. capillaris and B. pikei from other species in the complex; andspecimens without benzyldepsides.The latter were clustered in two well-supported groups, one with fumarprotocetraric acid as the main substance, and the other without it (including specimens with no detectable substances). If the structural relationships of the compounds were encoded in the presence/absence matrix (benzyldepsides vs depsidones), the same clustering was obtained.The analysis of combined morphological, geographical, and chemical characters resulted in the phenogram in Appendix 2b. Only terminal branches were supported, including few monophyletic morphospecies, although not well isolated from others. Neither accepted morphospecies nor an unequivocal number of phenotypic groups could be recognized. This ambiguity was largely attributable to phenotypically intermediate specimens, mainly from the Mediterranean Basin, and also by the presence of some shared characters amongst the morphospecies, such as the presence/absence of soralia, pseudocyphellae, extrolite composition, and thallus colour.
Phylogenetic tree
Due to the topological conflict between loci, three DNA matrices were used to generate three phylogenetic trees:a concatenated matrix including ITS, IGS, and GAPDH with 134 individuals consisting of 1 774 unambiguously aligned nucleotide position characters, with 83 parsimony informative (Pi) sites;FRBi15 with 93 individuals contained 569 unambiguously aligned nucleotide position characters, with 44 Pi sites; andFRBi16 with 80 individuals had 632 unambiguously aligned nucleotide position characters, with 160 Pi sites.No evidence of recombination events was detected in the concatenated matrix. The resulting tree (Fig. 6) had four well-supported main clades, G (Glabra, yellow), Ko (Kockiana, magenta), NA (North American, blue), and WD (Widely Distributed, red + green + brown). Clade G included only material of B. glabra, appearing as an isolated taxon sister to the other three clades, which showed an uncertain topology between them. Clade Ko included material named B. kockiana and two unidentified specimens, all collected in Alaska (USA). Clade NA comprised the previously recognized North American morphospecies group (the ‘North American endemic species’, Velmala et al. 2014) named as B. friabilis, B. inactiva, B. pikei, and B. pseudofuscescens. While these species were mixed in the tree, the group as a whole was resolved as monophyletic. The WD clade included specimens widely distributed but mainly from Europe (the ‘European and globally distributed species’ group, Velmala et al. 2014) under the names B. capillaris, B. fuscescens (syn. B. chalybeiformis and B. lanestris), B. implexa, B. kuemmerleana, and B. vrangiana. None of these previously recognised species formed a monophyletic group.The phylograms produced using the FRBi15 and FRBi16 markers had a different tree topology, not congruent among them or with that from the preceding concatenated dataset. In the FRBi15 reconstruction (Fig. 1), B. glabra was not represented due to the lack of primer annealing in the PCR process, and the tree could not be rooted. Several well-supported groups were produced, but did not follow any evident geographic, morphological, chemical, or SSR frequency pattern. Bryoria pikei L376 had a sequence with a putative recombination fragment with the B. vrangiana S10 clade in c. 50 % of the total length. This insertion placed the specimen out of the main parental group and it appears as its sister. Although marker FRBi16 (Fig. 1) also produced a well-resolved tree with supported nodes, the clades do not show any phenotypic and/or geographic structure. In this tree reconstruction, B. glabra was an isolated taxon and served to root the tree. In FRBi16 sequences, many putative recombination events were detected, suggesting a reticulate evolution. In both trees in Fig. 1, clade Ko (magenta) was recovered as monophyletic, but embedded between other named morphospecies.
Fig. 1
Phylogenetic relationships in Bryoria sect. Implexae based on FRBi15 and FRBi16 loci. Tree topology depicts the result of the Bayesian Markov chain Monte Carlo (B/MCMC) analysis. Posterior probabilities and bootstrap analysis for the supported nodes (≥ 0.95 and ≥ 70 %) are indicated at the main nodes. Lines connecting clades indicate putative recombination events, with main parents (continuous lines) and minor parents (discontinuous lines). Because the clade insertion in the trees is influenced by the recombination, clades with recombination are depicted with a discontinuous branch line. Note that clades with recombination appear as sister or close to the main parent but tending to be deviated towards the minor parent. — The coloured bar corresponds to the SSR genepool from Fig. 2, with specimens intermediate between two or more genepools in grey. The clades obtained, although well-supported, do not follow any evident geographic, morphologic, chemical, or microsatellite pattern.
STRUCTURE clustering
Of the 18 microsatellite markers, the nine that showed more than 95 % successful amplification across the samples were used (number of haplotypes shown in brackets, Appendix 3): Bi01 (17), Bi03 (6), Bi04 (8), Bi05 (5), Bi10 (8), Bi11 (9), Bi12 (12), Bi14 (6), and Bi19 (5). We allowed a maximum of three missing loci per specimen, a value reached only in seven samples. STRUCTURE was allowed to run to K = 12, but from K = 6 the clustering process started to be uninformative (Fig. 2). The likelihood results of the ΔK analysis (Evanno et al. 2005) indicated three as the most probable number of clusters (likelihood = –1232, ΔK = 2.2), the clades G, NA, and WD (Fig. 2, K = 3). Clade Ko, which appeared isolated in the concatenated phylogenetic tree (Fig. 6), could not be accepted as distinct under a K = 3 hypothesis. However, B. glabra, a morphologically delimited taxon, was not isolated at first in STRUCTURE. This could be attributable to the clustering algorithms being influenced by unbalanced sampling sizes, masking clade Ko, which appeared isolated at K = 6. From K = 4 to K = 6, the new groups appeared mainly inside the WD clade, showing that the samples from Europe were much more diverse than those from North America. Indeed, the NA clade was not split into subgroups even at K = 10. Apart from B. glabra, no other named morphospecies formed an exclusive cluster even at high K values.
Fig. 2
Bayesian inference of population structuring using STRUCTURE v. 2.3.4 (Pritchard et al. 2000, Falush et al. 2003) and nine microsatellite loci in Bryoria sect. Implexae. – Left. Results from the hypothesis of 2–6 clusters. Vertical bars represent specimen assignment probability into a genetic cluster depicted by the colours. Morphospecies names given to the specimens appear at the top. – Right. Detailed columns of the K = 6 hypothesis, the numbers representing the specimens shown in Table 1 to provide a better understanding of the components of each individual. — G = Glabra clade; Ko = Kockiana clade; NA = North American clade; WD = Widely Distributed clade; * = Bryoria pikei specimen 49.
Principal coordinates analysis (PCoA)
The PCoA analysis has a three-dimensional output represented here in two graphs, one comparing axis 1 against 2, and the other 1 against 3 (Fig. 3). The information percentage of each axis was 44.47 %, 15.06 %, and 14.44 %, respectively. Clade G (Fig. 3, yellow) appeared isolated, whereas clade Ko (Fig. 3, magenta) was admixed with NA clade (Fig. 3, blue), forming a single cluster. Clade WD was isolated from the others, but divided into two clusters, one corresponding to the red and brown groups in the K = 6 STRUCTURE output (Fig. 2), and one for the green group. Apart from B. glabra, none of the currently accepted morphospecies formed a defined group. Four reasonably isolated clusters could be distinguished, corresponding to the groups G, WDr (Widely Distributed, Fig. 3, red), WDg (Widely Distributed, Fig. 3, green), and Ko together with NA forming a single cluster.
Fig. 3
Principal Coordinate Analysis (PCoA) of microsatellite data in Bryoria sect. Implexae. Species names according to Velmala et al. 2014 (shape and colours) and STRUCTURE clusters (colours) for each specimen are represented in the three main coordinates. Note that the Ko clade does not appear isolated from the NA clade in any coordinate axis. — G = Glabra clade; Ko = Kockiana clade; NA = North American clade; WDr = Widely Distributed red clade; WDg = Widely Distributed green clade.
Haplotype network
The haplotype network of the concatenated data matrix, coding gaps as missing data, produced 39 haplotypes. Bryoria glabra specimens (Appendix 4, yellow) formed two haplotypes not connected with other members of the network, indicating genetic isolation of this species. One of the haplotypes was composed exclusively of South American specimens, whereas the other contained European, North American, and South American samples. Clade Ko (Appendix 4, magenta) fell into two haplotypes, one including specimens with psoromic acid and identified as B. kockiana, and the other clustering unidentified samples with no substances detected. This group was connected to the NA clade (Appendix 4, blue) by a long branch with 13 mutation steps. The NA clade was separated by nine mutations from the WD clade (Appendix 4, green, red, and brown). The WD green, red, and brown groups split by STRUCTURE (Fig. 2) formed a unique cluster. Four isolated clusters could be distinguished, corresponding to the groups G, Ko, NA, and WD.
Species delimitation programs
The ABGD, PTP, GMYC, and DISSECT programs (Table 3) use different algorithms, and consequently different numbers of putative species may be predicted. The genetic distance method (ABGD) gave the smallest number of putative species, whereas the coalescence methods (especially GMYC) the largest. Analyses also revealed the contribution of each locus to the postulated species delimitation, GAPDH being the most informative and constant marker. DISSECT analysis (Fig. 4) predicted five species corresponding to G, Ko, NA, and WD clades, and specimen B. pikei 5. Although the GMYC analysis also showed the B. pikei 5 specimen as a separate species, it was grouped in the NA clade in the other analyses. DISSECT showed two internal greyish square groups in WD, but they did not correspond exactly to the WDr and WDg groups in Fig. 2, 3 (STRUCTURE and PCoA analyses).
Table 3
Species delimitation analysis results for loci ITS, IGS, GAPDH and the concatenated data matrix in Bryoria sect. Implexae. Brackets indicate groups predicted as conspecific. — G = Glabra clade; Ko = Kockiana clade; NA = North American clade; WD = Wide Distributed clade; WDr = Wide Distributed red clade; WDg = Wide Distributed green clade; pik5 = Specimen Bryoria pikei 5.
Method
ITS
IGS
GAPDH
Concatenated
ABGD
2 = G + (Ko, NA, WD)
2 = G + (Ko, NA, WD)
4 = G + Ko + NA + WD
4 = G + Ko + NA + WD
PTP
2 = G + (Ko, NA, WD)
2 = G + (Ko, NA, WD)
4 = G + Ko + NA + WD
5 = G + Ko + NA + WDr + WDg
GMYCs
4 = G + (Ko, NA, WDg) + WDr + WDr
3 = G + (Ko, WD) + NA
4 = G + Ko + NA + WD
6 = G + Ko + NA + pik5 + WDr + WDg
GMYCm
4 = G + (Ko, NA, WDg) + WDr + WDr
4 = G + (Ko, WD) + NA1 + NA2
4 = G + Ko + NA + WD
5 = G + Ko + NA + WDr + WDg
DISSECT
–
–
–
5 = G + Ko + NA + pik5 + WD
Fig. 4
Similarity matrix from DISSECT analysis performed after clone correction in Bryoria sect. Implexae. Squares represent posterior probability (white = 0, black = 1) of pairs of specimens to belong to the same species. Resulting major groups are delimited by lines, which indicate the clade on the collapsed phylogenetic tree.
Node dates and demographic history
The calibrated maximum clade credibility chronogram for the concatenated data matrix is shown in Fig. 6. As only the ITS mutation rate is estimated in previous studies (Leavitt et al. 2012a, b), a second chronogram was prepared using this locus alone. Results from this analysis have to be treated with caution, as the species tree is not strictly clock-like (B. glabra has a shorter branch), and the ITS mutation rate has been taken from Melanohalea, a lichen-forming genus in the same family. Both analyses produced similar values, and the divergence of the B. glabra lineage was estimated at 6.9 Mya (95 % HPD = 3.5–10.8) in the concatenated matrix analysis, and 6.5 Mya (95 % HPD = 2.2–11.4) in the ITS data alone. The Ko, NA, and WD split was estimated at 1.0 Mya (95 % HPD = 0.3–2.2) from the concatenated matrix and 0.6 Mya (95 % HPD = 0.2–1.5) from the ITS data alone.Bayesian Skyline Plots (Fig. 5, left) indicate a recent population increase in the NA and WD clades. However, the sequences contained few informative mutations and the deepest coalescence was reached in around 700 000 yr, with no population changes detectable further back from this period. Tests of neutrality (Fig. 5, right) are commonly used to support inferences from Bayesian Skyline Plots. As indicated by non-significant Tajima’s D and Fs results, all sampled groups seem in mutation-drift equilibrium, with the exception of the GAPDH locus of the NA clade which had a significant negative D value (Fig. 5
bold). This could indicate a recent expansion or ‘purifying’ selection, as seen in the concatenated Bayesian Skyline analysis, but other loci did not support this hypothesis.
Fig. 5
Bryoria sect. Implexae. – Left. Bayesian Skyline Plots for each clade predicted by the ITS marker and the concatenated loci matrix. The X-axis of each graph represent time (in Myr), and the Y-axis represents the value for the log of the effective population size as relative changes, because generation times in Bryoria species are unknown. Grey shadows indicate the upper and lower 95 % credible intervals. – Right. Results from neutrality tests for each marker and clade, indicating (in bold) any statistically significant deviation from neutrality. — h = number of haplotypes; Fs = Fu’s Fs; D = Tajima’s D.
Markers ITS, IGS, and GAPDH indicate population stability over the recent past for clades NA and WD, with putative even more recent small population expansions. Due to the low variability of the loci, and the putative loss of demographic signals, this hypothesis is not confirmed by this analysis.
Integrated assessment of datasets
Depending on the analysis, different numbers of putative species were suggested, ranging from four to six (Table 4). All analyses, however, confirmed that the combination of morphological and chemical characters generally used for species circumscription in the complex was inadequate. GAPDH, despite its low variability, was the only marker tested that supported species-rank assignations for the clades G, Ko, NA, and WD (Table 3). ITS, one of the most used loci for DNA barcoding in lichen-forming fungi, did not unambiguously distinguish those clades. The new markers FRBi15 and FRBi16, despite their higher variability, showed inconclusive results and putative recombination events. The microsatellite data (Fig. 2) supported the DNA sequences results and reflected internal variability not revealed in our sequence data, showing that the WD cluster was much more diverse than NA, which had a particularly low diversity.
Table 4
Summary of the number of putative species suggested by the different methods used for each dataset in Bryoria sect. Implexae.
Method
Data
Figure / reference
Number of putative species
Traditional concept
DNA sequences and phenotypic
Velmala et al. (2014)
12
Chemical
Phenotypic
Appendix 2 Left
c. 4
Morpho-chemical
Phenotypic
Appendix 2 Right
Not conclusive
Phylogeny
DNA sequences of ITS, IGS, and GAPDH
Fig. 6
4 = G + Ko + NA + WD
Phylogeny
DNA sequences of FRBi15
Fig. 1
Not conclusive
Phylogeny
DNA sequences of FRBi16
Fig. 1
Not conclusive
STRUCTURE
Microsatellites
Fig. 2
5 = G + Ko + NA + WDr + WDg
PCoA
Microsatellites
Fig. 3
4 = G + (Ko, NA) + WDr + WDg
Haplotype Network
DNA sequences of ITS, IGS, and GAPDH
Appendix 4
4 = G + Ko + NA + WD
ABGD
DNA sequences of ITS, IGS, and GAPDH
Table 3
4 = G + Ko + NA + WD
PTP
DNA sequences of ITS, IGS, and GAPDH
Table 3
5 = G + Ko + NA + WDr + WDg
GMYCs
DNA sequences of ITS, IGS, and GAPDH
Table 3
6 = G + Ko + NA + pik5 + WDr + WDg
GMYCm
DNA sequences of ITS, IGS, and GAPDH
Table 3
5 = G + Ko + NA + WDr + WDg
DISSECT
DNA sequences of ITS, IGS, and GAPDH
Fig. 4
5 = G + Ko + NA + pik5 + WD
TAXONOMY
sect. (Gyeln.) Brodo & D. Hawksw., Opera Bot. 42: 114. 1977Basionym. Bryopogon sect. Implexae Gyeln., Feddes Repert. Spec. Nov. Regni Veg. 38: 223, 238. 1935.Type species. Bryoria implexa (Hoffm.) Brodo & D. Hawksw. 1977. ≡ Usnea [unranked] implexa Hoffm. 1796. = Bryoria fuscescens (Gyeln.) Brodo & D. Hawksw. 1977; but see below.Species with a fruticose, hair-like, subpendent to mainly pendent thallus, lateral spinules or spinulose branches absent, whitish grey to brown or black, often paler in the basal parts. Angles between branches variable, acute to obtuse or even rounded. Pseudocyphellae absent or present, then frequently inconspicuous, ± fusiform, concolorous or whitish. Soralia absent or present, tuberculate or fissural, white to dark. Isidia or isidioid spinules absent. Apothecia mainly absent, if present, usually afunctional. Chemistry varied, with no detectable or with one or a combination of major substances, including alectorialic, barbatolic, connorstictic, fumarprotocetraric, gyrophoric, norstictic, protocetraric, psoromic and possibly salazinic acids, atranorin, and chloroatranorin. Photobiont Trebouxia ‘hypogymniae’ (Lindgren et al. 2014).Notes — Most species included in Brodo & Hawksworth (1977) under Bryoria sect. Implexae were transferred to other sections in Myllys et al. (2011). In the light of our results (but see the Discussion later), Bryoria sect. Implexae includes the four species treated below. Comments on particular morphological or chemical traits that may be helpful for distinguishing these taxa are given under each species. Nevertheless, nearly all cited characters are shared by different taxa, so they can be interpreted as ‘cryptic’. The species names adopted here are epitypified by sequenced material here in order to fix their identities at the molecular level. This epitypification is essential to fix the application of these names as no DNA sequences are available and cannot be obtained from old type material of most species names. The old types cannot therefore be critically identified for purposes of the precise application of the names and so epitypes may be designated (Turland et al. 2018: Art. 9.9 and Ex 9).As molecular data are necessary for unambiguous species level identification in the taxonomy proposed here, we recommend using the collective ‘Bryoria fuscescens complex’ when referring to material lacking molecular data.(Gyeln.) Brodo & D. Hawksw., Opera Bot. 42: 83. 1977Basionym. Alectoria fuscescens Gyeln., Nytt Mag. Naturvidensk. 70: 55. 1932, nom. cons. (cf. Hawksworth & Jørgensen 2013).Synonyms. Lichen chalybeiformis L., Sp. Pl. 2: 1155. 1753, (nom. cons.) nom. rej. against Bryoria fuscescens (cf. Hawksworth & Jørgensen 2013).Bryoria chalybeiformis (L.) Brodo & D. Hawksw., Opera Bot. 42: 81. 1977.Usnea [unranked] implexa Hoffm., Deutschl. Fl., Zweiter Teil: 134. 1796.Bryoria implexa (Hoffm.) Brodo & D. Hawksw., Opera Bot. 42: 121. 1977.Parmelia jubata β. [var.] capillaris Ach., Methodus, Sectio post.: 273. 1803.Bryoria capillaris (Ach.) Brodo & D. Hawksw., Opera Bot. 42: 115. 1977.Alectoria jubata χ. [var.] lanestris Ach., Lichenogr. Universalis: 593. 1810.Bryoria lanestris (Ach.) Brodo & D. Hawksw., Opera Bot. 70: 88 1977.Alectoria kuemmerleana Gyeln., Nytt Mag. Naturvidensk. 70: 49. 1932.Bryoria kuemmerleana (Gyeln.) Brodo & D. Hawksw., Opera Bot. 42: 155. 1977.Alectoria prolixa var. subcana Nyl. ex Stizenb., Ann. Naturhist. Mus. Wien 7: 129. 1892, nom. rej. against Bryoria fuscescens (cf. Hawksworth & Jørgensen 2013).Bryoria subcana (Nyl. ex Stizenb.) Brodo & D. Hawksw., Opera Bot. 42: 91. 1977.Alectoria vrangiana Gyeln., Magyar Bot. Lapok 31: 46. 1932.Bryoria vrangiana (Gyeln.) Brodo & D. Hawksw., Opera Bot. 42: 97. 1977.Type specimens. FINLAND, Tavastia austr., Hollola, ad truncos Pini locis apricioribus in silva, Sept. 1882, J.P. Norrlin (Norrlin, Herb. Lich. Fenn. No. 46) (BP 33947 – lectotype designated by Hawksworth 1972: 217). — FINLAND, Etelä-Savo, Savitaipale, 600 m NW of Mustapää, 61, N1721° E27,6900°, 2005, L. Myllys 464 (HA.H9209715 (L139)) – epitype designated here, MycoBank MBT381730.Nomenclature — A large number of species rank names belong to this group, and are synonymised, but these have not been epitypified with sequenced material. Further information on synonyms and type materials can be seen in Hawksworth (1972), Brodo & Hawksworth (1977) and Velmala et al. (2014). Although no samples of Bryoria austromontana have been studied, the published description and images (Jørgensen & Galloway 1983) suggest this taxon also belong here.The earliest species rank epithets amongst these are chalybeiformis dating from 1753 (Lichen chalybeiformis), and implexa dating from 1799 (Usnea implexa). The former has been rejected against Bryoria fuscescens, but not against other species names apart from B. subcana (Hawksworth & Jørgensen 2013). A proposal to add the four earlier names Alectoria capillaris, Usnea implexa, A. kuemmerleana, and A. lanestris to the two against which Alectoria fuscescens is already conserved is being prepared separately. Protection against A. vrangiana is not required as it appeared in the same work as A. fuscescens. While the proposal is under discussion, the name B. fuscescens should be adopted in accordance with Rec. 14A.1 (Turland et al. 2018).We refrained from epitypifying and taking up any of the earlier and potentially competing names by epitypification primarily as the name B. fuscescens is the most commonly used species name in the complex, is well established, the most widely used* and is already conserved over two earlier species names in the complex. In addition, all the other names have been associated with particular morphotypes or chemotypes since the 1970s, and so their use might be mistaken as applying to a taxon with those particular traits.If the proposal for rejection of the previously mentioned competing synonyms is not accepted, the principle of priority would rule the use of the earliest and not already rejected, validly published name at the species rank, i.e., Usnea implexa (and then the combination Bryoria implexa), which would require epitypification by sequenced material in order to fix the precise application of that name. The species was first described from Germany but with no named locality, and neotypified by an unlocalised and undated specimen in Hoffmann’s herbarium in Moscow (Hoffmann 8569, MW) which may be part of the original material from Germany or have been collected later and perhaps in Russia (Hawksworth 1969a). As the neotype contained psoromic acid, and the epithet has therefore been applied to that chemotype, a potential sequenced epitype should represent that chemotype and ideally also have been collected in Germany. No such specimen was available to us during this study.Distribution — Widely distributed, known from cool temperate to boreal and arctic areas of Europe, Asia, North America, and Africa. There are also reports from Antarctica, Oceania, and South America, but material from those regions has not yet been studied molecularly and so we cannot confirm that they belong to this complex.Notes — Bryoria fuscescens is highly variable in morphology and chemistry, and many of the analysed specimens develop soralia. Further, atranorin, which is not normally detectable in the other three species accepted here, is frequently found in concentrated extracts from both sorediate and esorediate morphs.(Motyka) Brodo & D. Hawksw., Opera Bot. 42: 86. 1977Basionym. Alectoria glabra Motyka, Fragm. Florist. Geobot. 6: 448. 1960.Type specimens. USA, Washington, Olympic Peninsula, Clallam Co., Hurricane Ridge, 5800 ft, on trunk of Abies lasiocarpa, 24 July 1950, B.I. Brown & W.C. Muenscher 129 (US – holotype). — USA, Alaska, Mainland, Valley between the Bucher and Gilkey Glaciers, southern end of subalpine valley, on east side of creek running through valley, subalpine forest, N58°47’20.12" W134°30’0.10", 773 m elevation, on Tsuga mertensiana twigs, 4 Aug. 2011, K.L. Dillman 4Aug11:1 (UBC (L406)) – epitype designated here, MycoBank MBT381731.Distribution — Known from northern Europe (Scandinavia), and North and South America. In North America, it is most abundant in the Pacific North-West.Notes — Distinguishing features in well-developed specimens are the brownish thallus with a regular branching pattern, generally with obtuse and rounded angles between the branches, and broad oval and usually white soralia. It is, however, difficult to separate poorly developed or small specimens conclusively, so molecular sequences are recommended for unambiguous identifications. Only fumarprotocetraric and protocetraric acids have been detected in this species, and these are characteristically produced in the soralia.The Alaskan specimen is selected as the epitype here as sequences are available from all loci, whereas the material we have sequenced from Washington state (type locality) only has data on the ITS locus.Velmala, Myllys & Goward in Velmala et al., Ann. Bot. Fenn. 51: 361. 2014Type specimen. USA, Alaska, Fairbanks, North Star Borough, 26 July 2011, D. Nossov 20019-1 (UBC (L394) – holotype).Distribution — Known only from Alaska (USA) and British Columbia (Canada), on conifer branches.Notes — Few specimens of this species have so far been studied, and these are characterised by the absence of any whitish grey tone in the thallus, the lack of soralia, and greyish to brown branches with conspicuous, white to concolourous, broad, elongate-fusiform, sometimes slightly raised pseudocyphellae. It lacks TLC-detectable substances or produces psoromic acid. The not validly published designation Alectoria krogii D. Hawksw. 1972 may be synonymised here.(Gyeln.) Brodo & D. Hawksw., Opera Bot. 42: 127. 1977Basionym. Alectoria pseudofuscescens Gyeln., Ann. Hist.-Nat. Mus. Natl. Hung. 28: 283. 1934, and Rev. Bryol. Lichénol. 7: 51. 1934.Synonyms. Bryoria friabilis Brodo & D. Hawksw., Opera Bot 42: 118. 1977.Bryoria pikei Brodo & D. Hawksw., Opera Bot 42: 125. 1977.Bryoria inactiva Goward et al., Ann. Bot. Fenn. 51: 360. 2014.Type specimens. USA, Oregon, Benton County, Corvallis, on old apple trees, Dec. 1931, F.P. Sipe 669 (BP 33958 – holotype of Alectoria pseudofuscescens). — CANADA, British Columbia, 25 Sept. 2006, T. Goward 07-02-2011 (UBC (S222) – epitype selected here, MycoBank MBT381732; British Columbia, Clearwater Valley, 0.5 km S of Philip Creek, ‘Edgewood West’, 715 m, 9 Nov. 2011, T. Goward 11-61 (UBC (L347) – holotype of Bryoria inactiva).Nomenclature — A number of species rank names are synonymised to this taxon, but these have not been epitypified with sequenced material. All these names, however, are later in date than pseudofuscescens, and so could not have priority over that name. Further information on type materials can be seen in Brodo & Hawksworth (1977) and Velmala et al. (2014). Although no samples of Bryoria salazinica have been studied at molecular level, the published description and images (Brodo & Hawksworth 1977) suggests this taxon also belong here.Distribution — Only known from North America, growing on bark, branches, rock or soil.Notes — Characterised by the absence of soralia and detectable atranorin.
DISCUSSION
Phylogenetic relationships
Species concepts in Bryoria sect. Implexae have previously been based primarily on well-characterised northern European and North American specimens (Hawksworth 1972, Brodo & Hawksworth 1977, Velmala et al. 2014). Velmala et al. (2014) recognised 11 species on the basis of morphological and chemical characters, but many of these were not supported by molecular data, and different species names were accepted for taxa that could not be distinguished molecularly. We discovered that these demarcations broke down when specimens from more southern European populations were incorporated. This is shown in a phenetic analysis using only phenotypically diagnostic characters (Appendix 2), where the resulting groups are not resolved as clear-cut morphospecies. Indeed, any character previously used in the group could be used to define the three lineages of the Bryoria fuscescens complex (Fig. 6).Sexual structures are of major importance in species identification in fungi, but here the rarity of apothecial production has hampered their study in most Bryoria species. Any such features found would in any case be of limited practical value in identification as nearly all samples lack apothecia. Extrolite composition has been accorded a major role in species delimitation in the complex, but many of the substances that were considered to be of diagnostic value are biosynthetically closely related, being produced by the same gene cluster (pks genes; Keller & Hohn 1997), and may be environmentally influenced (Myllys et al. 2016, Lutsak et al. 2017).Integrative taxonomy, rather than phylogenies based only on neutral markers, are increasingly being used to resolve complex taxonomical groups (e.g., Dayrat 2005, Will et al. 2005, Lumley & Sperling 2011, Zamora et al. 2015, Caparrós et al. 2016). Microsatellites are also now widely used in intraspecific population studies because of their high variability (Widmer et al. 2012, Dal Grande et al. 2014), and in species complexes with diffuse genetic barriers, microsatellite data can improve DNA sequence resolution (Lumley & Sperling 2011, Vanhaecke et al. 2012). It is generally assumed that DNA sequence data reflect the evolution of the species, but these data only reflect the history of the studied loci, which may sometimes be different from the species history overall (Nichols 2001). In this case, we demonstrated that traditionally used loci (ITS, IGS, and GAPDH) and microsatellite data reveal similar clades, whereas other intergenic loci (FRBi15 and FRBi16) produced discrepant but statistically supported lineages. These incongruences may be due to recombination, hybridisation, or incomplete lineage sorting, as documented in many other species groups (e.g., Jakob & Blattner 2006, McGuire et al. 2007, Edwards et al. 2008, Stewart et al. 2014). In lichen-forming fungi, outcrossing and recombination have been demonstrated, for example, in Lobaria pulmonaria (Zoller et al. 1999, Singh et al. 2012, Keller & Scheidegger 2016), Letharia (Kroken & Taylor 2001a, b, Altermann et al. 2014), and Cladonia (Steinová et al. 2013).Apothecia are usually absent in Bryoria sect. Implexae, and even when present may not contain mature spores. If cryptic sexuality is not occurring, hybridization is unlikely to provide an explanation of our data. In the absence of sexual reproduction, any recombination is improbable, although some fungi lacking sexual structures show recombination events attributable to parasexual cycles (Schoustra et al. 2007). We did detect signals suggesting putative recombination in the FRBi loci, but not in the standard three loci used in the taxonomy adopted here. Recombination signals may reflect some mitotic recombination, actual or ancient sexual reproduction (Douhan et al. 2007) or be merely false positives produced by chance production of similar sequences. In any case, recombination alone is insufficient to explain all the discordances found. For example, only one putative recombination event was detected in FRBi15, and disentangling the FRBi16 recombination points is insufficient to obtain the topology of the three-locus phylogeny. Incongruences may also be caused by the analysis of different paralogs of FRBi15 and FRBi16 amplified with the new primers, but this seems improbable, as no paralogs have been detected in the SSRs of these loci (Boluda et al. unpubl.). However, our results indicate recent diversification and large effective population sizes in this lichenised complex. Thus, incongruences amongst loci seem rather attributable to incomplete lineage sorting.The different putative species concepts generated by the species delimitation programs (Table 3) may not be only due to the different algorithms applied, but also because some of the programs were designed for use in single-locus phylogenies (e.g., ABGD, PTP, or GMYC). Nevertheless, all the clustering analyses showed a tendency to distinguish four groups, G, Ko, NA, and WD (Table 4; Fig. 6). STRUCTURE was unable to define these groups until reaching the K = 6 hypothesis, which can be attributed to the highly unbalanced sampling sizes; the analyses shows that the WD cluster is much more variable than NA, which was not divided into subgroups until K = 10.Specimen 49 (identified as B. pikei, Fig. 2 asterisk), probably reflects the impossibility of unequivocally distinguishing that species from B. capillaris; however, as sequence amplification of this specimen failed, we cannot determine if this mismatch was due to misidentification or DNA contamination. Haplotype network analyses have been extensively used in infraspecific population and less frequently closely related species groups (Houbraken et al. 2012, Pino-Bodas et al. 2016). Even this type of nested clade phylogeographic analysis has some critics (e.g., Knowles 2008, Templeton 2009). The resulting network in the present case is concordant with that obtained from other analyses. If two DNA barcoding standard marker networks are obtained from a single analysis with a 95 % parsimony connection limit, members of each network might be considered as different species (Hart & Sunday 2007), showing the clear isolation of B. glabra from the other taxa in the complex. In the case of connections with several mutation steps, as between clades Ko, NA, and WD, taxon delimitation is below the species level, but in any case indicates some degree of genetic isolation.The relationships amongst the Ko, NA, and WD clades remains unresolved, indicating that the evolutionary history may be too complex to be adequately captured by dichotomous phylogenies based on a few neutral markers. Moreover, the putative presence of shared ancestral polymorphisms amongst the clades may be producing incompatible topologies, which result in clades with low support.We also performed analyses to estimate changes in past population sizes, which may have affected current clade diversity. Genealogies of most plant and animal species coalesce between 2.58–0.01 Mya (Grant 2015), and our estimated intervals are within this range. However, in our case the dates are relatively recent, with the oldest coalescences at 0.7 Mya. A flat graphic is generally interpreted as population stability but can also be due to a lack of detection power produced by small sample sizes. Moreover, a small rise in the curve near the present, seen in Fig. 5, can be a consequence of the random sampling of the MCMC haplotype trees (Grant 2015); this result must therefore be treated with caution. Our sequences do not bear imprints of ancient population history, but rather more recent population growth, for example by extensions northwards in post-glacial periods. The loss of information may also arise from bottlenecks (i.e., a marked reduction in population size), local extinctions, and subsequent recolonization. Additionally, the use of genes with low levels of polymorphisms, may impede a robust reconstruction of population sizes through time.
Species concept
Some species delimitation analyses, such as STRUCTURE, GMYC, PTP, can overestimate the number of taxa meriting formal recognition, particularly when sampling is uneven or in species with a strong intraspecific genetic structure (Altermann et al. 2014, Modica et al. 2014, Alors et al. 2016, Del-Prado et al. 2016). ABGD, in contrast, has been considered as rather conservative, less prone to species overestimation and less sensitive to unbalanced sampling. While that method only detects discontinuities in DNA sequence variation (Puillandre et al. 2011), it can also be expected to fail in species with strong population genetics structures, for example ones containing exclusive haplotypes. All species delimitation programs will provide a number of reasonably discrete groups (‘evolving lineages’) that should be evaluated for consideration as meriting species rank, but the decision has to be by taxonomists with experience in the group concerned. Some of our analyses suggest that groups such as WDr, WDg, or pik5 might merit species rank, but our experience, together with the results from other analyses shown here, leads us to reject this hypothesis. We conclude that the most pragmatic solution, supported by the general trend of the results from the different analyses we performed is to consider clades G, Ko, NA, and WD as the species Bryoria glabra, B. kockiana, B. pseudofuscescens, and B. fuscescens, respectively.Clade age can contribute to decisions as to species limits. Divergence time estimates can be robust if the analyses are performed with well-resolved phylogenies and can incorporate fossil calibration points, as in some vertebrates (Perelman et al. 2011). Contrarily, in lichenised fungi, fossils are rare and in many cases enigmatic or with ambiguous relationships (Thomas et al. 2014, Hawksworth 2015, Kaasalainen et al. 2015). In addition, generation times can be expected to be different among species, as is the case with nuITS locus mutation rates between herbaceous and woody plants, or even a difference of almost an order of magnitude between different plant genera (Kay et al. 2006). Here we used a nuITS mutation rate estimated from Melanohalea, a genus in the same family (Leavitt et al. 2012a), species of which frequently grow with Bryoria and reproduce asexually as well as sexually. The split of B. glabra from the other taxa in the B. fuscescens complex is estimated at c. 6.9 Mya, and clearly separated from the much later divergence of the other three species estimated at c. 1 Mya (0.6 Mya if only the nuITS locus is used). This contrasts with other lichenised species considered of recent origin, estimated around 2.5–5 Mya (Pliocene), with any Pleistocene speciation event rare and always older than 1 Mya (Amo de Paz et al. 2012, Leavitt et al. 2012a, b, Molina et al. 2016). As the three B. fuscescens complex clades seem to have diverged more recently, the extent of their reproductive isolation is uncertain, and the discovery of intermediate lineages amongst other named species from unsampled geographical regions, such as continental Asia remains possible.In the absence of supporting phenotypic, geographic, or ecological differences, the recent divergence, and the possibility of incomplete lineage sorting, clades Ko, NA, and WD may be considered as conspecific evolving lineages. It is, however, important to recognise the lineages formally in order to facilitate their conservation by enabling their threat status to be assessed by IUCN criteria. We decided not to adopt the rank of subspecies as that is now hardly used in mycology, and then not in any consistent way; traditionally this rank was used in plants for morphologically distinguishable populations with geographical differences and where intermediates occurred where they were sympatric (Stuessy 2009).The formal recognition of cryptic lineages at species level, as suggested by our analyses, emerges as the most appropriate solution. Cryptic speciation is now recognised as a common phenomenon in Parmeliaceae, and our results are in accord with other studies in which molecular markers in combination with statistical tools revealed genetically distinct lineages previously hidden under a single taxon name in this family (e.g., Singh et al. 2015, Alors et al. 2016, Del-Prado et al. 2016, Divakar et al. 2016, Leavitt et al. 2016). Further, this solution is in line with the increasing need to formally recognize cryptic species-level lineages in all fungi (Hibbett 2016); indeed, cryptic speciation may mean that there are on average ten or more fungi masked in formally named species (Hawksworth & Lücking 2017).Of the lineages recognised here, only the WD clade emerged as cosmopolitan, occurring in Europe, Asia, North America, and Africa (Appendix 5). NA and Ko have been collected so far only in North America, despite our extensive sampling in Europe (Appendix 5). Further sampling, especially in South America, Asia, and Africa, is needed before any finer-scale biogeographic patterns might be detected.The practical issue of naming older museum specimens and material in ecological surveys could be resolved by recognising the three groups as species within a broad concept, such as an aggregate, complex, or adding ‘s.lat.’. We considered commending the adoption of the suffix ‘agg.’ for material when precise molecular species identifications cannot be made. While this has been done in a few other groups of fungi (e.g., Parnmen et al. 2013, Pažoutová et al. 2015), ‘complex’ has come to be used more widely and was strongly favoured at the Cryptic Speciation in Fungi symposium in Utrecht in September 2017 (report awaited). We therefore suggest the use of ‘complex’ here but recognise some may prefer to use ‘agg.’ or ‘s.lat.’. Where DNA samples can be obtained and analysed, we recommend use of the GAPDH locus, as all the other tested markers are not able to distinguish with confidence the three species we recognize in the complex.
Infraspecific phenotypic variation
While our results support rejection of the morphospecies concept in this group of lichens, two main phenotypes can nevertheless often be distinguished by the naked eye in the field:the pale grey ‘capillaris’ morphotype (including B. capillaris and B. pikei, Fig. 6b); andthe fuscous brown to dark brown ‘fuscescens’ morphotype, in which most other species names are placed (Fig. 6a, c).The chemical characters are not always checked by field workers, and while the ‘capillaris’ morphotype typically has benzyldepsides, the ‘fuscescens’ morphotype lacks those compounds and has fumarprotocetraric acid or various depsidones. However, there are dark morphs with benzyldepsides (once called f. fuscidula), and pale grey ones with fumarprotocetraric acid (e.g., B. subcana) or other depsidones (e.g., B. kuemmerleana). It is conceivable that the two morphotypes originated before the separation of B. pseudofuscescens and B. fuscescens, as both colour variants and chemistries appear in both taxa. This phenomenon cannot be explained by a simple ongoing speciation event in which one lineage has originated new adaptations, but is still not isolated from the parental lineages, as neither are monophyletic in a paraphyletic clade.The difference in phenotype cannot be attributed to different algal partners as all material in the complex shares the same species and even in many cases the same nuITS haplotypes of Trebouxia (Lindgren et al. 2014, Boluda et al. unpubl.). Further, as we used neutral markers to detect gene-flow gaps between lineages, the phenotypes are also not the result of genetic isolation, and other possibilities must be considered.It has recently been reported that yeast morphs of the lichenicolous and gall-forming basidiomycete genus Cyphobasidium can be abundant in or on the outermost cortical tissues of Bryoria species (Spribille et al. 2016). Spribille et al. (2016) reported a possible relation between Cyphobasidium yeast abundance and vulpinic acid production in two other species of Bryoria, B. fremontii and B. tortuosa, and also visualised these yeasts in material identified as B. capillaris phenotypes. Contrary to the claims of Spribille et al. (2016), these fungi do not appear to be an integral part of the mutualism (Oberwinkler 2017). It is, however, feasible that the yeasts cells are able to develop to a greater extent in ‘capillaris’ morphotypes as the cortices can have lumpier polysaccharide deposits than do those of ‘fuscescens’ (Hawksworth 1969b, Boluda et al. 2014, Esseen et al. 2017). How the occurrence of yeast morphs of these lichenicolous fungi in the surface of the cortex could possibly determine colour morphotypes is obscure.Material referred to the ‘capillaris’ and ‘fuscescens’ phenotypes has been reported to show slight differences in water holding capacity (Esseen et al. 2015), and also the pigments may provide protection against excesses of light (Färber et al. 2014). Further, in southern Europe particularly, the ‘capillaris’ phenotype tends to be restricted to humid, shaded, and protected or undisturbed environments than the ‘fuscescens’ one, something already recognised by Motyka (1964). Additionally, in northern Europe, dark specimens containing barbatolic acid are much more common than in southern Europe, where they are extremely rare (Myllys et al. 2016, Esseen et al. 2017). As both phenotypes can grow side by side and even intermixed on the same trees, where environmental conditions must be identical, ecological plasticity has to be discounted. However, some unknown epigenetic modification could perhaps have a role in that process, as once a metabolic pathway is activated or silenced, it may be hardly modifiable under more or less neutral environmental conditions, transferring the phenotypes to the clonal offspring. Specimens with dark thalli, containing barbatolic acid, or with pale thalli with traces of barbatolic and also containing other extrolites, could represent transitional specimens.Molecular and morphological rates of divergence may sometimes be uncoupled. Incomplete lineage sorting arises when an ancestral polymorphism persists through a speciation event and each polymorphism can lead to different alleles being carried among descendants (Maddison 1997, Hartl & Clark 2007). Consequently, different tree topologies may be obtained depending which specimens or loci are used. Rosenberg (2003) has shown that 5.3N generations are needed for a species to acquire monophyly at 99 % of its loci given that all loci in the sister species are also monophyletic. That indicates that for a species of 1 M individuals with a generation time of 10 yr, the full monophyly will only be reached 50 M years after speciation, whereas only around 1 000 yr may be needed for species with small populations (Naciri & Linder 2015). Incomplete lineage sorting may be frequent in closely related taxa or during a speciation process (Hobolth et al. 2011, Blanco-Pastor et al. 2012, Saag et al. 2014, Naciri & Linder 2015), as may be considered in our case. The topological incongruence observed among the standard loci, FRBi15 and FRBi16, supports the incomplete lineage sorting hypothesis as one of the main reasons explaining why morphospecies are not monophyletic. While neutral markers are useful for understanding gene-flow patterns, adaptive markers provide the evolutionary pressure that contributes to speciation (Emelianov et al. 2004, Hey 2006, Holderegger et al. 2006). As adaptive markers are under natural selection, certain alleles can be present in some morphospecies and absent in others, even if there is gene flow amongst them (Lumley & Sperling 2011). The use of phylogenomic datasets may provide a more accurate and supported phylogenetic reconstruction, especially if the appropriate scale of loci variability is selected from all the genome (Leavitt et al. 2016). However, if there are high levels of incomplete lineage sorting, it might not be expected that morphospecies would appear forming supported clades. Nevertheless, genomic data may reveal few mutations linked to certain morphospecies, which would be producing adaptive traits. Darwin’s finches are an iconic example of a rapid speciation process, in which there is a mismatch between the phylogenetic species concept and phenotype-based taxonomy; in that case, genomic studies have detected specific loci subjected to selection pressure, which are directly related with the development of taxon-specific phenotypes (Lamichhaney et al. 2015). In Bryoria, supposed adaptive traits may be influenced by the genes involved in the production of certain extrolites or in the epicortical substances (Boluda et al. 2014), which may produce differential selection pressure for each morphotype, at least in some environments. The process might be similar to that of natural selection of the pale and melanic morphs of the Peppered Moth (Biston betularia) in Europe (Majerus 2009), impeding the fixation of a single morphotype in all populations. In our case also, high levels of incomplete lineage sorting mixed with a few phenotypically important genes under variable degrees of selection in different environments, may explain the mismatch observed between phenotypes and genotypes.
Appendix 1
Chemical, geographical and morphological characters used of Bryoria sect. Implexae samples in the phenogram reconstruction (Fig. 1).
Fig. 6
Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5