The Microbiomes of Seven Lichen Genera Reveal Host Specificity, a Reduced Core Community and Potential as Source of Antimicrobials.

The High Andean Paramo ecosystem is a unique neotropical mountain biome considered a diversity and evolutionary hotspot. Lichens, which are complex symbiotic structures that contain diverse commensal microbial communities, are prevalent in Paramos. There they play vital roles in soil formation and mineral fixation. In this study we analyzed the microbiomes of seven lichen genera in Colombian Paramos using 16S rRNA gene amplicon sequencing and provide the first description of the bacterial communities associated with Cora and Hypotrachyna lichens. Paramo lichen microbiomes varied in diversity indexes and number of OTUs, but were composed predominantly by the phyla Acidobacteria, Actinobacteria, Bacteroidetes, Cyanobacteria, Proteobacteria, and Verrucomicrobia. In the case of Cora and Cladonia, the microbiomes were distinguished based on the identity of the lichen host. While the majority of the lichen-associated microorganisms were not present in all lichens sampled, sixteen taxa shared among this diverse group of lichens suggest a core lichen microbiome that broadens our concept of these symbiotic structures. Additionally, we identified strains producing compounds active against clinically relevant microbial strains. These results indicate that lichen microbiomes from the Paramo ecosystem are diverse and host-specific but share a taxonomic core and can be a source of new bacterial taxa and antimicrobials.

Introduction

Symbiotic relationships between eukaryotes and microorganisms are ubiquitous (Gilbert et al., 2012), and often essential for the function and survival of the host, fulfilling roles that range from stress tolerance (de Zelicourt et al., 2013) and nutrient supply (Douglas, 2009; Ceja-Navarro et al., 2019) to defense against pathogens (Berendsen et al., 2012; Berg and Koskella, 2018). The composition of the microbial community associated with a particular host is defined by factors such as temperature and pH (Thomas et al., 2016), host genotype, nutrients (Shapira, 2016) and microbe-microbe interactions (Hacquard et al., 2015). Recent evidence indicates that hosts of the same species (Shapira, 2016), as well as evolutionarily related hosts (Ochman et al., 2010), harbor similar microbial communities. Core microbiomes, which are microbes consistently associated with a given host or found in a large fraction of samples from a particular environment (Shade and Handelsman, 2012), have been identified and cataloged in sponges (Thomas et al., 2016), corals (Ainsworth et al., 2015), insects (Sanders et al., 2014), plant roots (Yeoh et al., 2017) and mammals (Ley et al., 2008; Ochman et al., 2010).

Lichens represent some of the oldest and most diverse symbioses on Earth (Yuan et al., 2005). Lichens consist of a photobiont (cyanobacterium/alga) and a mycobiont (fungus), which together form a unique structure called the thallus (Nash, 2008). Lichens play a vital role in ecosystems as they are essential in soil formation, naked soil colonization, and nutrient uptake and release for plants (Nieboer et al., 1978; Seneviratne and Indrasena, 2006). Lichens can colonize a wide range of substrates, from natural surfaces to man-made materials such as plastic, rubber, metals and glass (Thomas and Nash, 2008). They can also tolerate extreme environmental conditions and offer a niche for diverse microorganisms (Shimizu, 2004; Garvie et al., 2008). The diversity of these lichen-associated microbial communities is not yet well characterized, and has only recently been investigated using high-throughput techniques (Cardinale et al., 2006; Grube et al., 2015). These studies indicate substantial microbial and functional diversity (Grube et al., 2015; Spribille et al., 2016; Cernava et al., 2017) that has been suggested to help protect the thalli against pathogens through the production of antimicrobials (Davies and Davies, 2010; Cernava et al., 2015a). The process of community establishment within lichens is poorly understood. It has been proposed that microbial composition is driven either by the photobiont (Hodkinson et al., 2012) or by geography and habitat (Aschenbrenner et al., 2014). In addition to differing from neighboring moss or bark habitats, it has also been suggested that microbial members present in adjacent mosses may facilitate lichen establishment (Aschenbrenner et al., 2017).

Recent studies have described the microbial communities associated with lichens using culture-independent strategies (Grube et al., 2015; Cernava et al., 2017; West et al., 2018). However, comparisons between studies are hindered by differences in sampling methods, data analyses, and the use of different lichen genera or species. Given the complexity of the lichen symbiotic structures, with recent evidence indicating that some lichens may be composed of multiple bacteria and more than one fungus (Spribille et al., 2016), it is important to study lichen microbiomes in order to understand their ecological role in the symbiosis.

