Variable retention harvesting influences belowground plant-fungal interactions of Nothofagus pumilio seedlings in forests of southern Patagonia

Study site

This study was carried out in the Argentine portion of Tierra del Fuego Island, where working timberlands constitute c.a. 220,000 ha of monospecific lenga stands (Collado, 2001). We studied three old-growth N. pumilio stands and adjacent areas of VR timber management (with AR and DR treatments) at Los Cerros Ranch (54°18′S, 67°49′W). Primary forest stands are 5–10 ha with uneven age structure and older trees aged up to 400 years old. The PF stands presented a closed canopy with a sparse understory of relatively few woody and herbaceous plants (Lencinas et al., 2008). Managed stands with the VR system include areas of AR, intact forest “islands” (60 m radius, one per hectare), surrounded by DR. In DR areas evenly dispersed seed trees maintain a basal area of dominant trees at 10–15 m² ha⁻¹ (Martínez Pastur et al., 2009). These sites are part of the Parcelas de Ecología y Biodiversidad de Ambientes Naturales en Patagonia Austral (PEBANPA), a long-term monitoring network in southern Patagonia (Peri et al., 2016b), and have been monitored for seedling recruitment dynamics since 2004.

Seedling sampling and regeneration measurements

To characterize seedling EMF communities across timber management treatments in relation to seedling biomass, we sampled seedlings in three replicate sites of each PF, AR, and DR treatment. Logging occurred in the growing season of 2004–2005, and in 2012 we focused our sampling on young seedlings, most of which were under seven years old and established after logging. At each site we sampled 10 lenga seedlings along three randomly positioned 50 m transects in close proximity (10–50 m) to permanent PEBANPA network regeneration plots, collecting one seedling approximately every five meters, for a total of 30 seedlings per site and 90 seedlings per treatment. In all sites our sampling did not account for distance to the rooting zone of lenga trees and therefore represents the full variation of EMF-seedling dynamics. At one of PF sites we sampled an additional transect. In total, we sampled 280 seedlings across 28 transects and nine sites.

Harvested seedlings were transported back to the Agroforestry Lab at the Centro Austral de Investigaciones Científicas (CADIC-CONICET) in Ushuaia, Argentina. Shoots were separated from roots. Aboveground biomass was dried at 70 °C for 48 h and weighed. The age of the seedling was verified by counting leaf whirls (Martínez Pastur et al., 2011). Thirteen of 280 seedlings established before 2005, e.g., prior to timber harvest, across all three treatments (three in AR, seven in DR, and three in PF). We chose not to exclude these seedlings because they were relatively well distributed across treatments and comparative bi- and multi-variate analyses indicated that they did not impact statistical results.

We obtained lenga regeneration metrics including sapling (>1 year-old plants) and seedling (<1 year old plants) density in the permanent plots of the PEBANPA network at our sites (Peri et al., 2016b). Survey methods for lenga regeneration are fully described in Martínez Pastur et al. (2011). In brief, seedling and sapling density was surveyed in 18 regeneration plots of 1 m² (500 × 20 cm) in each treatment (PF, AR, DR). We also measured volumetric soil water content (VSW) in the upper 10 cm of soil using a MP406 moisture probe (ICT, Armidale, Australia) at the beginning and end of the summer.

EMF root colonization

The root system of each seedling was rinsed gently to remove soil, and then was cut into 5 cm segments. Root segments were dispersed in water in a gridded Petri dish under a dissecting microscope at 10–40× magnification. We used the gridline intercept method (Brundrett et al., 1996) to count the proportion of fine roots colonized by EMF in relation to the total number of root-gridline intersections. We randomly selected five healthy EMF-colonized root tips per seedling and preserved them in RNAlater (Ambion, Inc., Austin, Texas, USA).

Molecular analysis of EMF composition

