Front. Microbiol. Frontiers in Microbiology Front. Microbiol. 1664-302X Frontiers Media S.A. 10.3389/fmicb.2022.1009588 Microbiology Original Research Giant pandas’ staple food bamboo phyllosphere fungal community and its influencing factors KangLiwen12 LuoWei3 DaiQinglong3 ZhouHong12 WeiWei12 TangJunfeng12 HanHan12 YuanYuan12 LongJuejie12 ZhangZejun12* HongMingsheng12* 1Liziping Giant Panda’s Ecology and Conservation Observation and Research Station of Sichuan Province (Science and Technology Department of Sichuan Province), China West Normal University, Nanchong, China 2Key Laboratory of Southwest China Wildlife Resources Conservation (Ministry of Education), China West Normal University, Nanchong, China 3Liziping National Nature Reserve Administration, Ya’an, China

Edited by: Lifeng Zhu, Nanjing Normal University, China

Reviewed by: Wei Zhang, Guizhou University, China; Jun Zhao, Nanjing Normal University, China; Dunwu Qi, Chengdu Research Base of Giant Panda Breeding, China

 *Correspondence: Zejun Zhang, zhangzejun66@163.com Mingsheng Hong, mingshenghong@cqu.edu.cn

This article was submitted to Microbial Symbioses, a section of the journal Frontiers in Microbiology

 30 09 2022 2022 13 1009588 02 08 2022 06 09 2022 Copyright © 2022 Kang, Luo, Dai, Zhou, Wei, Tang, Han, Yuan, Long, Zhang and Hong. 2022 Kang, Luo, Dai, Zhou, Wei, Tang, Han, Yuan, Long, Zhang and Hong

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Giant pandas have developed a series of foraging strategies to adapt to their special bamboo diets. Although bamboo is an important food resource for giant pandas in Liziping National Nature Reserve (Liziping NR), China, there are relatively few studies on their phyllosphere fungal community and its influencing factors. Herein, we used ITS1 amplification and metagenomic sequencing to analyze the phyllosphere fungi diversity and functions (KEGG, CAZyme, and antibiotic resistance gene) and explore the influencing factors for the three giant pandas foraging bamboo species (Arundinaria spanostachya, AS; Yushania lineolate, YL; and Fargesia ferax, FF) over different seasons (spring vs. autumn) in Liziping NR, China. We found that Ascomycota and Basidiomycota were the most dominant phyla in the bamboo phyllosphere. The alpha diversity (e.g., the Sobs index and Shannon index) was relatively higher in autumn samples than in spring samples, and the community structure differed significantly between the three bamboo species in spring and autumn. Some biotic and abiotic variables (e.g., the elevation and mean base diameter of bamboo) significantly influenced the abundance, diversity, and community structure of the bamboo phyllosphere fungal community. Moreover, the functional analysis showed the differences in the glycoside hydrolase community and antibiotic resistance gene (ARG) profile between spring and autumn samples. Co–occurrence network modeling suggested that AS phyllosphere fungal communities in autumn employed a much more complex network than that in spring, and the abundance of multidrug, tetracycline, and glycopeptide resistance genes was high and closely correlated with other ARGs. These results indicate that fungal community’s abundance, diversity, and community structure are mainly affected by the season, host species, and elevation. The season and host species are major factors affecting the biological functions (KEGG and CAZyme), ARGs, and interactions between sympatric bacterial and fungal communities in bamboo phyllosphere. This integrated study can provide a reference basis for the seasonal management of bamboo resources foraged by wild giant pandas, and predict the risk of antibiotic resistance in bamboo phyllosphere fungal flora in Liziping NR (Xiaoxiangling mountains), China.

 giant panda staple food bamboo phyllosphere fungi metagenome ARGs Liziping NR 31900337 2019QZKK05010502 KCXTD2022-7 National Natural Science Foundation of China10.13039/501100001809 The Second Tibetan Plateau Scientific Expedition Innovation Team Funds of China West Normal University Introduction

The giant panda (Ailuropoda melanoleuca) belongs to the order carnivore and has a typical carnivore gastrointestinal tract. Fossil evidence has shown that the giant panda at least 7 Myr ago started to forage bamboo (Jin et al., 2007; Wang X. M. et al., 2022). As the umami taste receptor gene Tas1r1 had lost its function and other reasons by about 4.2 Myr ago (Zhao et al., 2010), it is clear that giant pandas had already completed their dietary switch by about 2.0–2.4 Myr ago (Jin et al., 2007).

To adapt to the special bamboo diet, giant pandas have developed a series of foraging strategies, including the evolution of pseudo-thumbs, seasonal foraging, feeding point selection, habitat selection, etc. (Endo et al., 1999; Wei et al., 1999; Hong et al., 2015, 2016; Hu et al., 2017). However, in different seasons and mountain regions, giant pandas feed on different parts of bamboo and different kinds of bamboo; however, bamboo leaves are still their predominant food throughout the year (Wei et al., 2015). For example, giant pandas forage bamboo leaves most of the year except May–July in the Qinling Mountains. (Wu et al., 2017). Despite the fact that giant pandas only feed on bamboo, they have not evolved any genome-encoding enzymes specific for cellulose digestion (Hu et al., 2017). Thus, bamboo leaves are an important reason driving the composition and changes of intestinal flora of giant pandas; imbalances in this intestinal flora can lead to gastrointestinal disease, which is the most common cause of death in captive and wild giant pandas (Tun et al., 2014).

As an important food source for giant pandas, bamboo leaves have different nutrient and microbial compositions as compared with bamboo shoots, bamboo stems, and branches (Wei et al., 2015; Wu et al., 2017; Long et al., 2021). The health status of bamboo leaves, including the differences in seasonal nutrition and microorganisms in bamboo phyllosphere, could be the reason for the variations in giant pandas’ seasonal activities and physiology. Using the traditional pathology method, Xu et al. (2014) found that the microbial community composition of phyllosphere differed among different bamboo species, and microbial groups showed seasonal differences. There were also differences in species composition and frequency of endophytic bacteria and fungi in bamboo leaves (Helander et al., 2013), while endophytic fungi in Phyllostachys heteroclada differ in the leaf and branch (Zhou et al., 2017). Using high-throughput amplicon sequencing, Jin et al. (2020) found that the richness and diversity of bacteria and fungi differed significantly between bamboo species.

The phyllosphere is one of the most diverse ecosystems on earth (Lindow and Leveau, 2002; Corinne et al., 2016), inhabited by bacteria, archaea, fungi, and protists (Arnold et al., 2003; Lindow and Brandl, 2003). Bacteria are the most abundant colonizers in leaves. However, phyllosphere fungal colonizers are incredibly diverse (Lindow and Brandl, 2003; Finkel et al., 2011; U’Ren et al., 2012). During the growing season, bacteria dominate, yeasts followed, and filamentous fungi accounted for the least (John and Robin, 2000; Redford and Fierer, 2009). The phyllosphere community composition is significantly impacted by biological and abiotic factors, such as the host species, season, temperature, humidity, host attributes, etc. (John and Robin, 2000; Lindow and Brandl, 2003; Paulino-Lima et al., 2013; Kembel et al., 2014; Copeland et al., 2015).

As an important part of the ecosystem, phyllosphere microorganisms affect plant surface properties (Lindow and Brandl, 2003), which have a key role in the degradation of organic matter, nutrient conversion, and energy flow (Rastogi et al., 2013; Zhao et al., 2013; Horton et al., 2014). Phyllosphere fungi, which can prevent herbivores from eating plants, protect plants from all kinds of pathogens and also make the plant tolerant against abiotic stress (Rehman et al., 2021). The ability to degrade cellulose is widespread in fungi, and is particularly well represented in Ascomycota and Basidiomycota (Boer et al., 2005). In a paddy environment, Fan et al. (2019) found an indigenous fungus, Aspergillus cvjetkovicii, which exerted an effect of strong antagonism on Magnaporthe oryzae and a promoting effect on rice. Zahn and Amend (2019) found that phyllosphere fungi alter the timing and allocation of reproduction in Arabidopsis, resulting in later flowering and higher seed quality. However, with the abuse of antibiotics and the rapid spread and development of antibiotic resistance, antibiotic resistance has become one of the world’s greatest public health dangers (Chen et al., 2019; Pei et al., 2019). Unfortunately, phyllosphere is increasingly known as reservoirs of antibiotic resistance genes (Sun et al., 2021).