A large and unexplored diversity of lichens is located in the Paramo ecosystem (Lücking et al., 2014). At least 264 species of lichens are described to date in Paramos and new species are still being reported (Celis et al., 2016; González et al., 2019). The Paramo ecosystem is considered a unique biodiversity hotspot that consists of high-elevation regions distinguished by extreme daily temperature variations, nutrient-poor and acidic soils, heavy rains and high UV radiation (Myers et al., 2000). This Andean ecosystem, as many in Colombia, has been understudied for decades due to armed conflict (Sabater et al., 2017). Here we characterize and compare the microbial communities in seven lichen genera from two Colombian Paramos. These lichen genera are widely distributed in these Paramos and are among the most taxonomically diverse (Vargas and Pedraza, 2004). Using 16S rRNA gene amplicon sequencing, we describe these microbiomes and identify members common to all seven genera, expanding our understanding of the complexity underlying these symbiotic structures. In addition, we isolated bacteria producing antibacterial and antifungal compounds and propose that Paramo lichen microbial communities should be further studied as a valuable source of antimicrobials. The diversity and prevalence of lichen-associated microbial communities underscore the need to further understand their ecological roles in lichen function.

Materials and Methods

Sampling and Sample Processing

Samples were taken at Los Nevados and Chingaza National Natural Parks in Colombia at altitudes ranging from 3,600 to 4,160 meters above the sea level, using sterile forceps and immediately placed in sterile plastic bags at environmental temperature. Several individual thalli were taken in order to have representative samples of different lichen genera in both localities. Samples were processed within 48 h of collection. Metadata such as GPS location and types of substrate where the lichen was collected (Corticolous: Wood, Terricolous: Soil and Saxicolous: Rock) were taken (Supplementary Table S1). Lichen morphological identification at genus level was carried out through herbarium specimen comparison.

Genomic DNA of lichens was extracted using the PowerSoil DNA Isolation Kit (Qiagen, Germantown, MD, United States) with some modifications: 20mg of lichen were homogenized in a FastPrep-24 (MP Biomedicals, Irvine, CA, United States) for two 20 s cycles at 4 m/s, and then processed according to the manufacturer instructions. The 16S rRNA gene V3-V4 region spanning a read length of 300bp, was amplified with primers V3F (5′-CCTACGGGAGGCAGCAG-3′) and V4R (5′- GGACTACHVGGGTWTCTAAT-3′) with barcoded Illumina adapters as describe in the standard procedures of the Earth Microbiome Project[1]. Blank controls were also included in amplification for quality assurance. Each 20 μL PCR reaction was prepared with 4 μL 5x HOT FIREPol master mix (Solis BioDyne, Tartu, Estonia), 2 μL of each primer (10 μM), 2 μL of sample DNA and 12 μL PCR-grade water. The amplicons were pooled in equimolar concentrations using SequalPrep plate normalization kit (Invitrogen, Carlsbad, CA, United States) and then purified with AMPure XP beads (Beckman Coulter, Atlanta, GA, United States). Amplicons were sequenced on the Illumina MiSeq platform (2 × 250 bp paired ends) at the Microbiology Department, Harvard Medical School in Boston, United States.

Analyses of Sequence Data

Illumina reads were quality checked with FastQC and edited with Trimmomatic (Bolger et al., 2014) to remove adapter and low-quality sequences that included reads with ambiguous nucleotides (Q value < 25) and short reads (<200 bp). Edited reads were processed in Mothur (v1.40) (Schloss et al., 2009), by first removing sequences longer than 430pb (screen.seqs: maxambig = 0, maxlength = 430). Files were reduced to non-identical sequences (unique.seqs and count.seqs) to minimize computational effort. Non-redundant sequences were aligned (align.seqs) to a trimmed SILVA (v132) bacteria database (pcr.seqs: start = 7 697, end = 23,444, keepdots = F) provided by Mothur (Quast et al., 2012). Only sequences that were aligned to the expected position were kept (screen.seqs start = 2, end = 15,747, maxhomop = 8; filter.seqs: vertical = T, trump = .). Aligned sequences were again reduced to non-redundant sequences and de-noised (unique.seq; pre.cluster), checked for chimeras using the VSEARCH algorithm (chimera.vsearch: dereplicate = t), which were then filtered out (remove.seqs). Sequences were classified (classify.seqs) based on the Greengenes database provided by Mothur (McDonald et al., 2012). Possible undesirable misclassified lineages were removed (remove.lineage taxon = Chloroplast-Mitochondria-unknown-Archaea-Eukarya). Sequences were then clustered (cluster.split: splitmethod = classify, taxlevel = 4, cutoff = 0.03) and converted to shared file format (make.shared: label = 0.03) assigning taxonomy to each OTU (classify.otu: label = 0.03, relabund = t). For alpha-diversity analysis reads were normalized to 20,623. Representative sequences of OTUs were retrieved based on the distance among the clustered sequences (get.oturep). The non-normalized shared file with OTU counts was used for differential abundance analysis in beta-diversity with ALDEx2 (Gloor, 2015).

Diversity Comparisons and Statistical Analyses

