Population genomic analyses of RAD sequences resolves the phylogenetic relationship of the lichen-forming fungal species Usnea antarctica and Usnea aurantiacoatra
expand article infoFelix Grewe, Elisa Lagostina§, Huini Wu|, Christian Printzen§, H. Thorsten Lumbsch
‡ Field Museum of Natural History, Chicago, United States of America
§ Senckenberg Research Institute and Natural History Museum Frankfurt, Frankfurt, Germany
| Loyola University Chicago, Chicago, United States of America
Open Access


Neuropogonoid species in the lichen-forming fungal genus Usnea exhibit great morphological variation that can be misleading for delimitation of species. We specifically focused on the species delimitation of two closely-related, predominantly Antarctic species differing in the reproductive mode and representing a so-called species pair: the asexual U. antarctica and the sexual U. aurantiacoatra. Previous studies have revealed contradicting results. While multi-locus studies based on DNA sequence data provided evidence that these two taxa might be conspecific, microsatellite data suggested they represent distinct lineages. By using RADseq, we generated thousands of homologous markers to build a robust phylogeny of the two species. Furthermore, we successfully implemented these data in fine-scale population genomic analyses such as DAPC and fineRADstructure. Both Usnea species are readily delimited in phylogenetic inferences and, therefore, the hypothesis that both species are conspecific was rejected. Population genomic analyses also strongly confirmed separated genomes and, additionally, showed different levels of co-ancestry and substructure within each species. Lower co-ancestry in the asexual U. antarctica than in the sexual U. aurantiacoatra may be derived from a wider distributional range of the former species. Our results demonstrate the utility of this RADseq method in tracing population dynamics of lichens in future analyses.


Antarctica, Ascomycota, lichens, Parmeliaceae, phylogeny, RADseq


Over the last decades, the use of DNA sequence data to delimit species and reconstruct phylogenetic relationships has become standard (Barraclough and Nee 2001; de Queiroz 2007; Holder and Lewis 2003; Huelsenbeck et al. 2001; Taylor et al. 2000; Wiens and Penkrot 2002). In groups with high morphological plasticity and homoplasy in phenotypical data sets, such as fungi, molecular data have dramatically changed our understanding of evolution and coinciding taxonomic interpretations (Hibbett et al. 2007; James et al. 2006; Lutzoni et al. 2004; McLaughlin et al. 2009; Robbertse et al. 2006; Schoch et al. 2009; Spatafora et al. 2017; Spatafora and Robbertse 2010; Stajich et al. 2009).

The general lineage species concept (de Queiroz 2007) allows researchers to use different empirical data to test the hypothesis of lineage separation, including phenotypical characters and molecular data. The latter dataset often provides strong evidence if analysed within a rigorous statistical framework (Rannala 2015). With regards to species delimitation, numerous studies of lichen-forming fungi detected distinct lineages lacking obvious distinguishing phenotypical characters, the so-called cryptic species (Bickford et al. 2007; Crespo and Lumbsch 2010; Crespo and Pérez-Ortega 2009; Lumbsch and Leavitt 2011). However, some studies also demonstrated that morphologically distinct populations could not be separated using single- or multi-locus genetic data. These results have been interpreted either as an indication of recent diversification and incomplete lineage sorting (Leavitt et al. 2016a; Zhao et al. 2017) or that the phenotypes represented populations of the same species (Articus et al. 2002; Buschbom and Mueller 2006; Kotelko and Piercey-Normore 2010; Lohtander et al. 1998; Myllys et al. 2001; Velmala et al. 2009). The latter result was often found in so-called species pairs. These are lichens that differ in forming either ascomata and reproducing sexually or forming asexual diaspores (soredia), which propagate the fungal and photosynthetic partner simultaneously (Mattsson and Lumbsch 1989; Poelt 1970; Tehler 1982). Otherwise, these species are morphologically identical, but were traditionally regarded as distinct species due to their different reproductive modes (Poelt 1972). The Parmotrema perforatum group was used as a model system of species delimitation based on the reproductive mode and secondary metabolites (Culberson and Culberson 1973). However, a recent study suggests that the phylogenetic relationships between sexual and asexual populations might be more complex (Widhelm et al. 2016).

We here focus on a complex of two morphologically similar species that differ in their reproductive mode and are considered a species pair in the genus Usnea: U. aurantiacoatra and U. antarctica, the latter reproducing by asexual soredia (Walker 1985). Within the genus, there is group of species predominantly occurring in Antarctica and adjacent cool-temperate to polar regions with a thallus that consists of yellow (containing usnic acid) and blackish areas (caused by melanins). Species in this group, which is also called neuropogonoid (Lumbsch and Wirtz 2011; Wirtz et al. 2008), can be difficult to distinguish by their general appearance and hence, molecular data, such as DNA marker sequencing, can be helpful in delimiting lineages (Articus 2004; Lumbsch and Wirtz 2011; Seymour et al. 2007; Wirtz et al. 2008; Wirtz et al. 2012). Earlier studies based on morphological and chemical data considered the neuropogonoid species as a subgenus Neuropogon in Usnea (Lamb 1964; Walker 1985) or as a distinct genus Neuropogon (Krog 1976; 1982; Lamb 1939). Molecular studies confirmed Usnea (including Neuropogon) as a monophyletic genus within Parmeliaceae (Crespo et al. 2007); however, the relationship of Neuropogon and Usnea remained ambiguous. A two-marker DNA analysis of Usnea elevated Neuropogon to a generic rank (Articus 2004), but subsequent studies provided evidence that Neuropogon is polyphyletic with a core group nested within Usnea (Seymour et al. 2007; Wirtz et al. 2008; Wirtz et al. 2012; Wirtz et al. 2006). Multi-locus DNA sequence data could not delimit individuals of the species U. antarctica and U. aurantiacoatra (Seymour et al. 2007; Wirtz et al. 2012) suggesting that they might be conspecific. In contrast, a recent microsatellite study provided evidence that the two species represent isolated lineages (Lagostina et al. 2018). Given the contradicting results of multi-locus and microsatellite data, we decided to employ a reduced genomic dataset to revisit the species delimitation of U. antarctica and U. aurantiacoatra.