Most microorganisms do not exist in isolation, but live in the same space and share common substrates; these common substrates lead to either synergistic or antagonistic interactions, which directly or indirectly affect the growth of host plants and the absorption and utilization of nutrients (Mille-Lindblom et al., 2006; Schmidt et al., 2015; Ma et al., 2020). The interaction of bacteria and fungi can change soil structure and fertility, improve plant adaptability, and enhance the absorption and utilization of plants for nutrients and water (Porras-Soriano et al., 2009; Tang et al., 2020). Bacteria and fungi also work together to promote the rumen digestion of grass (Lee et al., 2000).

As mentioned above, only three studies have been conducted on the phyllosphere microflora of giant pandas’ staple food, bamboo (Helander et al., 2013; Jin et al., 2020; Long et al., 2021). Even fewer have been conducted on the phyllosphere microorganisms of bamboo leaves. Our previous research has focused on environmental factors influencing phyllosphere bacterial communities in bamboo in Liziping National Nature Reserve (Liziping NR), Sichuan, China (Long et al., 2021). According to the results of the Fourth National Giant Panda Survey, there are 22 wild giant pandas in the reserve (Sichuan Forestry Bureau, 2015). As of 2019, 1 rescued giant panda and 9 captive pandas were released into the reserve, alongside 8 surviving giant pandas (Hong et al., 2019). In this study, therefore, we focus on the phyllosphere fungal communities in giant pandas’ staple food bamboo and their influencing factors using next-generation sequencing (NGS) technology (ITS1 region amplicon sequencing and whole-genome shotgun sequencing). This study aimed to (i) compare interspecies differences and seasonal variations in the composition, diversity, and community structure of phyllosphere fungal communities in three bamboo species (Arundinaria spanostachya, AS; Yushania lineolate, YL; and Fargesia ferax, FF); (ii) explore the ecological factors that influence the phyllosphere fungal community changes; (iii) investigate the interactions between fungi and bacteria in bamboo phyllospheres using co-occurrence network analysis; (iv) reveal the differences in biological functions (KEGG and carbohydrate-active enzyme functions, CAZyme) of AS phyllosphere fungi between spring and autumn; and (v) obtain antibiotic resistance gene (ARG) abundances of the AS phyllosphere fungal community and their differences between spring and autumn, determining their interactive relationship. Our research results could provide a reference basis for the restoration, protection, and management of bamboo resources in wild giant panda habitats.

 Materials and methods Study area

Liziping NR (28°51′02″–29°08′42″N, 102°10′33″–102°29′07″E) is situated in Shimian County, Sichuan Province, China (Figure 1), and lies in the southwest region of the Sichuan Basin; the largest tributary of the Minjiang River; and the southeast region of the Gongga Mountains. The reserve covers an area of about 47,885 hm2; the elevation of the reserve is from 1,330 m to 4,551 m. The annual mean temperature and annual mean precipitation are 11.7°C–14.4°C and 800–1,250 mm, respectively. As the altitude increases, the vegetation types of the reserve transition from evergreen broadleaved forest and evergreen broadleaved mixed forest to coniferous and deciduous broadleaved mixed forest and evergreen coniferous forest, and subsequently, to alpine shrub and alpine meadow (Hong et al., 2015). AS, YL, and FF are the three main species foraged by giant pandas in the reserve. The most widely distributed species is AS, accounting for 38.08% of the total area of bamboo for giant pandas in this mountainous region (Sichuan Forestry Bureau, 2015). YL and FF account for 28.02% and 12.48% of giant pandas’ diets, respectively (Sichuan Forestry Bureau, 2015). AS is mainly distributed in regions higher than 2,500 m above sea level, and YL and FF are mainly distributed below 2,800 m above sea level. AS was the main food year-round, YL was partial food for winter, and FF was occasional food (Xie and Hong, 2020; Long et al., 2021).

Study area in Liziping National Nature Reserve, China.

 Experimental design

Field surveys and sampling were conducted in May (spring) and October (autumn) in 2020 (Figure 1). First of all, 4 transect lines were set up in AS, YL, and FF forests, respectively. The transect lines were set along the elevation gradient, and each transect line was separated by more than 200 m. Second, 3–5 survey plots (20 m × 20 m) were set on each transect line, and adjacent survey plots with an altitude distance of more than 50 m were arranged on the same transect. The bamboo species, latitude and longitude, altitude, tree layer, shrub layer, and other related variables were recorded and measured in each survey plot. Then, 1 bamboo plot (1 m × 1 m) was set at the central point of each survey plot and 5 m east and south of the central point, and the related variable data of bamboo layers were measured and recorded (Supplementary Table S1).

 Sample collection and DNA extraction

For each survey plot, no <200 g of mixed bamboo leaves was collected using sterile gloves. Samples were immediately transported to the laboratory within 2 h after collection and stored in a −20°C freezer. Every 200 g sample was sterile transferred to a plastic bag (24 cm × 35 cm) containing 200 ml of sterile precooled TE-buffer (10 mM Tris, 1 mM EDTA, and pH 7.5); subsequently, 0.05% Tween-80 was added (Hong et al., 2017). The samples were cleaned to collect microbiota; each sample was shaken, eddy current, and ultrasonically stirred in a TE-buffer for 5 min, and the plastic bags were then stored in ice water (at ~4°C) for each treatment step. An ultrasonic scrubber with a frequency of 40 kHz (Shanghai Kudos Instrument Co., Shanghai, China) was used to remove microorganisms from bamboo leaves. The leaves were filtered through three layers of sterile nylon mesh for cell suspension. After filtration, the cell suspension was placed in four 50 ml tubes and centrifuged at 4°C for 15 min with the rate of 2,200 × g to remove the supernatant and obtain the fungus slime. Multi-tube mycelium was collected in a 2.0 ml reaction tube and washed twice with TE-buffer. The mud was immediately frozen at −80°C until DNA was extracted.

The DNA was extracted using the E.Z.N.A.TM Soil DNA Kit (Omega, Norcross, GA, United States) using the manufacturer’s specifications; however, the following modifications were made: the frozen cell pellets were suspended in a 1 ml Kit-supplied SLX Mlus buffer containing 500 mg glass beads and ground at 65 Hz for 90 s on Tissue Lyser-24 (Jingxin, China). Cell fragment suspensions were immediately disposed of in accordance with instructions in the kit manual. Finally, the total DNA was obtained using an elution buffer of 50 μl.

 ITS1 amplification, quantification, and sequencing

The first internal transcribed spacer (ITS1) region of the fungus was amplified. Polymerase chain reaction (PCR) amplified 35 cycles of the target marker gene. Error-corrected 12 bp barcode primers were used for each sample to allow sample multiplexing. PCR products from all samples were quantified using PicoGreen dsDNA analysis and were aggregated at equimolar concentrations. Each library was submitted to Mega Biotech on the Illumina MiSeq PE300 platform. The raw reads have been deposited in the NCBI Sequence Read Archive (SRA) database, with accession number PRJNA718425.

 Sequence data analysis

FASTP (V0.20.0) was used to perform multiplex isolation and quality screening of the original ITS1 region sequencing reads (Chen et al., 2018) and using FLASH (V1.2.7) for merging. UPARSE (V7.1) was used to cluster the operational taxonomic units (OTU) of the sequences based on 97% similarity, and single sequences and chimera were removed. Each OTU representative sequence was classified and analyzed using the RDP Classifier (V2.2) and ITS database (e.g., Unite V8.0), with a confidence threshold of 0.7. To better convey the biological information of these samples, bar charts were used to visualize the average relative abundance of fungal communities at the phylum and genus levels.

 Statistical analysis

First, the Kruskal–Wallis H test and Wilcoxon rank-sum test were used to detect differences in the abundance of dominant fungi at the phylum and genus levels, respectively. This was conducted among different bamboo species in the same season and the same bamboo species in different seasons. We used the linear discriminant analysis (LDA) effect size (LEfSe; Segata et al., 2011) to identify the taxa that characterized the differences among three bamboo species (LDA score ≥ 3.5). The Sobs index and Shannon index were calculated for each sample using Motherur (V1.30.1) to estimate the species abundance and diversity of different bamboo species in spring and autumn, and the Wilcoxon rank-sum test was used to determine the significant level. Then, a Principlal Coordinate Analysis (PCoA) distance matrix based on the (un)weighted UniFrac method was used to evaluate the fungal community structure. Permutational MANOVA (PERMANOVA) was used to analyze the effects of different bamboo species and seasons on the fungal community structure based on the Bray–Curtis distance matrix.

