Population genomic analyses of RAD sequences resolves the phylogenetic relationship of the lichen-forming fungal species Usneaantarctica and Usneaaurantiacoatra.

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.

Introduction

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 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 : and , 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 in (Lamb 1964; Walker 1985) or as a distinct genus (Krog 1976; 1982; Lamb 1939). Molecular studies confirmed (including ) as a monophyletic genus within (Crespo et al. 2007); however, the relationship of and remained ambiguous. A two-marker DNA analysis of elevated to a generic rank (Articus 2004), but subsequent studies provided evidence that is polyphyletic with a core group nested within (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 and (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 and .

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 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 and sexual 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.

Methods

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 and 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 and 47 specimens were identified as based on their phenotypical characters (Walker 1985). As a reference sequence to filter for lichen-fungal loci of and during the RADseq processing, we sequenced a specimen of 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 and , but only the central axis in 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 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 . 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 and 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 and 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 (https://github.com/dereneaton/ipyrad/blob/master/docs/index.rst) 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 (http://tree.bio.ed.ac.uk/software/figtree/).

Analysis of population structure

To calculate differences in the population structure between and , 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 and 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 finerad_input.py included in fineRADstructure-tools (https://github.com/edgardomortiz/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.

Reproducibility

The reference sequence and all scripts that were used in this study are available online (https://github.com/felixgrewe/Usnea). All RAD sequences were deposited in the NCBI Sequence Read Archive (SRA) with accession number PRJNA505526.

Results

Reference assembly and RADseq results

We assembled a draft reference genome of to filter for fungal RAD loci from and . The Illumina NextSeq sequencing of the whole 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 had a total size of 24.1 Mbp and an estimated completeness of 72.18%.

We included 105 specimens of and 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 (EL0059, EL0281) and two samples of (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 and (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 clade showed higher internal support values than the clade. Within the clade, three sister relationships of (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 and 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. ADAPC plot of the densities of (blue) and (green) on the first retained discriminant function B Bar plot of group membership probabilities.

The separation of and 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, showed a higher estimated co-ancestry than (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 than (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 from Potter Peninsula, King George Island (EL0022, EL0034 and EL0051) and for six specimens of from Fildes Peninsula, King George Island (EL0444, EL0435, EL0445, EL450, EL0455 and EL0453). Moreover, two 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 specimens were close sister taxa.

Figure 4.

Clustered fineRADstructure co-ancestry matrix. A Full dataset including collected on the Antarctic Peninsula in addition to and collected on King George Island and Elephant Island B Reduced dataset with all and collected on King George Island and Elephant Island. Two major clades are corresponding to the two species (top-left) and (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.

Discussion

In this study, we used RAD sequencing for evaluating the delimitation of two predominantly Antarctic 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 and 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 and , 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 . 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 and 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 , compared to the asexually-reproducing 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 than 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 to include a much wider geographical range (Antarctic Peninsula in addition to King George and Elephant Island) compared to the sampling of (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 from the Antarctic Peninsula for a direct comparison, it should be noted that covers a wider geographical range than (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 . The exclusively sexual reproduces via the dispersal of fungal spores which are required to meet with an appropriate photobiont after germination. The asexual 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 . Finally, it was predicted that a small number of sexual individuals per generation — and rarely can be found with apothecia — is sufficient to make an apparently asexual population highly variable (Bengtsson 2003).

Despite the lower co-ancestry of compared to , we detected three pairs of very close relatives with high co-ancestry of (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 expressed similarly close relationships. However, we could detect substructure for a group of six individuals of collected at the same location and three specimens of 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 may, in future, provide a comprehensive picture of population structure and diversification.

Conclusion

We successfully used RADseq for phylogenetic and population genomic studies on two species of the lichen-fungal genus . Phylogenetic inference using RAD data clearly delimited the species and 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 than in the sexual , potentially derived from a wider geographical distribution of 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.