# Habitat- and soil-related drivers of the root-associated fungal community of Quercus suber in the Northern Moroccan forest

Fatima Zahra Maghnia, Younes Abbas, Frédéric Mahé, Benaissa Kerdouh, Estelle Tournier, Mohamed Ouadji, Pierre Tisseyre, Yves Prin, Naima El Ghachtouli, Salah Eddine Bakkali Yakhlef, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Fatima Zahra Maghnia, Younes Abbas, Frédéric Mahé, Benaissa Kerdouh, Estelle Tournier, et al. Habitat- and soil-related drivers of the root-associated fungal community of Quercus suber in the Northern Moroccan forest. PLoS ONE, 2017, 12 (11), 17 p. 10.1371/journal.pone.0187758. [http://dx.doi.org/10.1371/journal.pone.0187758](http://dx.doi.org/10.1371/journal.pone.0187758)

HAL Id: hal-02627010

https://hal.inrae.fr/hal-02627010

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

RESEARCH ARTICLE

# Habitat- and soil-related drivers of the root-associated fungal community of *Quercus suber* in the Northern Moroccan forest

Fatima Zahra Maghnia¹,²,³,⁴, Younes Abbas⁵, Frédéric Mahé²,³, Benaissa Kerdouh¹, Estelle Tournier²,³, Mohamed Ouadji¹, Pierre Tisseyre³,⁴, Yves Prin²,³, Naïma El Ghachtouli⁶, Salah Eddine Bakkali Yakhlef⁷, Robin Duponnois³,⁴, Hervé Sanguin²,³*

1 Forestry research center, Rabat, Morocco, 2 CIRAD, UMR LSTM, Montpellier, France, 3 LSTM, Univ Montpellier, CIRAD, IRD, INRA, Montpellier SupAgro, Montpellier, France, 4 IRD, UMR LSTM, Montpellier, France, 5 Polyvalent Laboratory, Multidisciplinary Faculty, University of Sultan Moulay Slimane, Béni Mellal, Morocco, 6 Laboratory of Microbial Biotechnology, Faculty of Sciences and Technology, University of Sidi Mohamed Ben Abdellah, Fez, Morocco, 7 Directorat for Higher Education, Training and Research, Ministry of Agriculture and Fisheries, Morocco

* herve.sanguin@cirad.fr

<figure><img src="image_4.png" /></figure>

## OPEN ACCESS

**Citation:** Maghnia FZ, Abbas Y, Mahé F, Kerdouh B, Tournier E, Ouadji M, et al. (2017) Habitat- and soil-related drivers of the root-associated fungal community of *Quercus suber* in the Northern Moroccan forest. PLoS ONE 12(11): e0187758.
https://doi.org/10.1371/journal.pone.0187758

**Editor:** Ricardo Aroca, Estacion Experimental del Zaidin, SPAIN

**Received:** January 27, 2017

**Accepted:** October 25, 2017

**Published:** November 20, 2017

**Copyright:** © 2017 Maghnia et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All codes and representative sequences of OTUs can be found in HTML format (S1 File) and raw data are available under the BioPproject ID PRJNA378471 (https://www.ncbi.nlm.nih.gov/bioproject).

**Funding:** This study was supported by the "Programme de Recherche Agronomique pour le Développement" (PRAD 04/10, Gestion et valorisation de la strate arbustive ectotrophe pour la conservation et la productivité durable des subéraies) and the Cirad PhD program "Doctorant

## Abstract

Soil fungi associated with plant roots, notably ectomycorrhizal (EcM) fungi, are central in above- and below-ground interactions in Mediterranean forests. They are a key component in soil nutrient cycling and plant productivity. Yet, major disturbances of Mediterranean forests, particularly in the Southern Mediterranean basin, are observed due to the greater human pressures and climate changes. These disturbances highly impact forest cover, soil properties and consequently the root-associated fungal communities. The implementation of efficient conservation strategies of Mediterranean forests is thus closely tied to our understanding of root-associated fungal biodiversity and environmental rules driving its diversity and structure. In our study, the root-associated fungal community of *Q. suber* was analyzed using high-throughput sequencing across three major Moroccan cork oak habitats. Significant differences in root-associated fungal community structures of *Q. suber* were observed among Moroccan cork oak habitats (Maâmora, Benslimane, Chefchaoun) subjected to different human disturbance levels (high to low disturbances, respectively). The fungal community structure changes correlated with a wide range of soil properties, notably with pH, C: N ratio (*P* = 0.0002), and available phosphorus levels (*P* = 0.0001). More than 90 below-ground fungal indicators (*P* < 0.01)—either of a type of habitat and/or a soil property—were revealed. The results shed light on the ecological significance of ubiquitous ectomycorrhiza (*Tomentella*, *Russula*, *Cenococcum*), and putative sclerotia-associated/ericoid mycorrhizal fungal taxa (*Cladophialophora*, *Oidiodendron*) in the Moroccan cork oak forest, and their intraspecific variability regarding their response to land use and soil characteristics.
du Sud". We thank Morris Traduction (http://andrewmorris.fr/) for English language editing.

**Competing interests:** The authors have declared that no competing interests exist.

# Introduction

Soil fungi are one of the most diverse groups of organisms on Earth [1] colonizing a wide range of ecological niches [2] and playing a central role in major ecological and biogeochemical processes, notably in forests [3,4]. Yet, forest ecosystems are highly threatened by global changes [5,6]. Recent extreme droughts in Southern Mediterranean forests have increased tree mortality [5]. This negative impact has been exacerbated by increasing human pressures such as deforestation and overharvesting [7]. In this context, the Mediterranean cork oak (Quercus suber) forest (covering more than 20,000 square kilometers) is particularly at risk due to bark exploitation, intensive agro-silvo-pastoral management, and oak decline [7]. The degradation of cork oak forests strongly affects soil fungi, notably ectomycorrhizal (EcM) fungi [8–10], a key component of nutrient cycling and plant productivity [11]. Conservation strategies of these forests are dependent on our understanding of the soil fungal biodiversity, notably the fungal components directly associated with plant roots. The conservation of soil fungi themselves has seldom been considered compared to other Eukaryotes [12] and initiatives are relatively recent (http://iucn.ekoo.se/en/iucn/welcome). For instance, in December 2016, 34 fungal species mostly distributed in Europe, were present in the IUCN Red List of threatened species (http://www.iucnredlist.org). The diversity and dynamic of soil and root-associated fungal communities in cork oak ecosystems has been mainly investigated in the Northern Mediterranean basin [8–10,13,14]. Hitherto, estimations of soil fungal diversity in Southern cork oak forests have been almost exclusively based on fungal sporocarp surveys [15], but this approach has been shown to provide a partial view of true below-ground soil fungal diversity, even for the EcM community [16]. To address the challenges of cork oak forest conservation, soil fungal diversity in Southern cork oak ecosystems must be more extensively explored and environmental drivers affecting fungal diversity identified. The influence of environmental factors has been shown to vary with the type of habitats, the spatial scale considered and the fungal taxa analyzed [2,17,18]. For instance, variations in phosphate (P) and nitrogen (N) soil content have been suggested as important drivers of intraspecific variability of Pisolithus spp. in the Maâmora habitat, the largest cork oak forest in Morocco [19].

The main goals of the present study are (i) to characterize in depth the molecular diversity of the root-associated fungal community of *Q. suber* (fungal endophytes, EcM fungi and fungi in adherent soil); (ii) to determine the main environmental factors (habitats, soil characteristics) driving this fungal community; and (iii) to identify fungal indicators associated with a type of habitat and/or a soil property.

# Material and methods

## Study site and sampling

The study was conducted in three habitats of the Moroccan cork oak forest, located in the Moroccan Northern Mountains known as “Chefchaoun” (35°15’5.14”N 005°30’6.68”W, 1534 m elevation), and in the lowland bordering the Atlantic Ocean (North-West of Morocco) known as “Maâmora” (34°17’06.186”N 6°28’30.792”W, 27 m elevation) and “Benslimane” (33°41’9.85””N, 6°54’7.26”W; 326 m elevation). The three habitats are under a Mediterranean-type climate characterized by hot and dry summers, mild and wet winters, and mean annual rainfall levels of 453 mm (Benslimane), 570 mm (Maâmora) and 880 mm (Chefchaoun). They are characterized by an abundant understory, notably *Cistus salviifolius*, *Lavandula stoechas*, and *Thymeleae lythroides* for Maâmora, and *Arbutus unedo* and *Pistacia lentiscus* for Benslimane, and *Erica arborea* and *Arbutus unedo* for Chefchaoun. Twenty seven cork oak trees were sampled between February and June 2013. The sampling design was based on the
selection of three plots per forest spaced 100 meters apart, each composed of three trees 20–30 meters away from any other. Roots with soil were sampled under the crown of each tree and stored at +4°C. The roots were rinsed under tap water to remove the non-adherent soil and observed under a binocular microscope to select root zones rich in ectomycorrhizal (EcM) fungi, dried and stored at -20°C. The product of this sampling process resulted in a fungal community including fungal endophytes, EcM fungi and fungi in adherent soil. This fungal community is hereafter named “root-associated fungal community”.

Soil physico-chemical parameters were measured at the LAMA Laboratory (Dakar, Senegal): pH, total nitrogen (N), total carbon (C), Carbon:Nitrogen ratio (C:N ratio), total and available phosphate (P), K+, Mg2+, Na+, Cation exchange capacity (CEC). The habitats are managed by the High Commission for Water and Forests and Combatting Desertification. The permissions for root and soil sampling were provided by the Forestry Research Centre of Rabat (Morocco).

## DNA extraction, ITS amplification and Illumina MiSeq sequencing

For each of the 27 cork oak tree samples, all root pieces and the adherent soil were subjected to liquid nitrogen grinding for homogeneization. The total DNA was extracted from a sub-sample (70–80 mg) using a FastPrep-24 homogenizer (MP biomedicals Europe, Illkirch, France) and the FastDNA® SPIN kit (MP biomedicals Europe) according to the manufacturer's instructions. The purity of DNA extracts was improved by adding 20–30 mg Polyvinylpolypyrrolidon (PVPP) during the first step of DNA extraction in order to avoid the presence of PCR inhibitors.

The Internal transcribed spacer ITS1 of the nuclear ribosomal RNA was amplified using the primers ITS1FI2 (5' -GAACCWGCGGARGGATCA-3') and ITS2 (5' -GCTGCGTTCTTCATCGATGC-3') [20]. The amplification reaction was performed in a final volume of 25 µl with the primers ITS1FI2 and ITS2 (0.6 µM each), 2 µl of DNA extract, 200 µM of each dNTP, 200 ng/ml BSA, GoTaq® DNA Polymerase (2 units) and 1X Green GoTaq® Reaction Buffer (Promega, Charbonnieres, France), with the following cycling conditions: 95°C for 15 min; 30 cycles of 95°C for 30 s, 58°C for 30 s, 72°C for 30 s; a final elongation step at 72°C for 5 min. To increase richness recovery and to limit PCR biases, three PCR replicates per sample were pooled and purified using an illustra GFX PCR DNA and Gel Band Purification Kit (GE Healthcare Life Sciences, Velizy-Villacoublay, France) following manufacturer's guidelines. All amplicon products were subjected to paired-end Illumina MiSeq sequencing (2x300 bp) by Molecular Research LP (MR DNA, TX, USA).

## Bioinformatic data processing

Paired Illumina MiSeq reads were assembled with vsearch v1.11.1 [21] using the command fastq_mergepairs and the option fastq_allowmergestagger. Demultiplexing and primer clipping was performed with cutadapt v1.9 [22], enforcing a full-length match for sample tags and allowing a 2/3-length partial match for forward and reverse primers. Only reads containing both primers were retained. For each trimmed read, the expected error was estimated with vsearch's command fastq_filter and the option eeut. Each sample was then dereplicated, i.e. strictly identical reads were merged, using vsearch's command derep_fulllength, and converted to FASTA format.

To prepare the clustering, the samples were pooled and submitted to another round of dereplication with vsearch. Files containing expected error estimations were also dereplicated to retain only the lowest expected error for each unique sequence. Clustering was performed with swarm v2.1.8 [23], using a local threshold of one difference and the fastidious option.
Molecular operational taxonomic unit (OTU) representative sequences were then searched for chimeras with vsearch's command uchime_denovo [24]. In parallel, representative sequences received taxonomical assignments using the stampa pipeline (https://github.com/frederic-mahe/stampa) and a custom version of the fungal reference database UNITE v7 (https://unite.ee/; [25]). In brief, the stampa pipeline requires the reference sequences to be trimmed with cutadapt, using the same primers as those used for the amplification of the environmental sequences. Using vsearch's exact global pairwise comparisons, each environmental sequence is compared to all reference sequences and is assigned to the closest hit. In case of a tie, the environmental sequence is assigned to the last-common ancestor of the co-best hits. The abbreviation "cf." is used throughout the text for interpretation at the species level because of potential taxonomic biases relative to length and variability of the amplified region (ITS1 region). As the UNITE database contains only fungal sequences, 20 sequences were added to our custom reference database to identify the plants also amplified by our ITS primers (see S1 File).

Clustering results, expected error values, taxonomic assignments and chimera detection
results were used to build a raw OTU table (script available in S1 File). Up to that point, reads
that could not be merged, reads without tags or primers, reads shorter than 32 nucleotides and
reads with uncalled bases ("N") were eliminated. To create the "cleaned" OTU table, additional
filters were applied to retain: non-chimeric OTUs, OTUs with an expected error divided by
length below 0.0002, OTUs containing more than 3 reads or seen in 2 samples, OTUs assigned
to plant or fungal taxa with at least 80% similarity or containing more than 10,000 reads. All
codes and representative sequences of OTUs can be found in HTML format (S1 File) and raw
data are available under the BioPproject ID PRJNA378471 (https://www.ncbi.nlm.nih.gov/
bioproject).

Statistics

Diversity (Shannon, inverse Simpson [1/D]), richness (number of MOTUs, Chao1) and even-
ness (Pielou) indexes were estimated using R [26] and the R package vegan version 2.4-3 [27],
and differences among cork oak habitats were assessed by non-parametric permutational mul-
tivariate analysis of variance (PERMANOVA), as implemented in the *perm.anova()* function
from the R package RVAideMemoire version 0.9-65 [28].

Fungal community membership was assessed using venn diagram analysis with the R pack-
age VennDiagram version 1.6.17 [29]. The differences in fungal community structure among
the three habitats were displayed with nonmetric multi-dimensional scaling (NMDS) imple-
mented in the *metaMDS()* function. Significance in fungal community structure variation was
also assessed using PERMANOVA in the *adonis()* function. Multivariate dispersion was esti-
mated using the *betadisper()* function and *permutest()* as it can affect PERMANOVA results.
Soil parameters were fitted to the NMDS using the *envfit()* function (9,999 permutations). Cor-
relation among soil parameters was assessed using the Pearson correlation coefficient, as
implemented in the *cor.test()* function. All functions are available in the R package vegan.
Table transformations in R were performed with the tidyverse packages version 1.1.1 [30], and
plots were visualized with the packages ggplot2 version 2.2.1 [31] for NMDS and ggtern ver-
sion 2.2.1 [32] for ternary plots.

The presence of fungal indicator species of a specific type of habitat (Maâmora, Benslimane, Chefchaoun) was determined using the corrected Pearson's phi coefficient of association (r.g; 9,999 permutations) implemented in the *multipatt()* function from the R package indicspecies [33]. Fungal indicator species with respect to soil properties (pH, C:N ratio, and available P) were determined using the indicator value (IndVal) index, as implemented in indicspecies' *multipatt()*. Two different probabilities were calculated, *i.e*. A (specificity), representing the
probability of a site to be defined by a given soil property, given that the species have been detected, and B (sensibility) representing the probability of finding the species in different sites characterized by a given soil property. We considered as valid indicators the OTUs showing both A (specificity) and B (fidelity) superior to 0.8 and 0.6 respectively, as recommended in [18].

# Results

## Composition of cork oak root-associated fungal community

A filtered dataset of 792,931 sequences was obtained from a raw dataset of 1,129,145 trimmed- and paired-sequences (see material and methods for more details). Sixty percent of sequences were affiliated to ITS sequences from plants (Streptophyta). Overall, a dataset of 315,597 fungal sequences (27 samples) was rarefied down to 3,447 sequences per sample to improve the robustness of fungal community comparison among the three habitats [34].

Analysis of taxonomic fungal community composition in the Moroccan cork oak forest detected 1,768 OTUs belonging to 4 known fungal phyla, 39 orders, 78 families, and 127 genera (S1 Table). Two percent of OTUs were uncharacterized at the phylum level (about 0.2% of sequences), and 43% at the genus level (about 30% of sequences). The root-associated fungal community of *Q. suber* was mostly composed of Basidiomycota and Ascomycota, 60% and 39%, respectively. However, a higher number of genera was found for Ascomycota compared to Basidiomycota, 94 and 31, respectively. The ten most abundant fungal orders (89% of sequences) were Thelephorales (17,084 sequences), Russulales (16,275), Agaricales (14,017), Hysteriales (12,815), Sebacinales (5,902), Chaetothyriales (5,622), Heliotales (4,796), Hypocreales (2,575), Capnodiales (1,895), and Eurotiales (1,690) (S1 Table). The ten most abundant fungal genera (55% of sequences) were *Cenococcum* (12,459 sequences), *Tomentella* (11,609), *Russula* (7,019), *Inocybe* (4,744), *Cortinarius* (4,568), *Cladophialophora* (3,046), *Hygrophorus* (2,643), *Lactarius* (2,196), *Sebacina* (1,667), and *Cryptosporiopsis* (1,317) (S1 Table).

## Habitat-related root-associated fungal indicators

No habitat-related impact (Maâmora, Benslimane, Chefchaoun) was detected on root-associated fungal community richness, diversity and equitability (Table 1), except for Chao1 estimation between Benslimane and Chefchaoun (*P* < 0.015).

By contrast, analysis of fungal OTU distribution in ternary plots revealed strong fungal community patterns (abundance and membership) among habitats (Fig 1). Russulales (40% of sequences) was the most abundant fungal order in Maâmora, Thelephorales (41%) in Benslimane, and Hysteriales (20% of sequences) and Agaricales (19%) in Chefchaoun. At least

<table><caption>Table 1. Alpha diversity of root-associated fungal communities of *Q. suber* in Moroccan cork oak forests.</caption><thead><tr><th></th><th>MOTUs number (Richness)</th><th>Chao1 (Richness)</th><th>Shannon's index (diversity)</th><th>Inverse Simpson's index (diversity)</th><th>Pielou's index (Evenness)</th></tr></thead><tbody><tr><td>Maâmora</td><td>211±30 a<sup>2</sup></td><td>476±67 ab</td><td>2.74±0.22 a</td><td>6.17±1.50 a</td><td>0.31±0.10 a</td></tr><tr><td>Benslimane</td><td>183±69 a</td><td>365±119 a</td><td>2.34±0.90 a</td><td>6.27±4.55 a</td><td>0.27±0.06 a</td></tr><tr><td>Chefchaoun</td><td>239±47 a</td><td>509±101 b</td><td>2.97±0.55 a</td><td>8.42±6.52 a</td><td>0.34±0.02 a</td></tr><tr><td>Forest</td><td>ns<sup>1</sup></td><td>*</td><td>ns</td><td>ns</td><td>ns</td></tr></tbody></table>

<footnote n="1">¹ Statistics were performed using PERMANOVA (Habitat type as factor). Data in the same column followed by the same letter are not significantly different according to PERMANOVA (*P* < 0.05).</footnote>

<footnote n="2">² *P* < 0.05; ‘ns’ non significant.</footnote>

<footnote n="2">² Values indicate mean ± standard deviation.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig 1. Distribution of fungal OTUs among habitats (Maâmora, Benslimane, Chefchaoun) in a ternary plot.** A circle represents an OTU, with the circle size corresponding to the abundance (number of sequences) of the OTU. The color indicates the taxonomic affiliation of OTUs at order level. The category "others" corresponds to all OTUs, except the one belonging to the eight most abundant orders).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0187758.g001

65% of OTUs were habitat-specific (318 OTUs for Maâmora, 378 for Benslimane, 461 for Chefchaoun), but these accounted for only 6% of total sequences (S1 Fig). The most abundant taxonomic orders for each habitat (> 50% sequences) composed of habitat-specific OTUs were Trechisporales, Agaricales, Russulales in Maâmora, Thelephorales in Benslimane, and Agaricales, Cantharellales, Heliotales, Sebacinales in Chefchaoun. Meanwhile, 17% of OTUs were shared among the three habitats, representing 84% of total sequences (S1 Fig). However, strong differences were observed for shared OTUs among habitats in terms of relative abundance. For instance, few sequences belonging to Agaricales OTUs were present in Maâmora (). A high variability in OTU abundance and distribution was also observed among habitats for the most abundant genera (Fig 2). Among the three most abundant fungal genera, *Cenococcum* and
<figure><img src="image_6.png" /><figcaption>**Fig 2.** Distribution of fungal OTUs among habitats (Maâmora, Benslimane, Chefchaoun) at genus level in ternary plots. The three most abundant genera are indicated, (A) *Cenococcum* (45 OTUs, 12,516 sequences), (B) *Tomentella* (140 OTUs; 11,605 sequences), (C) *Russula* (64 OTUs, 7,027 sequences). A circle represents an OTU, with the circle size corresponding to the abundance (number of sequences) of the OTU.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0187758.g002

*Russula* were predominant in Maâmora and Chefchaoun, whereas *Tomentella* was more predominant in Benslimane and Chefchaoun (Fig 2).

Permutational multivariate analysis of variance (PERMANOVA) based on the Bray-Curtis dissimilarity matrix confirmed that cork-oak habitats are significant drivers of root-associated fungal community structures (Table 2). Non-significant multivariate dispersion of data was detected, emphasizing the robustness of PERMANOVA results (Table 2). Pairwise PERMANOVA comparisons of cork oak habitats showed significant differences between root-associated fungal communities (S2 Table). Indicator species analysis revealed 40 OTUs significantly associated with Maâmora, 21 to Benslimane and 49 to Chefchaoun (S3 Table). The most significant OTUs (r.g > 0.5; *P* < 0.01) (Fig 3) belonged almost exclusively to Ascomycota, and were affiliated to Archaeorhizomycetaceae, Herpotrichiellaceae, Mycosphaerellaceae, Myxotrichaceae, Trichocomaceae, and unidentified Sordariomycetes for Maâmora; unidentified Dothideomycetes and Sordariomycetes for Benslimane; Dermateaceae, Geoglossaceae, Gloniaceae, Herpotrichiellaceae, Dothideomycetes, and Eurotiomycetes for Chefchaoun (S1 Table).

## Soil-related drivers of root-associated fungal community structures

Soil characteristics in the Moroccan cork oak forest were investigated (S4 Table) in order to identify the main abiotic soil parameters driving root-associated fungal community structures. Correlation analysis between each of the soil parameters (S5 Table) revealed strong positive

<table><caption>**Table 2.** Impact of forest habitat on cork oak root-associated fungal community structures.</caption><thead><tr><th>Model/Factors</th><th>Df</th><th>SS</th><th>MS</th><th>F. Model</th><th>R2</th><th>P-value<sup>1</sup></th></tr></thead><tbody><tr><td>PERMANOVA</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Habitat</td><td>2</td><td>1.701</td><td>0.85356</td><td>2.3665</td><td>0.16473</td><td>0.001</td></tr><tr><td>Residuals</td><td>24</td><td>8.6563</td><td>0.36068</td><td></td><td>0.83527</td><td></td></tr><tr><td>Total</td><td>26</td><td>10.3634</td><td></td><td></td><td>1</td><td></td></tr><tr><td>BETADISPER</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Habitat</td><td>2</td><td>0.033575</td><td>0.0167876</td><td>2.9099</td><td></td><td>0.077</td></tr><tr><td>Residuals</td><td>24</td><td>0.138459</td><td>0.0057691</td><td></td><td></td><td></td></tr></tbody></table>

Df, degrees of freedom; SS, sum of squares; MS, mean squares; R2, R-squared; F.Model, F value by permutation.

https://doi.org/10.1371/journal.pone.0187758.t002

<footnote n="1">¹ The significance of multivariate analysis of variance and dispersion was assessed with a permutational test (Iterations = 999).</footnote>
<table><caption>**Fig 3.** Major fungal indicators of habitat and soil property. Only fungal indicators with a significant association with at least one ecological condition (black box) are shown, i.e. fungal OTUs associated with a habitat with r.g > 0.5 and *P* < 0.01 and/or a soil property with A > 0.8, B > 0.6 and *P* < 0.01. For certain major fungal indicators, a low but significant association with a habitat (r.g > 0.5 and 0.01 < *P* < 0.05) and/or a soil property (A > 0.8, B > 0.6 and 0.01 < *P* < 0.05) was also revealed (gray box). See S3 and S6 Tables for details. The taxonomic assignment is provided until genus level, except for those with an unidentified genus (where higher taxonomic level is indicated).</caption><thead><tr><th rowspan="2">Taxonomic assignment</th><th rowspan="2"></th><th rowspan="2">Fungal OTU (name)</th><th rowspan="2">Fungal abundance (number of sequences)</th><th colspan="5">Habitat</th><th colspan="2">Soil property</th></tr><tr><th>Maâmora</th><th>Benslimane</th><th>Chefchaoun</th><th>Soil acidity</th><th>SOM decomposition rate</th><th>Soil available P content</th></tr></thead><tbody><tr><td rowspan="3">Archaeorhizomycetes</td><td>Archaeorhizomyces</td><td>236</td><td>360</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Mycosphaerellaceae</td><td>493</td><td>271</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Capnodiales</td><td>75</td><td>1270</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="10">Dothideomycetes</td><td></td><td>292</td><td>92</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>32</td><td>1117</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>56</td><td>151</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>362</td><td>79</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>471</td><td>91</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>1008</td><td>18</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Lophiostoma</td><td>45</td><td>540</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Saccharicola</td><td>28</td><td>294</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Dothideomycetes</td><td>78</td><td>170</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>1190</td><td>13</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="12">Eurotiomycetes</td><td></td><td>1137</td><td>17</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>33</td><td>2014</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>112</td><td>221</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>246</td><td>118</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>1194</td><td>10</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>1241</td><td>11</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>3363</td><td>5</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Herpotrichiellaceae</td><td>1188</td><td>19</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Chaetothyriales</td><td>492</td><td>31</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Aspergillus</td><td>82</td><td>314</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Paecilomyces</td><td>19070</td><td>3</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">Geoglossomycetes</td><td>Sarcoleotia</td><td>610</td><td>24</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Trichoglossum</td><td>1664</td><td>81</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">Lecanoromycetes</td><td>Degelia</td><td>703</td><td>28</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Cryptosporiopsis</td><td>95</td><td>548</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="4">Leotiomycetes</td><td></td><td>254</td><td>23</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Phialocephala</td><td>545</td><td>3</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>453</td><td>54</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Oldiodendron</td><td>105</td><td>472</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Pezizomycetes</td><td>Pezizales sp</td><td>551</td><td>5</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="10">Sordariomycetes</td><td></td><td>69</td><td>311</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>759</td><td>27</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Sordariales</td><td>928</td><td>5</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Sordariomycetes</td><td>1937</td><td>9</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>366</td><td>54</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>407</td><td>57</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Hygrophorus</td><td>51</td><td>2592</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Agaricales</td><td>376</td><td>10</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Inocybe</td><td>315</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Lactarius</td><td>60</td><td>1260</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Russula</td><td>72</td><td>37</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Russulaceae</td><td>16</td><td>3230</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="8">Agaricomycetes</td><td></td><td>705</td><td>197</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>4242</td><td>18</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>Sebacina</td><td>572</td><td>56</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>17</td><td>5448</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>780</td><td>70</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>66</td><td>347</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>85</td><td>736</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>230</td><td>24</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>358</td><td>393</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="2">unidentified</td><td>Thelephoraceae sp</td><td>38</td><td>1024</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Fungi sp</td><td>445</td><td>16</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0187758.g003

correlations among almost all of them, with the notable exceptions of (i) the C:N ratio, which is not correlated with pH, total N, total C, available P, Na⁺ and K⁺; and (ii) the available P, which is not correlated with pH and total N. Available P, Mg²⁺, K⁺ and Cation exchange capacity were the soil parameters the most strongly correlated (*r* > 0.75) with others.

NMDS ordination showed the habitat effect on root-associated fungal community structure (Fig 4), which was confirmed by a significant correlation (R² = 0.5512, *P* = 0.0001). All ten soil parameters fitted in NMDS were also significantly correlated with the root-associated fungal
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Nonmetric multidimensional scaling analysis of *Q. suber* root-associated fungal community structures (OTU level) in Moroccan cork oak habitats (Maâmora, Benslimane and Chefchaoun) and soil parameter fitting. The stress of the ordination is 0.1584. All significant soil parameters are shown by arrows (the length is proportional to the strength of the correlation. Factor fitting (Habitat) is indicated below the soil parameter. The most significant soil parameters and factor (*P* < 0.001) are indicated in bold in the table.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0187758.g004

<table><thead><tr><th>Soil parameter</th><th>NMDS1</th><th>NMDS2</th><th>R<sup>2</sup></th><th>Pr(&gt;r)</th></tr></thead><tbody><tr><td>pH</td><td>0.94034</td><td>0.34024</td><td>0.5238</td><td>0.0002</td></tr><tr><td>Total N</td><td>0.75175</td><td>0.65945</td><td>0.2645</td><td>0.0248</td></tr><tr><td>Total C</td><td>0.46470</td><td>0.88547</td><td>0.2676</td><td>0.0217</td></tr><tr><td>C:N ratio</td><td>-0.60758</td><td>0.79426</td><td>0.5592</td><td>0.0002</td></tr><tr><td>Total P</td><td>0.25916</td><td>0.96584</td><td>0.4598</td><td>0.0003</td></tr><tr><td>Available P</td><td>-0.20510</td><td>0.97874</td><td>0.5559</td><td>0.0001</td></tr><tr><td>Mg<sup>2+</sup></td><td>0.07739</td><td>0.99700</td><td>0.4554</td><td>0.0003</td></tr><tr><td>Na<sup>+</sup></td><td>0.22786</td><td>0.97369</td><td>0.3761</td><td>0.0030</td></tr><tr><td>K<sup>+</sup></td><td>0.35212</td><td>0.93596</td><td>0.3847</td><td>0.0020</td></tr><tr><td>CEC</td><td>0.23164</td><td>0.97280</td><td>0.3734</td><td>0.0025</td></tr><tr><td>Factor</td><td>-</td><td>-</td><td>R<sup>2</sup></td><td>Pr(&gt;r)</td></tr><tr><td>Habitat</td><td>-</td><td>-</td><td>0.5512</td><td>0.0001</td></tr></tbody></table>

community structure (Fig 4). The Maâmora habitat appeared as the most different in terms of soil characteristics compared to Benslimane and Chefchaoun habitats. The three most significant drivers of root-associated fungal community structures were pH, C:N ratio (*P* = 0.0002), and available P (*P* = 0.0001). Analysis of fungal indicators with respect to these three main soil properties revealed the significant association (A [specificity] > 0.8; B [sensitivity] > 0.6; *P* < 0.05) of 27 OTUs with soil acidity levels (pH), 28 OTUs with soil organic decomposition rates (C:N ratio), and 60 OTUs with soil available P contents (S6 Table). As observed for the fungal indicators of habitat, a higher diversity of Ascomycota than Basidiomycota (only Agaricomycetes) was significantly associated with soil properties. Among the most significant OTUs (*P* < 0.01), seven OTUs affiliated to *Cenococcum*, *Cladophialophora*, *Aspergillus*, *Crytoposporiopsis*, and *Ilyonectria*, were also indicators of a habitat (Fig 3).

## Discussion

We conducted an in-depth investigation of the molecular diversity of the root-associated fungal community of *Q. suber* in the Moroccan cork oak forest to decipher the main environmental (habitats, soil parameters) drivers of this fungal community and to characterize fungal indicators associated with a type of habitat and/or a soil property.

## Ectomycorrhiza but also ericoid mycorrhiza are predominant in Moroccan cork oak root-associated fungal community

EcM fungi (notably *Cenococcum*, *Tomentella*, *Russula*, *Inocybe*, *Cortinarius*) were the most abundant fungal groups associated with the roots, as previously observed in most oak forests [13,18,35-39]. Fungi belonging to *Tomentella*, *Russula*, and *Cenococcum* were proposed as major actors of forest ecosystems under drought conditions [8], and the predominance of *Cenococcum* and *Tomentella* confirmed previous EcM surveys based on morphological identification in the Maâmora habitat [40]. *Cenococcum* is one of the most common and abundant fungi
in forest ecosystems [41,42], but over 40 different OTUs were unequally distributed in the three forest habitats, confirming reports of intraspecific diversity at local and large geographical scales [43–45]. Remarkably, potential parasites of *Cenococcum sclerotia* [46], *Cladosphialophora* and *Oidiodendron* (10% of Ascomycota sequences) were present along with *Cenococcum*. An in-depth investigation of ecological interactions between sclerotia-associated fungi and *Cenococcum* may thus be important for a better understanding of fungal assembly rules affecting forest ecosystem functioning, particularly in Mediterranean drought scenarios.

By contrast, well-known oak-associated EcM fungal taxa, i.e. *Tuber* [47] and *Pisolithus* [15], represented a minor part of the root-associated fungal community of *Q. suber* (0.4% and < 0.00001% of total sequence number, respectively). The recently-characterized *Tuber* cf. *cistophilum*, isolated from Spanish acidic argillaceous soils [48], constituted 7% of *Tuber* abundance. Similarly, *Terfezia*, a genus of truffle-like fungi that is an important edible fungus for local Moroccan inhabitants, was also detected. The main species detected was *Terfezia cf. pini*, but at low abundance (0.0005% of total sequence number). Whereas, the low abundance of *Tuber* had been previously observed in the Sardinian cork oak forest [10], *Pisolithus* had been described as predominant in different Moroccan cork oak habitats [15]. Discrepancy in EcM community composition and abundance may be due to method-dependent biases as previously highlighted in other fungal surveys [9, 10, 38, 49].

Surprisingly, ericoid mycorrhizal (ErM) fungi, *Oidiodendron* [50] and *Cryptosporiopsis* [51], were relatively abundant in cork oak roots analyzed in the current study. The presence of *Oidiodendron* members had been also evidenced in EcM roots of *Quercus ilex* [50]. Two principal ecological implications of EcM-ErM interactions in forest ecosystem functioning have been proposed by [50], (i) nutrient exchanges among EcM and ErM plants through hyphal links of shared mycorrhizal fungi, and (ii) EcM plants acting as an ErM fungal reservoir for the efficient land recolonization of ErM plants. The significance of the second hypothesis is of particular importance for cork oak ecosystems, strongly affected by fire events in the Mediterranean basin [52]. It also raises the question: could other fungi, such as arbuscular mycorrhizal (AM) fungi, play a role in EcM plants? Evidence of AM fungal colonization of EcM plants and their functional role has been demonstrated [53–55]. However, AM fungi generally constitute a minor part of the total fungal community in EcM-dominated ecosystems [14]. In the current work, AM fungi were almost absent from the raw data using the ITS-based approach. A more AM specific approach based on the 18S rRNA gene could address AM fungal community structures in EcM ecosystems [56,57].

## Ecological specificity of Moroccan cork oak root-associated fungal community

A worldwide fungal survey highlighted differences between Mediterranean forests and temperate deciduous forests regarding EcM fungal communities [2]. The current results strengthened the specificity of the Mediterranean oak-associated EcM fungal community compared to Northern European oak forests. Indeed, the predominance of the three main EcM families, Thelephoraceae, Russulaceae and Gloniaceae in different Mediterranean oak forests [8,9,58] were confirmed in the current study. These accounted for almost half of total root-associated fungal abundance, whereas Cortinariaceae constituted one of the three predominant EcM families in Northern European oak forests [18,59]. The stronger predominance of Gloniaceae (*Cenococcum spp.*) in Mediterranean forests may be linked to its resistance to drought conditions [60,61] and thus also to precipitation rates [35].

At the Mediterranean scale, results also highlighted the ecological specificity of Moroccan cork oak root-associated fungal community. *Lactarius chrysorrheus*, a predominant
*Russulaceae* species in the Sardinian cork oak forest [10], widely distributed in Mediterranean oak forests [8,9,13,35,38], was not detected in the present study. Nevertheless, other known species of *Lactarius* were detected, notably *Lactarius cf. quietus*, which is considered as an oak specialist [18]. However, *L. quietus* distribution appears more related to Northern Europe oak ecosystems [18,59] than Mediterranean oak ecosystems [8,38]. A wide range of *Russula* (*Russula cf. heterophylla*, *Russula cf. olivobrunnea*, *Russula cf. violeipes*) and *Inocybe* (*Inocybe cf. bresado-lae*, *Inocybe cf. glabripes*, *Inocybe cf. posterula*, *Inocybe cf. subporospora*) species were detected specifically in Morocco and not in other Mediterranean cork oak forests [9,10,13,15,62]. By contrast, *Russula cf. decipiens* was also described in Corsica [62], *Russula cf. foetens* in Portugal [8], and *Russula cf. odorata* in Corsica and Sardinia [10,14,62], and only *Inocybe cf. asterospora* was described in the Portuguese cork oak forest [8]. Surprisingly, *Laccaria* and *Tricholoma*, two EcM Agaricales frequently observed at above and below-ground levels in Mediterranean oak ecosystems [9,10,13,14,38,39] were not detected in the current work. These two fungal genera had been however observed in the Moroccan cork oak forest at above-ground level, reflecting, as for *Pisolithus*, the discrepancy between above- and below-ground fungal abundance. Remarkably, *Archaeorhizomyces*, a genus of ubiquitous root endophyte fungi [63], was relatively abundant (> 3% of Ascomycota) in Moroccan cork oak roots. The presence of this genus in Mediterranean ecosystem surveys is rarely described and few data are available regarding its ecological role, but a continuum from root endophytic to free-living saprophytic life strategies has been proposed [64].

## Land use and soil parameters are major drivers of Moroccan cork oak root-associated fungal community

More than 90 fungal OTUs were significantly (*P* < 0.01) associated with habitat and/or soil properties, notably EcM and ErM fungi, underlining the value of fungi as below-ground indicators of forest status and environmental conditions [18,65]. The high number of habitat-associated fungal OTUs for Maâmaora and Chefchaoun suggested their strong ecological specificity. Differences in human pressure levels may explain a part of these specificities since Maâmaora is more disturbed (intensive cork exploitation, overgrazing, aging of tree population) than Chefchaoun (lower human pressure, natural tree regeneration). Remarkably, an intraspecific *Cenococcum* ecological preference related to land use and succession stages (association with Chefchaoun, low disturbances), and not only to drought resistance, was shown in the Moroccan cork oak forest. However, *Cenococcum* is usually described as a multi-stage fungus due to its worldwide predominance and stability in forest ecosystems [41,42,58]. Certain *Russula* OTUs might also be seen as indicators of disturbance levels (association with Maâmaora, high disturbances). *Russula* is generally considered as a late-successional fungal genus [62,66,67], which may explain its high sensitivity to land use as observed in the current study and in the Portuguese cork oak forest (Azul et al., 2009), but also its reliability as an indicator of ecosystems characterized by the aging of cork oak populations (Maâmaora being an open forest with low regeneration rate). *Russula* was also shown as highly sensitive to soil nitrogen levels, seasons, and cork oak decline [18,40,58,68], and whereas no decline symptoms were observed, N total content was significantly correlated to differences in the root-associated fungal community structure among the three cork oak habitats. The potential status of *Tomentella* as an indicator of land use [40] was unclear in the current study. *Tomentella* was associated with the habitat characterized by the lowest rainfall level, but its high sensitivity to seasonal fluctuations [58] makes it difficult to relate its resistance to environmental stresses such as drought.

ErM fungi were also among the most abundant fungal groups significantly associated with habitat disturbances, *Oidiodendron maius* [69,70] and *Cryptosporiopsis* [51], and were
respectively associated with high and low disturbances. This observation strengthens the hypothesis that EcM plants can act as reservoirs of ErM fungi [50], but also suggests that certain ErM fungal members such as *Oidiodendron* might play a role for efficient Ericaceae re-colonization of habitats characterized by strong Ericaceae layer degradation such as found in Maâmora.

A wide range of soil characteristics, whether or not related to land use, were also described as strongly affecting soil fungal communities in forest ecosystems [2,17,18]. In the current work, pH, C:N ratio, and available P appeared as the most significant soil parameters correlated to differences in the root-associated fungal community structure. Soil pH is described as the most influential driver of fungal OTU richness and composition [2]. While a narrow soil pH range was observed in the current study for the Moroccan cork oak habitats, a wide range of fungal OTUs associated with high and low acidity were identified. The impact of phosphorus and C:N ratio was shown as more taxa-specific [2]. Some EcM fungi notably able to adopt a saprophytic lifestyle depending on environmental conditions [71] were associated with C and P cycling in the current study. *Tomentella* species were notably associated with a high SOM decomposition rate and low available soil P content. By contrast, the EcM *Cenococcum* appeared strongly associated with high P soils.

The present study sheds new light on the below-ground interactions taking place in the Moroccan cork oak forest, a severely threatened ecosystem. A highly diverse soil fungi community, both in Ascomycota and Basidiomycota, strongly structured by habitat type (related to disturbance levels) and the soil properties has been observed. However, the characterization of robust fungal indicators of environmental conditions remains a critical point depending of the taxonomic level considered. Indeed, high intraspecific variability (OTU level) were observed for certain fungal genera or species, reflecting not only a taxonomic variability but a probable functional variability since their association with environmental conditions differed. This has been pointed out for several EcM fungi (e.g *Tomentella, Russula*) [66], and its extent at the whole community level remains a major issue in fungal-based functional ecology [72]. At the biogeographic level, the results suggest a Southern Mediterranean fungal pattern, reinforcing the need to investigate the Southern Mediterranean fungal diversity more extensively. In addition, the strong geographical structure of *Q. suber* in the Mediterranean basin [73] may also affect fungal diversity associated with cork oak, probably contributing to the differences between fungal community structures in Morocco, as compared to Sardinia and Corsica. Indeed, intra-specific plant host diversity has been described as a soil fungal diversity driver [74] but remains poorly investigated. The success of future conservation strategies thus depends on joint initiatives on a Mediterranean scale, associating plant population geneticists, botanists, naturalists, microbiologists and pedologists to picture in a broader framework of the tremendous complexity of plant-fungal and fungal-fungal assembly rules in forest ecosystems [16,62,75,76].

**Supporting information**

S1 File. Data and code snippets used to produce the results presented in this study.
(ZIP)

S1 Table. Taxonomic affiliation of root-associated fungal OTUs of *Q. suber* in the Moroccan cork oak forest (Maâmora, Benslimane, Chefchaoun).
(XLSX)

S2 Table. Variation in *Q. Suber* root-associated fungal community structures among Moroccan cork oak habitats.
(DOCX)
**S3 Table.** Fungal indicators of Moroccan cork oak habitats (Maâmora, Benslimane and Chefchaoun).
(DOCX)

**S4 Table.** Soil parameters in the Moroccan cork oak forest.
(XLSX)

**S5 Table.** Correlation coefficient between soil physical-chemical parameters among Moroccan cork oak habitats.
(XLSX)

**S6 Table.** Fungal OTU indicators with respect to soil properties (pH, C:N ratio, available P).
(DOCX)

**S1 Fig.** Venn diagram representing shared and unique EcM-related fungal OTUs among habitats (Maâmora, Benslimane, Chefchaoun).
(TIFF)

**Acknowledgments**

This study was supported by the "Programme de Recherche Agronomique pour le Développement" (PRAD 04/10, Gestion et valorisation de la strate arbustive ectotrophe pour la conservation et la productivité durable des subéraies) and the Cirad PhD program "Doctorant du Sud".
We thank Morris Traduction (http://andrewmorris.fr/) for English language editing.

**Author Contributions**

**Conceptualization:** Younes Abbas, Yves Prin, Salah Eddine Bakkali Yakhlef, Robin Duponnois, Hervé Sanguin.

**Data curation:** Fatima Zahra Maghnia, Frédéric Mahé, Benaissa Kerdouh, Estelle Tournier, Mohamed Ouadji, Pierre Tisseyre, Hervé Sanguin.

**Formal analysis:** Fatima Zahra Maghnia, Frédéric Mahé, Hervé Sanguin.

**Funding acquisition:** Yves Prin, Robin Duponnois, Hervé Sanguin.

**Investigation:** Fatima Zahra Maghnia, Younes Abbas, Yves Prin, Robin Duponnois, Hervé Sanguin.

**Methodology:** Fatima Zahra Maghnia, Younes Abbas, Frédéric Mahé, Benaissa Kerdouh, Estelle Tournier, Mohamed Ouadji, Pierre Tisseyre, Hervé Sanguin.

**Project administration:** Naïma El Ghachtouli, Robin Duponnois, Hervé Sanguin.

**Supervision:** Younes Abbas, Yves Prin, Naïma El Ghachtouli, Salah Eddine Bakkali Yakhlef, Robin Duponnois, Hervé Sanguin.

**Writing – original draft:** Fatima Zahra Maghnia, Frédéric Mahé, Hervé Sanguin.

**Writing – review & editing:** Fatima Zahra Maghnia, Younes Abbas, Frédéric Mahé, Yves Prin, Naïma El Ghachtouli, Salah Eddine Bakkali Yakhlef, Robin Duponnois, Hervé Sanguin.

**References**

1. Blackwell M. The Fungi: 1, 2, 3... 5.1 Million Species? Am J Bot. 2011; 98: 426-438. https://doi.org/10.3732/ajb.1000298 PMID: 21613136
2. Tedersoo L, Bahram M, Polme S, Koljalg U, Yorou NS, Wijesundera R, et al. Global diversity and geog- raphy of soil fungi. Science. 2014; 346: 1256688-1256688. https://doi.org/10.1126/science.1256688 PMID: 25430773

3. Courty PE, Buée M, Diedhiou AG, Frey-Klett P, Le Tacon F, Rineau F, et al. The role of ectomycorrhizal communities in forest ecosystem processes: New perspectives and emerging concepts. Soil Biol Bio- chem. 2010; 42: 679-698.

4. Gessner MO, Swan CM, Dang CK, McKie BG, Bardgett RD, Wall DH, et al. Diversity meets decomposi- tion. Trends Ecol Evol. 2010; 25: 372-80. https://doi.org/10.1016/j.tree.2010.01.010 PMID: 20189677

5. Allen CD, Macalady AK, Chenchouni H, Bachelet D, McDowell N, Vennetier M, et al. A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests. For Ecol Manag. 2010; 259: 660-684. https://doi.org/10.1016/j.foreco.2009.09.001

6. Sala OE, Chapin FS 3rd, Armesto JJ, Berlow E, Bloomfield J, Dirzo R, et al. Global biodiversity scenar- ios for the year 2100. Science. 2000; 287: 1770-4. PMID: 10710299

7. Gauquelin T, Michon G, Joffre R, Duponnois R, Génin D, Fady B, et al. Mediterranean forests, land use and climate change: a social-ecological perspective. Reg Environ Change. 2016; https://doi.org/10. 1007/s10113-016-0994-3

8. Azul AM, Sousa JP, Agerer R, Martin MP, Freitas H. Land use practices and ectomycorrhizal fungal communities from oak woodlands dominated by Quercus suber L. considering drought scenarios. Mycorrhiza. 2010; 20: 73-88. https://doi.org/10.1007/s00572-009-0261-2 PMID: 19575241

9. Azul AM, Castro P, Sousa JP, Freitas H. Diversity and fruiting patterns of ectomycorrhizal and saprobic fungi as indicators of land-use severity in managed woodlands dominated by Quercus suber —a case study from southern Portugal. Can J For Res. 2009; 39: 2404-2417. https://doi.org/10.1139/X09-148

10. Lancellotti E, Franceschini A. Studies on the ectomycorrhizal community in a declining Quercus suber L. stand. Mycorrhiza. 2013; 23: 533-542. https://doi.org/10.1007/s00572-013-0493-z PMID: 23503869

11. Smith SE, Read DJ. Mycorrhizal symbiosis. 3. ed., Repr. Amsterdam: Elsevier/Acad. Press; 2009.

12. Heilmann-Clausen J, Barron ES, Boddy L, Dahlberg A, Griffith GW, Nordén J, et al. A Fungal perspec- tive on conservation biology: Fungi and Conservation Biology. Conserv Biol. 2015; 29: 61-68. https:// doi.org/10.1111/cobi.12388 PMID: 25185751

13. Barrico L, Rodríguez-Echeverría S, Freitas H. Diversity of soil basidiomycete communities associated with Quercus suber L. in Portuguese montados. Eur J Soil Biol. 2010; 46: 280-287. https://doi.org/10. 1016/j.ejsobi.2010.05.001

14. Orgiazzi A, Lumini E, Nilsson RH, Girlanda M, Vizzini A, Bonfante P, et al. Unravelling soil fungal com- munities from different mediterranean land-use backgrounds. PLoS ONE. 2012; 7: e34847. https://doi. org/10.1371/journal.pone.0034847 PMID: 22536336

15. Bakkali Yakhlef SB, Kerdouh B, Mousain D, Ducousso M, Duponnois R, Abourouh M. Phylogenetic diversity of Moroccan cork oak woodlands fungi. Biotechnol Agron Soc Environ. 2009; 13: 521-528.

16. Baptista P, Reis F, Pereira E, Tavares RM, Santos PM, Richard F, et al. Soil DNA pyrosequencing and fruitbody surveys reveal contrasting diversity for various fungal ecological guilds in chestnut orchards: Fungal diversity by NGS and fruiting surveys. Environ Microbiol Rep. 2015; 7: 946-954. https://doi.org/ 10.1111/1758-2229.12336 PMID: 26391727

17. Creamer RE, Hannula SE, Leeuwen JPV, Stone D, Rutgers M, Schmelz RM, et al. Ecological network analysis reveals the inter-connection between soil biodiversity and ecosystem function as affected by land use across Europe. Appl Soil Ecol. 2016; 97: 112-124. https://doi.org/10.1016/j.apsoil.2015.08. 006

18. Suz LM, Barsoum N, Benham S, Dietrich H-P, Fetzer KD, Fischer R, et al. Environmental drivers of ectomycorrhizal communities in Europe's temperate oak forests. Mol Ecol. 2014; 23: 5628-5644. https://doi.org/10.1111/mec.12947 PMID: 25277863

19. Bakkali Yakhlef SE, Abourouh M, Ducousso M, Duponnois R, Delaruelle C, Mousain D. Intraspecific variability of Pisolithus spp. as a response to changes in soil characteristics in a Moroccan cork oak plantation. Mycology. 2011; 2: 283-290. https://doi.org/10.1080/21501203.2011.592155

20. Schmidt P, Balint M, Greshake B, Bandow C, Rombke J, Schmitt I. Illumina metabarcoding of a soil fun- gal community. Soil Biol Biochem. 2013; 65: 128-132.

21. Rognes T, Flouri T, Nichols B, Quince C, Mahé F. VSEARCH: a versatile open source tool for metage- nomics. PeerJ. 2016; 4: e2584. https://doi.org/10.7717/peerj.2584 PMID: 27781170

22. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.jour- nal. 2011; 17: 10. https://doi.org/10.14806/ej.17.1.200

23. Mahé F, Rognes T, Quince C, de Vargas C, Dunthorn M. Swarm: robust and fast clustering method for amplicon-based studies. PeerJ. 2014; 2: e593. https://doi.org/10.7717/peerj.593 PMID: 25276506
24. Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R. UCHIME improves sensitivity and speed of chimera detection. Bioinformatics. 2011; 27: 2194–200. https://doi.org/10.1093/bioinformatics/btr381
PMID: 21700674

25. Köljalg U, Nilsson RH, Abarenkov K, Tedersoo L, Taylor AFS, Bahram M, et al. Towards a unified paradigm for sequence-based identification of fungi. Mol Ecol. 2013; 22: 5271–5277. https://doi.org/10.1111/mec.12481
PMID: 24112409

26. R Core Team. R: A Language and Environment for Statistical Computing [Internet]. R Foundation for Statistical Computing, Vienna, Austria; 2017. Available: https://www.R-project.org

27. Oksanen J, Blanchet, FG, Friendly M, Kindt R, Legendre P, McGlinn D, et al. vegan: Community Ecology Package. R package version 2.4–1 [Internet]. 2016. Available: https://CRAN.R-project.org/package=vegan

28. Hervé M. RVAideMemoire: Diverse Basic Statistical and Graphical Functions. R package version 0.9–65 [Internet]. 2017. Available: https://CRAN.R-project.org/package=RVAideMemoire

29. Chen H. VennDiagram: Generate High-Resolution Venn and Euler Plots. R package version 1.6.17 [Internet]. 2016. Available: https://CRAN.R-project.org/package=VennDiagram

30. Wickham H. tidyverse: Easily Install and Load "Tidyverse" Packages. R package version 1.1.1. [Internet]. 2017. Available: https://CRAN.R-project.org/package=tidyverse

31. Wickham H. Ggplot2: elegant graphics for data analysis. New York: Springer; 2009.

32. Hamilton N. ggtern: An Extension to "ggplot2", for the Creation of Ternary Diagrams. R package version 2.2.1 [Internet]. 2017. Available: https://CRAN.R-project.org/package=ggtern

33. De Cáceres M, Legendre P. Associations between species and groups of sites: indices and statistical inference. Ecology. 2009; 90: 3566–3574. https://doi.org/10.1890/08-1823.1 PMID: 20120823

34. Gihring TM, Green SJ, Schadt CW. Massively parallel rRNA gene sequencing exacerbates the potential for biased community diversity comparisons due to variable library sizes. Env Microbiol. 2012; 14: 285–90.

35. Aponte C, García LV, Marañón T, Gardes M. Indirect host effect on ectomycorrhizal fungi: Leaf fall and litter quality explain changes in fungal communities on the roots of co-occurring Mediterranean oaks. Soil Biol Biochem. 2010; 42: 788–796. https://doi.org/10.1016/j.soilbio.2010.01.014

36. Morris MH, Smith ME, Rizzo DM, Rejmánek M, Bledsoe CS. Contrasting ectomycorrhizal fungal communities on the roots of co-occurring oaks (*Quercus* spp.) in a California woodland. New Phytol. 2008; 178: 167–176. https://doi.org/10.1111/j.1469-8137.2007.02348.x PMID: 18194145

37. Smith ME, Douhan GW, Rizzo DM. Ectomycorrhizal community structure in a xeric *Quercus* woodland based on rDNA sequence analysis of sporocarps and pooled roots. New Phytol. 2007; 174: 847–63. https://doi.org/10.1111/j.1469-8137.2007.02040.x PMID: 17504467

38. Richard F, Millot S, Gardes M, Selosse MA. Diversity and specificity of ectomycorrhizal fungi retrieved from an old-growth Mediterranean forest dominated by *Quercus ilex*. New Phytol. 2005; 166: 1011–23. https://doi.org/10.1111/j.1469-8137.2005.01382.x PMID: 15869659

39. Richard F, Moreau P-A, Selosse M-A, Gardes M. Diversity and fruiting patterns of ectomycorrhizal and saprobic fungi in an old-growth Mediterranean forest dominated by *Quercus ilex* L. Can J Bot. 2004; 82: 1711–1729. https://doi.org/10.1139/b04-128

40. Maghnia FZ, Sanguin H, Abbas Y, Verdinelli M, Kerdouh B, El Ghachtouli N, et al. Impact du mode de gestion de la subéraie de la Maâmora (Maroc) sur la diversité des champignons ectomycorhiziens associés à *Quercus suber*. C R Biol. 2017; 340: 298–305. https://doi.org/10.1016/j.crvi.2017.04.001
PMID: 28506468

41. Obase K, Douhan GW, Matsuda Y, Smith ME. Revisiting phylogenetic diversity and cryptic species of *Cenococcum geophilum* sensu lato. Mycorrhiza. 2016; 26: 529–540. https://doi.org/10.1007/s00572-016-0690-7
PMID: 26968743

42. LoBuglio KF. Cenococcum. In: Cairney JWG, Chambers SM, editors. Ectomycorrhizal fungi key genera in profile. Berlin, Heidelberg: Springer Berlin Heidelberg; 1999. pp. 287–309. Available: http://www.springerlink.com/index/10.1007/978-3-662-06827-4_12

43. Bahram M, Pölme S, Köljalg U, Tedersoo L. A single European aspen (Populus tremula) tree individual may potentially harbour dozens of *Cenococcum geophilum* ITS genotypes and hundreds of species of ectomycorrhizal fungi: Ectomycorrhizal fungi of a single aspen tree. FEMS Microbiol Ecol. 2011; 75: 313–320. https://doi.org/10.1111/j.1574-6941.2010.01000.x
PMID: 21114502

44. Douhan GW, Huryn KL, Douhan LI. Significant diversity and potential problems associated with inferring population structure within the *Cenococcum geophilum* species complex. Mycologia. 2007; 99: 812–819. PMID: 18333505
45. Jany J-L, Garbaye J, Martin F. Cenococcum geophilum populations show a high degree of genetic diversity in beech forests. New Phytol. 2002; 154: 651–659. https://doi.org/10.1046/j.1469-8137.2002.00408.x

46. Obase K, Douhan GW, Matsuda Y, Smith ME. Culturable fungal assemblages growing within *Cenococcum sclerotia* in forest soils. FEMS Microbiol Ecol. 2014; 90: 708–717. https://doi.org/10.1111/1574-6941.12428 PMID: 25229424

47. Bonito GM, Gryganskyi AP, Trappe JM, Vilgalys R. A global meta-analysis of Tuber ITS rDNA sequences: species diversity, host associations and long-distance dispersal. Mol Ecol. 2010; 19: 4994–5008. https://doi.org/10.1111/j.1365-294X.2010.04855.x PMID: 21040049

48. Alvarado P, Moreno G, Manjon JL. Comparison between *Tuber gennadii* and *T. oligospermum* lineages reveals the existence of the new species *T. cistophilum* (Tuberaceae, Pezizales). Mycologia. 2012; 104: 894–910. https://doi.org/10.3852/11-254 PMID: 22314593

49. Richard F, Roy M, Shahin O, Sthultz C, Duchemin M, Joffre R, et al. Ectomycorrhizal communities in a Mediterranean forest ecosystem dominated by *Quercus ilex*: seasonal dynamics and response to drought in the surface organic horizon. Ann For Sci. 2011; 68: 57–68. https://doi.org/10.1007/s13595-010-0007-5

50. Bergero R, Perotto S, Girlanda M, Vidano G, Luppi AM. Ericoid mycorrhizal fungi are common root associates of a Mediterranean ectomycorrhizal plant (*Quercus ilex*). Mol Ecol. 2000; 9: 1639–1649. https://doi.org/10.1046/j.1365-294x.2000.01059.x PMID: 11050558

51. Leopold DR. Ericoid fungal diversity: Challenges and opportunities for mycorrhizal research. Fungal Ecol. 2016; https://doi.org/10.1016/j.funeco.2016.07.004

52. Schaffhauser A, Curt T, Véla E, Tatoni T. Recurrent fires and environment shape the vegetation in *Quercus suber* L. woodlands and maquis. C R Biol. 2012; 335: 424–434. https://doi.org/10.1016/j.crvi.2012.04.006 PMID: 22721564

53. Egerton-Warburton LM, Querejeta JI, Allen MF. Common mycorrhizal networks provide a potential pathway for the transfer of hydraulically lifted water between plants. J Exp Bot. 2007; 58: 1473–1483. https://doi.org/10.1093/jxb/erm009 PMID: 17350936

54. Egerton-Warburton L, Allen MF. Endo- and ectomycorrhizas in *Quercus agrifolia* Nee. (Fagaceae): patterns of root colonization and effects on seedling growth. Mycorrhiza. 2001; 11: 283–290. https://doi.org/10.1007/s005720100134 PMID: 24549348

55. Toju H, Yamamoto S, Sato H, Tanabe AS, Gilbert GS, Kadowaki K. Community composition of root-associated fungi in a *Quercus*-dominated temperate forest: "codominance" of mycorrhizal and root-endophytic fungi. Ecol Evol. 2013; 3: 1281–93. https://doi.org/10.1002/ece3.546 PMID: 23762515

56. Bagella S, Filigheddu R, Caria MC, Girlanda M, Roggero PP. Contrasting land uses in Mediterranean agro-silvo-pastoral systems generated patchy diversity patterns of vascular plants and below-ground microorganisms. C R Biol. 2014; 337: 717–724. https://doi.org/10.1016/j.crvi.2014.09.005 PMID: 25433564

57. Lumini E, Orgiazzi A, Borriello R, Bonfante P, Bianciotto V. Disclosing arbuscular mycorrhizal fungal biodiversity in soil through a land-use gradient using a pyrosequencing approach. Env Microbiol. 2010; 12: 2165–2179.

58. Corcobado T, Moreno G, Azul AM, Solla A. Seasonal variations of ectomycorrhizal communities in declining *Quercus ilex* forests: interactions with topography, tree health status and Phytophthora cinna-momi infections. Forestry. 2015; 88: 257–266. https://doi.org/10.1093/forestry/cpu056

59. Courty P-E, Franc A, Pierrat J-C, Garbaye J. Temporal changes in the ectomycorrhizal community in two soil horizons of a temperate oak forest. Appl Environ Microbiol. 2008; 74: 5792–5801. https://doi.org/10.1128/AEM.01592-08 PMID: 18658284

60. Kerner R, Delgado-Eckert E, Del Castillo E, Müller-Starck G, Peter M, Kuster B, et al. Comprehensive proteome analysis in *Cenococcum geophilum* Fr. as a tool to discover drought-related proteins. J Prote-omics. 2012; 75: 3707–3719. https://doi.org/10.1016/j.jprot.2012.04.039 PMID: 22579754

61. Peter M, Kohler A, Ohm RA, Kuo A, Krützmann J, Morin E, et al. Ectomycorrhizal ecology is imprinted in the genome of the dominant symbiotic fungus *Cenococcum geophilum*. Nat Commun. 2016; 7: 12662. https://doi.org/10.1038/ncomms12662 PMID: 27601008

62. Taudiere A, Munoz F, Lesne A, Monnet A-C, Bellanger J-M, Selosse M-A, et al. Beyond ectomycorrhizal bipartite networks: projected networks demonstrate contrasted patterns between early- and late-successional plants in Corsica. Front Plant Sci. 2015; 6. https://doi.org/10.3389/fpls.2015.00006

63. Rosling A, Cox F, Cruz-Martinez K, Ihrmark K, Grelet GA, Lindahl BD, et al. *Archaeorhizomycetes*: unearthing an ancient class of ubiquitous soil fungi. Science. 2011; 333: 876–9. https://doi.org/10.1126/science.1206958 PMID: 21836015
64. Rosling A, Timling I, Taylor DL. *Archaeorhizomycetes*: Patterns of distribution and abundance in soil. In: Horwitz BA, Mukherjee PK, Mukherjee M, Kubicek CP, editors. Genomics of Soil- and Plant-Associated Fungi. Berlin, Heidelberg: Springer Berlin Heidelberg; 2013. pp. 333–349. Available: http://link.springer.com/10.1007/978-3-642-39339-6_14

65. Cudlin P, Kieliszewska-Rokicka B, Rudawska M, Grebenc T, Alberton O, Lehto T, et al. Fine roots and ectomycorrhizas as indicators of environmental change. Plant Biosyst—Int J Deal Asp Plant Biol. 2007; 141: 406–425. https://doi.org/10.1080/11263500701626028

66. Gao C, Zhang Y, Shi N-N, Zheng Y, Chen L, Wubet T, et al. Community assembly of ectomycorrhizal fungi along a subtropical secondary forest succession. New Phytol. 2015; 205: 771–785. https://doi.org/10.1111/nph.13068 PMID: 25303438

67. Twieg BD, Durall DM, Simard SW. Ectomycorrhizal fungal succession in mixed temperate forests. New Phytol. 2007; 176: 437–447. https://doi.org/10.1111/j.1469-8137.2007.02173.x PMID: 17888121

68. Avis PG, McLaughlin DJ, Dentinger BC, Reich PB. Long-term increase in nitrogen supply alters above- and below-ground ectomycorrhizal communities and increases the dominance of *Russula* spp. in a temperate oak savanna. New Phytol. 2003; 160: 239–253. https://doi.org/10.1046/j.1469-8137.2003.00865.x

69. Kohler A, Kuo A, Nagy LG, Morin E, Barry KW, Buscot F, et al. Convergent losses of decay mechanisms and rapid turnover of symbiosis genes in mycorrhizal mutualists. Nat Genet. 2015; 47: 410–415. https://doi.org/10.1038/ng.3223 PMID: 25706625

70. Rice AV, Currah RS. *Oidiodendron maius*: saprobe in sphagnum peat, mutualist in Ericaceous roots? In: Schulz BJE, Boyle CJC, Sieber TN, editors. Microbial Root Endophytes. Berlin, Heidelberg: Springer Berlin Heidelberg; 2006. pp. 227–246. Available: http://link.springer.com/10.1007/3-540-33526-9_13

71. Churchland C, Grayston SJ. Specificity of plant-microbe interactions in the tree mycorrhizosphere biome and consequences for soil C cycling. Front Microbiol. 2014;5. https://doi.org/10.3389/fmicb.2014.00005 PMID: 24523715

72. Johnson D, Martin F, Cairney JW, Anderson IC. The importance of individuals: intraspecific diversity of mycorrhizal plants and fungi in ecosystems. New Phytol. 2012; 194: 614–28. https://doi.org/10.1111/j.1469-8137.2012.04087.x PMID: 22489902

73. Magri D, Fineschi S, Bellarosa R, Buonamici A, Sebastiani F, Schirone B, et al. The distribution of *Quercus suber* chloroplast haplotypes matches the palaeogeographical history of the western Mediterranean. Mol Ecol. 2007; 16: 5259–5266. https://doi.org/10.1111/j.1365-294X.2007.03587.x PMID: 17995923

74. Tedersoo L, Bahram M, Cajthaml T, Põlme S, Hiiesalu I, Anslan S, et al. Tree diversity and species identity effects on soil fungi, protists and animals are context dependent. ISME J. 2016; 10: 346–362. https://doi.org/10.1038/ismej.2015.116 PMID: 26172210

75. Clemmensen KE, Finlay RD, Dahlberg A, Stenlid J, Wardle DA, Lindahl BD. Carbon sequestration is related to mycorrhizal fungal community shifts during long-term succession in boreal forests. New Phytol. 2015; 205: 1525–1536. https://doi.org/10.1111/nph.13208

76. Toju H, Sato H, Yamamoto S, Kadowaki K, Tanabe AS, Yazawa S, et al. How are plant and fungal communities linked to each other in belowground ecosystems? A massively parallel pyrosequencing analysis of the association specificity of root-associated fungi and their host plants. Ecol Evol. 2013; 3: 3112–24. https://doi.org/10.1002/ece3.706 PMID: 24101998