(C) PLOS One This story was originally published by PLOS One and is unaltered. . . . . . . . . . . Host genotype controls ecological change in the leaf fungal microbiome [1] ['Acer Vanwallendael', 'Department Of Plant Biology', 'Michigan State University', 'East Lansing', 'Michigan', 'United States Of America', 'Ecology', 'Evolution', 'Behavior', 'Great Lakes Bioenergy Research Center'] Date: 2022-08 Leaf fungal microbiomes can be fundamental drivers of host plant success, as they contain pathogens that devastate crop plants and taxa that enhance nutrient uptake, discourage herbivory, and antagonize pathogens. We measured leaf fungal diversity with amplicon sequencing across an entire growing season in a diversity panel of switchgrass (Panicum virgatum). We also sampled a replicated subset of genotypes across 3 additional sites to compare the importance of time, space, ecology, and genetics. We found a strong successional pattern in the microbiome shaped both by host genetics and environmental factors. Further, we used genome-wide association (GWA) mapping and RNA sequencing to show that 3 cysteine-rich receptor-like kinases (crRLKs) were linked to a genetic locus associated with microbiome structure. We confirmed GWAS results in an independent set of genotypes for both the internal transcribed spacer (ITS) and large subunit (LSU) ribosomal DNA markers. Fungal pathogens were central to microbial covariance networks, and genotypes susceptible to pathogens differed in their expression of the 3 crRLKs, suggesting that host immune genes are a principal means of controlling the entire leaf microbiome. Funding: The plant material is based upon work supported in part by the Great Lakes Bioenergy Research Center (glbrc.org), U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research under Award Number DE-SC0018409. Support for this research was provided by the National Science Foundation Long-term Ecological Research Program (DEB 1832042; https://beta.nsf.gov/funding/opportunities/long-term-ecological-research-lter ) at the Kellogg Biological Station and by Michigan State University AgBioResearch. We received further funding from the U.S. Department of Energy (energy.gov) through grants DE-SC0014156 to TEJ and DE-SC0017883 to DBL and National Science PGRP Award IOS1402393 to J.T.L. The work conducted by the U.S. Department of Energy Joint Genome Institute ( https://jgi.doe.gov/ ) is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. We hypothesized that the phyllosphere fungal microbiome develops seasonally as a successional community controlled by environmental factors, host genetics, and interspecific fungal–fungal associations. We used amplicon sequencing to compare the relative importance of these factors in the phyllosphere fungi of a replicated diversity panel of switchgrass (Panicum virgatum [ 31 ]). We tested whether communities change directionally and whether the trajectory of succession differed across switchgrass genetic subpopulations and across different growing sites. Additionally, we sought to uncover whether specific genetic loci underlie host control of the microbiome through genome-wide association study (GWAS) and RNA sequencing analyses. Finally, we investigated the roles of specific fungal taxa in the microbiome through network analysis. Specifically, we aimed to determine whether known switchgrass leaf pathogens [ 32 ] covary with nonpathogenic symbionts, or are peripheral to microbial communities. The phyllosphere microbiome, consisting of the microbes on and inside the plant leaf, comprises diverse taxa that impact plant health and productivity [ 15 – 18 ]. Leaf fungi in particular are common plant pathogens [ 19 ], but nonpathogenic taxa may perform beneficial functions for the host, including nutrient uptake and pathogen antagonism [ 20 – 25 ]. Since the phyllosphere microbiome of perennial plants is reassembled at the start of each growing season in freshly sprouted tissues, [ 26 , 27 ] it may show similar patterns to macro-scale secondary successional communities. Recent research has shown that host control of the leaf microbiome is often governed by numerous loci of small effect directly impacting relatively few microbes [ 28 – 30 ]. In the case of microbiomes, host factors governing microbial succession must also be considered. Since the composition of the microbiome can greatly impact host fitness, it can be evolutionarily beneficial for the host to play a role in the successional process, encouraging mutualist colonization while dispelling pathogens as the community assembles. Hosts express genes that influence colonizing microbes through several means, including immunity, morphological adaptations [ 12 ], and chemical exudation [ 13 ]. While the immune system is often effective at preventing detrimental infections, immune receptors may recognize and exclude beneficial microbes if elicitors are structurally similar to a pathogen, so specific immunity can have wider impacts on the microbiome [ 14 ]. Hosts require finely calibrated mechanisms for attracting beneficial microbes without attracting pathogens in a constant coevolutionary push and pull. Communities that colonize available niches in the process of succession follow certain predictable ecological patterns. Early-arriving species are typically those with effective long-range dispersal, while the climax community is dominated by species that can more effectively use resources under competition [ 6 ]. While these broad patterns are generalizable, the composition of any particular successional community depends greatly on both the habitat colonized and interspecific interactions such as priority effects, where the order of arrival of taxa governs the success of later arrivals [ 7 , 8 ]. While most successional theory is based on studies in macro-scale organisms, the principles of succession are evident in microbial communities as well, but on a more rapid timescale [ 9 – 11 ]. Microbial communities perform essential functions for their host organisms in all branches of life. In some systems, hosts can tightly control the microbes with which they form symbioses [ 1 , 2 ]. In others, the composition of the microbiome is more governed by ecological interactions such as the order of species arrival or abiotic conditions during colonization [ 3 , 4 ]. A key goal of microbial evolutionary ecology is to determine how both host and nonhost factors influence microbiome assembly [ 5 ], particularly in natural settings where host influence is more challenging to study. Our analyses mostly identified the taxa that were abundant across samples. Rare taxa can be important in microbial community functioning [ 42 ], but their role in overall ecological patterns is less clear and more challenging to study. Therefore, we only examined rare taxa that we expected a priori to play an important ecological role. Claviceps species were present in 119/760 samples, and were more highly abundant in the early season. Claviceps species produce alkaloid compounds that deter grazing [ 43 ], so this endophyte may play a role in protecting young grass shoots. Metarhizium, a related genus, was present at low abundances in 43/760 samples in the Columbia, Missouri and KBS, Michigan sites. Metarhizium species are insect-pathogenic fungi [ 44 ], so may provide a similar protective role. While patterns in this core group reflected major changes in the fungal microbiome, we used several alternate methods to identify important OTUs. In addition to the core, we modified trajectory analyses ( Fig 1B ) by computationally removing each OTU from the analysis and calculating the change in the overall community trajectory [ 34 ]. Nineteen OTUs significantly impacted trajectories when they were removed, all of which overlapped with the core group ( Fig 5C ). To examine priority effects, we used microbial temporal variability mixed linear models (MTV-LMMs), which identify taxa for which variation in earlier time points explains variation in later points [ 41 ]. Of the 153 OTUs we found in this analysis, 49 overlapped with the core, and 14 with both the core and trajectory analysis ( Fig 5C ). The 14 OTUs that were identified as important using all 3 methods ( Fig 5C and 5D ) included taxa from several putative functional guilds, including yeasts, pathogens, and mycoparasites. Network connections confirmed mycoparasitic interactions; we found a negative relationship between putative plant pathogens Mycosphaerella tassiana and Microdochium seminicola versus mycoparasitic Epicoccum dendrobii. In addition, we used indicator species analysis to identify OTUs that were overrepresented in leaves with fungal disease symptoms ( Fig 5A ). Although the major fungal pathogen of switchgrass, Puccinia novopanici, was not identified as a core taxon, indicator species analysis showed that putative mycoparasite Sphaerellopsis filum is present in the core and significantly associated with fungal infection symptoms (OTU_4; Fig 5A ). In addition to varying among functional groups, OTU covariance also significantly changed over time as supported by the bootstrap-permutation based network comparisons between each sampling point ( Fig 5A and S6 and S2 Tables). To identify positive or negative covariance temporal patterns within network members, we generated a Class-level heatmap showing the proportion of edges linking OTUs within or between each Class at each time point ( S8 Fig ). Due to the high proportions of Dothideomycetes (mixed guilds) and Tremellomycetes (yeast) in the core, the majority of edges at every time point were within (38.3% to 48.6%) and between (9.9% to 14.5%) OTUs in these classes. While the proportion of positive edges maintained more or less stable with time between OTUs in the Dothideomycetes (from 20.0% to 23.9%) and Tremellomycetes (from 23.0% to 17.4%) or within the 2 classes (from 3.6% to 2.8%), negative edges between classes increased from DOY 158 (6.7%) to DOY 233 (9.7%) and DOY 286 (11.7%). This may indicate competition for host resources between these 2 classes of fungi, resulting in more spatially heterogeneous distributions in the late season. (A) Covariance networks of core OTUs over time. Nodes are colored by each OTU’s relative abundance in infected leaves with visible symptoms. The shape of the node denotes network position, defined by Zi-Pi ratio. Edges are colored by the covariance sign. (B) Infection indicator taxa, including best taxonomic match and z-score for indicator analysis. (C) Number of OTUs identified as important by several methods: MTV-LMM analyses that indicate time-dependent OTUs, OTUs that impact the successional trajectory, and core OTUs with high occupancy-abundance. (D) Taxonomic information for the 14 OTUs identified in all 3 analyses in (C). Best match denotes the lowest taxonomic level confidently identified for each OTU using BLAST. Guilds were estimated based on published studies, references are in S1 Table . (E) Network statistics for fungal guilds, calculated as mean values across all time points, with (SD. Data underlying this figure can be found in S5 Data . DOY, day of year; MTV-LMM, microbial temporal variability mixed linear model; OTU, operational taxonomic unit; SD, standard deviation. To investigate how these functional guilds associate, we built covariance networks using OTU relative abundances at each time point ( Fig 5A ). We summarized network statistics across functional guilds to show that known grass pathogens are central to covariance networks, with high betweenness centrality (extent to which a node lies on paths connecting other nodes) and degree (overall number of connections; Fig 5E ). Standard deviation was high within this group, however, reflecting seasonal and within-group differences. Yeasts, in contrast, showed higher modularity (compartmentalization; Fig 5E ). This indicates that, while yeasts are overall more speciose in the core microbiome, they covary less with the rest of the microbial community than pathogens. Since yeasts are thought to be mostly commensal inhabitants of the outer leaf surface [ 40 ], this difference may reflect their ecological or spatial niche. Given the large differences in leaf pathogen susceptibility across switchgrass subpopulations, we sought to determine the influence of pathogenic fungi on other members of the fungal microbiome. We examined the taxonomic relationships of the 7392 OTUs in our dataset through a hybrid method that compares matches across fungal databases and BLAST (Basic Local Alignment Search Tool) hits [ 37 ]. We identified 6,756 OTUs as fungi, and excluded 633 plant, and 3 metazoan OTUs. We performed NMDS and PERMANOVA analyses using the full fungal community, but focused our taxon-specific analyses on OTUs at the focal site that were present at high occupancy across time and showed relatively high abundance, often referred to as the “core” microbiome [ 38 ]. This group consisted of 128 OTUs, the majority of which were Dothideomycetes (43.5%) and Tremellomycetes (28.7%, S1 Table ). We assigned each of the core OTUs to a functional guild when possible using published literature ( S1 Table ). Of the core group, 23 OTUs were grass pathogens, and 9 were documented pathogens of other plants. Four were known mycoparasites, fungi that prey upon other fungi. Three were generalist decomposers or had an unclear functional guild, and the remaining 52 were yeasts or yeast-like fungi. Compared to fungal species in soil, these taxa were especially enriched for grass pathogens and yeasts and contained much fewer saprophytes [ 39 ]. As further confirmation of the importance of the outlier SNP, we selected several genotypes that were not in the original study (n = 20) containing differing alleles of the outlier SNP. We used the same protocols as the first round of sampling, but only sampled at one location at 2 time points, KBS, Michigan. Fungal microbiomes in these samples conformed closely to our predictions, with allelic variation at the Chr02N_57831909 locus influencing microbiome structure (PERMANOVA F = 1.84, p = 0.016, R 2 = 0.064; S7 Fig ). Although this was a smaller subset of samples taken 3 years after the original samples, we were also able to confirm a strong temporal influence on microbiome structure between the 2 time points (PERMANOVA F = 16.80, p < 0.001, R 2 = 0.293; S7 Fig ). In addition, we sequenced the fungal large subunit (LSU) in these samples to exclude any primer bias against any taxonomic group. LSU sequences showed a similar temporal pattern to fungal ITS and were structured with the Chr02N_57831909 locus in a strikingly similar pattern to ITS (PERMANOVA F = 2.66, p = 0.001, R 2 = 0.086; S7 Fig ). For both amplicons, population structure as indicated by PCs 1 and 2 of the genomic SVD explain some variation in the microbiome, but not as much as allelic variation at Chr02N_57831909. Together, these results suggest that the Chr02N_57831909 locus has a stable and deep influence on the microbial community. We corroborated the importance of these candidate genes by comparing their expression levels in divergent genotypes at the 3 of the 4 sites where phyllosphere experiments were conducted, KBS, Columbia, and Austin. In each site, we sequenced leaf tissue RNA from multiple biological replicates (n ≥ 3) from 4 genotypes: 2 that are typically more susceptible to leaf fungal pathogens (Midwest upland VS16 and DAC) and 2 that typically are more resistant (Gulf lowland WBC and AP13) [ 32 ]. Consistent with host-gene driven variation in fungal community assembly, all 3 candidate genes were much more highly expressed in susceptible than resistant genotypes (Wald tests; Table 2 ). These differential genotype-specific patterns of expression were very similar across planting sites (likelihood ratio test for ecotype ✕ site interaction, p = 0.354). (A) GWAs of microbiome structure (NMDS2 at DOY 260). The lower solid line shows the 5% FDR threshold, and the upper dashed line shows the Bonferroni-adjusted alpha threshold for SNPs associated with the microbiome. (B) Outlier region on Chromosome 2N with nearby genes shown in red. (C) Expression-level differences for genes shown in (B). Leaf tissue for these samples was collected as part of a different study, performed at 3 of the same sites we used. FPKM are scaled differently in each gene facet. Data underlying this figure can be found in S4 Data . DOY, day of year; FDR, false discovery rate; FPKM, fragments per kilobase of transcript per million mapped reads; GWA, genome-wide association; NMDS, nonmetric multidimensional scaling; SNP, single nucleotide polymorphism. Such tight host–microbiome genetic diversity associations imply a genetic basis of influence on fungal community dynamics by their plant hosts. To identify the genetic loci that might underlie this pattern, we calculated genome-wide associations (GWAs) for microbiome structure. We used the second NMDS axis at DOY 260 from the above analysis (Figs 1 and S4 ) to represent microbiome structure, since it showed the greatest clustering with subpopulation (Figs 1A and 3C , additional time points in S5 Fig ) and controlled for large-scale host genetic structure by including a single variate decomposition of pairwise genetic distance as a covariate in the linear models. We found several loci associated with the phenotype at a 5% false discovery rate (FDR), but the GWA showed an excess of low p-values (quantile–quantile plot: S6 Fig ) so we used a more conservative Bonferroni-corrected threshold to identify significant SNPs ( Fig 4A ). This threshold revealed only 1 SNP on chromosome 2N significantly associated with microbiome structure, Chr02N_57831909. This SNP is closely linked to several genes in the switchgrass v5.1 genome annotation ( Fig 4B and 4C ). The 3 closest genes are homologous to receptor-like kinases (RLKs) annotated in the closely related Panicum hallii (2 copies of cysteine-rich receptor-like protein kinase 6; XP_025800480.1 and XP_025800481.1, and 1 copy of cysteine-rich receptor-like protein kinase 10; XP_025801715.1). This class of RLKs is diverse in plants, but is known to contain many immune receptors [ 36 ], indicating a potential role for these genes in host control of fungi. (A) Pairwise genetic distance (π) for all samples. Samples are ordered by hierarchical clustering. (B, C) Pairwise community distance (Bray–Curtis) for all samples, shown in the same order as genetic distances, for 2 sampling times, DOY 158 and DOY 260. Other sampling times shown in Supp. (D) Values of Mantel’s r shown indicate correlation between distance matrices for genetics and fungal communities at each time point, with subsetted subpopulations shown as faint lines. p < 0.01 for all tests in the combined populations. Data underlying this figure can be found in S3 Data . DOY, day of year; KBS, Kellogg Biological Station. Beyond differences at the level of subpopulations, we expected that within-subpopulation genetic differences would impact fungal diversity. To further examine genetic differences over time, we compared host genetic distances to fungal community differences between plants at the focal site, KBS. Genetic distances, calculated as Nei’s distance using 10.2 million single nucleotide polymorphisms (SNPs) [ 35 ], revealed that host population genetic structure largely matched the 3 major switchgrass genetic groups observed previously: “Gulf,” “Atlantic,” and “Midwest” [ 31 ] ( Fig 3A ). These 3 subpopulations are deeply diverged and serve as discrete gene pools within which we tested for host-driven fungal community divergence. Fungal community distances, calculated as Bray–Curtis community dissimilarity, varied across sampling dates, but largely recapitulated the genetic structure of switchgrass ( Fig 3B and 3C ). Notably, Mantel tests showed that fungal community structure was most closely correlated with host genetic structure at DOY 260, when most plants had set seed (r = 0.453), but declined as senescence progressed ( Fig 3D ). While we anticipated some degree of genetic influence, subpopulations were even more highly structured than expected, with almost half of the variation in fungal community distance explained by genetic distance when plants are setting seed (DOY 260). To confirm this pattern, we estimated pseudo-heritability values for community structure using the kinship matrix as a random effect in mixed models. Overall, there was weak heritability for variation on the second NMDS axis (H 2 = 0.132 ± 0.043), but this may be attributable to high variation across time points (DOY 158: 0.120 ± 0.188, DOY 212: 0.164 ± 0.146, DOY 233: 0.184 ± 0.733, DOY 260: 0.950 ± 0.176, DOY 286: 0.913 ± 0.352). Genotypes were sampled at 4 sites. From north to south: KBS, Michigan; Columbia, Missouri; Austin, Texas; and Kingsville, Texas. Northern sites are shown by symbols, and southern by open shapes. Color indicates phenological stage sampled, “Early” samples were taken just after emergence, “Mid” samples were taken during seed development, and “Late” samples were taken after senescence began. Data underlying this figure can be found in S2 Data . NMDS, nonmetric multidimensional scaling. Fungal microbiomes can be greatly influenced by environmental factors in addition to host factors. Therefore, we compared succession across environments by selecting a subset of 8 plant genotypes replicated in 3 additional sites across a latitudinal gradient in the USA. From north to south, these field sites were Columbia, Missouri; Austin, Texas; and Kingsville, Texas ( S2 Fig ). We sampled each site at 3 time points, standardized by phenology to account for seasonal differences across sites ( S3 Fig ). At most sites, collection date correlated with both NMDS1 and NMDS2 ( Fig 2 ; stress = 0.103). However, the northern and southern sites were divided on a diagonal line orthogonal to collection date. The northern sites KBS and Columbia, Missouri formed one cluster, while the southern sites, Austin, Texas and Kingsville, Texas formed another ( Fig 2 ). Differences across sites accounted for 29.6% of the variation in community dissimilarity across sites, but sampling date and subpopulation also structured the community to a lesser extent ( Table 1 ). While succession may show temporal patterns in southern sites, the composition of fungal communities on leaves is largely distinct. A significant site:sampling date interaction term indicates that succession differs across sites, but additional samples may be needed to fully understand this pattern. In order to directly test the differences in succession across subpopulations, we modeled changes in the multidimensional representation of fungal communities as directional trajectories [ 34 ]. Across the season, fungal communities on individual plants showed parallel changes over time, with almost no reversals to earlier community states ( Fig 1B ), strongly indicating a successional pattern. Switchgrass genetic subpopulations differed in both mean trajectory length (df = 3, F = 2.786; p = 0.0453) and mean overall direction (df = 3, F = 3.677; p = 0.0151). While little subpopulation difference is evident at the beginning of the season, climax fungal communities were markedly different in the Midwestern population, which showed the greatest divergence from others in trajectory direction ( Fig 1B , Midwest-Atlantic; Tukey honest significant difference [HSD] = 0.013, p = 0.050). This provided initial evidence that, while fungal dispersal is similar across plant subpopulations, host plants influence the climax state of fungal communities. Further, we tested for differences between genotypes that were presenting rust-associated symptoms and genotypes that were not presenting symptoms. Infected and symptomless plants changed along parallel trajectories that differed in length (df = 1, F = 5.274, p = 0.0238), but not in directionality (df = 1, F = 0.652, p = 0.421). Each point represents an individual plant sampled from the experimental plot at KBS, Michigan. (A) NMDS. Dates are shown as DOY. Points are colored by DOY, and switchgrass subpopulations as shapes. (B) Trajectory plots of principal coordinates of community distances. Transparent arrows represent individual switchgrass genotypes sampled over the 5 dates shown in (A), and colors show genetic subpopulations. Solid colored arrows show mean subpopulation trajectories. Data underlying this figure can be found in S1 Data . DOY, day of year; KBS, Kellogg Biological Station; NMDS, nonmetric multidimensional scaling. To determine the directionality of successional changes in the microbiome, we visualized community differences with nonmetric multidimensional scaling (NMDS). NMDS accurately preserved sample distances in reduced dimensions (Stress = 0.102; S1 Fig ) and revealed clear temporal community structure. NMDS1 clustered closely with the date of collection, while NMDS2 clustered more with host genetic subpopulation ( Fig 1A ). Notably, the first sampling date was highly distinct from the later time points, showing greater variation within that time point, as well as divergence from later time points ( Fig 1 ). To explore the statistical significance of patterns of succession, we used permutational multivariate analysis of variance (PERMANOVA). Both sampling date (day of year, DOY) and subpopulation had significant effects on community structure, but differed greatly in their explanatory power ( Table 1 ). At the focal site, KBS, collection date (DOY) explained the greatest amount of variation (19.4%), followed by genetic subpopulation (5.7%), and there was a significant date-by-subpopulation interaction (2%). To assess the impact of disease on leaf microbiome, we performed a separate test with the subset of samples for which we were able to collect both infected and symptomless leaves. We restricted permutations within individual plants to perform the equivalent of a paired test of infection effects, resulting in a significant infection term that explained 7.79% of the variation in community distance (p < 0.001). Switchgrass is a highly genetically diverse perennial grass native to North America, and both plant traits and switchgrass–microbe interactions vary across its range [ 31 – 33 ]. We leveraged this diversity to assess the difference in microbial communities across the 3 main switchgrass subpopulations by randomly selecting 106 genotypes from a diversity panel [ 31 ] planted at our focal site, the Kellogg Biological Station (KBS), Michigan, United States of America. Of these, 28 genotypes were from the Midwestern subpopulation, 38 from the Atlantic, 31 from the Gulf, and 9 showed signs of admixture between groups (Intermediate). These subpopulations differ in morphological and ecological characteristics, so we expected that fungal succession would differ as well across subpopulations. Switchgrass subpopulations correspond roughly to 3 morphological ecotypes: Midwest genotypes are mostly Upland, Gulf are Lowland, and Atlantic mostly Coastal [ 31 ]. Since the samples used in this study largely followed this pattern, subpopulation differences can also be considered ecotypic differences. We examined succession over time by sampling leaf tissue from each plant at 5 time points, then sequencing the internal transcribed spacer (ITS) region of the phyllosphere-associated fungi in and on the leaf. After quality filtering, we clustered 47.8 million ITS reads to 6,756 fungal operational taxonomic units (OTUs) that varied across genotypes and over time. Discussion Our results show strong support for the importance of time, geographic location and host genetics in influencing the switchgrass phyllosphere microbial succession over the growing season. We found evidence for clear successional dynamics that were consistent in direction across growing sites, but were distinct in community composition. Fungal communities were different across host genetic subpopulations, a pattern that may be driven by variation at 3 linked immune receptors. Leaf fungal communities are taxonomically diverse, but a few highly abundant pathogens and yeast species play a disproportionate role in shaping community progression. Viewing the switchgrass leaf microbial community through the lens of succession allowed us to delineate ecological patterns in these communities. Multidimensional scaling representations of the leaf communities at the focal site revealed a clear clustering by date of collection on the first NMDS axis (Fig 1). This indicates that, as we predicted, date of collection is an important source of variation in the switchgrass leaf fungal community. Further, measuring the trajectories of these communities showed that succession is both directional and deterministic, since no samples showed negative trajectories (reversals of succession) by the end of the season, and most samples followed a similar trajectory (Fig 1B). We observed similar patterns to other studies that show early-season leaves as highly distinct from later time points, perhaps owing to greater influence from soil microbes [16,45]. While the overall shape of trajectories was similar among samples, the Midwestern population deviated from others, particularly in the late season. The Midwestern population is notable since we have previously shown that it is more susceptible to several fungal pathogens such as leaf rust (Puccinia novopanici) and leaf spot (Bipolaris spp.; [32]; also see [46]) and has on average an earlier phenology than the other population groups [33]. Leaf microbiome relationships are consistently distinct in this population and may be linked to other traits such as cold tolerance that also differ [31,33]. In addition to temporal differences across subpopulations, the composition of fungal leaf communities differed markedly across geographic locations. This may be partially due to seasonality differences across the region we examined. The Kingsville, Texas site did not experience freezing temperatures between 1989 and 2020 (NOAA weather service), so perennial grasses in the region may have living aboveground tissue year-round. Growing season length has been shown as an important factor in governing the abundance and diversity of endophytic fungi [47], so it is unsurprising that we saw large differences across this latitudinal gradient. However, many other factors that influence fungal communities also differ across these sites, including precipitation regime, soil type, and surrounding vegetation, so further work is needed to determine if the growing season is truly the causal factor. We predicted that fungal communities would be impacted by host genetics as well as location. We found several lines of evidence for genetic control of the leaf microbiome. In addition to examining differing successional trajectories across subpopulations, we tested the covariance of genetic distance and fungal community differences using Mantel correlations. Genetic-fungal community correlations increased until DOY 260, then declined as host senescence began. Mantel tests are inappropriate for some ecological tests and often underestimate p-values, but can be useful for exploratory analysis of distance matrices [48]. While there was high variation in our pseudo-heritability estimates, the fact that they mirror temporal patterns in the Mantel tests strengthens the general trend of greater genetic associations in the late season. Previous studies have found similarly high variation in microbiome heritability estimates across time [29], so it is not surprising to see this in our case. Deng and colleagues calculated H2 for individual OTUs, which ranged from 0 to 0.66 and a Mantel’s correlation of r = 0.13 between genetic and microbiome composition in sorghum rhizosphere. This value of r is lower than we saw in our study, possibly since it was based on a relatively small subset of samples. When selecting samples for this study, we randomly chose equal numbers of samples from the 2 major switchgrass morphological ecotypes, upland and lowland switchgrass [49] (S2 Fig). Lowland switchgrass, which is more highly represented in Gulf and Atlantic subpopulations, is more resistant to several leaf fungal pathogens [32], so subpopulation differences may be at least partially driven by differences in immunity across these genotypes. Since pathogens such as Microdochium and Alternaria were among the most abundant taxa in our samples, their differences across subpopulations may have driven overall community differences. In addition to immunity, however, subpopulations differ in other traits that may contribute to fungal colonization differences, such as leaf wax content [50], exudate concentration [51], and phenology [33,49], so microbiome differences may be responding to multiple host plant traits. A replicated receptor-like kinase is associated with fungal differences We found one outlier SNP associated with microbiome structure. While there were several peaks in the Manhattan plot (Fig 5A), our analysis showed a strongly skewed distribution of observed versus expected p-value (S6 Fig), indicating a risk of Type I errors. This is probably attributable to the low sample size in this GWAS. The influence of the identified locus is fairly strong, contributing to a clear decrease on NMDS axis 2 when the minor allele is present (MAF = 0.083; S9 Fig). This SNP is not in Hardy–Weinberg equilibrium in switchgrass; we found only one minor-allele homozygote among our samples. This abnormal pattern may be attributable to structural variation at this locus. Switchgrass subpopulations vary widely in genome structure, which may result in alignment mismatches that resemble SNPs, particularly in regions with multiple gene copies [52]. Indeed, this region shows an elevated number of insertions and deletions compared to nearby sections of the 2N chromosome (S10 Fig, data from [31]) and is adjacent to a region dense with repetitive long terminal repeat retroelements (positions 60960000–60980000). Given the confirmatory results for this locus as well as the RNA sequencing results, however, we expect that there is a true phenotypic association with the locus, but that it may be with a structural variant rather than a true SNP. The 3 nearby genes we identified were replicated variants of a cysteine-rich RLK whose function has not been experimentally verified in Panicum. RLKs are one of the largest plant gene families, including over 600 members in Arabidopsis [36]. The best studied of these is FLS2, which detects the bacterial flagellin protein and initiates an immune response cascade [53]. The 3 RLKs we identified show high sequence similarity to immune-related cysteine-rich RLKs in Arabidopsis and Oryza, and contain the “stress-antifungal domain” PF01657, which has been linked to salt stress as well as fungal responses when present in several proteins [54,55]. Arabidopsis CRK5, for example, alters defense responses either through resistance to infection or programmed cell death, depending on how the gene is expressed [56], and CRK6 and CRK 14 are involved in the general non–self-response [57]. Related Arabidopsis genes may be the targets of immune repression by bacterial strains [58]. Similarly, the Oryza gene LIL1 (Os07g0488400) improves fungal rice blast resistance when overexpressed [59]. The pattern of these receptors being more highly expressed in pathogen-susceptible plants may seem counterintuitive given that many RLKs are immune receptors. However, this can often occur when pathogens produce effector proteins that target immune receptors [60]. Necrotrophic fungi in particular can benefit by over-inducing plant immune receptors to initiate programmed cell death, [61,62] then feeding on dead plant tissue. Since allelic variants at this locus have now been associated with variation in the fungal microbiome across several years in natural populations, it may represent a useful target for future research into genetic control of the leaf microbiome. Previous research has shown that microbiome control is often polygenic, with many contributing loci of small effect [28–30]. Uncovering only a single causal locus in this study may be a product of the relatively low sample size; there are more loci associated with microbiome community structure that did not meet the GWAS cutoff, but may contribute to a polygenic architecture for this trait. [END] --- [1] Url: https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3001681 Published and (C) by PLOS One Content appears here under this condition or license: Creative Commons - Attribution BY 4.0. via Magical.Fish Gopher News Feeds: gopher://magical.fish/1/feeds/news/plosone/