The advent of next-generation sequencing (NGS), also referred to as high-throughput sequencing, drastically changed the scale of molecular datasets for systematic analyses and revolutionised our ability to assess evolutionary histories of organisms (Kraus and Wink 2015; Wachi et al. 2018; Zimmer and Wen 2015). Many molecular studies, such as the former species delimitation efforts for neuropogonoid Usnea spp., were limited to, at most, a dozen markers because their production would require tedious lab work and costly Sanger-sequencing (Hoffman and Lendemer 2018; Wilkinson et al. 2017). Population genomics of closely related organisms often relied on the descriptive power of microsatellite markers (Hodel et al. 2016). Compared to these traditional lab methods, NGS techniques allow a relatively straight-forward production of genome-scale datasets. Direct sequencing NGS methods, such as de-novo genome sequencing (Ellegren 2014), re-sequencing (Stratton 2008) or RNAseq of expressed genes (Ozsolak and Milos 2011; Wickett et al. 2014), can provide whole genome-scale data but may still be limited by high sequencing costs. Therefore, these methods are rarely applied to population studies which require the sequencing of many individuals. However, often a subset of the genome entails sufficient polymorphisms to answer questions of phylogenetic or population genomic studies. Hence, many NGS methods for systematic analyses are designed to be economical and generate reduced genome representation datasets (Allendorf 2017; Davey et al. 2011). One of these reduced representation methods is genotype-by-sequencing and its altered approach, which is known as restriction associated DNA sequencing (RADseq) (Baird et al. 2008). We recently designed a RADseq approach for metagenomic data derived from symbiotic lichen genomes, which allows reduced representation genomic analyses of numerous individuals for population-scale studies (Grewe et al. 2017).

By using genome-wide single nucleotide polymorphisms (SNPs) produced by restriction site-associated DNA sequencing (RADseq) of predominantly Antarctic lichen-forming fungi, our main aim in this study was to clarify the taxonomy of asexual Usnea antarctica and sexual Usnea aurantiacoatra and test their hypothesised species pair relationship. To further support our findings, we applied population genomic methods to measure the degree of genomic divergence and infer the levels of co-ancestry for each species.


Sample collection and site description

Samples were collected in Antarctica between December 2015 and January 2017. From these collected specimens, we chose to compare 105 representative specimens of the species U. antarctica and U. aurantiacoatra for this study (see details of specimens in Suppl. material 1). All selected specimens were either collected on King George Island (65) and Elephant Island (19) of the South Shetland Islands or in the Northern part of the Antarctic Peninsula (21) near the research bases Esperanza and Primavera. Fifty-eight of the 105 selected specimens were identified as U. antarctica and 47 specimens were identified as U. aurantiacoatra based on their phenotypical characters (Walker 1985). As a reference sequence to filter for lichen-fungal loci of U. antarctica and U. aurantiacoatra during the RADseq processing, we sequenced a specimen of U. strigosa that was collected in Arkansas, U.S.A. (Suppl. material 1).

DNA extraction

Total metagenomic DNA was extracted either by following a cetyltrimethylammonium bromide (CTAB) protocol as modified by Cubero et al. (1999) or by using the ZR Fungal/Bacterial DNA MiniPrep Kit (Zymo Research, Irvine, CA, USA) as recommended by the manufacturer. We used the whole lichen thalli for DNA extraction from U. antarctica and U. aurantiacoatra, but only the central axis in U. strigosa to preferentially extract DNA from the lichen fungus (to avoid the photobiont). DNA concentrations of all samples were quantified with a Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA, USA).

Reference Sequencing and Assembly

We first deep-sequenced and assembled a reference sequence of an Usnea strigosa specimen to aid in mapping lichen-fungal loci during the processing of metagenomic RADseq data. A paired-end Illumina sequencing library of 150 bp read length was constructed from total DNA with the Nextera XT DNA Library Prep Kit (Illumina, San Diego, CA, USA) and sequenced on a NextSeq platform at the University of Illinois Chicago’s Sequencing Core Facility (Chicago, IL, USA). The resulting reads were quality trimmed using the programme Trimmomatic v0.36 (Bolger et al. 2014). Bases were trimmed when the average quality of 4-base sliding windows was below 15 and bases at the start and end of reads had a quality below 10. Subsequently, all trimmed reads, shorter than 25 bp, were filtered out (LEADING:10 TRAILING:10 SLIDINGWINDOW:4:15 MINLEN:25). The trimmed reads were used for a genome assembly with the programme SPAdes v3.5 (Bankevich et al. 2012) with default parameter settings. The assembled metagenomic scaffolds were loaded into the programme MetaWatt v3.5.3 (Strous et al. 2012) for a binning based on tetranucleotide frequencies. Scaffolds of fungal origin that clustered together were separated from the remaining scaffolds. All selected scaffolds that were larger than 10 kb were then included into the final reference sequence of U. strigosa. We used the Core Eukaryotic Gene Mapping Approach (CEGMA) to estimate the genomic completeness of the assembly (Parra et al. 2009). Finally, we created a Bowtie2 (Langmead and Salzberg 2012) database from the selected scaffolds for the mapping approach to filter for fungal RAD loci.

RADseq Library Preparation and Sequencing

RADseq libraries for Usnea antarctica and U. aurantiacoatra were prepared as described previously (Grewe et al. 2017). In short, for the RADseq library production, DNA isolations were pooled with sequence adapters (Rubin and Moreau 2016), subsequently digested with the restriction enzyme ApeKI (New England Biolabs, Ipswich, MA, USA) and ligated using T4 ligase (New England Biolabs). Up to 48 samples with compatible barcodes were pooled and selected for fragments of sizes between 300 and 500 bp using the BluePippin DNA size selection system (Sage Science, Beverly, MA, USA). The pooled libraries were amplified using the REDTaq ReadyMix (Sigma-Aldrich, St. Louis, MO, USA) prior to sequencing on an Illumina MiSeq using the MiSeq Reagent Kit v3 for 150 cycles (Illumina, San Diego, CA, USA) to produce single-end sequences with a length of 150 bp.

Assembly of RADseq datasets

