Mycorrhizal Fungi Associated With Juniper and Oak Seedlings Along a Disturbance Gradient in Central Mexico

Competition for resources between arbuscular mycorrhizal (AM) and ectomycorrhizal (ECM) plants can alter belowground mycorrhizal communities, but few studies have investigated host effects on both AM and ECM communities. In Central Mexico, the AM plant Juniperus deppeana is frequently used for reforesting areas affected by soil erosion, while the surrounding native forests are dominated by ECM oak trees. Oaks are capable of associating with both AM and ECM fungi during part of their life cycle (a feature known as dual mycorrhization) but it is unclear whether junipers possess such ability. To assess how juniper planting may affect belowground fungal interactions with oaks, we investigated mycorrhizal associations in J. deppeana and Quercus rugosa seedlings along a disturbance gradient: a native oak forest, a mixed Juniperus-Quercus population in secondary vegetation and a juniper site severely degraded by mining extraction. We measured root colonization and identified fungal communities using soil and root meta-barcoding of the ITS2 rDNA region. ECM fungal community composition was strongly affected by disturbance (regardless of host), while the community composition of AM fungi was mostly host-dependent, with a higher AM fungal richness in J. deppeana. Importantly, the fungal communities associated with Q. rugosa seedlings significantly changed in the vicinity of juniper trees, while those of J. deppeana seedlings were not affected by the presence of oak trees. Even though ECM fungal richness was higher in Q. rugosa and in the native forest, we detected a variety of ECM fungi associated exclusively with J. deppeana seedlings, suggesting that this plant species may be colonized by ECM fungi. Our results indicate that J. deppeana can alter ECM native fungal communities, with implications for its use in reforestation of mixed oak forests.


