Lichen holobionts show compositional structure along elevation

Holobionts are dynamic ecosystems that may respond to abiotic drivers with compositional changes. Uncovering elevational diversity patterns within these microecosystems can further our understanding of community‐environment interactions. Here, we assess how the major components of lichen holobionts—fungal hosts, green algal symbionts, and the bacterial community—collectively respond to an elevational gradient. We analyse populations of two lichen symbioses, Umbilicaria pustulata and U. hispanica, along an elevational gradient spanning 2100 altitudinal metres and covering three major biomes. Our study shows (i) discontinuous genomic variation in fungal hosts with one abrupt genomic differentiation within each of the two host species, (ii) altitudinally structured bacterial communities with pronounced turnover within and between hosts, and (iii) altitude‐specific presence of algal symbionts. Alpha diversity of bacterial communities decreased with increasing elevation. A marked turnover in holobiont diversity occurred across two altitudinal belts: at 11°C–13°C average annual temperature (here: 800–1200 m a.s.l.), and at 7°C–9°C average annual temperature (here: 1500–1800 m a.s.l.). The two observed zones mark a clustering of distribution limits and community shifts. The three ensuing altitudinal classes, that is, the most frequent combinations of species in holobionts, approximately correspond to the Mediterranean, cool‐temperate, and alpine climate zones. We conclude that multitrophic microecosystems, such as lichen holobionts, respond with concerted compositional changes to climatic factors that also structure communities of macroorganisms, for example, vascular plants.