Secondly, using the Spearman correlation analysis method, analysis of the relative abundance of fungal community and environmental factors (Supplementary Table S1). To investigate the spatial correlation between environmental factors and fungal communities, the Mantel test (based on the OUT level) was used to analyze the spatial correlation between the distance matrix of environmental factors and the UniFrac distance matrix of fungi. The canonical correspondence analysis (CCA) was used to determine samples, environmental factors, and the relationship between the fungi flora. The relationship between Alpha (the Sobs index and Shannon index) and Beta diversity (Bray-Curtis distance) and environmental factors was analyzed by the linear regression method.

Finally, the FUNGuild (Fungi Functional Guild) prediction method was used to predict the functional composition of fungi (Nguyen et al., 2016).

 Metagenome sequencing and analysis

Whole-genome shotgun sequencing using the Illumina sequencing HiSeq platform. The raw sequencing data were filtered for low-quality and host genomic bases using MetaGene (Sun et al., 2018). Through the CD-HIT (V4.6.1) for the sample prediction of gene sequence clustering (parameters: 90% identity and 90% coverage), non-redundant gene sets were constructed. SOAPaligner (V2.21) was used (Li et al., 2008; Fu et al., 2012), mapped the quality read data of each sample to the non-redundant gene set (parameter: 95% identity), and counted the abundance of genes in the sample. Non-redundant genes sets were compared with the NCBI NR database using DIAMOND (V0.8.35; Buchfink et al., 2015; parameter: Blastp. E-value ≤ 1E-5) to obtain the taxonomic annotation from the NR library. DIAMOND (V0.8.35) was used to map the non-redundant gene sets to Kyoto Encyclopedia of Genes and Genomes database to annotate the KEGG function. Non-redundant gene sets were compared with the CAZy database using HMMER (V3.1b2; parameter: Hmmscan. E-value ≤ 1E-5) to annotate carbohydrate-active enzyme. DIAMOND (V0.8.35) was used to map the non-redundant gene sets into CARD to obtain the annotation results for fungal antibiotic resistance genes (ARGs) in AS phyllospheres.

LEfSe was used to determine the relative abundance characteristics of the KEGG Level 2 (LDA score ≥ 2), CAZyme (LDA score ≥ 3), and CARD databases (LDA score ≥ 3.5) of the AS phyllosphere fungal community between spring and autumn (Segata et al., 2011).

 Co-occurrence network analysis

Using the OTU data of bacteria and fungi obtained by sequencing, OTUs with relative abundances >0.1% were selected. Network analyses were performed in R. Gephi (Version 0.9.2) used to construct microbial co-occurrence networks among AS, YL, and FF in spring and autumn. Antibiotic resistance genes with relative abundances >0.1% were selected, and network analyses were also performed in R. Gephi (Version 0.9.2) used to construct microbial co-occurrence network in AS. If the Spearman’s correlation coefficient |ρ| is not <0.8, and p values are no more than 0.01, suggests that a correlation between two items was considered statistically robust (Quintela-Baluja et al., 2019).

 Results Phyllosphere fungal alpha and beta diversity

A maximum of 19 samples were collected from each bamboo species in each study season (Supplementary Table S2). A total of 101 bamboo samples and 6,749,062 reads were obtained. The average target ITS1 reads of each sample were 66,822 ± 6,323. A total of 8,769 OTUs were obtained by clustering 97% similarity, among which 3,278 OTUs were shared among the three bamboo species (Supplementary Figure S1A). The fungal OTUs in the three bamboo species in autumn were significantly higher than those in spring (Supplementary Figures S1B, S2). Among the three bamboo species, FF had the highest number of OTUs in spring and autumn (Supplementary Figure S1A).

Taxonomic analysis classified fungal OTUs into 9 phyla and 1,053 genera. Most of the phyllosphere fungal community belonged to the phyla Ascomycota in all samples (Figure 2A). The rest comprised Basidiomycota and other phyla. The relative abundance of the phyla Ascomycota and Basidiomycota were not significantly different among the three bamboo species (Figure 2A, Supplementary Figures S4A,B; Supplementary Table S3). The relative abundance of the phylum Basidiomycota in YL phyllosphere was higher in spring than in autumn, and the relative abundance of the phylum Ascomycota in YL phyllosphere was lower in spring than in autumn (Figure 2A; Supplementary Table S3).

Taxonomic profiles of fungal communities among AS, YL, and FF phyllosphere. (A) Relative abundance of different phyllosphere fungi phyla among AS, YL, and FF in spring and autumn. (B) LEfSe analysis of phyllosphere fungi among AS, YL, and FF (LDA score ≥ 3.5). The Sobs index (C) and Shannon index (D) of AS, YL, and FF phyllosphere fungi in spring and autumn.

At the genus level, unclassified_p__Ascomycota, Trichomerium, unclassified_c__Sordariomycetes, unclassified_o__Helotiales, unclassified_o__Pleosporales, and Lapidomyces were the dominant phyllosphere fungi (Supplementary Figure S3). In spring, the relative abundance of unclassifed_p__Ascomycota was significantly higher in the YL phyllosphere (Supplementary Figure S4C; Supplementary Table S3). The relative abundance of Trichomerium was significantly lower in the AS phyllosphere, but a contrasting pattern was observed for Lapidomyces (Supplementary Figure S4C; Supplementary Table S3). In autumn, the relative abundance of unclassified_p__Ascomycota and unclassified_o__Helotiales were significantly higher in the YL phyllosphere (Supplementary Figure S4D; Supplementary Table S3). The highest relative abundance of unclassified_o__Pleosporales existed in the AS phyllosphere, but a contrasting pattern was observed for unclassified_c__Sordariomycetes (Supplementary Figure S4D; Supplementary Table S3). For all three bamboo species, the relative abundance of unclassified_c__Sordariomycetes in autumn was significantly lower than that in spring (Supplementary Table S3).

LEfSe analysis revealed that the LDA scores of microbial taxa differed significantly among AS, YL, and FF (Figure 2B). Four fungi groups (c__Dothideomycetes, g__Lapidomyces, f__Teratosphaeriaceae, and g__unclassified_f__Didymellaceae) in the AS phyllosphere were significantly enriched in comparison to the other two bamboo species. In YL phyllosphere, only o__Lecanorales was significantly enriched in comparison to the other two bamboo species. In FF phyllosphere, six groups of fungi (including f__Mycosphaerellaceae, g__unclassified_o__Capnodiales, f__unclassified_o__Capnodiales, g__Apiospora, g__Pezizella, and f__Bulleribasidiaceae) were significantly enriched in comparison to the other two bamboo species (Figure 2B).

The Sobs index and Shannon index were found to decrease from autumn to spring in all bamboo species (Figures 2C,D; Supplementary Table S4). AS phyllosphere fungal community had significantly lowest the Sobs index and Shannon index in spring (Figures 2C,D; Supplementary Table S4). In autumn, the Sobs index and Shannon index in FF phyllosphere were significantly higher than those of AS and YL (Figures 2C,D; Supplementary Table S4).

The PCoA results showed that the samples were clustered by season and species (Figure 3). Based on weighted UniFrac, three bamboo species clustered together in the same season, FF and YL clustered closer together. Based on unweighted UniFrac, the three bamboo species are easier to distinguish than the weighted UniFrac (Figure 3). Further, the PERMANOVA based on the Bray-Curtis distance matrix showed that there were significant differences in the phyllosphere fungal community structure between the three bamboo species in spring and autumn (Table 1).

PCoA of phyllosphere fungal communities based on the (un)weighted UniFrac distance matrixes (A,B) among AS, YL, and FF in spring and autumn.

PERMANOVA exploring the effects of three bamboo species and two seasons based on Bray–Curtis distance matrixes.

