Metabarcoding of Antarctic Lichens from Areas with Different Deglaciation Times Reveals a High Diversity of Lichen-Associated Communities

Lichens have developed numerous adaptations to optimise their survival under harsh abiotic stress, colonise different substrates, and reach substantial population sizes and high coverage in ice-free Antarctic areas, benefiting from a symbiotic lifestyle. As lichen thalli represent consortia with an unknown number of participants, it is important to know about the accessory organisms and their relationships with various environmental conditions. To this end, we analysed lichen-associated communities from Himantormia lugubris, Placopsis antarctica, P. contortuplicata, and Ramalina terebrata, collected from soils with differing deglaciation times, using a metabarcoding approach. In general, many more Ascomycete taxa are associated with the investigated lichens compared to Basidiomycota. Given our sampling, a consistently higher number of lichen-associated eukaryotes are estimated to be present in areas with deglaciation times of longer than 5000 years compared to more recently deglaciated areas. Thus far, members of Dothideomycetes, Leotiomycetes, and Arthoniomycetes have been restricted to the Placopsis specimens from areas with deglaciation times longer than 5000 years. Striking differences between the associated organisms of R. terebrata and H. lugubris have also been discovered. Thus, a species-specific basidiomycete, Tremella, was revealed for R. terebrata, as was a member of Capnodiales for H. lugubris. Our study provides further understanding of the complex terricolous lichen-associated mycobiome using the metabarcoding approach. It also illustrates the necessity to extend our knowledge of complex lichen symbiosis and further improve the coverage of microbial eukaryotes in DNA barcode libraries, including more extended sampling.


Introduction
Antarctica is the only continent where the dominant vegetation comprises cryptogamic species. Among these, lichens are the most diverse group, encompassing more than 400 described species [1] that mostly grow on soils or rocks. Thus, they are involved in soil formation processes and nutrient cycling and are key elements of Antarctic terrestrial ecosystems [2,3]. Being poikilohydric, lichens are particularly well-adapted to withstanding adverse environmental conditions and reach large population sizes and high coverage in Antarctic ice-free areas [4,5]. Lichens are a fascinating example of obligate fungal symbiosis and are well known to be composed of at least two main organisms: the fungus (mycobiont) and the algal (photobiont) partner [6,7]. In addition, many more or less closely associated additional players, such as accessory fungi, algae, and bacteria, are known to occur [8][9][10]. Thus, lichen thalli represent not only dual or triple symbiosis but also consortia with an unknown number of participants [11] or complex ecosystems [8,12] comprising lichenicolous and endolichenic fungi, as well as bacterial epi-and endobionts with various potential roles in the whole symbiotic relationship [13].
Until the end of the twentieth century, the taxonomy of lichen-forming fungi and their photobionts was based on morphology, anatomy, and chemistry. However, in recent years, molecular tools, bioinformatics, and growing DNA sequence databases have enabled the study of the inter-and intraspecific diversity of lichen-forming fungi, which harbour additional fungi, algae, and bacteria [14][15][16][17]. The complexity of lichen thalli is part of the reason that lichens are still poorly investigated via omics approaches [15,18]. To reduce problems with chimeric contigs, it is preferred that the sequences of lichen symbionts are obtained by cultivating individual partners axenically [16,17,19], but then, the full diversity of symbiosis is not revealed. Moreover, the isolation of lichen symbionts can be challenging because of the high risk of contamination or the time required for spore discharge (mycobiont). Hence, it might take months to one year before enough material can be obtained for further experiments and sequencing [20][21][22], which makes the metabarcoding and metagenomic approaches beneficial.
Next-generation sequencing studies have discovered lichen-specific bacterial communities and various secondary fungal lineages inhabiting lichens [13,15,23,24]. Several investigators have found basidiomycetous yeasts as symptomless endolichenic fungi on and in lichen thalli [10]. A specific group of basidiomycetous yeasts, Cyphobasidiales (Pucciniomycotina, Basidiomycota), was detected as a component of the lichen microbiome [13]. In addition to a phenotypic relationship, the authors associated the abundance of these basidiomycetous yeasts with the production of lichen secondary substances in Bryoria fremontii and B. tortuosa. These two aspects were fundamental to their hypothesis of basidiomycetous yeast being the third mutualistic partner within the lichen complex. However, a subsequent broad metagenomic study of 339 lichen species collected from the Appalachian Mountains in North America failed to detect basidiomycete yeasts in over 97% of the sampled species, leading the authors to question their ubiquity in lichens [25]. Next to secondary fungi, concomitant photobionts may alter the physiological properties of lichens. For example, Casano et al. [26] reported two Trebouxia algae with different physiological performances, which are always present in the lichen thalli of Ramalina farinacea.
The recent application of metabarcoding approaches using high-throughput sequencing (HTS) has highlighted that the microbiome diversity of lichens may be greater than we think (e.g., [27][28][29]) and can help detect cryptic lineages of lichen-associated communities in recently deglaciated soils. Thus, in a recent study, Zhang et al. [30] revealed a high diversity of fungal communities in eleven ice-free Antarctic habitats and elucidated the ecological traits of fungal communities in this unique ice-free area of maritime Antarctica. In addition, Coleine et al. [31] demonstrated the importance of sun exposure on the presence of functional groups of fungi in cryptoendolithic Antarctic communities.
Among the ribosomal cistron regions, the nuclear internal transcribed spacer (nrITS) has been proven to have the most clearly defined barcode gap between inter-and intraspecific variation. For this reason, it was successfully applied in the barcoding of various lichen genera (e.g., [32][33][34][35]) and lichen mycobiomes [36]. Thus, for example, to investigate cryptic diversity in Parmelia, Divakar et al. [34] used a DNA barcoding tool for rapid, accurate sample identification and detection and sampled 202 specimens of Parmelia s.s. worldwide, including Antarctica, using nrITS. In our study, we used next-generation sequencing to target the nrITS2 region. We used a metabarcoding approach to assess and characterise the diversity of the lichen-associated communities in Himantormia lugubris, Placopsis antarctica, P. contortuplicata, and R. terebrata. As the two Placopsis species occur in both recently deglaciated and areas that have been ice-free for longer on the South Shetland Islands, we used the genus Placopsis to investigate whether or not the associated communities differ between areas with different deglaciation times. Furthermore, specimens of this genus from the volcanic, and thus, warmer Deception Island [Olsacher 1956] were included to compare samples from this peculiar habitat with the ones from the other South Shetland Islands.

