Metagenomic analysis of lichen-associated bacterial community profiling in Roccella montagnei

A lichen is a composite organism formed of algae or cyanobacteria that live in a mutually advantageous symbiotic relationship with the filaments (hyphae) of fungus. Three lichen samples were obtained from diverse sites at a terrestrial habitat located in Coimbatore and coastal habitats located in Kanyakumari and Nagapattinam districts of Tamil Nadu. Amplification and sequencing of 16S rRNA V3–V4 regions were used for metagenomic study. Aside from the Next-Generation Sequencing data (NGS), distinct types of lichen microbiome profiles were clearly revealed. The bacterial diversity in the lichen genera of Roccella montagnei growing in coastal and terrestrial environments was further investigated using common and unique operational taxonomic units (OTUs) and the QIIME pipeline (1.9.1). Using similarity clustering, the heat map analysis depicts the abundance information of chosen OTUs as well as the similarity and difference between OTUs and lichen samples. Using multiple methods, the alpha and beta diversity analysis revealed that there were differences in all of the samples. However, UPGMA tree inference of comparable bacterial community in coastal habitat lichen samples compared to terrestrial habitat validates their evolutionary lineage. As a result, the bacterial population associated with corticolous lichen is dependent on geographic locations, growth substrate, and climatic circumstances of similar lichen genera produced in different habitats and tree substrates.


Introduction
Lichens are symbiotic organisms, which were formed by the mycobiont and photobiont partner. Lichens have recently been documented as advanced ecosystems that contain a distinct range of microbial lineages from all around the tree, rock and soil of life (Arnold et al. 2009;Hodkinson and Lutzoni 2009;Cardinale et al. 2012;Bates et al. 2011a;Hodkinson et al. 2012). Lichens have high tolerance to ecological stresses such as surviving in extremes of temperature, low nutrient or water availability and UV light intensities and, therefore, can be found in a wide range of habitats, where other organisms cannot survive (Beckett et al. 2008;Boddy et al. 2010).
The organic (process, biological process) process of a lichen structure begins when it comes into contact with a suitable photobiont. In most situations, the mycobiont dictates the appearance of a lichen structure; however, there are few exceptions where the photobiont controls the environment. However, the photobiont has an impact on lichen Communicated by Erko Stackebrandt. ontogeny, and the lichen structure is not fully grown until each partner contributes (Budel and Scheidegger 2008). The lichen structure is generally a visibly overlaid structure that includes the cortex, affiliated protista layer, medulla, and rhizinae in certain circumstances (Jahns 1973). The dependent partners make contact at intervals in the protista layer, but the relationship between the protista and the fungal hyphae is different. The top cortex is made up of densely packed blood cells that cover the lichen structure's highest point. This lichen component layer can be hundreds of micrometers thick, yet it can also appear as a thin coating (Budel and Scheidegger 2008). Apart from offering protection against high radiation such as ultraviolet (UV) (Solhaug et al. 2003), actinic radiation, and actinic ray light, the upper cortex also has numerous roles such as altering energy budgets and defending against lichenivorous animals (Budel and Scheidegger 2008).
Though the fungus hyphae surround the protoctist cells, they are gently interlinking and leave enough area for optimum intensity. All the same, photobionts are rarely found deep inside the lichen structure, and there are no unattached protoctists within the lichen structure (Jahns 1973). In primarily leafy or fruticose lichens, the medullary layer is an enlarged neighborhood of the inner lichen structure; this layer comprises long and loosely packed fungal hyphae, providing a soft layer with high internal airspace, permitting gas exchange at intervals the lichen structure (Honegger 1993). In some lichens, notably fruticose lichens, supporting tissue is typically designed at intervals of the medulla, making irregular forms that might be discovered in genera like genus Cladonia and Usnea (Budel and Scheidegger 2008).
The bacterial population found in lichens was first reported in the early twentieth century (Uphof 1925;Henkel and Yuzhakova 1936;Iskina 1938). Various bacterial species were identified to be linked with lichens in early investigations based on culture-dependent methods, including classes such as Gammaproteobacteria, Alphaproteobacteria and Firmicutes (Iskina 1938;Panosyan and Nikogosyan 1966;Henkel and Plotnikova 1973). The study of lichen-associated bacterial community was analyzed from three lichen samples such as Cladonia arbuscula, Lecanora polytropa and Umbilicaria cylindrical (Grube et al. 2009).
The 16S rRNA gene includes 1500 base pairs and nine hypervariable regions with different levels of conservation (V1-V9). More conservative zones are useful for recognizing higher-ranking taxa, whereas more rapidly growing zones can help identify genus or species. In the vast majority of cases, the V3/V4 region is chosen because it has the most nucleotide diversity and discriminating capability. It is important to note, however, that no one place is capable of distinguishing between all bacteria. More conservative zones are valuable for identifying higher-ranking taxa, whilst more rapidly developing ones can aid in genus or species identification. Metabarcoding with the 16S rRNA marker is commonly used in research of diverse microbial populations (Babilonia et al. 2018). Metagenomic analysis involves within the application of bioinformatics tools to review the genetic material from environmental, uncultured microorganisms. The next-generation sequencing of 16S rRNA lets the evaluation of bacterial diversity and distinguish thousands of organisms. There are three primary phases in metagenomic data processing: (1) assembly, (2) annotation, and (3) statistical analysis (Thomas et al. 2012).
Lichen is a symbiotic association of fungi and algae occurs and distributed in tree, rock and soil habitat. The symbiotic association of lichen establishment in tree, rock, soil is influenced by several bacteria, actinobacteria and cyanobacteria community. The type, kind of bacterial community in the lichen growing substrate like tree, rock, soil is unique and distinct in geographical locations and substrate chemical composition. Nevertheless, the report on corticolous lichen-associated bacterial community analysis in Roccella montagnei in Tamil Nadu, India is very scanty. Keeping above points, and to accomplish the knowledge gap, therefore, we are concerned to investigate bacterial community profiling of Roccella montagnei in corticolous lichens, which are grown under different and same tree substrate in distinct geographical locations by next-generation sequencing approach of Illumina Sequencing by 16S rRNA amplicon analysis.

