Bacteriomic Profiles of Rock-Dwelling Lichens from the Venezuelan Guiana Shield and the South African Highveld Plateau

Lichens are not only fungal–algal symbiotic associations but also matrices for association with bacteria, and the bacterial diversity linked to lichens has been receiving more attention in studies. This study compares the diversity and possible metabolism of lichen-associated bacteria from saxicolous foliose and fruticose taxa Alectoria, Canoparmelia, Crocodia, Menegazzia, Usnea, and Xanthoparmelia from the Venezuelan Guiana Shield and the South African Highveld Plateau. We used DNA extractions from the lichen thalli to amplify the eukaryotic 18S rRNA gene (rDNA) and the V3–V4 region of the bacterial 16S rDNA, of which amplicons were then Sanger- and MiSeq-sequenced, respectively. The V3–V4 sequences of the associated bacteria were grouped into operational taxonomic units (OTUs) ascribed to twelve bacterial phyla previously found in the rock tripe Umbilicaria lichens. The bacterial OTUs emphasized the uniqueness of each region, while, at the species and higher ranks, the regional microbiomes were shown to be somewhat similar. Nevertheless, regional biomarker OTUs were screened to predict relevant metabolic pathways, which implicated different regional metabolic features.


Introduction
Lichens are the symbiotic associations of fungal mycobionts and algal/cyanobacterial photobionts and play essential ecological roles in many ecosystems.They can be primary colonizers of bare rock, soil, or wood, initiating the process of soil formation and nutrient cycling.They can also contribute to ecosystems via nitrogen fixation when cyanobacteria participate in symbiosis and serve as food sources for various animals.In addition, lichens are indicators of air pollution, climate change, and a variety of other environmental changes [1,2].
As pioneer organisms, lichens are often the first to colonize bare rocks.They can do this because they are able to survive in extreme environments such as deserts, tundra, and bare rock surfaces [3].Of the estimated 5 million fungal species [4], only 156,287 have been included in the Species Fungorum (as of 12 December 2023) [5], and 19,387 are lichen-forming species [6], of which ca. 10-20% are regarded as rock-dwelling or epilithic lichens [1,7,8].
Bacterial communities associated with lichens have been studied from the holobiont viewpoint [9], especially with the advent of multi-omics and high-throughput sequencing techniques [10,11].Studies on lichen-associated microbiomes often include or target epilithic lichens and report the bacterial families of Acetobacteraceae, Acidobacteriaceae and Table 1.Sampling sites of epilithic lichens inhabiting the rocks in the Venezuelan Chiuri Tepui and the South African Highveld Plateau.Coordinates (latitudes and longitudes) and elevations (altitudes) were determined with GPSMAP62S (Garmin, Olathe, KS, USA).

Region Sample Code Latitude Longitude Altitude
Churi Tepui, Guiana Shield, Venezuela G01 05 The lichen samples from the South African Highveld Plateau were collected in October 2018 from rocks and cliffs in grassland and bush along a stream in the Golden Gate Highlands National Park (ca.28 • 52 ′ S, 28 • 60 ′ E; Figure 1, Table 1).The park is situated in the foothills of the Maluti Mountains and is characterized by high-altitude grasslands, rolling hills, valleys, and sandstone cliffs.The park has a high-altitude climate, with cool temperatures and low humidity.The average temperature ranges from 13 • C to 26 • C in summer and from 1 • C to 15 • C in winter.Rainfall is concentrated in the summer months from November to February, with an average annual rainfall of approximately 650 mm to 760 mm [28,29].
All Guiana Shield lichens (G01 to G12) and half of the South African lichens (SA01, SA03, SA05 and SA07) were collected from bare rocks in grassland, but the other half (SA02, SA04, SA06 and SA08) were collected from bare rocks in bush.Thalli of epilithic lichens were cut with a flamed field knife and put into Whirl-Pak bags (Nasco, Fort Atkinson, WI, USA).The lichen thalli samples were air-dried, stored in the dark on site, transferred to the laboratory at Hiroshima University, and frozen at −25 The lichen samples from the South African Highveld Plateau were collected in October 2018 from rocks and cliffs in grassland and bush along a stream in the Golden Gate Highlands National Park (ca.28°52′ S, 28°60′ E; Figure 1, Table 1).The park is situated in the foothills of the Maluti Mountains and is characterized by high-altitude grasslands, rolling hills, valleys, and sandstone cliffs.The park has a high-altitude climate, with cool temperatures and low humidity.The average temperature ranges from 13 °C to 26 °C in summer and from 1 °C to 15 °C in winter.Rainfall is concentrated in the summer months from November to February, with an average annual rainfall of approximately 650 mm to 760 mm [28,29].
All Guiana Shield lichens (G01 to G12) and half of the South African lichens (SA01, SA03, SA05 and SA07) were collected from bare rocks in grassland, but the other half (SA02, SA04, SA06 and SA08) were collected from bare rocks in bush.Thalli of epilithic lichens were cut with a flamed field knife and put into Whirl-Pak bags (Nasco, Fort Atkinson, WI, USA).The lichen thalli samples were air-dried, stored in the dark on site, transferred to the laboratory at Hiroshima University, and frozen at −25 °C in preparation for bulk DNA extraction.Relief Model [30] is the source of the map image.