Grouping factors Sums of Sqs Mean Sqs F model R2 Value of p
AS (Spring:Autumn) 2.298 2.298 11.947 0.260 0.001
YL (Spring:Autumn) 2.017 2.017 10.813 0.258 0.001
FF (Spring:Autumn) 1.611 1.611 7.132 0.192 0.001
Spring:Autumn 12.109 2.422 12.047 0.388 0.001
Autumn (AS:YL:FF) 3.568 1.784 8.845 0.269 0.001
Spring (AS:YL:FF) 4.843 2.422 12.087 0.340 0.001
Ecological factors influencing bamboo phyllosphere fungal communities

Mantel analysis results showed that elevation had the greatest effect on the phyllosphere fungal community (R = 0.617, p = 0.001; Supplementary Table S5). Furthermore, CCA showed that elevation (E), mean base diameter of bamboo (MBDB), mean height of bamboo (MHB), bamboo deaths (BD), tree diameter at breast height (TDBH), shrub number (SN), total number of live bamboo (TNLB), number of trees (NT), shrub coverage (SC), slope (S), tree height (TH), canopy density (CD), and water source distance (DW) also influenced the fungal community (Figure 4; Supplementary Table S6).

CCA analysis of environmental factors and phyllosphere fungal community among AS, YL, and FF. E, elevation; DW, water source distance; TDBH, trees diameter at breast height; TNLB, total number of live bamboo; SC, shrub coverage; BD, bamboo deaths; SN, shrubs numbers; TH, trees height; SH, Shrub height; CD, canopy density; MBDB, mean base diameter of bamboo; NT, number of trees; MHB, mean height of bamboo; S, Slope; BC, bamboo coverage.

The Spearman’s correlation heatmap showed that Basidiomycota abundance was significantly negatively correlated to shrub number (SN; Figure 5A). Ascomycota abundance was significantly positively correlated to tree height (TH), tree diameter at breast height (TDBH), shrub coverage (SC), shrub number (SN), total number of live bamboo (TNLB), and bamboo deaths (BD; Figure 5A). At the genus level, Lapidomyces exhibited significantly positive correlations with elevation (E) and bamboo coverage (BC), but significantly negative correlations with tree height (TH), number of trees (NT), mean height of bamboo (MHB), and mean base diameter of bamboo (MBDB; Figure 5B). unclassified_o__Helotiales and unclassified_o__Pleosporales were both positively correlated with tree height (TH), tree diameter at breast height (TDBH), shrub coverage (SC), shrub number (SN), total number of live bamboo (TNLB), and bamboo deaths (BD). unclassified_c__Sordariomycetes exhibited significantly negative correlations with shrub coverage (SC), shrub number (SN), total number of live bamboo (TNLB), and bamboo deaths (BD), but significantly positive correlations with mean height of bamboo (MHB). Trichomerium exhibited significant positive correlations with bamboo coverage (BC), mean height of bamboo (MHB), and mean base diameter of bamboo (MBDB), but negative correlations with elevation (E), tree diameter at breast height (TDBH), and shrub coverage (SC). unclassified_p__Ascomycota was significantly positively correlated to water source distance (DW) and total number of live bamboo (TNLB), but negatively correlated to number of trees (NT; Figure 5B).

Correlation heatmap showed the relative abundance of fungi phylum (A) and genus (B) in relation to environmental factors. R is shown in different colors, and the legend on the right is the color range of different R values. E, elevation; DW, water source distance; TDBH, trees diameter at breast height; TNLB, total number of live bamboo; SC, shrub coverage; BD, bamboo deaths; SN, shrubs numbers; TH, trees height; SH, Shrub height; CD, canopy density; MBDB, mean base diameter of bamboo; NT, number of trees; MHB, mean height of bamboo; S, Slope; BC, bamboo coverage. (*0.01 < p < 0.05, **0.001 < p < 0.01, ***p < 0.001).

Further linear regression revealed that elevation had the highest effect for the Bray–Curtis distance, but the highest ecological influence factor on the Sobs index and Shannon index was mean base diameter of bamboo (Supplementary Table S7).

 Co-occurrence patterns of fungal interactions with bacteria community

The results of co-occurrence network analysis showed that the network complexity, node connectivity, and average degree of the F. ferax microbial community were higher than those of A. spanostachya and Y. lineolate (Figure 6; Supplementary Table S8). In autumn, A. spanostachya communities employed a much more complex network than that in spring (Figures 5A,B), which was reflected in the node and edge number, as well as average degree, which was higher in autumn than in spring (Supplementary Table S8). However, the network complexity of F. ferax and Y. lineolate had no obvious differences in the two seasons (Figures 6C–F). F. Ferax microbial networks were mainly reflected in the interactions within the bacterial domain. A. spanostach and Y. lineolate microbial networks were mainly reflected in the interactions within the fungal domain (Figure 6; Supplementary Table S8).

Co-occurrence network analysis showing the fungi interactions with bacteria among AS, YL, and FF phyllosphere in spring and autumn (A–F).

 Function of the fungi (e.g., KEGG, CAZy, and ARGs)

The functional classification of fungi was performed based on the Fungi Functional Guild (Nguyen et al., 2016). Three ecological guilds – pathotroph, saprotroph, and symbiotroph – for ~40% of all fungal OTUs were detected (Supplementary Figure S5).

In order to further explore the variation of functional potential in AS phyllosphere, we screened 8 samples (spring: 4 samples; autumn: 4 samples) for whole-genome shotgun sequencing. A total of 3.91 million contigs and 2.47 Gb of assembly sequence were obtained, and the average contig N50 was 655 bp (Supplementary Table S9). The metagenomic analysis confirmed 402 KOs, including 46 KEGG Level 2 categories in fungi. In the KEGG level 2 categories, the relative abundance of global and overview maps (20%) was the highest in A. spanostachya, followed by neurodegenerative disease (7%), signal transduction (6%), and carbohydrate metabolism (4%; Supplementary Figure S6A). In spring, AS phyllosphere fungi showed high abundances in KEGG Level 2 categories of neurodegenerative disease, energy metabolism, environmental adaptation, and drug resistance: antimicrobial, endocrine, and metabolic disease (Figure 7A), and the functions of AS phyllosphere fungi are primarily a consequence of the fungi of the Beauveria and Cladosporium genus (Supplementary Table S10). In contrast, in autumn, translation; folding; sorting and degradation; endocrine systems; transcription; metabolism of cofactors and vitamins; immune systems; nucleotide metabolism; infectious diseases: bacterial, signaling molecules, and interactions; immune disease; cellular community-eukaryotes; metabolism of terpenoids and polyketides; and infectious diseases: parasitic, exhibited higher abundance (Figure 7A), and the functions of AS phyllosphere fungi are mainly a consequence of fungi of the Beauveria and Lachnellula genus (Supplementary Table S10).

LEfSe analysis of KEGG level 2 (A) (LDA score ≥ 2) and CAZy (B) (LDA score ≥ 3) of AS phyllosphere fungal community between spring and autumn. GH, glycoside hydrolases; GT, glycosyl transferases; CE, carbohydrate esterases; AA, auxiliary activities.

The relative abundance of CAZy in glycoside hydrolases (43%) was the highest in AS phyllosphere fungi which contained some cellulases (GH5, GH6, GH7, GH28, etc.), followed by glycosyl transferases (27%), auxiliary activities (16%), and carbohydrate esterases (11%; Supplementary Figure S6B; Supplementary Table S10). According to the LEfse results for the CAZyme, three families had significantly higher relative abundances in spring, including glycoside hydrolases (GHs), glycosyl transferases (GTs), and carbohydrate esterases (CEs; Figure 7B), and the functions of AS phyllosphere fungi are mainly a consequence of fungi of the Chaetothyriales and Hypocreales order (Supplementary Table S11). In autumn, four families had significantly higher relative abundances: GHs, CEs, GTs, and AAs, and the functions of AS phyllosphere fungi are mainly a consequence of fungi of the Chaetothyriales and Pleosporales orders (Supplementary Table S10). Interestingly, GHs had the highest differences between spring and autumn (Figure 7B).