Diversity within samples (alpha-diversity) was analyzed with the Shannon-Weaver (Shannon, 1997) and Simpson Index (Simpson, 1949). Richness of microbial communities was assessed based on the observed number of OTUs and the rarefaction curves using the R package Phyloseq (McMurdie and Holmes, 2013). Multiple comparisons of richness and diversity measures were performed by one-way ANOVA, including Tukey’s (equal SD) or Tamhane T2 (non-equal SD) corrections. P values of <0.05 were considered to be statistically significant. Microbial community comparisons (beta-diversity) were first assessed with a similarity tree of samples based on the Bray-Curtis distance similarity matrix and the WPGMA hierarchical clustering method. We used ALDEx2 analysis (ANOVA-Like Differential Expression tool for compositional data) (Gloor et al., 2014) to find OTUs that define the differences between lichen microbiomes. The ALDEx2 R package decomposes sample-to-sample variation into four parts (within-condition variation, between-condition variation, sampling variation, and general unexplained error) using Monte-Carlo sampling from a Dirichlet distribution (aldex.clr: denom = “all”) (Urbaniak et al., 2014; Freitas et al., 2018). The statistical significance of each OTUs was determined by the general lineal model and Kruskal-Wallis Test (aldex.kw) for one-way ANOVA to determine OTUs significantly different for the seven lichen genera under study. The significantly differentially abundant OTUs were used to generate a Principal Coordinate Analysis (PCoA) based on the Bray-Curtis index and a prevalence matrix based on presence/absence. A Neighbor-Joining tree with differentially abundant OTUs and their abundances was built with OTU sequences aligned by an iterative refinement method (FFT-NS-i) (Katoh et al., 2002, 2017).

To display the taxonomy of OTUs present in each lichen microbiome, sequences were aligned in MAFFT v.7 with default settings (Katoh et al., 2002), and the cladogram for each microbiome was constructed using the average linkage method (UPGMA) (Sokal, 1958).

Core Microbiome

OTU prevalence (20,174 OTUs) was calculated based on the count mean of each OTU in every sample and cataloged as core (prevalence ≥0.9), pan (prevalence ≥0.25 and <0.9) or peripheral (<0.25). Core OTU sequences were aligned by an iterative refinement method (FFT-NS-i) and clustered by Neighbor-Joining (Jukes-Cantor Model) on MAFFT v.7 (Katoh et al., 2002). Core OTU relative abundances (CLR-transformed) in each lichen genus were displayed on a violin plot from Prism8 (GraphPad_Software, 2019). Core OTUs sequences were aligned to sequences in NCBI using Blastn optimized for highly similar sequences. Reference sequences were chosen based on >98% identity value. Both reference and core sequences were aligned and clustered with the same parameters mentioned above.

Bacterial Isolation and Screen for Antimicrobial Activity

Lichens were briefly washed with sterile water to remove sediment and loosely attached microorganisms (González et al., 2005; Parrot et al., 2015). Samples were aseptically divided into small pieces (∼0.5 cm) using sterile scalpels. The pieces were homogenized in phosphate saline solution with glass beads (4 mm-Marienfield) using a vortex. 100 μL of 10-fold serial dilutions in phosphate saline solution were plated on four culture media mainly for isolation of bacteria: Actinomycete Isolation Agar (AIA, Difco, BD), International Streptomyces Project medium-2 (ISP2), both supplemented with nalidixic acid (150 mg/L) and nystatin (50 mg/L) (González et al., 2005; Chevrette et al., 2019), Gause Synthetic medium, and Gause Oligotrophic supplemented with potassium dichromate (80 mg/L) (Wang et al., 2014).

Plates were incubated at room temperature (∼25°C) until no more new colonies appeared (up to 20 days). Colonies were isolated and purified based on morphological characteristics [color, surface (smooth or rough), shape (circular, filamentous, irregular or punctiform) and edge format (regular or irregular)]. DNA was extracted from purified colonies using the phenol chloroform extraction method with the following modifications: glass beads (0.4 mm) were used to lyse cells in a FastPrep-24 homogenizer (MP Biomedicals) with two cycles of 20 s at 4 m/s, adding 500 μL Tris–HCl buffer p.H 8.0, 200 μL NaCl 2.8M and 34 μL SDS 0.8% in 2 mL tubes. Bacterial strains were stored in the same medium in which they were isolated with 20% glycerol at −80°C.

Bacterial DNA was used to amplify the 16S rRNA using 27F and 1492R universal primers (Lane, 1991; Hayashi et al., 2002): 27F 5′-AGAGTTTGATCCTGGCTCAG-3′ and 1492R: 5′-ACGGTTACCTTGTTACGACTT-3′. Each 25 μL PCR reaction contained 12.5 μL CorpoGen PCR Master mix, 0.5 μL of each primer (25 μM), 9.5 μL PCR-grade water. PCR amplification was done by 3 min denaturation at 94°C; 35 cycles of 30 s at 94°C, 45 s at 55°C, and 60 s at 72°C; and 6 min elongation at 72°C. Isolates negative for 16S rRNA gene were corroborated as fungi by amplifying the ITS region with primers ITS5 5′-GGAAGTAAAAGTCGTAACAAGG-3′ and ITS4 5′-TCCTCCGCTTATTGATATGC-3′ (White et al., 1990).

