We studied 197 specimens of the R. decipiens group representing the 15 species-level lineages (Blázquez et al. in preparation). Seven of those species correspond to previously known species (Krog and Østhagen, 1980a,b; Krog, 1990). Another six represent new species which are going to be described in Blázquez et al. (in preparation). They are referred here with the following provisional names: R. delicata nom. Prov., R. fortunata nom. Prov., R. gomerana nom. Prov., R. papyracea nom. Prov., R. sabinosae nom. Prov., and R. sampaiana nom. Prov. Two further species-level lineages will remain unnamed until further material is available and are referred in the present study as Ramalina sp. 1 and Ramalina sp. 2. The specimens were collected throughout 66 localities spanning the seven Canary Islands, the islands of Madeira and Porto Santo of the Madeira archipelago, and the islands of Sal and São Vicente of the Cape Verde archipelago (Supplementary Table 1).

A small thallus fragment from each specimen (c. 10 mm²) was excised from the tip of a single laciniae with the help of razor blade and forceps under a Nikon SMZ800 stereomicroscope. First, we inspected the thalli under a x80 magnification to check for areas free of epiphytic microalgae and fungi. Thallus tips were the only zones free of epiphytes in all samples. The choice of thallus tips is also supported by the finding of Molins et al. (2021), who reported that there were no differences in photobiont composition among different parts of the thallus in specimens of a species of the same genus (R. farinacea) from the same biogeographic region as the present study. Samples were washed with acetone to remove secondary metabolites and stored at −80°C. After 1 h of freezing, they were pulverized using TissueLyser II (Qiagen) with two crystal beads. Genomic DNA was extracted using E.Z.N.A.® Forensic DNA Kit (Omega Bio-Tek), following the instructions of the manufacturer. The second region of the ribosomal internal transcriber spacer (ITS2) was used as barcode to prospect photobiont diversity. We used primers FDGITS2-f and FDGITS2-r (Dal Grande et al., 2018) for amplification. These primers included Fluidigm CS1 and CS2 universal oligomer sequences at their 5′ ends. PCR reactions were carried out in a total volume of 15 μl, containing 3 μl of template DNA, 0.3 μl of each primer (10 μM), 0.6 μl of MyFi DNA Polymerase (Bioline, Sydney, Australia), 3 μl of 5x MyFi Reaction Buffer, and 7.8 μl of distilled water. PCR settings consisted in an initial denaturation at 95°C for 1 min; 35 cycles of 95°C for 15 s, 54°C for 15 s, and 72°C for 15 s; with a final extension at 72°C for 5 min. PCR products were checked in 1% agarose gels stained with SYBR™ Safe DNA Gel Stain (Thermo Fisher Scientific). PCR products were quantified using the Qubit dsDNA HS (High Sensitivity) Assay Kit (Thermo Fisher Scientific) and pooled in equimolar concentrations for sequencing on a single MiSeq run (Illumina, United States) at the RTSF Genomics Core at Michigan State University (East Lansing, Michigan). The Genomics Core completed library preparation by secondary PCR targeting the Fluidigm CS1/CS2 oligos at the ends of the primary PCR products. Primers used for the second PCR included dual indexed, Illumina library compatible sequences. Finished libraries were bulk normalized using Invitrogen SequalPrep DNA Normalization Plates (Thermo Fisher Scientific) and the recovered products pooled. The pool was QC’d and quantified using a combination of Qubit dsDNA HS (Thermo Fisher Scientific), Agilent 4,200 TapeStation HS DNA1000 (Agilent Technologies), and Kapa Illumina Library Quantification qPCR assays (Kapa Biosystems, Inc.). This pool was loaded onto an Illumina MiSeq v2 Standard flow cell and sequencing performed in a 2x250bp paired-end format using a MiSeq v2 500 cycle reagent cartridge. Custom sequencing and index primers complementary to the Fluidigm CS1/CS2 oligos were added to appropriate wells of the reagent cartridge. Base calling was done by Illumina Real-Time Analysis (RTA) v1.18.54, and output of RTA was demultiplexed and converted to FastQ format with Illumina Bcl2fastq v2.19.1.

