Host specialization and spatial divergence of bacteria associated with Peltigera lichens promote landscape gamma diversity

Abstract

Background

Lichens are micro-ecosystems relying on diverse microorganisms for nutrient cycling, environmental adaptation, and structural support. We investigated the spatial-scale dependency of factors shaping the ecological processes that govern lichen-associated bacteria. We hypothesize that lichens function as island-like habitats hosting divergent microbiomes and promoting landscape gamma-diversity. Three microenvironments —thalli, substrates, and neighboring soils— were sampled from four geographically overlapping species of Peltigera cyanolichens, spanning three bioclimatic zones in the Chilean Patagonia, to determine how bacterial diversity, assembly processes, ecological drivers, interaction patterns, and niche breadth vary among Peltigera microenvironments on a broad geographical scale.

Results

The hosts’ phylogeny, especially that of the cyanobiont, alongside climate as a secondary factor, impose a strong ecological filtering of bacterial communities within Peltigera thalli. This results in deterministically assembled, low diverse, and phylogenetically convergent yet structurally divergent bacterial communities. Host evolutionary and geographic distances accentuate the divergence in bacterial community composition of Peltigera thalli. Compared to soil and substrate, Peltigera thalli harbor specialized and locally adapted bacterial taxa, conforming sparse and weak ecological networks.

Conclusions

The findings suggest that Petigera thalli create fragmented habitats that foster landscape bacterial gamma-diversity. This underscores the importance of preserving lichens for maintaining a potential reservoir of specialized bacteria.

Background

Ecological processes governing the assembly of microbial communities have been investigated in both free-living and host-associated communities across a range of ecosystems, such as soils and aquatic environments, plant rhizospheres and leaf surfaces, and animal guts and surfaces [1]. The collective findings from these studies suggest that deterministic mechanisms tend to be predominant in free-living and plant-associated environments, generally associated with either abiotic conditions or biotic interactions [1, 2]. Contrastingly, stochastic processes take on a central role in ecosystems associated with animals, disconnected from selection based on differential symbionts’ fitness [1, 3]. Within host-associated communities, the host selection introduces additional complexity and feedback mechanisms into the assembly of its microbiome. Comprehending the ecological assembly processes of host-associated microbial communities could be particularly valuable due to their crucial role in providing essential functions for their hosts [4]. This understanding becomes pivotal in predicting the disturbance consequences of symbiotic complexes that fulfill critical services, thereby ensuring the stability and functionality of the ecosystems [5].

Underexplored and intriguing symbioses in this regard are lichens, which nowadays are considered self-sustaining miniature ecosystems resulting from the interaction between a resident fungus, known as the ‘mycobiont’, forming a symbiosis with extracellular microbial photosynthetic partners called ‘photobiont/s’, and an unspecified number of other microscopic organisms [6, 7]. These complex relationships weave throughout the evolutionary history of life on Earth, connecting diverse kingdoms and global ecosystems across space and time. Nonetheless, lichens maintain a unique status among all the cellular multi-kingdom symbioses. In most symbioses, a single organism acts as a structural scaffold, whereas lichens, which contain at least a filamentous fungus and a photobiont, have no a priori scaffold [8]. This association results in a unique joint structure, the lichen thallus, which represents various micro-niches for microorganisms. Microbial communities associated with lichens, renowned as classic models of successful and sustainable symbiosis, play a crucial role in the overall functionality of the meta-organism. They provide essential functions, including nutrient supply, resistance against biotic and abiotic stress factors, the production of vitamins and hormones, metabolite detoxification, and the degradation of aged lichen thallus components. As a result, lichen-associated bacteria engage in intricate interactions with the fungal, algal, or cyanobacterial partners, contributing to their hosts’ health, growth, and fitness [9].

Several studies have explored the bacterial species composition of terricolous, corticolous and saxicolous lichens within Lecanoromycetes class, encompassing genera such as Cladonia, Lecanora, Lobaria, Parmelia, Peltigera and Rhizoplaca [10,11,12,13,14]. Alphaproteobacteria emerge as the predominant and metabolically most active bacterial class within these lichens, alongside other notable bacterial lineages such as Firmicutes, Bacteroidetes, Verrucomicrobia, Acidobacteriaceae, Acetobacteraceae, Brucellaceae and Chloroflexi [15, 16]. Nevertheless, the specific bacterial community composition varied among different lichen species. Extrinsic factors, such as geographical distribution and environmental properties, along with intrinsic factors, like the phylogenetic relatedness of the mycobiont and photobiont, play pivotal roles in shaping the microbial communities within lichen [15, 17]. Besides, the acquisition of lichen microbiomes can occur through the local dispersal of vegetative propagules [15] and/or the recruitment of additional microbial strains from surrounding environmental resources such as the substrate [12]. All these variables ultimately culminate in establishing a distinct community structure in mature lichen thalli [15].

The prospect that the composition and functions of the bacterial microbiome may respond to ecological and climatic variations suggests potential enhancements in the adaptability of the holobiont. As lichen populations adapt to novel habitats, changes likely occur in their bacterial communities [16]. However, a remaining knowledge gap in the lichen microbiome study lies in comprehensively unraveling the critical assembly mechanisms of microbial communities under varying conditions. Bridging this knowledge gap is essential since community assembly is a dynamic process. Few studies have focused on the bacterial community assembly processes in biological soil crusts (BSCs), a life form found in the topsoil layer of drylands composed of lichens and other organisms such as cyanobacteria, algae, and bryophytes [18]. The bacterial community assembly in BSCs exhibits different drivers and patterns across successional stages [19] and across various lichen microenvironments [20]. Despite growing research interest in lichens, there remains a significant gap in understanding the spatial-scale dependency of assembly processes, ecological drivers, interaction patterns, and niche breadth of lichen-associated bacteria [21]. We hypothesize that lichens function as island-like habitats hosting divergent microbiomes and promoting landscape gamma-diversity due to (i) local stochastic colonization involving host ecological filtering of specialist taxa and (ii) regional dispersal limitations of these taxa influenced by geographic distance.

