Unravelling the Symbiotic Microalgal Diversity in Buellia zoharyi (Lichenized Ascomycota) from the Iberian Peninsula and Balearic Islands Using DNA Metabarcoding

Abstract

Buellia zoharyi is a crustose placodioid lichen, usually occurring on biocrusts of semiarid ecosystems in circum-Mediterranean/Macaronesian areas. In previous work, we found that this lichenized fungus was flexible in its phycobiont choice in the Canary Islands. Here we test whether geography and habitat influence phycobiont diversity in populations of this lichen from the Iberian Peninsula and Balearic Islands using Sanger and high throughput sequencing (HTS). Additionally, three thallus section categories (central, middle and periphery) were analyzed to explore diversity of microalgal communities in each part. We found that B. zoharyi populations hosted at least three different Trebouxia spp., and this lichen can associate with distinct phycobiont strains in different habitats and geographic regions. This study also revealed that the Trebouxia composition of this lichen showed significant differences when comparing the Iberian Peninsula with the Balearics thalli. No support for differences in microalgal communities was found among thallus sections; however, several thalli showed different predominant Trebouxia spp. at each section. This result corroborate that thallus parts selected for DNA extraction in metabarcoding analyses are key to not bias the total phycobiont diversity detected. This study highlights that inclusion of HTS analysis is crucial to understand lichen symbiotic microalgal diversity.

1. Introduction

Lichen holobionts are complex systems individualized from cyclical symbiotic associations [1], which represent one of the most paradigmatic examples of ecologically obligate mutualistic interactions between at least two organismal groups [2]: a fungus (the mycobiont) and a population of photosynthetic green or blue–green algae (the photobionts). The definition of lichens has lately been expanded to include other fungi [3,4], basidiomycete yeasts [5,6,7,8], non-photosynthetic bacteria [9,10,11,12,13] and additional populations of green microalgae [14,15,16,17,18].

Focusing only on the two main lichen partners, the myco- and phycobionts, Muggia et al. [19] proposed that these associations vary depending on the selectivity that the two partners show towards each other. Although the concept of partner selectivity is usually referred to from the mycobiont towards the phycobionts [20,21], selectivity and specificity are currently regarded as multidirectional [19]. High flexibility in the photobiont choice is considered as a strategy to survive under variable selective pressures, widening the ecological niche and therefore enabling the occurrence of lichens even in nutrient-poor and climatically-harsh environments [22,23,24]. Furthermore, the association with a wide range of symbionts, including locally adapted photobiont strains, has been hypothesized as a driver of lichenized fungal evolution [25,26]. Although the photobiont adaptive hypothesis has only been occasionally experimentally tested in lichens [24,27], several studies have shown a strong correlation between the photobiont genotype and the geographical and ecological distribution of the lichen [23,28,29,30,31,32].

Intrathalline microalgal diversity, i.e., multiple photobiont species coexisting within a single lichen thallus, has previously been observed in a number of lichen symbioses [14,15,26,27,33]. For the vast majority of lichens, however, photobiont diversity has been only investigated by Sanger sequencing. Indeed, traditional Sanger sequencing may fail to detect the less represented photobionts and lead to an underestimation of the whole photobiont diversity, while this drawback is overcome by the implementation of high throughput sequencing (HTS) [16,34,35,36]. Currently, an improved assessment of organismal diversity in lichens has been possible by applying environmental DNA (eDNA) metabarcoding analyses, based on high throughput sequencing of entire holobiome communities within a thallus [4,37,38,39,40]. Recently lichen-associated phycobiont communities, have been explored by eDNA metabarcoding approaches within the Mediterranean basin [18,29,34], and in other lichens with diverse geographic origins and growth forms [33,41,42]. The target in most of these studies has been lichens with a foliaceous or fruticose growth form, whereas the number of works investigating coexistence of microalgae in crustose lichen thalli has been comparatively few [14,35,43]. Only Moya et al. [17] have analyzed the microalgal diversity by Sanger sequencing and 454-pyrosequencing in crustose and lichenicolous lichens growing on gypsum soil crusts (biocrusts). In this study, the intrathalline coexistence of diverse microalgal lineages was confirmed in all the analyzed samples, and it was hypothesized that the range of associated phycobionts could be influenced by thallus morphology.

