Introduction
The magnitude of fungal diversity has been the subject of much attention, with estimates varying by at least two orders of magnitude, but seemingly converging on numbers from two to four million (Hawksworth & Lücking Reference Hawksworth and Lücking2017). Among the most important factors that have recently driven figures upwards have been discoveries of cryptic diversity in fungi, which have revealed a remarkable amount of overlooked morphological divergence (Bensch et al. Reference Bensch, Groenewald, Dijksterhuis, Starink-Willemse, Andersen, Summerell, Shin, Dugan, Schroers and Braun2010; Lücking et al. Reference Lücking, Dal-Forno, Sikaroodi, Gillevet, Bungartz, Moncada, Yánez-Ayabaca, Chaves, Coca and Lawrey2014; Moncada et al. Reference Moncada, Lucking and Suarez2014). In lichen fungi, the Parmeliaceae has been the subject of a great deal of research attention, with a recent focus on checking morphological concepts with barcoding (e.g. Crespo & Lumbsch Reference Crespo and Lumbsch2010; Molina et al. Reference Molina, Del-Prado, Divakar, Sánchez-Mata and Crespo2011; Divakar et al. Reference Divakar, Leavitt, Molina, Del-Prado, Lumbsch and Crespo2015). Barcoding for specimen identification in fungi using the nuclear ribosomal internal transcribed spacer (ITS) is effective across nearly 75% of investigated ascomycetes (Schoch et al. Reference Schoch, Seifert, Huhndorf, Robert, Spouge, Levesque, Chen, Bolchacova, Voigt and Crous2012), including the Parmeliaceae (Crespo & Lumbsch Reference Crespo and Lumbsch2010). The widely distributed species P. saxatilis s. lat. has been the subject of several studies showing that the broad species concept in wide use refers to a suite of morphologically similar species (a species complex or aggregate), that differ genetically and also appear to subtly differ in morphology and chemistry (Divakar et al. Reference Divakar, Leavitt, Molina, Del-Prado, Lumbsch and Crespo2015; Thell et al. Reference Thell, Tsurykau, Persson, Hansson, Åsegård, Kärnefelt and Seaward2017; Ossowska et al. Reference Ossowska, Guzow-Krzemińska, Dudek, Oset and Kukwa2018).
Parmelia saxatilis s. lat. is one of the most common and widely distributed macro-lichen species in the world, existing in both hemispheres including Antarctica (Molina et al. Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004). Crespo et al. (Reference Crespo, Molina, Blanco, Schroeter, Sancho and Hawksworth2002) discovered that a genetically and geographically distinct lineage exists within P. saxatilis s. lat. in Mediterranean regions. Subsequently, further evidence supported a separate species, named P. serrana Crespo et al. (Molina et al. Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004). Another monophyletic clade (Clade IV, sensu Crespo et al. Reference Crespo, Molina, Blanco, Schroeter, Sancho and Hawksworth2002) within the so-called Atlantic population of P. saxatilis has been described as P. ernstiae (Feuerer & Thell Reference Feuerer and Thell2002). Material studied by Feuerer & Thell (Reference Feuerer and Thell2002) of P. ernstiae can be distinguished from P. saxatilis and P. serrana based on its nucleotide sequence and also its morphology since studied specimens had a thick layer of crystals (pruina) on the upper surface and a mixture of both isidia and lobulae. The remaining members of P. saxatilis in Europe (referred to as the Atlantic population by Crespo et al. Reference Crespo, Molina, Blanco, Schroeter, Sancho and Hawksworth2002) are now known as P. saxatilis (L.) Ach. s. str. (henceforth referred to as P. saxatilis).
Further study of additional material shows that the degree of chemical and morphological similarity between the cryptic members of the P. saxatilis complex makes them very hard to distinguish without molecular identification. Although they are believed to have different climatic preferences which could assist with field identification (Molina et al. Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004), the true extent of their habitat ranges are not yet fully understood, particularly outside the original study area (Spain) (Thell et al. Reference Thell, Elix, Feuerer, Hansen, Kärnefelt, Schüler and Westberg2008). Recent work in the Nordic countries suggests that P. ernstiae is widespread but that P. serrana is more rarely encountered (Thell et al. Reference Thell, Tsurykau, Persson, Hansson, Åsegård, Kärnefelt and Seaward2017). It is suspected that there are further cryptic species hidden within P. saxatilis s. lat. (Molina et al. Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004; Thell et al. Reference Thell, Elix, Feuerer, Hansen, Kärnefelt, Schüler and Westberg2008; Divakar et al. Reference Divakar, Leavitt, Molina, Del-Prado, Lumbsch and Crespo2015) and that these will continue to be revealed as we develop our understanding of the genetic diversity within the complex. Indeed, additional cryptic species have already been discovered but their distributions are believed to be restricted to regions not including the UK. For example, P. mayi Divakar et al., is thought to be restricted to the northern Appalachian mountain territories of North America (Molina et al. Reference Molina, Del-Prado, Divakar, Sánchez-Mata and Crespo2011). The discovery of cryptic species within the P. saxatilis complex demonstrates our poor understanding of their ecology and of the extent of morphological diversity within and among these species.
Lichen diversity can differ greatly over many scales (Ellis Reference Ellis2012). To investigate niche preferences and incorporate the different habitat preferences at different scales, lichen samples were collected from five different climatic regions across Scotland, from rocks and different tree species, and from different microhabitats on each tree. The aims of this project were to: 1) identify which members of the P. saxatilis aggregate are present in Scotland; 2) identify any distinguishing morphological characteristics of the P. saxatilis s. lat. aggregate species, based on hypotheses proposed by previous work; 3) investigate the habitat and climatic preferences of each species. Molecular techniques were used to identify the collections to species level a priori. Habitat preferences and morphological characters were then investigated to characterize the species a posteriori. We use multivariate statistics to assess whether ecological, morphological and/or chemical correlates can be used to help confidently and consistently identify barcoded specimens.
Methods
Lichen specimens were collected systematically from five locations across a steep climatic gradient in Scotland: Glen Creran and Benmore (Argyll and Bute), Loch Earn and Dunkeld (Perthshire) and Ballater (Aberdeenshire) (Fig. 1). Locations differ in climate variables including average annual minimum temperature and average annual rainfall (Table 1). Specimens collected opportunistically from additional sites in Stirlingshire, near the middle of the climate gradient, were also examined.
Sampling locations in Scotland for lichen specimens from the Parmelia saxatilis aggregate. GC = Glen Creran, BM = Benmore, LE = Loch Earn, ST = Stirlingshire, DK = Dunkeld and BL = Ballater.
Climate data, location and sample information for each collection site. Climate averages from 1981–2010 were taken from the weather station closest to each study site (Met Office, 2016). All samples are numbered according to the Royal Botanic Garden Edinburgh DNA bank numbers and so are preceded by EDNA16-0044XXX. All species were found in sites where more than 20 samples were collected.
ern: corresponding to sequences from Parmelia ernstiae
sax: corresponding to sequences from Parmelia saxatilis s. str.
ser: corresponding to sequences from Parmelia serrana
At each location, 20 specimens were collected from several different woodlands to increase the chance of collecting representative genetic diversity within the area and to decrease the chance of collecting several specimens from one widely dispersed, asexually reproduced individual. Specimens were collected from oak, birch, pine, ash, willow and rocks in each location. Where possible, two specimens with differing positions were collected from each sampled tree and rock. The following were recorded for each specimen: GPS grid reference of the tree or rock, slope aspect, diameter at breast height (DBH), height of lichen from ground and lichen aspect. If collected from a branch then branch height, branch aspect and branch circumference were also recorded. After collection, specimens were dried at room temperature for 4 days, frozen for at least 7 days and then stored at −20 °C. Slope and lichen aspects were converted to folded aspects to account for the circular nature of the measurements, and were zero centred at 225° to reflect the highest heat load corresponding to SW aspects (McCune Reference McCune2007).
For some of the statistical analyses, locations were split into areas with low (750–1050 mm per year) and high (1400–1700 mm per year) average annual rainfall, corresponding to a natural break in epiphyte distributions according to rainfall (Ellis Reference Ellis2018). Glen Creran, Benmore and Loch Earn fall into the high average annual rainfall category and Stirlingshire, Dunkeld and Ballater fall into the low average annual rainfall category (Table 1).
DNA extraction, amplification and restriction digestion
Small sections (c. 4 mm2) of thallus were removed from each specimen and any visible non-target growth was removed to reduce the chance of contamination. Lichens were homogenized using a TissueLyser II® (QIAGEN, Valencia, USA). The tubes were shaken for 1·5 min at 20 Hz, 4–6 times until specimens were converted to a fine powder. DNeasy® Plant Mini Kits from QIAGEN (Valencia, USA) were used for DNA extraction according to the manufacturer's instructions with a few modifications: 420 µl of lysis/digestion buffer was added to each specimen after homogenization; samples were then placed in a thermo mixer for 1 h at 65 °C at 800 rpm; specimens were centrifuged and 110 µl of the uppermost clear liquid taken for use in the extraction process. A QIAxtractor® robot was used for the majority of specimens with an additional 16 individual specimens extracted by hand. A selection of the extractions were run on a 1% agarose gel and stained with SYBR Safe (Life Technologies Corporation, New York, USA) to check for DNA presence and integrity, and visualized with GeneSys image acquisition software. The DNA was deposited, for long-term storage, at the Royal Botanic Garden Edinburgh DNA bank (EDNA) and sample numbers (Table 1) are derived from EDNA numbers and so are preceded by EDNA16-0044XXX.
PCR amplification of the nrITS region was carried out using a Bio-Rad Tetrad® 2 thermal cycler (California, USA) and the fungal specific primers ITS4 (White et al. Reference White, Bruns, Lee, Taylor, Innis, Gelfand, Sninsky and White1990) and ITS1F (Gardes & Bruns Reference Gardes and Bruns1993). Reactions contained 1× TBT-PAR (Samarakoon et al. Reference Samarakoon, Wang and Alford2013), 1× NH4 Buffer (Bioline, London), 2·5 mM MgCl2 (Bioline, London), 0·2 mM dNTPs, 0·4 µM forward and reverse primers, 0·5 U BIOTAQ™ DNA polymerase (Bioline, London) and 1 µl of template DNA, raised to a final volume of 25 µl with autoclaved, distilled H2O. The thermal cycler program used was as follows: 94 °C for 3 min, followed by 30 cycles of 94 °C for 1 min, 55 °C for 1 min and 72 °C for 1·5 min, with a final extension for 5 min at 72 °C.
Restriction digestion has been used in several studies as a method of species delimitation (DePriest Reference DePriest1993; Horton & Bruns Reference Horton and Bruns2001). Based on sequences that were previously collected, it was speculated that P. saxatilis could be distinguished from P. ernstiae and P. serrana because previously analyzed P. saxatilis s. str. sequences consistently contained a Nsp1 cut site within the nrITS2 sequence (GGCATG^C) that was not present in the other species (C. Zuñiga & R. Yahr, unpublished data). Restriction digestion also allowed identification of immature Parmelia sulcata Taylor specimens that were misidentified as a species from the P. saxatilis complex; P. sulcata PCR products are cut twice by Nsp1 revealing three bands after gel electrophoresis. Each digestion reaction contained 17 µl of distilled H2O, 2·5 µl of CutSmart® Buffer (NEB, Massachusetts, USA), 0·38 µl Nsp1 (NEB, Massachusetts, USA) and 5 µl of PCR product. After the digestion mixes were prepared, they were run in a Bio-Rad Tetrad® 2 thermal cycler (California, USA). The thermal cycler program used was 37 °C for 60 min and 65 °C for 20 min. Restriction products were then stained with SYBR Safe (Life Technologies Corporation, New York, USA), run on a 2% agarose gel and visualized with GeneSys image acquisition software.
PCR product purification and sequencing
Fungal barcode nrITS sequencing was carried out on most of the specimens to verify their species identity based on their phylogenetic groupings and sequence similarity to GenBank and RefSeq reference sequences (O'Leary et al. Reference O'Leary, Wright, Brister, Ciufo, Haddad, McVeigh, Rajput, Robbertse, Smith-White and Ako-Adjei2016). Sequencing was also used to verify whether restriction digestion allows reliable identification of P. saxatilis. For PCR product purification, 2 µl of ExoSAP-IT® (Affymetrix, California, USA) was added to 5 µl of PCR product, followed by a thermal cycler programme of 37 °C for 15 min, then 80 °C for 15 min.
The sequencing reaction mix was prepared using the BigDye® Terminator v3.1 cycle sequencing kit (ThermoFisher Scientific, Massachusetts, USA) and the ITS4 and ITS1F primers as used for PCR. The sequencing reaction was performed in the Bio-Rad Tetrad® 2 thermal cycler (California, USA) with the following program: 25 cycles of 95 °C for 30 s, 50 °C for 20 s and 60 °C for 4 min. Sequencing was conducted by Edinburgh Genomics (Edinburgh, UK) on an ABI 3730XL capillary sequencer (Applied Biosystems, California, USA).
Bioinformatics
Sequences were quality checked, trimmed, aligned and edited using Sequencher 5.1 (Gene Codes Corporation, Michigan, USA). They were then aligned using the MAFFT sequence aligner (Katoh & Toh Reference Katoh and Toh2005) using the default settings. Sequences were trimmed again to remove the upstream small subunit and the downstream large subunit, leaving only the nrITS/fungal barcode region.
The nrITS barcodes were grouped into identical unique sequence groups using CD-HIT Suite (Huang et al. Reference Huang, Niu, Gao, Fu and Li2010) to remove identical sequence types. Two maximum likelihood phylogenies were produced on PhyML 3.0 (Guindon et al. Reference Guindon, Dufayard, Lefort, Anisimova, Hordijk and Gascuel2010), one with all sequences and one with only the unique sequence types. The phylogenies were made with default settings except from the following criterion modifications: automatic substitution model selection (beta version), BioNJ starting tree, nearest-neighbour interchange and subtree pruning and regrafting tree improvements. To find well-supported monophyletic groups, 100 bootstrap replicates were used. A sequence of Parmelia adaugescens Nyl. from Japan was used as the outgroup (GenBank ID: AY036992.1). Parmelia ernstiae, P. saxatilis s. str. and P. serrana reference sequences were included in the phylogenetic trees to allow identification of the monophyletic groups. The megaBLAST search function in GenBank was used to check the sequence identity determined using the phylogenetic trees, with default settings (NCBI 2016).
Morphological and chemical coding
Lichen specimens were observed under a dissection microscope in order to code a set of morphological characteristics corresponding to those used for species discrimination in prior studies of the group. The characters coded were: presence or absence and extent of pruinose layer, lobe tip colour, lobe form, shape of lobe tips (rounded or squared), predominant lobe tip roll direction (up or down), presence or absence of lobules, presence or absence of isidia, isidium shape and distribution.
The first two of these variables were ordinal: extent of pruinose layer was categorized as “0” when not present, “+” when found at the lobe tips of three or fewer lobes, “++” when found on all lobe tips and “+++” when found as a thick layer on all lobe tips and also on the main thallus. Lobe tip colour was recorded on a scale of brown shades fading inwards from the margin to the thallus centre in the following order: “very light brown” corresponding to pale or light brown gradient from tip, “light brown” corresponding to dark brown tip outline and light brown gradient from tip, and “brown” corresponding to a dark brown tip outline and dark brown gradient from tip (Supplementary Material Fig. S1, available online). Lobe characteristics were always based on the outermost growing tips of each thallus. The predominant lobe tip roll was obtained by observing the curvature of the lobes within 1 mm of the lobe edges. A thallus was recorded to have lobules if more than 50% of the asexual reproductive structures on the mature region of the thallus were lobules (as opposed to isidia). This prevented thalli with a small number of outgrown isidia from being included in the category of thalli with lobules.
Isidium shape was determined from the older parts of the thallus where isidia are best developed. They were categorized as “cylindrical” if isidia were predominantly longer than twice their diameter or “rounded” if they were predominantly shorter than twice their diameter. Distribution of isidia was determined by observing the distribution of newly forming isidia towards growing tips of lobes. It was categorized as “edge” if new isidia mostly formed around lobe margins, or “throughout” if new isidia mostly formed across the thallus surface and margins.
Major chemical constituents were examined using thin-layer chromatography for all specimens following standard protocols (Orange et al. Reference Orange, James and White2001). Relative Rf values were calculated by comparing with fumarprotocetraric, norstictic and atranorin-containing controls. Main chemical constituents were coded as either present or absent only, or as present, trace and absent. For distinguishing Rf values of fatty acids, a modified version of Solvent B from Arup et al. (Reference Arup, Ekman, Lindblom and Mattsson1993) was used at ratios of 70:36:9.
Morphological hypothesis testing
Hypotheses to be tested originated from prior work describing the segregate species in other geographical regions (Feuerer & Thell Reference Feuerer and Thell2002; Molina et al. Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004; Mattsson et al. Reference Mattsson, Lattman, Divakar and Crespo2013). Morphological or chemical characters were summarized by species (as determined by ITS sequence) and analyzed using chi-squared tests for those characters previously described as helping to discriminate species: lobe tip shape, lobe tip colour, pruinosity, shape of isidia, direction of lobe tip (e.g. rolled up or down), and presence/absence of lobaric, galbinic and fatty acids. Additional characters were also examined: terminal sinus size, distribution of isidia, position on tree (branch versus bole) and distribution with regard to average annual rainfall. In all cases, the null hypothesis was that all three species could not be differentiated based on the character in question. These were tested by tabulating contingency tables for each tested character by species (identified by sequences) and performing chi-squared tests in R (R Development Core Team 2013), and then calculating adjusted residuals following Sharpe (Reference Sharpe2015). All data were tested and a subset of data including only epiphytes was also tested. Critical values were adjusted for the number of contrasts using the Bonferroni correction, and the greater-than-two rule of thumb was applied to determine which cells contributed to the overall chi-squared test (MacDonald & Gardner Reference MacDonald and Gardner2000).
Multivariate analysis
Data were analyzed using classification trees implemented in the ‘rpart’ package (Maindonald & Braun Reference Maindonald and Braun2010; Crawley Reference Crawley2013) for R (R Development Core Team 2013). This methodology is used to classify specimens into groups by splitting their associated data (e.g. morphological characters, ecological measurements) at points which maximize the exclusive representation of a particular group, in our case species identified by barcoding. Splitting occurs in a hierarchical, iterative fashion and group representation is measured using cross-validation to determine an error rate between training and test data. Furthermore, a cost-complexity parameter (cp) is used to measure the balance between complexity (trees with more branches) and fit (error rate), showing how an increase in complexity is balanced by improvement in fit (Zuur et al. Reference Zuur, Ieno and Smith2007). Only epiphyte samples were used in this analysis and samples with more than 50% missing data were excluded. Classification trees were initiated at a high level of complexity (many branches, though up to a minimum number of specimens per node of five) and pruned at a point that minimized the cross-validated error. However, fitting by cross-validation is based on randomized sub-selections of data and the final result can be sensitive to this selection. Trees were therefore fitted 10 000 times and the final solution was determined at the cost-complexity value (cp) that most frequently returned an optimum tree (i.e. minimum cross-validated error).
Trees were fitted in three ways:
1. Using (a) morphological characters, including lobe tip colour and direction of margin roll, pruinosity, thallus colour, first and second sinus sizes, first and second lobe tip shapes, presence of lobules and distribution of isidia, and (b) the status of the metabolites consalazinic, salazinic, protocetraric, fumarprotocetraric, lobaric and galbinic acids, atranorin and the fatty acid compounds lichesterinic and protolichesterinic acids.
2. Using ecological measurements only, including climate (mean maximum and minimum annual temperatures and annual precipitation), substratum (tree species or saxicolous) and, if on a tree, DBH (whether on bole or branch) and height above the ground.
3. Using a combination of morphological characters (excluding chemistry) and ecological measurements.
Results
Across six sites spanning a strong gradient in temperature and rainfall in Scotland, 193 specimens were sampled, 191 were sequenced using ITS barcodes, and 189 were identified using both restriction digestion (RD) and sequencing. Parmelia saxatilis could be confidently identified in all our specimens because each had a group I intron c. 270 bp and each was cut with NspI, leaving two fragments of c. 160 bp and 600 bp. Three specimens could not be sequenced but were cut by NspI, and we infer these are also P. saxatilis s. str. Results of RD were confirmed by sequencing for 65 P. saxatilis specimens. Restriction digestion could not separate P. ernstiae from P. serrana.
As found previously using ITS sequencing, almost all sampled material of the P. saxatilis aggregate could be attributed to one of three clades, P. saxatilis s. str., P. ernstiae or P. serrana (Supplementary Material Fig. S2, available online). Of the molecularly-identified specimens, 83 belonged to the clade containing reference specimens of P. ernstiae and all had intraspecific pairwise uncorrected p-distances of ≤0·008 sequence similarity; 80 specimens belonged to the clade containing the lectotype of P. saxatilis s. str (MAF 6882) with pairwise uncorrected p-distances of ≤0·004, and all 27 P. serrana specimens had 100% similarity to the type (MAF 9756). One specimen (Dk16b, EDNA16-0044220; MH039745) was outside the main sequence group, with a pairwise uncorrected p-distance of 0·020–0·028 from other P. saxatilis specimens in the sample (but with 100% similarity with other specimens from Europe and North America). This specimen was excluded from further testing as its identity did not reliably match reference material. For each taxon, a single sequence type made up the majority of specimens found across the gradient of sites. Of all the specimens, 10 P. saxatilis s. str. were collected from rock, one specimen was collected from a fencepost and all other specimens were collected from trees.
Morphological traits
For all barcoded specimens (n = 190), morphological, chemical and ecological characters were summarized by species (Tables 2, 3 & 4). In all cases, the null hypothesis was that all three species could not be differentiated based on the character in question. The frequencies among character states for tip shape, lobe tip colour, and presence or absence of lobaric, galbinic, lichesterinic and protolichesterinic acids were all significantly different than expected among species, but in no case was any character state or combination of character states exclusive for a given species. The effect size (measured as standardized residuals of chi-squared tests) was very high (over 3) for a subset of features, which we interpret as potentially particularly useful for species hypotheses (Table 5). Similarly, position of thalli on trees (aspect) and relation to average annual rainfall differed quantitatively among species, but without constancy (data not shown). The degree of pruina cover was only marginally significant among species. Within P. ernstiae, differences were only marginally significant between pruina cover and rainfall categories, with significantly more than expected heavily pruinose thalli (+++) in the low rainfall parts of the sample (χ2 = 9·86, P = 0·02) and fewer than expected epruinose thalli, though that effect was not as strong (standardized residual of −1·10 vs 2·57 for heavily pruinose). No differences were found between analysis results of the full barcoded dataset and those 179 samples which were epiphytic only.
Frequency of observation and chi-squared tests of morphological character states among species of the Parmelia saxatilis aggregate identified by ITS barcoding (not including specimen EDNA16-0044220, which fell outside the main sequence group for P. saxatilis s. str.). For each species bold values are significantly higher than expected; italicized values are significantly lower than expected.
*: 0 = none; + = few tips with pruina; ++ = many tips with pruina; +++ = pruina throughout
**: 0 = very light brown or pale; 1 = dark brown at edge only; 2 = brown at edge and fading proximally
Presence/absence (1/0) and chi-squared tests of chemical characters among species of the Parmelia saxatilis aggregate identified by ITS barcoding (not including specimen EDNA16-0044220, which fell outside the main sequence group for P. saxatilis s. str.). Salazinic acid is constant among all specimens. For each species bold values are significantly higher than expected; italicized values are significantly lower than expected. Columns for % represent the percentage of all specimens in which the compound was found.
Frequency and chi-squared tests of ecological characters among species of the Parmelia saxatilis aggregate identified by ITS barcoding (not including specimen EDNA16-0044220, which fell outside the main sequence group for P. saxatilis s. str.). For each species bold values are significantly higher than expected; italicized values are significantly lower than expected.
Characters of species within the Parmelia saxatilis aggregate, with strong effect size in chi-squared analyses (residuals >|2|). Blank cells indicate effect sizes <|2|.
Multivariate analysis
Epiphytic specimens identified by sequencing were classified into groups using classification trees, first investigating the role of morphology and chemistry for accurate species discrimination. Using three character states for chemical characters (presence, trace quantity, absence), the most frequent optimum solution based on these characters occurred over all 10 000 bootstrapped model runs, cp = 1e–11, and the absolute cross-validated error rate was estimated at 12·4% (Fig. 2). In this case, specimens were sorted into species groups based on their lobe tip colour and chemistry only, without discrimination based on other morphological features.
Classification tree using morphological and chemical characters for 161 specimens of the P. saxatilis aggregate sampled across Scotland, excluding cases with missing data. The overall error rate was 12·4%. The terms ern, sax and ser correspond to specimens of Parmelia ernstiae, P. saxatilis s. str. and P. serrana respectively.
The optimal classification tree revealed morphological and chemical variability in terminal groups, although clear trends are apparent. For example, there was relatively high accuracy in the case of P. saxatilis with 2/3 of specimens with the darkest “brown” lobe tip colour, but greater uncertainty in making a positive identification with paler lobe tips, with very high percentages of both P. ernstiae and P. serrana and 1/3 of P. saxatilis specimens. Chemical constituents were helpful in determining these paler specimens, with lobaric acid absent from 85% of P. serrana specimens with pale lobe tips (vs 78% across the entire dataset; 5 with trace lobaric, and 1 with a full-sized spot in TLC) but present in 90% of specimens of P. ernstiae with paler tips.
Using presence/absence information only for chemical data increased the cross-validated error rate (to 16·7%). Following the initial split using lobe tip colour to determine the identification, specimens were further grouped based on presence/absence of lobaric and galbinic acids, terminal lobe tip shape and abundance of pruina (see Supplementary Material Fig. S3, available online). In that case, the terminal groups were less informative, with poor identification accuracy and a higher overall error rate. Using morphology alone, the cross-validated error rate was much higher, at 25·9%, suggesting a poor fit for the model overall.
Analysis based on ecological measurements appeared to be less reliable for specimen determination than the use of morphology and chemistry (data not shown). The optimum solution occurred over 98% of 10 000 bootstrapped model runs, with cp = 0·0175. However, the absolute cross-validated error was estimated at 37·4% using ecological data alone. In the first split, P. saxatilis was identified by a preference for cooler minimum temperatures. In warmer climates, a combination of microhabitat preferences was selected as having the best explanatory power but with very high ecological variability within species.
Finally, a combination of morphological (excluding chemistry) and ecological measurements was tested as a route to the field-based identification of specimens. A single best solution was identified over 91% of 10 000 bootstrapped model runs, with cp = 0·02298 and an absolute cross-validated error of 19·8% (Fig. 3). In terms of identification accuracy, the combination of morphology and ecology improves significantly compared to ecology alone but is less effective than morphology when including chemical characters determined by TLC under laboratory conditions.
Classification tree using morphological and ecological characters for 154 specimens of the Parmelia saxatilis aggregate sampled across Scotland, excluding cases with missing data. The overall error rate was 19·8%. The terms ern, sax and ser correspond to specimens of P. ernstiae, P. saxatilis s. str. and P. serrana respectively, and DBH is the diameter breast height of the bole of the tree and corresponds to tree age.
Discussion
Based on this study, we found that all members of the Parmelia saxatilis aggregate are common across a strong environmental gradient in Scotland, although P. serrana is less abundant than the others. For practical purposes, species in this complex could not be definitively distinguished based on morphology, chemistry, measured ecological preferences or any combination of these characters within the studied area. Despite clear trends which can often discriminate among species, considerable variation exists within all three species and there is overlap of characters and ecological correlates between species. We conclude that in this area, these species are cryptic and can only be confidently distinguished from one another by molecular methods. PCR length assays and restriction digestion analysis of the nrITS region with Nsp1 was found to be a reliable identification tool in this dataset for distinguishing P. saxatilis samples from P. ernstiae and P. serrana; however, for the latter pair of species, comparison to verified reference sequences is the only reliable method of discrimination. Based on classification tree analysis, the most reliable non-molecular dataset for distinguishing species within this geographical area appears to be chemical measurements by TLC, particularly in combination with the pigmentation of lobe tips, though this still had an error rate of over 10% in the current work.
Morphology
Consistent with other morphological studies of cryptic species (Crespo & Lumbsch Reference Crespo and Lumbsch2010), members of the P. saxatilis complex apparently lack diagnostic features. In this as in other studies, morphological characters are very subtle and hard to describe, and categorization can be subjective or “hard to use” (e.g. the degree of shininess on lobe margins; Arup & Berlin Reference Arup and Berlin2011). In this case, the hard-to-use character is the colour of the lobe tip, which is clearly correlated with specimen identity, with a high probability of the darkest tips being P. saxatilis s. str. However, the converse is not true, with c. 1/3 of P. saxatilis s. str. specimens having paler tips, along with most of the P. ernstiae and P. serrana specimens. In these paler-tipped specimens, the best characters for classification are compounds detectable only using thin-layer chromatography. Even using these specialist data, there is still a 12% error rate, where specimens held back from the model for cross-validation are not identified correctly. Using classification trees for prediction of diagnostic characters is necessarily reliant on the dataset included, so a sample of specimens from another geographical location, for example, could provide a different set of discriminating characters. As TLC is unfeasible for the majority of those interested in practical identification, we also tested easily observed ecological characters. This dataset produced a more error-prone model, suggesting morphology and chemistry together offer the best discriminating power. Identification based on morphology alone will be especially challenging if the identifier does not have prior knowledge of the morphological ranges that exist within each species. Patterns recently reported by Ossowska et al. (Reference Ossowska, Guzow-Krzemińska, Dudek, Oset and Kukwa2018) of predominantly laminal isidia and non-overlapping lobes in P. ernstiae, and overlapping lobes with isidia laminal and marginal in P. serrana were not consistently found when we checked our barcoded material, though coding of this character could differ between studies.
We tested the role of environmental plasticity displayed by species in contributing to the high within and low between species diversity by coding environmental factors with each collection made. Although such an explanation has been suggested in other studies of fungi in the Lecanorales (Articus et al. Reference Articus, Mattsson, Tibell, Grube and Wedin2002; Kotelko & Piercey-Normore Reference Kotelko and Piercey-Normore2010), we found no evidence to support this hypothesis. For example, we found little or no relationship between pruinosity level and environmental factors (tested as aspect; P. saxatilis, Kruskal-Wallis chi-squared = 6·6991, df = 3, P = 0·08213; P. ernstiae, Kruskal-Wallis chi-squared = 8·2816, df = 3, P = 0·04054), despite a correlation with water absorption capacity previously suggested by Hauck et al. (Reference Hauck, Jürgens, Brinkmann and Herminghaus2007). More heavily pruinose specimens (P. ernstiae category +++) than expected were found in climatic regions that have a lower average annual rainfall and average annual minimum temperature (χ2 = 9·62, P = 0·001). No such relationship was found for P. saxatilis and the chi-squared test was not conducted to investigate this relationship within P. serrana because of the small sample size. An intriguing result from Hauck et al.'s (Reference Hauck, Jürgens, Brinkmann and Herminghaus2007) study of saxicolous P. saxatilis thalli versus corticolous P. ernstiae showed that P. saxatilis was hydrophobic, associated with lack of pruinosity and a horizontal mode of growth, contrasting with the hydrophilic and pruinose P. ernstiae found on vertical trunks. That study used only a small number of samples, but it would be an interesting avenue to test more explicitly in these variable species.
A larger sample size for P. serrana would be needed to infer which morphological traits are most frequent in this species. The chi-squared analysis showed that (based on the small number of specimens collected) P. serrana is most frequently associated with rounded lobe tips, isidia across the thallus and light brown lobe tips. Supporting these observations, Molina et al. (Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004) described P. serrana as having lobe tips paler than P. saxatilis and a pale greenish grey to whitish grey thallus with apically rounded short lobes that are continuous to overlapping (large second order lobe sinuses) and centrally clustered isidia. In our material, the majority of specimens did have rounded lobes and over half of studied specimens did have centrally concentrated isidia, which contrasts with our sample of P. ernstiae specimens in which the majority of thalli did not have centrally concentrated isidia. It is apparent from examination of the material post-hoc (knowing the species sequence identity) that there appears to be a distinctly different pattern of isidia, including a combination of shape, colour and distribution, among the species and that our isidia distribution category “edge/throughout” mostly seems to capture this variation. However, a separate ecological study of propagule morphology relative to microsite characteristics might be informative. Propagule morphology in other genera has been shown to respond to environmental factors (e.g. Pseudocyphellaria; K. Bogomazova, unpublished data).
In other studies, P. ernstiae has been distinguished from P. saxatilis based on its thick pruinose layer, deep cut lobes and mixture of lobules and isidia, and lack of lobaric acid (Feuerer & Thell Reference Feuerer and Thell2002; Molina et al. Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004; Thell et al. Reference Thell, Elix, Feuerer, Hansen, Kärnefelt, Schüler and Westberg2008). While some P. ernstiae specimens which conform to this description were collected in the present work, the majority did not express all of these character states. For example, a thick layer of pruinosity (category +++) was not found to be significantly associated with P. ernstiae specimens, and these most frequently had no pruinose layer (44·8%). Lobules were found only in six P. ernstiae specimens and they were also found in four P. saxatilis samples, as noted by Feuerer & Thell (Reference Feuerer and Thell2002). Parmelia serrana is recognized by Molina et al. (Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004) and Mattson et al. (Reference Mattsson, Lattman, Divakar and Crespo2013) by having shiny, non-pruinose lobes, something also found in this study (92% lacked pruina). Taken together, this suggests that the distinguishing morphological traits for P. ernstiae that have been indicated in the literature are not useful for identification of Scottish material.
Molina et al. (Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004) identified P. saxatilis by narrow, non-overlapping lobes (with large, open sinuses). The majority of P. saxatilis samples in this study did have wide and open sinuses but we found a large variation, and many P. saxatilis specimens had lobes which did overlap. Molina et al. (Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004) also distinguished P. saxatilis as being “generally darker in colour, especially at the lobe margins” when compared to P. serrana. This conforms with the findings from our study, since P. serrana most frequently had light brown or very light brown lobe tips and P. saxatilis was most frequently darker with more extensively pigmented lobe tips.
The secondary chemistries of species in this aggregate vary among reports (Molina et al. Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004; Thell et al. Reference Thell, Elix, Feuerer, Hansen, Kärnefelt, Schüler and Westberg2008, Reference Thell, Tsurykau, Persson, Hansson, Åsegård, Kärnefelt and Seaward2017; Ossowska et al. Reference Ossowska, Guzow-Krzemińska, Dudek, Oset and Kukwa2018). For example, we found protocetraric acid only in P. saxatilis, whereas Molina et al. (Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004) report it in both P. saxatilis and P. serrana. Chemical analysis using TLC and HPLC for a small number of specimens by Thell et al. (Reference Thell, Elix, Feuerer, Hansen, Kärnefelt, Schüler and Westberg2008, Reference Thell, Tsurykau, Persson, Hansson, Åsegård, Kärnefelt and Seaward2017) found lobaric acid and a suite of fatty acids present in both P. serrana and P. ernstiae. In Scotland, there do appear to be quantitative differences in frequency of these compounds. Lobaric acid was found only infrequently, in six of 27 specimens of P. serrana in this study, corresponding with the more recently analyzed material from the Nordic countries (Thell et al. Reference Thell, Tsurykau, Persson, Hansson, Åsegård, Kärnefelt and Seaward2017). Also, in contrast to Ossowska et al. (Reference Ossowska, Guzow-Krzemińska, Dudek, Oset and Kukwa2018), 10–20% of specimens of P. ernstiae and around half of our examined specimens of P. serrana in the present study have a suite of fatty acids including those normally referred to as lichesterinic and protolichesterinic (TLC using multiple solvent systems and specimens reported to have these compounds shows that there is likely to be some confusion in the naming of these spots; P. Harrold, unpublished data).
Niche preferences
In this study, P. ernstiae and P. serrana were found to have a preference for wetter and warmer regions (in the hyperoceanic west), whereas P. saxatilis was found more frequently in the more continental, cooler and drier sampling sites. In contrast, Molina et al. (Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004) state that from their collections it appears as though P. serrana prefers more continental climatic conditions and P. saxatilis prefers more oceanic climatic conditions, though the actual climate data were not specified.
Parmelia serrana was found on branches much more frequently than the other two species, with 44% of specimens being found on tree branches rather than the main bole. A similar preference has been found in other lichen species; for example, Melanelixia subaurifera is more frequently found on twigs and branches than M. fuliginosa (Louwhoff & Esslinger Reference Louwhoff, Esslinger, Smith, Aptroot, Coppins, Fletcher, Gilbert, James and Wolseley2009). Very few specimens were collected from branches for this study because the branches were often too high to access. This might explain why so few P. serrana specimens were found. In future studies, tree branches should be targeted in specimen collection to investigate this relationship and discover whether P. serrana is more abundant than this study suggests. It is possible that the more frequent occurrence on branches may ameliorate the overall wetter Scottish conditions compared with more continental material.
The three lichen species were not found to have any preference for height on bole (though only up to 2 m was sampled, apart from branches), substratum tree DBH (suggesting that they have no preference for tree age), tree species, or aspect on the bole. Molina et al. (Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004) found that P. serrana and P. saxatilis are similarly distributed on deciduous and coniferous trees, and mossy rocks, although P. saxatilis is found on rocks much more frequently. In our study, both P. saxatilis and P. ernstiae were found on rocks but P. serrana was not. However, we sampled only up to two saxicolous specimens per site. Several authors report that P. ernstiae is infrequently found on rocks (Feuerer & Thell Reference Feuerer and Thell2002; Molina et al. Reference Molina, Crespo, Blanco, Lumbsch and Hawksworth2004; Thell et al. Reference Thell, Elix, Feuerer, Hansen, Kärnefelt, Schüler and Westberg2008).
Conclusions
Parmelia ernstiae, P. saxatilis s. str. and P. serrana populations are present in Scotland and cover a wide range of macroclimates. Inspection of morphological characters of these populations revealed that the three species are indeed cryptic and not consistently diagnosable by any combination of studied characters. High morphological diversity exists within all three species, which is partly attributable to environmental determination of morphological traits. For example, P. ernstiae was found to have a thicker pruinose layer in climatic regions with higher rainfall, supporting the idea that calcium oxalate crystals have a physiological function and that production of this thick pruinose layer is an adaptive response to the environment. We recommend that future studies quantify incident light by measuring canopy conditions at each collection and angle of substratum to further investigate the degree of hydrophobicity and pruina production.
By testing previous morphological and chemical hypotheses to distinguish a set of closely related species, we found that we could not consistently and accurately discriminate all specimens to sequence-determined species level using traditional characters, or even sets of characters, in multivariate analyses. However, from the data considered, it is evident that there are clear morphological, chemical and ecological trends that can be brought to bear to increase the likelihood of accurate identification consistent with ITS barcodes. In the present dataset, the combination of morphology and chemical analysis by TLC could discriminate among over 85% of specimens. In contrast, the measured ecological and microhabitat characteristics are less useful for species discrimination, resulting in a greater than 1/3 error rate. Using a combination of morphological and ecological observations, this error rate improves to c. 1/5.
Andrei Tsurykau checked our determinations of fatty acids, Chris Ellis contributed valuable discussions about analysis and assisting with R code, and Antje Ahrends provided the code for geographical coordinate conversion. All are warmly thanked for their help.
Supplementary Material
For supplementary material accompanying this paper visit https://doi.org/10.1017/S0024282919000069.