Bulk DNA Extraction from Lichen Thalli
Thalli of a lichen sample were cleaned with autoclaved Milli-Q ultrapure water, cut into pieces, ground to finer fragments, and homogenized, one gram of which was used for DNA extraction by the method detailed in previous studies [20,21].Although they were not examined under a microscope, cephalodia were not clearly visible.The extracted DNA was maintained at −20 °C until PCR amplification.

Amplification and Sequencing of Fungal/Algal 18S rDNA
Near-full-length 18S rDNA of lichen-forming fungi and algae were amplified on two TaKaRa Thermal Cyclers with the primer sets shown in Table 2 and using the thermal cycling described in the previous study [21].The PCR amplicons were purified and Sanger-sequenced at the Department of Gene Science, Natural Science Center for Basic Research and Development (N-BARD), Hiroshima University [21].Relief Model [30] is the source of the map image.

Bulk DNA Extraction from Lichen Thalli
Thalli of a lichen sample were cleaned with autoclaved Milli-Q ultrapure water, cut into pieces, ground to finer fragments, and homogenized, one gram of which was used for DNA extraction by the method detailed in previous studies [20,21].Although they were not examined under a microscope, cephalodia were not clearly visible.The extracted DNA was maintained at −20 • C until PCR amplification.

Amplification and Sequencing of Fungal/Algal 18S rDNA
Near-full-length 18S rDNA of lichen-forming fungi and algae were amplified on two TaKaRa Thermal Cyclers with the primer sets shown in Table 2 and using the thermal cycling described in the previous study [21].The PCR amplicons were purified and Sangersequenced at the Department of Gene Science, Natural Science Center for Basic Research and Development (N-BARD), Hiroshima University [21].

Amplification and Sequencing of V3-V4 Region of Bacterial 16S rDNA
Using the same prepared DNA, PCR amplification of the V3-V4 region of 16S rDNA was carried out with the specific primers 341F and 806R (Table 2).The thermal cycling conditions for the PCR were as follows: 95  C for 5 min.
To construct the sequence library and conduct paired-end 300 bp sequencing, the molecular diagnostic company, Environmental Research and Solutions Co. Ltd. (Kyoto, Japan), was employed [21].

Sequence Data Analysis and OTU Determination
The 18S rDNA sequences produced by Sanger sequencing were aligned using ClustalW and BioEdit to remove poor-quality sequences [31,32].Then, the remaining sequences were manually assembled and checked for chimeras using tree topology analysis [33].Finally, the resulting sequences were searched using BLAST to identify the lichen-forming fungi and algae.
The sequences of 18S rDNA for lichen-forming fungi and algae have been deposited in the DDBJ/ENA/GenBank database with accession numbers ranging from LC761218 to LC761237 and LC761244 to LC761263, respectively.The V3-V4 reads are accessible at the DDBJ Sequence Read Archive (DRA) and can be found under the accession num-berDRA015994.The associated BioProject and BioSample numbers are PRJDB15406 and SAMD00585845 to SAMD00585864, respectively, and the sample-to-number correspondence can be found in Table S1.

Diversity Indices and Bioinformatic Analyses of OTUs
The MTP pipeline of EzBioCloud was used to analyze the rarefaction curves.Specifically, using the alpha diversity indices (Chao1, Shannon, and Simpson indices), the richness and evenness of bacterial OTUs associated with the lichen samples were estimated by EzBioCloud MTP.Chao1 corresponds to a rarefaction curve asymptote as an estimator of species richness or OTU richness [36].Shannon and Simpson indices were also used to calculate an "effective number of species" (ENS) [37] or an effective number of OTUs of a sample.It is important to note that the Chao1 index takes singletons into account.
For beta diversity, the OTUs were clustered based on the UniFrac distance matrix [38], and biomarker OTUs were screened using the linear discriminant analysis (LDA) [39] and LDA-effect size method (LEfSe) [40] with LDA scores of 4.0 and 4.5 as the thresholds to screen, respectively.Differential abundance was analyzed using the analysis of compositions of microbiomes with bias correction (ANCOM-BC) [41].