Taxon Sampling
Sampling covered different successional stages of lichen communities, from the pioneer phase with the crustose lichens P. antarctica D. J. Galloway, R. I. L. Sm. & Quilhot ( Figure 1A), and P. contortuplicata I. M. Lamb to mature lichen communities, including the fruticose lichens R. terebrata Hook. f. & Taylor and H. lugubris (Hue) I. M. Lamb ( Figure 1B,C). The samples collected were identified using standard identification procedures [1,37], with voucher specimens deposited in the herbarium of SNSB-Botanische Staatssammlung München (M). were included to compare samples from this peculiar habitat with the ones from the other South Shetland Islands.

Taxon Sampling
Sampling covered different successional stages of lichen communities, from the pioneer phase with the crustose lichens Placopsis antarctica D. J. Galloway, R. I. L. Sm. & Quilhot ( Figure 1A), and Placopsis contortuplicata I. M. Lamb to mature lichen communities, including the fruticose lichens Ramalina terebrata Hook. f. & Taylor and Himantormia lugubris (Hue) I. M. Lamb ( Figure 1B,C). The samples collected were identified using standard identification procedures [1,37], with voucher specimens deposited in the herbarium of SNSB-Botanische Staatssammlung München (M). This study is based on specimens collected on the South Shetland Islands, Antarctica: Fildes and Potter Peninsula of King George Island, Ardley Island (close to the southwest end of King George Island), Coppermine Peninsula of Roberts Island, Byers Peninsula of Livingston Island, and four localities on Deception Island (Table 1). In total, 32, 10, and 11 thalli of the genus Placopsis were included in our study: 28 thalli of P. antarctica from longdeglaciation-time regions, seven from short-deglaciation-time regions, and ten from Deception Island, as well as four, three, and one thallus of P. contortuplicata from the respective areas. In addition, 14 specimens of Ramalina terebrata (7) and Himantormia lugubris (7) were included. Due to its unique geology, the deglaciation times of the Deception Island cannot be accurately given.   (Table 1). In total, 32, 10, and 11 thalli of the genus Placopsis were included in our study: 28 thalli of P. antarctica from long-deglaciation-time regions, seven from short-deglaciation-time regions, and ten from Deception Island, as well as four, three, and one thallus of P. contortuplicata from the respective areas. In addition, 14 specimens of R. terebrata (7) and H. lugubris (7) were included. Due to its unique geology, the deglaciation times of the Deception Island cannot be accurately given.