EMF root tips were rinsed with distilled water to remove excess RNAlater and then were lyophilized. Root tips were ground in lysis buffer with a motorized sterile pestle (Kontes, Rockwood, TN, USA). DNA was extracted from each root tip using the DNEasy Plant Mini Kit (QIAGEN Inc., Valencia, CA, USA) according to the manufacturer’s instructions with minor modifications suggested in Bent & Taylor (2010). Genomic DNA was amplified in 25 µl PCR reactions containing 25 mM MgCl2, 10 mM dNTPs, 50 µM forward primer ITS1-F (CTTGGTCATTTAGAGGAAGTAA (Gardes & Bruns, 1993)), 50 µM reverse primer ITS4 (TCCTCCGCTTATTGATATGC (White et al., 1990)), 10 mg/ml bovine serum albumin, Promega Go-Taq (Sigma-Aldrich, St. Louis, MO, USA), 5X Promega Green Taq buffer, and 5 µl of genomic DNA. Reaction mixes were prepared in 0.2 ml tubes and thermocycled in an MJ Research PTC-225 thermal cycler as follows: 96 °C for 3 min, 35 cycles of 94 °C for 30 s, 52 °C for 30 s, then 72 °C for 3 min, followed by 72 °C for 10 min. PCR products were gel checked and Sanger sequenced by Functional Biosciences (Madison, WI).

Bioinformatics

Sanger sequences were trimmed, quality filtered, and grouped into operational taxonomic units (OTUs), synonymous with EMF taxa, as follows. Initial end-trimming (error rate below 0.1 in a 25 base window for both ends) and assembly of paired end reads (where available) were carried out in Codoncode Aligner 4.03. Sequences with fewer than 200 high quality base calls (phred > 20) were deleted and all remaining sequences were exported in fastq format. Next, a more stringent end-trimming was carried out using FASTQ Quality Trim program in Galaxy (Goecks, Nekrutenko & Taylor, 2010), with window sizes of 5, 20 and 40 and a phred threshold of 20. Sequences were converted to fasta format in Galaxy then aligned using Muscle (Edgar, 2004) in AliView (Larsson, 2014) to identify and manually remove primer sequences. We then used QIIME 1.9 (Caporaso et al., 2010) to cluster sequences into OTUs via open-reference with the USEARCH method (Edgar, 2010), a 95% sequence identity threshold, and the UNITE sh_refs_qiime_ver7_97_s_20.11.2016 database (Kõljalg et al., 2013). Singletons were kept at this step (min_otu_size 1). Representative sequences for each OTU were submitted to the RDP Naïve Bayesian classifier using the Warcup 2 ITS training set to estimate OTU taxonomies (Deshpande et al., 2016). We also searched GenBank using discontiguous mega-BLAST, with uncultured sequences excluded, to confirm, and in some cases improve identifications produced by RDP. OTU trophic guilds were inferred based on investigator knowledge and literature searches. All sequences are deposited in GenBank under accession numbers MH019772 –MH019959.

Following the taxonomic identification of EMF, OTUs were eliminated if they were saprotrophs, yeasts, wood decay fungi, pathogens, or if identification was uncertain. For further ordination analysis singletons (OTUs occurring once) were also eliminated to reduce the disproportionate influence of extremely rare taxa (McCune & Grace, 2001), resulting in a matrix of 49 OTUs × 228 seedlings.

Data analysis

We conducted comparisons for stand metrics including basal area, soil moisture, and the regeneration metrics including seedling and sapling density for each VR treatment (Table 1). These comparisons were tested using one-way ANOVA and further pairwise comparisons using the Tukey test of honest significant difference (p < 0.05). All statistical analyses were performed in R (R Core Team, 2016) with the exception of EMF ordination analysis and rarefaction curves described below. The seedling biomass and colonization dataset and ectomycorrhizal community composition dataset are archived with Bonanza Creek LTER: doi: 10.6073/pasta/350b701beb4ac2ebd2ef3f29893e6ed5, http://www.lter.uaf.edu/data/data-detail/id/687; doi: 10.6073/pasta/df801ec9b7f3493b6eb3bf191124e41f, http://www.lter.uaf.edu/data/data-detail/id/686. Seedling recruitment and density data are available on PEBANPA databases through CADIC-CONICET Ushuaia.

Analysis of EMF root colonization