| INTRODUC TI ON
The holobiont concept holds that a multicellular host organism and its associated microorganisms form a discrete, yet dynamic biological unit (Simon et al., 2019).Members of the holobiont collectively ensure functioning of the entire consortium.An intact holobiont can therefore be considered an ecosystem that is in a state of dynamic equilibrium (Pita et al., 2018).Like other ecosystems, holobionts will thus respond compositionally to environmental change, for example along climatic and chemical gradients (Calegario et al., 2021;van der Loos et al., 2019;Peay et al., 2017).Moreover, biotic factors within the holobiont, such as host ecology or genotype, will govern its internal community structure (Bálint et al., 2013;Dal Dal Grande, Rolshausen, et al., 2018;Glasl et al., 2019).Understanding the environmental conditions and internal biotic interactions, that trigger compositional change in holobionts can help us evaluate vulnerabilities and adaptive capacities of such microecosystems, and possibly predict approaching tipping points under ongoing changing conditions (Pita et al., 2018).It can also deepen our knowledge of the roles which microorganisms play in their host's biology and ecology.
To address these questions, we here assess the community composition of lichen holobionts-integrated symbiotic units that consist of a fungal host, a suite of microalgal symbionts, and a bacterial community-along a climatic gradient.
The internal diversity and community composition of holobionts is shaped by a complex interplay between stochastic aggregation dynamics on the one hand, and selective factors mediated by the host and/or the environment on the other hand (Bordenstein & Theis, 2015;McFall-Ngai et al., 2013;Sullam et al., 2012;Theis et al., 2016).Correspondingly, functional and ecological integrity of holobionts throughout their distributional range will in part be contingent on internal compositional changes under distinct environmental conditions or along environmental gradients that show pronounced climatic transitions (Brener-Raffalli et al., 2018;Pantos et al., 2015;Theis et al., 2016;Zilber-Rosenberg & Rosenberg, 2008).
It remains contentious, however, if and under which circumstances selection on holobiont integrity might indeed supersede individual selection on hosts and microorganisms (Douglas & Werren, 2016;Kopac & Klassen, 2016;Moran & Sloan, 2015).To conclusively disentangle different selection regimes that potentially govern holobiont community composition, an important first step is thus to distinguish particular fractions of the overall diversity that do (or do not) co-vary with environmental and/or host-genetic variation along environmental gradients (Brener-Raffalli et al., 2018).
The degree of holobiont coherence can be evaluated based on the extent to which hosts and their associated microorganisms show concerted patterns of β-diversity, such as parallel community turnover or clustering, in response to changing selection regimes.
However, consistent holobiont patterns in natural populations are inherently difficult to characterize due to individual variation among hosts from the same population, as well as complex interactions between microbiomes and host ecology (e.g., Foster et al., 2017;Lozupone et al., 2012;Näpflin & Schmid-Hempel, 2018).One promising approach therefore involves the comparison of holobiont communities from simple, genetically distinct host organisms that only experience negligible ecological interferences other than hostsymbiont relations (Blankenchip et al., 2018;Goodrich et al., 2014;Reese & Dunn, 2018).
Importantly, lichens also typically lack confounding factors common to other host systems, such as effects of diet variation, herbivory, foodweb competition, or-in the case of rock-dwelling lichenscomplex and variable soil substrates (Goodrich et al., 2014;Näpflin & Schmid-Hempel, 2018).Diversity and distribution of individual components of the lichen holobiont are driven by the abiotic environment, as well as biotic interactions within the lichen thallus.
Fungal partners, for instance, can show intraspecific genomic differentiation along altitudinal clines (Dal Grande et al., 2017), and algal symbiont communities experience turnovers along latitude and altitude (Rolshausen et al., 2018(Rolshausen et al., , 2020;;Vargas Castillo & Beck, 2012;Werth & Sork, 2014).How the microbial communities associated with lichens respond to elevational gradients has not explicitly been addressed, but an effect of biogeographic region has been reported (Printzen et al., 2012;Sierra et al., 2020).Moreover, biotic factors that play a role in community assembly of a lichen holobiont include fungal host traits and photobiont identity.For example, the geographic distribution of algal symbionts is in part driven by fungal host identity (e.g., Dal Grande, Rolshausen, et al., 2018;Leavitt et al., 2015), and also the community structure of lichen-associated bacteria differs between closely related fungal host species, and between different types of photobionts associated with the host (Bates et al., 2011;Grube et al., 2009).While it is clear that composition and diversity of the lichen holobiont are driven by various biotic and abiotic factors, no study has yet addressed the combined responses of all major components of the lichen holobiont to an environmental gradient in a single system.
Open questions include whether different components of the holobiont respond to environmental gradients independently or in a concerted fashion, whether turnovers of holobiont components are gradual or abrupt, and whether explicit turnover conditions can be identified and geographically located.To address these knowledge gaps, we studied lichen holobiont composition (i.e., genomic differentiation of fungal hosts, abundance distributions of symbiotic green algae, and diversity patterns of bacterial communities) of two closely related host species along an elevational gradient spanning from Mediterranean to alpine climate regimes.We investigate to what extent host genomic differentiation and climatic transitions along the gradient coincide with algal symbiont distributions and compositonal changes in the bacterial communities.

| Sampling sites, host species, and algal symbionts
The two host species for which holobiont compositions were investigated are the rock-inhabiting lichen-forming fungi Umbilicaria pustulata and U. hispanica.Both taxa formerly belonged to the genus Lasallia (Davydov et al., 2017).These lichens form foliose, umbilicate thalli, and host algae of the genus Trebouxia as primary photosynthetic symbionts.Sampling was conducted along an elevational gradient on the south-facing slope of the Sierra de Gredos mountain range (Sistema Central, Spain), and included nine sampling sites ranging from c. 700 to 2100 m asl (Dal Grande, Rolshausen, et al., 2018).Both host species occur sympatrically in the mid-mountain zone (900-1800 m), but only U. pustulata extends into low-mountain regions (<900 m), and only U. hispanica extends into the highmountain regions (>1800 m); our sampling protocol thus included (i) a "mid-elevation overlap zone" where both host species co-occur at three sampling sites, (ii) a "low-elevation zone" where only U. pustulata occurs at three sampling sites, and (iii) a "high-elevation zone" where only U. hispanica occurs at three sampling sites (Figure 1).The sampled gradient encompasses three bioclimatic zones, spanning across Mediterranean ("mesomediterranean"), cool-temperate ("supramediterranean"), and alpine ("oromediterranean") with a range of annual precipitation from 400 to 850 mm and an annual mean temperature difference of 7.3°C between the lowest and highest sites.We summarized the climatic profile of the gradient based on 19 bioclimatic variables (Hijmans et al., 2005).A more detailed description of the sampling sites is given in Dal Grande, Rolshausen, et al. (2018).At each sampling site, we chose lichen thalli that were located on horizontal or gently sloping, fully sun-exposed rock faces.
We collected small samples (c. 8 mm in diameter) from 20 different lichen thalli per population over an area of c. 50 m² and with a minimal distance of 0.5 m between thalli, resulting in a total of 120 thallus fragments (individuals) across six sampling sites for each species (Figure 1).To minimize environmental contamination of thallus surfaces, all samples were taken with sterile tools, immediately transferred into sterile tubes, and thoroughly washed twice in sterile water before DNA extraction.Genomic DNA was then extracted from each thallus separately using the acetyl trimethylammonium bromide (CTAB)-based method described in Cubero and Crespo (2002).

| SNP genotyping of host species and characterization of algal symbionts
Genomic differentiation within each fungal host species was assessed based on a different sampling design including 100 sampling units (i.e., thallus fragments/individuals) per host per given altitude.Each population (100 individuals) was pooled resulting in F I G U R E 1 Schematic elevation profile showing the distribution of sampling sites (I to XI), host populations (h1: U. pustulata; h2: U. hispanica), microbial community cluster, and algal symbiont occurrence along the gradient (see also Dal Grande.Rolshausen et al., 2018).The basic climatic profile of the gradient-in reference to the bioclimatic framework of (Metzger et al., 2013)-reveals three main bioclimatic regions.The main climatic transition zones are located between 800 and 1000 m a.s.l.(Mediterranean to cool-temperate transition), and between 1800 and 1900 m a.s.l.(cool-temperate to alpine transition).°C refers to annual mean temperature and mm refers to annual precipitation.Indices I to IX refer to sampling sites six Pool-seq samples per host (Figure 1).Each pool contained equal amounts of DNA-measured with a Qubit fluorometer (dsDNA BR, Invitrogen)-from the 100 sampling units.We included six contamination controls, that were sequenced as well.Library preparation and sequencing on Illumina HiSeq2000 at ~90× coverage per pool were performed by GenXPro GmbH (Frankfurt am Main, Germany).Annotated reference genomes for each host were available for U. pustulata from (Dal Grande et al., 2017) and for U. hispanica from (Dal Grande, Meiser, et al., 2018).Trimmed and quality filtered reads (i.e., read length >80 bp and read quality >26) were mapped to the respective reference genome using BWA-MEM (Li & Durbin, 2009) with default parameter settings.We used SAMtools (Li, Handsaker, et al., 2009) to further process reads with a minimum mapping quality of 20 and single nucleotide polymorphisms (SNPs) were called using mpileup (Li, Li, et al., 2009).The software packages PoPoolation and PoPoolation2 (Kofler, Orozco-terWengel, et al., 2011;Kofler et al., 2011) were used to detect and mask indels (min-count 2, indelwindow 5) and to calculate pairwise genetic differentiation (F ST ) at coding regions (CDS) within hosts at a uniform coverage of 30 per population.For each host pairwise F ST comparisons were based only on coding SNPs with a minimum count of 4, a minimum quality of 20, and falling into the upper 0.5% tail of the overall F ST distribution.To visualize SNP-based F ST differentiation among populations for each host, we obtained a reduced set of pairwise distance matrices based on sample quantiles (0.975, 0.75, 0.5, 0.25, 0.025) of the full set of distances across all polymorphic SNPs.The resulting set of pairwise quantile distances was then jointly analysed using the DISTATIS generalization of classical multidimensional scaling (Abdi et al., 2005(Abdi et al., , 2009)).More details on genome annotation, SNP calling, and F ST analyses are given in Dal Grande et al. (2017).
Algal symbionts of the genus Trebouxia associated with the two host species were delineated via Illumina high-throughput metabarcoding of the ITS2 region (232 bp) of the rRNA operon (Coleman, 2009).Detailed descriptions of all molecular methods, bioinformatic analyses, and phylogenetic placements of the Trebouxia OTUs in our samples, along with an in-depth analysis of their community structure at all nine sampling sites, is given in Dal Grande, Rolshausen, et al. (2018).In the current study, we therefore did not attempt to reproduce the already published results, but rather to incorporate the results of Dal Grande, Rolshausen, et al. (2018) into a comprehensive discussion of elevational structuring of holobionts.Here, algal symbionts were thus described based on presence-absence patterns alone.Detailed analyses of algal community structure in U. pustulata and U. hispanica along the "Sierra de Gredos" gradient are provided by Dal Grande, Rolshausen, et al. (2018).

| Characterization of bacterial communities
Lichen-associated bacterial communities were characterized individually for each of the 240 thalli samples (i.e., 120 per host species).For this, we amplified the V4 domain of the 16S ribosomal RNA genes using the Illumina primer constructs of (Caporaso et al., 2011): 515f (GTGYCAGCMGCCGCGGTAA) and 806r (GGACTACNVGGGTWTCTAAT).Libraries of barcoded and normalized PCR products were constructed using the Fasteris PCR-free Metafast protocol and sequenced in three replicates on Illumina HiSeq in Rapid Run mode (2 × 300 bp) at FASTERIS, Switzerland (www.fasteris.com).Initial adapter and quality trimming was carried out with Trimmomatic (V0.32;Bolger et al., 2014), and adapter-based combinatorial demultiplexing was carried out using cutadapt (Martin, 2011).We used the high-resolution DADA2 method on pooled reads from all samples to infer amplicon sequence variants (ASV) for each sample (Callahan et al., 2016).Additional post-clustering curation was then carried out based on a conservative 97% similarity threshold, yielding operational taxonomic units (OTUs).Lastly, the LULU pipepline was applied with default parameters of similarity and cooccurence thresholds (Frøslev et al., 2017), to further remove highly similar artifacts among OTUs.The final curated set of OTUs was then taxonomically assigned based on the Silva SSU 138 reference database.

| Statistical analyses
To depict elevational patterns in bacterial α-diversity among sites, we calculated three indices-species richness, Shannon diversity, and Simpson diversity-along the full gradient, as well as within each host species.Sample completeness was evaluated via rarefaction and extrapolation (Chao et al., 2014;Hsieh et al., 2016; Figure S1), however, samples were not rarefied as recommended by McMurdie and Holmes (2014) and instead treated as compositional count data (Gloor et al., 2017).For each index we then used local estimated scatterplot smoothing (LOESS) to determine the shape of the responses.Changes in microbial community composition (β-diversity) were analysed (i) along the full gradient, comprising both host species at a total of nine sites, (ii) within each host species, and (iii) between host species, restricting analyses to the overlap zone where the two hosts co-occur in sympatry (i.e, between 1258 and 1699 m; Figure 1).For each set of turnover analyses (i.e, β-diversity for microbes and pairwise F ST within fungal host species), we primarily focused on detecting regions of significant directional turnover along the gradient and between the host species respectively, rather than describing turnover patterns between all possible pairwise combinations of sites (Legendre, 2014).This yielded a total of eight focal comparisons along the gradient (site I vs. site II; site II vs. site III, site III vs. site IV; site IV vs. site V, site V vs. site VI, site VI vs. site VII, site VII vs. site VIII, and site VIII vs. site XI), together with the hostspecific comparisons, and the comparison of pooled sites within the overlap zone (host 1 vs. host 2; Figure 1).Statistical associations between β-diversity and increasing elevation along the full gradient were tested via homogenous dispersion models between and within sampling sites based on distances to the respective group centroids (Anderson et al., 2006).Distances between centroids were tested using PERMANOVA (adonis in vegan; Oksanen et al., 2020).The latter was also used to depict effect sizes (η²) of β-diversity turnover for each of the aforementioned focal comparisons along the gradient.
Moreover, we tested altitude as a predictor of diversity differences along the full gradient; and the interaction of altitude and host species as predictors of diversity differences within the overlap zone.All distance-based models used Jaccard dissimilarity matrices.Analyses of β-diversity were conducted (i) for nonpartitioned Jaccard dissimilarity matrices, and in order to account for potential alpha-diversity effects, (ii) for the turnover component of the Jaccard dissimilarity only (Baselga, 2012).
We further scrutinized host-and climate-specific differences in bacterial communities by comparing differential relative OTU abundances (i) between hosts in the overlap zone and (ii) between Mediterranean low-altitude populations (sampling sites I-III) and alpine high-altitude populations (sampling sites (VII-IX).In this approach, an estimate of sampling error is obtained by taking Monte-Carlo samples from a Dirichlet distribution and used to transform OTU abundances based on centred log-ration transformation, yielding individual abundances that are relative to the mean abundance of all taxa in the respective sample (Fernandes et al., 2013(Fernandes et al., , 2014;;Gloor et al., 2016).Differences between and within hosts were then tested using Welch's t test with a Benjamini-Hochberg false discovery rate (FDR) correction.OTUs were considered as differentially abundant between hosts, or between genetic clusters within hosts, if their FDR adjusted p-values were significant (p < .05).We report differentially abundant OTUs at the taxonomic levels of "Order", "Family", and "Genus" (if yielded from taxonomic assignment); overall taxonomic profiles for all samples pooled for each sampling site along the gradient are presented at the level of "Order" and "Family" respectively.Finally, we use OTU prevalence at given detection thresholds to describe "core" communities (Salonen et al., 2012) that are present in the majority of samples across the entire gradient for each host species.Core communities are reported at the levels of "Order", "Family", and "Genus" (if yielded from taxonomic assignment).All the above statistical analyses were conducted in R 4.0.1 (R Core Team, 2020), using the following core packages: ALDEx2 (Fernandes et al., 2013), betapart (Baselga, 2012), DistatisR (Beaton et al., 2019), phyloseq (McMurdie & Holmes, 2013), and vegan (Oksanen et al., 2020).

| Genome-wide fungal host patterns
Both fungal host species showed a pronounced genomic differentiation that separated low elevation populations from high elevation populations.In the six U. pustulata (h1) populations we identified 165,623 polymorphic coding loci with a mean pairwise F ST of 0.148, in the six U. hispanica (h2) populations we found 204,313 polymorphic coding loci with a mean pairwise F ST of 0.
2).Likewise, richness estimates as well as α-diversity metrics (i.e., Simpson and Shannon diversity) showed decreasing trends along the gradient, both, within host species and for both host species combined (Figure 3).Across the region of overlapping occurrence (between c. 1,258 and 1,699 m a.s.l), we found a transitory diversity plateau (Figure 3), indicative of the concatenated diversity across both hosts.The analyses of "core" communities were based on a range of relative within-sample-abundance between 0.5 and 15% and a minimum prevalence across respective sample pools of 60% (Salonen et al., 2012).These parameters yielded a set of 13 bacterial OTUs in the core of U. pustulata (h1) and 10 bacterial OTUs in U.
hispanica (h2).Most of the OTUs overlapped between hosts, but we found five OTUs that were specific to the "core microbiome" of U.

pustulata (2 further unclassified Families Acetobacteraceae, 1 Genus
Granulicella, 1 Family Beijerinckiaceae, and 1 Genus Bryobacter), and two OTUs that were specific for the "core microbiome" of U. hispanica (1 further unclassified Family Acetobacteraceae, and 1 Genus Acidiphilium).See Figure S3 for a detailed description of both "core microbiomes".4).Analysis of bacterial community composition between host species could only be carried out between sampling sites IV, V, and VI-that is, where the two host ranges overlap.Here, we found significant compositional F I G U R E 3 α-diversity profile for host-associated bacterial communities along the elevational gradient, depicted as Hill numbers (i.e., Species richness, Shannon, and Simpson diversity).Dots indicate local sampling sites and LOESS smoothing curves determine the shape of the responses within host species (solid, blue and orange) and across hosts for the entire gradient (dashed, black).Orange =U.pustulata (h1), blue =U.hispanica (h2) differences between bacterial communities of U. pustulata versus U.

| Microbial β-diversity along the gradient and between host species
hispanica, with a between-host effect size that was comparable to the host-independent community turnover effects between climatic regions (F = 3.035, p < .050;multivariate ANOVA of homogenic dispersion between hosts: n h1 = 60, n h2 = 60; Figure 5).Moreover, multivariate variance partitioning within in the overlap zone revealed that the effect of host identity (F = 5.874, p < .001) was stronger compared to the effect of altitude (F = 2.057, p < .050),or the effect of the host × altitude interaction (F = 2.211, p < .050).Analyses of differential OTU abundances in communities relative to the mean OTU abundance in the respective sample pools yielded a set of 10 OTUs that showed significantly different abundance between hosts Marked turnover in host genetics and host-associated bacteria is not perfectly spatially synchronized, but the fungal genomic differentiation occurs one population below the bacterial community turnover F I G U R E 4 NMDS ordinations of Jaccard dissimilarity depicting β-diversity for host-associated bacterial communities (a) along the full gradient, comprising both host species at a total of nine sites, and (b) between host species (h1 and h2), restricting analyses to the overlap zone and pooling sampling sites where the two hosts co-occur in sympatry (i.e., sites IV, V, and VI between 1 258 and 1,699 m; Figure 1).Sampling sites are summarized with 95% confidence ellipses and the altitude-dependent climatic profile of the respective ordination space is shown as a heatmap insert with precipitation (mm) and temperature (°C) gradients depicted as colour gradients (see Figure S6 for a detailed climatic resolution of the ordination space).Effect sizes of directional turnover (η²), based on permutational ANOVAs of focal comparisons along the gradient and between hosts, are shown in (c).Colours and roman numerals correspond to Figure 1 (a) F I G U R E 5 β-diversity for hostassociated bacterial communities within hosts (h1 and h2) along the gradient.
Left panel shows NMDS ordinations of Jaccard dissimilarity with sampling sites summarized with 95% confidence ellipses.Right panel depicts effect sizes of directional turnover (η²), based on permutational ANOVAs of focal comparisons within hosts.Colours and roman numerals correspond to Figure 1 in each of the two hosts (Figures 1 and 4).Second, a pronounced community difference between the parapatric and sympatric (overlap) occurrence regions (Figure 5).Notably, the above results did not change when analyses were conducted based on the turnoverpartition of the Jaccard dissimilarity matrix, accounting for potential effects of α-diversity (Baselga, 2012; see also Figure S5).

| Green algal symbiont distributions
We detected a total of five green algal OTUs, belonging to the genus Trebouxia, along the gradient (Figure 1).OTU1 is a climate and host generalist occurring in all populations in both hosts.OTU2 is a host generalist, occurring in h1 and h2 regularly, but only below 1700 m.
OTU3 is found solely in the two bottommost populations of the gradient (Mediterranean climate), and associates only with the low elevation genotype of h1.OTU4 and OTU5 are high elevation OTUs and occur only with one of the hosts (h2) in the present study, but they may also associate with h1 in cold climates, as shown in a rangewide investigation of h1 (Rolshausen et al., 2018).Rolshausen et al.
(2018) identified a Trebouxia ecotype preferring warm and dry climate, which is identical to the present OTU3, and an ecotype preferring cool and wet climate, which corresponds to the present OTU4.

| DISCUSS ION
Holobionts are dynamic biological entities.Like ecosystems, they exhibit species networks, trophic levels, and responses to abiotic factors.The functional, ecological, and evolutionary consequences of habitat-specific beta-diversity among members of holobionts are not well understood.In part, this is due to a lack of suitable model systems that allow the comparison of holobiont communities from simple, genetically distinct host organisms for which ecological interferences other than host-symbiont relations can be neglected (Blankenchip et al., 2018;Goodrich et al., 2014;Reese & Dunn, 2018).We here contribute to the discussion by showing that lichen holobiont diversity is structured along an elevational gradient.In particular, we identify two major transition zones for two lichen holobionts along an elevational gradient, which are not detectable in morphological or chemical traits of the studied host species (see Singh et al., 2021;IS personal observation, October 2018).We show for the first time that concerted turnover occurs at different organizational levels of lichen holobionts: in the genomic differentiation of fungal hosts, in the presence/absence of certain symbiotic green algae, and in bacterial communities.
The first biogeographic transition zone is located between 800 and 1200 m a.s.l.(average annual temperature 11°C-13°C), enveloping the changeover from Mediterranean to cool-temperate climate.
This first zone marks a significant genome-wide differentiation in the fungal host U. pustulata (h1) as well as the lower distribution limit of fungal host U. hispanica (h2).Together with these boundaries within and between host species, we detect a marked low elevation cluster of bacterial communities associated with U. pustulata (Figure 4).The annual average temperatures (~1214°C) (Vančurová et al., 2021).
The second transition zone is located between 1500 and 1800 m a.s.l., enveloping the changeover from cool temperate to alpine climate.Here, we detect significant genetic differentiation within fungal host h2, the upper distribution limit of h1, and a compositional shift in the bacterial communities of h2.Analogously to the transition zone at lower altitude, bacterial turnover and host genomic turnover do not exactly coincide, but are offset by one population, suggesting that these two lichen components respond independently to conditions along the gradient.The second transition zone also marks the upper distribution limit of algal OTU2.Since this taxon is widespread in both hosts at lower altitudes, climatic factors-rather than host effects-probably prevent this alga from occurring in alpine habitats.
Similarly to vascular plant species that cannot cross the tree line, this algal photobiont apparently cannot persist under alpine conditions.In summary, the two transition zones identified here mark a clustering of distribution limits and abrupt community shifts within the lichen holobiont.Structuring takes place both within hosts and across the full gradient.Our results therefore suggest that the two transitions zones constitute ecological boundaries, which divide the lichen holobiont into three altitudinal diversity classes.If we view holobionts as ecosystems, that is, healthy when in equilibrium (Pita et al., 2018), each altitudinal class with its specific constellation of host genotypes and associated microorganisms (algae and bacteria), can be considered a "stable state" of the lichen ecosystem, a dynamic equilibrium that arises under the given environmental conditions and allows proper functioning of the holobiont.
The observation that the composition of microbial communities (and potentially their functionality) changes according to ecological and climatic variations, suggests a link to the adaptive potential of the lichen holobiont.Likewise, the adaptation of lichen populations to new habitats is presumably accompanied by changes in the bacterial communities-much similar to previously observed switches of the photobiont strains that correlate with different ecological conditions (Blaha et al., 2006;Grube & Berg, 2009;Rolshausen et al., 2018Rolshausen et al., , 2020)).Moreover, as also suggested by our results (Figures S3 and S4), comparative analyses of lichens from different environmental backgrounds have described the presence of an environment-independent "core" microbiome that coexists together with distinct environment-specific fractions of the total microbiome (Aschenbrenner et al., 2016;Grimm et al., 2021;Hawksworth & Grube, 2020).Further evidence suggesting a strong link between the adaptive potential of lichens and the variation in their microbiomes comes from studies that examined microbiome responses to long-term climate warming scenarios (Klarenberg et al., 2020).Interestingly, the transition zones that we find to structure the individual components of lichen holobionts in our study coincide with major ecological zones (biomes) known to structure communities of macroorganisms.It is thus possible that these areas encompass environmental conditions that act as universal environmental boundaries for macroorganisms and microorganisms alike.
Contrary to many horizontally acquired symbioses in which host genotype has only limited to no effect on microbiome composition (e.g., luminescent bacteria in the Hawaiian bobtail squid [Rader & Nyholm, 2012], the gut microbiota of laboratory mice [Hufeldt et al., 2010], insects [Kikuchi et al., 2007], a montane amphibian [Griffiths et al., 2018], and the human gut microbiome [Rothschild et al., 2018;Spor et al., 2011]), our results indicate the presence of a significant, albeit small fraction of host-specific bacteria, detected in the overlap zone of both host species.Recent studies in Arabidopsis thaliana showed that only a minor fraction of the phyllosphere microbiome is driven by host genotype and abiotic factors, but that this small fraction is ecologically important, because it may trigger cascades of microbe-microbe interactions that are critical for maintaining microbiome community structure (Agler et al., 2016;Banerjee et al., 2018).Hence, although the host-specific fractions probably comprised a rather minor portion of the lichen-associated bacteria, their ecological roles in the symbiosis need not to be underestimated.One promising way forward is to integrate current results with future reciprocal transplant experiments to further disentangle host genetics from potential abiotic effects and to identify specific fractions of the bacterial community involved in host ecology and evolution.
To conclude, we here show for the first time that lichen holobionts behave as ecological units in similar ways as macroecological communities along elevational gradients in that the holobiont subunits all exhibit clusters of diversity separated by two regional zones that broadly correspond to biogeographic transitions from Mediterranean to cool climates, and from cool to alpine climates.
Locating such discontinuities in microbial as well as host diversity is a promising step towards more integrated concepts of gradient biogeography, and can further help to detect and precisely locate climate sensitive areas (Angeler et al., 2016).Our study contributes to understanding ecological drivers of the lichen symbiosis (Vančurová et al., 2021), and specifying climatic niches of individual lichen components, important prerequisites for projecting future distributions of lichens (Nelsen et al., 2022).
156.The strongest directional pairwise differentiation (see Methods) along the gradient in U. pustulata was found between populations II (887 m) and III (1082 m), corresponding to the transition zone between Mediterranean and cool-temperate climate regions.The four uppermost U. pustulata populations (III-VI, 1082-1699 m) clustered together with no significant pairwise differentiation (Figure 2a-b), and the two lowermost U. pustulata populations (I and II; Figure 1) grouped together showing only minor differentiation.A similar pattern was found for U. hispanica with the strongest directional differentiation located between populations V (1480 m) and VI (1699 m), dividing U. hispanica populations into a high elevation group (populations VI-IX, 1699-2096 m), and a low elevation group (populations IV-V, 12,581,480 m).A pronounced genomic differentiation in U. hispanica thus resides at the upper end of the cool-temperate climate region (Figure 2c-d).