Identification of Mycobionts and Photobionts of the Epilithic Lichens
All of the mycobionts of the epilithic lichens from the Guiana Shield and Highveld Plateau were attributed to genera of the class Lecanoromycetes (Table 3), the largest class among the lichen-forming fungi [45].Similarity values were at least 98.58% (Table S2).The samples G07-G12 from the Guiana Shield and SA01-SA02 from the South African Highveld Plateau were affiliated with the same lichen species, Canoparmelia caroliniana (Nyl.)Elix and Hale, known to occur in the Americas and East Africa [46].The other samples were affiliated with different species by region.Fruticose lichens from the Guiana Shield were affiliated with the species Alectoria sarmentosa (Ach.)Ach.(G01, G02) and Usnea florida (L.) Weber ex F.H.Wigg.(G03, G04).The samples G05 and G06 were affiliated with Menegazzia terebrata (Hoffm.)A.Massal., a sub-cosmopolitan (excluding Antarctica) lichen [47].The South African samples SA03-SA05 were affiliated with sub-cosmopolitan species Xanthoparmelia conspersa (Ehrh.ex Ach.) Hale [48].The samples SA06-SA08 were affiliated with the only peltigeracean species in this study, Crocodia aurata (Ach.)Link (also a cosmopolitan species) [49].
All of the algal photobionts were most closely related to the green algal species ascribed to the genus Trebouxia, the most prevalent photobiont among lichens [50].Similarity values were 98.44% or higher (Table S3).However, paraphyly in the genus Trebouxia and delineation of new genera Asterochloris and Vulcanochloris leaves room for considering members of the new genera as photobionts of the studied lichens [51,52].
Cyanobacterial OTUs were present in all of the samples.A total of 107 cyanobacterial OTUs were detected, of which 82 were affiliated with species or species-rank taxa (Table S4).The ratio of cyanobacterial reads to the total read number of a sample ranged from 0.11% in G01 with 5 OTUs to 13.26% and 11.01% in G05 and G04 with 28 and 19 OTUs, respectively.The highest ratio for a single OTU was 5.80% of "PAC002560_g_uc" (genusrank, uncultured, details unknown) in G12, followed by 5.31% of "JN023297_s" (speciesrank) [53] in G04.The same OTU "JN023297_s" was present at a ratio as high as 2.51% in G04 and in other samples at <1%.Other OTUs occurring at >1% were "PAC000112_s" in G07 and G11 at 3.09% and 3.06%, respectively; "DQ914863_g_uc" [54] in G05 and G08 at 2.17% and 1.10%, respectively; "FJ465967_s" [55] in G05 at 2.14%; "Stigonema ocellatum" in G05 and G03 at 1.41% and 1.16%, respectively; and "AY326529_s" [56] in G06 at 1.10%.The OTU affiliated with Stigonema ocellatum was also detected in the Antarctic epilithic lichen Umbilicaria [20].The cyanobacteria represented by the OTUs with read frequencies > 1% may function as additional photobionts in tripartite lichens or multipartite lichens, as in the case of G05 with its two > 1% cyanobacterial OTUs; however, inclusion of cyanobacteria that are not associated with lichens but are nearby may not be ruled out (discussed later).
Table 3. Taxonomic classification and the closest species based on near-full-length 18S rDNA sequences of mycobiont fungi of the epilithic lichens from the Venezuelan Guiana Shield (G01 to G12) and the South African Highveld Plateau (SA01 to SA08).All of the listed taxa belong to the class Lecanoromycetes.

MiSeq-Generated V3-V4 Sequences and OTUs
A total of 1,010,008 raw reads from the 20 lichen samples were generated using Illumina MiSeq sequencing, which were then filtered to 848,814 valid paired reads to be grouped into OTUs.Based on the analysis records in the EzBioCloud database [34], the mean length of all valid reads was 403.8 bp.For Guiana and South Africa, the mean length of valid reads was 401.8 bp and 406.8 bp, respectively.
Rarefaction curves were generated using the read and OTU counts (Figure S1).The coverage of rarefaction analysis expressed the ratio of obtained OTUs (Table 4) against the estimated total OTUs (a rarefaction asymptote in Figure S1), the latter of which is equivalent to the alpha diversity index, Chao1 (shown later).The mean, minimum, and maximum coverage ratios were 95.42%, 87.24% (in G02), and 99.67% (in G06), respectively.The coverage ratios suggest that the valid reads generated in this study are sufficient for further statistical and bioinformatic analyses.Table 5 shows the regional numbers of taxa (OTU, species, genus, family, order, class, and phylum) detected only in the Guiana lichen samples, the South African samples, and both regions' samples.The observed bacterial OTUs showed higher percentages of regionspecific features.However, regional traits were somewhat ambiguous at the species and higher ranks, with the regions' common OTUs being more than half of total OTUs at the order, class, and phylum ranks.The results emphasize the uniqueness of each region at the OTU rank and the similarity of the two areas at higher ranks.

Taxonomic Composition of Lichen-Associated Bacterial Community
Compositions of the OTU-derived bacterial phyla in 20 lichen samples are shown in Figure 2. A total of 12 bacterial phyla are portrayed as the standard features in all 20 samples.Each lichen sample contained 10 to 23 bacterial phyla (Table 4), including the 4-13 phyla consisting of less than 1% (of total) reads in each sample.The most common were Acidobacteriota, Actinomycota, Armatimonadota, Bacteroidota, Chloroflexota, Cyanobacteria, Deinococcota, Gemmatimonadota, Planctomycetota, Pseudomonadota, Saccharibacteria_TM7 and Verrucomicrobiota.