A total of 125 antibiotic resistance genes were detected in eight bamboo samples from the AS phyllosphere fungi based on metagenomic sequencing (Supplementary Table S11). At the genus level, g__unclassified_o__Helotiales, g__Epicoccum, g__Mortierella, and g__Paraphaeosphaeria had multidrug resistance genes; g__unclassified_k__Fungi had multidrug and tetracycline resistance genes; g__Periconia had aminocoumarin resistance genes; g__Lachnellula had multidrug, tetracycline, and rifamycin resistance genes; g__Beauveria had multidrug, tetracycline, sulfonamide, fluoroquinolone, pleuromutilin, and beta-lactam resistance genes (Supplementary Table S11). Multidrug (42%) and tetracycline (24%) resistance genes were the most common types of ARGs in all samples (Figure 8A). The abundance of multidrug, tetracycline, and glycopeptide resistance genes was high and closely correlated with other ARGs (Figure 9). Beta-lactam resistance gene had only one edge to multidrug resistance genes (Figure 9). The remaining ARGs were glycopeptide, aminocoumarin, mupirocin, fluoroquinolone, rifamycin, etc. According to the LEfse results of the CARD, in spring, AS fungi displayed higher abundances of beta-lactam and phenicol resistance genes than in autumn; however, glycopeptide and bicyclomycin resistance genes exhibited higher abundance in autumn (Figure 8B).

Fungal antibiotic resistance genes (ARGs) in AS phyllosphere. Relative abundance of fungal antibiotic resistance genes (ARGs) in AS phyllosphere between spring and autumn (A). LEfSe analysis of fungal antibiotic resistance genes (ARGs) in AS phyllosphere between spring and autumn (B). (LDA score ≥ 3.5).

ARGs co-occurrence network of fungal resistance genes in AS phyllosphere. Nodes represent genes and are colored to indicate antimicrobial classes of ARGs.

 Discussion

In our study, Ascomycota and Basidiomycota were found to be the dominant phyla in AS, YL, and FF phyllosphere fungi microflora, which is similar to several previous studies (Helander et al., 2013; Xu et al., 2014; Jin et al., 2020). The fungal OTUs richness and diversity (e.g., Sobs index and Shannon index) of the three bamboo species in spring were significantly lower than those in autumn. This result matches that found by Zheng (2011), who determined that the diversity of Pinus tabulaeformis phyllosphere microbial communities exhibited the greatest diversity in autumn, followed by summer and spring. Further, Thompson et al. (1993) also found that the diversity of phyllosphere fungal communities of Beta vulgaris in spring was lower than that in autumn. The season is an important factor affecting phyllosphere microorganisms, whose diversity is higher in autumn than in spring because of higher temperature and humidity in autumn (Laforest et al., 2016).

The season is also a stronger driver of the phyllosphere fungal community structure than the host bamboo species based on the (un)weighted UniFrac distance matrix (Figure 3). However, this result differs from that of Laforest et al. (2016), who found that host species are a more important factor than the site or time in terms of determining the community of temperate tree phyllospheres.

The Mantel test found that elevation, water source distance, tree height, etc. significantly affected the phyllosphere fungal community, and effect of elevation is the most obvious (Supplementary Table S5), which is similar to the result of Long et al. (2021), who found that elevation, water source distance, etc. significantly impact the phyllosphere bacterial community. Similarly, Cordier et al. (2012) found that elevation affects the relative abundance and community composition of Fagus sylvatica fungal communities. Furthermore, CCA found that slope significantly affected the phyllosphere fungal community. Similarly, Bari et al. (2019) explored the effects of elevation, the field slope, site direction, and fungal spore fruit position height on trees, finding that the field slope was the main factor influencing the fungal spores abundance.

To explore the effect of the micro-ecological environment (biotic factor) on the bamboo phyllosphere fungal community, a co-occurrence network analysis was performed, which showed that FF had higher network complexity than AS and YL. This may be related to the higher Sobs and Shannon indices in FF phyllosphere microbial communities, because the higher Sobs and Shannon indices would increase the possibility of interactions between different microbial communities (Wang B. et al., 2022). This result was similar to those of Wang B. et al. (2022), who found that the high alpha diversity of microbial community in a given group provides additional potential interactions between groups. AS had higher network complexity in autumn than in spring. However, networks of FF and YL slightly differed in spring and autumn (Figure 6). It is possible that AS occurs at high altitudes and has more distinct climatic environments, in terms of both temperature and light. There are complex links between bacteria and fungi, bacteria and bacteria, and fungi and fungi in AS, YL, and FF. Bacteria and fungi usually have symbiotic relationships, with both positive and negative effects. For example, some bacteria stimulate mycorrhizal formation, and fungi and bacteria grow faster when they are present together (Bengtsson, 1992; Founoune et al., 2002; Mille-Lindblom et al., 2006); bacterial volatiles inhibit the growth of fungal spores, and fungal volatiles inhibit not only the growth of bacteria, but also the growth of other fungi (Schmidt et al., 2015).

In this study, the functions of AS phyllosphere fungi were mainly dominated by Beauveria (Supplementary Table S10), a genus of entomogenic fungi that plays an important role in pest control and has no pathogenic effect on humans and animals (Mascarin and Jaronski, 2016). Beauveria not only protects giant panda’s food resources but can also kill giant panda’s ectoparasites. Carbohydrate metabolism and GH enzymes are mainly a consequence of Eurotiomycetes and Dothideomycetes (Supplementary Table S10). It has been demonstrated that Eurotiomycetes and Dothideomycetes consume lignin when it is the only carbon source (Ferrari et al., 2021). The majority of vitamins come from an animal’s diet. Both the giant panda and red panda forage on AS as a staple food in Liziping NR (Wei et al., 1999). The abundant metabolism of cofactors and vitamins may be beneficial to the physical and mental health of giant pandas and red pandas (Badawy, 2014; Rudzki et al., 2021). Further, the abundant energy metabolism and environmental adaptation in spring indicated that spring was not a suitable season for phyllosphere fungi (Laforest et al., 2016). More abundant translation, folding, sorting and degradation, transcription, nucleotide metabolism, signaling molecules and interactions, and cellular community-eukaryotes indicated fungal cell viability and cell proliferation in autumn.

The LEfse results of the CAZyme show that glycoside hydrolases (GHs) exhibit the main differences, which supports the fact that phyllosphere microorganisms rely primarily on carbohydrate metabolism for nutrients (Kawaguchi et al., 2011; Ryffel et al., 2016). GHs, which can decompose complex carbohydrates, are a highly diverse group of key enzymes related to gut microbiota and their metabolic function (Lee et al., 2013). Zhu et al. (2011) found that the gut microbiome of giant pandas contains a small number of cellulases, and giant pandas and red pandas may obtain cellulose-digesting flora from phyllosphere microorganisms by eating bamboo leaves in large quantities (Jin et al., 2020).

In the AS phyllosphere, multidrug and tetracycline resistance genes were the most dominant ARGs. ARG co-occurrence network analysis showed that beta-lactam resistance genes were rarely correlated with other ARGs (Figure 9), reducing the possibility of co-selection of other antimicrobials (Yang et al., 2022). Many phyllosphere fungi have multiple resistance genes, proving that phyllospheres have become reservoirs for resistance genes (Sun et al., 2021). Multidrug resistance genes have been shown to be the most abundant ARGs in some vegetables’ phyllosphere (Blau et al., 2018; Yin et al., 2022). Multidrug resistance to many different classes of antibiotics, which impedes antibiotic treatment, is a huge risk to human life (Zhu et al., 2013; Ouyang et al., 2015), and is also a potential threat to giant panda health. Several recent studies have shown that diverse ARGs exist in the gut of wild and captive giant pandas, and multidrug resistance genes were the dominant type of ARGs in the giant panda gut (Guo et al., 2019; Hu et al., 2020; Mustafa et al., 2021; Zhu et al., 2021). A survey based on ARGs of human gut microbiota in three countries (China, Denmark, and Spain) found that tetracycline resistance genes were the most abundant group of ARGs, followed by MLSs and beta_lactam resistance genes (Hu et al., 2013). Given the abundance and diversity of ARGs in the bamboo phyllosphere, the diet has proven to be the main driving factor for ARG variation in the gut microbiome of giant pandas (Zhu et al., 2021).

 Conclusion

