Leaf-Associated Epiphytic Fungi of Gingko biloba, Pinus bungeana and Sabina chinensis Exhibit Delicate Seasonal Variations

Plant-leaf surface on Earth harbors complex microbial communities that influence plant productivity and health. To gain a detailed understanding of the assembly and key drivers of leaf microbial communities, especially for leaf-associated fungi, we investigated leaf-associated fungal communities in two seasons for three plant species at two sites by high-throughput sequencing. The results reveal a strong impact of growing season and plant species on fungal community composition, exhibiting clear temporal patterns in abundance and diversity. For the deciduous tree Gingko biloba, the number of enriched genera in May was much higher than that in October. The number of enriched genera in the two evergreen trees Pinus bungeana and Sabina chinensis was slightly higher in October than in May. Among the genus-level biomarkers, the abundances of Alternaria, Cladosporium and Filobasidium were significantly higher in October than in May in the three tree species. Additionally, network correlations between the leaf-associated fungi of G. biloba were more complex in May than those in October, containing extra negative associations, which was more obvious than the network correlation changes of leaf-associated fungi of the two evergreen plant species. Overall, the fungal diversity and community composition varied significantly between different growing seasons and host plant species.


Introduction
The above-ground surface of plants, also known as the phyllosphere, covers an estimated area of 6.4 × 10 8 km 2 , which is approximately twice as large as the land surface area [1,2], providing a unique habitat for the colonization, interaction and evolution of microorganisms with plants. Studies have reported that there are more than 10 6 microbial cells/cm 2 of leaf surfaces [3]. These microbial communities include epiphytes that colonize the leaf surfaces and endophytes that occupy the interior spaces of the leaves [4,5], and have significant impacts on host plant growth, development and health [6]. There may be positive, neutral or negative interactions between microbes and plants [3,7]. Studies on phyllosphere microbiology have mainly focused on leaves, which are the plant structures