Lichen sample preparation
The lichens were thoroughly washed with distilled water, surface disinfected by immersion in absolute ethanol for 5-10 s, followed by immersion in hydrogen peroxide (4%) for 90 s, and then rinsed three times for 60 s each with sterile demineralized water for surface sterilization.

Extraction of lichen metagenomic DNA
Following the modified procedures of Silva et al. (2014) and Hassan et al. (2018), total metagenomic DNA was extracted from the lichen samples. 2 g of the lichen samples was transferred to a sterile 15 mL centrifuge tube after crushing them into a paste form with a mortar and pestle using lysis buffer (3%(w/v) CTAB, 100 mM Tris HCl (pH 8), 2 M NaCl, 25 mM EDTA (pH 8), 4% β-Mercaptoethanol (v/v) and 5% PVP (w/v)). Heat shock method was administrated by keeping the samples at 95 ℃ water bath for 1 h and incubated in -80 ℃ freezer for 1 h. Immediately thawing samples by re-immersing them in a 95 °C water bath. Absolute ethanol was added double the sample volume, vortexed the samples, were centrifuged at 5000 rpm for 10 min. Then, supernatant was discarded from each samples and add 2 ml of lysis buffer, 10 μl lysozyme, 10 μl proteinase K and 10 μl RNase and vortexed for 15 min. The set up was incubated for 1.5 h at 60 °C with intermittent vortexing and a buffersaturated phenol: chloroform: isoamyl alcohol (1:1:1) was added to the mixture at an equal volume. The aqueous layer was transferred to a sterile centrifuge tube after 20 min of centrifugation at 8000 rpm at 4 °C, and an equal amount of chloroform: isoamyl alcohol (1:1) was added to the aqueous layer for the elimination of phenol traces and repeated twice. The supernatant was transferred to a sterile centrifuge tube with an equivalent amount of isopropanol and incubated on ice for 20 min after centrifugation. The supernatant was removed after centrifugation, and the pellet was rinsed with 100 μL of 80% ethanol and centrifuged for 5 min in 4 °C at 8000 rpm. The supernatant was removed, and the DNA pellets were dried completely at room temperature. The pellet was resuspended in 50 μL of TE buffer and kept at − 20 °C until needed. A total of nine lichen metagenomics DNA (three samples from each location) was extracted from Roccella montagnei and mixed together and used for analysis.