To compare mean root colonization for each retention type we used the Kruskal–Wallis test with a post hoc Dunn Test of multiple comparisons with the Bonferroni correction with the R package dunn.test (Dinno, 2017). We used a non-parametric test because EMF colonization levels were strongly skewed due the large number of root systems colonized at high levels (mean colonization across treatments (97.2% ± 0.01 S.E.)), and regardless of data transformation the data did not meet the assumption of normality.

Analysis of EMF community structure

For EMF taxa richness and evenness, we counted the EMF taxa occurring across all seedlings within each treatment to compute a total richness estimation for PF, AR, and DR areas. At the treatment level we also summed the number of species shared between each two-way and the full three-way comparison of treatments. To compute EMF α diversity at the seedling level, we counted the number of EMF taxa associated with each seedling from our sequence dataset. We tested for significant treatment differences in seedling-level α diversity using one-way ANOVA, and tested for differences between two-way comparisons using the Tukey test of honest significant difference (p < 0.05). We used the same model structure to test for treatment differences in α diversity at the site level. Furthermore, we computed sample-based abundance rarefaction curves for each treatment to assess variation in richness across treatments. Diversity statistics for rarefaction curves were calculated using EstimateS (Colwell, 2013), where we used 100 randomization runs and randomized the samples without replacement.

To evaluate patterns in seedling- and transect-level EMF composition (Kruskal & Wish, 1978; McCune & Mefford, 2011) we performed Nonmetric Multidimensional Scaling (NMDS) ordinations. We used the Sorensen (Bray–Curtis) distance measure to represent compositional dissimilarity. We identified the solution with the lowest stress relative to the number of axes (McCune & Grace, 2002). Finally, we calculated a post hoc proportion of variance represented by each axis by calculating the squared Pearson correlation (r²) between distances in the original distance metric and pair-wise distances of objects in the ordination space (McCune & Grace, 2002).

Preliminary analysis of composition associated with individual seedlings suggested a 1-dimensional solution, which would be highly unstable and not significantly different than random. We therefore excluded all OTUs that occurred on less than four seedlings. This resulted in a final matrix of 26 OTUs × 206 seedlings. Our final matrix was transformed using Beal’s smoothing to relieve the “zero truncation problem” (Beals, 1984), which is common with heterogeneous data with a large number of zeros. We exported three ordination axes and used those as representations of EMF composition in Random Forest and liner mixed-effects model analyses described below (Fig. S1).

We combined EMF composition on individual seedlings across each transect for ordination analysis of EMF composition in relation to environmental factors. Pooling EMF composition data at the transect level allowed us to evaluate variation in composition within and across sites. We correlated this main OTU matrix with a secondary environmental matrix and expressed high correlations (r ≥ 0.20) with a biplot overlay to indicate the direction and magnitude of significantly correlated environmental variables. The secondary matrix was made up of seven environmental variables: treatment (PF, AR, DR), site (n = 9), volumetric soil water content, average seedling biomass (g), average % root colonization, and average age of seedlings across the transect. We tested whether EMF composition varied by VR treatments using Multiple-Response Permutation Procedure (MRPP) (Berry, Kvamme & Mielke, 1983) and tested whether the abundance of any taxon was an indicator of the VR treatments with Indicator Species Analysis (Dufrene & Legendre, 1997). All ordination analyses of EMF communities were performed in PC-ORD 6.0 (MJM Software Design, Gleneden Beach, OR, USA).

Analysis of relationships between seedling biomass and EMF metrics

To evaluate the importance of retention treatments, seedling age, and EMF variables (colonization, richness, and composition) in predicting seedling biomass, we used Random Forest (RF) regression trees in the randomForest package in R (Liaw & Wiener, 2015). RF is an ensemble decision tree method that uses an algorithmic approach to make predictions based on the input variables (Breiman, 2001b). The RF approach allowed us to include all our measured variables, including untransformed root colonization and biomass data because RF can accurately predict variable response with data that is not normally distributed. We used RF to calculate variable importance scores for treatment, seedling, and fungal variables when predicting seedling biomass. We built 500 regression trees using random samples of the 206 observations. At each node in the tree, two predictor variables were chosen at random from the seven explanatory variables: management treatment, seedling age, and EMF composition represented by three NMDS axes, EMF richness, and % root colonization. For each predictor, RF computes the prediction error (mean square error). The difference between the classification of observed bootstrap data and the classification of permuted data is averaged across all trees and divided by the standard error (Breiman, 2001a; Cutler et al., 2007).