Antimicrobial activity was assayed using the double agar layer assay (Hockett and Baltrus, 2017), against seven microbes (tester strains) of medical importance, some of which can present serious risk due to their antimicrobial resistance profiles (World Health Organization, 2014). Salmonella enterica subsp. enteritidis (CG100, poultry isolate), Escherichia coli ATCC 25922 (CMPUJ060), Klebsiella pneumoniae (CGKp16, clinical isolate), Pseudomonas aeruginosa PAO1 (CG1159), Acinetobacter baumannii (CG577), Staphylococcus aureus (ATCC 25923, CG378) and Candida albicans (clinical isolate, CG303). Isolated strains were grown on solid medium for 10 days, covered with Mueller-Hinton agar containing 100 μL of an overnight culture of each tester strain (Saffari et al., 2016), and incubated at 37°C for 24 h. Isolates that displayed a growth inhibition halo of the tester strains were considered as antimicrobial producers.

Results

Microbial Diversity Varies Depending on the Lichen Host

To study the structure of bacterial communities in lichens, we collected samples from different lichen genera at Paramo ecosystems within two national parks in Colombia. The 57 lichen samples were classified into eleven genera (Supplementary Table S1), but only seven genera (Cora, Hypotrachyna, Usnea, Cladonia, Peltigera, Stereocaulon and Sticta), which corresponded to 47 individual samples, were found in both locations. Samples with three or more biological replicates from the seven genera were used for microbial community analyses. DNA was isolated from individual lichen samples to identify microbial community profiles by 16S rRNA sequencing, which resulted in a total of 3,412,279 reads (mean per sample: 72,601).

A total of 20,174 operational taxonomic units (OTUs) were identified across all samples, which ranged from 100 to 1,955 OTUs per sample. Rarefaction curves indicated that this richness was adequately sampled as many samples reached saturation (Supplementary Figure S1). Simpson and Shannon diversity indices were calculated after randomly subsampling to the lowest number of reads (20,623), showing a broad distribution among samples. Diversity was significantly different between Usnea and Hypotrachyna lichens, Figure 1A and Supplementary Table S2, using ANOVA (p = 0.037 and p = 0.026 for Simpson and Shannon indices, respectively), and between lichens Usnea and Sticta, which had the smallest and largest number of OTUs, respectively (ANOVA test p = 0.010).

FIGURE 1

Microbial diversity varies across seven lichen genera. (A) Diversity and richness values for samples per lichen genus, measured by the Shannon and Simpson indices and the number of Observed OTUs, respectively. Multiple comparisons of richness and diversity measures were performed by one-way ANOVA with P values of <0.05 considered to be statistically significant. Asterisks indicate significant differences between Usnea and Hypotrachyna for Simpson (p = 0.037) and Shannon (p = 0.026) indices, and between Usnea and Sticta (p = 0.01) for Observed of OTUs. (B) Relative abundance of the five most abundant phyla and classes. Data are presented as mean ± standard deviation (SD). (C) Taxonomic diversity of lichens with the total number of OTUs for a given genus (shown in the center of cladogram).

Taxonomic assignment of OTUs showed that lichen microbiomes were predominantly composed by members of the phyla Acidobacteria, Actinobacteria, Bacteroidetes, Cyanobacteria, Proteobacteria, and Verrucomicrobia (Supplementary Figure S2). In general, Proteobacteria and Cyanobacteria were the most abundant phyla (Figure 1B). Most strikingly, there were clear differences in microbial composition for the communities associated with Hypotrachyna and Cladonia which had less Cyanobacteria than the communities associated with other lichens. At the class level, Gammaproteobacteria, Alphaproteobacteria, and Nostocophycideae were the most abundant taxa (Figure 1B). Again, there were taxonomic differences at the class level among lichen genera, such as the low abundance observed for Gammaproteobacteria in Hypotrachyna. Cladonia, Usnea and Hypotrachyna also had very low mean abundance of the Nostocophycideae within their microbiomes, which was more abundant in Cora. To obtain an overview of the similarities and differences in taxonomy of the lichen microbiomes sampled, we generated cladograms with the OTU sequences present in all samples from a given lichen genus (Figure 1C). This taxonomy depicts the predominance of taxa from the phylum Proteobacteria in all lichens and the similarity in taxonomic composition of these seven lichen genera at the phylum level, despite differences in the number of identified OTUs, ranging from the lowest number in Peltigera to the highest in Hypotrachyna.

Lichens Can Define Microbiomes and Share Core Members