16S rRNA V3/V4 gene amplification by NGS
The following procedures were followed for the creation of next-generation sequencing libraries and Illumina MiSeq sequencing. A Qubit 2.0 Fluorometer was used to measure DNA samples (Invitrogen, Carlsbad, CA, USA). A MetaVxTM Library Preparation kit was used to produce amplicons from 30 to 50 ng of DNA (GENEWIZ, Inc., South Plainfield, NJ, USA). The hypervariable V3 and V4 regions of bacterial 16S rRNA were chosen for amplicon sequencing, followed by taxonomic analysis. GENEWIZ created a set of proprietary primers aiming at the relatively conserved areas of the V3 and V4 hypervariable regions of the 16S rRNA of bacteria.

Sequencing and library preparation
The Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) was used to confirm library quality control, and the Qubit 2.0 Fluorometer was used to quantify it. Following the manufacturer's recommendations, DNA libraries were multiplexed and put into an Illumina MiSeq device (Illumina, San Diego, CA, USA). The MiSeq control software package (MCS) integrated in the MiSeq instrument was used to sequence using a 2300 paired-end (PE) configuration; image analysis and base calling were accomplished using the MiSeq control software package (MCS). The metagenome sequence of lichen-associated bacterial community in Roccella montagnei was submitted in NCBI with the accession number of bioproject PRJNA741221.

Sampling and process of lichen Roccella montagnei
The Roccella montagnei was collected from two different habitats of sampling locations like coastal and terrestrial regions. In the present study, we have examined the bacterial community profiling of the Roccella montagnei collected from the different habitat, corticolous lichen grown under different and same tree substrate in distinct geographical locations (Table S3). A total of three Roccella montagnei lichen samples were used for metagenomic DNA extraction. The lichen samples were abbreviated as PR lichen from terrestrial habitat and RR, KR lichen from coastal habitat.
Lichen microbiome profiling QIIME pipeline (1.9.1) was used to analyze the bacterial diversity in the lichen genera of Roccella montagnei obtained from the coastal and terrestrial habitats. In the 16S rRNA dataset, 2.1 lakhs sequences number were obtained having a combined length of 390-470 base pairs (bp) ( Table S4).

Common and unique operational taxonomic units (OTUs)
16S rRNA V3-V4 gene profiles were generated from the lichen samples to compare and evaluate the diversity and abundance of bacterial association in the lichen Roccella montagnei. The circles with different colors in the Venn diagram ( Fig. 1), each represents a different sample of corticolous lichen grown under same tree substrate in distinct geographical locations and numbers within them represents the number OTUs which are either common or unique to the samples, the number inside the overlapping part denotes the common OTUs, whereas in the non-overlapping ones are unique to the samples.
A sum total of 424 OTUs were observed from three samples, out of which 205 OTUs were in common with all three lichen samples, 35 OTUs are unique to KR lichen sample from coastal habitat, 34 for PR sample from terrestrial habitat and about 42 for RR lichen from coastal habitat. It is also observed that the highest OTUs were shared between KR and RR lichen sample which are both from coastal habitats; on the other hand, PR and RR lichen share the lowest of 20 OTUs.
The heatmap analysis illustrates the abundance of bacterial community information of selected OTU as well as the similarity and difference across OTUs and lichen samples by similarity clustering. The figure showed the top 30 OTUs with the highest abundance of bacteriome.