We also created linear mixed-effects models using the nlme package (Pinheiro et al., 2018) to explore the fixed-effects of timber management treatment, seedling age, EMF composition represented by three NMDS axes, and EMF richness on seedling biomass with seedling nested in transect and site as a random factor. We evaluated data distributions and correlations between all predictor variables (habitat, age, root colonization, EMF richness, Axis1, Axis2, and Axis3) with a threshold spearman correlation of 0.4. None were excluded based on this criterion. Seedling biomass was log-transformed.

To further evaluate whether the presence of specific EMF taxa affected seedling biomass we correlated the transect-level species matrix with itself yielding a correlation matrix of OTUs that strongly correlated (>0.20) with each axis (Hollingsworth et al., 2013). We then used a simple ANOVA to evaluate differences in seedling biomass for seedlings that were colonized (yes/no) by each of these strongly correlated taxa.

VR impacts on seedling recruitment

Forest structure significantly varied among treatments. Basal area was greatest in PF, followed by AR, and then DR sites (F = 248.14, p < 0.01, Table 1) due to the removal of trees during timber harvesting. Harvesting treatment greatly influenced soil moisture (F = 15.66, p < 0.01). Dispersed retention sites had greater soil moisture compared to PF and AR sites (Table 1). We found marginal differences in seedling density by treatment (F = 2.73, p = 0.08, Table 1). Primary forest and AR sites had higher seedling densities compared to DR sites. Seedling densities ranged from 0 to 166 thousand ha⁻¹ in DR sites, 0 to 786 thousand ha⁻¹ in AR, and 0 to 980 thousand ha⁻¹ in PF sites. Sapling density followed the same trend, but with significant differences among treatments (F = 6.27, p < 0.01, Table 1). Sapling densities ranged from 35 to 203 thousand ha⁻¹ in DR sites, 78–620 thousand ha⁻¹ in AR, and 6–477 thousand ha⁻¹ in PF sites.

VR impacts on EMF root colonization

On average, colonization rates were high (97.2% ± 0.01 S.E.) across all treatments. EMF colonization varied from 38% to 100% among seedlings. Colonization levels varied with timber management treatment (Kruskal–Wallis χ² = 11.63, p < 0.01). On average AR sites had the highest mean colonization, which was significantly higher than colonization levels in the DR sites, but similar to those in PF sites (Table 1).

EMF taxonomic diversity

Our sequencing effort yielded 188 OTUs, 89 of which were EMF or likely EMF (Tables S1 and S2). The most abundant OTUs in the sequencing dataset were taxa in the EMF and likely EMF trophic niche with 17 EMF taxa occurring in the Phylum Ascomycota (Subphylum Pezizomycotina) and 72 in the Basidiomycota (Subphylum Agaricomycotina). All of the fungi in the Basidiomycota were in class Agaricomycetes and belonged to the orders Agaricales, Boletales, Cantharellales, Russulales, Sebacinales, and Thelephorales. The families Cortinariaceae (n = 29), Inocybaceae (n = 16), and Thelephoraceae (n = 8) were the most species rich. Despite the Basidiomycota showing higher EMF taxonomic diversity, the most abundant OTUs were ascomycetes (Table S1).

VR impacts on EMF taxa richness and occurrence

Summaries of our sequence data indicate that of the EMF taxa, 52 occurred in PF, compared to 53 in AR, and only 30 in DR. Twenty-eight taxa observed in PF stands also occurred in AR, while only 13 PF taxa occurred in DR. Seedlings from AR and DR areas shared 15 taxa. Overall, 10 of 89 taxa occurred in all treatments (Table S2). Despite some of the dominant taxa occurring across all VR timber treatments, the abundances of many of these dominants varied with VR treatment. For example, the most abundant taxon, OTU 1 Peziza depressa, occurred primarily in the DR areas. Conversely, some dominant Cortinarius, Tomentella, and Inocybe species had low abundances or were not observed in the DR sites (Table S2).