Sample Collection
The sampling location, sampling protocol and sample processing was described in detail by Bao et al. [27]. Briefly, the experimental site was chosen in the city of Beijing, China, including two sampling sites, a forest park (40 •  In order to minimize the influence of other factors, all the leaf samples were collected on sunny days under similar weather conditions, and the week before sample collection was also sunny. About 300 g of composite samples collected from each plant species were placed in a Labplas on ice, and quickly transported to the laboratory. In addition, all the leaf samples from the same plant were collected during two different sampling seasons (October and May). The outside temperature of sampling was as follows: 22 October 2017 (12 • C) and 10 May 2018 (27 • C).

DNA Extraction
The DNA extraction method for microbiomes in leaf samples was described in our previous study [27]. Briefly, for each composite sample, 30 g of leaves were randomly selected and placed into a 1 L sterile Erlenmeyer flask containing 500 mL of sterile PBS buffer (1 × phosphate-buffered saline buffer, pH 7.4). In order to obtain the microbial cells of leaves, the Erlenmeyer flask was first sonicated (frequency 40 kHz) in an ultrasonic cleaning bath for 6 min, then shaken at 200 r/min at 30 • C for 20 min and finally sonicated (frequency 40 kHz) for 3 min. To separate the microbial cells from the leaves, we filtered the cell suspension through a 0.22 µm sterile nylon membrane. Genomic DNA on each membrane could be extracted directly by the Fast ® DNA SPIN Kit (MP Biomedicals, Santa Ana, CA, USA), according to the manufacturer's instructions, after the membrane was minced with sterile scissors and stored at −80 • C.

Quantitative PCR Analysis
A primer set ITS1 (CTTGGTCATTTAGAGGAAGTAA)/ITS2 (GCTGCGTTCTTCATC-GATGC) was used to amplify the fungal internal transcribed spacer (ITS) region with the annealing temperature of 55 • C. Quantitative PCR assay was run by a CFX96 Optical Real-Time Detection System (Bio-Rad, Laboratories, Inc. Hercules, CA, USA). The total volume of each reaction was 20 µL, which contained 0.5 µM of each primer (10 mM), 1 µL of DNA template, 10.0 µL of SYBR ® Premix Ex Taq (Takara, Biotech, Dalian, China) and sterile double-distilled water (ddH 2 O). To obtain a standard curve of the target fragment, an ITS region fragment obtained by PCR was cloned into the pMD19-T vector (Takara, Biotech, Dalian, China) and then transformed into Escherichia coli JM109 competent cells. Plasmid DNA in the competent cells was extracted by the Plasmid Purification Kit (Takara, Biotech, Dalian, China). We selected and verified the plasmid DNA containing the correct fragment lengths. The selected plasmid DNAs were then used as templates for generating standard curves. Sterile water was run as the template, instead of the DNA, as blank controls. Amplification of the ITS region was mainly verified by melting curve analysis and agarose gel electrophoresis. Each DNA sample was assayed in a minimum of triplicate with R 2 values greater than 0.990 and amplification efficiencies ranging from 90% to 110%. The cycling condition was as follows: 40 cycles for 30 s at 95 • C, annealing for 30 s at 55 • C, extension at 72 • C for 30 s and a final extension at 72 • C for 8 min.

High-Throughput Sequencing
The extracted DNA was used as a template for the amplification of a 300-bp internal fragment of the ITS region, with a barcode primer set ITS1/ITS2. PCR analysis was performed in a reaction volume of 50 µL, containing 25 µL Premix Taq DNA polymerase, 1 µL DNA template (20 ng total DNA), 0.5 µL each forward/reverse primer (20 mM) and 23 µL ddH 2 O. Purified amplicons were pooled at equimolar concentrations and paired-end sequenced on an Illumina MiSeq PE250 platform (Majorbio, Shanghai, China), according to a standard protocol. The raw sequence data of leaf-associated fungi in October and May were deposited in the NCBI Sequence Read Archive (SRA) with the Accession Numbers PRJNA562843 and PRJNA562965, respectively. Raw reads from all fastq files were quality-filtered through Trimmomatic software (version 0.38) and barcodes were removed while assigning to the respective samples based on unique barcodes [28]. For each sample, paired-end reads were merged into full-length sequences using the FLASH software (version 1.2.11) [29]. These processed sequences were clustered by UPARSE (version 7.0.1090) [30] into operational taxonomic units (OTUs) with a sequence threshold of 97% similarity, while extracting representative sequences of the OTUs and filtering chimeras and singletons.

Statistical Analysis
Sequencing data were analyzed on an online platform using the Majorbio Cloud Platform at www.majorbio.com, accessed on 14 May 2019. The UNITE database (version 8.0) was used for the taxonomic identification of fungi. Based on the Bray-Curtis distance matrix, OTU-based community similarities were calculated by the vegan package of R software (version 3.3.1) using principal coordinate analysis (PCoA) and PerMANOVA analysis. Heatmap was created by R software (version 3.3.1) with vegan package and the scale represented the logarithm of the rarefied-reads number. Linear discriminant analysis (LDA) effect sizes were used to characterize the differences in fungal community compositions in different seasons, with an LDA score ≥ 2 indicating significant contributions (p < 0.05) to the model [31]. All the twelve-correlation network analyses were constructed at the genus level using the Maslov-Sneppen procedure [32,33], and the top 10 genera with Spearman's correlation coefficient (p < 0.01 and |SpearmanCoef| > 0.75) were visualized by Cytoscape (version 3.5.1). The functional modes of the leaf-associated fungal communities were predicted by FUNGuild (version 1.0) [34]. Origin (version 2016, OriginLab, Northampton, UK) was used to create the figures. One-way analysis of variance (ANOVA) of SPSS software (version 16, IBM, New York, NY, USA) was used to calculate the significant differences (p < 0.05) between different groups. The statistical significance level of all the analyses in this study was 0.05.

High Abundance and Diversity of Leaf-Associated Fungi
The abundances of fungi among the leaf samples of G. biloba, P. bungeana and S. chinensis from two locations over two growing seasons were detected by performing qPCR assays and were found to be different ( Figure 1). The results show that ITS region copies among the leaf samples range from 2.80 × 10 6 to 2.34 × 10 8 copies per gram of leaf. Seasonal variation was found to be the most important factor affecting the number of ITS region copies in different leaf samples of the three plant species. The abundance of leaf samples obtained in October 2017 was significantly higher than that of the leaf samples collected in May of the following year ( Figure 1A).
The fungal community structures among the leaf samples of G. biloba, P. bungeana and S. chinensis from two locations were assessed using high-throughput sequencing over two growing seasons. After sequences filtering and assembly, a total of 3,958,099 high-quality ITS reads were recorded from 60 detected leaf samples. Of these samples, the valid reads ranged from 41,275 to 84,327 for each sample and were rarefied to 41,275 for comparisons in leaf fungal communities in subsequent analyses (Table 1). Good's coverage of all samples was more than 99%, indicating that the sequencing depth included nearly all fungal communities in the detected leaf samples, which was sufficient to saturate fungal diversity. The results of fungal α-diversity (Shannon and Chao1 indices) in leaf samples of G. biloba show significant differences between two seasons, with the two indices in May being significantly higher than those in October. However, the difference in α-diversity between October and May was not significant in the leaf samples of P. bungeana and S. chinensis, except the Chao1 index of S-C in May ( Figure 1A). These results suggest that the leaf-associated fungal α-diversity of deciduous and evergreen plants may have different responses to seasonal changes. The PCoA analysis results of fungal β-diversity at the different locations and seasons reveal the obvious separation of the communities between two seasons in different samples ( Figure 1B). PerMANOVA analysis further showed that changes in seasons, plant species and locations all had significant effects on the leaf-associated fungal community structures (p = 0.001). Moreover, seasons (R 2 = 0.294) The fungal community structures among the leaf samples of G. biloba, P. bungeana and S. chinensis from two locations were assessed using high-throughput sequencing over two growing seasons. After sequences filtering and assembly, a total of 3,958,099 highquality ITS reads were recorded from 60 detected leaf samples. Of these samples, the valid reads ranged from 41,275 to 84,327 for each sample and were rarefied to 41,275 for comparisons in leaf fungal communities in subsequent analyses (Table 1). Good's coverage of all samples was more than 99%, indicating that the sequencing depth included nearly all fungal communities in the detected leaf samples, which was sufficient to saturate fungal diversity. The results of fungal α-diversity (Shannon and Chao1 indices) in leaf samples of G. biloba show significant differences between two seasons, with the two indices in May being significantly higher than those in October. However, the difference in α-diversity between October and May was not significant in the leaf samples of P. bungeana and S. chinensis, except the Chao1 index of S-C in May ( Figure 1A). These results suggest that the leaf-associated fungal α-diversity of deciduous and evergreen plants may have different responses to seasonal changes. The PCoA analysis results of fungal β-diversity at the different locations and seasons reveal the obvious separation of the communities between two seasons in different samples ( Figure 1B). PerMANOVA analysis further showed that changes in seasons, plant species and locations all had significant effects on the leaf-associated fungal community structures (p = 0.001). Moreover, seasons (R 2 = 0.294) and tree species (R 2 = 0.284) had a greater effect on the leaf-associated fungal communities than locations (R 2 = 0.086).

Variations in Community Composition of Leaf-Associated Fungi during October and May
The rarefied dataset represented fungal communities from 8 phyla, 35 classes, 103 orders, 255 families, 531 genera and 2401 OTUs. The dominant phyla of all the detected leaf samples were Ascomycota and Basidiomycota ( Figure A1). The relative abundance of the leaf-associated fungal community compositions at class level for the three plant species at two sites over two different seasons is shown in Figure 2A by performing a community bar plot analysis. Dothideomycetes was the most abundant class among all the detected leaf samples, except for S-F in May with a higher abundance of Cystobasidiomycetes than Dothideomycetes. Among these samples, S-C in October had the highest abundance of Sordariomycetes, S-F in October had the highest abundance of Eurotiomycetes, Taphrinomycetes and Ustilaginomycetes, and P-F in October had the highest abundance of Microbotryomycetes. Heatmap data displayed in Figure 2B show that the relative abundances of fungal community compositions among the three tree species over different seasons and different locations greatly vary at the genus level. The fungal community biomarkers of the three plant species in two different seasons were studied in detail by LDA analysis, and the graph shows the genus-level fungal biomarkers with LDA scores >2.5, p < 0.05 for each leaf sample ( Figure 3). These results further demonstrate that the compositions of leaf-associated fungi at the genus level in October can be clearly distinguished from those in May among the three plant species. In detail, for deciduous tree G. biloba, the number of enriched fungal genera was clearly greater in May than that in October. For evergreen trees P. bungeana and S. chinensis, the number of enriched genera was higher in October than that in May, especially the samples obtained from the forest. For the three plant species, Alternaria, Cladosporium and Filobasidium were mostly enriched in October compared to May.   Variations in fungal communities may be closely related to changes in functions [35,36], so revealing the modes of functional transitions is important for understanding the ecological processes of leaf-associated fungal communities in different seasons or locations. The relative abundances of specific trophic modes showed that the functional modes of the same plant in different seasons or locations were different, especially the distinct shift in the functional modes of S. chinensis (Figure 4). It also revealed that different plant species had fungal communities with different functional modes. The results demonstrate that 93.60% of the total OTUs of leaf-associated fungi are classified into different trophic modes. Among them, saprotrophs, pathotrophs and symbiotrophs accounted for 16.25%, 6.86% and 0.51%, respectively. The remaining 69.98% of trophic modes were classified into the multiple trophic modes. For G. biloba, the relative abundances of fungal pathotrophs on campus decreased from 6.53% in May to 0.50% in October, and the relative abundances of fungal pathotrophs in the forest decreased from 4.85% in May to 0.81% in October. For P. bungeana, the relative abundances of fungal pathotrophs on campus tended to decrease from 0.71% in May to 0.22% in October, and the relative abundances of fungal pathotrophs in the forest decreased from 6.00% in May to 4.04% in October, with little difference between the different seasons. For S. chinensis, the relative abundances of fungal pathotrophs on campus increased from 3.05% in May to 7.12% in October, and the relative abundances of fungal pathotrophs in the forest increased from 9.15% in May to 39.31% in October.   pathotrophs on campus tended to decrease from 0.71% in May to 0.22% in October, and the relative abundances of fungal pathotrophs in the forest decreased from 6.00% in May to 4.04% in October, with little difference between the different seasons. For S. chinensis, the relative abundances of fungal pathotrophs on campus increased from 3.05% in May to 7.12% in October, and the relative abundances of fungal pathotrophs in the forest increased from 9.15% in May to 39.31% in October.

Variations in Fungal Interactions of Leaf-Associated Fungi during October and May
In order to compare the differences in the interactions between the leaf-associated fungal communities of the three plant species in May and October at campus sampling sites, six correlation networks were constructed based on the genus level ITS dataset with significant (p < 0.01) and robust (q > 0.75 or q < −0.75) correlations (Figure 5). For G. biloba leaf-associated fungal communities, the correlation networks in October and May consisted of 123 pairs with 82 nodes, and 192 pairs with 123 nodes, respectively. For P. bungeana leaf-associated fungal communities, the correlation networks in October and May consisted of 84 pairs with 38 nodes, and 97 pairs with 53 nodes, respectively. For S. chinensis leaf-associated fungal communities, the correlation networks in October and May consisted of 125 pairs with 80 nodes, and 117 pairs with 85 nodes, respectively. Among the tested three plants, the changes in fungal network structures of deciduous tree G. biloba in different seasons were greater than those of evergreen trees P. bungeana and S. chinensis.

Variations in Fungal Interactions of Leaf-Associated Fungi during October and May
In order to compare the differences in the interactions between the leaf-associated fungal communities of the three plant species in May and October at campus sampling sites, six correlation networks were constructed based on the genus level ITS dataset with significant (p < 0.01) and robust (q > 0.75 or q < −0.75) correlations ( Figure 5). For G. biloba leaf-associated fungal communities, the correlation networks in October and May consisted of 123 pairs with 82 nodes, and 192 pairs with 123 nodes, respectively. For P. bungeana leaf-associated fungal communities, the correlation networks in October and May consisted of 84 pairs with 38 nodes, and 97 pairs with 53 nodes, respectively. For S. chinensis leaf-associated fungal communities, the correlation networks in October and May consisted of 125 pairs with 80 nodes, and 117 pairs with 85 nodes, respectively. Among the tested three plants, the changes in fungal network structures of deciduous tree G. biloba in different seasons were greater than those of evergreen trees P. bungeana and S. chinensis.
Six correlation networks with significant (p < 0.01) and robust (q > 0.75 or q < −0.75) correlations were constructed to understand the differences in the interactions of the leaf-associated fungi in different seasons on forest sampling sites ( Figure 6). The leaf fungal networks of G. biloba consisted of 64 pairs with 47 nodes, and 180 pairs with 107 nodes in October and May, respectively. The leaf fungal networks of P. bungeana consisted of 98 pairs with 57 nodes, and 77 pairs with 63 nodes in October and May, respectively. The leaf fungal networks of S. chinensis consisted of 111 pairs with 85 nodes, and 104 pairs with 64 nodes in October and May, respectively. Among the tested three plants, the fungal network structures of G. biloba in May were more complicated than that in October both at two sampling sites, especially in forest sampling sites. These results indicate that seasonal factors have a great influence on the interactions between the fungal community compositions of deciduous tree, and the influence is greater than that on evergreen trees P. bungeana and S. chinensis. Moreover, the number of negative interactions between the G. biloba leaf-associated fungi in May were much higher than that in October. It is inferred that competitive interaction might be one of the reasons why the total abundance of fungi in May was lower than that in October.  . Fungal ITS region genera-based correlation network for G. biloba, P. bungeana and S. chinensis on campus in October and May. A node represents a genus. A connection stands for a strong (Spearman's q > 0.75 or q < −0.75) and significant (p < 0.01) correlation. Edge widths are scaled according to their weights and edge colors indicate a positive (blue) or negative (red) correlation for the nodes they connect. The numbers on columns indicate the ratio of increase in May compared to October. G-C, G. biloba from campus; P-C, P. bungeana from campus; S-C, S. chinensis from campus. Figure 6. Fungal ITS region genera-based correlation network for G. biloba, P. bungeana and S. chinensis at forest in October and May. A node represents a genus. A connection stands for a strong (Spear-man´s q > 0.75 or q < −0.75) and significant (p < 0.01) correlation. Edge widths are scaled according to their weights and edge colors indicate a positive (blue) or negative (red) correlation for the nodes they connect. The numbers on columns indicate the ratio of increase in May compared to October. G-F, G. biloba from forest; P-F, P. bungeana from forest; S-F, S. chinensis from forest.

Discussion
Diverse microorganisms, including fungi, live on plant leaves, which can affect plant growth, development and even evolution [3,6]. In spite of a growing recognition of the importance of these microorganisms to plants and even ecosystems and their potential for crop improvement, the factors that influence phyllosphere microbial communities and their functional consequences require more research to thoroughly explore them. As a step towards a more comprehensive understanding of the dynamics of leaf-associated fungal community in different plant species, we investigated the leaf-associated fungi of three plant species during two growing seasons at different locations. Our study demonstrated clear differences in leaf-associated fungal communities between different growing seasons and plant species in the northern region, whereas small geographic distance seemed to have little effect on the leaf-associated fungal communities.
Ascomycetes (Dothideomycetes, Eurotiomycetes and Sordariomycetes) and Basidiomycetes (Cystobasidiomycetes, Tremellomycetes and Microbotryomycetes) dominated the epiphytic fungal communities in the studied plant species, in line with the results of previous studies [37,38]. Leaf-associated fungal community structure varied strongly throughout the different growing seasons; the core fungal taxa, especially, showed distinct seasonal abundance patterns at higher taxonomic levels. Katsoula et al. [39] observed Figure 6. Fungal ITS region genera-based correlation network for G. biloba, P. bungeana and S. chinensis at forest in October and May. A node represents a genus. A connection stands for a strong (Spearman's q > 0.75 or q < −0.75) and significant (p < 0.01) correlation. Edge widths are scaled according to their weights and edge colors indicate a positive (blue) or negative (red) correlation for the nodes they connect. The numbers on columns indicate the ratio of increase in May compared to October. G-F, G. biloba from forest; P-F, P. bungeana from forest; S-F, S. chinensis from forest.

Discussion
Diverse microorganisms, including fungi, live on plant leaves, which can affect plant growth, development and even evolution [3,6]. In spite of a growing recognition of the importance of these microorganisms to plants and even ecosystems and their potential for crop improvement, the factors that influence phyllosphere microbial communities and their functional consequences require more research to thoroughly explore them. As a step towards a more comprehensive understanding of the dynamics of leaf-associated fungal community in different plant species, we investigated the leaf-associated fungi of three plant species during two growing seasons at different locations. Our study demonstrated clear differences in leaf-associated fungal communities between different growing seasons and plant species in the northern region, whereas small geographic distance seemed to have little effect on the leaf-associated fungal communities.
Ascomycetes (Dothideomycetes, Eurotiomycetes and Sordariomycetes) and Basidiomycetes (Cystobasidiomycetes, Tremellomycetes and Microbotryomycetes) dominated the epiphytic fungal communities in the studied plant species, in line with the results of previous studies [37,38]. Leaf-associated fungal community structure varied strongly throughout the different growing seasons; the core fungal taxa, especially, showed distinct seasonal abundance patterns at higher taxonomic levels. Katsoula et al. [39] observed a strong seasonal effect on the composition of the leaf epiphytic fungal communities in a semi-arid Mediterranean ecosystem (summer vs. winter). The abundance of the total epiphytic fungal community decreased significantly from autumn to spring, which was contrary to the findings of olive trees in Mediterranean ecosystems [40]. Members of Dothideomycetes (Ascomycota) were at their highest abundance among all the detected leaf samples in the growing seasons of October and May. Dothideomycetes is the largest and most ecologically and functionally diverse class of fungi, containing human and plant pathogens, endophytes and epiphytes [41,42], which has been widely reported as a dominant leaf-associated taxa globally [38,43,44]. The relative abundance of Cystobasidiomycetes (Basidiomycota) in May was much higher than that in October, only in evergreen trees P. bungeana and S. chinensis. At the genus level, the core fungal microorganisms in this study included Phoma, Cladosporium, Epicoccum and Alternaria, which were found to be dominant in previous studies of other plant leaf-associated fungal communities [45][46][47][48]. The previous study using culture-dependent methods found that Alternaria, Cladosporium, Fusraium and Phoma were the main epiphytic fungi on Hornbeams [49], which were both included in our study with culture-independent methods. Furthermore, studies have shown that Epicoccum is an important endophytic fungus involved in the production of secondary metabolites, and some species, such as E. nigrum, can act as biocontrol agents to promote plant growth and resist brown-rot pathogen [50,51]. Therefore, the presence of this genus is of potential interest for future studies. Interestingly, the enrichment of Filobasidium in October might be related to a novel radiation resistance of Filobasidium sp. [52].
In addition to seasonal effects, host plant species have been observed to explain shifts in leaf-associated microbial communities [53,54], which follows the variations in fungal communities among different plant species in our study. The leaf-associated fungal abundance of Ginkgo biloba was significantly higher in May than in October, which may be related to the different maturation of leaves. Studies have shown that significant differences in epiphytic fungal richness were observed among olives obtained in different production systems and maturation stages [55]; whereas, in the same production system, there was no statistical significance in richness regarding olive cultivars, and this result is consistent with the reports of epiphytic fungal communities in olives [56] and mango fruits [57] from different cultivars studied using culture-dependent methods. The screening role of host plant species on the epiphytic fungal communities [38,43] has been fully confirmed and attributed to different ecological strategies as well as the chemical and functional properties of host plants. Different plant species have several different plant traits, including average plant height, leaf nitrogen concentration and plant genotypes. Previous studies have shown that plant height is closely related to changes in leaf-associated bacterial communities [58,59]. Studies of leaf-associated microbiomes in a variety of plant types showed that variations in the leaf-associated fungal communities were significantly related to the changes in leaf mass and nitrogen concentration per leaf area [25,37]. Genotypes play a critical effect in determining the colonization and establishment of leaf-associated microbial communities in plant species [43,60]. Many phenotypic properties derived from the plant genetic repertoire, such as leaf morphology, physiology and chemistry, may affect the interaction of plants with herbivores, pathogens, competitors and symbionts, including the leaf-associated microorganisms [61]. Leaf fungal communities may be subject to selective pressures derived from specific dynamics arising from plant phenotypes, resulting in the adaptation of local communities to individual host genotypes [62]. A previous study that sampled 32 Arabidopsis thaliana individuals and collected 21 air samples over a 73 day period confirmed that leaf microbial community compositions initially reflected the airborne community, but then formed a unique community, revealing great host selectively for community members [63]. A plausible mechanism for the effect of genotype on leafassociated microbial communities may be related to leaf-surface properties, such as leaf wax [53], ethylene perception [53], gamma-aminobutyrate pathway [64], or salicylic acid and jasmonic acid signaling defense pathways [65]. Genome-wide association studies revealed that the variation in the leaf-associated microbial community was significantly associated with host plant loci responsible for defense and cell wall integrity [66]. It can be seen that the mechanisms of deciphering the associations of host plants to form microbial communities is highly complex, and more extensive and in-depth research is required.
Epiphytes inhabiting the phyllosphere are under a set of selective pressures, which are distinct from those that endophytes are facing [67]. The epiphytic community was found to be significantly richer and more abundant than the endophytic community in some woody plant systems [40,68,69]. There was a great difference in the effects of seasons on the total abundance of phyllosphere epiphytic fungi and endophytic fungi. Gomes et al. [40] found that seasons had a significant effect on the total abundance of epiphytic fungi, but had no significant effect on endophytic fungi. It is significant to compare the phyllosphere epiphytic and endophytic fungal communities in the future studies.
Functional mode analysis found that, among all the samples, only Sabina chinensis at the forest sampling site in October had the highest attributes of leaf fungal pathotrophs, indicating that the host plants may be infected by pathogenic fungi, which requires us to pay more attention to the disease and health of S. chinensis. The network among leaf-associated fungi of the deciduous tree G. biloba in May was more complex than that in October, showing more negative correlations, which was consistent with the previous studies on bacteria, but not with changes in the fungal network of the evergreen trees P. bungeana and S. chinensis. The unstable state of competition for space and nutrients in leaf-associated fungi of the deciduous tree G. biloba in May could have resulted in lower overall fungal numbers than in October. Positive correlations in fungal-fungal network correlations were more common than negative correlations in leaf-associated core taxa associations in the microbial cooccurrence networks of the evergreen trees P. bungeana and S. chinensis, in agreement with the recent reports [25]. Among the core taxa associations in the leaf-associated microbial co-occurrence networks, fungal-bacterial associations were more likely to be negative than positive [25,70,71]. Studies of Arabidopsis root microorganisms have shown that these negative correlations reflect direct or indirect bacterial biocontrol of detrimental fungi, which have significant implications for plant survival and health [71]. In future studies, we need to explore more comprehensive microbial network correlations, including fungalfungal network correlations, bacterial-bacterial network correlations and fungal-bacterial network correlations, to protect plants against fungi and oomycetes through biological control methods.

Conclusions
We extended the knowledge about the leaf-associated epiphytic fungal microbiomes of three municipal greening plants, Pinus bungeana, Sabina chinensis and Gingko biloba, in northern cities of China. The leaf-associated epiphytic fungal communities from two different sites over two seasons were compared, revealing the strong impact of growing seasons and plant species on fungal community composition, exhibiting clear temporal patterns in abundance and diversity. Further investigations will elucidate the mechanism of this effect. Among the genus-level biomarkers, the abundances of Alternaria, Cladosporium and Filobasidium were significantly higher in October than in May in the three tree species. Network correlations between the leaf-associated fungi of G. biloba were more complex in May than in October, containing extra-negative associations, which was more obvious than the network correlation changes in leaf-associated fungi of the two evergreen plant species. Further studies will be conducted to investigate biomarkers as indicators of total microbial network correlations for protecting plant health.