Raw sequence data were processed using the DADA2 pipeline (Callahan et al., 2016) in R 4.0.3 (R Core Team, 2020). DADA2 takes a set of demultiplexed paired-end fastq files, filters the sequences based on their length and quality, assembles them into error-corrected amplicon sequence variants (ASVs), and removes chimeric ASVs. We chose this pipeline because ASVs make possible to work at a finer resolution than other approaches based on the creation of operational taxonomic units (OTUs). Besides, ASVs allow for greater reproducibility, reusability, and comprehensiveness than OTUs (Callahan et al., 2017). Following Dal Grande et al. (2018) we considered that ASVs that represented <0.005% (i.e., n = 100) of the total reads were most probably exogenous to the symbiosis and did not include them in downstream analyses. We then assigned taxonomy to the ASVs using all sequences assembled by Muggia et al. (2020) as a taxonomic reference. This was done using the assignTaxonomy function in DADA2. Given that microalgae can show ecophysiological differences even when they are phylogenetically close (Sadowsky et al., 2016) we chose to conduct all analyses in this study at the ASV level. All sequences obtained in this study are available in the SRA (NCBI) under BioProject PRJNA764073.

ASVs were aligned using MAFFT v.7.308 (Katoh et al., 2002) as implemented in Geneious® v.9.1.8. Sequences of Asterochloris glomerata (MW043487, Pino-Bodas and Stenroos, 2021) and Vulcanochloris guanchorum (KR952330, Vančurová et al., 2015) were included as outgroup. We inferred phylogenetic relationships among ASVs using Maximum Likelihood using RAxMLHPC2 v. 8.2.4 (Stamatakis, 2014) as implemented in CIPRES Science Gateway (Miller et al., 2011). We selected the GTRGAMMA substitution model and carried out 1,000 rapid bootstrap pseudoreplicates to evaluate nodal support. We considered nodes with bootstrap values equal or higher than 70% to be significantly supported. The resulting phylogenetic tree was visualized in FigTree 1.4.4¹ and edited using Adobe Illustrator CS5 (Adobe Systems Inc., San José, United States). Additionally, we generated haplotype networks for the Trebouxia species represented by more than three ASVs. We used the TCS method (Templeton et al., 1992) as implemented in PopART v.1.7 (Leigh and Bryant, 2015) for network construction. These networks were edited in Adobe Illustrator CS5 to show the prevalence of the different recovered ASV in each sampled Ramalina species and island.

To explore whether the frequency of ASVs was different among islands and Ramalina species we adjusted a two-way ANOVA for each ASV. ASV frequency was used as the response variable and island and Ramalina species as explanatory variables. We used Levene’s test to check whether data for each ASV met the normality assumption necessary to carry out a two-way ANOVA. In the cases in which our data did not meet the assumptions, we employed the non-parametric Kruskal-Wallis test. This was done using the functions aov, leveneTest, and kruskal.test in base R 4.0.3 (R Core Team, 2020). ASVs occurring only in one specimen were not included in this analysis.