At the seedling level, EMF richness was low, ranging from one taxon to three taxa per seedling (Table 1). These richness estimates include EMF that were observed once (e.g., singletons). Seedlings in PF stands had significantly greater EMF richness than seedlings in AR and DR, and seedlings in AR showed greater richness than seedlings in DR areas (F = 15.00, p < 0.001, Table 1). These patterns of α diversity shifted slightly when we assessed the richness at the site level. Here again richness varied by treatment (F = 11.16, p = 0.01) with PF and AR sites showing greater richness than DR sites (Table 1), but no significant difference in richness was detected between AR and PF. Richness ranged from 25 to 29 taxa detected in PF sites, 24–29 taxa in AR sites, and 10–19 in DR sites. This was supported by the species rarefaction curves, which indicated that EMF species richness was lower in DR sites (Fig. 1). The shape of the curves also suggest that we did not saturate EMF diversity with our sampling in any of the treatment categories, but we likely came closer in the DR sites.

VR impacts on EMF composition

To determine the effects on EMF composition of VR management treatments (AR and DR) in comparison to PF stands we used NMDS ordination to visualize the variation in taxon abundance and occurrence within and across our study sites. Three major gradients captured ∼75% of the variance in EMF composition across transects (Axis 1 = 21.5%, Axis 2 = 18.4%, and Axis 3 = 35.3% of the variation). The NMDS ordination had a final stress of 14.45 with a final instability of 0.00046 after 278 iterations (Fig. 2). Soil moisture and seedling biomass were correlated with fungal composition along Axis 2 (seedling biomass r² = 0.47; soil moisture r² = 0.45). Additionally, there was a relatively strong visual pattern of differentiation in composition between treatments (Fig. 2). This was confirmed by MRPP analysis (T =  − 7.43, A = 0.08, p < 0.01) with the strongest differentiation between DR and AR (T =  − 6.42, A = 0.08, p < 0.01) or PF (T =  − 6.72, A = 0.10, p < 0.01). There was no significant difference in species composition between PF and AR sites (T =  − 1.57, A = 0.02, p = 0.07). Indicator species analysis showed that the abundance of OTU14, Timgrovea ferruginea, was an indicator of DR (IV = 42.70, p = 0.03), and OTU34, Cortinarius cycneus, was an indicator of AR (IV = 58.30, p < 0.01).

Relationships between seedling biomass and EMF metrics

Random Forest analysis of individual seedlings indicated that treatment (PF, AR, or DR) was the most important variable predicting seedling biomass, followed by seedling age, EMF composition, EMF richness, and % root colonization (Fig. 3). Supporting these results, the mixed-effects model analysis also indicated that seedling biomass was affected by treatment (χ² = 144.73, p < 0.01) and seedling age (χ² = 33.07, p < 0.01), while EMF composition and richness were not significantly important to explaining variation in seedling biomass (Table 2). The mixed-effects model accounted for 70.7% of the variance explaining log-transformed seedling biomass ($R_{\mathrm{marginal}}^2=0.71$, $R_{\mathrm{conditional}}^2=0.76$). Mixed-effects model parameter estimates indicate that seedlings from the AR and DR sites had greater biomass than seedlings harvested from the PF sites (PF 0.06 g ± 0.00 S.E., AR 0.11 g ± 0.01 S.E., DR 0.56 g ± 0.04 S.E., Table 2). The mean age of seedlings in PF was 4.30 years ± 0.12 S.E., slightly younger than in DR 4.58 years ± 0.16 S.E. and significantly older than in AR (3.76 years ± 0.16 S.E.) (F = 8.07, p < 0.01, where D-A diff 0.82 p < 0.01, PF-A diff 0.54 p = 0.02, PF-D diff −0.28 p = 0.36). Although seedling biomass was higher in the DR, seedling and sapling density were lower in DR sites, than AR or PF (Table 1).

Taxon-specific relationships between EMF and seedling biomass