To address these issues, we selected the lichen genus Peltigera considering its wide distribution in Chilean Patagonia. This region provides an extensive latitudinal gradient encompassing diverse environmental contexts, bioclimatic features [22], and a rich abundance and diversity of Peltigera lichens [23]. Peltigera is a cosmopolitan genus of foliose and terricolous lichens forming symbiotic associations with cyanobacteria from the genus Nostoc (i.e., cyanobiont), although a few species associate with a green alga of the genus Coccomyxa as their main photobiont and Nostoc as their secondary photobiont [24]. More than ten bipartite species have been documented in Chile, including wide or narrow distributions [23,24,25]. To complement our understanding of how lichen microbial communities assemble, we focused on closely related species, which provide a more accurate mechanistic insight [26,27,28,29,30]. Specifically, we centered our study on four geographically overlapping Peltigera species spanning a 1,100 km transect in the Chilean Patagonian. Our objectives were to elucidate: (i) the balance between stochastic and deterministic ecological processes; (ii) the contribution of the hosts (mycobiont and cyanobiont), the environment (pedologic and bioclimatic properties), and the geography; (iii) the structure of potential co-interaction networks; and (iv) the habitat niche breadth and specialization degree, within bacterial communities of three micro-environments (i.e., soil, substrate, and thallus) associated with these lichens at local and regional geographical scales.

Methods

Study sites and samples

Specimens of terricolous Peltigera cyanolichens were obtained from the Peltigera lichen collection of the Microbial Ecology Laboratory of Universidad de Chile [23]. Specifically, we selected individuals of Peltigera frigida, Peltigera fuscopraetextata, Peltigera ponojensis/monticola 11, and Peltigera rufescens, due to their predominance in different habitats and sites from southern Chile. Eight specimens of each Peltigera species were collected from four sites in southern Chile: (i) Coyhaique National Reserve (Aysén Region, Chile; hereafter referred to as Coyhaique); (ii) Patagonia National Park (Aysén Region, Chile; hereafter referred to as Tamango); (iii) Karukinka Natural Park (Tierra del Fuego Island, Chile; hereafter referred to as Karukinka); (iv) Puerto Williams, (Navarino Island, Chile; hereafter referred to as Navarino). P. frigida and P. fuscopraetextata were collected from Nothofagus pumilio forests, while Peltigera ponojensis/monticola 11 and P. rufescens were collected from grasslands without forest cover (Fig. 1). At each site, eight samples of thalli and substrates of each Peltigera species and soils from each habitat (i.e., forest and grassland) were collected from approximately 100-meter-long transects. The whole sample set consisted of 126 thalli, 126 substrates, and 64 soils (Additional File 1: Table S1). Lichen thalli were collected at least 1 m apart to avoid resampling the same genetic individual. Substrate samples consisted of the soil beneath each thallus. Soil samples were collected at least 1 m away from any lichen to prevent its influence. Lichen thalli, along with their underlying substrates, were collected with clean tools and stored in new paper bags. Soil samples were collected at a depth of 10 cm, after removing the leaf litter with clean tools, and stored in sterile Falcon tubes. All samples were transported to the laboratory in cooled containers. In the laboratory, lichen thalli were separated from the attached substrates with a sterile brush and spatula. The thalli fragments were stored in sterile Eppendorf tubes, while the substrate and soil samples were sieved and stored in sterile Falcon tubes. All samples were stored at -20 °C.

Fig. 1

Localization of the sampling sites and examples of site topology and Peltigera spp. in Navarino Island, Chile (NAV). The pictures show the grassland (upper left) and the forest (bottom left) habitats and one of their two associated Peltigera species (P. fuscopraetextata, upper right, and P. rufescens, bottom right). White arrows point at the Peltigera thallus microenvironment. Sites are designated as follows: COY; Coyhaique, TAM; Tamango, KAR; Karukinka, NAV; Navarino

Genomic DNA extraction

Genomic DNA from thallus, substrate, and soil samples were extracted using the DNeasy PowerSoil Kit (Qiagen) with 0.15 g of thalli and 0.25 g of substrates and soils, with the following modifications: an overnight incubation at 4 °C after adding the precipitation solution C2 (step 8) and an increase in time from 30 s to 1 min in the centrifugation steps after bead beating (step 6), ethanol washing (step 16), and elution (step 21) [12, 31].

Molecular identification of lichen symbionts

The Peltigera genus presents several cryptic species and Nostoc cyanobionts that require multiple molecular markers to be identified [24, 25, 32]. Ribosomal barcode markers were then sequenced to confirm mycobiont and cyanobiont identification.