Taxonomic Composition of Lichen-Associated Bacterial Community
Compositions of the OTU-derived bacterial phyla in 20 lichen samples are shown in Figure 2. A total of 12 bacterial phyla are portrayed as the standard features in all 20 samples.Each lichen sample contained 10 to 23 bacterial phyla (Table 4), including the 4-13 phyla consisting of less than 1% (of total) reads in each sample.The most common were Acidobacteriota, Actinomycota, Armatimonadota, Bacteroidota, Chloroflexota, Cyanobacteria, Deinococcota, Gemmatimonadota, Planctomycetota, Pseudomonadota, Saccharibacteria_TM7 and Verrucomicrobiota.At the family level, the overall top five families in this study were acidobacterial Acidobacteriaceae and Bryobacteraceae, and alphaproteobacterial Acetobacteraceae, Beijerinckiaceae, and Sphingomonadaceae (Figure S4).These are compared with microbiomes of Thai tropical lichens, whose top five families are Beijerinckiaceae, Chthoniobacteraceae At the family level, the overall top five families in this study were acidobacterial Acidobacteriaceae and Bryobacteraceae, and alphaproteobacterial Acetobacteraceae, Beijerinckiaceae, and Sphingomonadaceae (Figure S4).These are compared with microbiomes of Thai tropical lichens, whose top five families are Beijerinckiaceae, Chthoniobacteraceae (phylum Verrucomicrobiota), Acetobacteraceae, Gemmataceae (phylum Planctomycetota), and an unidentified family in the order Tepidisphaerales (phylum Planctomycetota) [57].The difference in the top five families can be ascribed to biogeography as well as to host lichen species.Of ten Thai lichens, three and one belong to the families Parmeliaceae and Peltigeraceae, respectively, but all belong to genera that are different from those in this study.

