Phylogenetic insight into the Lecidea atrobrunnea complex – evidence of narrow geographic endemics and the pressing need for integrative taxonomic revisions

Abstract

Species of lichen-forming fungi (LFF) display an array of geographical distribution patterns. Among the broadly distributed lichen-forming fungal species, the degree of reproductive isolation and genetic substructure among populations varies widely, in some cases masking unrecognized diversity or meaningful biogeographical patterns. Lecidea atrobrunnea (Raymond ex Lam. & DC.) Schaer. s. lat. (Lecideaceae) is a widespread species complex that has been studied for over two centuries since its initial description. The diversity of the L. atrobrunnea group is highest in western North America, where a dizzying array of morphologies and chemistry can occur at local scales. Here we investigate whether the assumed cosmopolitan distribution of L. atrobrunnea s. lat. is an artifact of taxonomic limitations and masks biogeographical patterns in this species complex. To address these questions, we compiled sequence data from the standard fungal barcoding marker (ITS) for over 100 specimens within this complex, in addition to genome-scale data from a subset of these representing over 1600 single-copy nuclear genes spanning over 3 Mb of the genome. Our study corroborates the perspective that the morphologically and chemically variable Lecidea atrobrunnea group reflects a complex of distinct species-level lineages, with 42–83 candidate species inferred from the ITS region and high levels of diversity inferred from a subset of specimens using genome-scale data. However, both phenotype- and molecular-based species boundaries remained unsettled, with the most common nominal taxa recovered as highly polyphyletic and with conflict among different molecular species delimitation approaches. Our study also highlights the potential for geographically restricted species, with fascinating biogeographical patterns, challenging, in part, the assumed cosmopolitan distribution of L. atrobrunnea s. lat. This study provides valuable direction for future research that will be crucial in understanding diversification and establishing a robust taxonomy for this well-known species complex.

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. 2018; Moncada et al. 2018; Simon et al. 2018) to others with disjunct populations in similar habitats separated by incredible geographical distances (Culberson 1972; Bailey & James 1979; Lindblom & Søchting 2008) to some species that occur in diverse habitats across multiple continents (Werth 2011; Onuţ-Brännström et al. 2017). 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 1972). Among the LFF species studied with molecular sequence data, some show little evidence of strong population subdivision, even among intercontinental populations (Piercey-Normore 2006; Geml et al. 2010; Park et al. 2012; Garrido-Benavent et al. 2018). 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. 2014; Leavitt et al. 2018; Zhang et al. 2022). Rather than broad, intercontinental distributions, many LFF species, particularly epiphytes, have geographical ranges comparable to those of vascular plants (Culberson 1972). Factors driving varying levels of reproductive isolation among populations continue to be explored (Printzen & Ekman 2003; Printzen et al. 2013; Jorna et al. 2021; Werth et al. 2021), 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 2003; Garrido-Benavent et al. 2018), and these have been a general focus of lichenological research (Geml et al. 2010; Fernández-Mendoza et al. 2011; Leavitt et al. 2015; Onuţ-Brännström et al. 2017). 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. 2018; Magain et al. 2018; Lutsak et al. 2020).

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 1805). 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. 2020). This lichen is both morphologically and chemically polymorphic (Leuckert & Hertel 2003). Thallus morphology varies from highly dispersed to contiguous, with colours ranging from pale yellowish brown to olivaceous green to dark brown (Leuckert & Hertel 2003). 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. 2020), 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 2004; McCune et al. 2017). 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 2003). 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. (2020) 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 2003), where a dizzying array of morphologies and chemistry can occur at local scales (McCune 2017) (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 2003; McCune 2017). 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 2003; Hertel & Printzen 2004). 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 2004). Other less frequently encountered taxa include L. deplanaica (Hertel & Leuckert) McCune with deplanaic acid (McCune et al. 2017) and L. truckeei Herre with schizopeltic acid (Herre 1911). 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 2010; Lücking et al. 2021). 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 2011; Schneider et al. 2016). 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. 2011; Pino-Bodas et al. 2013; Boluda et al. 2019; Spjut et al. 2020). 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. 2017).

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. 2015; Zhang et al. 2022), 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. (2020), 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. (2020)); L. confluens and L. swartzioidea were used as outgroups (Ruprecht et al. 2020).

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 1972; Orange et al. 2001). Specimen identification followed McCune (2017) 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. 2012), 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. 2004; Lücking et al. 2021). For new specimens, we extracted total genomic DNA using the Wizard Genomic DNA Purification Kit (Promega, USA). Amplification and sequencing followed Hale et al. (2019). 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. 2012). 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. 2022). 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. (2020)) 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 2008; Rozewicki et al. 2017), 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. 2020), 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. (2020)) 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. 2021). 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. 2013). A maximum likelihood (ML) tree was reconstructed from the ITS MSA with IQ-TREE 2 (Minh et al. 2020), and the tree file was adjusted into an ultrametric tree using the ‘chronos’ function from the R package ape (Paradis et al. 2004). 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. 2015) 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. 2012), 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. 2021) 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. 2016), 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 2008; Rozewicki et al. 2017), implementing the default parameters. Ambiguously aligned regions in the individual BUSCO MSAs were removed using Gblocks (Talavera & Castresana 2007), implementing the default parameters. ML phylogenetic reconstructions of individual MSAs were inferred using IQ-TREE 2 (Minh et al. 2020), 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 2015). We also used the quartet-based species delimitation method SODA (Species bOundry Delimitation using ASTRAL) v. 1.0.1 (Rabiee & Mirarab 2021), implemented with ASTRAL using the default parameters. Subsequently, we used FASconCAT-G v. 1.05 (Kück & Longo 2014) to concatenate all BUSCO alignments into a single supermatrix (Tonini et al. 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. 2012), and contigs > 4000 bps were retained as targets for the HybPiper v. 1.2 pipeline (Johnson et al. 2016) as described above. Topologies were inferred from individual mitochondrial contigs using IQ-TREE 2 (Minh et al. 2020), with 1000 ultrafast bootstrap replicates. Alignments representing individual mitochondrial contigs were also concatenated using FASconCAT-G v. 1.05 (Kück & Longo 2014) and phylogenetic relationships were inferred using IQ-TREE 2 (Minh et al. 2020) and ASTRAL v. 5.7.8 (Mirarab & Warnow 2015). We used ASAP (Puillandre et al. 2021) 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. (2020)) 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. (2020) 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. (2020), 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. (2020) 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 2004; McCune 2017), 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. 2022).

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. 2012; Lücking et al. 2020). 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. 2013; Spjut et al. 2020), 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 2022). 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 2017). 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 2017). 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. 2016; Boluda et al. 2019; LaGreca et al. 2020) and others showing secondary metabolites coinciding with distinct lineages (Schmitt & Lumbsch 2004; Fehrer et al. 2008).

Specimen misidentification may also play a role in the interpretation of our results (Hertel & Printzen 2004; McCune 2017). 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. 2020). 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. 2020). 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. 2022) and L. deplanaica (McCune et al. 2017), 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. 2015; Zhang et al. 2022). 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. 2022). Other LFF species complexes with a mixture of geographically restricted species and others with broad distributions have been observed (Leavitt et al. 2013, 2018; Divakar et al. 2019; Zhang et al. 2022). 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. 2010; Fernández-Mendoza & Printzen 2013; Widhelm et al. 2019), isolation by distance (Walser et al. 2005; Widhelm et al. 2021), niche differentiation (Guttová et al. 2019; Allen et al. 2021; Moncada et al. 2021) 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. 2018; Jorna et al. 2021; Widhelm et al. 2021; Keuler et al. 2022), and similar approaches will be crucial for understanding diversification and phylogeography in this group.