Overall, our results indicated that the fungal community’s abundance, diversity, and community structure were mainly affected by the season, host species, and elevation. The season and host species are major factors affecting biological functions (KEGG and CAZyme), ARGs, and interactions between sympatric bacterial and fungal communities in bamboo phyllospheres. As the most important food resource for giant pandas in Liziping NR (Xiaoxiangling mountains; Hong et al., 2015; Sichuan Forestry Bureau, 2015; Xie and Hong, 2020), AS phyllosphere fungi contain many glycoside hydrolases (e.g., cellulases CH5, GH6, GH7, etc.), which may enter the giant pandas’ intestine through ingestion, helping them to digest cellulose and obtain nutrition. The functions of AS phyllosphere fungi were mainly dominated by Beauveria, which can be used for biological control of pests and diseases (Mascarin and Jaronski, 2016). Therefore, food microbes (bamboo phyllosphere microorganisms) may be an important technique for giant pandas’ gut microbes to degrade cellulose and also protect food resources of giant pandas (Zhu et al., 2011). Although various ARGs were found in AS phyllosphere fungi flora, these ARGs could naturally occur in bamboo phyllosphere fungi flora (D’Costa et al., 2011; Sun et al., 2021). Therefore, numerous comparative studies on the effects of human activities on bamboo phyllosphere microorganism ARGs should be carried out in giant panda habitats to predict the risk of antibiotics for wild pandas.

 Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

 Author contributions

LK: conceptualization, methodology, investigation, software, resources, formal analysis, and writing original draft. WL, QD, HZ, WW, JT, HH, YY, and JL: investigation, validation and data curation. ZZ: conceptualization, supervision, and writing–review and editing. MH: conceptualization, supervision, funding acquisition, writing–original draft, and review and editing. All authors contributed to the article and approved the submitted version.

 Funding

This study was supported by the National Natural Science Foundation of China (grant no. 31900337), The Second Tibetan Plateau Scientific Expedition (grant no. 2019QZKK05010502), and the Innovation Team Funds of China West Normal University (grant no. KCXTD2022-7).

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