Mycobionts were identified by analyzing the fungal 28 S rRNA gene, amplified with primers LIC24R [33] and LR7 [34], and the ITS region, amplified with primers ITS1F [35] and ITS4 [36], following previously published methodologies [37]. The DNA sequences of 28 S rRNA and ITS markers were checked using SnapGene Viewer software (http://www.snapgene.com). The sequences were then concatenated in MEGA 11 [38] and aligned using MAFFT (v.7) [39]. Operational taxonomic units (OTUs) were defined based on 100% identity. Mycobiont identity was confirmed through de novo phylogenetic analysis of 28 S and ITS sequences, using the T-BAS web platform, with the Peltigera reference tree [40], by the implementations of the phylogenetic placement algorithm EPA (Evolutionary Placement Algorithm) in RAxML phylogeny program (Substitution model: GTRGAMMA, and Generic parameter: HKY85). The resulting phylogenetic tree was edited in iTOL [41] to group samples by color and to collapse certain sections and clades of the genus. Among our 126 concatenated sequences, 14 OTUs belonging to the section Peltigera were established, and their closest associations were determined, including P. frigida, P. fuscopraetextata, P. ponojensis/monticola 11, and P. rufescens (Additional File 2: Fig. S1). Next, the ITS1 spacer of the hypervariable region (ITS1-HR) of the samples was compared with previously published reference sequences [24, 42] to complement the mycobionts identifications. This analysis confirmed that OTUs match with P. frigida, P. fuscopraetextata, and P. rufescens species (Additional File 2: Fig. S2). Differentiation between P. ponojensis/monticola 11 and 1a was achieved by analyzing the complete ITS region. Four nucleotides (positions 47, 133, 221, and 535) identify individuals as P. ponojensis/monticola 11 (Additional File 2: Fig. S3).

Nostoc cyanobionts’ haplotypes were identified through a phylogenetic analysis conducted on the 16S rRNA gene, amplified with the primers PCR1 and PCR18 [43]. Amplification and sequence analyses were carried out in accordance with the protocols described previously [37]. The DNA sequences of 16S rRNA were checked using SnapGene Viewer software (http://www.snapgene.com), aligned using MAFFT (v.7) [39], and operational taxonomic units (OTUs) were defined based on 100% identity. The evolutionary history was inferred using the Neighbour-Joining method with Kimura’s 2-parameter distance using the MEGA11 software [38]. Inter-site rate variation was accommodated through a gamma distribution model (gamma parameter = 0.42). The dataset, comprising reference sequences obtained from cyanolichens and free-living cyanobacteria, resulted in a final alignment of 670 nucleotide positions. The resulting phylogenetic tree was edited in iTOL [41] to assign colors to the samples for better visualization. Among the 126 samples, five distinct Nostoc haplotypes were identified: C01, C02, C03, C07, and C14 (Additional File 2: Fig. S4), as previously designated [44].

Specimens of P. frigida and P. fuscopraetextata were associated with Nosctoc sp. C01 or C02, while P. ponojensis/monticola 11 and P. rufescens were associated with Nosctoc sp. C03 or C07. Only one individual of P. ponojensis/monticola 11 was found to be associated with Nostoc sp. C14 (Additional File 1: Table S1).

16S rRNA gene amplicon sequencing

The V4-V5 region of the microbial 16S rRNA gene was sequenced using the Illumina MiSeq platform (300 bp × 2) (Macrogen, South Korea) to characterize the microbial community in thalli, substrates, and soil. To minimize the amplification of the cyanobiont and enhance the sequencing depth of other bacteria, primers 515 F [45] and 904RM (5’-CCCCGTCAATTCCTTTGAGTTTYAR-3’) were used, with the corresponding sequences of Illumina adapters. The 904RM primer is a modified version of 904R from Hodkinson & Lutzoni [46] and was shown through in silico analysis to effectively reduce the amplification of representatives of Archaea and Eukarya domains, as well as the Cyanobacteria phylum. All amplifications were carried out in a 25 µl reaction volume containing GoTaq^® Green Master Mix (Promega, WI, USA), 0.2 mM of each primer, 0.5 µg/µl of BSA, and 10 ng of genomic DNA in a Maxygene thermocycler (Axygen, CA, USA) using a ramp rate of 0.1 °C/s. PCR conditions were: initial denaturation of 3 min at 94 °C, followed by 30 cycles of 30 s at 94 °C, 30 s at 60 °C and 30 s at 72 °C, with a final extension of 10 min at 72 °C. The quality and size of the amplicons were visualized in 2% (w/v) agarose gels in TAE 1x buffer (40 mM Tris-acetate, 1mM EDTA [pH 8.0]) stained with GelRed™ (Biotium, CA, USA). Additional PCR negative controls (template-free samples) were included and showed no amplification. The band of interest was cut from the agarose gel and purified using the Wizard^®SV Gel and PCR Clean-Up System kit (Promega, WI, USA) to decrease the presence of non-specific amplicons in sequencing. Sequences obtained from the microbial communities associated with lichen thalli, substrates, and soils were deposited in the Sequence Read Archive (SRA) database of the National Center for Biotechnology Information (NCBI) under the accession number PRJNA931999.

The forward and reverse reads acquired from the 316 samples were processed using Mothur (v.1.48.0) [47], following the procedure and criteria previously detailed for pairing, trimming, and alignment steps [48]. Sequences were clustered into operational taxonomic units (OTUs) at 97% identity, discarding the OTUs conformed by a single sequence. Exceedingly rare OTUs with relative abundance < 0.005% were removed [49], resulting in 1,578 OTUs that still accounted for more than 75% of the total sequences. The OTU table was rarefied at 14,592 sequences, leading to the removal of a single sample (soil sample from Navarino locality; NAV20_PS08_sue). Rarefaction curves of OTU richness were calculated for each sample with the get_rarecurve from the microbiotaprocess R package [50]. Rarefaction curves demonstrated that most of the bacterial diversity within the samples was adequately captured, particularly in the thallus microenvironment (Additional File 2: Fig. S5).

Alpha-diversity measures and beta-diversity partitioning

Shannon and Faith’s phylogenetic alpha-diversity indices were computed for soil, substrate, and thallus datasets, using plot_richness and estimate_pd from phyloseq and btools R packages, respectively [51, 52]. Statistical comparisons among microenvironments were performed across the entire dataset, sampling sites, and Peltigera species using a pairwise Wilcoxon test implemented under the geom_signif function from the ggpubr R package [53]. The beta-diversity was partitioned into turnover and nestedness components based on Jaccard dissimilarities of bacterial communities associated with Peltigera spp. microenvironments, using the beta.pair function from the R package betapart [54].

Ecological assembly processes

The neutral abundance-frequency distribution model [55] was applied independently to soil, substrate, and thallus datasets, both individually for each Peltigera species and by combining them, using an available R code [56] to evaluate the contribution of stochastic ecological processes in community assembly. The OTUs either above or below the 95% confidence interval of the neutral model were considered as over- or underrepresented, respectively. Further, the deviance from the neutral model predictions was calculated for the OTUs conforming to the 40 most abundant families in Peltigera spp. thalli, which were classified as part of the ‘core microbiome’ if present in 100% of the thallus samples, while the remaining OTUs were allocated to the ‘variable microbiome’. The Quantitative Process Estimates (QPE) method [57] was applied to assess the contribution of stochastic (i.e., homogenizing dispersal, ecological drift, and dispersal limitation) and deterministic (i.e., variable and homogeneous selection) community assembly mechanisms depending on Peltigera habitat and microenvironment. Briefly, homogenizing dispersal results from high dispersal rates in space and leads to more similar communities, while dispersal limitation increases community dissimilarity due to low dispersal rates that, in turn, drive ecological drift. The variable selection leads to greater divergence in microbial communities, both in terms of compositional and phylogenetic dissimilarities, driven by inconsistent selective pressures. In contrast, homogeneous selection promotes convergence of bacterial communities [57]. The QPE method relies on the Raup-Crick (RC) and the β-nearest taxon (βNTI) indices, calculated using null models with 999 bootstraps. The relative strength of each ecological process is inferred from RC and βNTI indices [58]. A previously published script deriving from the original implementation of the QPE method was then applied to the OTU table of each Peltigera habitat and microenvironment [59]. Permutation tests were employed to assess the statistical significance of changes in process proportions resulting from both the site and the host species. These tests were conducted using the perm.test available in the exactranktests R package [60].

Ecological drivers

The variation in microbial community compositions according to the microenvironment was partitioned using the varpart function from the vegan R package to assess the relative contribution of cyanobiont and mycobiont hosts’ phylogenies, pedology, climate, and geography. Phylogenies of cyanobiont and mycobiont hosts were reconstructed with the web service NGPhylogeny.fr, using PhyML with the Smart Model Selection (SMS) algorithm [61, 62]. Pairwise distances were extracted from both reconstructed phylogenetic trees using the cophenetic function from the ape R package [63]. As the varpart function is designed to accept a maximum of four explanatory matrices, the respective correlation of host distance matrices with the microbial community was tested with Multiple Regression on Matrix distance (MRM), implemented under the MRM function of the ecodist R package [64]. Subsequently, mycobiont and cyanobiont distance matrices were combined by weighting them with their respective MRM coefficients. Pedologic variables of substrate and soil samples, including organic carbon (OC), phosphorus (Olsen.P), nitrate (N.NO3), ammonium (N.NH4), water content (WC), and pH, were obtained from a public repository [65]. Climatic conditions for each site were retrieved from the WorldClim database [66] using the raster R package [67]. The altitude of sampling points was inferred from their GPS coordinates using Google Earth Pro (v.7.3.6.9796). Pedologic and climatic data were scaled and centered using the scale function in R and were analyzed through principal components analysis (PCA) (Additional File 2: Fig. S6 and S7). To prevent collinearity among explanatory factors, the scores from PC1 to PC3 axis, capturing 97.6% and 81.3% of the variation in the climatic and pedologic datasets, respectively, were transformed into Euclidean distances using the vegdist function from the vegan R package [68]. The geographic distances in kilometers were calculated from the longitude and latitude coordinates using the earth.dist function from the fossil R package [69], followed by Hellinger transformation with the decostand function from the vegan R package. Microbial OTU tables were converted in Bray-Curtis distances with the vegdist function from the vegan R package. The significance of hosts, pedology, climate, and geography fractions was tested with separate distance-based redundancy analysis (db-RDA), and their contributions in variance decomposition were represented under Venn diagrams.

Co-occurrence networks

Microbial co-occurrence networks were constructed for each soil, substrate, and thallus sample, considering their site, habitat, and/or Peltigera species of origin (i.e., 40 networks in total, including 8, 16, and 16 for soil, substrate, and thallus, respectively). To ensure networks’ comparability, OTUs that were detected in less than 50% of the samples within a specific condition and had a relative abundance lower than 0.1% were discarded from the respective rarefied dataset, ensuring that each condition comprised a minimum of 215 OTUs. These filters led to the removal of a single network (thallus of P. fuscopraetextata in Coyhaique). The construction process utilized the trans_network function, using the SParse InversE Covariance Estimation for Ecological Association Inference (SPIEC-EASI) approach [70], with a default number of 100 SparCC simulations [71]. Individual network properties (i.e., clustering coefficient, centralization, average path length, and natural connectivity) and shared nodes and edges among networks were determined using the cal_network_attr and trans_venn functions, respectively. The significance of network property differences was assessed using pairwise Wilcoxon tests. The relative abundances and order affiliations of the linked nodes of positive edges in the co-occurrence networks were determined with the edge_tax_comp function. All the functions used were implemented in the microeco R package [72].

Ecological niche breadth

The ecological niche breadth (B_N) was calculated individually for each community member (i.e., individual OTU). Subsequently, the mean community-level niche breadth was calculated for each soil, substrate, and thallus sample. The analysis relied on the Levins’ method, implemented in the microniche R package [73]. To categorize OTUs as generalists, specialists, or neutralists, we compared their Levins’ B_N values to a null distribution derived from 999 randomly generated OTU distributions. Specifically, OTUs with Levin’s B_N falling below the fifth quantile were classified as specialists, while those with values exceeding the 95th quantile of the null distribution were designated as generalists. OTUs with Levins’ B_N falling within the range of the fifth and 95th quantiles were categorized as neutralists. Dunn tests were employed to assess the statistical significance of the differences in proportions of neutralists, generalists, and specialists among Peltigera microenvironments using the dunn.test R package [74].

Results and discussion

Ecological filtering of non-neutral OTUs constrains diversity in thallus bacterial communities

Alpha-diversity levels (i.e., Shannon and Faith indices) within thalli were relatively homogeneous and consistently lower compared to soils and substrates, indicating a consistent ecological filtering across the four sampling sites and the four Peltigera host species (Additional File 2: Fig. S8, S9, and S10). This finding indicates that the colonization of Peltigera thalli is limited to a narrower range of bacterial taxa than soils and substrates, as previously observed in P. frigida and other lichen genera [12, 20]. In specific conditions, including Faith phylogenetic diversity in Coyhaique and Navarino, Shannon and Faith phylogenetic diversity in P. ponojensis/monticola 11 and P. rufescens, and Faith phylogenetic diversity in P. frigida, alpha-diversity within substrates was higher than in soil and thalli samples (Additional File 2: Fig. S9 and S10). Additionally, regardless of the Peltigera species considered, thallus shared more OTUs with the substrate than with the soil, soil shared more OTUs with the substrate than with the thallus, and substrate have low proportion of unique OTUs (Additional File 2: Fig. S11). These patterns suggest that in certain ecological contexts, substrates may act as transitional microenvironment, or ecotones, characterized by a higher level of biodiversity resulting from the coexistence of bacterial taxa from adjacent but contrasting microenvironments (i.e., soil and thallus) [75].

The frequency and abundance of the OTUs conforming soil, substrate, and thallus bacterial communities were relatively well described by the neutral model, which accounted for 54% (soil), 61% (substrate), and 72% (thallus) of the bacterial community variance (Fig. 2a, b and c). The strong fit observed for thallus bacterial communities suggests that deterministic processes play a limited impact on these communities. This result seems contrary to the classical expectation that hosts tightly control their resident microbial communities, reducing the relative importance of stochastic processes, as previously observed in other host systems [48, 76,77,78,79]. However, we suspect that this pattern might arise from thallus OTUs having globally lower disparities in abundance and frequency with the model predictions compared to soil OTUs, yet still falling outside the model confidence intervals, thereby inflating the R² value. The estimated migration rates (m values) were higher in soils and substrates (m = 0.15 and m = 0.19, respectively) compared to the thalli (m = 0.12), indicating fewer migrators that could easily disperse among thallus communities [77].

Fig. 2

Neutral model fitting to (a) soil, (b) substrate, (c) thallus microbiomes at the OTU level, and (d) relative abundances of neutrally distributed, overrepresented and underrepresented OTUs among the three microenvironments. Each point represents an OTU colored according to its deviance pattern from the neutral model predictions. The 95% confidence intervals (dotted lines) are drawn around the predicted occurrence frequency (solid line). OTUs within, above and below the 95% confidence limits of the neutral model are considered neutral, over- and underrepresented, respectively. Values of R² and m indicate the fit to the neutral model and the estimate of migration rate between communities, respectively

The proportion of bacterial OTUs deviating from neutral predictions substantially differed across microenvironments. While soils were dominated by neutral OTUs, comprising 43.1% of their communities, substrates and thalli were dominated by overrepresented OTUs (39.8 and 44.8%, respectively) (Fig. 2d). The dominance of overrepresented OTUs further confirms that Peltigera thalli are mainly composed of bacterial taxa that are either preadapted to (or selected by) the hosts, or more able to disperse and pass through host ecological filtering [56, 80, 81]. Performing the neutral model analysis separately for each Peltigera species revealed a marked decrease in the numbers of over- and underrepresented OTUs in favor of neutral ones, indicating the presence of species-specific bacterial OTUs and illustrating the host’s importance in shaping the thallus microbiome (Additional File 2: Fig. S12).

When examining thalli samples, the overrepresented OTUs were distributed across 135 families, with the highest prevalence in Xanthobacteraceae, Solirubrobacteraceae, WD2101, and Gemmataceae (Fig. 3, Additional File 1: Table S2). These families have been previously documented in association with different lichen species, further supporting their potential significance in lichen symbiosis. For instance, representatives of the aerobic chemoheterotroph family Xanthobacteriaceae have been previously reported in various lichen species, including the epiphytic lichen Lobaria pulmonaria growing on maple trees in the Alps and the saxicolous lichen Endocarpon in arid desert crusts, and were suggested as possible nitrogen fixers [20, 82]. In contrast to a previous study where Solirubrobacter was predominant in the substrate of a local population of P. frigida [12], we found that over two-thirds of the OTUs assigned to Solirubrobacteraceae were overrepresented in thalli. Solirubrobacteraceae was also previously detected in other lichens, Lepraria yunnaniana and Punctelia borreri from the Yunnan region in Southwest China [83], suggesting that this family may be specialist of lichen microenvironments and gathers potential species-specific strains. Lastly, WD2101 and Gemmataceae families have been suggested to be potentially involved in the degradation of lichen exopolysaccharides such as lichenin; however, their activity within the lichen ecosystem remains to be investigated [84, 85].

Fig. 3

Deviance from the neutral model predictions of the 40 most abundant families in Peltigera spp. thalli. These families represent 71.7% of the total sequences from the rarefied dataset. Colors indicate the deviation of the OTUs from the neutral model predictions. Shapes refer either to the variable or core status of the OTUs. The numbers and mean deviance values of under- and over-represented OTUs confirming the selected families are provided in Additional File 13: Table S2. “Acido.” and “Plancto.” refer to Acidobacteria and Planctomycetes phyla, respectively

Remarkably, more than 90% of the 32 OTUs identified as core microbiome within the bacterial community associated with Peltigera thallus were overrepresented (Fig. 3, Additional File 2: Fig. S13). The widespread detection of these core OTUs within the Peltigera genus across a large geographic scale and their potential host-driven selection in the thallus microenvironment further emphasized their potential beneficial functions for the host [86]. Two core OTUs of Peltigera thalli microbiome, including the most abundant one, were assigned to Psychroglaciecola. This genus has been originally isolated from an arctic glacial foreland [87] and detected in another Peltigera species (Peltigera ponojensis), however its ecological role in the lichen symbiosis remains unknown [88]. More recently, the Psychroglaciecola genus has been detected in atmospheric samples, indicating an airborne dispersal that could explain its ubiquity in Peltigera thalli over large geographical scale [89, 90]. Two other core OTUs were assigned to Burkholderiaceae. This family has been reported as particularly abundant in P. frigida thalli and suggested to contribute to the sulfur cycling in lichen [91]. These overrepresented core OTUs represent intriguing candidates for genomic exploration of their ecology in a symbiotic context.

The underrepresented OTUs were associated with only 65 families, primarily found in Beijerinckiaceae, Sphingomonadaceae, Acetobacteraceae, and Rhizobiaceae (Fig. 3, Additional File 1: Table S2). Moreover, these underrepresented OTUs displayed a narrower phylogenetic diversity than the overrepresented OTUs (mean pairwise phylogenetic distances 1.08 vs. 1.20, respectively, Wilcoxon test, p-value < 2.2e-16). Being underrepresented does not necessarily indicate that these taxa are being selected against by the host (e.g., pathogens); it may also indicate that these closely related OTUs could represent site- or host-specific taxa providing distinct crucial metabolic functions in certain ecological contexts [26, 56]. Consistently, these families were previously detected as dominant and consistent across the bacterial community associated with P. frigida, suggesting their ecological relevance for the host [12].

Local-scale variable selection drives divergence of thallus bacterial communities

When inferring the assembly processes governing microenvironment bacterial communities within the same habitat or species, dispersal limitation and ecological drift were detected as the two major assembly processes in soil and substrate communities (Fig. 4a, Additional File 1: Table S3). The habitat (i.e., grassland and forest) influenced the ecological processes governing bacterial community assembly in soils and substrates. Although, the contributions of deterministic processes (i.e., variable and homogeneous selection) were roughly the same between forest and grassland soils (0.9% and 14.3–19.9%, respectively), forest soils exhibited reduced dispersal limitation and higher levels of homogenizing dispersal and ecological drift in comparison to grassland soils (25.9% vs. 55.4%, 14.3% vs. 2.1%, and 44.6% vs. 29.8%, respectively) (Fig. 4a, Additional File 1: Table S3). These results markedly contrast with a previously published meta-analysis [1], which showed that in non-saline soils, processes are predominantly deterministic (i.e., variable selection), while stochastic processes, primarily represented by homogeneous dispersal, play a lesser role. However, the differences we reported between grassland and forest habitats suggest that the ecological processes cannot be generalized among soil samples. Previous studies consistently reported a diversity of factors, such as the plant aboveground diversity and density, and the belowground root biomass, among others, can modulate the ecological process shaping the soil bacterial community assembly [92, 93]. The high consistency in processes’ contributions in soils and substrates suggests that Peltigera hosts probably have a limited role in shaping the bacterial community of these microenvironments.

Fig. 4

Ecological processes affecting the assembly of bacterial communities associated with Peltigera microenvironments. The relative contributions of the different assembly processes in governing the bacterial communities of Peltigera microenvironments are provided (a) under individual sample distance (i.e., within the same microenvironment, habitat, site, and species) and averaged by species (for substrate and thallus) or habitat (for soil), and (b) under different combination of Peltigera species (habitat for soil) and sites (S. habitat: same habitat, D. habitat: different habitat, S. host: same host, D. host: different host, S. site: same site, D. site: different site). Details of the values associated with (a) and (b) are provided in Tables S3 and S4, respectively. The host- and site-induced changes in ecological process contributions (c) were evaluated by comparing mean contribution values of the processes from D. host S. site versus S. host S. site and S. host D. site versus S. host S. site, respectively. The significance of the changes was assessed using a permutation test, represented in the plots as follows: ***; p < 0.001, **; p < 0.01, *; p < 0.05. Error bars represent 95% confidence intervals. P. pon/mon 11 = P. ponojensis/monticola 11

In contrast to soils and substrates, variable selection dominates in thalli, contributing to 67.0–79.2% of the community assembly (Fig. 4a, Additional File 1: Table S3). A previous study on P. frigida reported that while substrates and thalli shared a large part of their bacterial communities, they remained distinguishable microenvironments [12]. Consistently, our results show that major shifts in ecological assembly processes drive differences in community composition between these two microenvironments. Surprisingly, the assembly processes driving convergent bacterial communities were negligible in thalli (homogenizing dispersal: 0.0%, homogeneous selection: 0–1.8%) (Fig. 4a, Additional File 1: Table S3). Since this part of the analysis exclusively considered interindividual variations without accounting for the host species and site factors, a consistent filtering of the surrounding bacterial community might have been expected among individuals of a given Peltigera species. Instead, our results suggest that individual heterogeneity among hosts and differential symbiont fitness may lead to bacterial community structure divergence among thalli [94]. This trend, commonly observed in plant surface samples [1], aligns with the notion that part of the lichen microbiome could be acquired from ‘bacterial rain’, as demonstrated in the phyllosphere [95]. However, the underlying mechanisms responsible for the ecological filtering in Peltigera thallus remain to be thoroughly investigated. One plausible scenario involves a ‘kill-the-winner’ dynamic within Peltigera thalli, where dominant community members are eliminated (e.g., by phages), followed by the stochastic colonization of the resulting niche space by specialist bacterial taxa [96]. While this mechanism has not been directly evidenced in lichens, a previous study reported a significant presence of viral genes associated with bacteriophages within Peltigera lichens [97]. Alternatively, early colonizers of the thallus could influence the stable state of the thallus bacterial community through the benefit of the ‘priority effect’ [98]. Both scenarios would result in a noisy community assembly across Peltigera populations, dominated by differentially selected and phylogenetically diverse dominant bacterial taxa [99]. This aligns with the observed low number and fluctuating relative abundances of the core bacterial taxa in Peltigera thallus and the variation in the most abundant OTUs among individuals within the same Peltigera species and site (Additional File 2: Fig. S13 and S14a). Additionally, the contribution of stochastic processes was marginally and positively correlated with the mean cumulative relative abundances of core bacterial OTUs in the thallus metacommunity (Additional File 2: Fig. S14b). This indicates that both fluctuations in core OTU abundances and colonization by new taxa contribute to the divergence of thallus bacterial communities. The high specificity of the Peltigera thallus microenvironment would cancel out the stochasticity of bacterial colonization, constraining proliferation to only specialist taxa that are less diverse yet phylogenetically overdispersed (i.e., less related than expected by chance).

Interhost and regional dispersal limitation depicts thalli as island-like habitats

When incorporating the site and Peltigera species (or habitat in the case of soils) factors, we observed that both factors had a relatively modest yet statistically significant impact on the contributions of assembly processes within microenvironment bacterial communities (Fig. 4b). Notably, there was an average increase in dispersal limitation in soils (+ 36.7%), substrates (+ 24.3%), and thalli (+ 6.8%), coupled with a decrease in ecological drift in soils (-29.0%), substrates (-20.5%), and thalli (-9.8%) across different sites or Peltigera species (or habitat) (Fig. 4c, Additional File 1: Table S4). The increase in dispersal limitation, and decrease in ecological drift, was particularly evident when both factors were distinct (Fig. 4b). While site and species (or habitat) switching significantly reduced a stochastic process leading to convergent bacterial communities in soils and substrates (i.e., homogenizing dispersal, -8.2% and − 2.1%, respectively), thalli remained primarily dominated by deterministic processes leading to bacterial communities’ divergence (Fig. 4c) [11, 13, 17]. This suggests that at least a portion of the Peltigera thallus-associated bacterial community may have its dispersal constrained, possibly due to its reliance on host phenotypic traits [100]. Both site and host species (or habitat) exert a similar influence on the relative importance of ecological processes that shape bacterial communities associated with Peltigera thallus. This amplifies the degree of interhost isolation and the beta-diversity divergence, contributing to higher gamma-diversity at regional scale. The higher bacterial communities’ turnover in thallus compared to substrate and soil samples further supports this pattern (Additional File 2: Fig. S15). Indeed, the bacterial taxa substitution, generally influenced by species dispersal capacity or the degree of specificity to different biotic or abiotic conditions [101], has been previously shown to enhance gamma-diversity of plant-associated bacterial communities [102]. Given these findings, we propose extending the ‘virtual island’ concept to Peltigera lichens [103], as previously observed with ectomycorrhizal fungi [104]. Thalli can be likened to fragmented habitats, akin to patches of host lichens in a non-host matrix. Within these ‘habitat patches’, the bacterial community structure diverges due to host-mediated ecological filtering and dispersal limitation, leading to nonequilibrium species composition, a phenomenon in line with predictions from island biogeography theory [105]. In other words, the variability among bacterial community thallus composition emphasizes the role of Peltigera lichens as enhancers of the landscape microbial gamma-diversity [106].

Major contribution of cyanobiont and climate on thallus bacterial communities

Using the variation partitioning analysis, we examined the relative contribution of each explanatory variable (hosts, pedology, climate, geography) on Peltigera-associated microenvironments. Pedology, climate, and geography fractions had the highest values in soils compared to substrates and thalli (R² = 0.16, p < 0.001; R² = 0.14, p < 0.001; R² = 0.10, p < 0.001, respectively, Table 1 and Additional File 2: Fig. S16). In contrast, pedology and geography fractions had the lowest contribution in thalli compared to soils and substrates (R² = 0.06, p < 0.001 and R² = 0.07, p < 0.001, respectively, Table 1 and Additional File 2: Fig. S16). The lower geography influence in thalli aligns with the predictions of dispersal limitation discussed earlier, indicating that thallus bacterial communities might benefit from higher dispersal capacity compared to soil and substrate communities. While the mechanisms behind the dispersal and recruitment of lichen microbiomes remain elusive, some authors previously suggested that lichen propagules may facilitate the co-dispersal of associated bacteria [107]. Furthermore, as demonstrated in the past [108], bird-mediated airborne dispersal of vegetative lichen propagules could contribute to the reduced geographical structure of thallus bacterial communities observed in our study. Notwithstanding these findings, a slight but significant geographic influence was still evident in shaping thallus bacterial communities. This suggests that thalli from the same site tend to have more similar communities than those from distant populations, thus indicating that propagule dispersion is somewhat limited to a small geographic distance [107]. To further elucidate the terrestrial connectivity of bacterial symbionts across Peltigera populations, quantification at a finer diversity level (‘microdiversity’) is required, as previously conducted in diverse ecological contexts [109,110,111].

The hosts’ factor emerges as the dominant driver of Peltigera bacterial community composition, significantly contributing to the variance observed in thalli (R² = 0.19, p < 0.001) and substrates (R² = 0.17, p < 0.001) (Table 1 and Additional File 2: Fig. S16). This robust correlation between hosts’ phylogenetic distance and dissimilarities within the bacterial community represents the first evidence of a phylosymbiotic signal at the intragenus level within cyanolichens, underscoring the influence of phylogenetically conserved historical traits across the Peltigera genus in shaping the ecological filtering of the surrounding bacterial communities [12, 17, 112]. This observation further underscores Peltigera species as “smart harvesters”, evolving to enrich, in specialized structures, host-adapted symbionts fulfilling multiple functions within the symbiotic system [15].

Interestingly, the MRM analysis revealed that the regression coefficient between bacterial community and cyanobionts exhibited a more than ten-fold magnitude increase compared to mycobionts in both microenvironments (Additional File 1: Table S5). The pairwise phylogenetic distances were significantly greater in the mycobionts compared to cyanobionts, but an inverse pattern was observed for diversity, as evidenced by the density pattern of pairwise comparisons (Additional File 2: Fig. S17). Overall, this aligns with the assumption that cyanobiont partners would be more locally adapted than their fungal counterparts and drive regional lichen guild evolution [113]. Consistently, a previous study reported a major influence of the photobiont in the bacterial community associated with various lichen species distributed across North and Central America [100].

Climate emerged as the second most important factor driving thallus bacterial communities (R² = 0.11, p < 0.001) (Table 1 and Additional File 2: Fig. S16). Previous studies have extensively documented the impact of climate variations on lichen host distribution and diversity, including within the Peltigera genus, responding to temperature and moisture conditions [114,115,116]. However, the detection of climate’s influence on lichen microbiome is unprecedented, likely facilitated by our broad latitudinal sampling encompassing temperate, temperate/sub-Antarctic, sub-Antarctic bioclimatic regions [22]. This finding strongly suggests that future variations in climate will likely affect the bacterial symbionts of lichens, potentially disrupting symbiotic homeostasis and functionality [117, 118], therefore altering their numerous ecosystemic services [119]. Furthermore, the contribution of climate was more than twice as high in thallus than in the substrate microenvironment (Table 1 and Additional File 2: Fig. S16). In this context, thallus bacterial communities, more exposed to climatic variations than those in the belowground parts (e.g., substrate), exhibit greater sensitivity to sunlight, air temperature, and humidity [120], making them potential ideal sentinels for tracking climate variations.

Table 1 Variance partitioning for beta-diversity of bacterial communities associated with Peltigera microenvironments, constrained by pedology, geography, climate, and host variables.

Sparse ecological network and enhanced specialization in thallus bacterial communities

Upon analyzing the driving factors of microbial communities, around 70–80% of the variations remain unexplained (Additional File 2: Fig. S16). Thus, we further explored the biotic factors represented by interactions between bacterial species through ecological networks. Most of the interactions were positive co-occurrences, without significant differences among microenvironments (soils: 54.0%, substrates: 52.6% and thalli: 52.7%, Kruskal-Wallis, chi-squared = 1.32, p-values > 0.5). However, positive co-occurrences in thallus bacterial communities involved different bacterial taxa compared to those from soil and substrate, the former being dominated by taxa belonging to the Acetobacteraceae family, followed by the Planctomycetes WD2101 and Sphingomonas genera (Additional File 2: Fig. S18). Thallus ecological networks harbored lower values of average degree numbers, clustering coefficient, centralization, density, and natural connectivity (Fig. 5a, c, d, e, and f), and higher values of average path length compared with soil and substrate (Fig. 5b). This aligns with previous findings [121], suggesting less complexity and interconnectivity in thallus microbial networks compared to spatially close microenvironments. Network complexity, influenced by resource availability [122], likely ties to higher diversity and abundance of bacteria associated with nitrogen and phosphorus cycling in soils and substrates than in thalli [31, 123].

Fig. 5

Comparison of co-occurrence networks’ properties between soil and Peltigera microenvironments. Each point corresponds to a network calculated with the samples from either the same site and species in the case of lichens or the site and habitat in the case of soils. Colors refer either to the habitat or the Peltigera species. The significance of the comparisons was assessed using the Wilcoxon test, represented in the plots as follows: ***; p < 0.001, **; p < 0.01, *; p < 0.05, NS; no significant. P. pon/mon 11 = P. ponojensis/monticola 11

Despite the substantial overlap of nodes, ranging from 40.7–59.2% across sites among Peltigera species and from 41.9–57.8% across Peltigera species among sites, and the presence of core OTUs (Additional File 2: Fig. S13), we did not observe shared edges across thallus ecological networks (Additional File 1: Table S6). These findings indicate that each lichen individual is defined by a unique arrangement of ecological interactions within its bacterial communities, which further align with the ‘virtual island’ concept and contribute to the acknowledged adaptable nature of lichen symbioses [124]. Consequently, this emphasizes the need for further research to elucidate the specific ecological mechanisms driving these intricate network structures. Overall, our results indicate that thallus ecological networks are less complex, less resilient, and more vulnerable to environmental changes, and thus more likely to collapse compared with soil and substrate communities [122].

To delve into the underlying causes for these interaction patterns, we assessed the habitat niche breadth of microenvironment communities and their distribution among neutralist, generalist, and specialist OTUs. Consistently with their relatively weaker network structures, thallus communities were associated with a narrower habitat niche breadth (p-values < 0.001, Fig. 6a), and exhibited a higher proportion of specialist OTUs (> 60%, p-values < 0.001, Fig. 6b) compared to those in soils and substrates. Remarkably, substrate communities displayed a broader habitat niche breadth than those in soils and thalli, hosting the highest number of generalist OTUs (Fig. 6a and b). Given their limited habitat range, specialists might be restricted to relatively unique habitats with specific favorable environmental conditions, thus constraining their prevalence across a broad geographic scale [125, 126]. The reduced niche breadth observed in thallus communities further emphasizes that the lichen communities are mainly influenced by species sorting, where local environmental conditions act as filters [127]. In addition, lichen sub-compartments, varying in age, as well as external and internal surfaces, create chemically and physiologically unique micro-niches, fostering the emergence of diverse bacterial communities [15]. Therefore, lichen thalli stand out as exceptional and intricate micro-ecosystems that might serve as important reservoirs of specialist bacteria, which could play an essential role in sustaining the ecological functions of their environments.

Fig. 6

Habitat niche breadth (a) and proportions of habitat neutralists, generalists, and specialists (b) of the bacterial communities associated with Peltigera microenvironments. (a) Means of the Levins niche breadth index (Bj) were calculated for each OTU within a community (i.e., sample). The significance of the differences between microenvironments was assessed using the Wilcoxon test (***; p < 0.001). (b) Different lowercase letters represent significant differences assessed by the Dunn test (all p-values < 0.001) in the proportion of habitat specialists, neutralists, and generalists among Peltigera microenvironments

Conclusions

Our study leverages an intriguing ecological context of four geographically overlapping Peltigera species spanning a large-scale transect to offer a comprehensive understanding of the ecological assembly processes, driving factors, and species coexistence patterns within bacterial communities associated with lichens and their surrounding environments on a broad geographical scale. We find that deterministic ecological filtering, primarily influenced by the cyanobiont and secondarily by climate, plays a central role in shaping the bacterial communities within Peltigera thalli. This selection, which varies across sites and species, contributes significantly to enriching landscape gamma-microbial diversity. The increase in hosts’ phylogenetic and geographical distances accentuates the importance of dispersal limitation in assembling thallus bacterial communities. Furthermore, the thallus bacterial community establishes interaction networks with low complexity and stability and exhibits a narrow ecological niche, thereby serving as an important reservoir of specialized bacteria. Overall, these findings underscore lichens as fragmented habitats and their vulnerability to environmental perturbations, thus emphasizing the critical importance of preserving these remarkable and fragile micro-ecosystems.