INTRODUCTION
Temperate-subtropical forests occupy ca. 18% of the Mexican territory and harbor a high plant diversity, with over 160 oak species, including 109 endemics (Valencia-A, 2004). Although these forests are geographically distributed in the tropical zone, some authors consider them to be part of the temperate biome due to their cold mountain climate and the presence of dominant temperate elements, such as pines and oaks, in the vegetation (Rzedowski, 2006). Quercus rugosa Née is a dominant ectomycorrhizal (ECM) tree in pine-oak forests across the Mexican mountain ranges, with a distribution ranging from Honduras to the southern United States (Valencia-A, 2004). It possesses an efficient re-sprouting capacity that contributes to its regeneration following disturbance (Pausas et al., 2004;Cooper et al., 2018). Native forests in Mexico have been severely affected by land use changes related to agriculture, grazing, logging and mining, with an estimated loss of 662,000 ha of primary forests and a tree cover decrease of 8.1% between 2002 to 2020 (Global Forest Watch, 2014;Pérez-Suárez et al., 2014). In particular, Mexico is a major destination for mining investment in Latin America, with over 25,000 concessions granted from 2000 to 2010, equivalent to more than 13% of the national territory (Caìrdenas, 2013). These activities have strongly degraded temperate-subtropical forests in Central Mexico, calling for the urgent promotion of reforestation practices. Unfortunately, restoration plans are often conducted with plant species that are not native of the area and are, in some cases, exotic to the country (Tellez et al., 2020).
Juniperus deppeana Steud. is a drought-tolerant species native to the southwestern United States and Central Mexico (Little, 1968;Herrerías Mier and Nieto de Pascual Pola, 2020). Its capacity to colonize water-and nutrient-poor sites makes it an ideal candidate for reforestation, particularly to control soil erosion (Bush, 2008). The genus typically associates with arbuscular mycorrhizal (AM) fungi (Soudzilovskaia et al., 2020) although some studies have reported individuals forming ECM associations (Reinsvold and Reeves, 1986;Dean et al., 2015). Juniperus species have been largely re-introduced for reforesting pine-oak forests of Central Mexico that are naturally dominated by ECM associations (García-Guzmán et al., 2017;Argüelles-Moyao and Garibay-Orijel, 2018). Co-occurring AM and ECM plants are likely to spontaneously share mycorrhizal fungi in their native range (Kadowaki et al., 2018;Toju and Sato, 2018), but these aspects have scarcely been addressed for reforestation purposes. Some studies reported a loss in native plant species richness (including oaks) and alterations in nutrient acquisition following juniper encroachment Williams et al., 2013). Changes in ECM and AM plant densities may also impact plant demography and community structure through neighboring interactions mediated by beneficial and pathogenic soil fungi (Liang et al., 2020). On the other hand, access to an existing mycorrhizal network can favor the establishment and survival of emerging seedlings (Taudiere et al., 2015;O'Donnell et al., 2020). For example, AM fungi can facilitate the establishment of Juniperus seedlings in Quercus-dominated stands (Bush, 2008) as well as the survival of Q. rugosa seedlings (Olivera-Morales et al., 2011). These examples emphasize the need for more research on belowground interactions between AM and ECM communities for reforestation purposes.
Arbuscular mycorrhizal and ECM fungi often co-occur in temperate and subtropical forest soils, with contrasting strategies of nutrient acquisition that diversely affect their host performance as well as carbon (C) and nutrient cycling (Phillips et al., 2013;Williams et al., 2013;Brundrett and Tedersoo, 2018). For example, ECM fungi tend to promote conspecific and heterospecific seedling growth across co-occurring ECM plants, while AM fungi tend to negatively affect conspecific seedlings growth because of competition for root space and resources (Bennett et al., 2017;Liang et al., 2020). Nonetheless, other studies repeatedly found that AM associations had a neutral or positive effect on seedling growth (van der Heijden and Horton, 2009). Moreover, few studies have been conducted on plants that are colonized by both AM and ECM fungi during part of their life cycle, a feature known as dual mycorrhization (Queralt et al., 2019;Tedersoo and Bahram, 2019;Teste et al., 2019). The benefits and costs of dual mycorrhization for plants depend on a variety of factors including plant age and nutrient demands (Meinhardt and Gehring, 2012;Holste et al., 2017). Because of the high C cost of ECM associations, the colonization of ECM plants by AM fungi may be favored when plant growth is prioritized and nutrients are not limiting (Holste et al., 2017). Thus, seedlings of plant species that typically form ECM associations, such as oak, pine and fir, are occasionally colonized by AM fungi (Dickie et al., 2001;Wagg et al., 2008;Toju et al., 2013). Moreover, sharing AM fungi through common mycorrhizal networks likely contributes to ECM seedling survival until they are able to form more stable associations with ECM fungi (Selosse et al., 2006;O'Donnell et al., 2020). For example, laboratory experiments demonstrated that these early successional AM associations in ECM seedlings correlated with positive growth response and phosphorus (P) uptake (Lapeyrie and Chilvers, 1985;Van Der Heijden and Kuyper, 2001), but negative or neutral effects have also been reported (Egerton-Warburton and Allen, 2001;Kariman et al., 2012;Holste et al., 2017). On the other hand, primarily AM plants can gain access to new organic nitrogen (N) and P sources when additionally colonized by ECM fungi, but this comes at a higher C cost to their host (Jansa et al., 2019;Teste et al., 2019). Thus, the direction and strength of these AM-ECM interactions is context-dependent (Kadowaki et al., 2018;Teste et al., 2019;Grünfeld et al., 2020).
In Central Mexico, J. deppeana is frequently used to control soil erosion in sites that are naturally dominated by mixed oak forests. To assess how juniper planting may affect belowground fungal interactions with oaks, we investigated mycorrhizal interactions between J. deppeana (primarily AM) and Q. rugosa (primarily ECM) seedlings along a disturbance gradient: a native oak forest, a mixed site with Juniperus-Quercus populations in secondary vegetation, and a severely degraded J. deppeana site with high soil erosion from mining extraction. We hypothesized that (1) similar AM and ECM fungal species colonize J. deppeana and Q. rugosa seedlings growing side-by-side in the mixed site; (2) disturbance and the presence of J. deppeana will negatively affect ECM fungal richness, while AM fungal richness will not be affected along the gradient; (3) the fungal community structure will be strongly affected by disturbance and the presence of the co-occurring host in the mixed site.

Study Area and Vegetation
Our study took place in the mining area of Basaltex, municipality of Tepetlaoxtoc in the Patlachique mountain range, State of Mexico,Mexico (19.55 • N,98.73 • W, Figure 1). The native vegetation of the site was originally dominated by pine-oak forests that have been drastically degraded in the region due to logging, grazing and other human activities. During the last 10 years, the extraction of volcanic rocks (i.e., tezontle or andesite) used for construction work has decimated the native vegetation with few forest remains along steep canyons and mountain slopes. Severe soil erosion resulting from these mining activities also hinder forest regeneration.
A vegetation study was conducted prior to soil and root sampling to identify dominant plant components in sites with various degrees of degradation from mining and grazing. The vegetation of 32 plots (Figure 1) was investigated with a modified IFRI method (Escutia-Sánchez, 2004), consisting of surveys within circumscribed circles of 10 m radius for trees, 3 m radius for shrubs and 1 m radius for herbaceous plants. Plant species were identified according to Rzedowski and Rzedowski (2005). Vegetation structure of the tree layer was investigated in each plot by counting the number of individuals and recording the diameter at breast height (1.3 m) and the canopy cover based on foliage projection on the ground. The Importance Value Index (IVI) was used to evaluate the overall importance of each species in the plant community (Cottam and Curtis, 1956). The IVI of each taxon was calculated from the sum of its relative density, canopy cover and frequency per plot, hence varying from 0 to 300. The two focal species Q. rugosa and J. deppeana showed the highest IVI (Supplementary Table 1) in the area, with Q. rugosa dominating the tree layer over other oak and pine species in forest patches, while J. deppeana was the most frequent species in secondary vegetation resulting from logging, grazing and mining activities (López-Hernández et al., 2017).

Soil and Root Sampling
Based on the vegetation survey, we selected three study sites distanced by >1 km along a disturbance gradient (Figure 1): (i) a native site located inside the nearby forest, where Q. rugosa dominates the tree layer over other oak and pine species and J. deppeana is absent (Supplementary Table 1); (ii) a mixed site in the transition zone with native forest where J. deppeana and Q. rugosa co-occur (regeneration zone); and (iii) a disturbed site with a population of J. deppeana growing spontaneously on eroded soil in the absence of any other tree species. During the wet season in September 2018, we collected the whole root system of six J. deppeana (from mixed and disturbed sites) and six Q. rugosa (from native and mixed sites) seedlings per site, at a distance of > 10 m from each other, for a total of 24 samples. In the mixed site, seedlings were collected in the vicinity (< 3 m) of an adult tree of the other species. Samples were kept on ice and processed within 12 h. Terminal roots were carefully rinsed with 1% Tween 80 (Sigma-Aldrich) and sterilized water to reduce artifacts from the presence of hyphae or spores at the root surface, then cut into small segments. At least 10 segments of fine roots were conserved in DNeasy PowerSoil extraction buffer (QIAGEN), while another set of > 10 segments were preserved in 96% ethanol for anatomical observations. We also collected soil cores in the vicinity of half of the seedlings randomly selected in each site, for a total of 12 samples. Each sample was composed of four soil cores (5 cm diam. × 10 cm depth) taken < 50 cm from seedlings in cardinal directions. Samples were kept on ice and within 12 h, soils were sieved (2 mm mesh) and 0.25 g was stored in DNeasy PowerSoil extraction buffer (QIAGEN) while ca. 1 g of soil was dried at 60 • C for 48 h to determine soil moisture. The remaining soil from each sample was air dried for the characterization of the soil chemistry. The following edaphic parameters were measured at the Institute of Geology of the National Autonomous University of Mexico (UNAM): total C and N using an elemental analyzer CNHS/O Perkin Elmer 2400 series II; pH in 1:2.5 soil:CaCl 2 0.01 M solution; available P following Olsen (1954) and the cation-exchange capacity of calcium (Ca), aluminum (Al), hydrogen (H), magnesium (Mg) and potassium (K) by atomic absorption in 0.5 M NH 4 Cl solution using a PinAAcle 900 H (Perkin Elmer).

Root Colonization
The number of ECM root tips per cm in Q. rugosa samples was counted under a Leica MS5 stereomicroscope at 10-40 × on five root segments randomly selected per seedling. For anatomical observations of AM colonization, roots from J. deppeana preserved in 96% ethanol were depigmented with H 2 O 2 and 10% KOH at 60 • C and colored in a solution of 1% HCl with 1-2 drops of black Schaeffer ink (Vierheilig et al., 2005). Root colonization was evaluated by counting the number of AM (hyphae, arbuscules, vesicles) and endophytic (septate hyphae) structures using the "grid-line intersect" method (McGonigle et al., 1990). For each root sample, sections of five root segments were measured consecutively three times. Based on our limited sampling, we were not able to study AM colonization from oak roots nor ECM colonization in juniper samples, which will require further investigation.

Root and Soil Meta-Barcoding
DNA was extracted using the DNeasy PowerSoil kit (QIAGEN) following the manufacturer protocol. The ITS2 region of ribosomal DNA was amplified with the Platinum PCR SuperMix High Fidelity kit (Invitrogen) using fungal-specific primers ITS4ngsUNI and gITS7ngs that were designed to detect a broad spectrum of fungi (Tedersoo and Lindahl, 2016). PCR conditions were: initial denaturation at 94 • C for 1 min, followed by 35 cycles at 94 • C for 30 s, 52 • C for 30 s, 68 • C for 30 s, and a final elongation at 68 • C for 7 min. Amplicon libraries were constructed using barcodes with Nextera adapters (Illumina), normalized at equimolar concentrations, purified with Agencourt AMPure XP beads (Beckman) and sequenced using a MiSeq 300 × 2 paired-end protocol (Illumina) at the University of Arizona Genetics Core. Raw data are available at NCBI's Sequence Read Archive, Bioproject PRJNA744703.
Sequence quality filtering and clustering into operational taxonomic units (OTU) at 97% similarity was carried out in AMPtk (Palmer et al., 2018) with default parameters except for: quality-trimming -min_len 150 as the minimum length to keep a sequence, -m 10 as the minimum number of reads per OTU, -uchime_ref ITS for chimera filtering, and for filtering against index-bleed -min_reads_otu 10 as the minimum number of reads for an OTU to be retained as valid. OTU taxonomy was assigned in AMPtk by aligning reference sequences with the UNITE database (Kõljalg et al., 2005). OTUs that could not be taxonomically identified were manually checked by performing BLAST searches in UNITE (version 8.2) based on similarity thresholds for family, genus and species at >90, >95, and >97%, respectively (Tedersoo et al., 2015). The trophic mode of each OTU was assigned using FUNGuild (Nguyen et al., 2016) implemented within AMPtk. When FUNGuild failed to assign a guild as a result of taxonomic uncertainty, we assigned them manually when BLAST hits matched a UNITE "Species Hypothesis" with > 97% similarity and > 90% coverage. The OTU table was filtered based on positive and negative controls (Nguyen et al., 2015) and dubious OTUs with abundance < 0.05% of total reads per sample were filtered out.

Data Analyses
All tests were carried out in R 3.6.2 (R Core Team, 2019) using custom scripts and the phyloseq, vegan, mvabund, miceco, metacoder and ggplot2 packages (McMurdie and Holmes, 2013;Wickham, 2016;Foster et al., 2017;Oksanen et al., 2020;Wang et al., 2020;Russel, 2021). We transformed sequence read counts into relative abundances by averaging the number of reads per sample, multiplied by 10,000 and transformed to the next integer to be used as counts (Truong et al., 2019). Data were sub-divided into the following subsets prior to analyses: "All fungi" (ALL), "ECM fungi" (ECM), "AM fungi" (AM) and "saprophytic fungi" (SAP). Two samples were removed from the subset ECM and AM, respectively, because they did not contain any OTUs with these trophic modes.
Soil variables as well as root colonization and species richness (Fisher alpha index and observed richness, as the number of OTUs per sample) were compared between sites (native, mixed, disturbed), plant hosts (J. deppeana, Q. rugosa) and sample type (root, soil) using non-parametric Kruskal-Wallis tests and pairwise Mann-Whitney U-tests with Bonferroni corrections to account for multiple comparisons. We constructed a Raup-Crick dissimilarity matrix based on OTU presence/absence (after eliminating OTUs only present in one sample) and visualized fungal community composition by non-metric multidimensional scaling (NMDS). The effect of site, plant host and sample type on fungal community composition was tested with PERMANOVA using 999 random permutations (Adonis tests).
To predict how site and plant host affected individual fungal species in the community, we analyzed root samples separately by fitting generalized linear models (GLM) and performing univariate and multivariate analyses of deviance with non-normal error distribution for: (1) the 50 most frequent OTUs; (2) OTUs agglomerated per fungal families; and (3) the most species-rich ECM and AM families in our dataset (Thelephoraceae, Sebacinaceae, Russulaceae, Inocybaceae, and Glomeraceae, Supplementary Table 2). We fitted the GLM models with a binomial distribution for presence/absence data. Finally, we visualized shared and unique root mycorrhizal fungal OTUs between sites and plant hosts with a Venn diagram and we analyzed root samples using heat trees to visualize differences in OTU abundance between treatments (host/site combinations) at family level, for taxa with relative abundance higher than 0.1% of total read count per sample. Taxonomic groups were colored based on significant differences between the median proportion of reads per treatment as determined using a Wilcox rank-sum test followed by a Benjamini-Hochberg correction for multiple testing (Foster et al., 2017). The intensity of the color is relative to the log-2 ratio of difference in median proportions, analogously to up-regulated or down-regulated loci in differential gene expression analysis.

Root Colonization and Soil Variables Along the Disturbance Gradient
Arbuscular mycorrhizal colonization was found in 100% of the J. deppeana roots sampled, as well as ECM root tips in 100% of Q. rugosa roots. Endophytic root colonization was significantly lower (F = 260.80, P < 0.001) and AM colonization significantly higher (F = 143.47, P < 0.001) in the mixed site compared to the disturbed site (Table 1). Alternatively, ECM colonization in Q. rugosa did not vary significantly between the native and the mixed site. Regarding soil variables, available P was significantly higher in the native site, while C was significantly lower in the disturbed site and N significantly higher in the native site compared to the disturbed site ( Table 1). All other soil characteristics did not significantly differ across sites.
Fungal Diversity and Community Structure Between Sites and Plant Hosts 4,012,542 reads passed quality filtering and clustered into 2,047 OTUs. ECM fungi represented 42% of reads and 10% of OTUs; AM fungi 0.6 % of reads and 5% of OTUs; SAP fungi 12% of reads and 20% of OTUs; while 42% of reads and 60% of OTUs could not be assigned to a guild. Fungi that could not be taxonomically assigned to a species represented 84% of the total number of OTUs, while 61%, 47% and 18% of OTUs could not be assigned a genus, family and order respectively, emphasizing the paucity of information from rootassociated fungi from neotropical regions currently occurring in reference databases (Corrales et al., 2018).
Accumulation curves showed that our small sampling size captured the dominant fungal taxa in the communities but likely missed additional rare taxa (Supplementary Figure 1). The number of OTUs detected in roots and soil beneath Q. rugosa (1,356 OTUs) was similar to those detected beneath J. deppeana (1,409 OTUs), while OTU richness was higher in the mixed site (1,457 OTUs in total, 131 ECM, 67 AM and 321 SAP) than in the disturbed (889 OTUs in total, 52 ECM, 70 AM and 202 SAP) and in the native (772 OTUs in total, 99 ECM, 25 AM, 170 SAP) sites. These differences, however, remain non-significant when looking at total observed richness between sites, plant host or sample type. Nevertheless, ECM observed richness per sample was significantly higher in the native site, while SAP observed richness was significantly higher in the mixed site ( Table 2). ECM observed richness was significantly higher in Q. rugosa, while AM observed richness was significantly higher in J. deppeana ( Table 2). A majority of AM (51 OTUs) and ECM (89 OTUs) taxa were detected in both root and soil samples (Supplementary Figure 2) while ECM observed richness per sample was significantly higher in soil compared to roots ( Table 2). All other correlations remained non-significant and tests based on Fisher alpha indices gave highly similar results (data not shown). The global fungal community (ALL) varied significantly between sites (R 2 = 0.55, P = 0.001) and plant hosts (R 2 = 0.43, P = 0.001), but not sample type (Figure 2A and Supplementary  Table 3). When considering fungal guilds separately, we found significant differences between sites for ECM (R 2 = 0.23, P = 0.003) and SAP (R 2 = 0.31, P = 0.003), but not for AM fungi (Figures 2B-D and Supplementary Table 3). Changes in community composition between plant hosts were significant only for SAP (R 2 = 0.32, P = 0.001). When considering each plant host separately, ECM (R 2 = 0.26, P = 0.006) and SAP (R 2 = 0.42, P = 0.005) fungal community composition differed significantly between sites in Q. rugosa, while no significant differences were detected between sites for any fungal guild in J. deppeana samples.
Site significantly predicted the composition of the 50 most frequent OTUs of ALL, ECM and SAP fungi in root samples ( Table 3); it also significantly predicted ALL OTUs agglomerated by family, as well as the composition of three of the four most species-rich ECM families (i.e., Thelephoraceae, Sebacinaceae and Inocybaceae) and of the most species-rich AM family (Glomeraceae) ( Table 3). Univariate tests showed that site was a significant predictor for the distribution of an unknown species from the order of Helotiales (SAP), the family Hyaloscyphaceae (unknown fungal guild) and the genus Leohumicola (SAP). On the other hand, plant host was a significant predictor for the 50 most frequent OTUs of ALL, ECM and AM fungi in root samples (Table 3); it also significantly predicted ALL and AM OTUs agglomerated by family, as well as the distribution of the families Thelephoraceae (ECM), Sebacinaceae (ECM) and Glomeraceae (AM). Univariate tests showed that plant host was a significant predictor for the same OTUs as above in the family Hyaloscyphaceae (unknown fungal guild) and the genus Leohumicola (SAP), as well as one OTU from the subphylum Glomeromycotina (AM) and three unknown species from the phylum Ascomycota (unknown fungal guild), the order Acarosporales (unknown fungal guild) and a parasitic fungus from the genus Neonectria.

Host Specificity of Mycorrhizal Communities
When looking at fungal OTUs shared between plant hosts in root samples, we found that the majority of AM OTUs (72%, 52) were associated exclusively with J. deppeana, while 22% (16) were shared between Q. rugosa and J. deppeana and 6% (4) were only detected in Q. rugosa (Figure 3A). Most of the AM OTUs detected in Q. rugosa belonged to the Glomeraceae family, including OTUs from the genera Acaulospora, Claroideoglomus and Funneliformis found exclusively in Q. rugosa. In contrast, only 50% (56) of ECM OTUs were associated exclusively with Q. rugosa, while 24% (26) were shared between Q. rugosa and J. deppeana, and 26% (28) were detected only in J. deppeana ( Figure 3B); 9 ECM OTUs were exclusive to J. deppeana in the disturbed site. ECM OTUs detected in Q. rugosa belonged to a variety of fungal families, while ECM OTUs detected only in J. deppeana belonged to the genera Entoloma, Hygrophorus, Meliniomyces, Peziza, Ramaria, Tarzetta and Trichophaea.
Some fungal families detected in root samples displayed notable differences in relative abundance between sites and plant hosts. For example, Glomeraceae (AM) read abundance was significantly higher in J. deppeana than in Q. rugosa in general, while in the mixed site, families comprising ECM members such as Entolomataceae, Pyronemataceae and Clavariaceae were significantly higher in J. deppeana than in Q. rugosa (Figure 4). For comparisons between sites within each plant host, Tuberaceae (ECM) was significantly higher in the native site and Glomeromycetes (AM) significantly higher in the mixed site in Q. rugosa roots, while Tricholomataceae (ECM) were significantly lower in the mixed than in the disturbed site in J. deppeana roots.

DISCUSSION
Interactions between mycorrhizal fungi and their plant hosts shape below-and aboveground communities, but few studies have addressed the outcomes of ECM and AM fungi associated with the same plant host (Queralt et al., 2019;Tedersoo and Bahram, 2019;Teste et al., 2019). Likewise, the costs and benefits of associations with either ECM or AM fungi across environmental gradients remains largely unknown in the context of forest restoration (Phillips et al., 2013;Williams et al., 2013).

Ectomycorrhizal Fungal Networks as a Tool in Forest Restoration
Similarly to other studies on pine-juniper and mixed spruce forests (Hubert and Gehring, 2008;Otsing et al., 2021), disturbance and the presence of neighboring junipers FIGURE 2 | Non-metric multidimensional scaling (NMDS) ordination of fungal community composition of root and soil samples of Juniperus deppeana (in red) and Quercus rugosa (in blue) seedlings between sites (native = square, mixed = circle, disturbed = triangle) based on Raup-Crick dissimilarity distances. According to PERMANOVA, (A) the global fungal community (ALL) varied significantly between sites (R 2 = 0.55***) and plant hosts (R 2 = 0.43***), but not sample type (stress = 0.25); (B) Ectomycorrhizal (ECM) fungal communities varied significantly between sites (R 2 = 0.23**) but not plant hosts or sample type (stress = 0.25); (C) Arbuscular mycorrhizal (AM) fungal communities did not vary significantly between sites, plant hosts or sample type (stress = 0.12); (D) Saprophytic (SAP) fungal communities varied significantly between sites (R 2 = 0.31**), plant hosts (R 2 = 0.32***) and sample types (0.16683*) (stress = 0.26). Significance levels = *P ≤ 0.05, **P ≤ 0.01, and ***P ≤ 0.001. significantly affected the ECM fungal community composition of oak seedlings ( Figure 2B). On the other hand, the ECM fungal community structure associated with J. deppeana seedlings did not significantly change with disturbance nor with the presence of Q. rugosa. Dual mycorrhizal tree species have the ability to connect through shared ECM or AM mycorrhizal networks that can favor seedling establishment by facilitating access to otherwise unreachable nutrient sources (Nara, 2006;Teste et al., 2009). In a recent experiment, O'Donnell et al. (2020) suggested that juniper can facilitate oak seedling's survival by cooperating through a common network of AM fungi, until oak individuals develop more specific ECM associations. These benefits could nevertheless be context-dependent since other studies did not find any positive effects of AM colonization in oaks growing in the vicinity of junipers (Dickie et al., 2001;Bush, 2008). Our findings suggest that the establishment of oak seedlings could be sustained by a shared ECM fungal network with junipers, as evidenced by the ECM fungal OTUs shared by J. deppeana and Q. rugosa seedlings ( Figure 3B). Therefore, oak seedlings may benefit from the ECM fungal network of other plants, including junipers, to overcome nutrient limitations in highly disturbed sites where mature oak trees are absent and oak-specific ECM inoculum is not available (Teste et al., 2009;Kadowaki et al., 2018). Notwithstanding, these assumptions depend on the type and degree of disturbance and other factors, including the density of oak and juniper populations may additionally impact oak regeneration through neighboring interactions mediated by beneficial and pathogenic soil fungi (Liang et al., 2020). Further studies conducting inoculation experiments are needed to confirm the evidence of a shared ECM mycorrhizal network between oaks and junipers and their effect on seedling recruitment. Nevertheless, oaks are a source of diverse ECM inoculum and can form compatible ECM associations with the roots of other plant species, such as eucalypts (Santolamazza-Carbone et al., 2019). Our study sites are located in the transition between temperate and tropical biomes, likely favoring ECM fungal diversity and their association with a broad range of plant hosts, including junipers (Argüelles-Moyao et al., 2017). Moreover, the presence of other ECM trees species of Pinus and Quercus in the native and mixed sites can increase ECM diversity through neighboring effects (Molina and Horton, 2015). Two of the most species-rich ECM families that we detected across sites (Thelephoraceae and Sebacinaceae, Supplementary Table 2) associated with both plant hosts. Both families tend to be hostgeneralists that are abundant in subtropical and tropical oak forests (García-Guzmán et al., 2017). Some species from these families also have short-exploration type hyphae that allow them to regenerate faster compared to taxa with greater extraradical hyphae, thus recovering more easily from disturbances (Deslippe et al., 2011;Correia et al., 2021). To date a limited number of ECM fungal species, mostly from temperate and boreal regions, have been tested for reforestation purposes (Policelli et al., 2020). Given their diversity, abundance and apparent resistance to disturbance, further studies on these ECM fungal families are relevant for the restoration of less explored subtropical and tropical oak forests.

Mycorrhizal Associations Across Plant Hosts and Sites
Host preference is a key factor affecting ECM fungal diversity and community structure (Dickie, 2007;Otsing et al., 2021) and we expected that the presence of a primarily AM plant (Juniperus) would negatively affect ECM fungal diversity (Haskins and Gehring, 2004). Surprisingly, ECM fungal richness did not significantly change in the presence of juniper (Table 2), and J. deppeana seedlings in the mixed and disturbed sites associated with a diverse range of ECM fungi (Figures 3B, 4) including several OTUs that exclusively associated with juniper in the genera Entoloma, Hygrophorus, Meliniomyces, Peziza, Ramaria, Tarzetta and Trichophaea. The effect of disturbance may contribute to the presence of rare ECM taxa (Kałucka and Jagodziński, 2016). For example, Tarzetta is an early-successional ascomycete that abundantly colonizes seedlings across different forest ecosystems, particularly in sites that face environmental stress or disturbance (Correia et al., 2021). Similarly, species of Peziza are stress-tolerant taxa that can intensely colonize roots following fire disturbances (Pulido-Chavez et al., 2021), while the basidiomycete genus Entoloma is frequently found in fragmented young forest stands and is not very prone to dispersal limitation (Boeraeve, 2019). Because of its primarily AM status, studies that explore ECM fungal diversity in juniper are scarce (Bush, 2008;Hubert and Gehring, 2008). Based on our limited material, we were not able to perform comprehensive anatomical observations of the typical ECM root tip structures (mantle and Hartig net) from J. deppeana roots that are needed to consider the species as dual mycorrhizal (Teste et al., 2019). Nonetheless, some ECM taxa such as Peziza and other Pezizalean fungi have been reported to form ectendomycorrhiza with a thin or fragmented mantle and poorly developed Hartig net (Tedersoo et al., 2006). Ectendomycorrhizal morphologies have been frequently reported in tropical regions, colonizing seedlings or as early-successional species following disturbances (Fujimura et al., 2005;Salgado Salomón et al., 2013;Alvarez-Manjarrez et al., 2018). Additionally, juniper seedlings from the disturbed site were isolated from oaks, pines and any other ECM tree species (Supplementary Table 1), suggesting that these ECM associations are not the result of roots interwinding with other FIGURE 4 | Heat trees illustrating differences in OTU abundance between treatments (host/site combinations) at family level, for taxa with relative abundance higher than 0.1% of total read count per sample. The size of nodes in the large gray cladogram depicts the number of taxa identified at each taxonomic level. The smaller cladograms represent pairwise comparisons between treatments, with orange-brown nodes indicating a significantly higher abundance in the treatment displayed horizontally (based on Wilcox rank-sum tests with Benjamini-Hochberg corrections for multiple testing), while blue-green nodes indicate a significantly higher abundance in the treatment shown vertically.
Arbuscular mycorrhizal fungal communities were relatively uniformly distributed across sites, but associated more abundantly with J. deppeana than with Q. rugosa (Figures 3A,  4 and Table 2). AM fungal taxa are generally quick to colonize recently established vegetation in the vicinity of forested patches . Hence, the lack of significant change in AM fungal community composition across sites suggests that most of these taxa are generalist species able to colonize oak and juniper seedlings across various environments and disturbance gradients. Nevertheless, OTUs in the families Paraglomeraceae and Ambisporaceae were exclusively detected in juniper roots, suggesting that some level of host specialization may occur among AM fungi. It is to be noted that meta-barcoding studies on AM fungi often target the small subunit of ribosomal DNA (SSU), because of poor amplification with fungal ITS primers and the existence of accurate SSU sequence database for AM identification (Öpik et al., 2010;Tedersoo et al., 2018). Even though our approach may have overlooked some AM fungal diversity, it was demonstrated that SSU and ITS metabarcoding data respond similarly to environmental shifts (Lekberg et al., 2018).

The Effect of Disturbance on Mycorrhizal Fungi
The response of AM fungi to disturbance depends on the type and intensity of disturbance, the AM fungal species involved, as well as local environmental conditions (van der Heyde et al., 2017). For example, forest fragmentation can alter mycorrhizal colonization through changes in microsite nutrient conditions (Soudzilovskaia et al., 2015). Deforestation and mining activities gradually led to the disappearance of tree cover in the disturbed site (Figure 1) that can strongly reduce the soil inoculum of infective mycorrhizal fungi even though the plants are facing nutrient stress (Asmelash et al., 2016;Policelli et al., 2020). We found that AM colonization of J. deppeana roots was significantly lower in the disturbed site (Table 1), likely due to abiotic stresses linked to soil erosion. Accordingly, endophyte fungal colonization of J. deppeana roots was higher in the disturbed site as a consequence of abiotic stress linked to disturbance (Giauque et al., 2019;Yan et al., 2019). In contrast with studies showing that plants growing in soils with P shortage tend to be more heavily colonized by AM fungi (Bush, 2008;Kluber et al., 2012), available P was significantly lower in the disturbed site (Table 1). However, it has been demonstrated that plants do not rely on AM fungi for P acquisition in situations with very low P availability, but rather on the release of root exudates (Raven et al., 2018). In our study, soil available P was highest in the native site, in correlation with a higher ECM fungal richness (Tables 1, 2). Few studies have investigated the role of ECM fungi for P uptake (Köhler et al., 2018;Queralt et al., 2019), but evidence showed that ECM-dominated stands can have easier access to organic P than AM-dominated forests (Albornoz et al., 2016;Rosling et al., 2016). Moreover, plants that form both AM and ECM associations can allocate C preferably to the fungal symbiont that provides more benefits relative to the abundance of readily available nutrients (Albornoz et al., 2016). Therefore, P availability was likely not limiting in the native site, hence facilitating associations with a high diversity of ECM fungi for oak seedlings (Tedersoo and Bahram, 2019).
Disturbances can alter ECM fungal communities through changes in the soil C and N pools that, in turn, are likely to affect soil organic matter degradation and the rates of nutrient uptake for their host (Lindahl and Tunlid, 2015). In the native site, denser tree cover and subsequent litter accumulation contributed to soil C input (Table 1), while deforestation and soil erosion induced C losses in the disturbed site (Gómez-Guerrero and Doane, 2018). ECM fungal colonization did not significantly change between the native and the mixed sites, but we observed significant differences in both ECM fungal richness and community composition between sites ( Figure 2B and Table 2). In parallel, we also observed a decrease in SAP fungal richness in the native site (Table 2) in response to higher ECM fungal abundance and mycelial density that likely displaced SAP fungi (Fernandez and Kennedy, 2016). The lower ECM fungal richness in the disturbed site was likely due to a decrease in ECM propagules and hyphal network that can serve as inoculum for new seedlings (Kottke, 2002;Gebhardt et al., 2007). On the other hand, a dense population of mature oak and pine individuals in the native site (Supplementary Table 1) contributed to a diverse ECM fungal inoculum (Gebhardt et al., 2007;Peay et al., 2007Peay et al., , 2011. OTUs from the ECM families Clavulinaceae, Tricholomataceae and Tuberaceae were almost exclusively associated with oaks (Figure 4). Members of Tricholoma can abundantly colonize the roots of pine and oak species (Henke et al., 2015), while Tuber is commonly found in association with oaks (Morris et al., 2008). On the other hand, species of Clavulina are widely distributed in tropical latitudes in association with a variety of plant hosts (Argüelles-Moyao et al., 2017;Pérez-Pasos et al., 2019). Further studies on the ecology and host preference of these ECM fungal families would be relevant for the conservation of subtropical and tropical oak forests.

CONCLUSION
Our study suggests that J. deppeana and Q. rugosa seedlings share a common ECM fungal network and that these interactions are likely to favor oak recruitment in sites with high soil erosion. If its dual mycorrhizal status is confirmed with further anatomical studies, Juniperus may promote multifunctionality across temperate and subtropical forests through the association with generalist ECM fungi (Grünfeld et al., 2020). Further studies including plant response on a more extended sampling and over a longer timeframe are necessary to increase the robustness of our findings and assess the success of restoration strategies using dual mycorrhizal plants. Nevertheless, our findings suggest that J. deppeana is likely a good candidate for reforestation, considering its capacity to form mycorrhizal associations with a wide variety of fungi and share the mycorrhizal networks of other plant species. In particular, the potential ability of junipers to form ECM associations may enhance vegetation recovery and oak seedling establishment in highly disturbed sites where mature oak trees have been wiped out. Notwithstanding, these assumptions depend on the type and degree of disturbance and other factors, including the density of oak and juniper populations, need to be carefully evaluated in reforestation plans.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: https://www.ncbi.nlm.nih. gov/bioproject/, PRJNA744703.

AUTHOR CONTRIBUTIONS
CT designed the research. LP-L, JE-S, MDO-R, and CT conducted fieldwork. LP-L and JE-S carried out the vegetation survey. CM-G and CT performed the anatomical study. AB-C, MDO-R, and CT conducted the molecular work. AB-C and CT analyzed the data and wrote the manuscript. All authors contributed to the manuscript revision.

FUNDING
This work was funded by the Royal Botanic Gardens Victoria and a Problemas Nacionales grant #2015-218 from the Mexican Consejo Nacional de Ciencia y Tecnología (CONACYT).

ACKNOWLEDGMENTS
We would like to thank Saraí Montes Recinas and Nayeli Y. Orozco Salazar for their assistance during fieldwork and with the processing of samples. Susana Valencia Ávalos and colleagues from the Laboratory of Vascular Plants of the Faculty of Sciences at UNAM contributed with plant identifications. Lucy Mora and Christina D. Siebe Grabach conducted the soil analyses at the Institute of Geology of UNAM. Barbara B. Fransway and Jonathan R. Galina-Mehlman provided technical assistance with Illumina sequencing at the University of Arizona. Rosanne Healy revised the language and two reviewers provided valuable comments on the manuscript.