Ectomycorrhizal Fungal Communities and Their Functional Traits Mediate Plant–Soil Interactions in Trace Element Contaminated Soils

fpls-09-01682 November 16, 2018 Time: 20:55 # 1 ORIGINAL RESEARCH published: 20 November 2018 doi: 10.3389/fpls.2018.01682 Edited by: Stefano Castiglione, University of Salerno, Italy Reviewed by: Raffaella Balestrini, Consiglio Nazionale delle Ricerche (CNR), Italy Sabine Dagmar Zimmermann, Centre National de la Recherche Scientifique (CNRS), France *Correspondence: Marta Gil-Martínez [email protected] †These authors have contributed equally to this work Specialty section: This article was submitted to Plant Microbe Interactions, a section of the journal Frontiers in Plant Science Received: 14 June 2018 Accepted: 29 October 2018 Published: 20 November 2018 Citation: Gil-Martínez M, López-García Á, Domínguez MT, Navarro-Fernández CM, Kjøller R, Tibbett M and Marañón T (2018) Ectomycorrhizal Fungal Communities and Their Functional Traits Mediate Plant–Soil Interactions in Trace Element Contaminated Soils. Front. Plant Sci. 9:1682. doi: 10.3389/fpls.2018.01682 Ectomycorrhizal Fungal Communities and Their Functional Traits Mediate Plant–Soil Interactions in Trace Element Contaminated Soils Marta Gil-Martínez1*†, Álvaro López-García2†, María T. Domínguez3, Carmen M. Navarro-Fernández1, Rasmus Kjøller2, Mark Tibbett4and Teodoro Marañón1 1Department for Protection of the Soil, Plant and Water System, Institute of Natural Resources and Agrobiology of Seville, Spanish National Research Council, Seville, Spain, 2Department of Biology, University of Copenhagen, Copenhagen, Denmark, 3Área de Edafología y Química Agricola, Departamento de Cristalografía, Mineralogía y Química Agrícola, Universidad de Sevilla, Seville, Spain, 4Centre for Agri-Environmental Research and Soil Research Centre, School of Agriculture, Policy and Development, University of Reading, Reading, United Kingdom There is an increasing consensus that microbial communities have an important role in mediating ecosystem processes. Trait-based ecology predicts that the impact of the microbial communities on ecosystem functions will be mediated by the expression of their traits at community level. The link between the response of microbial community traits to environmental conditions and its effect on plant functioning is a gap in most current microbial ecology studies. In this study, we analyzed functional traits of ectomycorrhizal fungal species in order to understand the importance of their community assembly for the soil–plant relationships in holm oak trees (Quercus ilex subsp. ballota) growing in a gradient of exposure to anthropogenic trace element (TE) contamination after a metalliferous tailings spill. Particularly, we addressed how the ectomycorrhizal composition and morphological traits at community level mediate plant response to TE contamination and its capacity for phytoremediation. Ectomycorrhizal fungal taxonomy and functional diversity explained a high proportion of variance of tree functional traits, both in roots and leaves. Trees where ectomycorrhizal fungal communities were dominated by the abundant taxa Hebeloma cavipes and Thelephora terrestris showed a conservative root economics spectrum, while trees colonized by rare taxa presented a resource acquisition strategy. Conservative roots presented ectomycorrhizal functional traits characterized by high rhizomorphs formation and low melanization which may be driven by resource limitation. Soil-to-root transfer of TEs was explained substantially by the ectomycorrhizal fungal species composition, with the highest transfer found in trees whose roots were colonized by Hebeloma cavipes. Leaf phosphorus was related to ectomycorrhizal species composition, specifically higher leaf phosphorus was related to the root colonization by Thelephora terrestris. These findings support that ectomycorrhizal fungal community composition and their functional traits mediate plant performance in metal-contaminated soils, and have a high influence on plant capacity for phytoremediation of contaminants. The study also corroborates the Frontiers in Plant Science | www.frontiersin.org 1November 2018 | Volume 9 | Article 1682 fpls-09-01682 November 16, 2018 Time: 20:55 # 2 Gil-Martínez et al. Ectomycorrhizal Fungi Mediate Plant–Soil Interactions overall effects of ectomycorrhizal fungi on ecosystem functioning through their mediation over the plant economics spectrum. Keywords: ecosystem processes, heavy metal, microbiome, phytoremediation, Quercus ilex subsp. ballota (holm oak), root economics spectrum, symbiosis, trace element transfer INTRODUCTION There is an increasing consensus that microbial communities have an important role in mediating ecosystem processes. In recent years, and thanks to the development of molecular approaches, several studies have focused on the interaction between plants and soil microbial communities to reveal the potential of microbes to drive vegetation diversity and dynamics (Bever, 2003;Wardle et al., 2004;van der Heijden et al., 2015; Erktan et al., 2018;Rutten and Gómez-Aparicio, 2018). As vegetation determines how ecosystems function to a large extent, plant microbiomes indirectly affect the provision of multiple ecosystem services (Friesen et al., 2011;Van der Putten et al., 2013). In addition, some studies have highlighted the existence of feedback processes between plants and soil organisms (Bever et al., 2010;Brinkman et al., 2010), suggesting not only the potential of microbes to modify plant communities but also the role of plant communities and their traits at structuring microbial community compositions (de Vries et al., 2012; Aponte et al., 2013;Bauman et al., 2016;López-García et al., 2017). Although the effect of plant hosts on their microbiomes has often been studied from a taxonomic point of view (Aponte et al., 2010;de Vries et al., 2012;Kurm et al., 2018), little is known about how soil microbial functional traits are affecting the functioning of plant species. It is debatable whether the features of microbes associated to individual plants (species composition and trait distribution) can be actually defined as plant traits, as they are not heritable features, according to the definition of Garnier et al. (2016). Often, microbial traits in the root microbiome are referred as “biotic root traits” (Bardgett et al., 2014). Recently, some authors have considered the use of traits in the root microbiome as an extension of the plant species phenotype for explaining functional changes in plant communities along environmental gradients and it has been included in a multidimensional root trait framework (NavarroFernández et al., 2016;Weemstra et al., 2016). The influence of the plant microbiome from a trait-based perspective usually requires assessment of individual species in communities (Díaz et al., 2007), and this has been proven to be very challenging when working with microbes (see Crowther et al., 2014). According to the current thinking on ecological assembly, recording traits at individual species level will allow to differentiate between response and effect traits (Zirbel et al., 2017). The links between the response of microbial communitylevel traits to environmental conditions and the effects of these microbial traits changes on plant functioning is an important knowledge gap to be filled in current microbial ecology studies, although the existence of these links have been predicted previously (see Koide et al., 2014). Mycorrhizal fungi are recognized for their importance for plant foraging of soil resources (Tibbett and Sanders, 2002;van der Heijden et al., 2015;Köhler et al., 2018), particularly in plant species with relatively thick absorptive roots (Eissenstat et al., 2015;Liu et al., 2015). Coevolution of plant and fungal partners has been recently suggested by Chen et al. (2018), based on their description of a root-fungal functional complementarity in nutrient foraging. However, how mycorrhizal and plant traits are interrelated, for example aligned into the common root economics spectrum framework, and how mycorrhizal traits mediates soil–plant relationships are still open questions that need to be addressed (Weemstra et al., 2016). This mycorrhiza–root association improves plant health by enhancing resistance to diverse stresses like drought, salinity, heavy metals and pathogens, among others (van der Heijden et al., 2015). Therefore, mycorrhizal mediation on plant performance might be especially important in highly stressful environments, such as trace element (TE) contaminated soils. In these soils, mycorrhizal fungi enhance plant nutrition, stress tolerance and soil structure and, consequently, promote the recovery of the functions in the degraded soil (van der Heijden and Scheublin, 2007;Firmin et al., 2015). Association with mycorrhizal fungi can also play an important role in the transfer of TEs through the soil–root continuum, an issue of special relevance for the management of TE contaminated sites. For instance, the phytostabilization approach is a phytoremediation technology that combines the use of soil amendments and plants to immobilize pollutants into the soil, thus reducing the risks of transfer of these pollutants through the aboveground food web (Mendez and Maier, 2008). A prerequisite to apply this approach to large contaminated areas is that the plants used to remediate the soil can retain TEs at the rhizosphere level, and do not accumulate them into their aboveground biomass (Bolan et al., 2011;Madejón et al., 2018a). In relation to this, ectomycorrhizal (ECM) fungi may provide protection against metal toxicity through avoidance (i.e., extracellular precipitation, biosorption to cell walls, reduced uptake) and sequestration (i.e., intracellular chelation, compartmentation into fungi vacuoles) (Hartley et al., 1997;Jentschke and Godbold, 2000;Bellion et al., 2006). Therefore, phytoremediation of TE polluted soils can be facilitated by ECM fungi as they adapt to TE stress promoting the host growth (Wen et al., 2017). In this study, we aimed to elucidate the role of ECM community in the plant nutritional status and the transfer of TEs through the soil-root-leaf continuum in a largescale phytoremediation case study. Holm oaks (Quercus ilex subsp. ballota) root and leaf functional traits were analyzed in trees growing on remediated soils exhibiting a gradient of anthropogenic TE contamination. Relationships between Frontiers in Plant Science | www.frontiersin.org 2November 2018 | Volume 9 | Article 1682 fpls-09-01682 November 16, 2018 Time: 20:55 # 3 Gil-Martínez et al. Ectomycorrhizal Fungi Mediate Plant–Soil Interactions nutrient/TE concentrations in plants and the structure on the ECM communities were evaluated. Ectomycorrhizal community composition and morphological traits along the same gradient of soil contamination were previously reported by López-García et al. (2018). Here, relationships between ECM, soil, root, and leaf variables were explored in order to understand the importance of the ECM community assembly in the soil–plant relationships in holm oak trees. We hypothesized that (i) plant traits (i.e., morphological and chemical) of holm oak would change along TE gradient; (ii) ECM fungal communities, would partly mediate plant response to TEs, and thus a significant fraction of the plant nutrient status and transfer of TEs from soils to leaves will be explained by ECM variables (either species composition or functional traits) (iii) ECM fungal communities lead the intraspecific variation of root functional traits. MATERIALS AND METHODS Study Area The study was conducted at the Guadiamar Green Corridor (SW Spain), an area affected by a large mining accident in 1998 (the Aznalcóllar mine spill; Madejón et al., 2018b). The failure of a large tailings storage facility was one of the largest mining accidents in Europe to the date, which provoked the release of ca. 6 hm3of metalliferous tailings (water and sludge) over 55 km2of the Guadiamar River basin. As a result, soils were severely polluted with several TEs, mainly As, Cd, Cu, Pb, Tl, and Zn (Cabrera et al., 1999). After the accident, a large scale soil remediation program was conducted, which included the removal of the deposited sludge and the soil surface using heavy machinery, followed by the application of organic matter and calcium-rich amendments to immobilize TEs into the soil. The affected lands, mostly under agricultural production until the mining accident, were purchased by the Regional Administration, and then afforested using native tree and shrub species (Domínguez et al., 2008). Despite these remediation operations, contamination levels were highly variable across short distances in the Corridor and some patches are still highly degraded, due to acid drainage of the remnants of the sludge that lead to soil acidification and to a high solubility of toxic TEs (Domínguez et al., 2016). The climate of the study area is typically Mediterranean, with mild rainy winters and warm dry summers. Average annual temperature is 19◦C (minimum monthly mean of 9◦C in January, and maximum of 27◦C in July) and annual average rainfall is 484 mm. The study area harbors soils with different geology adding additional variation to the patchily distributed levels of TEs. Typical bedrock types at the North of the Corridor are slate and schist, and derived soils are naturally acidic. In the South (further than 15 km away from the mine) geological substrate tends to be dominated by calcarenite and marls originating neutral to calcareous loam soils. Potential vegetation is dominated by sclerophyllous Mediterranean forests, in particular by ECM holm oak in the alluvial terraces. Sampling Design The study was conducted in April 2016, 16 years after the application of soil amendments and the plantation of the former agricultural lands with native trees and shrubs. Holm oak was the target species of the study, given that it was intensively used to afforest the alluvial terraces of the affected area. Four sites were selected along a gradient of soil pollution across the affected area. A site location map and a general description of these soils as well as their classification is provided in López-García et al. (2018). Site 1 and Site 3 were located at the North of the corridor, while Site 2 and site 4 were located at the South of the Corridor. At each site, 10 holm oak trees were randomly selected (N= 40 trees). All these trees were planted at the same time (Autumn 2000) and with similar seed provenance. For each tree, roots (and their associated ECM fungi) were sampled by carefully tracing from the stems of the tree to the roots belowground in the four cardinal directions. Around 200 g of root material was collected from each direction, i.e., subsamples. Root samples were used to characterize the main root functional traits and the ECM community (see López-García et al., 2018 for ECM characterization). Soil samples (0–20 cm depth) were taken with an auger from the four directions under each tree canopy and were pooled to a total of 500 g to make a composite sample per tree. Likewise, fully expanded leaf samples were taken from the four cardinal directions of the tree canopies to obtain a composite sample of leaves for each tree. Soil Chemical Analyses Soil chemical analyses were conducted for the study reported in López-García et al. (2018). Soil samples were air-dried and sieved to <2 mm for chemical analysis. Soil pH, Ca, K, P, NH4, NO3, total C, total N and total TEs were measured, following the methodologies described in that paper. For this study, available concentrations of S and TEs were also analyzed. Sulfur and TEs were extracted from samples (<60 µm) with a 0.01 M CaCl2 solution (Houba et al., 2000) and analyzed by inductively coupled plasma spectrophotometry (ICP-OES) using a Varian ICP 720-ES (simultaneous ICP-OES with axially viewed plasma). Soil Enzyme Activities The activity of three extracellular enzymes involved respectively in C, N, and P cycling [β-glucosidase (BGL), N-acetylglucosaminidase (NAG), and acid phosphatase (ACP)] were measured as indicators of microbial activity in the collected soils. These enzymes were analyzed colorimetrically by incubation with p-nitrophenyl-linked substrates at 37◦C for 1 h, according to Tabatabai and Bremner (1969); Tabatabai (1982), Parham and Deng (2000), respective methods. ECM Species Composition and Functional Traits Molecular analysis of ECM in root samples, as well as quantification of ECM functional traits (abundance of rhizomorphs, emanating hyphae, and melanin content) were conducted by López-García et al. (2018). Briefly, a composite sample of 28 root fragments per tree was obtained by selecting Frontiers in Plant Science | www.frontiersin.org 3November 2018 | Volume 9 | Article 1682 fpls-09-01682 November 16, 2018 Time: 20:55 # 4 Gil-Martínez et al. Ectomycorrhizal Fungi Mediate Plant–Soil Interactions the seven longest root fragments in each of the four root subsamples collected from each tree. A random individual root tip per root fragment was photographed for posterior trait quantification (presence of emanating hyphae and rhizomorph and colorimetric estimation of melanization, see López-García et al., 2018, Appendix 1). Community weighted means (CWMs), i.e., the averaged value for these traits per tree, was calculated as the number of root tips exhibiting emanating hyphae or rhizomorphs divided by the total number of quantified root tips (Lepš et al., 2011). The color value was averaged between the 28 root tips of each tree for having an overall estimation of the ECM melanization of the community. The remaining material was used for the quantification of the percentage of root length colonized by ECM fungi, using the gridline intersect method (Brundrett et al., 1996;Navarro-Fernández et al., 2016). All these data was reported by López-García et al. (2018) and was included in the statistical analyses in order to evaluate the influence of ECM communities on holm oak status. A small portion of each root tip was cut and immersed separately into 10 µl of Extraction Solution (Extract-N-AmpTM Plant PCR Kit by Sigma-Aldrich) and the protocol of the manufacturer was followed to extract its DNA. PCR amplification was conducted using primers ITS1F (Gardes and Bruns, 1993) and ITS4 (White et al., 1990) following the procedure described in López-García et al. (2018), and Sanger sequenced. Sequences were blasted against the UNITE database (Kõljalg et al., 2013) and those found to correspond to ECM fungi were grouped by genera or family (see López-García et al., 2018 for details) and compared against the UNITE database (Kõljalg et al., 2013) for their taxonomic placement and Species Hypothesis determination. The number of root tips belonging to each root was used as abundance data. Plant Functional Traits Root and leaf functional traits were measured specifically for this study, following the protocol described in Pérez-Harguindeguy et al. (2013). Morphological root traits included specific root length (SRL), specific root area (SRA), and root dry matter content (RDMC) and were measured with WinRHIZO 2009 (Regent Instruments, Quebec, CA, United States). Specific leaf area (SLA) and leaf dry matter content (LDMC) were measured in a subsample of 10 leaves per tree: leaves were scanned and analyzed with Image-Pro 4.5 (Media Cybernetic, Rockville, MD, United States). After sampling, we selected 10 separated leaves from each tree and washed them with deionized water to determine the Chlorophyll Content Index (CCI) with a SPAD-502 chlorophyll meter (Minolta Camera, Co. Ltd., Osaka, Japan) taking three measurements per leaf. Subsamples of roots and leaves collected from each tree were used for chemical analysis. This root material can be considered as the symbiotic combination of plant and fungi tissues. These subsamples were washed with distilled water, dried at 70◦C for at least 48 h, and ground. Total C and N were determined by using a Flash 2000 HT elemental analyzer (Thermo Scientific, Bremen, Germany). Trace elements (As, Cd, Cu, Fe, Mn, Ni, Pb, and Zn) and macronutrients (S, P, K, Ca, and Mg) were determined by ICP-OES after digestion of plant tissues by wet oxidation with concentrated HNO3in a Digiprep MS block digester (Domínguez et al., 2008). Data Analysis In order to explore the relationships among ECM and plant variables we conducted a preliminary selection of key variables to be included in subsequent multivariate and modeling analyses. As the aim of the work was to evaluate whether plant performance (nutrient status and TE accumulation) is mediated by ECM communities in these soils, the subset of variables used as predictor variables included soil background properties and TEs, ECM species composition and ECM traits. The subset of response variables included TEs transfer from soil to root and leaves, root traits and leaf traits (Supplementary Figure 1). Soil variables and the characterization of ECM community, published by LópezGarcía et al., 2018, were used for the analysis of the present study. A preliminary analysis of variance (ANOVA) to compare differences in soil, root, and leaf variables among sampling sites was performed. We checked for normality and homoscedasticity of data, and when assumptions were not met data were log or square root transformed. When these assumptions were met a Tukey’s Honest post hoc test followed. Otherwise, a nonparametric Kruskal–Wallis test and a Dunn’s test corrected by Bonferroni post hoc were performed. Due to the dataset complexity, and in order to remove correlations and to reduce collinearity between soil variables, a principal component analysis (PCA) was performed to select a non-collinear subset of soil TEs to be used as predictors of plant traits in subsequent statistical analysis. Original data was logtransformed for normalization. Most correlated TEs with the first two axes of each PCA were selected for subsequent analyses. In order to reduce ECM fungal species composition into two dimensions, a principal coordinate analysis (PCoA) was performed with the operational taxonomic units (OTUs) matrix (Legendre and Gallagher, 2001). The first two PCoA axes were selected (Supplementary Table 1) (Pinheiro and Bates, 2000; Zuur, 2009). To evaluate the influence of soil and ECM variables on plant nutritional status and its functional traits we applied both correlational analysis and linear mixed models. In order to understand the relationships between the response and predictor variables, we first performed Pearson’s correlation tests, adjusted with Benjamini-Hochberg correction (Benjamini and Hochberg, 1995). Those soil and ECM variables showing a significant correlation with plant variables were considered as fixed effect factors in univariate linear mixed effect models, with sampling site as random factor. The significant variables from univariate models were included additively in multivariate models, however variance inflated factors (VIFs) were calculated and variables with VIF >3 were removed to avoid collinear predictors (Zuur et al., 2010). Models were compared against a null model, assuming no influence of any of these predictors on plant variables. The best and most parsimonious predictive models were selected based on the Akaike information criterion corrected for small sample sizes (AICc; Burnham and Anderson, 2002). Selected models were Frontiers in Plant Science | www.frontiersin.org 4November 2018 | Volume 9 | Article 1682 fpls-09-01682 November 16, 2018 Time: 20:55 # 5 Gil-Martínez et al. Ectomycorrhizal Fungi Mediate Plant–Soil Interactions FIGURE 1 | Principal component analysis (PCA) ordination of (A) soil total trace elements (TEs) and (B) soil available TEs, sampled at 0-20 cm depth under holm oak trees (N= 40) and classification by site. fitted, and marginal and conditional R2values were computed. Marginal R2(R2LMMm) is the variance explained by fixed factors, while conditional R2(R2LMMc) is variance explained by both fixed and random factors (Nakagawa and Schielzeth, 2013). Requirements for normality and homoscedasticity of residuals were fulfilled in all the selected models. All statistical analyses were carried out using the R software v.3.3.2 (R Core Team, 2016), using packages ggplot2 (Wickham, 2009), MuMIn (Barton, 2017), nlme (Pinheiro et al., 2016), psych (Revelle, 2017), and vegan (Oksanen et al., 2016). RESULTS Soil Characterization As reported by López-García et al. (2018), soil pH was significantly different among sites; sites 2 and 4 showed a significantly higher pH than site 3 and, specially, than site 1 (Supplementary Table 2). Available TEs levels decreased from Site 1 to Site 4. About soil nutrients, sites 1 and 3 showed significantly higher NH4, NO3, total N and organic C than sites 1 and 4. Calcium concentration was significantly higher at site 2 with respect to the other sites. Phosphorus contents were not significantly different among sites (Supplementary Table 2). All soil enzyme activities presented the highest activity at site 3 and NAG and ACP activities were found significantly lower at site 2 (Supplementary Table 2). Reduction of Trace Element and Community Composition Variables for Model Analysis Soil total TEs PC1 and PC2 ordination axis explained most of the total variance (86.68%) in the chemical composition of soils (Figure 1A). Axis 1 and 2 represented the variation of two clear groups of TEs which were orthogonal to each other. Axis 1 correlated well with total As, Cd, Cu, Pb, S and Zn, which tended to covariate. Axis 2 showed a high covariance between Mn and Ni. Likewise, the first two axes of available TEs explained most of the total variance (83.24%) (Figure 1B), being Zn and Mn the most correlated with axes PC1 and PC2 respectively. The final selected TEs included in the subsequent analyses were: total As, Fe and Mn, and available Mn and Zn concentrations. Available Cd was not chosen because some of the samples were below the detection limits. Lower guideline values (LGVs) for contaminated soils (Ministry of the Environment Finland, 2007) were exceeded for As, Cu, and Pb at site 1 (Supplementary Table 2). Fifty five OTUs were recorded belonging to ECM fungal species in 494 successfully sequenced root tips (published in López-García et al., 2018). In summary, these taxa comprised 14 families and 19 genera. The presence of rare species was common among the study: 19 of 55 OTUs were only identified in one root tip (Supplementary Figure 2). Two species, Hebeloma cavipes and Thelephora terrestris dominated the communities with 83 and 61 root tips, respectively (Supplementary Figure 2). The first two axes of PCoA of the ECM fungal communities explained a 25.08% of the variance in community composition. PCoA axis 1 (13.36%) showed a gradient from rare to abundant species (Figure 2). A clear pattern was also found in PCoA axis 2 (11.72% of explained variance) showing a transition of ECM fungal communities from Thelephora terrestris to Hebeloma cavipes. Relationships Between Soil and ECM Variables and Root Traits In general, nutritional root status was found to be more affected by biotic factors than by abiotic ones when univariate models Frontiers in Plant Science | www.frontiersin.org 5November 2018 | Volume 9 | Article 1682 fpls-09-01682 November 16, 2018 Time: 20:55 # 6 Gil-Martínez et al. Ectomycorrhizal Fungi Mediate Plant–Soil Interactions FIGURE 2 | Principal coordinate analysis ordination of ECM fungal species in symbioses with holm oak roots and classification by site. Arrows indicate signicant ECM species (p<0.05). were run. Root C was the variable that was best explained by the considered predictors (Table 1). Soil Ca and available Mn (Estimate = 0.10) were those variables explaining the greatest variation in root C (univariate models), followed closely by melanization (Estimate = −0.23) (Table 1). Species composition PCoA1 (Estimate = −1.13), PCoA2 (Estimate = −1.14), and rhizomorph formation (Estimate = 0.08) also presented an effect on root C but to a lesser extent. Soil Ca content presented a negative effect on C root concentration. ECM species composition (PCoA1 and PCoA2 scores) was also related to root C; a higher root C was observed where Hebeloma cavipes species was dominant. In terms of biotic CWM traits, a low melanin content and high rhizomorph presence were also affecting root C (univariate models). Root N was best explained by single total Zn, which exerted a positive effect (Table 1). In terms of biotic effects, species composition PCoA1 had a marginally significance influence (Estimate = 0.05). Root C:N ratio corroborated the role of total Zn on root N, as Zn was found negatively correlated to C:N ratio. The best model for root C:N included species composition PCoA1, exerting a negative effect, which confirmed that the presence of the most abundant species, Hebeloma cavipes and Thelephora terrestris, negatively influenced root N content. Root P was not explained by any abiotic factor but was positively affected by ECM species composition PCoA1 (Figure 3A and Table 1). As well as root N, root P was found to be lower when Hebeloma cavipes and Thelephora terrestris species were abundant in roots, therefore the symbioses with other species, here considered as rare due to their lower abundance, probably improved the nutritional status of holm oak roots in terms of P. Root N:P ratio was not significantly explained by any abiotic or biotic factor. Both morphological root traits, SRA and RDMC were better explained by species composition PCoA1 than by any abiotic TABLE 1 | Univariate and multivariate linear mixed models showing significant soil and ECM fungi fixed effects for each of the root traits and model explained variance. Response variable Individual effects of soil factors Individual effects of ECM fungal factors Combining significant effects into the best predictive model Nutrients and EA Trace elements PCoA axis 1 and 2 Fungal traits CWM Linear mixed effect models Variance Root Variable pVariable pVariable pVariable pModel SE t p R2LMMm R2LMMc C (%) Ca <0.001 Av. Mn <0.001 PCoA1 0.048 Melanization 0.003 C = 45.43 – 0.001 Ca 0.0003 −4.11 <0.001 0.31 0.31 PCoA2 0.044 Rhizomorph 0.047 N (%) – – Total Zn 0.006 (PCoA1) (0.053) – – N = 0.25 +0.0006 Total Zn 0.0002 2.90 0.006 0.19 0.19 C:N – – Total Zn 0.010 PCoA1 0.012 – – C:N = 141.70 – 29.64 PCoA1 11.07 −2.68 0.012 0.22 0.40 P (%) – – – – PCoA1 0.009 – – P = 0.09 +0.02 PCoA1 0.007 2.80 0.009 0.21 0.22 SRA (m2kg−1) – – – – PCoA1 0.045 – – SRA = 9.47 +0.62 PCoA1 0.30 2.10 0.045 0.12 0.12 RDMC (mg g−1) – – – – PCoA1 0.006 – – RDMC = 427.04 – 22.04 PCoA1 7.40 −2.98 0.006 0.22 0.22 SRA, specific root area; RDMC, root dry matter content; Av., available; CWM, community weighted mean; EA, enzyme activity; SE, standard error; R2LMMm, marginal variance; R2LMMc, conditional variance. Frontiers in Plant Science | www.frontiersin.org 6November 2018 | Volume 9 | Article 1682 fpls-09-01682 November 16, 2018 Time: 20:55 # 7 Gil-Martínez et al. Ectomycorrhizal Fungi Mediate Plant–Soil Interactions FIGURE 3 | Relationship between selected key ECM fungal species composition and traits and their effects on (A) root P, (B) root dry matter content, (C) leaf P, (D) translocation factor Fe, (E) soil-to-root Fe transfer, and (F) soil-to-root Zn transfer in studied sites. Frontiers in Plant Science | www.frontiersin.org 7November 2018 | Volume 9 | Article 1682 fpls-09-01682 November 16, 2018 Time: 20:55 # 8 Gil-Martínez et al. Ectomycorrhizal Fungi Mediate Plant–Soil Interactions factors, but their effects were opposite (Table 1). The presence of the most abundant species, Hebeloma cavipes and Thelephora terrestris reduced SRA but increased RDMC (Figure 3B), while SRL was not significantly explained by any of the measured soil or ECM fungal factors. In summary, a key effect of species composition was found for those variables related to root nutrition variables. The abundance of Thelephora terrestris and, in special, Hebeloma cavipes species seemed to be related to high C, C:N ratio and RDMC values, and low N, P and SRA values in holm oak roots. Marginal and conditional R2for all the response variables, except C:N ratio, were similar. Variance explained by conditional R2for the C:N ratio response almost doubled the marginal R2 (Table 1). Relationships Between Soil and ECM Variables and Transfer of Trace Elements to Roots Transfer of TEs from soil to root seemed to be mainly driven by biotic factors: species composition PCoA2 and melanization CWM (Table 2). The soil-to-root transfer of As (RS As) was related to species composition PCoA2 and ECM melanization (Table 2). A high abundance of Hebeloma cavipes species and high melanin content seemed to be associated to a high As transfer to roots. A negative relationship with soil As was also found (Estimate = −0.0004). The soil-to-root transfer of Fe (RS Fe) was positively explained by soil NO3and melanization (Figure 3E and Table 2). Species composition PCoA1 was also positively related (Estimate = 0.026), meaning that in those soils where rare species were abundant, Fe transfer to roots was higher. Rhizomorphs formation was negatively related (Estimate = −0.01). The soil-to-root transfer of Mn (RS Mn) was best explained by abiotic variables, namely soil Ca and soil Mn, which were negatively associated with this transfer (Table 2). Species composition PCoA2 showed an individual negative effect on Mn transfer (Estimate = −0.06), therefore the abundance of Hebeloma cavipes species in the soil was found to be positively related to Mn transfer from soils to roots. The soil-to-root transfer of Zn (RS Zn) was negatively affected by species composition PCoA2 (Figure 3F and Table 2). Therefore, as previously found for As, a higher abundance of Hebeloma cavipes species increased the soil-to-root transfer of Zn. In this case, the significant effect of melanization was negative (Estimate = −0.01), opposite to the As and Fe transfers. Soil Ca (Estimate = −0.00009) showed an individual negative effect on soil-to-root transfer of Zn. Marginal and conditional R2showed similar percentage of variances for Fe, Mn, and Zn transfer but transfer of As was more explained by the site random effect (conditional R2) than the biotic fixed effects (Table 2). Relations of Soil Factors With Translocation of Trace Elements to Leaves Translocation of TEs from roots to leaves were explained by different abiotic and biotic factors (Table 2), depending on the element. Translocation of As was not significantly explained by any individual abiotic or biotic factor. Due to the nonsignificant fixed effect of the model for the response variable translocation factor of As, a covariate Cu transfer was studied. Translocation factor of Cu was highly explained by soil Cu and species composition PCoA2 (R2LMMm = 0.55; R2LMMc = 0.76). Soil Cu contamination showed a significant negative effect on the Cu translocation (p<0.001) while PCoA2 showed a significant positive effect (p= 0.013), therefore Cu translocation was favored on Thelephora terrestris dominated soils. Iron translocation from roots to leaves was only significantly explained by the biotic emanating hyphae, showing a positive relationship (Figure 3D and Table 2). Translocation factor of Mn was significantly related to NAG enzyme activity and species composition PCoA1, being this last variable the most explicative, showing a positive effect (Estimate = 8.54) (Table 2). Translocation factor of Zn was only significantly explained by soil Zn, however, Zn showed a negative effect on Zn transfer (Table 2). High differences between marginal and conditional R2 variance were found for all TE translocation factors, except for Fe (Table 2). Relations of Soil Factors With Leaf Traits Nutritional status of holm oak leaves were, in general, highly affected by soil P and ECM fungal species composition (PCoA1 and PCoA2 factors; Table 3). Leaf C was highly explained by a combination of abiotic and biotic factors (soil As and species composition PCoA1 factors, Table 3). Both predictor variables showed a strong positive relationship with leaf C. Emanating hyphae was also found to influence leaf C content, but negatively (Estimate = −0.09), when univariate relationships were analyzed. Leaf N was significantly influenced by soil P (Table 3) which explained a high proportion of variance of leaf N. No biotic factor was identified as significant for leaf N. Leaf C:N ratio was also highly explained by soil P but a negative effect was observed, in coherence with leaf N effects. A biotic effect was significantly found in relation to species composition PCoA1. The positive effect (Estimate = 2.04) of PCoA1 on this ratio showed consistency with model effects on leaf C. In summary, the results from these models showed a higher leaf C content and, therefore a higher C:N ratio, in those sites with particular abiotic characteristics (high As contamination and low soil P), and associated with certain biotic features: low abundance of Hebeloma cavipes and Thelephora terrestris, and low emanating hyphae. Leaf P was best explained by species composition PCoA2 alone, which had a positive effect on this response variable (Figure 3C and Table 3). Soil P and available Mn had also a significant influence on leaf P, according to abiotic univariate models; soil P had a positive effect (Estimate = 0.0007) while Mn availability showed a negative effect (Estimate = −0.0005) on leaf P. Leaf N:P ratio was best predicted by species composition PCoA2 alone; the negative effect between PCoA2 and this the ratio corroborates the previous leaf P results. No abiotic variables were found to have a significant effect on leaf N:P ratio. To summarize, a higher leaf P and a lower N:P ratio were found Frontiers in Plant Science | www.frontiersin.org 8November 2018 | Volume 9 | Article 1682 fpls-09-01682 November 16, 2018 Time: 20:55 # 9 Gil-Martínez et al. Ectomycorrhizal Fungi Mediate Plant–Soil Interactions TABLE 2 | Univariate and multivariate linear mixed models showing significant soil and ECM fungi fixed effects for each of the soil-to-root (RS) transfer and translocation factor (TF) and model explained variance. Response variable Individual effects of soil factors Individual effects of ECM fungal factors Combining significant effects into the best predictive model Nutrients and EA Trace elements PCoA axis 1 and 2 Fungal traits CWM Linear mixed effect models Variance Transfer Variable pVariable pVariable pVariable pModel SE t p R2LMMm R2LMMc RS As – – Total As 0.036 PCoA2 0.008 Melanization 0.050 RS As = −0.32 −0.04 PCoA2 PCoA2 0.01 PCoA2 −3.54 PCoA2 0.001 0.27 0.63 +0.006 Melanization Mel 0.002 Mel 2.90 Mel 0.007 RS Fe NO3<0.001 – – PCoA1 0.002 Melanization 0.004 RS Fe = −0.12 + 0.009 NO3 NO30.002 NO34.39 NO3<0.001 0.48 0.51 Rhizomorph 0.021 +0.003 Melanization Mel 0.0009 Mel 2.99 Mel 0.005 RS Mn Ca 0.003 Total Mn 0.004 PCoA2 0.039 – – RS Mn = 0.61 − 0.00005 Ca −0.0003 Total Mn Ca 0.00001 Total Mn 0.00009 Ca −4.09 Total Mn −3.83 Ca <0.001 Total Mn <0.001 0.44 0.44 RS Zn Ca <0.001 – – PCoA2 <0.001 Melanization 0.041 RS Zn = 0.50 −0.16 PCoA2 0.03 −5.29 <0.001 0.46 0.51 TF Fe – – – – – – Hyphae 0.026 TF Fe = 0.009 + 0.0002 Hyphae 0.0001 2.32 0.026 0.13 0.13 TF Mn NAG <0.001 – – PCoA1 0.004 – – TF Mn = 3.62 +1.82 PCoA1 0.58 3.15 0.004 0.16 0.75 TF Zn – – Total Zn 0.007 – – – – TF Zn = 0.98 −0.002 Total Zn 0.0007 −2.88 0.007 0.22 0.61 NAG, N-acetyl-glucosaminidase; CWM, community weighted mean; EA, enzyme activity; SE, standard error; R2LMMm, marginal variance; R2LMMc, conditional variance. Frontiers in Plant Science | www.frontiersin.org 9November 2018 | Volume 9 | Article 1682