Sample Preparation, Lysis, and DNA Extraction
From each specimen, ten branch tips 5 mm long from the fruticose lichens R. terebrata and H. lugubris, or ten areoles of about 0.2 mm 2 from the crustose lichen Placopsis, without any visible contamination of aerophytic organisms, were selected and pooled in 2 mL sterile sample tubes. These tissue samples were subsequently homogenised using a FastPrep96 (MP Biomedicals, Eschwege, Germany), twice for 60 s, with the addition of ceramic beads. Genomic DNA was extracted as described in Beck and Mayr [39], but instead using the Invisorb ® Spin Plant Mini Kit (STRATEC Molecular GmbH, Berlin, Germany) according to the manufacturer's protocol.

DNA Amplification and Metabarcoding
DNA metabarcoding was conducted at the AIM Lab (AIM-Advanced Identification Methods GmbH, Leipzig, Germany). From each sample, 5 µL of extracted total DNA was used for PCR. Plant MyTAQ (Bioline, Luckenwalde, Germany) and high-throughput sequencing (HTS)-adapted mini-barcode primers were used for targeting the Internal Transcribed Sequence 2 (ITS2; primers (gITS7 and ITS4) and amplification following Ihrmark et al. [40]). Amplification success and fragment lengths were verified via gel electrophoresis. Amplified DNA was cleaned using a 1% sodium acetate and 70% ethanol precipitation method and resuspended in 50 µL purified water for each sample before proceeding. Illumina Nextera XT (Illumina Inc., San Diego, CA, USA) indices were ligated to the samples in a second PCR reaction, whereby we applied the same annealing temperature as for the first PCR reaction but with only seven cycles. Ligation success was confirmed via gel electrophoresis. DNA concentrations were measured using a Qubit fluorometer (ThermoFisher Scientific, Waltham, MA, USA) and adjusted to 40 µL pools, each containing equimolar concentrations of 100 ng/µL DNA template. Pools were purified using MagSi-NGSprep Plus beads (Steinbrenner Laborsysteme GmbH, Wiesenbach, Germany) with a final elution volume of 20 µL. High-throughput sequencing (HTS) was performed on an Illumina MiSeq (Illumina Inc., San Diego, CA, USA) using v3 chemistry (2 × 300 paired-end reads, 600 cycles, maximum of 25 million reads).

Barcode Sequence Analysis, Processing, and Amplicon Sequencing Variant (ASV) Identification
Sequence processing was performed with the VSEARCH v2.4.3 suite [41] and cutadapt v1.1439 [42]. Since not all of the sequenced samples yielded reverse reads of high enough quality to enable paired-end merging, only forward reads were utilised. Reads were removed in cases where the number of mismatches was >5, the alignment was shorter than 16 base pairs, or the identity percentage of the alignment was <90%. Forward primers and adaptors were removed with cutadapt [42]. Quality filtering was performed with the fastq_filter program in VSEARCH [41] (fastq_maxee 2, minimum length of 100 bp). Subsequently, sequences were dereplicated using VSEARCH [41] with derep_fulllength, first at the sample level, then concatenated into one fasta file, which was then dereplicated. Chimeric sequences were filtered from the large fasta file with VSEARCH [41] using uchime_denovo. The remaining sequences were clustered into amplicon sequencing variants (ASVs) at 97% identity with cluster_size, and an ASV table was created using VSEARCH [41] with usearch_global. To reduce false positives, a cleaning step was employed, which excluded read counts in the ASV table of less than 0.01% of the total. ASVs were blasted against a custom database downloaded from UNITE (February 2020), NCBI GenBank (January 2022), and BOLD (October 2021), including taxonomy and barcode index number (BIN) information (contained in the BOLD database), using Geneious (v.10.2.5; Biomatters, Auckland, New Zealand) and following the methods described by Moriniere et al. [43]. The resulting csv file, which included the ASV ID, BOLD process ID, BIN, Hit-%-ID value (percentage of overlap similarity, i.e., identical base pairs, of an ASV query sequence with its closest counterpart in the database), length of the top BLAST hit sequence, and phylum, class, order, family, genus, and species information for each detected ASV, was exported from Geneious and combined with the ASV table generated by the bioinformatic pipeline (Supplementary Data S1). The taxonomy of the most complete UNITE BLAST was used in subsequent analyses (Supplementary Data S1).