In order to explore the relative importance of mycobiont identity, macroclimate, and island in shaping photobiont communities we performed distance-based redundancy analysis (dbRDA, Legendre and Anderson, 1999). Following Shankar et al. (2017), we calculated an abundance-weighed phylogenetic UniFrac distance-based dissimilarity matrix. This was done with the R package phyloseq (McMurdie and Holmes, 2013). We transformed the ASV table into relative ASV abundances using the function decostand of the vegan package (Oksanen et al., 2012) with the method “total” and combined it with the taxonomy of the ASVs and the RAxML phylogeny previously generated into a phyloseq object using the merge_phyloseq function. The dissimilarity matrix was calculated from this object using the UniFrac function. The weighed UniFrac dissimilarity matrix was used as the response matrix in the dbRDA. We used three explanatory matrices: (1) first matrix corresponded to the Ramalina species, (2) second to the island where samples were collected, (3) and third matrix included five macroclimatic variables downloaded from WorldClim (Fick and Hijmans, 2017): average temperature (°C), precipitation (mm), solar radiation (kJ m⁻² day⁻¹), wind speed (m s⁻¹), and water vapor pressure (kPa). Because lichens are photosynthetic and poikilohydric symbioses, we chose to use these instead of the 19 widely used bioclimatic variables also available in WorldClim, that are derived from temperature and rainfall only. All layers of the five variables (data for each month of the year are available) were downloaded at a resolution of 30 s (~1 km²) and merged into year averages using the calc function of the raster R package (Hijmans et al., 2015). Then, variable values for each sample were obtained through their geographic coordinates using the extract function from the raster package. Before computing the dbRDA we checked for correlation between WorldClim variables in our samples. We calculated a Pearson correlation coefficient matrix and generated a correlation hierarchical cluster plot with the absolute correlation values. Precipitation and water vapor pressure were dropped as they were correlated (absolute correlation >0.5) with solar radiation and average temperature, respectively. Then we inspected the remaining variables to check their normality. To do so we generated a histogram for each variable with the hist function of base R. Solar radiation was found to be not normal and was log-transformed. dbRDA was then calculated with the dbrda function of the vegan R package (Oksanen et al., 2012). Variance partitioning between the three explanatory matrices was estimated using adjusted R². Adjusted R² significance was assessed for each fraction with a permutation-based ANOVA test with 2000 permutations. This was done in vegan using the functions varpart and anova.cca. Significance of the dbRDA as a whole was also assessed with a permutation-based ANOVA test with 2000 permutations, this time with the function anova. To visualize variance partitioning we generated a Venn diagram plotting the output of the varpart function.

Based on the ASVs abundance table we calculated for each Ramalina species: number of studied thalli, percentage of thalli that showed algal coexistence, minimum number of ASVs in one thallus, mean number of ASVs per thallus, and maximum number of ASVs in one thallus. For the species that showed algal coexistence, we also calculated: minimum relative abundance of the main photobiont found in one thallus, mean relative abundance of the main photobiont across all thalli showing algal coexistence and maximum relative abundance of the main photobiont found in one thallus.

We built two different adjacent matrices i × j, where i represents mycobionts and j ASVs. In the first matrix, interaction frequency between mycobiont i and ASV j was calculated as the sum of reads for ASV j found in all specimens of mycobiont i previously rarefied to the depth of the sample with the lower number of reads (Weiss et al., 2017). In the second matrix, number of read are not taken into account and the interaction frequency was based in the number of specimens of the mycobiont i where ASV j was found. The rarefaction was done with the rarefy_even_depth function of the phyloseq (McMurdie and Holmes, 2013) R package. The presence-absence conversion was done using the function decostand of the vegan package (Oksanen et al., 2012) with the method “pa.” For each network, we calculated the parameter d’ (Blüthgen et al., 2006) to estimate the specialization of the mycobionts toward the ASVs. This index ranges from 0 (no specialization) to 1 (high specialization). We also calculated species degree, which is the sum of links per species, and ND, which is species degree normalized by the number of possible partners. All analyses were carried out in the R package bipartite (Dormann et al., 2008), using the functions plotweb and specieslevel. Following the approach of Gómez et al. (2015), we used the quantitative modularity (Q) of the bipartite networks to explore the occurrence of different photobiont niches in the species of the R. decipiens group. Optimal modular configuration was calculated using the Beckett algorithm (Beckett, 2016) as implemented in the computeModules function of the bipartite package. The algorithm was run 20 times with 10⁸ MCMC steps and a tolerance level of 10^–10, retaining the iteration with the highest likelihood value as the optimal modular configuration. In order to check for significance, we built a null model based on 500 random networks calculated using the function nullmodel and the vaznull null model which implements the null model proposed by Vázquez et al. (2007) in which matrix connectance remains the same as in the original network but total marginals change. We calculated modularity for each of the random matrices and calculated a z-score as (X_observed − μ_null)/σ_null, being X_observed the actual value of the parameter, μ_null the mean of the parameter for the population of matrices in the null model, and σ_null their standard deviation. p-values for the z-score were calculated as the number of elements of the null model showing higher or lower values than the observed value divided by the total number of elements in the null model. To explore if the resulting modules were a consequence of the species geographic distribution we compared the geographic distances within and between modules. This was done with PERMANOVA tests using the adonis function of the vegan R package (Oksanen et al., 2012). The geographic distance matrix was obtained with the distm function of the geosphere R package (Hijmans et al., 2017).