F
Genome-wide differentiation in coding regions of both host species along the gradient.(a-b) U. pustulata, (c-d) U. hispanica.Left panels show Distatis compromise configurations based on pairwise F ST quantile distances between the populations of each host with ellipses depicting 95% confidence.Right panels show effect sizes of linear directional differentiation along the gradient, depicted as pairwise F ST differences between neighbouring populations of each host species For three replicated libraries (r1-r3), we analysed a total of 5.498.248(r1), 5.809.732(r2), and 6.221.131(r3) raw paired-end reads.After quality filtering, adapter-trimming, and demultiplexing, DADA2 detected 1351 ASVs in r1, 1,392 ASVs in r2, and 1,407 ASVs in r3, which were then summed up and merged to a total of 2001 ASVs that went into downstream analyses.After conservative 97% similarity clustering (yielding 1020 OTUs) and further merging (via LULU pipeline, yielding 558 OTUs) we assigned taxonomy identification based on the SILVA SSU 138 reference database.Another 162 OTUs were then excluded based on Genus-level taxonomy (i.e., chloroplasts, Eukaryotes, and unclassified), resulting in a total of 396 bacterial OTUs.Lastly, we subset our data based on a conservative OTU abundance threshold of 0.005% (Bokulich et al., 2013), resulting in a final data set of 217 bacterial OTUs.The taxonomic profiles for the most abundant (>1% OTU abundance) Orders and Families from the final OTU set are depicted in the Supporting Information (Figures S2.1, S2.2): U. pustulata (h1) and U. hispanica (h2) holobionts along the studied gradient were dominated by the following microbial Orders: Acetobacterales, Acidobacteriales, and Armantimonadales; and by the microbial Families Acetobacteraceae and Acidobacteriaceae.Correlative analyses between elevation and taxa abundance showed a pattern of decreasing relative abundances from low to high altitudes in both host species.Detailed results of abundance-elevation correlations are included in the Supporting Information (Figures S2.1, S2 Analysis of the bacterial OTU communities in holobionts-combined for both host species-revealed three prominent altitudinal classes that were separated by significant directional turnover of community composition along the gradient (F = 2.208, p < .001;multivariate ANOVA of homogenic dispersion among sampling sites).Compared to random effect size expectations, community turnover was most pronounced between sampling sites II and III, sampling sites III and IV, and sampling sites VI to VII (Figure4).The low-elevation class comprised sampling sites I to III between c. 700 and 1100 m a.s.l., spanning from the Mediterranean climatic region into the lower cool-temperate climatic region; the mid-elevation class comprised sampling sites IV to VI between c. 1200 and 1700 m a.s.l., spanning the cool-temperate climatic region; and the high-elevation class comprised sampling sites VII to IX between c. 1800 and 2100 m a.s.l., spanning from the upper cool-temperate climatic region into the alpine climatic region (Figures 1 and h1 and h2 in the overlap zone (i.e., sites IV-VI; Figure S4.1), and a set of 27 OTUs that showed significantly different abundance between low elevation sites (I-III) and high elevation sites (VII-IX).Taxonomic details and distinct relative abundances for all significant OTUs are shown in Figures S4.1 and S4.2.Microbial community turnover within host species (h1 and h2) revealed two main patterns: First, a pronounced compositional difference between the two genetically differentiated host-populations (low vs. high elevation) in U. pustulata and U. hispanica, respectively.
sharp transition in bacterial community composition between populations III (1080 m) and IV (1260 m) is geographically close, but not identical to the intraspecific genomic differentiation of the fungal host (h1) and the disappearance of algal OTU3 between populations II (800 m) and III (1080 m), suggesting abiotic drivers are shaping bacterial communities, rather than fungal genotype or dominant symbiotic alga (Figure 5).Notably, other species of lichenized fungi with distributions transgressing the Mediterranean-cold temperate biome boundary also show algal symbiont turnovers at comparable