Plant leaves are complex microecosystems abounding with microbial life (Vorholt 2012). Globally, total leaf area is estimated at 2 × 10 14 m 2 (Bar-On and Milo 2019), providing a ubiquitous habitat for microbial interaction, colonization, and evolution with plants. Indeed, studies report upward of 10 6 microbial cells/cm 2 of leaf surface (Lindow and Brandl 2003). These microbial communities, termed the phyllosphere microbiome, inhabit either the surface (epiphytes) or interior (endophytes) of the leaf (Carvalho and Castillo 2018) and can have important effects on plant health (Stone et al. 2018). Phyllosphere research has traditionally focused on agriculturally relevant pathogens; however, a growing number of studies indicate beneficial effects of the phyllosphere microbiome such as in pathogen (Ritpitakphong et al. 2016) or herbivore defense (Saleem et al. 2017). Therefore, management of phyllosphere communities provides numerous opportunities for improving crop productivity.

Plant-fungus interactions are an important contributor to ecosystem-level processes such as carbon and nutrient cycling (Tedersoo and Bahram 2019;Zhu and Miller 2003). Most studies to date have focused on fungal communities in the root-soil interface (the rhizosphere), given the importance of rhizosphere fungi in pathogenesis, decomposition, and plant nutrient uptake (Bonito et al. 2019;Deveau 2016;Shakya et al. 2013). Yet phyllosphere RESEARCH fungal communities are often extremely diverse (Izuno et al. 2016;Kirschner 2018) and can dictate plant productivity. For example, fungal foliar pathogens are responsible for a variety of devastating plant diseases such as leaf rusts, blights, and powdery or downy mildews in a variety of crop species (Jain et al. 2019). Additionally, common (nonpathogenic) phyllosphere fungi can either increase or decrease the severity of pathogenesis, as shown for Melampsora leaf rust disease in the bioenergy crop Populus (Busby et al. 2016). Finally, other phyllosphere fungi are hypothesized to be dormant saprotrophs, which play key roles in decomposition and nutrient cycling following leaf senescence (Osono 2006). Nevertheless, controlled microbiota management is a major challenge due, at least in part, to the wide variety of factors that can affect phyllosphere fungal community structure and assembly. For example, seasonality (Breeze and Dix 1981), host genotype (Qian et al. 2018), environmental gradients such as elevation and temperature (Cordier et al. 2012), location in the canopy (Izuno et al. 2016), and location in versus on the leaf (epiphytes versus endophytes, respectively) (Yao et al. 2019) can all influence fungal community structure in the phyllosphere. Additional work exploring phyllosphere fungal community dynamics is needed to improve understanding of phyllosphere ecology and function (Kirschner 2018).

Switchgrass (Panicum virgatum L.) is a warm-season perennial bunchgrass found throughout North America. The species has received considerable attention as a model biofuel crop, given its rapid growth rate, wide geographic distribution, and minimal management requirements (Wright and Turhollow 2010). Phyllosphere fungi play a major role in switchgrass ecology and agriculture. In particular, pathogens can significantly reduce yields (Gravert et al. 2000;Krupinsky et al. 2004). Recent culturedependent studies have reported that switchgrass fungal phyllosphere communities are environmentally structured (Whitaker et al. 2018) and are highly diverse, with numerous taxa that can enhance plant height and biomass in experimental manipulations (Xia et al. 2018). Given the potential for microbiome management to improve switchgrass productivity and stress resistance, further foundational work exploring fungal phyllosphere dynamics in switchgrass is needed. In addition, an important consideration of microbiome analysis is that phyllosphere fungi do not exist in isolation: recent reports on the abundance and diversity of phyllosphere bacterial communities in switchgrass (Gagne-Bourgue et al. 2013;Grady et al. 2019) suggest the potential for many thousands of complex interactions that can influence phyllosphere community composition and function. Although numerous studies have begun exploring intra-and interdomain phyllosphere co-occurrence in a variety of plant species and environments (Agler et al. 2016), our understanding of leaf microbial co-occurrences and their potential ecological relevance is still limited (Hassani et al. 2018).