We thank the staff of Liziping National Nature Reserve Administration for their assistance in data collection and experiments.

 Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2022.1009588/full#supplementary-material

 References Arnold A. E. Mejía L. C. Kyllo D. Rojas E. I. Maynard Z. Robbins N. . (2003). Fungal endophytes limit pathogen damage in a tropical tree. Proc. Natl. Acad. Sci. 100, 15649–15654. doi: 10.1073/pnas.2533483100, PMID: 14671327 Badawy A. A. (2014). Pellagra and alcoholism: a biochemical perspective. Alcohol Alcohol. 49, 238–250. doi: 10.1093/alcalc/agu010, PMID: 24627570 Bari E. Aghajani H. Ohno K. Shahi R. Hale M. Bahmani M. (2019). Ecology of wood-inhabiting fungi in northern forests of Iran. For. Path. 49:e12501. doi: 10.1111/efp.12501 Bengtsson G. (1992). Interactions between fungi, bacteria and beech leaves in a stream microcosm. Oecologia 89, 542–549. doi: 10.1007/BF00317161, PMID: 28311885 Blau K. Bettermann A. Jechalke S. Fornefeld E. Vanrobaeys Y. Stalder T. . (2018). The transferable Resistome of produce. mBio. 9, e01300–e01318. doi: 10.1128/mBio.01300-18, PMID: 30401772 Boer W. Folman L. B. Summerbell R. C. Boddy L. (2005). Living in a fungal world: impact of fungi on soil bacterial niche development. FEMS Microbiol. Rev. 29, 795–811. doi: 10.1016/j.femsre.2004.11.005, PMID: 16102603 Buchfink B. Xie C. Huson D. H. (2015). Fast and sensitive protein alignment using DIAMOND. Nat. Methods 12, 59–60. doi: 10.1038/nmeth.3176, PMID: 25402007 Chen Q. L. Cui H. L. Su J. Q. Penuelas J. Zhu Y. G. (2019). Antibiotic resistomes in plant microbiomes. Trends Plant Sci. 24, 530–541. doi: 10.1016/j.tplants.2019.02.010, PMID: 30890301 Chen S. Zhou Y. Chen Y. Gu J. (2018). Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, i884–i890. doi: 10.1093/bioinformatics/bty560, PMID: 30423086 Copeland J. K. Yuan L. Layeghifard M. Wang P. W. Guttman D. S. (2015). Seasonal community succession of the phyllosphere microbiome. Mol. Plant-Microbe Interact. 28, 274–285. doi: 10.1094/MPMI-10-14-0331-FI, PMID: 25679538 Cordier T. Robin C. Capdevielle X. Fabreguettes O. Desprez-Loustau M. L. Vacher C. (2012). The composition of Phyllosphere fungal assemblages of European beech (Fagus Sylvatica) varies significantly along an elevation gradient. New Phytol. 196, 510–519. doi: 10.1111/j.1469-8137.2012.04284.x, PMID: 22934891 Corinne V. Arndt H. Annabel J. P.. Ursula S. Stéphane C. Cindy E. (2016). The Phyllosphere: Microbial Jungle at the Plant–Climate Interface. Annu. Rev. Ecol. Evol. S. 47, 1–24. doi: 10.1146/annurev-ecolsys-121415-032238, PMID: 9950422 D’Costa V. M. King C. E. Kalan L. Morar M. Sung W. W. Schwarz C. . (2011). Antibiotic resistance is ancient. Nature 477, 457–461. doi: 10.1038/nature10388 Endo H. Yamagiwa D. Hayashi Y. Koie H. Yamaya Y. Kimura J. (1999). Role of the giant panda’s‘pseudo-thumb’. Nature 397, 309–310. doi: 10.1038/16830, PMID: 9950422 Fan X. Matsumoto H. Wang Y. Hu Y. Liu Y. Fang H. . (2019). Microenvironmental interplay predominated by beneficial Aspergillus abates fungal pathogen incidence in Paddy environment. Environ. Sci. Technol. 53, 13042–13052. doi: 10.1021/acs.est.9b04616, PMID: 31631659 Ferrari R. Gautier V. Silar P. (2021). Chapter three – lignin degradation by ascomycetes. Adv. Bot. Res. 99, 77–113. doi: 10.1016/bs.abr.2021.05.006 Finkel O. M. Burch A. Y. Lindow S. E. Post A. F. Belkin S. (2011). Geographical location determines the population structure in phyllosphere microbial communities of a salt-excreting desert tree. Appl. Environ. Microbiol. 77, 7647–7655. doi: 10.1128/AEM.05565-11, PMID: 21926212 Founoune H. Duponnois R. Bâ A. M. Sall S. Branget I. Lorquin J. . (2002). Mycorrhiza helper bacteria stimulate ectomycorrhizal symbiosis of Acacia holosericea with Pisolithus alba. New Phytol. 153, 81–89. doi: 10.1046/j.0028-646X.2001.00284.x Fu L. Niu B. Zhu Z. Wu S. Li W. (2012). CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics 28, 3150–3152. doi: 10.1093/bioinformatics/bts565, PMID: 23060610 Guo W. Mishra S. Wang C. D. Zhang H. M. Ning R. H. Kong F. . (2019). Comparative study of gut microbiota in wild and captive giant pandas (Ailuropoda melanoleuca). Genes 10:827. doi: 10.3390/genes10100827, PMID: 31635158 Helander M. Jia R. Huitu O. Sieber T. N. Jia J. Niemelä P. . (2013). Endophytic fungi and silica content of different bamboo species in giant panda diet. Symbiosis 61, 13–22. doi: 10.1007/s13199-013-0253-z Hong M. S. Peng G. X. Keyhani N. O. Xia Y. X. (2017). Application of the entomogenous fungus, Metarhizium anisopliae for leafroller (Cnaphalocrocis medinalis) control and its effect on rice phyllosphere microbial diversity. Appl. Microbiol. Biotechnol. 101, 6793–6807. doi: 10.1007/s00253-017-8390-6, PMID: 28695229 Hong M. S. Wei W. Yang Z. S. Yuan S. B. Yang X. Y. Gu X. D. . (2016). Effects of timber harvesting on Arundinaria spanostachya bamboo and feeding-site selection by giant pandas in Liziping nature reserve, China. For. Ecol. Manage. 373, 74–80. doi: 10.1016/j.foreco.2016.04.039 Hong M. S. Wei W. Zhou H. Tang J. F. Han H. Zhang Z. J. (2019). Creative conservation in China: releasing captive giant pandas into the wild. Environ. Sci. Pollut. Res. 26, 31548–31549. doi: 10.1007/s11356-019-06384-4, PMID: 31493077 Hong M. S. Yuan S. B. Yang Z. S. Yang X. Y. Gu X. D. Huang F. . (2015). Comparison of microhabitat selection and trace abundance of giant pandas between primary and secondary forests in Liziping nature reserve, China: effects of selective logging. Mamm. Biol. 80, 373–379. doi: 10.1016/j.mambio.2015.05.003 Horton M. W. Bodenhausen N. Beilsmith K. Meng D. Z. Muegge B. D. Subramanian S. . (2014). Genome-wide association study of Arabidopsis thaliana leaf microbial community. Nat. Commun. 5:5320. doi: 10.1038/ncomms6320, PMID: 25382143 Hu T. Dai Q. L. Chen H. Zhang Z. Dai Q. Gu X. D. . (2020). Geographic pattern of antibiotic resistance genes in the metagenomes of the giant panda. Microb. Biotechnol. 14, 186–197. doi: 10.1111/1751-7915.13655, PMID: 32812361 Hu Y. B. Wu Q. Ma S. Ma T. X. Shan L. Wang X. . (2017). Comparative genomics reveals convergent evolution between the bamboo-eating giant and red pandas. Proc. Nat. Acad. Sci. U. S. A. 114, 1081–1086. doi: 10.1073/pnas.1613870114, PMID: 28096377 Hu Y. F. Yang X. Qin J. J. Lu N. Cheng G. Wu N. . (2013). Metagenome-wide analysis of antibiotic resistance genes in a large cohort of human gut microbiota. Nat. Commun. 4:2151. doi: 10.1038/ncomms3151 Jin C. Z. Ciochon R. L. Dong W. Hunt R. M. Liu J. Y. Jaeger M. . (2007). The first skull of the earliest giant panda. Proc. Natl. Acad. Sci. U. S. A. 104, 10932–10937. doi: 10.1073/pnas.0704198104, PMID: 17578912 John H. A. Robin F. H. (2000). The Ecology and Biogeography of Microorganisms on Plant Surfaces. Annu. Rev. Phytopathol. 38, 145–1880. doi: 10.1146/annurev.phyto.38.1.145, PMID: 12676659 Jin L. Wu D. F. Li C. W. Zhang A. Y. Xiong Y. W. Wei R. P. . (2020). Bamboo nutrients and microbiome affect gut microbiome of giant panda. Symbiosis 80, 293–304. doi: 10.1007/s13199-020-00673-0 Kawaguchi K. Yurimoto H. Oku M. Sakai Y. (2011). Yeast methylotrophy and autophagy in a methanol-oscillating environment on growing Arabidopsis thaliana leaves. PLoS One 6:e25257. doi: 10.1371/journal.pone.0025257, PMID: 21966472 Kembel S. W. O’Connor T. K. Arnold H. K. Hubbell S. P. Wright S. J. Green J. L. (2014). Relationships between phyllosphere bacterial communities and plant functional traits in a neotropical forest. Proc. Natl. Acad. Sci. U. S. A. 111, 13715–13720. doi: 10.1073/pnas.1216057111, PMID: 25225376 Laforest I. Messier C. Kembel S. W. (2016). Host species identity, site and time drive temperate tree phyllosphere bacterial community structure. Microbiome 4:27. doi: 10.1186/s40168-016-0174-1, PMID: 27316353 Lee S. Cantarel B. Henrissat B. Gevers D. Birren B. W. Huttenhower C. . (2013). Gene-targeted metagenomic analysis of glucan-branching enzyme gene profiles among human and animal fecal microbiota. ISME J. 8, 493–503. doi: 10.1038/ismej.2013.167, PMID: 24108330 Lee S. S. Ha J. K. Cheng K.-J. (2000). Relative contributions of bacteria, protozoa, and fungi to in vitro degradation of orchard grass cell walls and their interactions. Appl. Environ. Microbiol. 66, 3807–3813. doi: 10.1128/aem.66.9.3807-3813.2000, PMID: 10966394 Li R. Li Y. Kristiansen K. Wang J. (2008). SOAP: short oligonucleotide alignment program. Bioinformatics 24, 713–714. doi: 10.1093/bioinformatics/btn025, PMID: 18227114 Lindow S. E. Brandl M. T. (2003). Microbiology of the phyllosphere. Appl. Environ. Microbiol. 69, 1875–1883. doi: 10.1128/AEM.69.4.1875-1883.2003, PMID: 12676659 Lindow S. E. Leveau J. H. J. (2002). Phyllosphere microbiology. Curr. Opin. Biotechnol. 13, 238–243. doi: 10.1016/s0958-1669(02)00313-0 Long J. J. Luo W. Xie J. M. Yuan Y. Wang J. Kang L. W. . (2021). Environmental factors influencing Phyllosphere bacterial communities in Giant pandas’ staple food bamboos. Front. Microbiol. 12:748141. doi: 10.3389/fmicb.2021.748141, PMID: 34803968 Ma B. Wang Y. Ye S. Liu S. Stirling E. Gilbert J. . (2020). Earth microbial co-occurrence network reveals interconnection pattern across microbiomes. Microbiome 8:82. doi: 10.1186/s40168-020-00857-2, PMID: 32498714 Mascarin G. M. Jaronski S. T. (2016). The production and uses of Beauveria bassiana as a microbial insecticide. World J. Microbiol. Biotechnol. 32:177. doi: 10.1007/s11274-016-2131-3, PMID: 27628337 Mille-Lindblom C. Fischer H. J. Tranvik L. (2006). Antagonism between bacteria and fungi: substrate competition and a possible tradeoff between fungal growth and tolerance towards bacteria. Oikos 113, 233–242. doi: 10.1111/j.2006.0030-1299.14337.x Mustafa G. R. Li C. W. Zhao S. Y. Jin L. He X. P. Shabbir M. Z. . (2021). Metagenomic analysis revealed a wide distribution of antibiotic resistance genes and biosynthesis of antibiotics in the gut of giant pandas. BMC Microbiol. 21:15. doi: 10.1186/s12866-020-02078-x, PMID: 33413128 Nguyen N. H. Song Z. Bates S. T. Branco S. Tedersoo L. Menke J. . (2016). FUNGuild: an open annotation tool for parsing fungal community datasets by ecological guild. Fungal Ecol. 20, 241–248. doi: 10.1016/j.funeco.2015.06.006 Ouyang W. Y. Huang F. Y. Zhao Y. Li H. Su J. Q. (2015). Increased levels of antibiotic resis-tance in urban stream of Jiulongjiang River, China. Appl. Microbiol. Biotechnol. 99, 5697–5707. doi: 10.1007/s00253-015-6416-5, PMID: 25661810 Paulino-Lima I. G. Azua-Bustos A. Vicuña R. González-Silva C. Salas L. Teixeira L. . (2013). Isolation of UVC-tolerant bacteria from the hyperarid Atacama Desert. Chile. Microb. Ecol. 65:325–335. doi: 10.1007/s00248-012-0121-z, PMID: 31351383 Pei M. Zhang B. He Y. Su J. Gin K. Lev O. . (2019). State of the art of tertiary treatment technologies for controlling antibiotic resistance in wastewater treatment plants. Environ. Int. 131:105026. doi: 10.1016/j.envint.2019.105026, PMID: 31351383 Porras-Soriano A. Soriano-Martín M. Porras-Piedra A. Azcón R. (2009). Arbuscular mycorrhizal fungi increased growth, nutrient uptake and tolerance to salinity in olive trees under nursery conditions. J. Plant Physiol. 166, 1350–1359. doi: 10.1016/j.jplph.2009.02.010, PMID: 19342122 Quintela-Baluja M. Abouelnaga M. Romalde J. Su J. Q. Yu Y. Gomez-Lopez M. . (2019). Spatial ecology of a wastewater network defines the antibiotic resistance genes in downstream receiving waters. Water Res. 162, 347–357. doi: 10.1016/j.watres.2019.06.075, PMID: 31295654 Rastogi G. Coaker G. L. Leveau J. H. (2013). New insights into the structure and function of phyllosphere microbiota through high-throughput molecular approaches. FEMS Microbiol. Lett. 348, 1–10. doi: 10.1111/1574-6968.12225, PMID: 23895412 Redford A. J. Fierer N. (2009). Bacterial succession on the leaf surface: a novel system for studying successional dynamics. Microb. Ecol. 58, 189–198. doi: 10.1007/s00248-009-9495-y, PMID: 19221834 Rehman F. U. Sultan A. Kalsoom M. Adnan M. Ilyas M. A. Hayat Q. . (2021). Chemistry of microbial life in Phyllosphere: a review. Green Rep. 2, 1–8. doi: 10.36686/Ariviyal.GR.2021.02.04.009 Rudzki L. Stone T. W. Maes M. Misiak B. Samochowiec J. Szulc A. (2021). Gut microbiota-derived vitamins-underrated powers of a multipotent ally in psychiatric health and disease. Prog. Neuro-Psychopharmacol. Biol. Psychiatry 107:110240. doi: 10.1016/j.pnpbp.2020.110240, PMID: 33428888 Ryffel F. Helfrich E. J. Kiefer P. Peyriga L. Portais J. Piel J. . (2016). Metabolic footprint of epiphytic bacteria on Arabidopsis thaliana leaves. ISME J. 10, 632–643. doi: 10.1038/ismej.2015.141, PMID: 26305156 Schmidt R. Cordovez V. de Boer W. Raaijmakers J. Garbeva P. (2015). Volatile affairs in microbial interactions. ISME J. 9, 2329–2335. doi: 10.1038/ismej.2015.42, PMID: 26023873 Segata N. Izard J. Waldron L. Gevers D. Miropolsky L. Garrett W. S. . (2011). Metagenomic biomarker discovery, and explanation. Genome Biol. 12:R60. doi: 10.1186/gb-2011-12-6-r60, PMID: 21702898 Sichuan Forestry Bureau (2015). The Giant Panda of Sichuan: Report of the Fourth National Giant Panda Survey. Sichuan Science and Technology Press, Chengdu. Sun Y. Y. Li M. Li Y. Y. Li L. X. Zhai W. Z. Wang P. . (2018). The effect of clostridium butyricum on symptoms and fecal microbiota in diarrhea-dominant irritable bowel syndrome: a randomized, double-blind, placebo-controlled trial. Sci. Rep. 8, 2964–2975. doi: 10.1038/s41598-018-21241-z, PMID: 29445178 Sun Y. Snow D. Walia H. Li X. (2021). Transmission routes of the microbiome and resistome from manure to soil and lettuce. Environ. Sci. Technol. 55, 11102–11112. doi: 10.1021/acs.est.1c02985, PMID: 34323079 Tang H. Zhang N. Sun B. Liang Y. (2020). Effect of interaction between arbuscular mycorrhizal fungi and rhizosphere bacteria in farmland soils on nutrients utilization. Acta Microbiol Sin. 60, 1117–1129. In Chinese. doi: 10.13343/j.cnki.wsxb.20200204 Thompson I. P. Bailey M. J. Fenlon J. S. Fermor T. R. Lilley A. K. Lynch J. M. . (1993). Quantitative and qualitative seasonal changes in the microbial community from the phyllosphere of sugar beet (Beta vulgaris). Plant Soil 150, 177–191. doi: 10.1007/BF00013015 Tun H. M. Mauroo N. F. Chan S. Y. Ho J. Wong M. T. Leung F. C. (2014). Microbial diversity and evidence of novel homoacetogens in the gut of both geriatric and adult giant pandas (Ailuropoda melanoleuca). PLoS One 9:e79902. doi: 10.1371/journal.pone.0079902, PMID: 24475017 U’Ren J. M. Lutzoni F. Miadlikowska J. Laetsch A. D. Arnold A. E. (2012). Host and geographic structure of endophytic and endolichenic fungi at a continental scale. Am. J. Bot. 99, 898–914. doi: 10.3732/ajb.1100459, PMID: 22539507 Wang B. Huang S. Y. Li Z. C. Zhou Z. C. Huang J. Y. Yu H. L. . (2022). Factors driving the assembly of prokaryotic communities in bulk soil and rhizosphere of Torreya grandis along a 900-year age gradient. Sci. Total Environ. 837:155573. doi: 10.1016/j.scitotenv.2022.155573, PMID: 35504392 Wang X. M. Su D. F. Jablonski N. G. Ji X. P. Kelley J. Flynn L. J. . (2022). Earliest giant panda false thumb suggests conflicting demands for locomotion and feeding. Sci. Rep. 12:10538. doi: 10.1038/s41598-022-13402-y, PMID: 35773284 Wei F. W. Feng Z. J. Wang Z. W. (1999). Habitat selection by giant pandas and red pands in XiangLing Mountains. Acta Zool. Sin. 45, 57–63. In Chinese Wei F. W. Hu Y. B. Yan L. Nie Y. G. Wu Q. Zhang Z. J. (2015). Giant pandas are not an evolutionary cul-de-sac: evidence from multidisciplinary research. Mol. Biol. Evol. 32, 4–12. doi: 10.1093/molbev/msu278 Wu Q. Wang X. Ding Y. Hu Y. B. Nie Y. G. Wei W. . (2017). Seasonal variation in nutrient utilization shapes gut microbiome structureand function in wild giant pandas. Proc. Biol. Sci. 284:20170955. doi: 10.1098/rspb.2017.0955, PMID: 28904137 Xie J. M. Hong M. S. (2020). Effects of ecological factors on growth of Arundinaria spanostachya shoots in Liziping National Nature Reserve, China. Glob. Ecol. Conserv. 23:e01121. doi: 10.1016/j.gecco.2020.e01121 Xu X. L. Bai Y. Liu H. Yang C. L. Ling Y. G. (2014). Microflora analysis on six kinds of bamboo leaves in Sichuan Province. J. Huazhong Agric. Univ. 33, 53–59. In Chinese Yang L. Shen Y. B. Jiang J. Y. Wang X. Y. Shao D. Y. Lam M. . (2022). Distinct increase in antimicrobial resistance genes among Escherichia coli during 50 years of antimicrobial use in livestock production in China. Nat. Food 3, 197–205. doi: 10.1038/s43016-022-00470-6 Yin Y. Zhu D. Yang G. Su J. Q., and Duan,G. L. (2022). Diverse antibiotic resistance genes and potential pathogens inhabit in the phyllosphere of fresh vegetables. Sci. Total Environ. 815:152851. doi: 10.1016/j.scitotenv.2021.152851, PMID: 34990692 Zahn G. Amend A. S. (2019). Foliar fungi alter reproductive timing and allocation in Arabidopsis under normal and water-stressed conditions. Fungal Ecol. 41, 101–106. doi: 10.1016/j.funeco.2019.04.002 Zhao D. Y. Liu P. Fang C. Sun Y. M. Zeng J. Wang J. Q. . (2013). Submerged macrophytes modify bacterial community composition in sediments in a large, shallow, freshwater lake. Can. J. Microbiol. 59, 237–244. doi: 10.1139/cjm-2012-0554, PMID: 23586747 Zhao H. B. Yang J. R. Xu H. L. Zhang J. Z. (2010). Pseudogenization of the umami taste receptor gene Tas1r1 in the giant panda coincided with its dietary switch to bamboo. Mol. Biol. Evol. 27, 2669–2673. doi: 10.1093/molbev/msq153, PMID: 20573776 Zheng Y. (2011). Study on microorganism and bacteriostasis of leaf circumference of Pinus tabulaeformis. Jilin Agric. 10:65. Zhou Y. K. Shen X. Y. Hou C. L. (2017). Diversity and antimicrobial activity of culturable fungi from fishscale bamboo (phyllostachys heteroclada) in China. World J. Microbiol. Biotechnol. 33:104. doi: 10.1007/s11274-017-2267-9, PMID: 28466300 Zhu Y. G. Johnson T. A. Su J. Q. Qiao M. Guo G. X. Stedtfeld R. D. . (2013). Diverse and abundant antibiotic resistance genes in Chinese swine farms. Proc. Natl. Acad. Sci. U. S. A. 110, 3435–3440. doi: 10.1073/pnas.122274311, PMID: 23401528 Zhu D. Lu L. Zhang Z. J. Qi D. W. Zhang M. C. O’Connor P. . (2021). Insights into the roles of fungi and protist in the giant panda gut microbiome and antibiotic resistome. Enviro. Int. 155:106703. doi: 10.1016/j.envint.2021.106703, PMID: 34139588 Zhu L. Wu Q. Dai J. Zhang S. Wei F. (2011). Evidence of cellulose metabolism by the giant panda gut microbiome. Proc. Natl. Acad. Sci. U. S. A. 108, 17714–17719. doi: 10.1073/pnas.1017956108, PMID: 22006317