Complex photobiont diversity in the marine lichen Lichina pygmaea

Abstract Lichens are a well-known symbiosis between a host mycobiont and eukaryote algal or cyanobacterial photobiont partner(s). Recent studies have indicated that terrestrial lichens can also contain other cryptic photobionts that increase the lichens’ ecological fitness in response to varying environmental conditions. Marine lichens live in distinct ecosystems compared with their terrestrial counterparts because of regular submersion in seawater and are much less studied. We performed bacteria 16S and eukaryote 18S rRNA gene metabarcoding surveys to assess total photobiont diversity within the marine lichen Lichina pygmaea (Lightf.) C. Agardh, which is widespread throughout the intertidal zone of Atlantic coastlines. We found that in addition to the established cyanobacterial photobiont Rivularia, L. pygmaea is also apparently host to a range of other marine and freshwater cyanobacteria, as well as marine eukaryote algae in the family Ulvophyceae (Chlorophyta). We propose that symbiosis with multiple freshwater and marine cyanobacteria and eukaryote photobionts may contribute to the ability of L. pygmaea to survive the harsh fluctuating environmental conditions of the intertidal zone.


Introduction
Lichen symbiosis is undergoing renewed scrutiny, in part because of the application of molecular ecology techniques such as metabarcoding (Hawksworth & Grube, 2020;Smith et al., 2020). Historic assumptions about the 'single fungus and single eukaryote alga or Cyanobacteria' nature of the symbiosis do not take into account tripartite lichens (i.e., those containing cyanobacterial and algal photobionts) (Rikkinen et al., 2002;Henskens et al., 2012), and have been further challenged by the more recent discovery of other fungi (Millanes et al., 2016;Spribille et al., 2016;Mark et al., 2020), eukaryote algae (chlorobionts) and Cyanobacteria (cyanobionts) (Henskens et al., 2012;Moya et al., 2017) within lichens in addition to the primary myco-and photobionts (Smith et al., 2020), alongside a complex diversity of heterotrophic bacteria (Grube et al., 2009). This increased diversity of organisms co-existing within the lichen thallus has added support to an ongoing reinterpretation of lichens as even more complex ecological associations comprising numerous interacting components (Honegger, 1991;Hawksworth & Grube, 2020).
Amongst this complex lichen community, photobionts are fundamental members responsible for providing sugars (and in the case of some cyanobacterial photobionts, available nitrogen) essential for the survival of the lichen complex (holobiont). Extended photobiont diversity is likely widespread in terrestrial lichens (Dal Forno et al., 2021) and has been hypothesized to be associated with increased ecological adaptive capability (Casano et al., 2011;del Campo et al., 2013;Muggia et al., 2013;Ruprecht et al., 2014). By recruiting a diverse range of photobionts, it has been proposed that lichens may benefit from enhanced ecological plasticity, facilitating the ability to survive a range of environments (Casano et al., 2011;del Campo et al., 2013;Muggia et al., 2013). For example, the soil crust lichen Psora decipiens supports multiple lineages of Trebouxia and Asterochloris allowing it to dominate soil crusts from deserts to alpine sites (Ruprecht et al., 2014). Additionally, photobiont 'switching' may allow for a further increase in fitness where multiple mycobiont lineages share locally adapted photobiont genotypes (Piercey-Normore & DePriest, 2001).
Changing environmental conditions are characteristic of the habitat occupied by marine lichens in the intertidal zone, with marine lichens being exposed to a uniquely dynamic range of pressures not experienced by terrestrial lichens including daily tidal cycles of seawater coverage. Lichina pygmaea (Lightf.) C. Agardh is a saxicolous fruticose lichen in the family Lichinaceae that inhabits the upper intertidal zone ( Figure 1A). It can be found throughout the North Atlantic coastline, including the area surrounding Plymouth, UK (Naylor, 1930). L. pygmaea was historically thought to form a symbiosis with Calothrix (Cyanobacteria) (Whitton, 2012), which was updated through molecular characterization to Rivularia (Ortiz-Álvarez et al., 2015). Marine lichens such as L. pygmaea remain relatively under-studied compared with terrestrial taxa, and the degree to which complex photobiont diversity occurs in marine lichens is yet to be determined. The potential for multiple cyanobionts within members of the Lichinaceae was proposed by Bubrick & Galun (1984), but the extent of this has not yet been investigated using molecular ecology techniques. Here, we perform a broad metabarcoding survey of the bacterial (16S rRNA gene) and eukaryote (18S rRNA gene) phototrophic taxa associated with the L. pygmaea thallus to establish the potential for complex photobiont diversity in a marine lichen.

