Introduction
Species of lichen-forming fungi (LFF) display a broad array of geographical distribution patterns and population substructure, ranging from narrow endemics (Allen et al. Reference Allen, McKenzie, Sleith and Alter2018; Moncada et al. Reference Moncada, Mercado-Dízz and Lücking2018; Simon et al. Reference Simon, Goffinet, Magain and Sérusiaux2018) to others with disjunct populations in similar habitats separated by incredible geographical distances (Culberson Reference Culberson1972; Bailey & James Reference Bailey and James1979; Lindblom & Søchting Reference Lindblom and Søchting2008) to some species that occur in diverse habitats across multiple continents (Werth Reference Werth and Fontaneto2011; Onuţ-Brännström et al. Reference Onuţ-Brännström, Tibell and Johannesson2017). Among the broadly distributed LFF, the degree of reproductive isolation and genetic substructure among populations also varies widely. Traditionally, non-vascular cryptogams, including LFF, were thought to be easily disseminated and without well-defined geographical ranges (Culberson Reference Culberson1972). Among the LFF species studied with molecular sequence data, some show little evidence of strong population subdivision, even among intercontinental populations (Piercey-Normore Reference Piercey-Normore2006; Geml et al. Reference Geml, Kauff, Brochmann and Taylor2010; Park et al. Reference Park, Jeong and Hong2012; Garrido-Benavent et al. Reference Garrido-Benavent, de los Ríos, Fernández-Mendoza and Pérez-Ortega2018). In other cases, genetic data have revealed unexpected population subdivisions and hidden species diversity masked within many assumed cosmopolitan or widespread lichens (Lücking et al. Reference Lücking, Dal-Forno, Sikaroodi, Gillevet, Bungartz, Moncada, Yánez-Ayabaca, Chaves, Coca and Lawrey2014; Leavitt et al. Reference Leavitt, Westberg, Nelsen, Elix, Timdal, Sohrabi, St. Clair, Williams, Wedin and Lumbsch2018; Zhang et al. Reference Zhang, Clancy, Jensen, McMullin, Wang and Leavitt2022). Rather than broad, intercontinental distributions, many LFF species, particularly epiphytes, have geographical ranges comparable to those of vascular plants (Culberson Reference Culberson1972). Factors driving varying levels of reproductive isolation among populations continue to be explored (Printzen & Ekman Reference Printzen and Ekman2003; Printzen et al. Reference Printzen, Domaschke, Fernández-Mendoza and Pérez-Ortega2013; Jorna et al. Reference Jorna, Linde, Searle, Jackson, Nielsen, Nate, Saxton, Grewe, Herrera-Campos and Spjut2021; Werth et al. Reference Werth, Meidl and Scheidegger2021), but recognizing that the distributions of lichens, vascular plants and other organismal groups are influenced by the same broad ecological and historical factors is inescapable.
High altitude/latitude habitats often support species of LFF that occur in similar habitats across the globe (Weber Reference Weber2003; Garrido-Benavent et al. Reference Garrido-Benavent, de los Ríos, Fernández-Mendoza and Pérez-Ortega2018), and these have been a general focus of lichenological research (Geml et al. Reference Geml, Kauff, Brochmann and Taylor2010; Fernández-Mendoza et al. Reference Fernández-Mendoza, Domaschke, García, Jordan, Martin and Printzen2011; Leavitt et al. Reference Leavitt, Divakar, Ohmura, Wang, Esslinger and Lumbsch2015; Onuţ-Brännström et al. Reference Onuţ-Brännström, Tibell and Johannesson2017). Exploring the evolutionary histories and genetic diversity of widespread LFF provides insight into the processes and factors that have shaped the diversity and distribution of species over time and space, in addition to helping resolve and stabilize taxonomy (Grewe et al. Reference Grewe, Lagostina, Wu, Printzen and Lumbsch2018; Magain et al. Reference Magain, Tniong, Goward, Niu, Goffinet, Sérusiaux, Vitikainen, Lutzoni and Miadlikowska2018; Lutsak et al. Reference Lutsak, Fernández-Mendoza, Kirika, Wondafrash and Printzen2020).
Lecidea atrobrunnea (Raymond ex Lam. & DC.) Schaer. s. lat. (Lecideaceae) has been studied for over two centuries since its initial description in 1805 (de Lamarck & de Candolle Reference de Lamarck and de Candolle1805). Since that time, it has been documented on all continents in montane and alpine belts of mountains, in addition to climatically similar habitats in subpolar and polar regions (Ruprecht et al. Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020). This lichen is both morphologically and chemically polymorphic (Leuckert & Hertel Reference Leuckert and Hertel2003). Thallus morphology varies from highly dispersed to contiguous, with colours ranging from pale yellowish brown to olivaceous green to dark brown (Leuckert & Hertel Reference Leuckert and Hertel2003). While the taxonomy of lecideoid lichens is typically based on a limited number of microscopic traits, such as spore size and septation or ascus-type (see Ruprecht et al. Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020), the taxonomy of the L. atrobrunnea group has largely been based on morphological and chemical characters.
Secondary metabolite variation has led, in part, to the description of a number of taxa at both species and infraspecific ranks (Hertel & Printzen Reference Hertel, Printzen, Nash, Ryan, Diederich, Gries and Bungartz2004; McCune et al. Reference McCune, Curtis and Di2017). Depsides of the orcinol-type with long aliphatic side chains, depsidones of the β-orcinol type and dibenzofurans are the typical substances used in the taxonomy of the L. atrobrunnea complex (Leuckert & Hertel Reference Leuckert and Hertel2003). However, the taxonomy of the group remains unsettled. To proceed towards a comprehensive taxonomic revision of the L. atrobrunnea complex, integrating inferences from molecular sequence data with phenotypic traits will be critical for a robust and stable taxonomy. Recently, Ruprecht et al. (Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020) provided evidence that the L. atrobrunnea group comprises a monophyletic clade within a well-supported lineage of lecideoid lichen (‘L01’) that also includes L. confluens, L. promiscens Nyl., L. swartzioidea Nyl. and an undescribed species. However, taxonomic sampling for members of the L. atrobrunnea group remains relatively sparse.
The diversity of the L. atrobrunnea group is highest in western North America (Leuckert & Hertel Reference Leuckert and Hertel2003), where a dizzying array of morphologies and chemistry can occur at local scales (McCune Reference McCune2017) (Fig. 1). Furthermore, this complex ranks among the most commonly encountered lichens in montane and alpine habitats in the western USA. Over 15 different chemotypes are recognized within this species complex in North America (Leuckert & Hertel Reference Leuckert and Hertel2003; McCune Reference McCune2017). Based on current taxonomy, L. atrobrunnea s. str. is characterized by the production of confluentic acid and/or 2ʹ-O-methylperlatolic acid, although acid-deficient and other chemical variations also occur (Leuckert & Hertel Reference Leuckert and Hertel2003; Hertel & Printzen Reference Hertel, Printzen, Nash, Ryan, Diederich, Gries and Bungartz2004). Two morphologically similar taxa occurring predominantly in alpine habitats are L. protabacina Nyl., characterized by the production of the stictic acid syndrome (=L. atrobrunnea subsp. stictica Hertel & Leuckert), and L. syncarpa Zahlbr., characterized by the production of the norstictic acid syndrome (=L. atrobrunnea subsp. saxosa Hertel & Leuckert). Lecidea perlatolica Hertel & Leuckert is less common and is distinguished by the production of perlatolic acid (Hertel & Printzen Reference Hertel, Printzen, Nash, Ryan, Diederich, Gries and Bungartz2004). Other less frequently encountered taxa include L. deplanaica (Hertel & Leuckert) McCune with deplanaic acid (McCune et al. Reference McCune, Curtis and Di2017) and L. truckeei Herre with schizopeltic acid (Herre Reference Herre1911). Challenges in interpreting the morphological and chemical variation has led to the contemporary practical treatment of members of the L. atrobrunnea complex as either an unusually variable species or a confusing complex of taxa.
Figure 1. Morphological variation in the Lecidea atrobrunnea complex from the La Sal Mountain Range in Utah, USA. A, Lecidea aff. promiscens (Leavitt 18565 [BRY-C]). B, L. atrobrunnea s. lat. (Leavitt 21360 [BRY-C]). C, ‘L. atrobrunnea s. lat. clade 1, 3.1’ (Leavitt 18299 [BRY-C]). D, ‘L. atrobrunnea s. lat. clade 1, 3.2’ (Leavitt 18646 [BRY-C]). E, ‘L. atrobrunnea s. lat. clade 1, 4.1’ (Leavitt 21175 [BRY-C]). F, ‘L. atrobrunnea s. lat. clade 3’ (Leavitt 18658 [BRY-C]). G, ‘L. protabacina 1’ (Leavitt 18650 [BRY-C]). H, ‘L. protabacina 2’ (Leavitt 18653 [BRY-C]). I, L. syncarpa (Leavitt 18488 [BRY-C]). J, L. syncarpa (Leavitt 18462 [BRY-C]). K, L. syncarpa (Leavitt 18582 [BRY-C]). L, L. aff. syncarpa (Leavitt 18642 [BRY-C]). Scales = 2 mm. In colour online.
Over recent decades, the use of DNA sequence data has foundationally transformed our understanding of lichen biogeography and taxonomy (Printzen Reference Printzen2010; Lücking et al. Reference Lücking, Leavitt and Hawksworth2021). Prior to the widespread availability of molecular techniques, lichen classification was primarily based on morphological, anatomical and chemical characters, which may not consistently circumscribe natural groups (Lumbsch & Leavitt Reference Lumbsch and Leavitt2011; Schneider et al. Reference Schneider, Resl and Spribille2016). However, the use of genetic data is not a panacea for resolving taxonomy and has frequently resulted in unsettled or contentious taxonomic conundrums (Leavitt et al. Reference Leavitt, Johnson, Goward and St. Clair2011; Pino-Bodas et al. Reference Pino-Bodas, Martín, Burgaz and Lumbsch2013; Boluda et al. Reference Boluda, Rico, Divakar, Nadyeina, Myllys, McMullin, Zamora, Scheidegger and Hawksworth2019; Spjut et al. Reference Spjut, Simon, Guissard, Magain and Sérusiaux2020). Despite the potential for introducing additional complexities into lichen taxonomy and our associated understanding of biogeographical patterns, molecular systematics are the foundation for advancing our understanding of Earth's fungal diversity (Spatafora et al. Reference Spatafora, Aime, Grigoriev, Martin, Stajich, Blackwell, Heitman, Howlett, Crous, Stukenbrock, James and Gow2017).
To add to the two centuries of study into the L. atrobrunnea s. lat. group, here we use DNA sequence data to assess species diversity within this complex. By focusing on recent collections made from western North America, we aim to test whether the observed morphological and chemical variation reflects a complex of jointly occurring similar taxa or simply a phenotypically polymorphic species. Is the assumed cosmopolitan distribution of L. atrobrunnea s. lat. an artifact of taxonomic limitations (Singh et al. Reference Singh, Dal Grande, Divakar, Otte, Leavitt, Szczepanska, Crespo, Rico, Aptroot and Cáceres2015; Zhang et al. Reference Zhang, Clancy, Jensen, McMullin, Wang and Leavitt2022), and what insight into the processes and factors that have shaped the diversity and distribution can be gained from utilizing molecular sequence data? To address these questions, we compiled sequence data from the standard fungal barcoding marker for over 100 specimens within this complex, in addition to genome-scale data from a subset of these. Our results provide a crucial perspective into the messy, unsettled taxonomy and biogeography of this well-known lichen, highlighting that well-designed phylogenetic studies will be required for any future taxonomic revision.
Materials and Methods
Taxon sampling
Sampling was focused on members of the Lecidea atrobrunnea complex in western North America. We attempted to compile currently available DNA sequence data, augmented by our own recent sampling in western North America (n = 120; see Supplementary Material File S1, available online). For this study, sampling included L. atrobrunnea s. lat (n = 49, total), including the subspecies L. atrobrunnea subsp. atrobrunnea (19) and L. atrobrunnea s. lat. (30; GenBank identifications or otherwise ambiguous determinations), L. ‘fuscoatra’ (L.) Ach. (2), L. glacierensis A. Abbas & R. Mamut (4), L. perlatolica (1), L. promiscens (16), L. protabacina (10), L. swartzioidea (1), L. syncarpa (24), Lecidea ‘sp. 1’ sensu Ruprecht et al. (Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020), and six unidentified specimens represented by sequences downloaded from GenBank. Samples originated from Antarctica (n = 6), Argentina (18), Austria (3), China (10), Greenland (Denmark; 1), Norway (1), Turkey (3) and the USA (78). Exploratory phylogenetic analyses were used to corroborate the placement of specimens within the Lecidea atrobrunnea complex and close relatives (clade ‘L01’ sensu Ruprecht et al. (Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020)); L. confluens and L. swartzioidea were used as outgroups (Ruprecht et al. Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020).
For specimens housed at the Herbarium of Non-Vascular Cryptogams, Brigham Young University, Provo, Utah, USA (BRY-C), morphological characters were assessed using an Olympus SZH dissecting microscope. Anatomical observations were made on material mounted in water with an Olympus BH-2 microscope. Secondary metabolites were identified using thin-layer chromatography (TLC), following standard methods with solvent system ‘G’ (Culberson Reference Culberson1972; Orange et al. Reference Orange, James and White2001). Specimen identification followed McCune (Reference McCune2017) using the TLC and morphological data generated here. For specimens represented by sequences downloaded from GenBank and not seen by us, determinations on GenBank were provided in quotation marks.
DNA extraction and sequencing
We compiled sequence data representing the standard fungal DNA barcode for fungi, the internal transcribed spacer region (ITS; Schoch et al. Reference Schoch, Seifert, Huhndorf, Robert, Spouge, Levesque, Chen and Fungal2012), for all specimens included here. The ITS has previously been shown to have high discriminatory power among species in Lecanoraceae and is a powerful tool for a ‘first pass’ assessment of species-level diversity (Grube et al. Reference Grube, Baloch and Arup2004; Lücking et al. Reference Lücking, Leavitt and Hawksworth2021). For new specimens, we extracted total genomic DNA using the Wizard Genomic DNA Purification Kit (Promega, USA). Amplification and sequencing followed Hale et al. (Reference Hale, Fisher, Keuler, Smith and Leavitt2019). We performed amplifications using Cytiva PuReTaq Ready-To-Go™ PCR Beads (Fisher Scientific, USA) following the manufacturer's protocol. Complementary strands were sequenced with the same primers used for amplifications, and sequencing reactions were performed using BigDye 3.1 (Applied Biosystems, Foster City, CA, USA). The PCR products were run on an ABI 3730 automated sequencer machine (Applied Biosystems) at the DNA Sequencing Center at Brigham Young University (Provo, Utah, USA).
Single marker approaches may be inadequate for delimiting species boundaries, particularly among closely related species (Dupuis et al. Reference Dupuis, Roe and Sperling2012). To explore the consistency between inferences made from the single marker ITS dataset and phylogenomic datasets, 27 specimens representing diversity in the L. atrobrunnea complex were selected for short-read shotgun sequencing for subsequent genome skimming (Zhang et al. Reference Zhang, Clancy, Jensen, McMullin, Wang and Leavitt2022). Specimens for Illumina sequencing were chosen based on genetic diversity initially observed from samples available at BRY-C (exclusively collected from western North America). For these specimens, total genomic DNA was extracted from lichen thalli, using the E.Z.N.A. Plant DNA DS Mini Kit (Omega Bio-Tek, Inc., USA) following the manufacturer's protocol. We prepared total genomic DNA following the standard Illumina whole genome sequencing (WGS) library preparation process with Adaptive Focused Acoustics for shearing (Covaris), followed by an AMPure cleanup process. The DNA was then processed with the NEBNext Ultra™ II End Repair/dA-Tailing Module end-repair, together with the NEBNext Ultra™ II Ligation Module (New England BioLabs), while using standard Illumina index primers. Libraries were pooled and sequenced using the HiSeq 2500 sequencer in high output mode, by the DNA Sequencing Center at Brigham Young University, USA, with 125 cycle paired-end (PE) reads.
Candidate species delimitation using the standard marker for fungal barcoding
We initially evaluated the monophyly of the L. atrobrunnea + L. promiscens aggregate (clade ‘L01’ sensu Ruprecht et al. (Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020)) based on the standard fungal DNA barcoding marker (ITS). An initial multiple sequence alignment (MSA) was inferred from a genus-wide sampling of ITS sequences using the program MAFFT v. 7 (Katoh & Toh Reference Katoh and Toh2008; Rozewicki et al. Reference Rozewicki, Yamada and Katoh2017), implementing the G-INS-i alignment algorithm and ‘1PAM/K = 2’ scoring matrix, with an offset value of 0.9, and the remaining parameters set to default values. To assess the monophyly of the L. atrobrunnea aggregate, we inferred a maximum likelihood (ML) ITS topology from the initial MSA, using IQ-TREE 2 (Minh et al. Reference Minh, Schmidt, Chernomor, Schrempf, Woodhams, von Haeseler and Lanfear2020), with 1000 ultrafast bootstrap replicates. A second MSA was inferred exclusively from sequences recovered in the clade representing the L. atrobrunnea + L. promiscens aggregate inferred from the broader MSA.
The ITS MSA representing the L. atrobrunnea + L. promiscens aggregate (clade ‘L01’ sensu Ruprecht et al. (Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020)) was subsequently used for two different sequence-based species delimitation analyses. For the genetic distance approach, the ITS MSA was analyzed using the species delimitation program ASAP (Assemble Species by Automatic Partitioning; Puillandre et al. Reference Puillandre, Brouillet and Achaz2021). ASAP uses pairwise genetic distances from a single-locus alignment to delimit candidate species partitions, and here we used the ASAP web server (https://bioinfo.mnhn.fr/abi/public/asap/) with a Jukes-Cantor (JC69) model. For the tree-based species delimitation analysis, we employed bPTP, a Bayesian implementation of the Poisson Tree Processes model (Zhang et al. Reference Zhang, Kapli, Pavlidis and Stamatakis2013). A maximum likelihood (ML) tree was reconstructed from the ITS MSA with IQ-TREE 2 (Minh et al. Reference Minh, Schmidt, Chernomor, Schrempf, Woodhams, von Haeseler and Lanfear2020), and the tree file was adjusted into an ultrametric tree using the ‘chronos’ function from the R package ape (Paradis et al. Reference Paradis, Claude and Strimmer2004). The ultrametric ML tree was analyzed to delimit the number and grouping of candidate species using the bPTP web server (https://species.h-its.org/ptp/), implementing 500 000 MCMC generations and a burn-in proportion of 0.1.
Assembling phylogenomic datasets
We prepared two genome-scale datasets for phylogenomic analyses, one comprising the Benchmarking Universal Single-Copy Orthologs (BUSCO) from the nuclear genome (Simão et al. Reference Simão, Waterhouse, Ioannidis, Kriventseva and Zdobnov2015) and the second representing the majority of the mitochondrial genome. BUSCO genes were extracted from a draft assembly of L. atrobrunnea s. lat. assembled for this study. We selected L. atrobrunnea specimen ‘sl18283’ (BRY-C) for the de novo draft assembly. Paired-end reads were assembled using SPAdes v. 3.14.1 (Bankevich et al. Reference Bankevich, Nurk, Antipov, Gurevich, Dvorkin, Kulikov, Lesin, Nikolenko, Pham and Prjibelski2012), with default parameters and assigned k-mer values at 21, 33, 55 and 77. Scaffolds from the SPAdes assembly were assessed to identify single-copy nuclear genes for downstream phylogenomic reconstructions using BUSCO v. 5.2.2 (Manni et al. Reference Manni, Berkeley, Seppey, Simão and Zdobnov2021) and the ‘ascomycota_odb10’ dataset for comparison.
BUSCO markers passing the quality filters were used as the targets in the HybPiper v. 1.2 pipeline (Johnson et al. Reference Johnson, Gardner, Liu, Medina, Goffinet, Shaw, Zerega and Wickett2016), implementing the ‘reads_first.py’ function. We drew the HybPiper results with the ‘get_seq_lengths.py’ function and visualized the coverage heatmap using R v. 4.1.2 (R Development Core Team 2012), with the ‘geom_raster function’ from the ggplot2 package. Genes with less than 50% coverage were excluded. The remaining contigs from all BUSCO genes across all samples were assembled using the ‘retrieve_sequences.py’ function, and MSAs were generated for individual BUSCO regions with MAFFT v. 7 (Katoh & Toh Reference Katoh and Toh2008; Rozewicki et al. Reference Rozewicki, Yamada and Katoh2017), implementing the default parameters. Ambiguously aligned regions in the individual BUSCO MSAs were removed using Gblocks (Talavera & Castresana Reference Talavera and Castresana2007), implementing the default parameters. ML phylogenetic reconstructions of individual MSAs were inferred using IQ-TREE 2 (Minh et al. Reference Minh, Schmidt, Chernomor, Schrempf, Woodhams, von Haeseler and Lanfear2020), with 1000 ultrafast bootstrap replications. The resulting individual gene trees were analyzed to estimate a species tree using the coalescent-based summary method ASTRAL v. 5.7.8 (Mirarab & Warnow Reference Mirarab and Warnow2015). We also used the quartet-based species delimitation method SODA (Species bOundry Delimitation using ASTRAL) v. 1.0.1 (Rabiee & Mirarab Reference Rabiee and Mirarab2021), implemented with ASTRAL using the default parameters. Subsequently, we used FASconCAT-G v. 1.05 (Kück & Longo Reference Kück and Longo2014) to concatenate all BUSCO alignments into a single supermatrix (Tonini et al. Reference Tonini, Moore, Stern, Shcheglovitova and Ortí2015). A maximum likelihood topology was inferred from this supermatrix using IQ-TREE 2, with 1000 ultrafast bootstrap replicates.
In addition to the single-copy nuclear genome data matrix, we also assembled mitochondrial phylogenomic datasets. Mitochondrial contigs were identified and extracted from a draft SPAdes assembly of L. atrobrunnea s. lat. generated for this study using a custom BLAST search. Mitochondrial contigs were assessed to investigate relative read coverage using the iterative Geneious read mapper (Kearse et al. Reference Kearse, Moir, Wilson, Stones-Havas, Cheung, Sturrock, Buxton, Cooper, Markowitz and Duran2012), and contigs > 4000 bps were retained as targets for the HybPiper v. 1.2 pipeline (Johnson et al. Reference Johnson, Gardner, Liu, Medina, Goffinet, Shaw, Zerega and Wickett2016) as described above. Topologies were inferred from individual mitochondrial contigs using IQ-TREE 2 (Minh et al. Reference Minh, Schmidt, Chernomor, Schrempf, Woodhams, von Haeseler and Lanfear2020), with 1000 ultrafast bootstrap replicates. Alignments representing individual mitochondrial contigs were also concatenated using FASconCAT-G v. 1.05 (Kück & Longo Reference Kück and Longo2014) and phylogenetic relationships were inferred using IQ-TREE 2 (Minh et al. Reference Minh, Schmidt, Chernomor, Schrempf, Woodhams, von Haeseler and Lanfear2020) and ASTRAL v. 5.7.8 (Mirarab & Warnow Reference Mirarab and Warnow2015). We used ASAP (Puillandre et al. Reference Puillandre, Brouillet and Achaz2021) to delimit candidate species from the mitochondrial alignments.
Except web server services, all analyses were performed at the supercomputational facility of the Faculty of Science, Kasetsart University (SciKU HPC) in Bangkok, Thailand.
Results
New ITS sequences generated for this study were deposited in GenBank under Accession nos. OR180026–OR180048; the ITS alignment of the L. atrobrunnea + L. promiscens aggregate (clade ‘L01’ sensu Ruprecht et al. (Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020)) spanned 554 aligned nucleotide position characters (Supplementary Material File S2, available online). Short-read data are available in the NCBI Sequence Read Archive under Accession PRJNA951751.
Phylogenetic circumscription of the L. atrobrunnea species complex
The ML topology inferred from the ITS alignment representing the ‘L01’ clade sensu Ruprecht et al. (Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020) recovered two main lineages: the ‘L. atrobrunnea clade’ and the ‘L. promiscens clade’ (Fig. 2; Supplementary Material File S3, available online). The L. atrobrunnea clade included specimens representing L. atrobrunnea, L. ‘fuscoatra’ (two specimens from Turkey, probably misidentifications), L. glacierensis, L. perlatolica, L. protabacina, L. aff. promiscens and L. syncarpa, and was recovered with 100% bootstrap support (BS). The L. promiscens clade included a number of specimens identified as L. promiscens, Lecidea ‘sp. 1’ sensu Ruprecht et al. (Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020), and a single specimen representing L. swartzioidea; this clade was recovered as sister to the L. atrobrunnea group with strong support.
Figure 2. Maximum likelihood topology of the Lecidea atrobrunnea and L. promiscens clades, clade ‘L01’ sensu Ruprecht et al. (Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020) inferred from ITS sequence data (see Supplementary Material File S3, available online, for the full tree without collapsed clades). Each clade includes a provisional name based on taxonomic identification, detected secondary metabolites are shown in bold (‘CON’ = confluentic acid, ‘MP’ = 2ʹ-O-methylperlatolic acid, ‘NOR’ = norstictic acid, ‘PER’ = perlatolic acid, ‘STI’ = stictic acid and ‘NA’ = not available; a ‘/’ indicates separate chemotypes found in distinct specimens recovered within the clade), and in parentheses, the geographical origin of specimens recovered within each clade, with the number of individual samples shown in brackets. The four subclades of L. atrobrunnea clade 1 are indicated with grey boxes with subclade numbers in quotes, and the major clade designation is shown in bold and in the tip label.
Our analyses of the ITS data revealed a high degree of phylogenetic substructure within both major clades inferred here (Fig. 2). In many cases, well-supported clades coincided with distinct geographical regions. Most samples within the L. promiscens clade originated from southern South American and Antarctica, although specimens from the Colorado Plateau and Rocky Mountain regions of the western USA were also recovered in two separate lineages within this clade. However, specimens identified as L. promiscens were not recovered as monophyletic, and several specimens collected from western USA were recovered within the L. atrobrunnea clade and provisionally called here the ‘L. aff. promiscens clade’ (Fig. 2).
In the L. atrobrunnea clade, specimens identified as L. atrobrunnea were not recovered as monophyletic and were recovered in four separate lineages, in addition to two clades representing L. protabacina (Fig. 2). Here, the four major L. atrobrunnea clades were arbitrarily called ‘L. atrobrunnea clade 1’, ‘L. atrobrunnea clade 2’, ‘L. atrobrunnea clade 3’ and ‘L. atrobrunnea clade 4’. Within ‘L. atrobrunnea clade 1’, four well-supported subclades were recovered. Specimens within ‘L. atrobrunnea clade 1’ were chemically polymorphic, and this clade also included specimens identified as L. perlatolica and L. syncarpa (Supplementary Material File S1, available online). Due to the limited thallus size of some specimens, secondary metabolites were not tested for specimens recovered in ‘L. atrobrunnea clade 2’. We did not detect 2ʹ-O-methylperlatolic or confluentic acids, secondary metabolites typical for L. atrobrunnea, in the single specimen representing ‘L. atrobrunnea clade 3’. Rather, this specimen contained several unidentified secondary metabolites. ‘L. atrobrunnea clade 4’ was represented exclusively by sequences downloaded from GenBank and secondary metabolite information was not available.
Specimens identified as L. protabacina (producing the stictic acid syndrome) were recovered in two separate clades, one represented exclusively by specimens collected from the Colorado Plateau (‘L. protabacina 1’), and the second clade comprised of a single specimen from Austria and another from the Colorado Plateau (‘L. protabacina 2’) (Fig. 2). ‘L. protabacina 1’ was recovered as sister to a recently described species, L. glacierensis, from the Tianshan Mountains, Xinjiang Province, China.
Most specimens representing L. syncarpa were recovered in a well-supported clade sister to L. ‘fuscoatra’, that were probably misidentifications (Fig. 2). Specimens recovered within the two major subclades in the ‘L. syncarpa complex’ produced either norstictic acid and accessory acids or had no detectable secondary metabolites (Fig. 2). We did not detect 2ʹ-O-methylperlatolic or confluentic acids in the norstictic acid-producing specimens. As previously noted, several norstictic acid-producing specimens were also recovered within ‘L. atrobrunnea clade 1’ (Fig. 2).
Candidate species delimitation using the standard fungal DNA barcode
From the 120 ITS sequences representing the ‘L01’ clade (L. atrobrunnea clade + L. promiscens clade), the ASAP species delimitation analysis provided the strongest support for 89 partitions, for example, ‘candidate species’ (40–89 partitions inferred in the ten top-scoring partitions; Supplementary Material File S4, available online). The tree-based species delimitation method, bPTP, delimited more than 100 species and was not considered further due to the assumed excessive splitting. Candidate species inferred using ASAP were subsequently considered within a phylogenetic context (see Fig. 2 and Supplementary Material File S3). In cases where ASAP species partitions were inferred to be closely related in the ITS topology and originated from the same geographical location, these were collapsed, and ultimately 37 species hypotheses were considered here (Fig. 2).
A total of 23 species hypotheses were circumscribed representing L. atrobrunnea s. lat.: 17 within ‘L. atrobrunnea clade 1’; two within ‘L. atrobrunnea clade 2’; one within ‘L. atrobrunnea clade 3’; one within ‘L. atrobrunnea clade 4’; and two representing L. protabacina (=L. atrobrunnea subsp. stictica). Generally, each species hypothesis was comprised of samples from the same geographical region, with only two exceptions: ‘L. atrobrunnea 1 2.2’ was represented by a specimen from China (GenBank) and another from southern Utah, USA; ‘L. atrobrunnea 1 4.1’ was represented by a specimen from southern Nevada and another from southern Utah, USA (Fig. 2). Candidate species ‘L. atrobrunnea 1 1.5’ represented the only sampled specimen identified as L. perlatolica.
The two specimens from Turkey, L. ‘fuscoatra’ (probably misidentifications), comprised a single species partition in the ASAP analysis. Similarly, specimens representing L. glacierensis (all from China) were delimited as a single species partition.
Of the two candidate species within the ‘L. syncarpa complex’, one was comprised of two specimens from the La Sal Mountains in southern Utah, USA (no secondary metabolites were detected in these) and the second of representative sequences from the Colorado Plateau and Great Basin (USA), in addition to a single specimen from China (Fig. 2). A third candidate species, ‘L. aff. syncarpa’, was represented by a single specimen from the La Sal Mountains in southern Utah, USA.
Specimens representing L. promiscens corresponded to six candidate species, five of which included specimens collected from southern Argentina or Antarctica (Fig. 2). Candidate species ‘L. promiscens agg. 1.5’ also included specimens from China and western USA, in addition to a specimen from southern Argentina. Candidate species ‘L. promiscens agg. 1.4’ was collected from alpine habitat in the Absaroka Range in southern Montana, USA. Furthermore, the candidate species ‘L. aff. promiscens’ (nested within the ‘L. atrobrunnea complex’) comprised one specimen from southern Arizona and two from the La Sal Mountains in southern Utah, USA.
The regionally focused sampling in the La Sal Mountains, Utah, USA (an alpine sky island on the Colorado Plateau) revealed high levels of diversity, with 14 of the 38 species hypotheses represented in this limited geographical area.
Inferences from genome-scale data
From the draft genome assembly of L. atrobrunnea s. lat. (sl18283 [BRY-C]; Supplementary Material File S5, available online), 1663 complete and single-copy BUSCO genes were recovered from a total of 1706 BUSCO genes searched, 97.5% of all BUSCO genes (Supplementary Material File S6, available online). Two additional BUSCO genes were removed that had less than 50% coverage across the 25 metagenomic samples (Supplementary Material File S7, available online). Attempts to generate Illumina sequencing libraries for representatives of the ‘L. protabacina 2’ and ‘L. aff. syncarpa’ clades failed. Alignments of 1661 BUSCO genes were included in the subsequent phylogenomic inferences (Supplementary Material File S8, available online). Concatenating the 1661 BUSCO alignments resulted in a supermatrix alignment spanning 3 211 186 bps. Concatenating the individual BUSCO alignments with ambiguously aligned regions removed resulted in a supermatrix spanning 3 180 693 bps (Supplementary Material File S9, available online).
Both ML and ASTRAL analyses of the 25 specimens selected for high-throughput sequencing inferred a consistent, well-resolved phylogeny (Fig. 3A). In both the nuclear and mitochondrial phylogenies, ‘L. atrobrunnea clade 1’, ‘L. atrobrunnea clade 2’, ‘L. atrobrunnea clade 3’, ‘L. protabacina 1’, ‘L. aff. promiscens’, ‘L. promiscens agg. 1.5’ and ‘L. syncarpa complex’ were recovered as monophyletic, although relationships among these clades differed among phylogenomic datasets (Fig. 3). The SODA species delimitation analysis of the 1661 individual gene trees resulted in 16 candidate species (Fig. 3A). The SODA species delimitations mostly overlapped with the ASAP partitions based on the ITS alignment.
Figure 3. Maximum likelihood phylogenomic inferences from 1661 single-copy nuclear genes (A) and a 74.9 Kb alignment from the mitochondrial genome (B). Bootstrap values < 100% are shown on branches with dashed lines displaying support values. Species delimited from 1661 single-copy nuclear gene trees using SODA are shown in the panel in ‘A’; species partitions inferred from five different mitochondrial fragments using ASAP are shown in the panel in ‘B’, each column representing results inferred from each individual mitochondrial fragments alignment. In the tip labels, numbers in bold indicate the major subclade designation inferred from phylogenetic analysis of ITS sequence data, followed by lineage codes (numbers with decimal points) and individual specimen number in brackets. In colour online.
The ML analysis of the concatenated mitochondrial data (spanning 74 906 bps) representing the 25 samples selected for high-throughput sequencing is reported in Fig. 3B. Most major clades inferred from the nuclear phylogeny were consistently recovered in the concatenated mitochondrial topology, although ‘L. protabacina 1’ was recovered as polyphyletic. The ASAP analyses of the multiple sequence alignments of the five mitochondrial fragments > 4000 bps resulted in two to nine candidate species, with variable specimen assignments based on the mitochondrial fragment analyzed (Fig. 3B).
Discussion
Our study corroborates the perspective that the morphologically and chemically variable Lecidea atrobrunnea group reflects a complex of distinct species-level lineages. This inference is supported by both single-marker candidate species delimitation analyses (Supplementary Material File S4, available online) and species delimitation analyses based on over 1600 single-copy nuclear loci (Fig. 3). However, both phenotype- and molecular-based species boundaries remained unsettled, due, in part, to limitations with current sampling. The most commonly occurring taxa in montane habitats in western North America, L. atrobrunnea, L. protabacina and L. syncarpa (Hertel & Printzen Reference Hertel, Printzen, Nash, Ryan, Diederich, Gries and Bungartz2004; McCune Reference McCune2017), are not monophyletic, with representatives of each taxon occurring in multiple, distinct clades (Fig. 2). Furthermore, we see evidence for narrow geographical distributions for many candidate species-level lineages, although current sampling is woefully incomplete. Below we discuss the taxonomic and evolutionary implications of our study, providing a crucial perspective for future research into this widespread species complex.
Our results suggest that species-level diversity within the L. atrobrunnea complex may be vastly underestimated, with both single marker and phylogenomic datasets delimiting unexpectedly high levels of species-level diversity (Fig. 3, Supplementary Material File S4). However, any interpretation of species boundaries and the total number of species in this complex would be preliminary given the incompleteness of the current sampling. Intensive sampling across taxonomic and geographical scales will be required before a formal taxonomic revision is implemented, and the impact of careful taxonomic sampling, including specimens collected throughout the candidate species’ distributions, on phylogenetic and taxonomic inferences can only be speculated. We predict that species-level lineages within the L. atrobrunnea clade will be found to comprise many narrow endemics and a more limited number of truly widespread species. Similar patterns have also been suggested for another montane/alpine lichen that commonly co-occurs with members of the L. atrobrunnea clade, Lecanora polytropa s. lat. (Zhang et al. Reference Zhang, Clancy, Jensen, McMullin, Wang and Leavitt2022).
Broader taxonomic and specimen sampling, coupled with genetic data and quantitative phenotypic data, will be required to robustly delimit species boundaries within an integrative framework (Fujita et al. Reference Fujita, Leaché, Burbrink, McGuire and Moritz2012; Lücking et al. Reference Lücking, Aime, Robbertse, Miller, Ariyawansa, Aoki, Cardinali, Crous, Druzhinina and Geiser2020). Inferences from the standard DNA fungal barcode (nrITS) will probably continue to provide a valuable perspective that can be integrated with geographical distributions, secondary metabolite variation, and other morphological characters to begin thorough revisions for this group. We highlight that conflicting molecular-based species circumscriptions are commonplace in empirical studies (Carstens et al. Reference Carstens, Pelletier, Reid and Satler2013; Spjut et al. Reference Spjut, Simon, Guissard, Magain and Sérusiaux2020), and careful biological interpretation of the results is critical for meaningful species delimitations. Best practices in the taxonomic revision of this clade should include, as far as possible: 1) representation of all species worldwide from the L. atrobrunnea complex; 2) representation of variation across the entire geographical ranges of species and potential species-level lineages; 3) examination and comparison of type specimens for all the included named taxa in order to apply names to the appropriate clades (Williams Reference Williams, Monro and Mayo2022). This is a tall order for this common, cosmopolitan species complex with an equally complex taxonomic history.
Secondary metabolites play a central role in the current taxonomy of the Lecidea atrobrunnea complex (McCune Reference McCune2017). However, chemical identification can be challenging. While thin-layer chromatography is often necessary for specimen determination, accurate metabolite identification is complicated by close R f values for a number of compounds, secondary metabolites in low concentrations, and/or scant medulla for testing (McCune Reference McCune2017). Our results suggest that, while secondary metabolites will probably continue to play a crucial role in the taxonomy of this complex, differences in secondary metabolite variation on their own will not be adequate to diagnose distinct species. For example, specimens producing the norstictic acid syndrome are found in multiple clades in both ‘L. atrobrunnea clade 1’, the ‘L. syncarpa complex’ and ‘L. aff. syncarpa’ (Fig. 2). These specimens typically have paler brown to medium brown thalli (Fig. 1L). However, these characters alone are insufficient to determine which of the ‘L. syncarpa’ clades the specimen belongs to. Similarly, specimens representing L. protabacina, producing the stictic acid syndrome, were recovered in two distinct clades (Fig. 2). In other cases, distinct DNA-based candidate species occasionally comprised specimens with different chemistries, for example L. atrobrunnea s. lat. 1 ‘2.1’, L. atrobrunnea 1 ‘2.3’, and L. atrobrunnea s. lat. 1 ‘4.1’ (Fig. 2). The mixed utility of secondary metabolites as diagnostic characters has been demonstrated in other lichens, with some cases showing chemically polymorphic species (Mark et al. Reference Mark, Saag, Leavitt, Will-Wolf, Nelsen, Tõrra, Saag, Randlane and Lumbsch2016; Boluda et al. Reference Boluda, Rico, Divakar, Nadyeina, Myllys, McMullin, Zamora, Scheidegger and Hawksworth2019; LaGreca et al. Reference LaGreca, Lumbsch, Kukwa, Wei, Han, Moon, Kashiwadani, Aptroot and Leavitt2020) and others showing secondary metabolites coinciding with distinct lineages (Schmitt & Lumbsch Reference Schmitt and Lumbsch2004; Fehrer et al. Reference Fehrer, Slavíková-Bayerová and Orange2008).
Specimen misidentification may also play a role in the interpretation of our results (Hertel & Printzen Reference Hertel, Printzen, Nash, Ryan, Diederich, Gries and Bungartz2004; McCune Reference McCune2017). For example, two GenBank accessions identified as L. fuscoatra (HQ605927, HQ605930) are probably misidentifications, since L. fuscoatra falls outside of the clade ‘L01’, which includes the L. atrobrunnea group (Ruprecht et al. Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020). In the present study, we found specimens identified as L. promiscens occurring in two distantly related clades, including one within the L. atrobrunnea complex (Figs 1A & 2). Future investigations into these specimens might reveal subtle differences that could lead to different determinations. Similarly, we relied here on data available from GenBank and the associated taxonomic identifications. GenBank sequences are prone to error due to misidentification and/or the limitations imposed by the lack of a polished species-level taxonomy for many groups, including LFF (Pentinsaari et al. Reference Pentinsaari, Ratnasingham, Miller and Hebert2020). While taxonomic revisions and reinterpretations of diagnostic taxonomic characters are beyond the scope of the present study, future work on the taxonomy of the L. atrobrunnea complex will rely on accessible vouchered specimens and a thorough reconsideration of phenotypic variation. Broader taxonomic sampling, including sampling multiple distinct populations representing different taxa, will be crucial to revising taxonomy in the group. Ultimately, exploring a broad range of potentially diagnostic phenotypic characters, including secondary metabolite variation, may help to create a stable taxonomy for this complex. For example, recent work has resulted in the description of new species within this complex, L. glacierensis (Mamut et al. Reference Mamut, Jiamahat and Abbas2022) and L. deplanaica (McCune et al. Reference McCune, Curtis and Di2017), based on combinations of distinctive chemistry, morphology, distribution and ecology.
Our study highlights the potential for geographically restricted species, with fascinating biogeographical patterns, challenging, in part, the assumed cosmopolitan distribution of L. atrobrunnea s. lat. (Singh et al. Reference Singh, Dal Grande, Divakar, Otte, Leavitt, Szczepanska, Crespo, Rico, Aptroot and Cáceres2015; Zhang et al. Reference Zhang, Clancy, Jensen, McMullin, Wang and Leavitt2022). Given the limitations with specimen sampling in the present study (sparse sampling in western North America, with extremely limited sampling on other continents), any interpretation of biogeographical patterns must be regarded as tentative. However, our results provide tantalizing evidence that some species-level lineages may have very narrow geographical distributions. Sampled L. atrobrunnea s. lat. populations in distinct mountain ranges in western North America rarely shared candidate species, and different mountain ranges often harboured multiple, distinct species-level lineages (Supplementary Material File S1, available online). Phylogeographic substructure was inferred from variation in ITS sequence data (Fig. 2), and it is likely that generating ITS sequences from specimens collected from landscape-level sampling will provide a critical insight into the spatial distributions of candidate species and limitations to dispersal (Zhang et al. Reference Zhang, Clancy, Jensen, McMullin, Wang and Leavitt2022). Other LFF species complexes with a mixture of geographically restricted species and others with broad distributions have been observed (Leavitt et al. Reference Leavitt, Fernández-Mendoza, Pérez-Ortega, Sohrabi, Divakar, Vondrák, Lumbsch and St. Clair2013, Reference Leavitt, Westberg, Nelsen, Elix, Timdal, Sohrabi, St. Clair, Williams, Wedin and Lumbsch2018; Divakar et al. Reference Divakar, Wei, McCune, Cubas, Boluda, Leavitt, Crespo, Tchabanenko and Lumbsch2019; Zhang et al. Reference Zhang, Clancy, Jensen, McMullin, Wang and Leavitt2022). Similarly, we predict that some species-level lineages inferred here will have narrow geographical distributions and others will have broad, intercontinental distributions.
Understanding the processes and factors that have shaped the diversity and distribution of the L. atrobrunnea complex remains enigmatic. How historical processes (Otálora et al. Reference Otálora, Martínez, Aragón and Molina2010; Fernández-Mendoza & Printzen Reference Fernández-Mendoza and Printzen2013; Widhelm et al. Reference Widhelm, Grewe, Huang, Mercado-Díaz, Goffinet, Lücking, Moncada, Mason-Gamer and Lumbsch2019), isolation by distance (Walser et al. Reference Walser, Holderegger, Gugerli, Hoebee and Scheidegger2005; Widhelm et al. Reference Widhelm, Grewe, Huang, Ramanauskas, Mason-Gamer and Lumbsch2021), niche differentiation (Guttová et al. Reference Guttová, Fačkovcová, Martellos, Paoli, Munzi, Pittao and Ongaro2019; Allen et al. Reference Allen, McMullin, Wiersma and Scheidegger2021; Moncada et al. Reference Moncada, Mercado-Díaz, Magain, Hodkinson, Smith, Bungartz, Pérez-Pérez, Gumboski, Sérusiaux and Lumbsch2021) and other factors interacted to give rise to the species diversity in the L. atrobrunnea complex remains to be explored. A robust, comprehensive reconstruction of the evolutionary history of the L. atrobrunnea complex is lacking, and our study highlights the potential of genome-scale datasets for resolving relationships within this complex (Fig. 3). A robust phylogeny and species delimitations in the L. atrobrunnea complex will provide the framework for investigating character evolution and phylogeography. The power of genome-scale approaches for understanding the process of speciation has recently been highlighted (e.g. Allen et al. Reference Allen, McKenzie, Sleith and Alter2018; Jorna et al. Reference Jorna, Linde, Searle, Jackson, Nielsen, Nate, Saxton, Grewe, Herrera-Campos and Spjut2021; Widhelm et al. Reference Widhelm, Grewe, Huang, Ramanauskas, Mason-Gamer and Lumbsch2021; Keuler et al. Reference Keuler, Jensen, Barcena-Peña, Grewe, Lumbsch, Huang and Leavitt2022), and similar approaches will be crucial for understanding diversification and phylogeography in this group.
Acknowledgements
We dedicate this study to Prof. Pier Luigi Nimis, whose research, including lichen diversity in montane and alpine habitats, has been transformative. This research was funded by the Canyonlands Natural History Association (Moab, UT, USA; Discovery Award 20-02-USFS), the United States Geological Survey (G19AC00400 and G19AC00269), and the Department of Biology, Brigham Young University. We thank Nastassja Noell, Barb Smith, Larry St. Clair, Griffin Leavitt and Mike Felix for invaluable help with fieldwork. The funders had no role in the design of the study, in the collection, analyses, or interpretation of data, in the writing of the manuscript, or in the decision to publish the results.
Author ORCIDs
Nopparat Anantaprayoon, 0000-0003-0834-4513; Jason Hollinger, 0000-0003-2465-2487; Ekaphan Kraichak, 0000-0002-8437-2180; Steven D. Leavitt, 0000-0002-5034-9724.
Supplementary Material
The Supplementary Material for this article can be found at https://doi.org/10.1017/S0024282923000270.
Supplementary Material File S1. List of specimens included in this study, including the NCBI GenBank Accession Nos.
Supplementary Material File S2. The multiple sequence alignment of the internal transcribed spacer (ITS) region from members of the Lecidea atrobrunnea and L. promiscens clades – clade ‘L01’ sensu Ruprecht et al. (Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020).
Supplementary Material File S3. Maximum likelihood topology of the Lecidea atrobrunnea and L. promiscens clades – clade “L01” sensu Ruprecht et al. (Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020) inferred from ITS sequence data – see Fig. 2 for the tree with collapsed clades. Specimen numbers can be matched to the list of specimens in supplementary file S1.
Supplementary Material File S4. Comparison of sample assignments to candidate species across the ten best ASAP partitions for the Lecidea atrobrunnea and L. promiscens clades – clade “L01” sensu Ruprecht et al. (Reference Ruprecht, Fernández-Mendoza, Türk and Fryday2020) inferred from a multiple sequence alignment of the internal transcribed spacer region (ITS, standard fungal barcode marker) and comprising 113 sequences.
Supplementary Material File S5. Draft metagenomic assembly of L. atrobrunnea s. lat. (Leavitt 18283 [BRY-C]). The de novo assembly was performed using SPAdes v3.14.1, and the resulting scaffolds represent metagenomic sequences, including the predominant lichen-forming fungi L. atrobrunnea s. lat.
Supplementary Material File S6. Benchmarking Universal Single-Copy Orthologs (BUSCO) from the nuclear genome L. atrobrunnea s. lat. (Leavitt 18283 [BRY-C]).
Supplementary Material File S7. Heatmaps showing coverage across Benchmarking Universal Single-Copy Orthologs (BUSCO) used for the phylogenomic analyses of members of the Lecidea atrobrunnea complex.
Supplementary Material File S8. Individual gene alignments of 1,661 Benchmarking Universal Single-Copy Orthologs (BUSCO) markers of members of the Lecidea atrobrunnea complex.
Supplementary Material File S9. The concatenated 1,661 Benchmarking Universal Single-Copy Orthologs (BUSCO) gene alignments filtered using GBlocks.