Some areas distributed throughout the central Iberian Peninsula show characteristic gypsum outcrops colonized by a predominant group of crustose lichens. Buellia zoharyi Galun is a circum-Mediterranean/Macaronesian lichen forming crustose placodioid thalli, which usually occurs on biocrusts in semiarid areas (Figure 1a,b) [44,45,46]. Specifically, this species predominantly grows on gypsum soils (Figure 1b) [47,48,49], occasionally on basaltic lava flows in the Canary Islands [50] and in calcareous substrate in the Balearic Islands (Figure 1c). Chiva et al. [51] analyzed mycobiont diversity in B. zoharyi in the Mediterranean region including populations in the Canary and Balearic Islands. Recently, Molins et al. [50] applied a multidisciplinary approach to describe microalgae diversity from B. zoharyi, covering the entire distribution range of the species in the Canary Islands. Different Trebouxia spp. were detected as the primary phycobiont in each island (Trebouxia cretacea—Fuerteventura, Trebouxia asymmetrica—Lanzarote and Trebouxia sp. arnoldoi—Tenerife), which indicated that B. zoharyi is flexible regarding photobiont choice and that it depended on geography. Moreover, coexistence of various Trebouxia spp. within a thallus were detected by using specific primers-PCR. Thus, the aim of this study is to analyze the microalgal diversity, both by Sanger and HTS, in B. zoharyi from the remaining known distribution, i.e., the Iberian Peninsula and Balearic Islands, and to test whether geography is further influencing microalgae diversity in this lichen.

2. Materials and Methods

2.1. Sampling and Experimental Design

In this study, 153 Buellia zoharyi thalli collected from 18 geographical populations were analyzed (Table S1). We included fresh (n = 120) and herbarium samples (n = 33). Fresh specimens were air-dried for one day after sampling and then stored at −20 °C. Lichen thalli were examined under a stereo microscope to remove surface contamination (i.e., sand, mosses and fragments of other lichen species) or infected areas by lichenicolous fungi. Selected fragments from different parts of each thallus were randomly excised and pooled together (termed MIX, Figure 1d). These fragments were cleaned with a sterile razor blade under a stereo microscope and were superficially sterilized following Arnold et al. [52] to ensure that sequenced microalgae have an intrathalline origin.

2.2. DNA Extraction, Primary Photobiont PCR Amplification and Sanger Sequencing