Sampling
Sampling and processing were carried out according to the protocol described in Parrot et al. (2015). In summary, L. pygmaea thalli were collected from three locations close to Plymouth in the south-west of the UK ( Figure 1B) in spring 2016: Rame Head (lat = 50.3228, long = −4.2194), West Wembury (lat = 50.3159, long = 4.0856), and Church Reef (lat = 50.3159, long = −4.0844). At each site, a single thallus was sampled using sterile forceps, transported to the laboratory at 4°C and processed within 1 h. Each of the three thalli were inspected using a dissecting microscope to remove any debris or invertebrates and separated into 1 g subsamples (N = 5). To separate the intra-and extrathalline communities, thalli were washed twice by placing in a universal tube with 20 mL sterile seawater and rotating for 10 min. Wash water was filtered through a 0.45 μm cellulose nitrate filter. Both the washed thalli (intrathalline community) and filter (extrathalline community) were homogenized using a Mini-Beadbeater machine (Biospec Products) at 42 rpm for 40 s.

DNA extraction
DNA was extracted from the washed thalli (intrathalline) and the wash water filters (extrathalline) using the PowerSoil DNA Isolation Kit (MoBio, CA, USA) and stored at −20°C. The V4 region of the bacterial 16S rRNA gene was amplified using primers 515F and 926R (Caporaso et al., 2011) using Phusion® High-Fidelity DNA Polymerase (New England Biolabs) and PCR was performed under the following conditions: 94°C for 5 min, followed by 30 cycles of 94°C for 30 s, 53°C for 30 s and 72°C for 1 min, and a final elongation step at 5°C for 10 min. The V4 region of the eukaryotic 18S gene was targeted with the primer pair E572F and E1009R (Comeau et al., 2011), and PCR was performed under the following conditions: 98°C for 30 s, followed by 30 cycles at 98°C for 10 s, 55°C for 30 s, 72°C for 30 s, with a final extension at 72°C for 5 min. Sequencing was performed on the Illumina MiSeq platform (Integrated Microbiome Resource, Dalhousie University, Canada).