In this study, we investigate seasonal dynamics of phyllosphere epiphytic fungal communities in switchgrass. We complement and extend a recent study of the switchgrass phyllosphere bacterial and archaeal microbiome (Grady et al. 2019) with an investigation of fungal community seasonal dynamics using a subset of the same switchgrass leaf samples. We also leverage the published bacterial community dataset of Grady et al. (2019) to explore intra-and interdomain network associations in the switchgrass phyllosphere. That study found consistent seasonal dynamics in bacterial communities, with similar groups of early-, mid-, and late-season bacterial taxa across two growing seasons (Grady et al. 2019). We hypothesized that fungal communities would exhibit strong seasonal shifts, similar to that observed for bacterial communities in the switchgrass phyllosphere, resulting in a hyperconnected fungal-bacterial correlation network.

Experimental design, sample collection, and processing. The field site, sampling scheme, and sample processing are described in detail by Grady et al. (2019). In that study, switchgrass phyllosphere and soil bacterial and archaeal communities were assessed using 16S V4 amplicon sequencing over two growing seasons. Here, we complemented that study by using a subset of those same samples from the 2017 growing season for fungal internal transcribed spacer 2 (ITS2) amplicon sequencing and community analysis.

Briefly, switchgrass leaf and bulk soil samples were collected from the switchgrass monoculture plots (G5) of the Great Lakes Bioenergy Research Center Biofuel Cropping System Experiment in Hickory Corners, MI (42°23941.60 N, 85°22923.10 W), where switchgrass has been grown continuously since 2008. Switchgrass (P. virgatum L. var. Cave-in-Rock) is grown in replicate 30-by-40m plots: the westernmost 4.5 m of each plot is maintained with no fertilizer (nitrogen free), while the remainder of the plot receives nitrogen at 50 lb. acre _ 1 in a single yearly application in May. Every 3 weeks across the 2017 growing season (April through September), we sampled from plots 1 to 4 in both the main and nitrogenfree subplots. Leaf samples were collected at each sampling date starting from 15 May 2017 (after switchgrass emergence). Leaves were collected and pooled at three flags along a standardized path in each plot. Leaves were collected using ethanol-sterilized gloves and stored in sterile Whirl-Pak bags. We also collected two bulk soil cores (10 cm in depth by 2 cm in diameter) at each of the three flags in each subplot, sieved them through 4-mm mesh, then pooled the three replicates per subplot for soil chemistry analysis. Samples were transported on wet ice before being stored at _ 80°C. At each sampling date, plant height and soil temperature were measured on a per-plot basis, while weather data were collected from the Michigan State University Weather Station Network for Kellogg Biological Station (https://mawn.geo.msu.edu).

Leaf dry matter content was assessed using 10 leaves per plot following (Cornelissen et al. 2003). The dried leaves were then ground to a fine powder using a Sampletek 200 vial rotator with iron roll bars (Mavco Industries, Lincoln, NE, U.S.A.). Leaf carbon and nitrogen were measured on a Costech ECS 4010 elemental analyzer (Costech Analytical Technologies Inc., Valencia, CA, U.S.A.). Soil pH, lime, P, K, Ca, Mg, organic matter, NO 3 , NH 4 , and percent moisture were measured by the Michigan State University Soil and Plant Nutrient Lab (http://www.spnl.msu.edu/). Phyllosphere surface microbial DNA was extracted from 5 to 10 leaves per subplot (approximately 5 g of leaf tissue) using a benzyl chloride liquidliquid extraction, followed by isopropanol precipitation, according to the methods described by Suda et al. (2008) and also applied by Grady et al. (2019). Following extractions, DNA concentrations were quantified using a Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, U.S.A.) and normalized. Paired-end amplicon sequencing was completed by the Department of Energy's Joint Genome Institute (JGI) on an Illumina MiSeq, using ITS9F (GAACGCAGCRAAIIGYGA) and ITS4R (TCCTCCGCTTATTGATATGC) primers to target the ITS2 ribosomal DNA (rDNA) region (Ihrmark et al. 2012;White et al. 1990). In total, 52 leaf samples (minimum of 4 and maximum of 8 per sampling date) were collected across the growing season (see Supplementary Table S1 for sampling effort and distribution across sampling dates, plots, and fertilization treatments). The DNA extractions used for ITS2 amplicon sequencing were aliquots of the same DNA extractions used for 16S ribosomal RNA gene V4 amplicon sequencing in the study by Grady et al. (2019).

Sequence processing and statistical analysis. Raw forward and reverse Illumina ITS reads were first quality evaluated with FastQC (Andrews 2010). Reads were then merged using USEARCH v11 (Edgar 2010), Illumina adapters and sequencing primers were removed using Cutadapt (Martin 2011), and quality filtered and trimmed to equal length (180 bp) in USEARCH v11 (Edgar 2016;Edgar and Flyvbjerg 2015). Highly conserved 5.8S (downstream) and 28S (upstream) regions of the ITS2 were removed to avoid biasing operational taxonomic unit (OTU) clustering and taxonomy assignments based on sequence similarity of conserved rDNA regions; then, reads were trimmed to an equal length of 180 bp. Sequences were dereplicated and singletons were removed, then clustered into OTUs, based on 97% similarity, using the UPARSE (Edgar 2013) algorithms. Taxonomy assignment of OTU representative sequences was performed using CONSTAX (Gdanetz et al. 2017) with the Eukaryotic UNITE+INSD database (version 8.2, released 2 April 2020) (Abarenkov et al. 2020). To filter nonfungal sequences as well as OTUs unidentified at the kingdom level out of the dataset, we ran CONSTAX twice. First, we used 0.6 confidence to assign sequences at the kingdom level, then filtered nonfungal and nonidentified sequences from the dataset. We then reclassified the fungal OTUs at 0.8 confidence to clean the taxonomic hierarchy at lower ranks.

All statistical analyses were performed in R version 3.6.1 (R Core Team 2019). Package RSQLite v2.1.2 (Muller et al. 2019) was used for data retrieval, and packages stringr v1.4.0 (Wickham 2019), tibble v2.1.3 (Muller and Wickham 2019), tidyr v0.8.3 (Wickham and Henry 2019), dplyr v0.8.3 (Wickham et al. 2019), and purrr v0.3.2 (Henry and Wickham 2019) were used for data manipulation, while ggplot2 v3.2.0 (Wickham 2016), gridExtra v2.3 (Baptiste 2017), scales v1.0.0 (Wickham 2018), and cowplot v1.0.0 (Wilke 2019) were used for data visualization. First, samples were rarefied to the minimum number of sequences observed per sample (n = 48,088) using vegan v2.5-5 (Oksanen et al. 2019). Fungal species accumulation curves were calculated as the cumulative number of unique fungal species detected across the growing season. Next, we performed a principal coordinates analysis (PCoA) based on Bray-Curtis dissimilarities of fungal phyllosphere communities. We assessed the effects of environmental and leaf characteristics on fungal community composition by fitting variables to the PCoA scores using envfit() in vegan. Parameters that had significant explanatory value (P < 0.05) were included as vectors in the ordinations. In addition, we tested for the effects of sampling date and fertilization and their interaction on community composition using permutational multivariate analysis of variance (Anderson 2001) implemented in vegan with 999 iterations.)

We generated abundance-occupancy plots, which illustrate the proportion of samples that a given taxon occurs in versus the relative abundance of that taxon (Shade and Stopnisek 2019). We designated a prioritized "core phyllosphere taxa" as those OTUs that were detected in all of the samples of at least one sampling date (n = 8 for all sampling dates, including four plots containing each two subplots of N+/N _ , with the exception of the September sampling date, where n = 4). We pooled the N+/N _ subplots to focus our next analyses on understanding seasonal dynamics. To assess whether groups of core taxa exhibit similar seasonal dynamics, we performed a hierarchical cluster analysis. Z-scored seasonal dynamics in relative abundance for each core OTU were used to build a complete linkage distance matrix. We then extracted k = 7 clusters of taxa to correspond to the 7 days of sampling and plotted shifts in relative abundance of each cluster over time to visualize their seasonal dynamics. We also visualized seasonal dynamics in core taxa abundances aggregated by fungal phylum, class, and genus. We calculated a diversity (species richness, Shannon diversity, and Pielou's Evenness), and also tested for differences in group dispersions with PERMDISP2 (Anderson 2006) implemented in vegan. We then conducted an analysis of variance (ANOVA) for each of these a diversity and dispersion measures. Pielou's evenness was square-root transformed and dispersions were log 10 transformed to meet ANOVA assumptions. The homogeneity of variance assumption was assessed using Bartlett's test (Snedecor and Cochran 1989) and Levene's test (Levene 1960) in the car package v3.0-3 (Fox and Weisberg 2019), and normality of residuals was assessed using the Shapiro-Wilk test (Shapiro and Wilk 1965).

Next, a switchgrass phyllosphere core bacterial OTU table was recreated from the study of Grady et al. (2019) using their publicly available data and analysis pipeline and filtered for the samples shared with the phyllosphere ITS dataset. The sequencing depth of that study was 1,000 reads per sample and was exhaustive to capture bacterial phyllosphere diversity, as shown in the rarefaction curve and supporting materials comparing rarefaction depths by Grady et al. (2019). We then conducted a network analysis on the core bacterial and fungal taxa using the Molecular Ecological Network Analysis Pipeline (Deng et al. 2012;Zhou et al. 2010Zhou et al. , 2011)). The majority parameter was set to 12, which indicates that an OTU must appear in at least 12 samples to be included in the network analysis. We note that core taxa were defined as taxa occurring in all samples of a given sampling date, and each sampling date had up to eight samples. Therefore, the majority parameter of 12 guaranteed that only core OTUs detected in multiple sampling dates were included in the network, allowing us to detect correlated changes in seasonal dynamics between taxa. For the network analysis, missing data were kept blank and rarefied read numbers were log transformed. Correlations were assessed using the Pearson correlation coefficient with a random matrix theory threshold of 0.800 (indicating P < 0.05). Calculated network properties were compared with 100 random networks generated using the same number of nodes and edges. Cytoscape version 3.7.1 was used for network visualization (Shannon et al. 2003).

Data availability and computing workflows. Raw reads are available on the JGI Genome Portal under Project ID 1191516 https:// genome.jgi.doe.gov/portal/SwiandphyiTagsII_FD/Swiandphyi-TagsII_FD.download.html. All plant and environmental data, as well as computational workflows and custom scripts, are available on GitHub (https://github.com/ShadeLab/PAPER_Bowsher_GLBRC_ITS_ 2020_Phytobiomes).

Sequencing summary. Sequencing efforts resulted in 54.0 million reads for 52 leaf samples (Supplementary Table S2). Subsequent leaf sample sequence processing resulted in an average of 245,007 reads per sample, with a minimum of 48,088 and a maximum of 690,267 after removing singletons. Rarefying the dataset to the minimum number of reads (48,088) resulted in an average of 344 fungal OTUs per sample (minimum of 179, maximum of 468) (Fig. 1).

Fungal community composition is associated with numerous plant and environmental factors. The first two principal coordinate axes captured 48.7% of the variation in the dataset (Fig. 2). Fungal communities were significantly affected by both sampling date (F = 15.37, P = 0.001) and fertilization (F = 2.32, P = 0.024), and their interaction (F = 2.22, P = 0.031). Communities were strongly distinguished by date of sampling (P = 0.002), as indicated by the separation between sampling week centroids and the clear trajectory from weeks one through seven. Variation in phyllosphere fungal communities was also associated with variation in minimum air temperature (P = 0.047), as well as leaf dry matter content, leaf nitrogen content, leaf carbon/nitrogen ratio, and plant height (all P = 0.001) (Fig. 2; Supplementary Table S3).

Seasonal effects on phyllosphere fungal communities. The cumulative number of unique phyllosphere fungal taxa sharply increased between the first and second sampling dates and continued to increase throughout the growing season, albeit at a progressively slower rate (Fig. 3A). Fungal species richness (number of taxa detected at each sampling date) increased after the first sampling date and remained high until late summer, when it decreased across the final two sampling dates (Fig. 3B). Shannon diversity and Pielou's evenness exhibited similar patterns, with the highest values found at the May sampling date (Fig. 3C andD). Dispersion of phyllosphere communities (i.e., within-group community variance) was high on the first sampling date, then was low and stable through the end of the growing season (Supplementary Fig. S1).

Dynamics of core phyllosphere fungi and interdomain networks with bacterial taxa. Across the 52 phyllosphere samples, we detected 227 core fungal OTUs, which we defined as OTUs which were present in all of the samples (subplots) on at least one sampling date (Fig. 4A). These OTUs ranged in occupancy from 15.4 to 100% (mean of 67.1%) across all phyllosphere samples. We next used hierarchical clustering to group core OTUs into clusters distinguished by changes in relative abundance across the growing season (Supplementary Fig. S2), and plotted shifts in relative abundance of OTUs in each cluster over time to visualize differences between clusters in their broad seasonal dynamics (Fig. 4B). We found several distinct patterns in seasonal dynamics (Fig. 4B). Two groups of OTUs exhibited the highest abundance in the early season, then steadily decreased throughout the remainder of the growing season (Fig. 4B, clusters 3 and 5). One group of OTUs generally increased throughout the growing season, followed by a decrease at the final two sampling dates (cluster 1), while others were at minimal abundance throughout the growing season followed by an increase at the end of the season (clusters 4 and 6). Still others exhibited the highest abundance in midseason (clusters 2 and 7). Seasonal dynamics differed between phylogenetic groups: Aureobasidium and Neoascochyta spp. peaked in relative abundance in early summer, then decreased in abundance across the growing season, while Epicoccum spp. were at high abundance early in the growing season yet steadily increased throughout much of the growing season (Fig. 4C). On the other hand, the smut fungus Tilletia peaked in relative abundance midsummer but was characterized by a very low relative abundance both early and late in the season (Fig. 4C). Taxa that were unidentified at the genus level (82 of the 227 core OTUs) increased dramatically at the final two sampling dates (Fig. 4C). Contrasting seasonal dynamics across phylogenetic groups was also apparent at higher taxonomic levels. For instance, class Dothideomycetes (Ascomycota) was at its highest abundance early and late in the growing season, whereas members of class Exobasidiomycetes (Basidiomycota) were in highest abundance midseason (Supplementary Fig. S3), mirroring the pattern seen in overall a diversity (Fig. 3B).

Considering the 43 phyllosphere samples shared between the current work and that of Grady et al. (2019), we constructed a fungal-bacterial co-occurrence network. The network's node degree distribution fit to the power law model (R 2 = 0.903), and differed from a random network in modularity, connectedness, and average clustering coefficient (Supplementary Table S4). The network analysis resulted in 147 nodes (124 fungi out of 227 fungal core OTUs and 23 bacteria out of 41 bacterial core OTUs) (Fig. 5). This indicates that approximately 55% of the total core OTUs in each domain were found in at least 12 samples (and, therefore, were detected on at least two sampling dates), and exhibited sufficiently strong correlations for inclusion in the network. These nodes were connected by 274 edges in total: 85 bacteria-fungi (29 positive and 56 negative), 183 fungi-fungi (115 positive and 68 negative), and 6 bacteria-bacteria (3 positive and 3 negative) correlations.

Plant leaves are inhabited by diverse microorganisms, which can influence plant health and productivity (Stone et al. 2018). Despite a growing understanding of the important roles these microbes play in plant and ecosystem ecology and their recognized potential in crop improvement programs, the factors influencing the phyllosphere microbiome and their functional consequences are not well understood. As a step toward a more holistic view of microbial community dynamics in the phyllosphere, we investigated fungal community structure on switchgrass leaves throughout a growing season, and related those seasonal dynamics to the dynamics of core bacterial community members that were previously reported from the same experiment (Grady et al. 2019).

Fungal communities exhibit lineage-specific temporal patterns in relative abundance. Switchgrass phyllosphere fungal community composition differed strongly across the growing season and in response to fertilization. It should be noted that edge effects can strongly influence the plant microbiome such as through increased incidence of foliar pathogens (Johnson and Haddad 2011). Therefore, edge effects could at least partially explain the apparent fertilization effect on fungal communities given the layout of the experimental plots. In any case, seasonal effects were stronger than fertilization effects on community composition, and the primary goal of our study was to establish broad patterns (i.e., regardless of fertilization levels) in phyllosphere core taxa. Therefore, fertilized and nonfertilized samples were pooled for the remainder of our analyses.

Interestingly, fungal community shifts across the growing season were at least partially phylogenetically conserved: core fungal taxa exhibited distinct temporal patterns of abundance at higher taxonomic levels. For example, members of Dothideomycetes (Ascomycota) were at highest abundance early and late in the growing season. Dothideomycetes has been widely as a dominant phyllosphere taxa globally (Abdelfattah et al. 2015;Izuno et al. 2016;Qian et al. 2018;Yao et al. 2019), is the largest and most diverse class of fungi, and includes numerous agriculturally relevant pathogens as its best-studied members but also includes endophytes and mutualists (Schoch et al. 2009). Other groups of fungi Fig. 3. Seasonal patterns in switchgrass phyllosphere fungal a diversity: A, unique fungal species accumulation, where each point is the mean of four to eight replicate leaves; B, species richness; C, Shannon diversity; and D, Pielou's evenness. Operational taxonomic units were defined at 97% amplicon sequence identity. Data are based on subsampling depth of 48,088 sequences per sample. Fig. 4. Dynamics of switchgrass core phyllosphere members. A, Abundance-occupancy of switchgrass phyllosphere taxa and their inclusion in the core. Taxa that were detected in all replicate plants in at least one sampling date were considered to be in the core (green). Noncore taxa (white or open circles) are also indicated. B, Patterns of core taxa that share similar temporal changes, as determined by hierarchical clustering of standardized dynamics. OTUs = operational taxonomic units. C, Patterns of core taxa relative abundances aggregated by fungal genus (i.e., each genus contains one or more OTUs) considering only genera where mean relative abundance was ³ 0.01 for at least one sampling date. For panels B and C, the y-axis indicates the mean abundance of core taxa relative to the total (both core and noncore) number of rarefied reads, while error bars indicate one standard deviation from the mean. Data are based on subsampling depth of 48,088 sequences per sample.

exhibiting seasonal patterns include members of Exobasidiomycetes (Basidiomycota), which often exhibit a yeast phase and were in highest abundance midseason, the converse of the pattern seen in Dothideomycetes. Indeed, a recent in vitro study of bean phyllosphere isolates reported that Dothideomycetes phyllosphere isolates possess effective antagonistic mechanisms toward other fungi, likely due to the secretion of inhibitory metabolites (Prior et al. 2017), which might explain, in part, the lower overall a diversity when Dothideomycetes isolates were in high abundance. Nevertheless, it is important to note that members of Dothideomycetes are extremely diverse, both ecologically and functionally (Ohm et al. 2012). At the genus level, the core fungal microbiome included Aureobasidium and Alternaria spp., which have been previously reported as belonging to the root core microbiome of switchgrass (Singer et al. 2019). Alternaria alternata, a known switchgrass pathogen (Vu et al. 2012), was a particularly dominant core taxon that comprised an average of 10.4% (± 4.9% standard deviation) of the rarefied fungal community dataset across all phyllosphere samples, while Tilletia, a genus of smut fungi that includes known switchgrass pathogens (Gravert et al. 2000), was highly abundant midseason. The core phyllosphere microbiome also included Aureobasidium pullulans, a known source of hemicellulosedegrading enzymes (Myburgh et al. 1991). Fungal species richness was highest midsummer, mirroring bacterial species richness, as previously reported (Grady et al. 2019). This likely reflects species accumulation through the bulk of the growing season, followed by microbial mortality and reduced diversity as leaves senesced in autumn.

Leaf traits associated with plant growth strategies explain variation in fungal communities. In addition to seasonal effects, several plant traits and seasonal environmental characteristics significantly explained variation in phyllosphere fungal communities, including leaf nitrogen concentration, leaf dry matter content, and average plant height. These findings are in partial agreement with previous work in this system, where plant height was associated with variation in bacterial communities (Grady et al. 2019). In addition, our results are similar to a recent study of 51 tropical tree species which found that variation in phyllosphere fungal communities was significantly associated with variation in leaf nitrogen concentration and leaf mass per area (a composite trait dependent upon leaf dry matter content and leaf thickness) (Kembel and Mueller 2014). However, it is important to note that these results could be confounded by the correlation between sampling date and plant traits. For example, later sampling dates (i.e., older plants) were associated with greater plant height, lower leaf nitrogen concentration, and higher leaf dry matter content, in agreement with recent work on ontogenetic shifts in leaf functional traits (Mason et al. 2013). Unfortunately, the design of our study does not allow us to distinguish between sampling date and factors correlated with sampling date as the causative agent shaping fungal phyllosphere communities. Nevertheless, the strong differences in communities across the growing season, with clear lineage-specific trajectories, suggests the potential for moving phyllosphere microbiome research toward a more predictive science.

Fungal-bacterial co-occurrences are more common than bacteria-bacteria in the switchgrass phyllosphere. In a network analysis considering only core fungal and bacterial OTUs in the phyllosphere, over half of the fungal and bacterial core phyllosphere OTUs merited inclusion in the intra-and interdomain microbial correlation network. The most highly connected fungal core OTU (20 connections) was identified as belonging to genus Taphrina (class Taphrinomycetes), while Sphaerellopsis filum (19 connections) and Pseudorobillarda phragmitis (16 connections) were also hub taxa. Interestingly, S. filum is a mycoparasite of Puccinia emaculata, one of the causative agents of leaf rust in switchgrass, and has been suggested as a potential biocontrol agent (Black 2012). The most highly connected bacterial core OTU (15 connections) was an unidentified Proteobacterium, while a member of genus Pseudomonas (14 connections) and a member of Sphingomonas, two common phyllosphere taxa in Arabidopsis (Bodenhausen et al. 2013;Delmotte et al. 2009), were also bacterial hub taxa.

Among phyllosphere core taxa associations in the microbial cooccurrence network, fungal-fungal network correlations were quite common, with positive correlations more common than negative (115 and 68, respectively, out of 183 total connections). On the other hand, fungal-bacterial associations were more likely to be negative than positive (56 and 29, respectively, out of 85 total connections), in agreement with recent reports (Agler et al. 2016;Durán et al. 2018). Work in the Arabidopsis root microbiome suggests that these negative associations reflect direct or indirect bacterial biocontrol of detrimental fungi, with important implications for plant health and survival (Durán et al. 2018). Bacteria were more frequently associated with fungi than with other bacteria (85 fungal-bacterial associations versus only 6 bacterial-bacterial associations), likely due, at least in part, to the sixfold greater number of fungal core taxa compared with bacterial core taxa (227 versus 41, respectively). The 85 fungal-bacterial connections detected here included 22 distinct bacterial OTUs, 6 of which were identified as Hymenobacter (phylum Bacteroidetes, family Cytophagaceae), which has previously been reported as a dominant endophyte in several plant species, including switchgrass (Ding and Melcher 2016), and 4 of which were identified as Sphingomonas (phylum Proteobacteria, family Sphingomonadaceae). Of the 54 distinct fungal core OTUs involved in fungal-bacterial associations, 24 belonged to class Dothideomycetes, with 16 of those belonging to order Pleosporales. Altogether, cross-domain connections comprised approximately 31% (85 of 274) of all network connections, roughly three times as many as found in the rhizosphere of the annual crop common bean, where cross-domain interactions comprised 9% of all interactions, and the same methods were used to determine core taxa (Stopnisek and Shade 2019). We speculate that the greater diversity of core fungal-bacterial interactions found here may be related to the life history of the plant. For example, perennial crops may provide more fungal-bacterial interaction opportunities than annual crops due to legacy effects of an established cropping system that has enriched for core microbiome members. However, these differences could also reflect fundamental differences between aboveground and belowground processes, and more work involving cross-plant comparisons must be done to test this hypothesis.

In summary, we found a strong impact of sampling date on phyllosphere epiphytic fungal community composition, with clear lineage-specific temporal patterns in relative abundance. We also found that leaf nitrogen concentration and dry matter content were associated with variation in fungal communities, potentially indicating the phyllosphere microbiome as a component of plant growth strategies. Finally, interdomain associations were substantially more than intradomain associations for bacteria, and negative fungalbacterial associations were common, suggesting the importance of high-level phylogenetic diversity and negative feedback processes in core microbiome assembly. Although our findings highlight the complexity of phyllosphere microbiome assembly, the clear temporal patterns in lineage-specific fungal abundances suggest the potential for accurately predicting shifts in fungal phyllosphere communities throughout the growing season, a key research priority for sustainable agriculture.