We next analyzed if microbiomes differed based on lichen genus. A hierarchical clustering of lichen sample community composition based on a pairwise Bray-Curtis dissimilarity matrix, indicated that the microbiomes of samples belonging to the same lichen genus were more similar to one another that to those present in different lichen genera (Supplementary Figure S3). To determine the specific taxa driving these difference we used ALDEx2 (Fernandes et al., 2013) software to identify OTUs that were significantly different in abundance among lichen genera. In total we identified 177 OTUs with significant differences as determined by the expected p-value of the Kruskal-Wallace test and the general lineal model-ANOVA (p < 0.05). To compare the various lichen microbiomes, we constructed a prevalence matrix based on the presence/absence of OTUs using these 177 differentially abundant taxa. As can be seen in Figure 2A, the microbiomes from samples belonging to the same lichen genus were more similar to one another than to those from other lichens. A Principal Coordinates Analysis (PCoA) also showed that the microbial communities associated with Cora lichens, and to a lesser extent with Cladonia, clustered close together (Figure 2B), whereas no such clustering was observed when samples were distinguished by growth surface (rock, soil or tree bark) or geographical location (Chingaza vs. Nevados Paramo) (Supplementary Figure S4). A Neighbor-Joining tree constructed with these 177 OTU sequences again revealed a preponderance of taxa from the Proteobacteria, although the relative abundance of taxonomic families within the Proteobacteria varied according to the lichen genus (Supplementary Figure S5). For example, the Acetobacteraceae family was more abundant in Cladonia and Usnea, unlike the family Sphingomonadaceae that was more abundant in Cora and Sticta lichens. Other phyla such as Bacteroidetes, Verrucomicrobia, Cyanobacteria, Acidobacteria and Actinobacteria were also well represented and showed differences in abundance across lichens. Some phyla, like Armatimonadetes, Firmicutes and TM7, were represented by a single OTU. These analyses indicate that the microbial community is predominantly defined by the lichen host rather than by location or growth substrate.

FIGURE 2

Members of the same lichen genus share similar microbial communities. Data based on the 177 OTUs with significant differences in abundance found by ALDEx2. (A) Prevalence matrix based on OTU presence/absence. A WPGMA hierarchical clustering method was used to group the significantly different OTUs (shown on bottom) on a dendrogram based on a Bray-Curtis dissimilarity matrix. Each row represents the presence/absence of OTUs in a given lichen sample (shown left). (B) Principal Coordinates Analysis (PCoA) based on the Bray-Curtis index shows microbiomes of lichens Cora and Cladonia differentiated from other lichens.

To detect if lichens harbored a core microbiome, shared OTUs were identified by registering both the presence (prevalence) and the total counts (abundance) for each of the originally identified 20,174 OTUs in all samples. An OTU was considered to be part of the core microbiome if it was present in at least 90% of samples (≥90%) (Figure 3A). OTUs with a prevalence <25% were cataloged as peripheral, taxa that might be an extension of the environment or substrate on which the lichen grows. OTUs with a prevalence between ≥25 and <90% represent pan taxa that might be occasionally present in lichens but are not widely distributed across samples and lichen genera. Sixteen OTUs were shared among all lichens sampled (Figure 3A). Their abundances ranged from 2,777 to 29,245 counts per OTU. Core OTUs corresponded to Proteobacteria (eleven OTUs), Acidobacteria (four OTUs), and Cyanobacteria (one OTU) (Figure 3B). The eleven Proteobacteria OTUs belonged to three orders, Rhizobiales, Rhodospirillales and Sphingomonadales, while the Acidobacteria OTUs corresponded only to the order Acidobacteriales. The Cyanobacteria OTU remained unclassified according to the taxonomic assignment with the Greengenes database. While these sixteen core taxa represented a minor part of the total number of OTUs (Figure 3A) they were among the most abundant in the data set. However, their abundance was variable among the different lichen genera, as can be seen for the Cyanobacteria OTU (Supplementary Figure S6). Thus, these diverse lichens appear to harbor a small core microbiome.

FIGURE 3

Lichen core microbiome of samples from two Paramos. (A) Prevalence of the total 20 174 OTUs across the 47 lichen samples and their total number of observations (counts). Taxa were defined as Core (prevalence ≥0.9), Pan (prevalence ≥0.25 and <0.9) and Peripheral (prevalence <0.25). OTUs are colored by phylum. (B) Neighbor-Joining tree of the lichen core microbiome.

To see if these sixteen core taxa were also present in lichens from different geographical sites, we aligned our sequences against 16S rRNA gene sequences available from the NCBI database (Identity values >98%). Of 34 sequences identified, eleven corresponded to uncultured bacteria from lichens (Supplementary Figure S7). However only in one case the lichen species was identified (Ramalina pollinaria, GenBank-ID MG996731.1) (Supplementary Table S3). Most of the 34 sequences retrieved from NCBI corresponded to samples found in cold environments such as Glaciers and Tundra, which have environmental conditions comparable to those of the Paramo. These results suggest that the OTUs of our core microbiome are similar to bacteria found in lichens described in other studies.

Paramo Lichens Are a Rich Source of Antimicrobials