Total genomic DNA was isolated and purified using the DNeasy Plant Mini kit (Qiagen, Hilden, 121 Germany) following the manufacturer’s instructions. Two algal loci were amplified to confirm the identity of the primary photobiont in the 153 lichen thalli: a region of the chloroplast LSU rDNA gene, using the algal specific primer pair 23SU1 and 23SU2 [53], and the nuclear ribosomal internal transcribed spacer (nrITS), using primers nr-SSU-1780 [54] and ITS4 [55]. PCR reactions were performed following Moya et al. [56]. PCR products were visualized on 2% agarose gels and purified using the Gel Band Purification Kit (GE Healthcare Life Science, Buckinghamshire, UK). Amplified PCR products were sequenced with an ABI 3730XL using the BigDye Terminator v.3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA). Sanger sequences were visualized and manually evaluated with Chromas v.2.6.6.0 (http://technelysium.com.au/wp/chromas/ (accessed on 20 May 2021)). Taxonomic assignment of Sanger-sequenced strains in the genus Trebouxia is often problematic and BLAST searches have shown high failure rates. Therefore, we double-checked the results obtained with an online BLAST search (https://blast.ncbi.nlm.nih.gov/Blast.cgi (accessed on 20 May 2021)) against the GenBank sequence database with the entry query filter “UTEX OR SAG”, and by conducting phylogenetic analyses using an accurate Trebouxia database (data not shown, Table S2).

2.3. Metabarcoding Assay

A total of 36 DNAs of fresh samples collected between 2012-2015, treated as MIX were selected from the Iberian Peninsula (ES-IP1, ES-IP2, ES-IP3, ES-IP4 and ES-IP5) and Balearic Islands’ (ES-BI1, ES-BI2 and ES-BI3) populations to be further analyzed using Illumina HTS. As a population of B. zoharyi is usually made up of different sized thalli, we selected eight thalli (one per locality) with a width of 4–6 cm. Adjusting the method of Moya et al. [34] to potentially corroborate the hypothesis of differential distribution of the microalgae throughout the specimens, each thallus was further divided into three parts: central (C, corresponding to the innermost area), middle (M, middle area) and periphery (P, the most external area of the thallus) (n = 24, Figure 1e). These C, M and P samples were cleaned and sterilized as previously explained. Total genomic DNA was isolated and purified as previously mentioned.

Chlorophyta microalgal communities associated with the thalli (n = 60, MIX = 36 + C-M-P= 24) were assayed using Illumina HTS of the first section of the ITS1 of the rRNA operon, proposed as a universal barcode across eukaryotic kingdoms [57]. Amplicons for Illumina MiSeq sequencing were generated by nested PCR: in the first PCR the forward nr-SSU-1780 and the reverse ITS4 primers were used and 27 amplification cycles were run; in the second PCR round, three replicates were generated using algal specific primers nr-SSU-1780/5.8S [18,34] modified with Illumina overhang adaptors (forward overhang: 5′-TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCTGCGGAAGGATCATTGATTC; reverse overhang: 5′-GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGATCCGTTGTTGAGAGTTGTC) and 22 amplification cycles were run. Nested PCRs (25 μL) contained 2.5 μL of 10X buffer, 0.4 μM nr-SSU-1780/5.8S 2R primers, 0.2 mM dNTPs and 0.6 u/μL of ExTaq (Takara, Shiga, Japan) and sterile Milli-Q water was used to bring to volume. The PCR conditions were 1 cycle of 95 °C for 2 min, (27/21) number of cycles of 94 °C for 30 s, 56 °C for 45 s and 72 °C for 1 min, and a final extension of 72 °C for 5 min. These three replicates were then pooled together. PCR reactions were performed as described in Moya et al. [34] and were purified using AMPure XP beads (Beckman Coulter, Brea CA, USA). Indexing PCR and the addition of Nextera sequence adapters were performed using Nextera XT Index kit (Illumina Inc., San Diego, CA, USA) following the protocol for Illumina L library preparation. Finally, a second purification step was carried out using AMPure XP beads (Beckman Coulter, Brea, CA, USA) and libraries were quantified and pooled together. The libraries were sequenced on an Illumina MiSeq platform using the MiSeq Reagent Kit v.3 (paired end 2 × 300 bp) (Genomics Core Facility at the University of Valencia, Spain). Negative control extraction and PCRs were performed and included in the run to check potential contaminants.

2.4. Bioinformatic Analyses

Quality control analysis of the Illumina MiSeq paired-end reads was performed using FastQC v.0.11.8. Raw reads were processed using Quantitative Insights into Microbial Ecology 2 (QIIME2 v.2018.11, [58]). Demultiplexed paired-end sequence reads were preprocessed using DADA2 [59], a package integrated into QIIME 2 that accounts for quality filtering, denoising, joining paired ends and removal of chimeric sequences. The first 20 bp were trimmed from forward and reverse reads before merging to remove adaptors. DADA2 denoises datasets to resolve sequence variants that differ by as little as one nucleotide and are referred to as amplicon sequence variants (ASVs, [59,60]). We performed further data filtering by removing any ASV present with less than 1000 reads in the entire dataset. The final ASVs table contained a total of 255 ASVs (Table S3). All the sequences obtained in this study have been submitted to NCBI/SRA under the bioproject identification number PRJNA720561.

2.5. Phylogenetic Assignment of Trebouxia ASVs

The identity of Trebouxia ASVs was further checked by phylogenetic assignment. The sequence alignment was prepared including: (i) 176 Trebouxia ASVs sequences, (ii) selected sequences published by Leavitt et al. [28] and Moya et al. [34] and (iii) a selection of Trebouxia type species available from the Culture Collection of Algae at Göttingen University (SAG), from the Culture Collection of Algae at the University of Texas (UTEX) and reference sequences of Trebouxia sp. TR9 (FJ418565). A multiple alignment was built using MAFFT v.7.0 [61] using default parameters. The substitution model GTR + G was the most appropriate for the nrITS region according to the Akaike information criterion (AIC) using JModelTest v.2.1.4 [62]. The phylogenetic trees were inferred by Bayesian inference (BI) and maximum likelihood (ML) approaches. ML analysis was implemented in RAxML 8 [63] using the GTRGAMMA substitution model. Bootstrap support (BS) was calculated based on 1000 pseudo replicates [64]. BI was carried out in MrBAYES v.3.2 [65]. Settings included two parallel runs with six chains over 20 million generations, starting with a random tree and sampling after every 200th step. We discarded the first 25% of the data as burn-in, and the corresponding posterior probabilities (PPs) were calculated from the remaining trees. The phylogenetic tree was visualized in FIGTREE v.1.4.2 (http://tree.bio.ed.ac.uk/software/Figuretree/ (accessed on 20 May 2021)) [66]. All analyses were run at the CIPRES Science Gateway v.3.3 webportal [67]. ASVs corresponding to other green microalgae genera such as Asterochloris, Coccomyxa, Coenochloris, Diplosphaera, Elliptochloris, Myrmecia, Pseudostichococcus and Vulcanochloris and uncultured Trebouxiophyceae were identified by BLAST searches.

2.6. Microalgal Diversity at Community Scale

A ß-diversity analysis was conducted to illustrate differences in the composition of Trebouxia microalgae across sample categories. First, we considered geography, and grouped the sampling localities into two regions, i.e., the Iberian Peninsula and the Balearic Islands, as these are separated by a stretch of the Mediterranean Sea and display slightly different thermotypic and ombrotypic ranges [68]. Second, we checked whether the central, middle and peripheral (C, M and P) sections of thalli hosted different Trebouxia communities. Non-metric multidimensional scaling (NMDS) ordinations were computed on the basis of Bray–Curtis dissimilarities calculated with an ASVs table normalized by the cumulative sum scaling method (CSS, [69]). To generate statistical support for the observed differences, we ran the analysis of similarities (ANOSIM, [70]) test with 999 permutations. Subsequently, we calculated the number of observed ASVs and Shannon diversity indices (α-diversity) for the above-mentioned geographic regions, and for each sampling locality within, using a non-rarefied dataset. Differences were statistically evaluated with the nonparametric Mann–Whitney and Kruskal–Wallis tests. All analyses were conducted in the MicrobiomeAnalyst [71] and all the tests regarded p-values below 0.05 as significant.

3. Results

3.1. Primary Trebouxia Detected by Sanger Sequencing (nrITS + LSU rDNA)

Three Trebouxia microalgae species were collectively detected in the 153 lichen thalli examined (Figure 2a). Based on phylogenetic reconstructions Trebouxia spp. tend to form five well-supported monophyletic clades-clade ‘A’, clade ‘C’, clade ‘I’, clade ‘S’ and clade ‘D’. According to this clade code introduced for Trebouxia by Leavitt et al. [28], and subsequently applied by Moya et al. [34], Muggia et al. [72] and Xu et al. [73], these Trebouxia species belong to clade ‘A’ (arboricola/gigantea), and most are precisely assigned to the taxa T. cretacea, T. asymmetrica and Trebouxia sp. OTUA25. Trebouxia cretacea occurred as the unique primary phycobiont in ten populations (ES-IP1, ES-IP5, ES-IP-6, MA1-3, GR, CY, IL and AZ). Trebouxia asymmetrica was exclusively detected in ES-BI1. In ES-IP4, eighteen thalli showed T. cretacea and two T. asymmetrica. IR, LY, ES-IP3 and ES-BI3 combined specimens with T. cretacea and others with Trebouxia sp. OTUA25. Higher flexibility in the Trebouxia choice was observed in ES-IP2 and ES-BI2, where T. cretacea, T. asymmetrica and Trebouxia sp. OTUA25 were detected as primary phycobionts.

3.2. Microalgal Metabarcoding

A total of 15,760,017 ITS1 raw reads were generated, of which 8,632,592 passed the demultiplexing step and quality filter. After excluding fungal sequences and filtering the ASVs represented by less than 1000 reads, the dataset was finally built of 255 ASVs (5,814,389 reads) (Table S3), with a minimum of 5686 (in sample ES-IP5_3_MIX) and a maximum of 261,120 (in ES-IB3_1_MIX) reads per sample. These ASVs were classified by BLAST searches into Trebouxia spp. (89.5% reads), and other green microalgae genera (10.5%, Table S3). In 52 samples, more than 70% of the reads belonged to members of Trebouxia.

3.3. Taxonomic Assignment of the Microalgal Taxa

A total of 21 Trebouxia lineages or taxonomic units were detected in the 153 samples using a combination of BLAST identity searches and phylogenetic analyses (Figure S1). The aligned Trebouxia ITS1-5.8S fragment was 220 bp in length. BI and ML phylogenetic trees were topologically congruent and also consistent with previous studies [19,28,41,72,73,74,75]. According to the Trebouxia clade nomenclature, ASVs of 12 Trebouxia species belonged to clade ‘A’, more particularly to T. asymmetrica, T. cretacea, T. jamesii, Trebouxia sp. TR9, Trebouxia sp. OTUA18 (KR913237), Trebouxia sp. OTUA12 (KR913203), Trebouxia sp. OTUA25 (KR913259), Trebouxia sp. OTUA55 (KT819993), Trebouxia sp. OTUA33 (FJ626728), Trebouxia sp. OTUA46 (FJ792797), Trebouxia sp. OTUA52 (MN622994) and Trebouxia sp. aff. OTUA56 (FJ92797). Of the remaining ASVs, three were placed in clade ‘I’: Trebouxia potteri, Trebouxia sp. OTUI52 (MH035938) and Trebouxia sp. aff. OTUI03 (KR913991); six in clade ‘S’: Trebouxia sp. OTUS52 (MH035971), Trebouxia sp. OTUS08 (GQ375345), Trebouxia sp. OTUS50 (RRJNA330765), Trebouxia sp. aff. OTUS13 (FJ626727), Trebouxia sp. aff. OTUS02 (GQ375315) and Trebouxia sp. aff. OTUS08 (GQ375345). Trebouxia cretacea and Trebouxia sp. OTUS52 showed the highest number of reads, represented by 61 and 47 ASVs, respectively. Trebouxia sp. OTUI52, T. asymmetrica, Trebouxia sp. OTUA52 and Trebouxia sp. OTUI03 were represented by 17, 11, 8 and 7 ASVs, respectively. The remaining Trebouxia spp. were represented only by one, two or three ASVs (Figure S1).

In the group “other green microalgae”, a total of eight genera were included, i.e., Asterochloris, Coccomyxa, Coenochloris, Diplosphaera, Elliptochloris, Myrmecia, Pseudostichococcus and Vulcanochloris, and a category named “uncultured Trebouxiophyceae” (Tables S3 and S4). Asterochloris, Coccomyxa, Coenochloris and Diplosphaera were detected in all populations; Elliptochloris was found in six, and the remaining four microalgae occurred in one or three populations (Table S4).

3.4. Relative abundance of Trebouxia by Locality and Region

The composition of Trebouxia species varied between locations and regions (Figure 2b,c). Trebouxia spp. under 60,000 reads, representing less than 10% of reads in all the localities, were grouped and named as “other Trebouxia spp.”. Trebouxia cretacea was predominant in thalli of ES-IP1, ES-IP3, ES-IP5 and ES-BI3 (more than 60% of reads), whereas six other Trebouxia spp. were detected with a relative abundance ranging from 10% to 30%. T. asymmetrica was predominant in ES-BI1 (80% of reads) and five other Trebouxia spp. represented 20% of reads. On the contrary, thalli from ES-IP2 and ES-IP4 showed the balanced co-presence of T. cretacea and T. asymmetrica and four other Trebouxia spp. represented 20–30% of the reads. Only 35% of reads in ES-BI2 belonged to T. cretacea and 45% of reads to five other Trebouxia spp.

In detail, all thalli treated as MIX analyzed in ES-IP1 and ES-BI3 showed a pattern in which more than 40% of the reads corresponded to T. cretacea (Figure S2). The primary photobiont identified by Sanger in these thalli fitted with the most abundant OTU sequenced by Illumina (T. cretacea). Illumina of some thalli from ES: IP3-IP4-IP5 revealed a heterogeneous pattern, and Sanger identified either T. asymmetrica or T. cretacea as the primary photobiont. In ES-IP2, the primary phycobiont identified by Sanger matched with the most abundant OTU sequenced by Illumina in three thalli (two T. asymmetrica and one T. cretacea), but two other thalli revealed the above-mentioned heterogeneous scenario. Three and two thalli from ES:BI1-BI2, respectively, showed the same most predominant phycobiont both by Sanger and Illumina T. asymmetrica/T cretacea; one showed the heterogenous pattern, and in the other one Sanger and Illumina did not identify the same Trebouxia spp.

3.5. Trebouxia Microalgal Community Differences between Regions and Thallus Sections

The NMDS ordination constructed on the basis of Bray–Curtis dissimilarities showed a weak differentiation in Trebouxia communities in lichen thalli collected in the Iberian Peninsula and the Balearic Islands (stress: 0.23428, Figure 3a). The ANOSIM provided statistical support for the observed differences (R: 0.17175), with a p-value < 0.015. Only ES-BI1_5 and ES-IP_1 maintained T. asymmetrica and T. cretacea respectively, as the most predominant phycobiont (more than 60% of reads) in each part. The remaining six thalli showed different predominant Trebouxia spp. and minor microalgae composition between each section (Figure S3). However, no support was found for differences in microalgal communities among the three thallus section categories (i.e., central, middle and periphery), which was reflected in an NMDS ordination with a high stress value (0.21282; Figures S3 and S4) and a corresponding ANOSIM R value of −0.042039 (p-value < 0.77).

Regarding α-diversity analyses, values for the number of observed ASVs and the Shannon diversity indices tended to be higher for the Iberian Peninsula than the Balearic Islands (Figure 3b,c), but the observed differences did not receive statistical support (Mann–Whitney’s U values of 161 and 102, respectively, and p-values > 0.05). When sampling localities were considered, the highest value ranges for the observed ASVs and Shannon indices were shown by ES-IP1, 2 and 4. Differences among localities were statistically supported for the number of observed ASVs (Kruskal–Wallis’ K of 17.256; p-value: 0.015816) but not for the Shannon index (Kruskal–Wallis’ K of 9.4707; p-value: 0.2206). Similar value ranges for both α-diversity indices and statistical tests results were obtained with a rarefied ASVs data matrix (data not shown).

4. Discussion

Numerous mycobionts have been shown to associate with identical species of Trebouxia, while others exhibit higher phycobiont flexibility [28,42,76,77,78,79,80]. Buellia zoharyi was described by Molins et al. [50] to be flexible in its photobiont choice, as this lichen-forming fungus associated with three Trebouxia species (T. cretacea, T. asymmetrica and Trebouxia sp. arnoldoi) in a restricted geographic area, the Canary Islands. The present study completed this sampling, covering the entire distribution range of this lichen, and reinforced the idea that this fungus can associate with different phycobiont strains in different geographic regions. Trebouxia cretacea associated with the fungus across the whole studied area as the primary photobiont. However, in the Balearic island of Mallorca (ES-BI1), the fungus switched from Trebouxia cretacea and associated instead with T. asymmetrica, and in Formentera (ES-BI2) B. zoharyi hosted at least three different Trebouxia spp. Flexibility in the phycobiont choice extends the ecological amplitude of lichens; therefore, it may constitute an important adaptive strategy for the colonization of extreme habitats [22,23,24].

This study revealed that in the case of B. zoharyi the Trebouxia composition showed significant differences when comparing the Iberian Peninsula with the Balearic Islands. Molins et al. [18] explored microalgal diversity in lichen thalli of the model species Ramalina farinacea from different habitats. Similarly, Trebouxia species diversity and composition inside the thalli of R. farinacea was observed to correlate with the geographic origins of the samples: continental Mediterranean vs insular Macaronesian. Both results revealed that intrathalline microalgal diversity strongly depends on the geographic area and the habitat type.

The dispersal of B. zoharyi over medium to long distances can be accomplished by both meiotic ascospores or thallus fragments [48,81]. Both types of reproductive strategy are suitable for long-distance dispersal, even across the Mediterranean Sea, based on evidence of other lichens showing disjunctive populations [82]. Photobiont diversity can be shaped by the reproductive and dispersal strategies of the mycobiont [83,84]. Molins et al. [50] suggested that geological events occurred during Pleistocene glaciations (when connectivity among the Canary Islands, Africa and the Iberian Peninsula was increased) could have influenced the Trebouxia patterns observed in the Canary Islands. The colonization of the Canary Islands by B. zoharyi must have originated from the Moroccan coast to Fuerteventura, evidenced as the maintenance of the symbiont pattern (T. cretacea), and then B. zoharyi must have colonized the Canary Islands subsequently and made use of ecologically adapted phycobionts (T. asymmetrica and Trebouxia sp. arnoldoi) in those islands. Similar explanation based on geological events could be suggested in the case of the Balearic Islands, however to corroborate this hypothesis more specimens and locations (i.e., Menorca) need to be analyzed. Differences in the growth substrate has been suggested by other authors [75,79,85] as a factor that may influence phycobiont diversity. Related with this idea, it should be pointed out that the Balearic Islands, Libya, Israel and Iran hold calcareous biocrusts and the remaining localities gypsum ones. The presence of sufficient Trebouxia spp. available to be acquired by the fungus could vary in each type of substrate, but to corroborate this hypothesis a complementary HTS approach from the substrate should be included.

Previous metabarcoding studies, Moya et al. [34] and Molins et al. [18], speculated that the diversity of microalgae in thalli of R. farinacea seemed to correlate with different parts of the thallus. Molins et al. [18] revealed no significant variation in microalgal diversity along the middle and apical parts, but, in contrast, a slightly greater microalgal diversity in the basal parts. In the case of B. zoharyi, no support for differences in microalgal communities was found among the three thallus section categories (i.e., central, middle and periphery). R. farinacea is a shrubby, fruticose lichen that grows epiphytically on all kind of trees [86], while B. zoharyi thrive always on soils as a crustose placodioid lichen dominant of the biocrusts. These differences in the thallus architecture and habitat type could influence the diversity of microalgae inside each thallus. No support for differences in microalgal communities was found among thallus sections; however, several thalli showed different predominant Trebouxia spp. at each section. This result corroborate that thallus parts selected for DNA extraction in metabarcoding analyses are key to not bias the total phycobiont diversity detected.

The diversity of phycobionts has only recently been explored by eDNA metabarcoding approaches, and has been focused on species within the Mediterranean basin to date [4,29,34]. In contrast to HTS approaches, traditional and largely applied DNA barcoding using Sanger sequencing was able to detect only the most abundant phycobiont in the thalli [17,87]. The HTS techniques have surpassed and almost substituted traditional Sanger sequencing because of their greater depth and resolution in the taxonomic assessment of multiple taxa at once. In this context, lichen symbioses have represented suitable microniches for the application and testing of the performance of HTS and metabarcoding to uncover the previously unknown species diversity [18,34,39,87,88,89,90]. The present study corroborates that inclusion of metabarcoding analyses is key to understand microalgal diversity in lichens. However, all analyses performed using the multicopy nrITS showed methodological limitations, which potentially bias the results presented, due to the variation in the copy numbers across the different microalgae species. Therefore, the relative abundance of algal groups could not accurately depict the true relative abundance of lichen-associated algae, given the potential for a very wide range of nrDNA copy numbers of these algal groups. Moreover, in the case of B. zoharyi, a crustose placodioid lichen with a characteristic lobed architecture and an intimate relationship with the soils, preparation of the thalli for metabarcoding analysis could be difficult as a percentage of the microalgae detected could originate from the substrate and not be intrathalline.

Paul et al. [87], compared Sanger sequencing and high-throughput metabarcoding for inferring photobiont diversity in lichens. They determined that if the second most abundant microalga exceeded 30% of the total HTS reads in a sample Sanger sequencing generally failed and generated ambiguous Sanger sequences showing double peaks. In this study, 25% of the samples analyzed by HTS and Sanger disagreed concerning the Trebouxia spp. detected, but in these samples the second phycobiont exceeded 30%. These controversial results reflect specimens showing algal multiplicity in which if two phycobionts are highly represented inside the thallus, Sanger sequencing could randomly amplify only one of them [35]. This study highlights that inclusion of HTS analysis is crucial to understand lichen symbiotic microalgal diversity.