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.

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 1–4) 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.