Multi-element fingerprinting and high throughput sequencing identify multiple elements that affect fungal communities in Quercus macrocarpa foliage

Diverse fungal mutualists, pathogens and saprobes colonize plant leaves. These fungi face a complex environment, in which stochastic dispersal interplays with abiotic and biotic filters. However, identification of the specific factors that drive the community assembly seems unattainable. We mined two broad data sets and identified chemical elements, to which dominant molecular operational taxonomic units (OTUs) in the foliage of a native tree respond most extremely. While many associations could be identified, potential complicating issues emerged. Those were related to unevenly distributed OTU frequency data, a large number of potentially explanatory variables and the disproportionate effects of outlier observations.Key words: community assembly, environmental filter, fungi, heavy metal enrichment, nutrient enrichment, oak, Quercus macrocarpaHyperdiverse fungal communities inhabit the foliage of most plants^1,2 and these fungal communities have been reported for virtually every plant that has been examined.³ Baas-Becking hypothesis states that environment selects microbial communities from the abundant and possibly globally distributed propagule pools.⁴ Although the foliage-associated communities—like other microbial communities—are suspected to be sensitive to environmental drivers, determination of the mechanisms that control the assembly of these foliar communities has remained difficult and elusive. Some of the proposed mechanisms include distance limitations to propagule dispersal,^5–7 volume limitations to propagule loads,⁷ or limitations set by the environmental conditions either on the scale of the site of fungal colonization⁸ or more broadly on a landscape level.^6,9 The forces that may control the fungal community assembly are overlaid by additional biotic controls that include compatibilities between the fungi and host species^10,11 or genotypes^6,12 and the competitive or facilitative interactions among the component fungal genotypes.^6,10–13 Although a variety of potential controls for the foliage-associated fungal communities have been speculated, very little consensus exists on the relative importance of the different drivers. For example, while macronutrient and heavy metal enrichment may have an influence on the composition fungal communities¹⁴ and populations,¹⁵ relative importance of various chemical elements in the foliage remains yet to be investigated.To evaluate the use of multi-element fingerprinting data produced by Inductively Coupled Plasma Mass Spectrometry (ICP-MS) in combination with high throughput 454-pyrosequencing for determining influential chemical elements in structuring of the leaf-associated fungal communities, we mined a recent dataset¹⁶ that explored the effects of urbanization on the diversity and composition of the fungal communities associated with a native tree Quercus macrocarpa. From a total list of more than 700 non-singleton fungal OTUs, we selected fifty with highest overall frequency to provide an observationrich dataset for elemental effect assessment; these OTUs accounted for 84.5% of all sequences. Even so, many of these OTUs had a number of zero frequencies (Fig. 1), highlighting one of the difficulties in the use of environmental sequencing data. We omitted one OTU (OTU630 with a likely affinity to Trimmatostroma cordae Mycosphaerellaceae]) that was strongly affected by the original land use design (urbanization; Wilcoxon rank sum test with a Bonferroni adjustment) and therefore unlikely to be representative for the present analyses of elemental drivers. This OTU was replaced with one with the next highest frequency. The frequencies of these 50 OTUs were investigated in the context of concentrations for 29 elements after the omission of five (Ag, Au, C, δ¹³C, δ¹⁵N) in the final analyses because of their strong association with the land use or the difficulty of finding a biological relevance. Of the remaining elements three (Fe, Cr and Ni) had pairwise correlations exceeding 0.98 between the three pairings; others showed no similar high correlations. To allow comparable evaluation across the broad array of elements, all concentrations were standardized to have a mean equal to zero and a standard deviation equal to one.Open in a separate windowFigure 1Rank-ordered distribution of observed frequencies for those OTU s whose frequency had an extreme slope when associated with the concentrations of one or more chemical elements in the mixed effects model. The asterisk denotes one extreme frequency for OTU 313 with a value 0.8636. Numbers in parentheses indicate the number of observations with a frequency equal to zero. The OTU s were assigned to approximate taxa using BLAST:²⁰ 425: Alternaria alternata (Pleosporaceae); 46: Phoma glomerata (Pleosporaceae); 686: Aureobasidium pullulans (Dothioraceae); 520: Davidiella tassiana (Davidiellaceae); 567: Cladosporioum tenuissimum (Davidiellaceae); 313 Oidium heveae (mitosporic Erysiphaceae); 586: Erysiphe hypogena (Erysiphaceae); 671: Mycosphaerella microsora (Mycosphaerellaceae); 555: Pestalotiopsis sp. (Amphisphaeriaceae); 607: Pleiochaeta setosa (incertae sedis).To rank elements according to their magnitude of association with the abundance of each OTU, a total of 1,450 models (50 OTUs times 29 elements) relating element concentration to OTU abundance were fit to the data. For each model, OTU frequency was the dependent variable, element concentration and time (a factor with three levels) were fixed effects, and—to account for the spatial arrangement of the experimental units—random effects associated with tree nested within site were included in the error structure. Time by element interactions were also investigated and tested using a likelihood ratio test. These mixed effect models were fit using R and the package lme4 (www.rproject.org).Statistical “tests of significance” that produce p-values can be sensitive to assumptions or outliers. Because of this and the fact that our analyses evaluated a total of 1,450 models, p-values themselves were not considered a reliable measure of importance when associating elements with OTU frequency. Instead, we emphasized metrics that highlight extraordinary findings rather than rely on tests of statistical significance. This approach facilitates finding few elements that have the strongest effect on OTU frequency. Note that the use of standardized element concentrations (above) provided slope coefficients that are comparable across all models. “Extreme slopes”, i.e., models where the OTU response to element concentration was strongest, were identified as those with estimated slope coefficients in the lower or upper 2.5 percentile, i.e., those farther than 1.8 standard deviations from the mean across all estimated slopes (Fig. 2). Using this approach, we identified a total of 69 models with extreme slopes (Open in a separate windowFigure 2Distribution of estimated slopes (i.e., the slope for element concentration) for a model relating OTU frequency to element concentration, time and a concentration by time interaction, including a tree-nested-within-site random effect. The mean across all 1,450 OTU s is approximately zero; the two vertical lines identify upper and lower 2.5 percentiles, beyond which the slopes were considered extreme (large black symbols). The horizontal line identifies the cut off maximum leverage (0.24), above which the slopes were considered to have observations with high leverage. Models with observations with a high leverage were tested for extreme slopes by refitting without those observations. Models are ranked from bottom to top in order of increasing leverage and the element for which the high-leverage observations and extreme slopes were recorded are identified on the right y-axis.

Table 1

Slopes identified as extreme in our analyses

Element	OTU 425	OTU 46	OTU 686	OTU 520	OTU 567	OTU 313	OTU 586	OTU 671	OTU 555	OTU 607
B	+*	+*					+*
Ba						−	−
Ca			−*			(−)*	−*	(+)*	+**
Cd	+					+	(+)
Ce	+	(+)						−
Co	+**						−*
Cr							−*
Cu	+*	−**					−*
Fe							−*
Hg	+**					−*
K	(−)		+	−		+		(−)		(+)
Li	(+)*						(+)*	−*
Mn		+*
Mo							−*
N	−*					+*	(+)*
Na	+
Ni							−*
P	−*									(+)*
Pb		+**					−*
Rb		+**	+*	−*			−*
S	(−)*	+*				+*	+*
Sc	−					(−)	−
Se		−
Sn	−					−	(−)
Sr		+*
Y		+*			−*	+*	(+)*
Zn	(−)*	+*				−**				(+)*

Open in a separate windowPositive slopes are indicated by +, negative by −. Parentheses indicate where a statistically significant (α = 0.05) interaction was observed (likelihood ratio test). Extreme slopes with observations with high leverage are identified by an asterisk (*) and those where omission of high-leverage observations lead to a non-extreme slopes are identified by two asterisks (**). Note that eight of the ten OTU s in the table had an extreme slope with at least one element concentration after accounting for high leverage and interactions in the model.Unfortunately, the models with extreme slopes were often affected by high leverage observations (outliers in the explanatory variables) that may have exerted substantial influence on the magnitudes of the slopes. We accounted for this by computing leverage values based on the fixed effect model matrix (element concentration and time) for each model. High leverage was defined as those observations with leverage approximately twice the mean leverage over all samples for a particular model as is considered conventional by some authors.¹⁷ This value was approximately 0.24 for our models. The models with high leverage and extreme slopes were re-evaluated by refitting the model to the data after omission of the influential observations. Of the 69 models with extreme slopes only 22 were void of influential observations by our metric (Fig. 1). Our analyses included the possibility of identifying those models that were affected by numerous low frequencies and a few high frequency observations. We argue that the few higher frequencies are most likely indicative of those elements that also have extreme concentrations in the same samples; we did not want to miss such findings. Second, no one element controls the occurrence of all or even majority, of the OTUs, but the OTUs appear to respond positively or negatively to different drivers. This is strongly visible even among the eight that remained through our rigorous evaluation of a vast number of models. This can be interpreted in the context of a niche. Foliage represents a complex abiotic physicochemical habitat within which organisms are sorted based by stochastic arrival parameters, but also by either environmental tolerances or nutritional preferences. Those fungi best able to colonize and invade the available substrate under any given combination of the complex physical and chemical environmental matrix will persist and be detected most frequently. Thirdly, even for one OTU, many elements may have strong and occasionally opposing effects. For example, for OTU425, B, Cd, Ce, Cu, Na, had positive effects, whereas N, P, Sc had negative effects (18^,19 it is tempting to speculate on species replacement or on tolerance to nutrient enrichment as a result of changes in the abiotic chemical environment. However, one must exercise caution: as we point out above, a number of other alternative factors come to play when a correlative relationship like this is considered across two discrete and complex datasets. Several heavy metal concentrations also showed either positive or negative associations with the fungal OTU frequencies. To exemplify, the frequencies of OTUs 313 and 425 were positively associated with the concentrations of Cd and OTU 46 was positively associated with Zn, whereas OTUs 313 and 586 were negatively associated Hg and Pb concentrations, respectively. Does this mean that these species differ in their sensitivities to these particular heavy metals? Not necessarily, but these observational data provide a starting point for more explicit hypothesis-driven experiments that allow for specific elucidation of the fungal responses to these elements and may guide future experimentation.We conducted a high-dimensional exploratory analysis to evaluate potential effects of element concentration on OTU frequencies. Using a repeated measures mixed effects model, we were able to compile a brief list of chemical elements with the most likely (based on these data) strongest effects on the abundances of the dominant components of the phyllosphere-associated fungal communities. Complicating the use of usual methods of statistical inference (i.e., use of p-values) was the sparseness in the occurrence of many OTUs across samples and outlying observations in the concentration of some elements. We chose the extreme slopes approach that allowed ranking associations between OTU frequency and element concentration with no assumptions regarding normality or equivariance that may be violated using traditional tools of inference (e.g., Analysis of Variance). Still, some of the observed associations may have been affected by extreme leverage points (outliers in the explanatory variables) and these were accounted for in the present analyses by model re-evaluation after omission of the high-leverage observations. While our analyses identified a number of biologically meaningful associations between chemical elements and molecular OTUs, rigorous experimentation is mandatory to establish causative relationships.