Alpha and Beta Diversity
Alpha diversity indices, Chao1, Shannon, and Simpson (Table 6), were used to calculate the effective number for species (ENS) [37].Chao1, Shannon, ENS values and observed OTU numbers have positive correlations with species and evenness, and Simpson index values negatively correlate with species and evenness.Therefore, higher Chao1, Shannon, ENS values, observed OTU numbers, and lower Simpson index values found in South African samples indicate higher species richness and evenness.Due to different calculation methods, Shannon and Simpson indices could not be used to estimate the richness of bacterial species.Comparatively, Chao1 values were close to estimated OTU numbers and may better represent species richness of a large sample size, as reported in other studies [25,58].
Beta diversity analysis demonstrated regional separation between the Guiana Shield (G01 to G12) and South Africa (SA01 to SA08) at the species rank (Figure 3) and the genus and class ranks (Figure S6).Regional separation is unclear at the family, order, and phylum ranks (Figure S6), which may be influenced by the "common" taxa at these ranks between the two regions (Table 5).Weak intra-regional separation or intra-regional variation between the bush samples (SA02, SA04, SA06 and SA08) and the grassland samples (SA03, SA05 and SA07) is implied; however, grouping of grassland SA01 with bush samples obscures the intra-regional variation.Nevertheless, monospecific intra-regional variation is also seen in the Guiana Shield samples, such as those between G05 and G06 of Menegazzia terebrata as well as those between G07 and G12 of Canoparmelia caroliniana (Figures 2 and 3).
(SA03, SA05 and SA07) is implied; however, grouping of grassland SA01 with bush samples obscures the intra-regional variation.Nevertheless, monospecific intra-regional variation is also seen in the Guiana Shield samples, such as those between G05 and G06 of Menegazzia terebrata as well as those between G07 and G12 of Canoparmelia caroliniana (   7, were selected by setting the LDA score threshold to 4.5.Substantial biomarkers for the Venezuelan Guiana Shield included two taxa, the family Acetobacteraceae and the order Rhodospirillales, both affiliated with the phylum Pseudomonadota.Biomarkers for the South Africa Highveld Plateau included seven taxa: the genus Sphingomonas, the family Sphingomonadaceae, the order Sphingomonadales, the genus-raked EU289441_g of the family Beijerinckiaceae, and the order Frankiales affiliated with the class-rank Actinomycetota_c of the phylum Actinomycetota.No biomarkers at the species rank were identified, with a threshold value of 4.5.The regional distinctions of OTUs are presented in the phylogenetic cladogram (Figure 4).Substantial biomarkers, listed in Table 7, were selected by setting the LDA score threshold to 4.5.Substantial biomarkers for the Venezuelan Guiana Shield included two taxa, the family Acetobacteraceae and the order Rhodospirillales, both affiliated with the phylum Pseudomonadota.Biomarkers for the South Africa Highveld Plateau included seven taxa: the genus Sphingomonas, the family Sphingomonadaceae, the order Sphingomonadales, the genus-raked EU289441_g of the family Beijerinckiaceae, and the order Frankiales affiliated with the class-rank Actinomycetota_c of the phylum Actinomycetota.No biomarkers at the species rank were identified, with a threshold value of 4.5.No biomarker OTU was identified at the species rank with an LDA score >4.5 (Table 7).However, by reducing the threshold to 4.0, six biomarkers, two from the Guiana Shield and four from the Highveld Plateau, were detected at the species rank and were subject to differential abundance analysis by ANCOM-BC.The most substantial biomarkers, nitrogen-fixing Beijerinckia mobilis (order Hyphomicrobiales, family Beijerinckiaceae) for Guiana Shield and EU289441_s affiliated with beta-carotene-producing Lichenibacterium (order Hyphomicrobiales, family Lichenibacteriaceae [59]) for South African Highveld Plateau, are shown in Figure 5, with other biomarkers shown in Figure S7.No biomarker OTU was identified at the species rank with an LDA score >4.5 (Table 7).However, by reducing the threshold to 4.0, six biomarkers, two from the Guiana Shield and four from the Highveld Plateau, were detected at the species rank and were subject to differential abundance analysis by ANCOM-BC.The most substantial biomarkers, nitrogenfixing Beijerinckia mobilis (order Hyphomicrobiales, family Beijerinckiaceae) for Guiana Shield and EU289441_s affiliated with beta-carotene-producing Lichenibacterium (order Hyphomicrobiales, family Lichenibacteriaceae [59]) for South African Highveld Plateau, are shown in Figure 5, with other biomarkers shown in Figure S7.
Six biomarkers at species rank with LDA > 4 were further predicted for metabolic pathways of lichen-associated bacteria.At the KEGG Level 1, i.e., the highest metabolic categories on the KEGG database, these biomarker OTUs were related by relative abundances to the following five significant pathways: "Metabolism", "Genetic information processing", "Unclassified", "Environmental information processing", and "Cellular processes".The "Metabolism" pathway showed the highest relative abundance, as high as over 50% in OTUs from both sampling regions (Figure S8).
At Level 2, i.e., sub-metabolic categories on the KEGG database, each biomarker OTU was related to the 25 pathways (Figure S9), of which the top five were carbohydrate metabolism (10.32% in the Guiana Shield samples and 10.71% in the South African samples), membrane transport (10.30% in Guiana and 11.13% in South Africa), amino acid metabolism (10.20% in Guiana and 10.58% in South Africa), replication and repair (7.55% in Guiana and 7.14% in South Africa) and energy metabolism (7.11% in Guiana and 6.39% in South Africa).The most substantial difference was identified in "membrane transport," which was more dominant in the South African OTUs.However, the greatest difference was only 0.83%.Six biomarkers at species rank with LDA > 4 were further predicted for metabolic pathways of lichen-associated bacteria.At the KEGG Level 1, i.e., the highest metabolic categories on the KEGG database, these biomarker OTUs were related by relative abundances to the following five significant pathways: "Metabolism", "Genetic information processing", "Unclassified", "Environmental information processing", and "Cellular processes".The ʺMetabolismʺ pathway showed the highest relative abundance, as high as over 50% in OTUs from both sampling regions (Figure S8).
At Level 2, i.e., sub-metabolic categories on the KEGG database, each biomarker OTU was related to the 25 pathways (Figure S9), of which the top five were carbohydrate metabolism (10.32% in the Guiana Shield samples and 10.71% in the South African samples), membrane transport (10.30% in Guiana and 11.13% in South Africa), amino acid metabolism (10.20% in Guiana and 10.58% in South Africa), replication and repair (7.55% in Guiana and 7.14% in South Africa) and energy metabolism (7.11% in Guiana and 6.39% in South Africa).The most substantial difference was identified in ʺmembrane transport,ʺ which was more dominant in the South African OTUs.However, the greatest difference was only 0.83%.
At Level 3, i.e., the most miniature metabolic category on the KEGG database, each biomarker OTU was related to 237 pathways (Figure S10), among which significant differences were identified in "transporters" (4.31% in Guiana and 4.91% in South Africa) and "bacterial motility proteins" (1.66% in Guiana and 2.04% in South Africa); however, the greatest difference was only 0.59%, identified in "transporters".In addition, only 17 KEGG Level 3 metabolic pathways were found in the biomarker OTUs from Guiana (red) and South Africa (green), suggesting that metabolic pathways predicted from the two regions were different (Figure 6).At Level 3, i.e., the most miniature metabolic category on the KEGG database, each biomarker OTU was related to 237 pathways (Figure S10), among which significant differences were identified in "transporters" (4.31% in Guiana and 4.91% in South Africa) and "bacterial motility proteins" (1.66% in Guiana and 2.04% in South Africa); however, the greatest difference was only 0.59%, identified in "transporters".In addition, only 17 KEGG Level 3 metabolic pathways were found in the biomarker OTUs from Guiana (red) and South Africa (green), suggesting that metabolic pathways predicted from the two regions were different (Figure 6).

Discussion
Five and three species of epilithic lichens were sampled from the Venezuelan Guiana Shield and the South African Highveld Plateau, respectively.All of these species but one showed no overlaps between the regions, Carolina shield lichen (Canoparmelia caroliniana) was the only lichen species common to the two regions: six specimens, G07 to G12, from the Guiana Shield and two specimens, SA01 and SA02, from South Africa (Table 3).The associated bacterial microbiomes of the monospecific samples showed apparent regional features (Figure 3), implying that the associated microbiomes are controlled more regionally by climate, particularly by rainfall [24][25][26][27], than host lichen species.In contrast, the monospecific (Canoparmelia caroliniana) samples from grassland and bush, SA01 and SA02, respectively, showed no or little intra-regional variation (Figure 3), which may implicate a case of control by host species rather than inter-regional habitat variety in the vicinity.It is suggested that different alphaproteobacterial families of the order Rhizobiales (a synonym of Hyphomicrobiales [60]) are distributed in lichens [61], which is not yet well explained but may contribute to monospecific intra-regional variation.Lichen-Associated Rhizobiales 1 (LAR1) is a lineage of previously uncultured and most frequent non-cyanolichen-associates within the order Rhizobiales [13].Nitrogen fixation was presumed as an eco-physiological function of LAR1 with which to rival cyanobacteria [13]; however, strains of Lichenihabitans psoromatis, the first cultured LAR1 species, possess no relevant genes in their genomes [62].This study detected LAR1 only at <1% read abundance.Instead, the most abundant biomarker for Guiana Shield was affiliated with nitrogen-fixing Beijerinckia mobilis (Figure 5).Interplay among the hosts, green algal/cyanobacterial

Discussion
Five and three species of epilithic lichens were sampled from the Venezuelan Guiana Shield and the South African Highveld Plateau, respectively.All of these species but one showed no overlaps between the regions, Carolina shield lichen (Canoparmelia caroliniana) was the only lichen species common to the two regions: six specimens, G07 to G12, from the Guiana Shield and two specimens, SA01 and SA02, from South Africa (Table 3).The associated bacterial microbiomes of the monospecific samples showed apparent regional features (Figure 3), implying that the associated microbiomes are controlled more regionally by climate, particularly by rainfall [24][25][26][27], than host lichen species.In contrast, the monospecific (Canoparmelia caroliniana) samples from grassland and bush, SA01 and SA02, respectively, showed no or little intra-regional variation (Figure 3), which may implicate a case of control by host species rather than inter-regional habitat variety in the vicinity.It is suggested that different alphaproteobacterial families of the order Rhizobiales (a synonym of Hyphomicrobiales [60]) are distributed in lichens [61], which is not yet well explained but may contribute to monospecific intra-regional variation.Lichen-Associated Rhizobiales 1 (LAR1) is a lineage of previously uncultured and most frequent non-cyanolichen-associates within the order Rhizobiales [13].Nitrogen fixation was presumed as an eco-physiological function of LAR1 with which to rival cyanobacteria [13]; however, strains of Lichenihabitans psoromatis, the first cultured LAR1 species, possess no relevant genes in their genomes [62].This study detected LAR1 only at <1% read abundance.Instead, the most abundant biomarker for Guiana Shield was affiliated with nitrogen-fixing Beijerinckia mobilis (Figure 5).Interplay among the hosts, green algal/cyanobacterial photobionts, LAR1, and other bacterial associates may influence "monospecific intra-regional variation" and would be a focus of the emerging multi-meta-omics [10].
Among the associated bacterial OTUs, cyanobacterial OTUs were present in all of the lichen samples at ratios ranging from 0.11% to 13.26% of the total OTUs (Table S4).High frequencies, >10%, were seen in Usnea florida (G04) and Menegazzia terebrata (G05) but not in samples of the same species, G03 and G06, respectively, indicating not speciesrelated but site-related frequencies of cyanobacterial OTUs.Site-specific inclusion of not lichen-associated but nearby cyanobacteria may not be ruled out.Crocodia aurata (SA06 to SA08), potentially representing a cephalolichen, showed low frequencies, 0.46% to 1.62%, of cyanobacterial OTUs, implying little involvement of cephalodial cyanobacteria in the studied specimens.
The contribution of the cyanobacterial partner to the overall functioning of the lichen can vary depending on the species of the lichen and the specific ecological conditions.Although a technical threshold of 1% was set in this study, the eco-physiologically significant threshold for defining a cyanobacterium as a photobiont in a lichen must be specified.This can vary depending on the particular circumstances and the ecological importance of the cyanobacterium to the lichen community.
The photobiont biomass can range from less than 1% to greater than 90% of the total biomass of the lichen.This variation results from the different types of photobionts and their roles in symbiosis and from environmental factors such as light intensity, temperature, and water availability.Moreover, this biomass distribution can profoundly affect the morphology and ecology of lichens [1,63].A recent comparison of the two sequencing methods revealed that most lichens have a single dominant photobiont genotype, which is representative of the vast majority of the thallus population [64].Cyanobacteria occur in ca.10% of the nearly 20,000 lichen species known as cyanolichens [65].
Cyanolichens are found in various terrestrial habitats, including tropical rainforests, semideserts, and arctic tundra.Their diversity and abundance are highest in humid environments [66].Epiphytic species thrive in the moist and cool conditions of higher elevations in tropical mountains and maritime regions at higher latitudes.They are frequently abundant in the epiphyte communities of old-growth boreal and temperate forests, where they intercept and help retain atmospheric moisture, sequester nutrients, and provide habitat and food for numerous invertebrates [67].Numerous epiphytic species flourish in microhabitats characterized by moderate light intensities, abundant moisture, and periodic drying events [68].
Tripartite cyanolichens host both green algal and cyanobacterial photobionts.In these lichens, the cyanobacteria, which are typically a minor component of the total photobiont biomass, are confined to structures known as cephalodia.Some green algal lichens frequently form ephemeral associations with free-living cyanobacteria, most likely to gain access to a fixed nitrogen source [24,69].
Field evidence implies that the identity of the Nostoc symbiont of bipartite and tripartite lichens depends on the degree of lichenization of the mycobiont than on the collection site [70].In addition, the same Nostoc can be modified by the fungal host to function as the primary photobiont in a bipartite association or as a partner in a tripartite association [71].The lichen-forming fungus can regulate the specific function of the cyanobacterium to optimize its fitness [72,73].
No cyanobacterial OTUs were screened as regional biomarkers.The (non-cyanobacterial) biomarkers were used to predict regional metabolic features such as "Energy metabolism", "Porphyrin and chlorophyll metabolism", and "Photosynthesis" for the Guiana Shield (Figure 6), which may be related to non-cyanobacterial photosynthetic activities, in addition to cyanobacterial photosynthesis, in the humid highland.In contrast, "Bacterial motility proteins" and "Flagellar assembly" were predicted for the South African Plateau (Figure 6), which may be linked to features that would help expand bacterial distribution and adhesion to a substrate in the dry highland.Geographical settings such as climate may thus have more impact on the microbiomes and relevant metabolisms associated with the epilithic lichens than host lichen species.This view, however, should be tested by comparing KEGG predictions derived from more diverse habitat features.
Cellular functions, as well as gene expressions of bacterial associates, would also be viewed as an integral part of a meta-organism, lichen, rather than as the sum of responses of individual species [61].

Figure 1 .
Figure 1.Sampling sites of epilithic lichens at Churi Tepui in the Venezuelan Guiana Shield and the Golden Gate Highlands National Park in the South African Highveld Plateau.The ETOPO1 GlobalRelief Model[30] is the source of the map image.

Figure 1 .
Figure 1.Sampling sites of epilithic lichens at Churi Tepui in the Venezuelan Guiana Shield and the Golden Gate Highlands National Park in the South African Highveld Plateau.The ETOPO1 GlobalRelief Model[30] is the source of the map image.

Figure 2 .
Figure 2. Compositions of the OTU-derived bacterial phyla of lichens from the Venezuelan Guiana Shield (G01 to G12) and the South African Highveld Plateau (SA01 to SA08).Twelve phyla were observed with >1% read abundances.Compositions of the OTU-derived bacterial classes, orders, families, and genera are shown in Figures S2-S5.

Figure 2 .
Figure 2. Compositions of the OTU-derived bacterial phyla of lichens from the Venezuelan Guiana Shield (G01 to G12) and the South African Highveld Plateau (SA01 to SA08).Twelve phyla were observed with >1% read abundances.Compositions of the OTU-derived bacterial classes, orders, families, and genera are shown in Figures S2-S5.

Figure 3 .
Figure 3. PCA plot (left) and hierarchical clustering dendrogram (right) of OTU-derived bacterial species of the lichens from the Venezuelan Guiana Shield (red) and the South African Highveld Plateau (green).PCA plots at higher taxa (genus, family, order, class, and phylum) are shown in Figure S6.The regional distinctions of OTUs are presented in the phylogenetic cladogram (Figure 4).Substantial biomarkers, listed in Table7, were selected by setting the LDA score threshold to 4.5.Substantial biomarkers for the Venezuelan Guiana Shield included two taxa, the family Acetobacteraceae and the order Rhodospirillales, both affiliated with the phylum Pseudomonadota.Biomarkers for the South Africa Highveld Plateau included seven taxa: the genus Sphingomonas, the family Sphingomonadaceae, the order Sphingomonadales, the genus-raked EU289441_g of the family Beijerinckiaceae, and the order Frankiales affiliated with the class-rank Actinomycetota_c of the phylum Actinomycetota.No biomarkers at the species rank were identified, with a threshold value of 4.5.

Figure 3 .
Figure 3. PCA plot (left) and hierarchical clustering dendrogram (right) of OTU-derived bacterial species of the lichens from the Venezuelan Guiana Shield (red) and the South African Highveld Plateau (green).PCA plots at higher taxa (genus, family, order, class, and phylum) are shown in Figure S6.

Figure 4 .
Figure 4. Microbiomic biomarkers of bacteria associated with the lichens from the Venezuelan Guiana Shield (red) and the South African Highveld Plateau (green), as shown in the cladogram generated by LEfSe [40].The concentrically arranged nodes correspond to the domain bacteria, phylum, class, order, family, genus, and species from innermost to outermost.Red and green nodes/shades indicate significantly higher relative abundances of taxa.The diameters of the node circles are proportional to the abundance of corresponding taxa.

Figure 4 .
Figure 4. Microbiomic biomarkers of bacteria associated with the lichens from the Venezuelan Guiana Shield (red) and the South African Highveld Plateau (green), as shown in the cladogram generated by LEfSe [40].The concentrically arranged nodes correspond to the domain bacteria, phylum, class, order, family, genus, and species from innermost to outermost.Red and green nodes/shades indicate significantly higher relative abundances of taxa.The diameters of the node circles are proportional to the abundance of corresponding taxa.

Figure 5 .
Figure 5. Significant differences (p < 0.05) in relative abundances of the most substantial biomarker OTUs from the Venezuelan Guiana Shield (red) and the South African Highveld Plateau (green) analyzed by ANCOM-BC.(Left) the most substantial biomarker for the Venezuelan Guiana Shield, Beijerinckia mobilis affiliated with the phylum Pseudomonadota.(Right) the most substantial biomarker for the South African Highveld Plateau, EU289441_s, affiliated with the genus Lichenibacterium of the same phylum Pseudomonadota.The bottoms and tops of boxes indicate the first and third quartiles, respectively; the bottoms and tops of whiskers indicate the 1.5 interquartile range beyond the lower and upper quartiles, respectively; the circles indicate the original data (including outliers); the crosses indicate the averages; and the horizontal lines indicate the medians.Other substantial biomarker OTUs are shown in Figure S7.

Figure 5 .
Figure 5. Significant differences (p < 0.05) in relative abundances of the most substantial biomarker OTUs from the Venezuelan Guiana Shield (red) and the South African Highveld Plateau (green) analyzed by ANCOM-BC.(Left) the most substantial biomarker for the Venezuelan Guiana Shield, Beijerinckia mobilis affiliated with the phylum Pseudomonadota.(Right) the most substantial biomarker for the South African Highveld Plateau, EU289441_s, affiliated with the genus Lichenibacterium of the same phylum Pseudomonadota.The bottoms and tops of boxes indicate the first and third quartiles, respectively; the bottoms and tops of whiskers indicate the 1.5 interquartile range beyond the lower and upper quartiles, respectively; the circles indicate the original data (including outliers); the crosses indicate the averages; and the horizontal lines indicate the medians.Other substantial biomarker OTUs are shown in Figure S7.

Figure 6 .
Figure 6.Excerpt of KEGG Level 3 pathways based on the biomarker OTUs of microbiomes associated with lichens from the Venezuelan Guiana Shield (red) and the South African Highveld Plateau (green).The horizontal axis indicates relative abundances (%) to be compared between the two regions.The pathways are screened by the cutoff values of >|0.0015| for the two regions' mean relative abundance distances.Significant differences evaluated at p < 0.05 are indicated on the right side.

Figure 6 .
Figure 6.Excerpt of KEGG Level 3 pathways based on the biomarker OTUs of microbiomes associated with lichens from the Venezuelan Guiana Shield (red) and the South African Highveld Plateau (green).The horizontal axis indicates relative abundances (%) to be compared between the two regions.The pathways are screened by the cutoff values of >|0.0015| for the two regions' mean relative abundance distances.Significant differences evaluated at p < 0.05 are indicated on the right side.

Table 2 .
Forward (F) and reverse (R) primers for PCR amplification of the target sequences.
• C for 3 min; 25 cycles of 95 • C for 30 s, 55 • C for 30 s, and 72 • C for 30 s; and finally, 72

Table 4 .
Numbers of MiSeq-generated V3-V4 region reads, 97% similarity-based OTUs, OTU-derived species, genera, families, orders, classes, and phyla in each sample.Due to overlaps among samples, the subtotal and total numbers of taxa are smaller than the simple sums.Mean lengths of valid reads are also listed.Samples G01 to G12 were collected in the Venezuelan Guiana Shield, and samples A01 to A08 were collected in the South African Highveld Plateau.

Table 5 .
Numbers of assigned OTUs and OTU-derived taxa (species, genera, families, orders, classes, and phyla) that were detected only in the Venezuelan Guiana Shield, only in the South African Highveld Plateau, and in both highland regions.The total numbers are the same as those in Table4.

Table 6 .
Alpha diversity indices (Chao1, Shannon, and Simpson) for the bacterial OTUs of 12 epilithic lichen samples from Venezuelan Guiana Shield (G01 to G12) and eight samples from South African Highveld Plateau (SA01 to SA08).The Shannon and Simpson indices calculated values of the effective numbers of species (ENS).

Table 7 .
Biomarker OTUs and the corresponding taxa with the LDA scores > 4.5 identified in the lichen-associated microbiomes from the Venezuelan Guiana Shield and the South African Highveld Plateau.

Table 7 .
Biomarker OTUs and the corresponding taxa with the LDA scores > 4.5 identified in the lichen-associated microbiomes from the Venezuelan Guiana Shield and the South African Highveld Plateau.