Taxonomic analysis of bacterial community abundance on Roccella montagnei
The bacterial community abundance was analyzed in Roccella montagnei for all taxonomic levels and the results were observed as follows, their respective stacked bar plots. The most abundant phyla in the 3-lichen Roccella montagnei samples, which were in coastal and terrestrial habitat were (i) Proteobacteria with 95.7% present in Nagercoil lichen sample, 96.5% in Pollachi lichen sample, and a 96.1% in lichen sample collected from Vedaranyam; (ii) Bacteroidetes having a presence of 1.3% in Nagercoil lichen sample, 1.87% in Pollachi lichen sample and 1.7% in Vedaranyam sample (  (Fig. 2).

Alpha diversity analysis of the lichen microbiome communities
Alpha diversity is mainly used to reflect the diversity of each lichen sample, which estimates the number of species in the microbial community as well as the abundance and diversity of species in environmental communities through a series of statistical indices. In which are summarized in terrestrial and coastal habitat lichens in distinct geographic locations of tree as growing substrate of Cocos nucifera and Maba buxifolia in Nagercoil, Vedaranyam, Pollachi of Tamil Nadu, India (Table S7). The indices for community richness calculation include: Chao and ACE index; the Shannon and Simpson index for community diversity calculations; the Goods coverage indices to measure sequencing depth (Coverage). Nagercoil lichen had the highest Shannon index (4.7) which indicates abundance of the lichen-associated bacterial community present in a significantly maximum in the trees of Cocos nucifera in coastal habitat in comparison with other two lichen samples, i.e., Pollachi sample (4.1) and Vedaranyam sample (4.5).
The Simpson's index shows that quantify the diversity of species richness and evenness is almost the same in Nagercoil lichen from Cocos nucifera in terrestrial locations (0.9) and Vedaranyam lichen from Maba buxifolia in Coastal locations (0.9) and both the samples were higher than the Pollachi Cocos nucifera in terrestrial location (0.8).
The Fig S8 shows the rarefaction curve, i.e., it tries to determine all of the bacterial diversity in a sample at a particular sequencing depth by sub-sampling the sequence to a minimum and plotting it against the observed OTUs which enables the comparison among the different lichen samples in of different and same tree substrate in distinct geographical locations. The curves represent the bacterial diversity as a function of sequencing depth. The X axis is the number of valid sequences extracted, and the Y axis is the number of OTUs (Observed OTUs). Each sample is represented by one curve with a unique color. The number of OTUs increases with the increase of extracted sequence count until reaching a plateau, which indicates the number of detected OTUs will not increase with the number of extracted sequences and reflects the reasonable sequence depth.
Rank-abundance curve is used to analyze bacterial diversity. Rank-abundance curve reflects both bacterial species abundance and species uniformity in lichen. The abundance of species is reflected by the length of the curve on the X axis. The more extended on the X axis, the more abundant the species is. Species uniformity is reflected by the shape of the curve. The smoother the curves, the higher the species uniformity in lichen Roccella montagnei.
Each curve in the figure above corresponds to an individual lichen Roccella montagnei samples from three locations. The X axis is the relative abundance of the OTU in descending order. The Y axis is the relative abundance of the OTU. '100' on the X axis indicates the OTU in the sample is ranked as the 100th abundant in descending order, and the corresponding value on the Y axis is the percentage of the sequence count in the OTU (the number of sequences of the OTU divided by the total number of sequences) (Fig. S9).

Beta diversity of the lichen microbiome communities
Beta diversity analysis is used to determine how different the lichen bacterial community in coastal samples is from that in terrestrial habitats, as well as to disclose elements of microbial ecology that are not apparent simply looking at the sample composition. Beta diversity analysis is frequently divided into two categories. They can be quantitative (e.g., Bray-Curtis or weighted UniFrac) or qualitative (e.g., Bray-Curtis or weighted UniFrac) (Goodrich et al. 2014). The NMDS, PCA, and PCoA figures were all created using beta diversity distance matrices.

PCA and PCoA plot analysis
PCA analysis (Principal Component Analysis): the difference and distance between samples can be reflected by the analysis of the gene functional distribution of different samples of Roccella montagnei. The first, second, and third main components are represented by PC1, PC2, and PC3, respectively (Fig. 3a). The percentage after the main component, which reflects how much information the primary component can extract from the original data, represents the contribution rate of this component to sample difference. The distance between coastal and terrestrial lichen demonstrates the similarity of the distribution of functional categories in the sample.
The PCoA (Principal Co-ordinates Analysis) is a method for analyzing and visualizing data differences and similarities. Illustrations from the same group are identified by the similar color and form. The PCoA plot for the first and second principal coordinates is shown in the PC1 vs PC2 image (Fig. 3b), with the first and second axes representing the first and second principal coordinates, respectively. The contribution of the related coordinate to the sample variations is shown by the percentage number in the axis label, which illustrates how much information is recovered from the original data. The distance between the sample points demonstrates the similarity of the bacterial population in the lichen.
According to the Principle Component Analysis (PCA) and Principal Co-ordinates Analysis (PCoA), the lichen samples from Nagercoil, Pollachi, and Vedaranyam were differed significantly. Mostly, the Nagercoil lichen from Cocos nucifera in coastal habitat and Vedaranyam lichen from Maba buxifolia in coastal locations, indicated similarity than the other sample comparisons.

UPGMA tree
Cluster analysis uses evolutionary data from sample sequencing to determine if samples in a well-known environment are significantly different from an evolutionary branch in bacterial communities. In UPGMA (Unweighted pair group method with arithmetic mean) tree, each branch represents a sample and difference of colors will represent different groups present. The results (Fig. 4) infer the all the 3-lichen bacterial community belongs to the same group and the evolutionary microbial lineage between Nagercoil lichen from Cocos nucifera from coastal habitat and Vedaranyam lichen from Maba buxifolia from Coastal habitat are similar and the highly different from the Pollachi lichen from Cocos nucifera from terrestrial habitat.

Discussion
The composition of bacterial microbiota has been well studied in lichens, nevertheless there are only limited reports in respect to comparing the bacteriome of the same lichen species from different types of locations in Tamil Nadu, India and Roccella montagnei is a species whose bacteriome had not been explored previously. In this study, we have collected Roccella montagnei from three different locations that were near coastal and terrestrial habitats and metagenomic analysis and their bacteriobiome records through NGS-Illumina sequencing platform. The predominance of Proteobacteria, Bacteroidetes, Actinobacteria, Acidobacteria and Verrucomicrobia in phylum level revealed and consistent with previous reports (Hodkinson et al. 2012;Leiva et al. 2021;Printzen et al. 2012;Sierra et al. 2020). In the lichen genera such as Cora, Hypotrachyna, Cladonia, Usnea, Sticta, Stereocaulon and Peltigera, the presence of these members of this phylum is considered to play important roles such as symbioses by supplying nutrients, iron and phosphate mobilization, and fixing nitrogen (Grube et al. 2015;Sigurbjornsdottir and Vilhelmsson 2016).
Due to chemical and physiological variations, lichenassociated bacteria colonize different thallus portions in varying abundances and patterns. Lobaria pulmonaria's upper and bottom surfaces are covered in alphaproteobacteria (Cardinale et al. 2012). Other dorsoventrally arranged lichen species, such as Umbilicaria sp., have the same pattern (Grube et al. 2015). Cardinale et al. 2008 andBates et al. (2011b) both found bacterial species that were present in all lichen samples. Actinobacteria, Firmicutes, and Gammaproteobacteria were found in four different foliose green algal lichen species (Parmelia sulcata, Rhizoplaca chrysoleuca, Umbilicaria americana, and Umbilicaria phaea) collected from granite rock outcrops in northern Colorado (40.01° N,105.47° W) in previous studies of lichen-associated bacteria. The dominance of Gammaproteobacteria over Alphaproteobacteria at class level was found in our results and this finding is supported by similar discovery reported by Sierra et al. (2020), while most of the reports majorly showed the dominance of alphaproteobacterial which is ubiquitous in all the lichen of different habitats (Almendras et al. 2018;Aschenbrenner et al. 2014;Cernava et al. 2015b;Hodkinson et al. 2012;Leiva et al. 2021;Printzen et al. 2012).
In the present study, microbiome analysis in lichen reveals that, among phylum level, Proteobacteria and Bacteroidetes are highly dominant and similar in all three different locations including the coastal and terrestrial lichen samples. Nevertheless, differences in all other phylum such as Acidobacteria, Actinobacteria, Firmicutes and Cyanobacteria. Among the 13 phylums, only two phylums were similar, whereas microbial community differs in 11 phyla. We observed high bacterial community differences in class, order, family and genus level in Roccella montagnei. The predominance of Pseudomonadales of Gammaproteobacteria indicates a different result from the previous studies that revealed typically the dominance of Rhizobiales, which is only the second most dominant order in our study, other orders such as Rhodospirillales, Sphigomonadales were observed in reports of Almendras et al. (2018) Sierra et al. (2020). Xanthomonadaceae and others like Sphingomonadaceae and Acetobacteraceae that are mostly similar to the reports. The genus level bacteriome exhibits the major dominance of unclassified genera; on the other hand, the classified genera were dominated by Pseudomonas Cernava et al. 2015a, b;Lee et al. 2014), Ochrobactrum which is surprisingly never before found in the lichen bacteriome and is followed by Stenotrophomonas in lung lichen Loboria pulmonaria were sampled from three different locations in Austria (Cernava et al. 2015a, b). The occurrence and presence of Stenotrophomonas indicate that the lichen can survive drought and high levels of salinity by producing protecting compounds such as spermidine and different osmolytes, which is relatable in our study as two of the samples where from coastal habitat, Sphingomonas Graham et al. 2018;Lee et al. 2014) that helps the lichen grow by producing growth-promoting hormones, Acinetobacter (Park et al. 2016) and others. The alpha and beta diversity analysis is done to identify the similarity of the microbial population of the lichen form the different habitats and all the analysis has led to the conclusion of similarity of coastal lichens and then distinctively different from the terrestrial lichen as proven by the UPGMA tree at an evolutionary lineage level in the microbial communities.