Bioinformatics
16S rRNA and 18S rRNA gene reads were processed through the DADA2 pipeline (Callahan et al., 2016), which resolves amplicon sequence variants (ASVs) at single-nucleotide resolution (Callahan et al., 2017). Reads were trimmed and truncated to remove primers and low-quality sequences respectively. Following error estimation, paired-end reads were merged, and chimeras were removed. Taxonomic classification of resolved ASVs was performed using the RDP naïve Bayesian classifier (Wang et al., 2007) using SILVA version 132 (Quast et al., 2013) for 16S rRNA gene ASVs and PR2 version 4.12.0 (Guillou et al., 2013) for 18S rRNA gene ASVs. Read counts, taxonomic classifications, and sample metadata were compiled as a phyloseq object prior to downstream analysis (McMurdie & Holmes, 2013). L. pygmaea reads were removed from analysis. A median of 15,115 16S rRNA gene reads and 1551 processed 18S rRNA gene reads per sample were recovered (Raw sequencing data are deposited in the European Read Archive under ENA study PRJEB43698). Alpha-(Observed, Shannon Index, Pielou's Evenness) and beta-diversity (Bray-Curtis dissimilarity) analyses were carried out using phyloseq and differential abundances of specific ASVs between the extra-and intrathalline samples were assessed using DESeq2 (Love et al., 2014), based on nonnormalized read count data and a significance threshold of α = 0.05, with a value of +1 assigned to the matrix to allow for statistical comparison between samples where no reads were returned. This analysis allowed for a direct test of which ASVs are enriched within the thalli and therefore most likely to be acting as photobionts. All analyses were performed in the R environment (R Core Team, 2017).

Phylogenetic trees
Identities of cyanobacterial (accession numbers MZ169867-MZ169909) and eukaryote algal (accession numbers MZ190002-MZ190020) ASVs identified by DESeq2 were confirmed by performing BLAST searches (megablast, excluding uncultured/environmental sample sequences, e-value = 0.05) against the NCBI nt database (9 November 2020) (see Supplementary material). Two 18S rRNA gene sequences (ASV_433 and ASV_251) had indeterminate identities and were removed from further analysis. To build phylogenetic trees, best BLAST hits for each ASV and an additional 2-5 close matches were obtained. Where possible, information about source material was recorded. L. pygmaea derived Rivularia sequences from Ortiz-Álvarez et al. (2015) were included in the cyanobacterial tree and photobionts of Hydropunctaria maura and Wahlenbergiella striatula from Thüs et al. (2011) in the algal tree. Sequences were aligned using MAFFT (using G-INS-i) and manually trimmed to remove gappy tails and poorly aligned regions. The 16S rRNA gene sequence ASV_170 was poorly aligned and excluded from further analysis. Phylogenies were reconstructed using the Cipres (Miller et al., 2011) implementation of RAXML-HPC (Stamatakis, 2014) using the GTR-CAT model and 1000 bootstrap replicates. Outgroups were specified for both Cyanobacteria (Pseudanabaena) and Ulvophyceae (Trentepohlia) trees. Resulting trees were first checked in FigTree (https://github.com/rambaut/figtree), before being annotated in R using ggtree (Yu et al., 2017). All figures were manually edited for aesthetics in Inkscape (https://inkscape.org/).

Results
Cyanobacteria dominated the bacterial community within all three thalli accounting for up to 49% of 16S rRNA gene reads (Figure 2A), while the green algae family Ulvophyceae (Chlorophyta) dominated the intrathalline eukaryotic community, accounting for up to 82% of the 18S rRNA gene reads (Figure 2A). Analyses of alpha-diversity indicated that despite increased abundance, Cyanobacteria exhibited significantly lower diversity (Shannon, P ≤ 0.001) and evenness (Pielou, P ≤ 0.001) inside the sampled thalli ( Figure 2B). Similarly, Ulvophyceae showed a non-significant (Shannon, P = 0.0781) reduction in diversity within the thalli and a significant reduction in evenness (Pielou, P ≤ 0.005) ( Figure 2B). Analyses of betadiversity revealed a clear distinction between the intrathalline and extrathalline communities of both Cyanobacteria (Adonis, P ≤ 0.001) and Ulvophyceae (Adonis, P ≤ 0.001) ( Figure 2C).
A total of six dominant (average normalized read count >20) cyanobacterial ASVs were significantly enriched (P ≤ 0.005) within at least one location ( Figure 3A). Three significantly enriched ASVs, ASV_6 (MZ169869), ASV_18 (MZ169873) and ASV_55 (MZ169883) accounting for up to 8.7, 4.6 and 1.5% of the intrathalline bacterial community respectively, and a fourth ASV_73 (MZ169887), which was only found at Rame Head where it comprised up to 6.1% of the community, all belonged to Rivularia, the established L. pygmaea photobiont (Ortiz-Álvarez et al., 2015). However, ASV_2 (MZ169867) and ASV_4 (MZ169868) (ASV_4 was only found at West Wembury) which made up to 20.6% and 18.1% of the intrathalline bacterial community respectively were both related to Pleurocapsa, a genus not previously associated with L. pygmaea. Other cyanobacterial ASVs were derived from the genera Acaryochloris, Phormidesmis and Synechocystis, some of which showed increased abundance in the extrathalline community. Within the phylogenetic tree (100 sequences, total alignment length = 1465 bp), all Rivularia ASVs occurred within a clade containing mostly marine and saline water derived sequences, while the highly abundant Pleurocapsa ASV_2 (MZ169867) fell within a freshwater dominated clade (Figure 4).

Discussion
Here we report a considerable diversity of cyanobacterial and algal lineages associated with the interior of the L. pygmaea thallus. Complex diversity of lichen cyanobionts such as this has not been widely reported, despite early indications from the Lichinaceae family (Bubrick & Galun, 1984). A previous molecular-based survey of L. pygmaea cyanobionts (Ortiz-Álvarez et al., 2015), which used conventional PCR and Sanger sequencing only amplified the dominant cyanobiont Rivularia. The whole community metabarcoding approach used here (which detects all sequence variants even at low abundance) suggests that the diversity of cyanobacteria contained within the thallus of L. pygmaea may be more complex than previously considered (Ortiz-Álvarez et al., 2015), and the existence of high diversity of photobionts within the lichen thallus such as shown here has likely been overlooked in many lichens. Our findings mirror those of Onuț-Bränström et al. (2018) who discovered multiple Trebouxia genotypes in Thamnolia and Cetratia thalli following Ion Torrent sequencing, compared with a single genotype recovered by Sanger sequencing. More widespread use of community metabarcoding approaches to investigate photobiont diversity may therefore reveal a greater diversity of photobionts and hitherto unrecognized lichen partnerships.
While many terrestrial tripartite lichens are well recognized (Rikkinen et al., 2002;Henskens et al., 2012), to the authors' knowledge the presence of Ulvophyceae within the L. pygmaea thallus is the first documented occurrence of a putative chlorobiont in this marine species. This potentially places L. pygmaea amongst several other cyanolichens shown to also contain chlorobionts (Henskens et al., 2012). In particular, L. pygmaea appears to associate with Blidingia, a blade-like multicellular alga that forms an unusual association with Turgidosculum ulvae (Verrucariaceae) whereby the mycobiont colonizes algal filaments (Pérez-Ortega et al., 2018). The exact nature of the relationship of L. pygmaea and Blidingia is yet to be established, but hints toward a unique interaction that has been previously overlooked. The other potential photobiont, Paulbroadya is known to form symbioses with several other freshwater and marine lichens in the family Verrucariaceae (Thüs et al., 2011;Gasulla et al., 2019) (Figure 5). The relationship here between the L. pygmaea-associated Paulbroadya and the same genus associated with the crustose marine lichen H. maura indicates the potential for sharing of photobionts between marine lichens. Photobiont 'switching' can occur between closely and distantly related lichens occupying the same habitat (Piercey-Normore & DePriest, 2001;Dal Forno et al., 2021), offering a fitness advantage through the acquisition of locally adapted photobionts, and is likely commonplace (Yahr et al., 2006(Yahr et al., , 2004Magain et al., 2017). The potential for photobiont sharing in L. pygmaea is supported by the repeated occurrence of the same Ulvophyceae ASVs between sample sites, and comparisons between photobiont communities in L. pygmaea and H. maura from the same site are now necessary to establish whether the same algal genotypes occur within both species of lichen. Interestingly, this is in contrast to several of the cyanobacterial lineages, including Rivularia, which exhibited variation between sites indicating either reduced dispersal and/or vertical transmission of cyanobionts leading to subsequent localized population structuring.  The phylogenetic placement of many of the potential L. pygmaea photobionts identified in this study within clades dominated by marine-derived lineages suggests that the lichen draws from a pool of marine photobionts, supporting the concept that photobiont selection is influenced by prevailing environmental conditions (Piercey-Normore & Deduke, 2011). Although the full ecological significance of the apparent complex photobiont diversity in L. pygmaea is yet to be determined, a selection of photobionts that enhance survival under different conditions could be a key factor allowing the lichen to occupy the marine environment. Other studies have proposed that the presence of multiple coexisting photobionts with different physiological properties may contribute to increased fitness of the lichen holobiont (Yahr et al., 2006;Casano et al., 2011). Since L. pygmaea is able to photosynthesize both in air and submerged in seawater (Raven et al., 1990), the ability of L. pygmaea to host a variety of photobionts may provide it with an ecological advantage with multiple marine-derived photobionts promoting survival in the intertidal zone by photosynthesizing when in seawater, while the presence of some freshwater-derived photobionts may allow for increased rates of photosynthesis at low tide and during periods of high freshwater input (e.g. rain events).
Further research into the relationship between L. pygmaea and associated photobionts is now necessary, including determining biogeographic factors in symbiont recruitment, establishing relatedness of photobionts between different marine lichen species and the extent of photobiont sharing, and photobiont localization within the thallus. Determining how much these photobionts contribute to carbon and, in the case of the cyanobacteria, nitrogen fixation will be important for understanding their contribution to biogeochemical processes, and relevant physiological experiments may be carried out using growth chambers and other experimental approaches (Henskens et al., 2012). RNA studies could help to quantify photobiont activity in relation to environmental parameters including seasonal, tidal and diel cycles. This work could contribute to our understanding of the