We calculated genetic diversity statistics at both the island and Ramalina species levels. We used the software DnaSP v.6.12.01 to calculate two indices: haplotype diversity (Hd, Nei, 1987), and nucleotide diversity (π, Nei, 1987). π was calculated using the Jukes and Cantor correction (Lynch and Crease, 1990). Instead of using the raw ASVs abundance table as input, we generated DNA alignments for each species and each island, which included a representative sequence of each ASV present in each sample. The two indices were then calculated for each alignment.

We used the phylogeny of the R. decipiens group which will be published in a forthcoming paper (Blázquez et al. unpublished) to calculate the phylogenetic signal of the specialization, species degree, ND, haplotype diversity and nucleotide diversity. We employed two different methods to calculate the phylogenetic signal: Pagel’s λ (Pagel, 1999) as implemented in the phylosig function of the phytools R package (Revell, 2012) and Blomberg’s K (Blomberg et al., 2003) as implemented in the phylosignal function of the picante R package (Kembel et al., 2010). Pagel’s λ can have values between 0 (no phylogenetic signal) and 1 (Brownian motion model, that is, correlation among trait values is proportional to the extent of shared ancestry). Blomberg’s K can be interpreted as follows: if K = 1 the trait under analysis follows a mode of evolution that is consistent with Brownian motion. If K > 1 close relatives are more similar than expected under Brownian motion. On the other hand, in K values <1 close relatives are less similar than expected (Blomberg et al., 2003).

We obtained a total of 2,949,259 raw reads, of which 2,175,408 passed DADA2 quality filter. After removing ASVs represented by less than 100 reads, 66 ASVs remained (Figure 1). The 66 ASVs clustered into the four main Trebouxia clades (Leavitt et al., 2015; Muggia et al., 2020), but not into the newly reported clade D (Xu et al., 2020). 62 ASVs belonged to four formally described species (i.e., Trebouxia jamesii, T. aggregata, T. decolorans, and T. australis) and to nine of the undescribed Trebouxia lineages reported by Muggia et al. (2020). Four ASVs did not match with any of the sequences assembled by Muggia et al. (2020). These ASVs represented three new lineages, two of them in clade A and one in clade I. Thirty-three of the 66 ASVs clustered into Trebouxia A39. This clade was identified as Trebouxia sp. TR9 through its reference sequence (FJ418565). The most common photobionts found in this study were ASVs 1 and 2, both belonging to Trebouxia sp. TR9 (Figure 2) and diverging by just one nucleotide. They were present on all islands and often coexist in the same localities. These two ASVs were the most common photobionts in the Canarian and Madeira archipelagos, together they represented the main photobionts in 78 and 85% of the Canarian and Madeiran thalli, respectively. Regarding the Cape Verde archipelago, the main photobionts from all samples corresponded to one of the five ASVs placed in the Trebouxia clade C25. On the other hand, 41 ASVs were only found in one thallus, although only 6 of them acted as the main photobiont.

ASVs belonging to Trebouxia clades A and C were found in all archipelagos. Clades I and S were not found in Cape Verde samples. Mycobiont identity had a limited effect on the frequency of ASVs (Table 1). Of the 25 ASVs present in more than one thallus only one, belonging to the species Trebouxia jamesii, was significantly associated to some Ramalina species. The island of origin had a greater influence on ASV frequency, with a significant effect in 8 of the 25 included ASVs. The greatest effect was observed in ASVs 5, 7, and 11. These three ASVs belong to Trebouxia clade C 25 and act as main photobionts in thalli collected in Cape Verde, but also occur in small numbers in El Hierro, La Gomera, and Porto Santo. An interaction between island an mycobiont identity was found in ASV 9. However, this ASV only occurs in three thalli: Two specimens of R. decipiens collected in Fuerteventura and one of R. sp. 2 collected in La Gomera.