Five OTUs were correlated with Axis 2 (Axis 2 = 18.7% of 75% total variance explained), the same axis that was correlated with seedling biomass: Peziza depressa, Clavulina cirrhata, Cenococcum geophilum, Inocybe fibrillosibrunnea, and Cortinarius amoenus. The presence of P. depressa and I. fibrillosibrunnea were significantly related to seedling biomass, while C. amoenus was marginally related to seedling biomass (Table 3). The presence of P. depressa was correlated with higher mean seedling biomass (Table 3). Peziza depressa was most abundant in DR areas (Table S2). Conversely, the presences of I. fibrillosibrunnea and C. amoenus were correlated to lower mean seedling biomass (Table 3) and these OTUs occurred in the PF and AR but not the DR sites (Table S2).

EMF composition in relation to other Southern Hemisphere studies

The most abundant EMF taxa observed in our study were ascomyceteous fungi of the order Pezizales. However, the basidiomycetes were more speciose than the ascomycetes. Within the Basidiomycota, the most species-rich order and family were the Agaricales and the Cortinariaceae, respectively. These patterns of taxonomic diversity mirror those in surveys of both root tips (Nouhra et al., 2013) and fungal fruiting bodies in the Nothofagus forests of South America (Garnica, Weiß& Oberwinkler, 2003; Truong et al., 2017). BLAST results for the five most abundant OTUs in our EMF dataset matched taxa observed in Chile and Argentina. Some of these taxa were previously thought to only occur in Australasia (e.g., OTU1, top BLAST description Ruhlandiella sp. (Truong et al., 2017)). Interestingly, the dominant families (Cortinariaceae, Inocybaceae, and Thelephoraceae) observed in our root tip survey are also abundant mycobionts in Australian forests (Tedersoo et al., 2009; Horton et al., 2017), and the most species-rich families observed at high latitudes in the Northern Hemisphere (Timling et al., 2012). Our sequence dataset supports recently described patterns of EMF diversity in South America (Truong et al., 2017) and draws attention again (Tedersoo et al., 2008) to parallels between EMF taxonomic diversity at high latitudes in both the Northern and Southern Hemispheres.

VR impacts on EMF

We found that N. pumilio seedling root systems were highly colonized by EMF (∼97% across all the studied treatments). These findings are similar to reports that both young and mature Nothofagus roots have high colonization levels of >70% (Diehl, Mazzarino & Fontenla, 2008; Fernández et al., 2015) and reinforce findings from the Northern Hemisphere where inoculum levels in logged areas remain high (Perry et al., 1982; Jones, Durall & Cairney, 2003) even when strong differences in fungal composition occur after harvest. Despite overall high colonization levels, mycorrhizal colonization was lowest in DR sites akin to observations of reduced EMF root presence in harvested New Zealand Nothofagus forests (Dickie, Richardson & Wiser, 2009). Nonetheless, even in DR sites, colonization levels were ∼95%, indicating that EMF inoculum was not severely limited for seedlings that recruited successfully.

Dispersed retention sites supported lower α diversity than PF and AR sites. Several studies from the Northern Hemisphere report reduced EMF richness with increasing distance from the forest edge or individual green trees for both seedlings that naturally established before logging (advanced regeneration) and experimentally established seedlings (Kranabetter, 1999; Hagerman, Sakakibara & Durall, 2001; Kranabetter & Friesen, 2002; Outerbridge & Trofymow, 2004; Teste, Simard & Durall, 2009). We did not directly measure the impact of seedling proximity to seed trees or the forest edge in the DR sites. Instead our results show that, on average, the DR areas have reduced EMF richness compared to the other treatments. This does not exclude the possibility that there are mycorrhizal hot-spots within DR areas due to proximity to mature tree rooting zones; however, our measurements of EMF attributes were not highly variable from seedling to seedling in the DR, which we would expect if there were great contrasts in EMF dynamics within versus beyond the rooting zone of seed trees. The low density of seed trees in the DR area may therefore account for the lower EMF richness on seedlings.