Sample Pooling, Data Exclusion, and Plausibility Control
Of the 429 ASVs retrieved, 359 ASVs could be identified at the genus level or below via metabarcoding (Supplementary Data Table S1). Obvious contaminants known not to occur in Antarctica were excluded from the analysis. Data on associated eukaryote composition were compiled and analysed, focussing on presence/absence because metabarcoding does not allow for accurate estimations of symbiont quantity [43,44]. Relative read percentages are thus intended to provide rough estimates.

KRONA Metagenomic Visualisation
Interactive HTML charts were created using KronaTools v2.7 [45] (https://github. com/marbl/Krona/wiki/KronaTools, accessed on 20 April 2021). Three sets of charts were created based on ASV table taxonomic annotations from NCBI, UNITE, and BOLD. In each case, a set of charts was created for each sample and for reads summed over all samples. Read counts of different ASVs annotated to the same species were combined in each chart. We used a custom script to extract the ASV counts and associated taxonomic annotations from the final Excel results table. Then, we created the intermediate sample count (.TAX) files for KronaTools using a bash script obtained from https://github.com/ GenomicaMicrob/OTUsamples2krona (accessed on 20 April 2021), where the charts were created using the command (ktImportText [SAMPLE.TAX] -n SAMPLE -o SAMPLE.html). The snapshot function of the HTML Krona charts was used to export figures that showed the relevant organism distribution.

Diversity Estimations
To estimate the expected diversity in the three communities, coverage and sample-sizebased rarefaction/extrapolation curves were computed using iNEXT [46]. Curves with 95% confidence intervals were calculated for Hill numbers q = 0 (species richness), q = 1 (the exponential of Shannon's entropy index), and q = 2 (the inverse of Simpson's concentration index). The size in the rarefaction/extrapolation curve was extrapolated to triple the minimum observed sample size, guided by an estimated asymptote (endpoint = 32), as recommended by Chao et al. [47,48].
The alpha diversity was calculated using the diversity function implemented in the vegan package [49] based on the number and relative abundance of taxa (species, or ASVs) to quantify the diversity in each of the three communities. The outliers were subsequently removed from this analysis (i.e., those whose alpha index was less than 0.4 and higher than 1.6). To test if at least one of the diversity means was different from the others, an analysis of variance (ANOVA) in the vegan package was employed [49]. Subsequently, the pairwise comparisons implemented in the honest significant difference test function (HSD.test) from the agricolae package [50] were used to analyse whether a set of communities' means were different from each other. All graphics were created using the ggplot2 R package [51]. The scripts for coverage, sample-size-based rarefaction/extrapolation curves and alphadiversity calculations are provided in the Supplementary Materials Scripts.

Eukaryotes Associated with the Lichens Himantormia, Placopsis, and Ramalina
Our results reveal a large number of lichen-associated organisms. Most reads expectedly corresponded to the mycobionts of the four lichens investigated, viz. P. antarctica, P. contortuplicata, R. terebrata, and H. lugubris, members of the class Lecanoromycetes. Given that P. contortuplicata is represented by a few specimens, we united the data for both Placopsis species in the Krona plots. The major photobionts of these lichens, Stichococcus in the two Placopsis species and Trebouxia in Ramalina and Himantormia comprised the next most abundant ASVs. Only 4% of the ASVs remained without a determined close match, primarily due to an unidentified Trebouxia symbiont of Himantormia (ASV_6; 2.8%) assigned to this genus in the NCBI search. The Ascomycota and Basidiomycota ASVs from different classes were recovered, namely from Dothideomycetes, Eurotiomycetes, Leotiomycetes, Arthoniomycetes, Sordariomycetes, Lichenomycetes, Saccharomycetes, Agaricomycetes, and Tremellomycetes, respectively. Given this large spectrum of coverage, the presented metabarcoding method is well suited to analysing lichen-associated eukaryotes, with an emphasis on lichen-associated fungi, as they comprise the largest group based on the number of reads and the number of ASVs (Figure 2). Dothideomycetes, Eurotiomycetes, Leotiomycetes, Arthoniomycetes, Sordariomycetes Lichenomycetes, Saccharomycetes, Agaricomycetes, and Tremellomycetes, respectively Given this large spectrum of coverage, the presented metabarcoding method is well suite to analysing lichen-associated eukaryotes, with an emphasis on lichen-associated fung as they comprise the largest group based on the number of reads and the number of ASV (Figure 2).

ASVs Associated with the Fruticose Lichens Ramalina terebrata and Himantormia lugubris
Our results show 168 ASVs associated with the seven specimens of R. terebrata (Figure 3). Most reads were related to the mycobiont (68%). Next to the major photobiont (Trebouxia, comprising 98% of all Chlorophyta reads), the algal genera Stichococcus (1%), as well as Coccomyxa, Desmococcus, Prasioloa, and Symbiochloris, were recorded in very small amounts (Supplementary Table S1). Associated basidiomycete fungi were present in relatively high amounts, in addition to the ascomycete fungi. Four percent of the reads belonged to the genus Tremella, and members of the recently described lichen-associated basidiomycete yeasts, the Cystobasidiomycetes (0.3%; four out of seven specimens; Supplementary Table S1). The associated Ascomycetes comprised only about one-third of the reads assigned to Basidiomycetes (Figure 3) and belonged to the classes Dothideomycetes (mainly Capnodiales) and Sordariomycetes (mainly Hypochreales). Compared to the rich spectrum of lichen-associated organisms in R. terebrata, only about one-third (54) ASVs, were recorded for the seven specimens of H. lugubris. Again, most reads belonged to the mycobiont (55%), and for algae, almost only the major photobiont, Trebouxia, was recovered. The relatively low number of mycobiont reads in H. lugubris compared to R. terebrata (55% versus 68%) may be due to its massive central chondroid axis with thick cell walls [52], resulting in comparatively fewer cells per sample volume, and consequently, lower DNA content. Stichococcus (0.03%), as well as Coccomyxa and Desmococcus (0.02% each), were also recovered, but in tiny amounts (Supplementary Table S1). In addition to the primary symbionts, most reads were assigned to Parateratosphaeria marasasii (Capnodiales, Dothideomycetes; ASV_10), but with 80.4% identity only (Supplementary  Table S1). Thus, this taxon, identified in all seven Himantormia specimens, was sequenced for the first time. The next most abundant ASVs, Gorgomyces and Chalara (Leotiomycetes), were recorded in only two of the seven specimens.
Taken together, the most striking difference between the associated organisms in R. terebrata and H. lugubris was the exclusive occurrence of Basidiomycetes in R. terebrata. On the other hand, taxon ASV_10 from Dothideomycetes, whose nrITS2 sequence was most similar to that of P. marasasii (Capnodiales), was the most frequently associated organism and present in all seven Himantormia specimens. At the same time, it was detected only in tiny amounts in one Ramalina specimen.

ASVs Associated with the Crustose Lichen Genus Placopsis Growing on Soils with Different Deglaciation Times
A total of 129 ASVs were detected in the 53 Placopsis specimens. These specimens were separated into three groups: two differing in deglaciation time and one from the volcanic Deception Island (due to its unique geology, the deglaciation times of the island cannot be accurately given). The specimens collected from soils that have been deglaciated for longer than 5000 years contained 122 of these ASVs, and those collected from Deception Island comprised 75 ASVs. The specimens collected from soils that were deglaciated less than 5000 years ago had 60 ASVs (Supplementary Table S1). These numbers were caused by the higher number of associated fungi growing on soils with longer deglaciation times. Members of the Dothideomycetes, Leotiomycetes, and Arthoniomycetes were restricted to the lichen specimens from soils that have been deglaciated for longer ( Figure 4).
To account for the uneven sample number between the three communities (namely, that with more than 5000 years of deglaciation time, that with less than 5000 years deglaciation time, and that from the volcanic island Deception Island (32, 10, and 11 specimens, respectively), coverage-based and sample-size-based rarefaction/extrapolation curves were calculated using iNEXT [46] for a base sample size of 32, as recommended by Chao et al. [47] ( Figure 5). Consistently, the estimate of ASV diversity (q = 0) was lowest for the samples from areas with a deglaciation time of less than 5000 years, followed by Deception Island, with the highest species number occurring in the regions with deglaciation times of longer than 5000 years. The 95% confidence intervals did not overlap for the samples from areas with different deglaciation times using the sample-size-based approach, indicating a significantly higher species number in the communities with longer deglaciation times. The species number estimates for the rarefied and extrapolated numbers of species were 100, 135, and 145, respectively; thus, the species number of lichen-associated ASVs was estimated to be about 50% higher on soils that have been deglaciated for more than 5000 years compared to those from soils that were deglaciated less than 5000 years. The extrapolation curves of less sampled communities tend to grow over the actual sampling size, which indicates that more comprehensive data are necessary to corroborate these results. To account for the uneven sample number between the three communities (namely, that with more than 5000 years of deglaciation time, that with than 5000 years deglaciation time, and that from the volcanic island Deception Island (32, 10, and 11 specimens, respectively), coverage-based and sample-size-based rarefaction/extrapolation curves were calculated using iNEXT [46] for a base sample size of 32, as recommended by Chao et al. [47] ( Figure 5). Consistently, the estimate of ASV diversity (q = 0) was lowest for the samples species number estimates for the rarefied and extrapolated numbers of species were 1 135, and 145, respectively; thus, the species number of lichen-associated ASVs was es mated to be about 50% higher on soils that have been deglaciated for more than 5000 yea compared to those from soils that were deglaciated less than 5000 years. The extrapolati curves of less sampled communities tend to grow over the actual sampling size, whi indicates that more comprehensive data are necessary to corroborate these results. Figure 5. A comparison of coverage-(left side) and sample-size-based (right side) rarefaction (so line) and extrapolation (dotted line) curves with 95% confidence intervals for Hill numbers q (species richness-upper panel), q = 1 (the exponential of Shannon's entropy index-middle pan and q = 2 (the inverse of Simpson's concentration index-lower panel) for taxa associated with Antarctic lichen genus Placopsis growing on soils with different deglaciation times. Sample-si based extrapolation was done up to the largest sample size, i.e. 32, guided by an estimated asym tote. Solid symbols denote reference samples. Samples collected from the volcanic soil of Decept Island (circles, orange colour), soils deglaciated less than 5000 years ago (triangles, blue colour), a soils deglaciated more than 5000 years ago (squares, purple colour).

Figure 5.
A comparison of coverage-(left side) and sample-size-based (right side) rarefaction (solid line) and extrapolation (dotted line) curves with 95% confidence intervals for Hill numbers q = 0 (species richness-upper panel), q = 1 (the exponential of Shannon's entropy index-middle panel), and q = 2 (the inverse of Simpson's concentration index-lower panel) for taxa associated with the Antarctic lichen genus Placopsis growing on soils with different deglaciation times. Sample-size-based extrapolation was done up to the largest sample size, i.e. 32, guided by an estimated asymptote. Solid symbols denote reference samples. Samples collected from the volcanic soil of Deception Island (circles, orange colour), soils deglaciated less than 5000 years ago (triangles, blue colour), and soils deglaciated more than 5000 years ago (squares, purple colour).
The alpha diversity of lichen-associated communities as expressed by the Shannon Index was also highest for the samples from soils with longer deglaciation times, followed by those collected on Deception Island, and was lowest in samples taken from soils with deglaciation times of less than 5000 years ( Figure 6). Despite these differences, the HSD test showed that the alpha diversity of all three communities did not differ significantly, assigning them to one group. The alpha diversity of lichen-associated communities as expressed by the Shannon Index was also highest for the samples from soils with longer deglaciation times, followed by those collected on Deception Island, and was lowest in samples taken from soils with deglaciation times of less than 5000 years ( Figure 6). Despite these differences, the HSD test showed that the alpha diversity of all three communities did not differ significantly, assigning them to one group. Figure 6. Box-and-whisker plot comparing alpha diversity of lichen-associated communities as expressed by the Shannon Index for taxa associated with the Antarctic lichen genus Placopsis growing on soils with different deglaciation times. Samples collected from the volcanic soil of Deception Island (left), soils deglaciated less than 5000 years ago (middle), and soils deglaciated more than 5000 years ago (right). Box-plot indicates median and lower and upper quartiles, with whiskers drawn within the 1.5 IQR; red dots indicate individual values.
Given the not-yet-complete sampling (see above) we only mention some interesting ASVs and give the currently known distribution pattern only for those taxa present in more than 20% of the samples from a given origin, as a guidance for further studies. Concerning the associated algae (i.e., not the main green algal symbiont), ten different ASVs of the genus Trebouxia were present in quite a low abundance, irrespective of the deglaciation time. Of these, Trebouxia sp. 'A04' (Trebouxiales, Chlorophyta; ASV_197) was unique to Deception Island and was present in three out of the eleven specimens. The Antarctic snow alga Chlorominima collina [53] was found in two lichen specimens from more recently deglaciated soils, while it occurred in only one specimen from soils with a longer deglaciation history. It did not appear in the specimens from Deception Island. In contrast, the Given the not-yet-complete sampling (see above) we only mention some interesting ASVs and give the currently known distribution pattern only for those taxa present in more than 20% of the samples from a given origin, as a guidance for further studies. Concerning the associated algae (i.e., not the main green algal symbiont), ten different ASVs of the genus Trebouxia were present in quite a low abundance, irrespective of the deglaciation time. Of these, Trebouxia sp. 'A04' (Trebouxiales, Chlorophyta; ASV_197) was unique to Deception Island and was present in three out of the eleven specimens. The Antarctic snow alga Chlorominima collina [53] was found in two lichen specimens from more recently deglaciated soils, while it occurred in only one specimen from soils with a longer deglaciation history. It did not appear in the specimens from Deception Island. In contrast, the genera Desmococcus and Symbiochloris were only present in soil specimens with a longer deglaciation time. The ASVs closely related to Stichococcus allas were mainly found in specimens from Deception Island (five out of eleven), while it occurred in only 10% of the specimens from the other South Shetland Islands.
Intriguingly, the Basidiomycete Sistotrema autumnale was found only in specimens from Deception Island (three out of eleven) and in one specimen from the areas that have been deglaciated for longer than 5000 years (Figure 4). This organism is already known from a soil sample in Antarctica (Unite: UDB07059034; [54]) and as a lichenicolous fungus from Physcia aipolia (Unite: MN902070; [55]). Dothideomycetes were exclusively found in areas that have been deglaciated for longer than 5000 years. A taxon related to the genus Bryochiton (Capnodiales; ASV_14) was the most frequent in seven out of the thirtytwo specimens. However, its closest relative, B. perpusillus, had only 89.5% sequence identity and, thus, clearly represented a distinct species based on sequence identity. The Arthoniomycetes were exclusively found in areas that have been deglaciated for over 5000 years. Of these, the genus Bryostigma is known to contain lichenicolous species [56].
A possible member of Herpotrichiellaceae (Chaetothyriales, Eurotiomycetes; ASV_207) was found to be unique to areas deglaciated less than 5000 years ago and was present in two out of the ten specimens, with 83.3% identity in the NCBI BLAST search. The next closest hit showed 82.4% identity to Dothideomycetes, but this is most likely due to a misidentification, as the next similar GenBank hits referred to Chaetothyriales as well. Thus, ASV_207 is most likely a member of Chaetothyriales, with no close relative sequenced so far. Moreover, the specimens from Deception Island contained a unique ASV related to a member of the genus Tetracladium sp. (Helotiales, Leotiomycetes; ASV_119). It was found to be present in three out of the eleven specimens.

Conclusions
Our study is the first to employ a DNA metabarcoding approach to lichens growing on soils with differing deglaciation times collected in Antarctica. This is a baseline study, that revealed the first distribution patterns of lichen-associated eukaryotes. However, given the small sample size of communities with short deglaciation times, the results still need further corroboration by extending this sample size. Nevertheless, some patterns already emerge: Overall, there are many more Ascomycete taxa associated with lichens compared to Basidiomycota. Significantly more lichen-associated eukaryotes are estimated to be present in areas with deglaciation times of longer than 5000 years compared to more recently deglaciated areas ( Figure 5). The occurrence of several lichen-associated organisms seems to be related to the deglaciation time of their substrate. An intriguing example is a member of the genus Bryochiton (Capnodiales), which was found exclusively in 20% of the Placopsis growing in ice-free areas with deglaciation times of longer than 5000 years, while a member of the Herpotrichiellaceae (Chaetothyriales) was unique to areas deglaciated less than 5000 years ago. As for Basidiomycetes, lichenicolous Tremella and Cystobasidiomycetes seem restricted to R. terebrata. Some taxa are found almost exclusively in one of the lichens studied and are present in all specimens analysed. Interesting examples are ASV_16, which is related to the Tremella genus and occurs almost exclusively in R. terebrata, and ASV_10, which is distantly related to P. marasasii (Dothideomycetes) and occurs almost exclusively in H. lugubris.
Although ten different parts of the lichen specimens were pooled to account for heterogeneity (viz., thallus tips and areoles, respectively), considerable variation in the associated organisms can still be observed. Thus, the Cystobasidiomycete representatives were only detected in about half of the specimens, possibly due to the patchy distribution or varying quantities of these associated eukaryotes. These results are in line with those obtained by Dal Grande et al. [57] and Rolshausen et al. [58], where the importance of environmental factors and host identity on the community structure of green algal symbionts and lichen holobionts was shown. Taken together, the environment seems to be the most crucial factor for the formation of lichen holobionts.
The ongoing growth of DNA barcode libraries, such as NCBI and UNITE, with the regular addition of new reference barcodes, will increase the capability of this metabarcoding approach in the future. Although many barcoding efforts have been undertaken in recent years (see Introduction), microbial eukaryotes are still not well covered. This fact can be exemplified by ASV_207, most closely related to Herpotrichiellaceae (Chaetothyriales, Eurotiomycetes), but with only 83.8% identity. In addition, caution still needs to be taken when interpreting results from one database alone. This can be exemplified by ASV_207 being assigned to different classes of Ascomycetes: to Dothideomycetes in NCBI BLAST (identity 84.2%) and Eurotiomycetes in the UNITE database (identity 83.8%). However, the BLAST result is most likely erroneneous, as the next similar GenBank hits also refer to Chaetothyriales. Thus, ASV_207 seems to be a member of Chaetothyriales, whose close close relative has not been sequenced so far. In the context of climate change, not only are temperatures expected to rise but also large areas of the previously ice-covered land surface are expected to be exposed, estimated to be up to 25% of the current ice-free area [59]. Connectivity between the previously relatively isolated ice-free patches will be increased, leading to an easier exchange of organisms, with an anticipated effect on biodiversity [60]. Consequently, it is essential to have baseline information on terrestrial biodiversity, as shown in this study. Yet, a significant fraction of this biodiversity is still not well analysed, if at all. This is, however, very important, as the successional dynamics of such communities may differ, as has recently been shown for bacterial and fungal communities in recently deglaciated soils of the maritime Antarctic [61]. In addition, it has been demonstrated that lichen-associated eukaryotes may differ between seasons [62]. In Antarctica, the seasons are significantly pronounced and considerably differ, e.g., increased snow cover and sub-zero temperatures. Still, it has been shown that lichens can perform photosynthesis below the snow, and even in a frozen state [63,64]. Thus, seasonal dynamics of lichen-associated organisms may occur even under the harsh conditions in Antarctica, which should be studied in detail in the future.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes14051019/s1, The details involved in the article can be found in Supplementary Data Table S1 and Supplementary Materials Scripts.
Author Contributions: A.B. designed the study, collected the specimens, guided the molecular work, prepared the photographs, participated in data analyses, and wrote the manuscript; J.G. performed part of the molecular work, conducted the diversity estimations, and assisted in manuscript writing and revision; A.C.-K. provided logistic support and organised travel to Antarctica, participated in specimen collection, and revised the MS. All authors have read and agreed to the published version of the manuscript.
Funding: The molecular part of this work was supported by the Staatliche Naturwissenschaftliche Sammlungen Bayerns (SNSBinnovativ to AB). Logistic support was provided by the National Fund for Scientific and Technological Development (ANID-FONDECYT 1118745 to AC-K) and the Instituto Antártico Chileno (INACH RT-2716 to AC-K). Financial support was also provided by the DFG to AB for attending the SPP 1158 Topic Workshop "Polar Genomics", held from the 16