In order to investigate if the studied lichens harbored bacteria that produce antimicrobials, the original 57 lichens collected in both parks were processed by plating on multiple media. 122 isolates were obtained from 37 samples belonging to eleven lichen genera (Supplementary Table S4). Of these isolated strains, 112 were Bacteria and 10 were Fungi, based on PCR amplification of the 16S rRNA gene and ITS region, respectively. Approximately 62% of these isolates (n = 76) were obtained from lichens collected in Chingaza, while 38% (n = 46) were from Nevados. All isolates were tested for antimicrobial activity using a double agar layer assay against seven microbial strains (S. aureus, A. baumannii, P. aeruginosa, C. albicans, S. enterica, E. coli, and K. pneumoniae). 28% of the bacterial strains (n = 32) and 30% of fungi (n = 3) displayed antimicrobial activity against at least one of the strains tested, the majority of which (n = 26) were recovered from Chingaza lichens. We found antimicrobial activity against all seven strains tested, but the most detected activities were against A. baumannii, S. aureus, and C. albicans. Of the 35 isolates showing antimicrobial activity, 21 strains were active against multiple strains. Additionally, two bacteria isolated from lichens Psoroma and Yoshimuriella exhibited activity against a multi-drug-resistant K. pneumoniae strain cataloged in our laboratory as resistant to β-lactams (Cephalosporins, β-lactamase Inhibitors, Carbapenems, Monobactams), Amikacin, Ciprofloxacin and Meropenem based on Antimicrobial Susceptibility Standards values (Clinical and Laboratory Standards Institute [CLSI], 2019).

Discussion

Despite the importance of lichens for ecosystems, there is still limited understanding of their biology and, especially, of the assembly and function of the associated microbial communities. The pervasiveness of microbial communities associated with lichens suggest that at least some of these microbes may be more than transient associates in these symbiotic structures. Here we describe the microbiomes of lichens from the Andean Paramo ecosystems, high mountain habitats that harbor endemic species and are an important reservoir of lichen diversity. Colombia is known to harbor at least 10% of the described lichen species in the world (Bernal et al., 2015; Gonzalez-Salazar et al., 2017), likely a low-end estimate (Lücking et al., 2014). This study extends previous observations by characterizing and comparing microbiomes present in seven lichen genera, Cora, Hypotrachyna, Cladonia, Usnea, Sticta, Stereocaulon and Peltigera and, to our knowledge the first assessment of the microbiome composition of Cora and Hypotrachyna lichens. Here we identified both transient members, expected due to location and exposure to variable environmental conditions, and more permanently associated taxa that indicate a tight relationship and can orient studies aimed at understanding how these microbes contribute to lichen function and ecology.

The 16S rRNA gene analyses of these Paramo lichen microbiomes showed diverse and complex microbial taxonomic profiles that varied among samples, as has been observed in other studies (Aschenbrenner et al., 2016). The seven lichen genera sampled harbored microbiomes composed mainly of the phyla Proteobacteria, Acidobacteria, Verrucomicrobia, Bacteroidetes and Actinobacteria, consistent with previous reports (Grube and Berg, 2009; Bates et al., 2011). There was, however, variation in the relative abundance of these groups across the various lichens and among individual samples of the same lichen genus, an indication of the heterogeneity that can be expected in environmental microbiomes. The observed variation in community structures, particularly within a given genus, could be due to the fact that not all lichens were classified to the species level and therefore the samples from a given genus could include multiple species. In fact, some of the genera sampled here are considered to be among the largest in terms of the number of lichenized species, with Cladonia containing the most known species among the lichens studied (Lücking et al., 2017).

Despite the observed variability, lichen microbiomes showed a high abundance of Proteobacteria, as has been documented in other studies (Cardinale et al., 2008; Aschenbrenner et al., 2014). Members of this phylum are thought to play important roles in lichen symbioses by providing nutrients, mobilizing iron and phosphate, and fixing nitrogen (Grube et al., 2015; Sigurbjörnsdóttir and Vilhelmsson, 2016). Among the Proteobacteria, Gammaproteobacteria and Alphaproteobacteria were the most abundant groups. Unlike other studies of lichen microbiomes, the Gammaproteobacteria were more abundant than the Alphaproteobacteria in many of our samples. Alphaproteobacteria has been consistently reported as an abundant member of the bacterial microbiome in lichens (Cardinale et al., 2008; Grube and Berg, 2009; Bates et al., 2011; Mushegian et al., 2011). The detection of Gammaproteobacteria and their predominance over Alphaproteobacteria in our dataset, differs from previous studies and could be due to the fact that here we sampled different lichen genera from a novel geographic region (Grube et al., 2009; Mushegian et al., 2011). It is also possible that this discrepancy is due to methodological differences such as sample collection, time of storage before processing, which in our case was up to 48 h and could alter community composition, extraction protocols and sequencing and analysis platforms, all of which can have an effect on the resulting community profiles and subsequent inferences.