Conclusions
The metagenomic analysis of lichen Roccella montagnei data has been generated; the findings illustrate that a total of 424 OTUs were obtained from three lichens. Among the total OTU Vedaranyam lichen coastal location sample showed the greatest number of unique OTUs (42) and lowest was Pollachi lichen from Cocos nucifera in terrestrial habitat. However, the majority OTUs (57) were shared between Nagercoil lichen from Cocos nucifera (Coastal); Vedaranyam lichen from Maba buxifolia in Coastal locations than by Pollachi sample the others. The taxonomic analysis of microbial community abundance on Roccella montagnei showed the abundance of Proteobacteria, Bacteroidetes and Actinobacterias in Phylum level; Gammaproteobacteria, Alphaproteobacteria and Betaproteobacteria in Class level; Pseudomonadales, Rhizobiales and Xanthomonadales in Family level and Unclassified genera, Pseudomonas, Ochrobactrum and Stenotrophomonas in the Genus level. All the alpha and beta diversity analyses using various tools have showed variations present in all samples; nevertheless, it has led to the conclusion of similarity of Nagercoil and Vedharanyam coastal lichen samples, the UPGMA tree confirms their evolutionary lineage. The outcome reveals that the bacterial community variation in lichen is due to the geographic origin, the biotic and abiotic factors surrounding the environment, the substrate in which they grow, the microhabitat conditions around the lichen sample. Since, chemical and physiological variations, lichen-associated bacteria colonize different thallus portions in varying abundances and patterns. The data were analyzed in terms of the lichenassociated bacterial community's stability under shifting abiotic settings. They add to our knowledge of a fascinating lifeform by demonstrating how the bacterial component of the holobiont maintains its integrity in the face of adverse environmental conditions. In conclusion, this study finding will be a contribution of bacterial community profile that were uncultivable and their diversity, similarity and dissimilarity have been revealed to this society.