The dbRDA analysis returned a F = 3.265 with a p = 0.001. The constrained axes explained 0.34 of the variance in ASV community structure (Figure 3). Island was the variable with the greatest simple effect. On the other hand, macroclimate was the variable with the smallest simple effect. However, the main driver of ASV community structure was the mixed effect of island and macroclimate, explaining a greater variance proportion than all other simple and mixed effects combined. Mycobiont identity was the least important variable, taking its simple and mixed effects into account.

Out of 197 thalli studied, 102 showed coexistence of several ASVs (51.77%). Specimens from the Canary Islands and Cape Verde showed a higher proportion of multiple ASVs within the thallus (56 and 58%, respectively) than specimens from the Madeira archipelago (38%). Algal coexistence was observed in all species (Table 2) except Ramalina sp. 1, which is most likely an artifact given that the species is represented in our dataset by only one specimen. The mean number of ASVs per thallus was 2.19, with a median of 2. Mean ASV count per thallus ranged from 1.30 in Ramalina hamulosa to 4.67 in Ramalina pluviariae. The sample with the highest ASV count was a Ramalina maderensis collected in Lanzarote (SPO 6337) containing 10 algae. However, on average the main photobiont amounted to 94% of the reads in thalli showing algal coexistence. The mean abundance of the main photobiont in thalli showing coexistence ranged from 83% in Ramalina erosa to 99% in Ramalina sp. 2. There were only six thalli, of varied origin and representing four Ramalina species, in which the main photobiont accounted for less than 70% of the reads. In most of them, the main and secondary photobionts were different ASVs belonging to the same Trebouxia species.

The bipartite interaction networks (Figure 4; Supplementary Figure 1) clearly depict ASVs 1 and 2 as the main photobionts of the Ramalina decipiens group. Ramalina sp. 1 is the only species not associated with these ASVs. This is most probably due to the poor representation of this species in the dataset. Excluding Ramalina sp. 1, whose only thallus associates with another ASV belonging to Trebouxia sp. TR9, specialization (d’) based on rarefied reads varied from 0.98 in Ramalina delicata to 0.51 in Ramalina maderensis. Species degree varied from 37 in Ramalina decipiens to 4 in Ramalina erosa and ND varied from 0.56 in Ramalina decipiens to 0.06 in Ramalina erosa. Species degree and ND were highest in the two most sampled species (R. decipiens and R. maderensis). Specialization parameters obtained from the data non considering read counts returned similar results for species degree and ND, but overall lower values of d’, that ranged from 0.59 in R. erosa to 0.27 in R. decipiens. All specialization parameters are summarized in Table 3. The network based on rarefied reads was significantly modular (Q = 0.18, z-score = 24.75, value of p <0.001), and detected the presence of six modules (Figure 5). The modularity analyses of the second network resulted in a non-modular pattern (Q = 0.21, z-score = −0.68, value of p = 0.241). The modular structure recovered in the first analysis seems to be a consequence of the mycobiont species distributions inasmuch as the geographic distances between and within modules were significantly different (F = 7.98, R² = 0.17, p < 0.001).

Nucleotide diversity (π) and haplotype diversity (Hd) were calculated at the island and species levels (Figure 6). At the island level, π ranged from 0.036 in Tenerife to 0.096 in Sal. Hd values ranged from 0.678 in Porto Santo to 0.928 in Fuerteventura. On average, Cape Verde had the highest π value (0.092) followed by Madeira (0.060) and the Canary Islands (0.052). Regarding Hd, Cape Verde again showed the highest value (0.810), followed by the Canary Islands (0.803) and Madeira (0.730). At the species level, π values ranged from 0.014 in Ramalina hamulosa to 0.077 in Ramalina papyracea. Hd values ranged from 0.709 in Ramalina erosa to 0.905 in Ramalina pluviariae.

The parameters d’, species degree, and ND derived from both networks were not phylogenetically structured according to Pagel’s λ and Blomberg’s K. Regarding genetic diversity, neither π nor Hd were phylogenetically structured either. Phylogenetic signal statistics are summarized in Table 4.