Taxonomic profiling revealed Cyanobacteria as abundant members of all microbiomes. In fact, Cyanobacteria were present in Cladonia, Hypotrachyna and Usnea, which are considered to be chlorolichens (Ahmadjian, 1993; Skaloud and Peksa, 2010; Anna et al., 2014; Rafat et al., 2015), that is to say, lichens with a green alga as its major or unique photobiont (Lange and Wagenitz, 2004). Tripartite lichens, which have both an algal and a cyanobacterial photobiont are known to make up a small number of lichens (about 3–4%) and introduce greater complexity to these structures since both photobionts can contribute to photosynthesis (Henskens et al., 2012; Rikkinen, 2017). While the presence of abundant Cyanobacteria could suggest an important role in these chlorolichens, such as nitrogen fixation and/or photosynthesis, further analyses would have to be done to confirm their precise functions. Given the limitation of Illumina sequencing, which provides information for only ∼300 bp of the V3-V4 portion of the 16S rRNA gene, these Cyanobacterial taxa could not be classified at higher phylogenetic levels with the available databases (Supplementary Figure S8). Additional metagenomic sequencing or the isolation of these microbial members would be needed to further determine if these are novel cyanobacterial species and to assess their possible roles within lichens.

Several studies have identified patterns in the structure of lichen microbial communities (Grube et al., 2009; Bates et al., 2011; Bjelland et al., 2011; Aschenbrenner et al., 2014). Specific bacterial taxa have been associated with some lichens (Hodkinson and Lutzoni, 2009), as well as predominance of particular groups, such as the Alphaproteobacteria in Cladonia arbuscula (Cardinale et al., 2008). These differences are thought to be driven by biotic and abiotic factors, of which the photobiont and large-scale geographical distance apparently determine the composition in different lichen types (Hodkinson et al., 2012). It has also been suggested that a lichen’s secondary metabolite production could drive microbiome structure (Hodkinson et al., 2012; Cernava et al., 2015b; Leiva et al., 2016). In this work, we used a differential abundance analysis approach (ALDEx2) to identify OTUs that varied in abundance among samples. ALDEx2 takes into account the compositional nature of microbiome data (Gloor et al., 2017), which means that the limited number of sequences obtained in any sequencing platform do not necessarily represent the number of sequences present in a given sample. This pipeline considers sample variation, which can be due to technical variations such as library preparation and sequencing output, to identify taxa that are significantly different between groups, reducing the false discovery rate frequently associated with other standard approaches for high-throughput sequencing data (Fernandes et al., 2013; Gloor and Reid, 2016; Gloor et al., 2017). This strategy also removes biases associated with standard data analysis that frequently defines bacterial community patterns based mostly on abundant taxa (Shade and Handelsman, 2012), and can overlook rare OTUs or low abundant taxa that may be important for host function (Reveillaud et al., 2014). This bias is evident in a variety of ecosystems where rare taxa have been seen to be essential for the dynamics of microbial communities (Dohrmann et al., 2013; Shade et al., 2014).

By using this differential abundance analysis, we identified a set of 177 OTUs with significant differences, from a total of 20,174 OTUs, that indicated that these microbiomes were not defined by geographical location or growth substrate (foliose, fruticose, or crustose), as reported in some cases (Aschenbrenner et al., 2014). In our work, lichen microbiomes appeared to be driven by the lichen genus, for the case of Cora and to a lesser extent for Cladonia samples. Interestingly, Cora is the only lichen sampled here that is known to have a Basidiomycete mycobiont (Lücking et al., 2017), suggesting that the fungal host could be important in shaping this microbiome. For Cladonia lichens, the community clustering was not as evident but might be further examined by identifying if lichen species play a role in defining community structure (Hodkinson et al., 2012), something that could not be done given that we did not classify the lichens to the species level (Grube et al., 2009). The 177 significantly different OTUs mainly belonged to the class Alphaproteobacteria (Supplementary Figure S5). However, the relative abundance of Alphaproteobacteria varied across the seven lichen genera. Interestingly, some taxa from the family Acetobacteraceae were more abundant in lichens Cladonia and Usnea, while the Sphingomonadaceae taxa were more prominent in lichens Cora and Sticta. Some authors have hypothesized that this variation of Alphaproteobacteria abundance might depend on the type of lichen photobiont (Hodkinson et al., 2012), with the order Rhodospirillales dominating in chlorolichens and Sphingomonadales in cyanolichens. However, further analyses to determine the type of photobiont in our samples are needed in order to explain the variation of Alphaproteobacteria taxa.