We observed that more taxa were shared between the PF and AR treatments than between PF and DR, supporting the idea that the AR patches act as refugia for species diversity (Soler et al., 2015) more effectively than the individual seed trees in the DR. Although, EMF composition was similar between PF and AR sites and differed in DR sites, there were only two taxa that showed strong gradients in abundance related to VR treatment (Timgrovea ferruginea and Cortinarius cycneus). Both of these fungi have been observed in environmental sampling in the Southern Hemisphere, but their individual ecologies are not well known (Danks, Lebel & Vernes, 2010; Johnston et al., 2010; Nouhra et al., 2013).

EMF effects on seedling performance

We detected species-level correlations between five individual EMF taxa and seedling biomass, corroborating findings from greenhouse and field studies demonstrating that the loss or gain of an EMF taxon can influence seedling performance (Jones, Durall & Cairney, 2003 and references therein). We also observed significant correlations between soil moisture, seedling biomass, and EMF community composition. Yet the causality of these relationships is unclear. Monitoring at the permanent PEBANPA network sites shows that DR areas have overall lower seedling and sapling density, which may release seedlings from competition resulting in overall larger seedlings (Martínez Pastur et al., 2009; Martínez Pastur et al., 2011). These larger seedlings may in turn select for particular fungi. Alternatively, the EMF in DR sites may promote greater seedling growth as has been observed in early-successional northern forests (Kranabetter, 2004). Outside of the intraspecific competition dynamics between seedlings and their EMF, the differences in soil moisture related to VR treatment could be a driving factor in explaining the variation in EMF composition (Boucher & Malajczuk, 1990), as well as the positive correlation between EMF and host plant biomass in sites with higher moisture (Kennedy & Peay, 2007). Yet a fourth explanation of these correlations between seedling size, recruitment density, and EMF may be the strong abiotic gradients across the VR treatments (Martínez Pastur et al., 2007; Martínez Pastur et al., 2011; Martínez Pastur et al., 2014). The high light and moisture levels in DR, compared to AR and PF sites may support relatively high biomass accrual when seedlings do recruit successfully. However, the consistently lower recruitment levels regardless of seed availability and good seedbed conditions (Martínez Pastur et al., 2011; Torres et al., 2015; Table 1) in DR sites may be due in part to reduced EMF metrics, and we therefore suggest that EMF may play a larger role in recruitment than we were able to measure with our biomass metric. Further experimentation is needed to disentangle the hierarchy of VR effects on the outcomes of seedling-EMF interactions (sensu Kranabetter, 2004; Nicholson & Jones, 2017).

VR impacts and conservation of forest biodiversity

Our study of the VR management regime can be used to infer how forests managed for multiple goals impact the diversity of EMF and plant-fungal interactions at the seedling phase in southern Nothofagus forests. Unlike other studies from our study region that have found that AR and DR sites increase the species richness of various taxa of insects, birds and plants (Soler et al., 2015; Soler et al., 2016), we found that AR sites maintained EMF richness and DR supported a lower number of taxa. The contradiction with other observations in our study system of VR effects on biodiversity likely has to do with the obligate symbiosis of EMF with N. pumilio host plants and is supported by findings from EMF studies in the Northern Hemisphere (Kranabetter, De Montigny & Ross, 2013). Within our study region, Soler et al. (2015) found that AR maintained the forest structure and microenvironment conditions found in PF stands. With variation in the physical environment more similar between PF and AR than AR and DR or PF and DR, it follows that treatment-level α diversity, evenness, and root colonization would vary more significantly between DR and the other treatments. The DR treatment negatively affected forest structure (Soler et al., 2015), resulting in a reduction in EMF host plants. Therefore, the shift in EMF composition between treatments is likely due to the presence and abundance of adult trees and may be further modulated by EMF dispersal capabilities (Peay et al., 2012; Horton, 2017), competitive abilities (Kennedy et al., 2007; Kennedy et al., 2011), and disturbance tolerance of EMF present in DR sites. A study of the impacts of dispersed green tree retention and AR showed that together these treatments maintained sporocarp production of EMF (Luoma et al., 2004), while in our study AR seem most important for the maintenance of EMF community dynamics.