The raw reads of U. antarctica and U. aurantiacoatra from the MiSeq sequencing were processed and assembled as described earlier for metagenomic datasets of lichens (Grewe et al. 2017). This process used a combination of the ipyRAD ( and pyRAD (Eaton and Ree 2013) pipelines with an additional mapping step that filtered for lichen-fungal loci with a reference sequence. Subsequently, we refer to the raw Illumina RAD sequences as ‘read’ and name the clustered reads per individual sample ‘loci’; the final matrices are alignments of homologous loci from multiple samples with nucleotide substitutions referred to as ‘SNP’. In pyRAD, we set the datatype to genotype-by-sequencing (gbs), ploidy to haploid (1), a similarity threshold for the clustering of reads within and between individuals to 90% (.90) and a minimum coverage of four individuals per locus (4). For the reference-based filtering of RAD loci, we used Bowtie2 with adjusted parameters to allow one permitted mismatch (−N 1), a seed length of 20 (−L 20), up to 20 seed extension attempts (−D 20) and a maximum “re-seeding” of 3 (−R 3). Following an initial round with all sequenced samples, we re-ran step 7 of pyRAD and excluded samples with less than 1000 recovered loci. We used the filtered pyRAD output files, such as unlinked_snps, alleles and vcf, for further analyses.

Phylogenetic reconstructions

Phylogenetic trees were calculated from all unlinked SNPs of the filtered RADseq dataset, i.e. a matrix that was limited to one SNP per RAD locus. This matrix was used for a RAxML v7.2.8 (Stamatakis 2006) maximum likelihood analysis using the GTR + G model. For each analysis, 100 bootstrap replicates were calculated using the fast bootstrapping option implemented in RAxML (Stamatakis et al. 2008). The resulting phylogenetic tree was midpoint rooted and drawn to scale with FigTree v1.4.3 (

Analysis of population structure

To calculate differences in the population structure between U. antarctica and U. aurantiacoatra, we created a reduced dataset that included all sites with a minor allele frequency (MAF) greater than or equal to 0.05 and greater than 50% coverage using vcftools v0.1.15 (Danecek et al. 2011). This reduced vcf file was converted into a genind object from the R package adegenet v2.0.2 (Jombart and Ahmed 2011; Jombart et al. 2010). The genind object was appended with additional information settings for haploid genomes and the population memberships for samples according to their initial identification based on morphological characters. With all information enclosed, the genind object became subject to population genetics analyses encoded in R.

To determine the degree to which both populations are subdivided, we estimated Gst (Nei 1973) and Hedrick’s standardised genetic differentiation measures G'st (Hedrick 2005) and Jost’s D (Jost 2008) by using the R package mmod v1.3.3 (Winter 2012). Gst is a good measure when the mutation rate is small relative to migration rate; contrarily, G’st and D fit to data with high mutation rates and two populations (Whitlock 2011). We used these multiple statistics to get a comprehensive measure of genomic differentiation. In a population pairwise comparison, we calculated these indices per site and plotted the values by frequency in separate histograms for Gst, G’st and D.

The genetic structure of samples of U. antarctica and U. aurantiacoatra was evaluated with the Discriminant Analysis of Principal Components (DAPC) implemented in the R package adegenet v2.02. This non-parametric method first transforms the data using a principle components analysis (PCA) and subsequently distinguishes between two or more groups using a discriminant analysis (DA). The DAPC was conducted by using the first 60 principal components and all (two) DA-eigenvalues. In addition to the display of the genetic variation in genomic space, the DAPC allows a prediction of the group membership probability for each sample which is visualised in a STRUCTURE-like plot.

In addition to the nonparametric approach with DAPC, we used a model-based method to detect population subdivision using the programme fineRADstructure (Malinsky et al. 2018). This software is specifically designed to measure population structure amongst haplotypes inferred from RADseq datasets. We used the script included in fineRADstructure-tools ( to convert the pyRAD alleles output into the input format for fineRADstructure. During the conversion, we also reduced the dataset to only unlinked loci (default parameter) with a minimum sample number of 4 (--minsample 4). As recommended by the authors, we then re-ordered the unsorted RAD loci with the script sampleLD.R which is part of the fineRADstructure package. Next, we used the scripts, RADpainter and fineSTRUCTURE, which are both implemented in fineRADstructure, to measure the population structure. First, we calculated the co-ancestry matrix with RADpainter for haploid datasets (-p 1). We then used fineSTRUCTURE for the Markov chain Monte Carlo (MCMC) clustering algorithm with the following arguments: -x 100,000, -z 100,000 and -y 1,000. We also started fineSTRUCTURE with the arguments -m T and -x 10,000 to run a simple tree-building algorithm on the data of the co-ancestry matrix. Finally, the co-ancestry matrix, MCMC output and the coalescence tree were loaded into the programme ‘Finestructure GUI’ for visualisation.


The U. strigosa reference sequence and all scripts that were used in this study are available online ( All RAD sequences were deposited in the NCBI Sequence Read Archive (SRA) with accession number PRJNA505526.


Reference assembly and RADseq results

We assembled a draft reference genome of U. strigosa to filter for fungal RAD loci from U. antarctica and U. aurantiacoatra. The Illumina NextSeq sequencing of the whole U. strigosa lichen resulted in 8,552,530 metagenomic paired-end reads. First, we trimmed these raw data which reduced the paired-end reads to 8,366,962 (97.78% of raw data). The trimmed read pairs were then assembled into 16,932 scaffolds (N50 = 12,750 bp) with a total size of 40.9 Mbp (including 1,187 scaffolds of sizes larger than 10 kb). Metagenomic binning identified 28.92 Mbp of the assembly as fungal derived from which we selected 1,100 scaffolds (N50 = 23,562 bp) with sizes larger than 10 kb; all but two of these scaffolds were continuous assemblies (contigs). The sorted draft genome of U. strigosa had a total size of 24.1 Mbp and an estimated completeness of 72.18%.

We included 105 specimens of U. antarctica and U. aurantiacoatra that were collected in Antarctica in four RADseq libraries (Suppl. material 2). The sequence read number of each sample varied widely from 13,659 for sample EL0059 to 1,942,819 for sample EL0074 with an average sequence read number of 488,468 (sd = 313,604). The number of loci (within sample clusters) that pyRAD generated from these sequences directly correlated with the initial number of sequences (R2 = 0.8017, Suppl. material 3). An average of 21.8% (sd = 2.9%) of all loci mapped to the lichen fungus reference genome and, of these, an average of 85.4% (sd = 5.5%) were included into the final pyRAD dataset. The numbers of loci before and after the mapping were directly correlated (R2 = 0.7598, Suppl. material 3); however, the number of mapped loci reached saturation at an average of 6,496 (sd = 801) for samples with more than 40,000 initial loci. In addition, the number of mapped loci were strongly correlated to the number of loci included in the final dataset (R2 = 0.9869, Suppl. material 3). Two samples of U. antarctica (EL0059, EL0281) and two samples of U. aurantiacoatra (EL0415, EL0437) had less than 1,000 loci in the final dataset and were removed from the analysis. All remaining 101 samples in the final dataset had on average of 4,143 (sd = 1,316) loci (Suppl. material 2).

Phylogenetic analysis of RADseq data

The phylogenetic analysis of the RADseq data showed two distinct and highly supported clades corresponding to the phenotypically circumscribed species U. antarctica and U. aurantiacoatra (Figure 1). The phylogenetic tree was calculated from a matrix with 7,087 positions and 53.24% gaps. Most internal relationships within each clade remained unresolved; however, the U. antarctica clade showed higher internal support values than the U. aurantiacoatra clade. Within the U. antarctica clade, three sister relationships of U. antarctica (EL0001 and EL0409, EL0382 and EL0390, EL0713 and EL0743) had a 100% bootstrap support and short branches, indicating low genomic divergence.

Figure 1. 

Phylogenetic tree inferred from the U. antarctica and U. aurantiacoatra RADseq data. The clades of each species are highlighted by brackets. Bootstrap values are indicated at the branches. The unit of branch length is substitutions per site. Note that branches leading to both major clades were abbreviated by 0.4 substitutions per site.

Population genomic analyses of RADseq data

We determined the degree to which both species complexes are subdivided by Gst, G’st and D measurements. For these analyses, we included only SNPs with a MAF greater than 0.05 and more than 50% coverage. This reduced the RAD dataset to a total of 4,132 SNPs. We plotted the frequency of Gst, G’st and D measures for each SNP (Figure 2). A strong tendency towards 1 for most SNPs in all three measures strongly indicated that genomes of both species were completely isolated. This was also supported by the average measures of Gst, G’st and D of 0.70, 0.93 and 0.60, respectively.

Figure 2. 

Pairwise Gst, G’st and D distribution. Pairwise values of Nei’s Gst (green), Hedrick’s G’st (blue) and Jost’s D (yellow) are plotted by their frequency.

The same reduced dataset of 4,132 SNPs was used to differentiate the genomes by their variation in a non-parametric approach with a DAPC (Figure 3A). The DAPC combines a PCA with a DA for a separation of genomes based on their variance between groups rather than the total variance of the sample. The resulting clusters of both species were clearly separated in genomic space and showed no evidence for admixture. In addition, the group membership probabilities indicated absolute discrimination of the two species by the DAPC assigning each individual with 100% probability to their respective species (Figure 3B).

Figure 3. 

Genomic variation by non-parametric DAPC. A DAPC plot of the densities of U. antarctica (blue) and U. aurantiacoatra (green) on the first retained discriminant function B Bar plot of group membership probabilities.

The separation of U. antarctica and U. aurantiacoatra was further supported by the results of a Bayesian model-based approach with the programme fineRADstructure. By converting the pyRAD allele output for fineRADstructure, we reduced the dataset to 3,803 unlinked SNPs with a minimum coverage of 4 samples. The resulting clustered co-ancestry matrix showed that both species shared more co-ancestry within each other than between species (Figure 4). By comparing both species clusters, U. aurantiacoatra showed a higher estimated co-ancestry than U. antarctica (Figure 4A). To avoid a sampling bias, we reduced the dataset for the fineRADstructure analysis to include only samples collected on King George Island and Elephant Island. This reduced the dataset to 80 samples and 3,652 unlinked SNPs with a minimum coverage of 4 samples. The resulting plot of the reduced dataset also showed higher shared co-ancestry within each species compared to that between species, but estimated higher co-ancestry of U. antarctica than U. aurantiacoatra (Figure 4B), opposite to the entire dataset. In addition, both matrices visualised different degrees of intraspecific co-ancestry and suggested substructure for a group of three specimens of U. antarctica from Potter Peninsula, King George Island (EL0022, EL0034 and EL0051) and for six specimens of U. aurantiacoatra from Fildes Peninsula, King George Island (EL0444, EL0435, EL0445, EL450, EL0455 and EL0453). Moreover, two U. antarctica specimen pairs collected on King George Island (EL0001 and EL0409, EL0382 and EL390) and one pair collected on the Antarctic Peninsula (EL0713 and EL0743) showed the highest degrees of co-ancestry demonstrating very close relatedness, such as sister or clonal relationships. These results agreed with the phylogenetic inference (see above) in which the same U. antarctica specimens were close sister taxa.

Figure 4. 

Clustered fineRADstructure co-ancestry matrix. A Full dataset including U. antarctica collected on the Antarctic Peninsula in addition to U. antarctica and U. aurantiacoatra collected on King George Island and Elephant Island B Reduced dataset with all U. antarctica and U. aurantiacoatra collected on King George Island and Elephant Island. Two major clades are corresponding to the two species U. antarctica (top-left) and U. aurantiacoatra (bottom-right). The top and left trees were calculated from the co-ancestry matrix to sort the individuals by their population structure. The matrix is diagonally split into the top-right half showing raw data and the bottom-left half displaying aggregated data.


In this study, we used RAD sequencing for evaluating the delimitation of two predominantly Antarctic Usnea species. Phylogenetic evidence and population genomic analyses of the RADseq data strongly supported that the two species represent independent lineages. Although both species showed no overlapping genomic structure in a DAPC, we could compare levels of co-ancestry and detect genomic substructure within each species in a fineRADstructure plot.

In previous studies using multi-locus approaches, including the ITS barcoding marker (Schoch et al. 2012), the relationship of U. antarctica and U. aurantiacoatra remained unresolved and, since specimens of both species did not separate as different clades, conspecificity of the species was not ruled out (Seymour et al. 2007; Wirtz et al. 2012). Our study using RADseq supports the results obtained using microsatellite data that suggested the two species are distinct lineages (Lagostina et al. 2018). In U. antarctica and U. aurantiacoatra, the taxonomic interpretation of species pairs as separate species (Poelt 1972) is supported.

We developed a RADseq method for lineages involved in intimate symbiotic associations (Grewe et al. 2017), which we here successfully implemented for the use of delimiting two species. Different to the previously described RADseq method that used a reference genome from a lichen-fungal culture, we successfully generated a reference genome from a metagenomic de-novo assembly of U. strigosa. The filtering of the metagenomic assembly for fungal derived content reduced the size and completeness of the fungal reference (28.92 Mbp, CEGMA: 72.18%) compared to the reference genome assembly from a lichen-fungal culture which was used in earlier studies (31.6 Mbp, CEGMA: 96.77%) (Grewe et al. 2017; Leavitt et al. 2016b). However, the saturation of successfully mapped loci to the reference (Suppl. material 3) suggested that the maximum number of possible mapped loci was reached for samples with many initial loci. Therefore, although using a smaller reference and less fungal derived loci than in our initial study (Grewe et al. 2017), this RADseq approach still was successful in mapping a large number of fungal loci sufficient for phylogenetic and population genomic methods. This widens the potential application of RADseq for intimate symbiotic organisms and includes studies where cultures of one symbiotic partner are not readily available.

RADseq data are extremely powerful, since the method generates a matrix of thousands of homologous loci derived from randomly distributed regions across the genome. Many studies have successfully used large RADseq datasets for phylogenetic analysis which were difficult to resolve due to insufficient signals in available markers (Eaton and Ree 2013; Escudero et al. 2014; Hipp et al. 2014; Vargas et al. 2017; Wagner et al. 2018). Our phylogenetic and population genomic results from the RADseq dataset clearly delimited U. antarctica and U. aurantiacoatra into two lineages (Figures 14) supporting the acceptance of two species. This confirms that closely related species are difficult to separate using sequence-based multi-locus approaches and great care should be taken when interpreting results from molecular studies when it comes to testing for conspecificity. On the other hand, the microsatellite-based multi-locus approach by Lagostina et al. (2018) rendered almost identical results, including nearly 100 % correct assignment of samples to their species.

The fineRADstructure matrix estimated lower co-ancestry (and hence higher genotypic variation) for the sexually-reproducing U. aurantiacoatra, compared to the asexually-reproducing U. antarctica when comparing samples that were collected in the same geographic range (Figure 4B). This result agrees with earlier observations that asexual populations have lower genotypic variation than sexual populations in modelling approaches (Balloux et al. 2003) and empirical measures (Delmotte et al. 2002). Moreover, Lagostina et al. (2018) inferred lower genetic variability for U. antarctica than U. aurantiacoatra using 23 microsatellite loci. These authors also used samples collected in mixed stands of both species from King George and Elephant Island. When we increased our sampling of U. antarctica to include a much wider geographical range (Antarctic Peninsula in addition to King George and Elephant Island) compared to the sampling of U. aurantiacoatra (King George and Elephant Island only), the matrix indicated increased levels of co-ancestry and a lower genotypic variation (Figure 4A). Although this comparative analysis is lacking collections of U. aurantiacoatra from the Antarctic Peninsula for a direct comparison, it should be noted that U. antarctica covers a wider geographical range than U. aurantiacoatra (Walker 1985) and this wider species distribution might increase genetic variability. The difference in distribution may result from the main form of reproductive units of both Usnea. The exclusively sexual U. aurantiacoatra reproduces via the dispersal of fungal spores which are required to meet with an appropriate photobiont after germination. The asexual U. antarctica on the other hand is in majority vegetatively reproducing via soredia, which already include the photobiont. Therefore, even if both reproductive units are dispersed over similar distances, the success rate of colonisation may be higher for soredia and explain the overall wider distribution and therefore genetic variability of U. antarctica. Finally, it was predicted that a small number of sexual individuals per generation — and U. antarctica rarely can be found with apothecia — is sufficient to make an apparently asexual population highly variable (Bengtsson 2003).

Despite the lower co-ancestry of U. antarctica compared to U. aurantiacoatra, we detected three pairs of very close relatives with high co-ancestry of U. antarctica (Figure 4). The three pairs were collected on Elephant Island, King George Island and on the Antarctic Peninsula, respectively and may indicate almost immediately related clones. On Elephant Island and the Antarctic Peninsula, the pairs were collected in the same locations with a greater chance to pick up clones. However, the clonal pair from King George Island must have dispersed between Fildes and Potter Peninsula over ice or water boundaries prior to our collection. Contrarily, none of the individuals of U. aurantiacoatra expressed similarly close relationships. However, we could detect substructure for a group of six individuals of U. aurantiacoatra collected at the same location and three specimens of U. antarctica collected at different locations on King George Island, which indicates the potential of this analysis to identify (sub)population structure. Using this detailed method to measure co-ancestry on a deeper sampling of individuals of Usnea may, in future, provide a comprehensive picture of population structure and diversification.


We successfully used RADseq for phylogenetic and population genomic studies on two species of the lichen-fungal genus Usnea. Phylogenetic inference using RAD data clearly delimited the species U. antarctica and U. aurantiacoatra into two lineages, which were irresolvable using multi-locus DNA sequence markers. Furthermore, the RADseq approach offered sufficient genotyping data for conclusive population genomic analyses. We used RADseq to measure lower co-ancestry in the asexual U. antarctica than in the sexual U. aurantiacoatra, potentially derived from a wider geographical distribution of U. antarctica in our sample. These results show that RADseq has much potential for future phylogenetic and population genomic studies on lichens, particularly for groups of organisms which remained unresolved by multi-locus markers.


The authors thank the Robert A. Pritzker Center for Meteoritics and Polar Studies and Tawani Foundation for financial support. We want to acknowledge the Lauer family for the generous gift of a MiSeq sequencer to the Pritzker Lab for Molecular Systematics at the Field Museum. Fieldwork by Elisa Lagostina was funded by the German Science Foundation through grant Pr 567/18-1. Mikhael Andreev, Aline Lorenz and Mayara Scur are gratefully acknowledged for contributing samples from Elephant and King George Island and the Antarctic Peninsula. We also thank Marta Ribeiro and Kevin Feldheim (Chicago) for their invaluable help in the laboratory and Sarah Cowles (Miami), Steve D. Leavitt (Provo) and Pradeep K. Divakar (Madrid) for comments on the manuscript.


  • Articus K, Mattsson JE, Tibell L, Grube M, Wedin M (2002) Ribosomal DNA and beta-tubulin data do not support the separation of the lichens Usnea florida and U. subfloridana as distinct species. Mycological Research 106: 412–418.
  • Baird NA, Etter PD, Atwood TS, Currey MC, Shiver AL, Lewis ZA, Selker EU, Cresko WA, Johnson EA (2008) Rapid SNP discovery and genetic mapping using sequenced RAD markers. Plos One 3(10): e3376.
  • Balloux F, Lehmann L, de Meeûs T (2003) The Population Genetics of Clonal and Partially Clonal Diploids. Genetics 164: 1635–1644.
  • Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, Lesin VM, Nikolenko SI, Son P, Prjibelski AD, Pyshkin AV, Sirotkin AV, Vyahhi N, Tesler G, Alekseyev MA, Pevzner PA (2012) SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing. Journal of Computational Biology 19: 455–477.
  • Bickford D, Lohman DJ, Sodhi NS, Ng PKL, Meier R, Winker K, Ingram KK, Das I (2007) Cryptic species as a window on diversity and conservation. Trends in Ecology & Evolution 22: 148–155.
  • Buschbom J, Mueller GM (2006) Testing “species pair” hypotheses: Evolutionary processes in the lichen-forming species complex Porpidia flavocoerulescens and Porpidia melinodes. Molecular Biology and Evolution 23: 574–586.
  • Crespo A, Lumbsch HT, Mattsson JE, Blanco O, Divakar PK, Articus K, Wiklund E, Bawingan PA, Wedin M (2007) Testing morphology-based hypotheses of phylogenetic relationships in Parmeliaceae (Ascomycota) using three ribosomal markers and the nuclear RPB1 gene. Molecular Phylogenetics and Evolution 44: 812–824.
  • Crespo A, Pérez-Ortega S (2009) Cryptic species and species pairs in lichens: A discussion on the relationship between molecular phylogenies and morphological characters. Anales Del Jardin Botanico De Madrid 66: 71–81.
  • Cubero OF, Crespo A, Fatehi J, Bridge PD (1999) DNA extraction and PCR amplification method suitable for fresh, herbarium-stored, lichenized, and other fungi. Plant Systematics and Evolution 216: 243–249.
  • Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, McVean G, Durbin R, 1000 Genomes Project Analysis Group (2011) The variant call format and VCFtools. Bioinformatics 27: 2156–2158.
  • Davey JW, Hohenlohe PA, Etter PD, Boone JQ, Catchen JM, Blaxter ML (2011) Genome-wide genetic marker discovery and genotyping using next-generation sequencing. Nature Reviews Genetics 12: 499–510.
  • Delmotte F, Leterme N, Gauthier JP, Rispe C, Simon JC (2002) Genetic architecture of sexual and asexual populations of the aphid Rhopalosiphum padi based on allozyme and microsatellite markers. Molecular Ecology 11: 711–723.
  • Eaton DAR, Ree RH (2013) Inferring phylogeny and introgression using RADseq data: An example from flowering plants (Pedicularis: Orobanchaceae). Systematic Biology 62: 689–706.
  • Escudero M, Eaton DAR, Hahn M, Hipp AL (2014) Genotyping-by-sequencing as a tool to infer phylogeny and ancestral hybridization: A case study in Carex (Cyperaceae). Molecular Phylogenetics and Evolution 79: 359–367.
  • Grewe F, Huang JP, Leavitt SD, Lumbsch HT (2017) Reference-based RADseq resolves robust relationships among closely related species of lichen-forming fungi using metagenomic DNA. Scientific Reports 7: 9884.
  • Hibbett DS, Binder M, Bischoff JF, Blackwell M, Cannon PF, Eriksson OE, Huhndorf S, James T, Kirk PM, Lücking R, Lumbsch HT, Lutzoni F, Matheny PB, McLaughlin DJ, Powell MJ, Redhead S, Schoch CL, Spatafora JW, Stalpers JA, Vilgalys R, Aime MC, Aptroot A, Bauer R, Begerow D, Benny GL, Castlebury LA, Crous PW, Dai YC, Gams W, Geiser DM, Griffith GW, Gueidan C, Hawksworth DL, Hestmark G, Hosaka K, Humber RA, Hyde KD, Ironside JE, Koljalg U, Kurtzman CP, Larsson KH, Lichtwardt R, Longcore J, Miadlikowska J, Miller A, Moncalvo JM, Mozley-Standridge S, Oberwinkler F, Parmasto E, Reeb V, Rogers JD, Roux C, Ryvarden L, Sampaio JP, Schussler A, Sugiyama J, Thorn RG, Tibell L, Untereiner WA, Walker C, Wang Z, Weir A, Weiss M, White MM, Winka K, Yao YJ, Zhang N (2007) A higher-level phylogenetic classification of the Fungi. Mycological Research 111: 509–547.
  • Hipp AL, Eaton DAR, Cavender-Bares J, Fitzek E, Nipper R, Manos PS (2014) A framework phylogeny of the American Oak clade based on sequenced RAD data. Plos One 9(4): e102272.
  • Hodel RDGJ, Segovia-Salcedo MC, Landis JB, Crowl AA, Sun M, Liu X, Gitzendanner MA, Douglas NNA, Germain-Aubrey CC, Chen S, Soltis DE, Soltis PS (2016) The report of my death was an exaggeration: a review for researchers using microsatellites in the 21st century. Applications in Plant Sciences 4.
  • Holder M, Lewis PO (2003) Phylogeny estimation: Traditional and Bayesian approaches. Nature Reviews Genetics 4: 275–284.
  • Huelsenbeck JP, Ronquist F, Nielsen R, Bollback JP (2001) Evolution – Bayesian inference of phylogeny and its impact on evolutionary biology. Science 294: 2310–2314.
  • James TY, Kauff F, Schoch CL, Matheny PB, Hofstetter V, Cox CJ, Celio G, Gueidan C, Fraker E, Miadlikowska J, Lumbsch HT, Rauhut A, Reeb V, Arnold AE, Amtoft A, Stajich JE, Hosaka K, Sung GH, Johnson D, O’Rourke B, Crockett M, Binder M, Curtis JM, Slot JC, Wang Z, Wilson AW, Schussler A, Longcore JE, O’Donnell K, Mozley-Standridge S, Porter D, Letcher PM, Powell MJ, Taylor JW, White MM, Griffith GW, Davies DR, Humber RA, Morton JB, Sugiyama J, Rossman AY, Rogers JD, Pfister DH, Hewitt D, Hansen K, Hambleton S, Shoemaker RA, Kohlmeyer J, Volkmann-Kohlmeyer B, Spotts RA, Serdani M, Crous PW, Hughes KW, Matsuura K, Langer E, Langer G, Untereiner WA, Lücking R, Büdel B, Geiser DM, Aptroot A, Diederich P, Schmitt I, Schultz M, Yahr R, Hibbett DS, Lutzoni F, McLaughlin DJ, Spatafora JW, Vilgalys R (2006) Reconstructing the early evolution of Fungi using a six-gene phylogeny. Nature 443: 818–822.
  • Kotelko R, Piercey-Normore MD (2010) Cladonia pyxidata and C. pocillum; genetic evidence to regard them as conspecific. Mycologia 102: 534–545.
  • Krog H (1976) Lethariella and Protousnea, two new lichen genera in the Parmeliaceae. Norwegian Journal of Botany 23: 83–106.
  • Krog H (1982) Evolutionay trends in foliose and fruticose lichens of the Parmeliaceae. Journal of the Hattori Botanical Laboratory 52: 303–311.
  • Lamb IM (1964) Antarctic lichens. I. The genera Usnea, Ramalina, Himantormia, Alectoria, Cornicularia. 38: 1–34.
  • Leavitt SD, Divakar PK, Crespo A, Lumbsch HT (2016a) A matter of time – understanding the limits of the power of molecular data for delimiting species boundaries. Herzogia 29: 479–492.
  • Leavitt SD, Grewe F, Widhelm T, Muggia L, Wray B, Lumbsch HT (2016b) Resolving evolutionary relationships in lichen-forming fungi using diverse phylogenomic datasets and analytical apporaches. Scientific Reports 6: 22262.
  • Lumbsch HT, Wirtz N (2011) Phylogenetic relationships of the neuropogonoid core group in the genus Usnea (Ascomycota: Parmeliaceae). Lichenologist 43: 553–559.
  • Lutzoni F, Kauff F, Cox C, McLaughlin D, Celio G, Dentinger B, Padamsee M, Hibbett D, James TY, Baloch E, Grube M, Reeb V, Hofstetter V, Schoch C, Arnold AE, Miadlikowska J, Spatafora J, Johnson D, Hambleton S, Crockett M, Shoemaker R, Sung G-H, Lücking R, Lumbsch HT, O’Donnell K, Binder M, Diederich P, Ertz D, Gueidan C, Hansen K, Harris RC, Hosaka K, Lim Y-W, Matheny B, Nishida H, Pfister D, Rogers J, Rossman A, Schmitt I, Sipman H, Stone J, Sugiyama J, Yahr R, Vilgalys R (2004) Assembling the fungal tree of life: progress, classification, and evolution of subcellular traits. American Journal of Botany 91: 1446–1480.
  • Malinsky M, Trucchi E, Lawson DJ, Falush D (2018) RADpainter and fineRADstructure Population Inference from RADseq Data. Molecular Biology and Evolution 35: 1284–1290.
  • Myllys L, Lohtander K, Tehler A (2001) beta-tubulin, ITS and group I intron sequences challenge the species pair concept in Physcia aipolia and P. caesia. Mycologia 93: 335–343.
  • Poelt J (1970) Das Konzept der Artenpaare bei den Flechten. Vorträge aus dem Gesamtgebiet der Botanik. Neue Folge 4: 187–198.
  • Poelt J (1972) Die taxonomische Behandlung von Artenpaaren bei den Flechten. Botaniska Notiser 125: 77–81.
  • Rubin BER, Moreau CS (2016) Comparative genomics reveals convergent rates of evolution in ant-plant mutualisms. Nature Communications 7: 12679.
  • Schoch CL, Seifert KA, Huhndorf S, Robert V, Spouge JL, Levesque CA, Chen W, Bolchacova E, Voigt K, Crous PW, Miller AN, Wingfield MJ, Aime MC, An KD, Bai FY, Barreto RW, Begerow D, Bergeron MJ, Blackwell M, Boekhout T, Bogale M, Boonyuen N, Burgaz AR, Buyck B, Cai L, Cai Q, Cardinali G, Chaverri P, Coppins BJ, Crespo A, Cubas P, Cummings C, Damm U, de Beer ZW, de Hoog GS, Del-Prado R, Dentinger B, Dieguez-Uribeondo J, Divakar PK, Douglas B, Duenas M, Duong TA, Eberhardt U, Edwards JE, Elshahed MS, Fliegerova K, Furtado M, Garcia MA, Ge ZW, Griffith GW, Griffiths K, Groenewald JZ, Groenewald M, Grube M, Gryzenhout M, Guo LD, Hagen F, Hambleton S, Hamelin RC, Hansen K, Harrold P, Heller G, Herrera G, Hirayama K, Hirooka Y, Ho HM, Hoffmann K, Hofstetter V, Högnabba F, Hollingsworth PM, Hong SB, Hosaka K, Houbraken J, Hughes K, Huhtinen S, Hyde KD, James T, Johnson EM, Johnson JE, Johnston PR, Jones EB, Kelly LJ, Kirk PM, Knapp DG, Koljalg U, Kovacs GM, Kurtzman CP, Landvik S, Leavitt SD, Liggenstoffer AS, Liimatainen K, Lombard L, Luangsa-Ard JJ, Lumbsch HT, Maganti H, Maharachchikumbura SS, Martin MP, May TW, McTaggart AR, Methven AS, Meyer W, Moncalvo JM, Mongkolsamrit S, Nagy LG, Nilsson RH, Niskanen T, Nyilasi I, Okada G, Okane I, Olariaga I, Otte J, Papp T, Park D, Petkovits T, Pino-Bodas R, Quaedvlieg W, Raja HA, Redecker D, Rintoul T, Ruibal C, Sarmiento-Ramirez JM, Schmitt I, Schussler A, Shearer C, Sotome K, Stefani FO, Stenroos S, Stielow B, Stockinger H, Suetrong S, Suh SO, Sung GH, Suzuki M, Tanaka K, Tedersoo L, Telleria MT, Tretter E, Untereiner WA, Urbina H, Vagvolgyi C, Vialle A, Vu TD, Walther G, Wang QM, Wang Y, Weir BS, Weiss M, White MM, Xu J, Yahr R, Yang ZL, Yurkov A, Zamora JC, Zhang N, Zhuang WY, Schindel D, Fungal Barcoding Consortium (2012) Nuclear ribosomal internal transcribed spacer (ITS) region as a universal DNA barcode marker for Fungi. Proceedings of the National Academy of Sciences of the United States of America 109: 6241–6246.
  • Schoch CL, Sung GH, Lopez-Giraldez F, Townsend JP, Miadlikowska J, Hofstetter V, Robbertse B, Matheny PB, Kauff F, Wang Z, Gueidan C, Andrie RM, Trippe K, Ciufetti LM, Wynns A, Fraker E, Hodkinson BP, Bonito G, Groenewald JZ, Arzanlou M, de Hoog GS, Crous PW, Hewitt D, Pfister DH, Peterson K, Gryzenhout M, Wingfield MJ, Aptroot A, Suh SO, Blackwell M, Hillis DM, Griffith GW, Castlebury LA, Rossman AY, Lumbsch HT, Lücking R, Büdel B, Rauhut A, Diederich P, Ertz D, Geiser DM, Hosaka K, Inderbitzin P, Kohlmeyer J, Volkmann-Kohlmeyer B, Mostert L, O’Donnell K, Sipman H, Rogers JD, Shoemaker RA, Sugiyama J, Summerbell RC, Untereiner W, Johnston PR, Stenroos S, Zuccaro A, Dyer PS, Crittenden PD, Cole MS, Hansen K, Trappe JM, Yahr R, Lutzoni F, Spatafora JW (2009) The Ascomycota tree of life: A phylum-wide phylogeny clarifies the origin and evolution of fundamental reproductive and ecological traits. Systematic Biology 58: 224–239.
  • Seymour FA, Crittenden PD, Wirtz N, Øvstedal DO, Dyer PS, Lumbsch HT (2007) Phylogenetic and morphological analysis of Antarctic lichen-forming Usnea species in the group Neuropogon. Antarctic Science 19: 71–82.
  • Spatafora JW, Aime MC, Grigoriev IV, Martin F, Stajich JE, Blackwell M (2017) The fungal tree of life: from molecular systematics to genome-scale phylogenies. Microbiology Spectrum 5(5): FUNK-0053-2016.
  • Spatafora JW, Robbertse B (2010) Chapter 4 : phylogenetics and phylogenomics of the Fungal Tree of Life. In: Borkovich KA, Ebbole DJ (Eds) Cellular and Molecular Biology of Filamentous Fungi. ASM Press, 36–49.
  • Strous M, Kraft B, Bisdorf R, Tegetmeyer HE (2012) The binning of metagenomic contigs for microbial physiology of mixed cultures. Frontiers in Microbiology 3: 410.
  • Taylor JW, Jacobson DJ, Kroken S, Kasuga T, Geiser DM, Hibbett DS, Fisher MC (2000) Phylogenetic species recognition and species concepts in fungi. Fungal Genetics and Biology 31: 21–32.
  • Vargas OM, Ortiz EM, Simpson BB (2017) Conflicting phylogenomic signals reveal a pattern of reticulate evolution in a recent high-Andean diversification (Asteraceae: Astereae: Diplostephium). New Phytologist 214: 1736–1750.
  • Velmala S, Myllys L, Halonen P, Goward T, Ahti T (2009) Molecular data show that Bryoria fremontii and B. tortuosa (Parmeliaceae) are conspecific. Lichenologist 41: 231–242.
  • Wachi N, Matsubayashi KW, Maeto K (2018) Application of next-generation sequencing to the study of non-model insects. Entomological Science 21: 3–11.
  • Wagner ND, Gramlich S, Hörandl E (2018) RAD sequencing resolved phylogenetic relationships in European shrub willows (Salix L. subg. Chamaetia and subg. Vetrix) and revealed multiple evolution of dwarf shrubs. Ecology and Evolution 2018(00): 1–13.
  • Walker FJ (1985) The lichen genus Usnea subgenus Neuropogon. Bulletin of the British Museum (Natural History), Botany Series 13: 1–130.
  • Wickett NJ, Mirarab S, Nam N, Warnow T, Carpenter E, Matasci N, Ayyampalayam S, Barker MS, Burleigh JG, Gitzendanner MA, Ruhfel BR, Wafula E, Der JP, Graham SW, Mathews S, Melkonian M, Soltis DE, Soltis PS, Miles NW, Rothfels CJ, Pokorny L, Shaw AJ, DeGironimo L, Stevenson DW, Surek B, Villarreal JC, Roure B, Philippe H, dePamphilis CW, Chen T, Deyholos MK, Baucom RS, Kutchan TM, Augustin MM, Wang J, Zhang Y, Tian Z, Yan Z, Wu X, Sun X, Wong GK-S, Leebens-Mack J (2014) Phylotranscriptomic analysis of the origin and early diversification of land plants. Proceedings of the National Academy of Sciences of the United States of America 111: 4859–4868.
  • Widhelm TJ, Egan RS, Bertoletti FR, Asztalos MJ, Kraichak E, Leavitt SD, Lumbsch HT (2016) Picking holes in traditional species delimitations: an integrative taxonomic reassessment of the Parmotrema perforatum group (Parmeliaceae, Ascomycota). Botanical Journal of the Linnean Society 182: 868–884.
  • Wiens JJ, Penkrot TA (2002) Delimiting species using DNA and morphological variation and discordant species limits in spiny lizards (Sceloporus). Systematic Biology 51: 69–91.
  • Wilkinson MJ, Szabo C, Ford CS, Yarom Y, Croxford AE, Camp A, Gooding P (2017) Replacing Sanger with Next Generation Sequencing to improve coverage and quality of reference DNA barcodes for plants. Scientific Reports 7.
  • Wirtz N, Printzen C, Lumbsch HT (2008) The delimitation of Antarctic and bipolar species of neuropogonoid Usnea (Ascomycota, Lecanorales): a cohesion approach of species recognition for the Usnea perpusilla complex. Mycological Research 112: 472–484.
  • Wirtz N, Printzen C, Lumbsch HT (2012) Using haplotype networks, estimation of gene flow and phenotypic characters to understand species delimitation in fungi of a predominantly Antarctic Usnea group (Ascomycota, Parmeliaceae). Organisms Diversity & Evolution 12: 17–37.
  • Wirtz N, Printzen C, Sancho LG, Lumbsch HT (2006) The phylogeny and classification of Neuropogon and Usnea (Parmeliaceae, Ascomycota) revisited. Taxon 55: 367–376.
  • Zhao X, Fernández-Brime S, Wedin M, Locke M, Leavitt SD, Lumbsch HT (2017) Using multi-locus sequence data for addressing species boundaries in commonly accepted lichen-forming fungal species. Organisms Diversity & Evolution 17: 351–363.
  • Zimmer EA, Wen J (2015) Using nuclear gene data for plant phylogenetics: Progress and prospects II. Next-gen approaches. Journal of Systematics and Evolution 53: 371–379.