A core microbiome of 16 OTUs was identified across the 47 lichens sampled from two distant Paramos. The limited number of shared OTUs reflects the complexity and diversity of these lichen microbiomes and the fact that seven different genera, and possibly many uncharacterized species, were analyzed. The identification of shared microbial taxa was also surprising and suggests that additional members may be relevant to the function and success of these symbiotic structures, allowing them to thrive in extreme and nutrient-poor locations. Some of these core taxa were found to be similar to 16S rRNA gene sequences found in other lichens from distant sites including extreme environments such as the Arctic. Further sampling and deeper sequencing efforts might help to determine if this core is in fact conserved in other lichens. In contrast, the same analysis within a more tightly defined phylogenetic group, such as a single species, could identify a more robust core community. The lichen core microbiome included representatives of three phyla, Proteobacteria, Acidobacteria and Cyanobacteria. The Proteobacteria core members belonged to the Alphaproteobacteria class, which functional omics studies have highlighted as essential for nutrient supply and lichen growth (Cardinale et al., 2008; Hodkinson and Lutzoni, 2009; Bates et al., 2011; Erlacher et al., 2015). These Alphaproteobacteria were represented by the orders Rhizobiales, Rhodospirillales and Sphingomonadales that have been previously reported as crucial for the maintenance of lichens (Hodkinson and Lutzoni, 2009; Erlacher et al., 2015). Finally, our core microbiome indicated the presence of a single cyanobacterial OTU, even though there were other highly abundant Cyanobacteria in >60% of our lichen samples (Pan microbiome). Previous reports have shown that cyanobacterial symbionts can be shared among different lichen types (Rikkinen et al., 2002), while in other cases the lichen mycobiont might be strongly selective in the choice of cyanobiont (O’Brien et al., 2005; Rikkinen, 2013). Future studies focused on cyanobacterial specificity within lichen microbiomes could disentangle the roles that these taxa are playing within lichen thalli.

Soil bacteria have traditionally been the major source of antimicrobials (Crits-Christoph et al., 2018; Chevrette et al., 2019), but most of these compounds are derived from relatively few culturable microbial taxa (Cragg and Newman, 2013). With the rapid and widespread increase of multi-drug-resistant bacteria, there is a pressing need for new antimicrobials (Payne, 2008) that has prompted exploration of different ecosystems. Bacteria producing bioactive compounds have been isolated from some lichens such as Lobaria (Cernava et al., 2015a), and Cladonia (Davies et al., 2005) and even from marine lichens (Parrot et al., 2015). These bacteria belong mainly to the phylum Actinobacteria, a group well-known for its biosynthetic capacity and antimicrobial production (Barka et al., 2016; Liu et al., 2017). Members of this phylum have been consistently reported as members of lichen microbiomes (Selbmann et al., 2010; Bates et al., 2011; Bjelland et al., 2011), and in situ analyses of C. arbuscula have shown that these bacteria are located within the thallus structure (Cardinale et al., 2008). Our lichen microbiomes had a high abundance of actinobacterial taxa, and antimicrobial activity screening showed that lichen bacterial isolates produced molecules active against diverse microorganisms, including the multi-drug-resistant pathogen K. pneumoniae. Taken together, these results suggest that lichen microbiomes from underexplored ecosystems such as the Paramo, could be an important source of novel bacteria and antimicrobials. These antimicrobial-producing bacteria could be crucial for the defense of lichen thalli against pathogens or for the maintenance of microbial community balance within the symbiosis. In addition to further analyses of potential bioactive compounds, metagenomic studies of our lichen isolates should help to identify bacterial species, biosynthetic gene clusters and their metabolic potential.

Conclusion

Here we described the microbiomes of seven lichen genera (Usnea, Cladonia, Peltigera, Stereocaulon, Sticta, Cora, and Hypotrachyna), including the first description of the bacterial communities from Cora and Hypotrachyna lichens, and the presence of a core lichen microbiome. These Paramo lichen microbiomes were dominated by the phyla Proteobacteria, Cyanobacteria, Acidobacteria, Verrucomicrobia, Bacteroidetes, and Actinobacteria. These microbiomes varied among lichens and were distinguished based on host identity rather than location or growth substrate. Importantly, we found a core community of sixteen OTUs present in all samples. The core community was composed of members from only three phyla, Proteobacteria, Acidobacteria and Cyanobacteria, suggesting that there is high selectivity regarding which bacteria can establish close associations across all lichens. Microbes isolated from these lichens produced antifungal and antibacterial compounds which suggests that these ecosystems could be further probed as a source of natural products.

Safety Statement

Bacterial tester strains were manipulated following institutional biosafety rules for BSL-2 pathogens.

Author’s Note

The preprint version of this manuscript has been deposited at bioRxiv-the preprint server for biology as: Sierra et al., 2019. The microbiomes of seven lichen general reveal host specificity, a reduced core community and potential as source of antimicrobials. doi: 10.1101/789032.

Data Availability Statement

Sequence data of lichen microbiomes are available in NCBI under accession number PRJNA558995. OTUs and taxonomy tables together with the figure scripts are available on GitHub: https://github.com/mariaasierra/Lichen_Microbiome.

Author Contributions

MS carried out sampling, laboratory work, data analysis and manuscript writing. TS, DD, and CM helped with data and statistical analyses and manuscript editing. GP helped with sequencing. RK helped with sampling and analysis of data. BM carry out identification of lichens. MZ conceived the study, supervised work, helped with sampling and with writing of the manuscript. All authors provided input on the manuscript, read and agreed to the contents of the final version.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.