Diversity of bacteria associated with lichens in Mt. Yunmeng in Beijing, China

Study area

Lichen samples were collected at Mt. Yunmeng (40°26′N–40°38′N, 116°30′E–116°50′E), which is located in the northwest of Miyun District, Beijing. The climate of this area is characterized by warm temperate semi-humid continental monsoon. The average annual temperature is 10.9 °C, the average annual precipitation ca. 700 mm, with 76% of the precipitation occurring from June to August. The vegetation was severely destroyed by human activities in the past, but the government classified Mt. Yunmeng and its adjacent areas as a nature reserve and a national forest park in the early 1990s. Since then, a secondary temperate deciduous broadleaf forest developed and became dominant, with planted coniferous forest interspersing some areas. Today about 90% of the Mt. Yunmeng area is covered by forest (http://www.yunmengshan.org.cn/ymsjqjs).

While most parts of Mt. Yunmeng were classified as nature reserve, its eastern part has developed into a series of scenic spots which attract large numbers of tourists. One of these spots is the Black Dragon Pool (160 to 380 m a.s.l), which is a 4.5 km long canyon, with a small creek and a riverine forest on both sides. Tourists’ activities may have destroyed or negatively affected the local lichen communities.

Collection of materials

Lichen samples were collected in October starting from the Black Dragon Pool, crossing the border between the scenic spot and the nature reserve, and ending at the heart of nature reserve at about 1,130 m a.s.l. Lichen samples were identified referring to Brodo, Sharnoff & Sharnoff (2001) and Nash et al. (2007). In total, twenty samples were obtained with three from the scenic spot (F1–F3) and seventeen from the nature reserve (F4–F10, C1–C10), most of which were saxicolous (growing on granites), while only two were epiphytic (growing on Quercus sp.). The number of crustose (C1–C10) and foliose (F1–F10) lichens were equivalent, while no fruticose lichens were found at all (for detailed information see in Table 1). The differences in the number of lichen samples were mainly due to the limitation of the number of locations available for collection in the field. Lichen samples (one entire thalli of each individual) were collected using sterile blade and forceps and placed in sterile polythene bags, sealed and labelled. Samples were transported to the laboratory using an ice box and then stored at −80 °C for further processing. Lichen samples were washed with ultra-pure laboratory grade water to remove dirt and debris. To avoid any possible bias depending on the extreme sensitivity of lichens to environmental changes, all sampling procedures were conducted in one day and under similar microenvironmental conditions.

DNA extraction, 16S rRNA gene amplification and sequencing

HiPure Soil DNA Kit was used to extract total genome DNA according to manufacturer’s protocols (Magen Biotechnology Co., Ltd, Guangzhou, China). DNA concentration was monitored by Qubit® dsDNA HS Assay Kit. A total of 20–30 ng of DNA was used to generate amplicons containing V3–V4 hypervariable regions using a MetaVX Library Preparation kit (GENEWIZ, Inc., South Plainfield, NJ, USA). The forward primer contained the sequence 5′–7CCTACGGRRBGCASCAGKVRVGAAT–3′ and the reverse primer, 5′–GGACTACNVGGGTWTCTAATCC–3′, which were designed by GENEWIZ (Suzhou, China). The polymerase chain reaction (PCR) system contained 2.5 μl of TransStart buffer, 2 μl of dNTPs, 1 μl of each primer, 0.5 μl of TransStart Taq DNA polymerase, 20 ng template DNA, and ddH2O which was added to obtain the total volume of 25 μl. The PCR ran as follows: 3 min of denaturation at 94 °C, 24 cycles of 5s at 95 °C, 90 s of annealing at 57 °C, 10 s of elongation at 72 °C, and a final extension at 72 °C for 5 min. Indexed adapters were added to the ends of the amplicons by limited cycle PCR. Finally, the library was purified with magnetic beads. Paired-end (PE) sequencing was carried out with the Illumina Miseq Platform (Illumina, San Diego, USA) at Genewiz, Inc.

Bioinformatic analysis

PRINSEQ (0.20.4) was used for sequence quality control. After quality filter (reads with length <200 bp) and removal of chimeric sequences, all sequences were grouped into operational taxonomic units (OTUs) using VSEARCH (1.9.6) (sequence similarity was set to 97%) against the SILVA v138 database. The OTU table was filtered to remove chloroplast sequences, and then OTUs were classified into different taxonomic levels using the Ribosomal Database Project (RDP) classifier (2.13). The names of annotated bacteria follow the List of Prokaryotic names with Standing in Nomenclature (LPSN) (Parte et al., 2020). Alpha diversity indices (Chao 1 and Shannon), which indicate bacterial community richness and diversity, were calculated using Mothur software (1.35.1). Principal Co-ordinates Analysis (PCoA) and Unweighted Pair Group Method with Arithmetic Mean (UPGMA) based on the Bray-Curtis distance matrices were employed to determine beta diversity and cluster samples, respectively. Analysis of Similarities (ANOSIM) was applied to test the differences between bacterial communities associated with foliose and crustose lichens. Linear discriminant analysis Effect Size (LEfSe) was conducted using the Galaxy framework online tool (https://huttenhower.sph.harvard.edu/galaxy/), which can find microorganisms with significant differences between groups by calculating the contribution of each microbial abundance to the overall difference. Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt2) was conducted to predict the function of lichen-associated microbiota based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) and the Clusters of Orthologous Groups (COG) databases. All statistical analyses were performed with SPSS 22.0 (IBM Corp., Armonk, NY). Differences in phylum, order and genus relative abundances are presented as Means ± SE. Alpha diversity indices were calculated by using the Independent–sample t–test. A P-value < 0.05 was considered statistically significant, while a P-value < 0.01 indicated that differences were highly significant. The raw data obtained in this study have been submitted to the NCBI Sequence Read Archive under Bioproject number PRJNA972522.

16S rRNA gene sequencing data statistics

After filtration of low-quality regions and removal of chimeras, a total of 1,232,513 effective sequences were obtained, with 36,494–79,179 for each sample and an average length of 442.04 bp (Table 2). Using cluster analysis, we identified 426 OTUs at 97% similarity level, including 16 phyla, 30 classes, 61 orders, 83 families and 125 genera.

The rarefaction curve (Fig. 1A) shows that the number of OTUs in the sample increases with sequencing depth, the curve finally flattens, indicating that the sequencing had covered most bacteria in the samples. Overall, the observed OTUs which stand for species richness in foliose lichens was higher than that of crustose lichens. The rank abundance curve (Fig. 1B) shows two aspects of diversity, namely richness and evenness. In general, the evenness and richness of the bacteria associated with foliose lichens was higher than that with crustose lichens.

Composition of bacteria associated with foliose and crustose lichens

Under the conditions of this study, at the phylum level, the most dominant phylum of bacteria associated with foliose and crustose lichens was Pseudomonadota, with a relative abundance of 55.88 ± 4.45% and 69.76 ± 5.48%, in foliose and crustose, respectively (Fig. 2). The phyla Actinomycetota (foliose: 13.44 ± 2.61%; crustose: 9.46 ± 2.70%) and Acidobacteriota (foliose: 13.63 ± 1.83%; crustose: 12.15 ± 2.72%) were also present with a relatively high abundance. At the order level (Fig. 3), the dominant orders of bacteria were Acetobacterales and Hyphomicrobiales, with relative abundance of 4.44 ± 1.02%, 6.13 ± 0.71% for foliose lichens and 6.62 ± 2.11%, 4.39 ± 1.08% for crustose lichens. At the genus level (Fig. 4), Sphingomonas (foliose: 5.44 ± 1.09%; crustose: 3.14 ± 1.16%) and Acidiphilium (foliose: 2.5 ± 0.54%; crustose:3.99 ± 1.14%) had higher relative abundance in both groups of lichens.

Diversity of bacterial communities associated with foliose and crustose lichens

Under the conditions of this study, Alpha diversity results showed that the Chao1 indices of foliose and crustose lichens (Fig. 5a) were 246.57 ± 20.42 and 151.40 ± 30.38, respectively (P = 0.02). The Shannon indices of bacterial communities associated with foliose and crustose lichens (Fig. 5b) were 3.69 ± 0.55 and 2.71 ± 0.56 respectively but did not show significant differences between the two groups (P = 0.43).

The UPGMA tree (Fig. 6) and the PCoA (Fig. 7) unraveled the microbial composition between different samples. Performing ANOSIM analysis, we found a significant difference between the community structures of bacteria associated with foliose and crustose lichens (P = 0.02).

Taxon differences of lichens-associated bacteria between foliose and crustose lichens

According to phylogenetic map (Fig. 8), there are nine bacterial taxa showed significant differences between foliose and crustose lichens in this study. Among the bacteria associated with foliose lichens, the relative abundance of Chloroflexi, Bacteroidia, Bacteroidota, Micromonosporae, Micromonosporales, Chitinophagae, Chitinophagales, Blastocellaceae, Blastocellales was significantly higher than that of bacteria associated with crustose lichens. No taxon was found to have a significantly higher relative abundance with crustose lichens than foliose lichens.

PICRUSt2 gene function predictions

Functional annotations were based on KEGG and COG databases. Most of the bacterial gene functions were related to metabolism (Fig. 9A), such as the metabolisms of carbohydrates, amino acids, vitamins, and cofactors. The abundance of genes related to amino acid transport and metabolism, function unknown, etc. were high (Fig. 9B). No significant difference in the functional abundance between associated bacteria of foliose and crustose lichens was found under the conditions of this study.