Soil Fungal Community Characteristics and Mycelial Production Across a Disturbance Gradient in Lowland Dipterocarp Rainforest in Borneo

Introduction

Tropical forests represent the most biodiverse terrestrial ecosystems on the planet (Myers et al., 2000) and provide a globally important carbon (C) sink (Pan et al., 2011). Increasingly rapid reductions in forest cover through land-use and land-conversion threatens the capacity of tropical forest to support biodiversity (Hansen et al., 2013; Powers and Jetz, 2019) and store C (Baccini et al., 2017). Significant forest degradation is occurring in Southeast Asia, with highest rates in Borneo (Bryan et al., 2013). Since the early 1970s, expansion of industrial oil palm plantation (OP) has resulted in the loss of >30% of forest cover (18.7 Mha of old-growth forest), while more than 70% of the remaining forest has been affected by selective logging for extraction of commercially valuable timber (Gaveau et al., 2014, 2016).

Ecosystem functions and biogeochemical cycles (e.g. C storage and nutrient dynamics) are underpinned by complex plant-soil interactions, mediated through reciprocal feedbacks between aboveground vegetation and soil microbial communities (Wardle et al., 2004; van der Heijden et al., 2008; Bever et al., 2010; van der Putten et al., 2013; Cortois et al., 2016). Changes in soil properties associated with forest disturbance have been shown to have significant and long-lasting effects on microbial communities and diversity across multiple biomes (Hartmann et al., 2013; Mcguire et al., 2015). However, despite the importance of plant-soil microbe interactions, understanding of relationships between above-below-ground communities alongside soil properties in tropical ecosystems remains limited (Barberán et al., 2015; Mueller et al., 2016). A small number of recent studies have highlighted the sensitivity of soil fungal communities to land-use change in Southeast Asia, including effects of selective logging and forest conversion to oil palm plantation on community structure and diversity (Kerfahi et al., 2014; Mcguire et al., 2015), although the direct impacts on functioning of fungal communities are unclear.

The mycelium produced by actively foraging soil fungi is a key component of C and mineral nutrient cycles in terrestrial ecosystems (Johnson et al., 2002; Finlay, 2008; Cairney, 2012). The extraradical mycelia of mycorrhizal fungi act as direct pathways for the reciprocal exchange of mineral nutrients from soils and photosynthetic C between plants and fungal symbionts. This mutualistic association benefits the host plant through increased surface area for the absorption of mineral nutrients from soil, in exchange for plant-derived photosynthates required for fungal growth and survival (Smith and Read, 2008; Itoo and Reshi, 2013; Chen et al., 2016). Mycorrhizal symbiosis is also known to improve plant resistance to disease and environmental perturbations including drought and exposure to harmful pollutants (Smith and Read, 2008). Plant productivity is enhanced through this mutually beneficial partnership (Lambers et al., 2008; Nasto et al., 2014; Wurzburger and Clemmensen, 2018), in turn leading to greater belowground C allocation by plants (Orwin et al., 2011) and C supply to the wider soil microbiome (Drigo et al., 2012; Nottingham et al., 2013). Recent research has evidenced nutrient transfer between plants of the same or even different species through common mycelial networks (CMNs) that physically connect two or more individuals, which can support interplant exchanges over long distances (Bever et al., 2010; Barto et al., 2012; Babikova et al., 2013; Gorzelak et al., 2015). While this process can influence plant survival and growth, it may also support establishment of young trees by providing existing infrastructure for seedlings to access (Nara, 2006; Gorzelak et al., 2015). While fungal hyphal structures themselves provide a C sink, they may also contribute to soil organic matter pools through mycelial turnover (Wallander et al., 2013) and protection of organic substrates through effects on soil aggregation (Wilson et al., 2009; Rillig et al., 2015). Although the precise role of mycelial production and activity in soil C dynamics is contentious and remains an area of active research (Zak et al., 2019), recent studies have highlighted links between mycorrhizal community structure and soil C storage. For example, dominance of ectomycorrhizal (EcM) vs. arbuscular mycorrhizal (AM) fungi may influence soil C accumulation potentially through competition with saprotrophic fungi for resources required for the decomposition of organic matter (Averill et al., 2014).

Mycelial production is regulated by a range of abiotic and biotic factors such as soil mineral nutrient availability (Hagerberg et al., 2003; Nilsson and Wallander, 2003; Potila et al., 2009; Ekblad et al., 2013), soil moisture (Majdi et al., 2008), climate (Bakker et al., 2015), phenological/seasonal changes in belowground C allocation (Ekblad et al., 2013), and direct grazing by soil mesofauna which may inhibit or stimulate mycelial production (Ek et al., 1994; Setäl et al., 1999). Plant and fungal community structures themselves are also critical for determining mycelial growth and turnover (Clemmensen et al., 2015), as shifts in vegetation and mycelial fungal communities related to forest management have been linked to patterns in mycelial production (Hagenbo et al., 2018). However, quantification of mycelial production rates has largely been confined to boreal and temperate biomes, mostly focussing on EcM-dominated systems in Scandinavia (Ekblad et al., 2013).

The lowland rainforests of Borneo are characterized by the high abundance and canopy dominance of tree species belonging to the obligate EcM-associating Dipterocarpaceae family (Whitmore, 1984; Taylor and Alexander, 2005; Brearley, 2012), with dipterocarp community assemblages potentially mediated by their mycorrhizal partners (Essene et al., 2017). Overall, the majority of tree species in these hyper-diverse forests also associate with AM fungi as in most tropical forests (Mcguire et al., 2008). Dipterocarps are directly targeted through selective logging due to their commercial value as timber (Appanah and Turnbull, 1998), making way for establishment of other tree species that are likely to be AM-forming (Mcguire et al., 2008). A recent study highlighted the effects of such tree community alterations on net primary productivity (NPP) allocation (Riutta et al., 2018), with shifts from canopy NPP to woody NPP fractions in selectively logged (SL) Bornean lowland dipterocarp rainforest. Corresponding shifts in tree functional traits have also been reported, indicating traits associated with C capture and growth to be more pronounced in SL forest compared to structural and persistence traits in old-growth (OG) forest (Both et al., 2019). These alterations to aboveground productivity patterns and plant functional traits along with modified soil properties due to logging may have strong implications for soil mycelial fungal community characteristics (i.e. the structure and diversity of actively foraging soil fungi), already indicated by shifts observed in bulk soil EcM fungal communities with SL (Mcguire et al., 2015), although these relationships and subsequent impacts on productivity to date have not been investigated. Establishment of OP involves the conversion of forest to expansive monoculture plantations of the AM-associating African species Elaeis guineensis (Phosri et al., 2010; Gaveau et al., 2014), and is associated with dramatic loss of forest-dwelling species, severely impacting on biodiversity (Fitzherbert et al., 2008; Sodhi et al., 2010; Wilcove et al., 2013). The soil microbiome is no exception, with previous studies documenting pronounced differences in both bacterial and fungal communities in OP relative to OG and SL forest, which may be relatively more similar (Lee-Cruz et al., 2013; Kerfahi et al., 2014; Mcguire et al., 2015). In particular, OP can greatly reduce richness and relative abundance of EcM fungi, which may be almost entirely absent from this land-use (Kerfahi et al., 2014; Mcguire et al., 2015). These alterations may have drastic consequences for mycelial productivity and related ecosystem functions, but impacts remain as yet unresolved.

The overarching aim of this study was to evaluate the impact of selective logging and forest conversion to oil palm plantation on mycelial fungal community characteristics and productivity, and explore the role of soil and vegetation properties as drivers of mycelial fungal community composition and productivity in tropical lowland dipterocarp rainforest. We constructed our experiment to address the following hypotheses:

H1. Mycelial fungal community attributes will differ across land-use types, with more distinct communities in OP relative to OG and SL forest (Mcguire et al., 2015). EcM relative abundance is expected to be highest in OG forest, followed by SL forest and OP (Kerfahi et al., 2014). This will be reflected in rates of mycelial production which will decline with disturbance intensity, along with removal of EcM-associating trees. Differences observed will correspond to alterations in soil and environmental properties, tree community composition and functional characteristics between land-use types resulting from disturbance.

H2. Across forest land-uses, mycelial fungal community attributes and production will be explained by a combination of soil and environmental properties, tree community composition and functional characteristics (see above).

H3. Across forest land-uses, shifts in mycelial fungal community attributes (particularly mycorrhizal) will correlate with rates of mycelial production (Clemmensen et al., 2015; Hagenbo et al., 2018).

Methods

Study Sites

This study was carried out in the Malaysian state of Sabah, northern Borneo. This region is characterized by a moist tropical climate, with average daily temperature of 27°C and annual precipitation of 2,600–2,700 mm (Walsh and Newbery, 1999; Kumagai and Porporato, 2012). Its central Southeast Asian location is considered aseasonal, with no distinct wet or dry season (Richards, 1996; Dykes, 1997; Itioka and Yamauti, 2004), although experiences irregular inter-annual dry periods which may average a total of ~1.4 months per year (Walsh and Newbery, 1999; Kumagai et al., 2005). Sampling was conducted in a total of nine 1 ha plots across OG and SL lowland dipterocarp rainforest and OP (see Supplementary Information 1 for map of sampling locations). Four plots were distributed across OG forest: two plots were situated in the Danum Valley Conservation Area (DVCA) (4.951°, 117.796° and 4.953°, 117.793°), and two in the Maliau Basin Conservation Area (MBCA) (4.747°, 116.970°, and 4.754°, 116.950°), coded as DAN-04, DAN-05, MLA-01 and MLA-02, respectively in the ForestPlots database (www.forestplots.net). These reserves have undergone little or no anthropogenic disturbance with legal protection from logging granted in 1976 and 1981, respectively (Marsh and Greer, 1992; Hazebroek et al., 2004). Four plots were located in SL forest in the Kalabakan Forest Reserve (SAF-01: 4.732°, 117.619°; SAF-02: 4.739°, 117.617°; SAF-03: 4.691°, 117.588° and SAF-04: 4.765°, 117.700°) and one in OP (4.641°, 117.452°; Benta Wawasan Sdn. Bdh.), in established research sites within the large-scale forest fragmentation study Stability of Altered Forest Ecosystems (SAFE) Project (Ewers et al., 2011). Two SL plots had been selectively logged twice and two four times, with forests undergoing timber extraction first during the mid-1970s (~113 m³ ha⁻¹ timber removed), followed by up to three more times between 1990 and 2008 (~37–66 m³ ha⁻¹ cumulative timber removed) (Fisher et al., 2011; Riutta et al., 2018). The OP plot was located in a stand aged ~7 years at the time of sampling, where fertilizer applications were roughly two applications of 3–4 kg bags palm⁻¹ year⁻¹ (comprising a mixture of diammonium phosphate, potassium chloride, ammonium sulfate, magnesium sulfate, and borax pentahydrate) prior to this study (J. Drewer, personal communication). All plots in OG and SL forest were previously established as part of the Global Ecosystem Monitoring (GEM) network (http://gem.tropicalforests.ox.ac.uk/) for the long-term evaluation of forest carbon cycling and productivity since 2011 (Malhi et al., 2015; Riutta et al., 2018), and have recently been characterized for tree community assemblage, plant traits, soil microbial communities and soil physicochemical properties (Elias et al., 2018; Both et al., 2019).

Sampling Design

Three 20 × 20 m subplots were randomly chosen per 1 ha plot for assessing mycelial community attributes and productivity. Ten hyphal in-growth bags were installed per subplot at a randomly chosen location. Bags were buried at 50 cm intervals along two parallel transects (5 bags per transect) spaced 1 m apart to ensure adequate sample recovery and account for spatial variability in fungal community composition and mycelial production.

Hyphal In-Growth Bags

Hyphal in-growth bags were constructed from two 5 cm × 5 cm squares of fine pore-size nylon mesh (Plastok, UK) sealed with a soldering iron. A mesh size of 41 μm was chosen to allow access to fungal hyphae, but prevent ingress of plant roots (Fisher et al., 2013; Wallander et al., 2013). Bags were filled with 25 g oven-dried, heat-sterilized quartz sand (150°C for 24 h). In-growth bags were installed between April-May 2016 by burying vertically at the soil surface covering a depth of ~0–5 cm. This was chosen to include the organic layer and interface between organic and mineral horizons to maximize fungal in-growth (Lindahl et al., 2007; Wallander et al., 2013). Canopy litterfall tends to be relatively stable throughout the year in these forests, reflecting little inter-annual climatic variation (Kumagai and Porporato, 2012). Therefore, in-growth bags were harvested after a period of 6 months to allow fungal colonization, rather than the duration of a growing season which is commonly used in temperate or boreal systems (Wallander et al., 2013). Bags were frozen in a field laboratory upon collection, and transported on ice to Lancaster University, UK where they were stored at −20°C prior to analysis. Bags were examined and discarded where there was clear damage or root ingress. Due to low recovery of undamaged bags, three in-growth bags were used per subplot to account for local variation in mycelial colonization (with the exception of five composite samples: one bulked from two in-growth bags, and four from one bag due to high damage or loss rates in some subplots). Bags were opened in the laboratory and sand was carefully bulked per subplot and hand-mixed, providing 27 composite samples (three per 1 ha plot) for quantification of mycelial abundance and molecular analysis of fungal communities. Only one 1 ha plot was established in the OP site, due to logistical challenges and as spatial variation in this land-use type was expected to be low (as suggested by low fungal beta diversity and variation in EcM relative abundances in OP relative to OG and SL in Borneo; Kerfahi et al., 2014). As such the three composite OP samples were treated as independent replicates, each representing an individual plot-level sample for all analyses.

Hyphal Extraction and Estimation of Mycelial Productivity

Hyphae were extracted from sand using an adapted floatation method (Bakker et al., 2015). Five-gram subsamples were transferred into 50 ml sample tubes, then filled up to 30 ml with 4 M KCl solution. This extractant was used to ensure floatation of all hyphae and avoid clumping and sinking of fungal filaments, which may be hydrophilic or hydrophobic (Ekblad et al., 2013), preventing adequate separation from sand particles (as found in prior laboratory tests using distilled water). Tubes were vortexed at full speed for 1 min, left to stand for 30 s to allow hyphal material to reach the surface of the solution, and 25 ml was decanted into a clean sample tube, avoiding the transfer of sand particles. This process was repeated once more on the same sand subsample. The extract was transferred into a sample pot and made up to 100 ml by rinsing the sample tube twice with 25 ml distilled water. Hyphal material was prepared for measurement using the membrane filtration technique (Hanssen et al., 1974). The extract was evenly mixed using a sample mixer at 700 rev min⁻¹, and 10 ml aliquots were transferred in two steps using a 5 ml pipette during mixing into a 15 ml glass filtration tube mounted on a nitrocellulose filter membrane (Merck Millipore, USA; 1.2 μm pore-size, 25 mm diameter). The suspension was filtered with a vacuum pump, and hyphae were stained in the tube using Lactophenol cotton blue under a fume hood for 20 min. The filter membrane was then rinsed with distilled water until filtrate ran clear, and left to air dry overnight. Filter membranes were mounted onto microscope slides with immersion oil for transparency and sealed under a coverslip.

Hyphae on filter membranes were photographed at high resolution (4,080 × 3,072 pixels) using a microscope-mounted camera at 100x magnification (Olympus BX51, Olympus DP71), allowing identification of hyphal structures with diameters ≥ 1 μm. Photographs were used for hyphal measurement to maximize membrane coverage and capture of within-sample variation not possible through direct microscope observation at higher magnification. Twenty-five photographs (each covering an area of ~1.53 mm²) were taken per membrane along four crossing transects positioned at 45° to each other, as described by Boddington et al. (1999). Hyphal length was estimated from photographs using the gridline intersect method (Tennant, 1975). A regular grid of grid size 50 × 50 μm was digitally placed on top of each photograph using ImageJ (Schneider et al., 2012; Shen et al., 2016). Image contrast and brightness were altered as necessary to improve visualization of hyphae. All intersections between hyphae and gridlines were counted, and total hyphal length was calculated using Tennant's formula. The total hyphal length in each subsample was estimated using hyphal length per area in photographs and area of membrane used for filtration. 10 g of fresh sand from each composite sample was oven-dried at 105°C to constant weight to calculate moisture content for standardization of hyphal length estimations. Hyphal length was then calculated in mm g⁻¹ dry sand for statistical analysis of mycelial production.

Soil, Environmental and Vegetation Characteristics

Soil physicochemical data for corresponding subplots were obtained from an existing dataset (Elias et al., 2018). Briefly, five soil samples were collected within each subplot (3 cm diameter gouge auger) March-April 2015. Organic soil layer depth was measured before separation from underlying mineral soil. Organic layer soil samples were bulked per subplot and analyzed for organic layer pH, total C, total nitrogen (N), total phosphorus (P), inorganic P and texture (% sand, silt and clay). pH in water was measured on fresh soils using a pH meter (1:2.5 soil to deionised water) after shaking overnight at 100 rev m⁻¹ on an orbital shaker and standing for 30 min (Landon, 1984). The remaining soils were air-dried at 40°C to constant weight and passed through a 2 mm sieve. Subsamples for total C and N analysis were dried at 65°C for 48 h and milled to a fine powder with a pestle and mortar. Total soil C and N contents were determined by dry combustion at 900°C using an Elementar Vario Max CN analyser (Elementar Analysensysteme, Hanau, Germany). Samples were digested using sulphuric acid-hydrogen peroxide (Allen, 1989) for soil total P. Inorganic P was extracted using a Bray No. 1 extractant (Bray and Kurtz, 1945). P contents of extracts and digests were determined using the molybdenum-blue method (Anderson and Ingram, 1993), read at 880 nm on a spectrophotometer (HITACHI-UV-VIS, Japan). Soil texture was measured by the pipette method (Miller and Miller, 1987). Soil bulk density was determined from one additional sample taken per subplot using a 7.5 cm diameter volumetric ring, dried at 105°C for 24 hours after removal of roots and stones (Emmett et al., 2008). Soil moisture content (top 12 cm), temperature (10 cm depth) and air temperature (20 cm above soil surface) values for corresponding subplots were accessed from datasets of continuous sampling of GEM plots, as described in Marthews et al. (2014). Slope measurements were taken upon in-growth bag harvesting using a clinometer at each subplot corner and center, and values were averaged at the subplot level. Altitude was recording using a GPS in each subplot center. Forest structural characteristics [stem density, basal area, and mean and maximum diameter of stems at breast height (DBH) ≥ 10 cm, pioneer tree proportion of basal area, and Leaf Area Index (LAI)] and tree productivity metrics [canopy net primary production (NPP), woody NPP, root NPP, total NPP] for corresponding subplots were also accessed from GEM plot datasets, using temporally-averaged data collected between 2011 and 2016 (Marthews et al., 2014; Riutta et al., 2018). Tree taxonomic community datasets were constructed at the 1 ha plot level from a previous survey undertaken July-December 2015, where all individual trees DBH ≥ 10 cm within three 20 × 20 m subplots were taxonomically identified (Both et al., 2019). Dipterocarp basal area was derived from tree community data at the plot level.

Molecular Analysis of Mycelial Fungal Communities and Data Pre-Processing

DNA was extracted from 0.2 g sand from in-growth bags using the PowerSoil^® DNA Isolation Kit and protocol (MoBio Laboratories). Amplicon libraries were constructed according to a dual indexing strategy with each primer consisting of the appropriate Illumina adapter, 8-nt index sequence, a 10-nt pad sequence, a 2-nt linker and the amplicon specific primer (Kozich et al., 2013). Fungi were targeted by amplifying the ITS2 region using primers GTGARTCATCGAATCTTTG and TCCTCCGCTTATTGATATGC (Ihrmark et al., 2012). Although the capability of detecting AM fungi using ITS primers is debated (Hart et al., 2015), recent studies have shown that patterns in diversity and community composition can be adequately identified within sample types such as soil (Berruti et al., 2017; Lekberg et al., 2018). Amplicons were generated using a high fidelity DNA polymerase (Q5 Taq, New England Biolabs). After an initial denaturation at 95°C for 2 min, PCR conditions were as follows: Denaturation at 95°C for 15 s; annealing at 52°C; annealing times were 30 s with extension at 72°C for 30 s; cycle numbers were 25; a final extension of 10 min at 72°C was included. Amplicon sizes were determined using an Agilent 2200 TapeStation system, samples were normalized using SequalPrep Normalization Plate Kit (Thermo Fisher Scientific) and pooled. The pooled library was quantified using a Qubit dsDNA HS kit (Thermo Fisher Scientific) prior to sequencing with an Illumina MiSeq using V3 600 cycle reagents at a concentration of 8 pM with a 5% PhiX Illumina control library. The sequencing run produced in excess of 18 million reads passing filter. Sequences were processed in R (R Core Team, 2013) using the DADA2 package to quality filter, merge, de-noise and assign taxonomies (Callahan et al., 2016). Sequence reads were trimmed to 225 and 160 bases, forward and reverse, respectively. Filtering settings were maximum number of Ns (maxN) = 0, maximum number of expected errors (maxEE) = 1. Sequences were dereplicated and the DADA2 core sequence variant inference algorithm applied. mergePairs and removeBimeraDenovo functions were used at default settings to merge forward and reverse reads and remove chimeric sequences. The amplicon sequence variants (ASVs) were subject to taxonomic assignment using assignTaxonomy and the training database UNITE version 7.2 (UNITE Community, 2017).

Fungal functional guild classifications were assigned to ASVs using the FUNGuild annotation tool (Nguyen et al., 2016). Only ASVs with unambiguous (non-multiple) classifications of “probable” or “highly-probable” confidence rankings were considered for analysis. These were used for calculating relative abundances of fungal guilds and sub-setting saprotrophic, mycorrhizal and pathogenic fungal datasets for assessment of diversity and community dissimilarity (EcM and AM fungal diversity and community dissimilarities could not be analyzed separately due to absence of read counts in a number of samples). Sequencing data were pre-processed (steps described below) and alpha diversity indices (ASV richness, Shannon index) and fungal guild relative abundances calculated in R version 3.5.1 (R Core Team, 2013) using the phyloseq package (Mcmurdie and Holmes, 2013). Only ASVs assigned to the kingdom of Fungi were retained for downstream analysis (99.34% of total reads), and all singleton ASVs were removed. Sub-setting by fungal guilds was conducted on the full unrarefied dataset to maximize the number of ASV reads available for analysis of functional groups. Sample sequencing depth was normalized for each group by rarefying to the minimum read counts of 4881, 151, 75, and 67 per sample for overall, saprotrophic, mycorrhizal and pathogenic fungal groups, respectively. Analyses of diversity metrics and community dissimilarities were repeated using unrarefied and whole-sample only rarefied datasets (before sub-setting fungal groups), which showed broadly consistent results.

Statistical Analyses

All statistical analyses were conducted in R version 3.5.1 (R Core Team, 2013) and significance was considered at the p ≤ 0.05 level. For univariate analyses, linear mixed effects regression models (LMMs) were constructed in the lme4 R package (Bates et al., 2015). This included testing (H1) the effect of land-use type on relative abundances of mycelial fungal guilds and alpha diversity, mycelial production and soil, environmental and vegetation characteristics, and (H2 & H3) relationships between these variables across OG and SL plots. Post-hoc pairwise comparisons were conducted with the emmeans R package (Lenth et al., 2019) with Bonferroni correction to identify statistically different variable means between OG, SL and OP land-use types. To control for potential within-plot pseudoreplication, plot ID was included as a random intercept term. Significance was evaluated using the Satterthwaite degrees of freedom approximation (Luke, 2017). Normality of model residuals were evaluated using Shapiro-Wilk tests, and variables were log- and exponentially-transformed where necessary to improve model fit. Kruskal-Wallis tests were conducted when residual normality could not be satisfactorily achieved using data averaged at the plot level.

Mycelial fungal community compositions across all land-use types were visualized with PCoA using Bray-Curtis dissimilarities via the phyloseq, vegan (Oksanen et al., 2019) and ggplot2 (Wickham, 2016) packages. To test differences in fungal community compositions between OG, SL and OP (H1), Bray-Curtis community dissimilarities were calculated from data averaged at the plot level (n = 11) using the merge_samples function in phyloseq. Differences between land-use types were tested with PERMANOVA using the adonis vegan function, and statistically different groups were identified using the pairwise.adonis function in the pairwiseAdonis R package (Martinez Arbizu, 2019) controlling for the False Discovery Rate (FDR). All permutational tests were run with 9,999 permutations with the exception of pairwise multiple comparisons, where full enumeration was used. Indicator analyses were conducted to identify specific fungal taxa associated with different land-use types using the labsdv R package (Dufrêne and Legendre, 1997; Roberts, 2016). Differences in tree community composition between OG and SL forest (H1) were tested with PERMANOVA using Bray-Curtis community dissimilarities calculated from plot-level community data. Hellinger-transformation was applied to all community data prior to multivariate analyses (Legendre and Borcard, 2018) to control for the effect of rare taxa. Homogeneity of multivariate dispersion across land-use types (an assumption of PERMANOVA) for fungal and tree community dissimilarities was evaluated using the betadisper vegan function for overall and pairwise tests (FDR-corrected).

Relationships between mycelial fungal community compositions and soil, environmental and vegetation characteristics (H2) across forest (OG and SL) plots were evaluated using distance-based redundancy analysis (db-RDA). Tree community composition was included in analysis as represented by the first three PCoA axes, explaining most of the variation (60.39% cumulative eigenvalues). Prior to analysis, highly correlated variables within 1) soil and environmental and 2) vegetation groups were identified with correlograms using the corrplot R package (Wei and Simko, 2017). Variables correlated with Pearson's r ≥ ~ 0.7 were removed. Soil and environmental variables were treated as one group as soil characteristics were considered to be of primary ecological importance, while driven by environmental factors. The best predictors of fungal Bray-Curtis dissimilarities were identified through forward-selection using the criteria of adjusted R² and significance level of p < 0.05 (Blanchet et al., 2008) with the ordiR2step vegan function. Here, Bray-Curtis dissimilarities were calculated at the subplot level and permutational tests were restricted by plot ID to control for the nested sampling design using the permute R package (Simpson et al., 2019). Variation in fungal community dissimilarities was partitioned by soil and environmental and vegetation components using the varpart function, and significance of components and individual predictors was tested with partial db-RDA in vegan. Relationships between mycelial fungal and vegetation communities (H2) were tested using plot-level Bray-Curtis dissimilarities with Mantel tests (Spearman's rank correlation) in vegan. Relationships between mycelial fungal community compositions and mycelial production (H3) were tested with PERMANOVA using Bray-Curtis dissimilarities and the same permutational scheme described above.

Results

Overall, 3,565 fungal ASVs from 10 phyla (Figure 1) and 374 genera were detected across all samples of OG forest, SL forest and OP. Of the ASVs that were assigned fungal functional guilds, the majority of reads comprised saprotrophic (39.9%), followed by mycorrhizal [27.0%: 83.1% of which EcM; 16.7% AM; 0.2% ericoid mycorrhizal (ErM)], pathogenic (21.0%: 51.4% of which animal pathogenic; 48.6% plant pathogenic), parasitic (7.1%), endophytic (2.2%), lichenised (2.0%) and epiphytic (0.7%) fungi.

H1 - Impact of Land-use Type on Mycelial Fungal Community Attributes, Mycelial Productivity and Soil, Environmental, and Vegetation Characteristics

Land-use type significantly affected community dissimilarities for overall, saprotrophic, mycorrhizal and pathogenic fungal groups (Figures 2A–D; see Table 1 for summary of statistics). Pairwise comparisons identified significant differences between OP and forest land-uses for overall and pathogenic fungal community dissimilarities. No significant differences were found between land-use types for saprotrophic or mycorrhizal community dissimilarity in pairwise comparisons (p > 0.05). However, further tests using only OG and SL plots showed no significant effect of land-use type on saprotrophic or mycorrhizal community dissimilarities (PERMANOVA: p = 0.223; p = 0.115, respectively), indicating OP to be the driver of the significant effect found across the full disturbance gradient for these fungal groups. Community dissimilarity dispersions were homogenous between all land-use types for all fungal groups (betadisper: p > 0.05). No significant differences were found in community dissimilarities or alpha diversity metrics (ASV richness and Shannon index) between OG and SL for any fungal group (p > 0.05), despite significant differences in soil, environmental and vegetation properties between forest land-uses (means by land-use type summarized in Table 2; test statistics for significant differences are reported in Supplementary Information 2). However, relative abundance of AM fungi was significantly higher in SL relative to OG (Figure 3) and ErM fungal relative abundance was higher in SL relative to all other land-use types (Figure 3) (see Table 3 for means of all mycelial fungal community attributes by land-use type; test statistics for significant differences are reported in Table 4). Overall pathogenic fungal relative abundance was marginally significantly higher in SL compared to OG when tested without OP samples (F = 4.25, R² = 0.16, p = 0.051). The vast majority of mycorrhizal reads were attributed to EcM fungi in forest plots (89.6%: 98.31% in OG; 80.90% in SL), but were mostly AM fungi in OP (69.4% AM; 30.6% EcM). Total mycorrhizal and EcM fungal relative abundances were significantly lower in OP compared to forest plots (Figure 3). Mean parasitic fungal relative abundance was an order of magnitude higher in OP relative to OG and SL, although no significant differences between land-use types were found due to the considerable variation in OP samples (post-hoc tests: p > 0.05). For alpha diversity metrics, only mycorrhizal fungal ASV richness was significantly affected by land-use type (F = 5.80, R² = 0.33, p = 0.025; Figure 4A), with lower values in OP compared to forest plots (OG—OP: p = 0.028; SL—OP: p = 0.018). Alpha diversity indices (ASV richness, Shannon index) did not differ between OG and SL for any fungal group (post-hoc tests: p > 0.05). Indicator analysis identified certain fungal taxa of the families Aspergillaceae, Sporocadaceae, and Nectriaceae to be most indicative of OG forest, while those of Ophiocordycipitaceae were characteristic of SL forest and Nectriaceae, Aspergillaceae and Ophiocordycipitaceae characteristic of OP (Supplementary Information 3).

Hyphal length values ranged from 0.24 to 8.50 m g⁻¹ of dry sand across the disturbance gradient, and were significantly affected by land-use type (F = 6.50, R² = 0.33, p = 0.006; Table 3; Figure 4B). Hyphal length values were significantly lower in OP relative to forest land-use types (OP–OG: p = 0.007; OP—SL: p = 0.008) but did not significantly differ between OG and SL (p = 0.999).

Relative to OG, significantly higher values were associated with SL for soil bulk density, altitude and proportion of pioneer species, and lower values for LAI and canopy, woody and total NPP (Table 2; Supplementary Information 2). Root NPP, soil clay content and organic layer depth in OP were significantly lower relative to SL, while soil silt content was significantly higher. OP soil bulk density was significantly higher relative to OG. Total and inorganic soil P were an order of magnitude higher in OP relative to both OG and SL. While stem density was significantly lower in OP relative to both forest types, basal area and mean stem diameter were significantly higher in OP. Maximum stem diameter was also higher in OP relative to SL, however all belonging to the same species of oil palm. Plot-level tree community dissimilarity significantly differed between OG and SL (PERMANOVA: F = 1.43, R² = 0.19, p = 0.028; betadisper: p > 0.05).

H2 and H3 - Relationships Between Forest Mycelial Fungal Community Attributes, Soil, Environmental and Vegetation Characteristics and Mycelial Production

Overall, saprotrophic, mycorrhizal and pathogenic fungal community dissimilarities were significantly related to soil and vegetation characteristics, although overall variance explained was low for each group (see Table 5 for summary of statistics). The best predictors selected for overall fungal community dissimilarity were soil pH and LAI. Of these, soil pH was negatively related to relative abundance of animal pathogenic fungi (F = 8.63, R² = 0.27, p = 0.008), while LAI was negatively related to relative abundance of AM fungi (F = 8.04, R² = 0.28, p = 0.015). The only predictor selected for saprotrophic community dissimilarity was soil inorganic P, while the best predictors for mycorrhizal fungal community dissimilarity were pH and dipterocarp basal area. Soil pH, tree community PCoA axis 2, tree basal area and soil inorganic P were all identified as predictors of pathogenic fungal community dissimilarity. Plot- level tree community dissimilarity was significantly related to plot-level community dissimilarity of mycorrhizal fungi (r = 0.55, p = 0.015; Figure 5), but no other fungal group (p > 0.05). Dipterocarp and mycorrhizal fungal community dissimilarities were not significantly related at the plot-level (p > 0.05). Additional tests showed no significant relationship between either EcM or AM fungal relative abundance and soil C content (p > 0.05).

Hyphal length was not significantly associated to any mycelial fungal community attributes (H3), although a negative relationship with EcM fungal relative abundance was marginally significant (F = 3.57, R² = 0.13, p = 0.072). There was no significant link between hyphal length and any of the vegetation properties measured (H2), including tree productivity metrics, but was significantly positively related to soil inorganic P (F = 5.25, R² = 0.19, p = 0.032; Figure 6). However, the association between hyphal length and soil inorganic P was primarily driven by one sampling point in SL forest with high values in both parameters (Figure 6), and was found to be non-significant when this data point was excluded from analysis (F = 0.83, R² = 0.04, p = 0.371).

Discussion

We did not find significant shifts in community structure of any mycelial fungal group between OG and SL forest (Figures 2A–D), despite clear differences in soil and vegetation properties (higher soil bulk density, proportion of pioneer trees, woody and total NPP, and lower LAI and canopy NPP in logged forest) and plot-level tree community composition, broadly in line with previous findings from these study sites (Riutta et al., 2018; Both et al., 2019). This finding contrasts with our hypothesis (H1) and recent surveys which found the compositions of bulk soil fungal communities in tropical Southeast Asian dipterocarp forest to be highly sensitive to logging (Kerfahi et al., 2014; Mcguire et al., 2015). Our results suggest that actively foraging fungi may be more resilient to disturbance than the wider fungal community in bulk soil. However, direct comparability with other studies may be affected by methodological differences in sampling related to community capture rates (Wallander et al., 2001; Ekblad et al., 2013; i.e. micro-environmental conditions inside in-growth bags potentially selecting for certain fungal species). We also highlight that the in-growth bag method recovers only active hyphae, therefore avoiding confounding effects of relic DNA on soil microbial community patterns that can be a concern in studies across environmental gradients (Lennon et al., 2018). The total number of ASVs detected was lower than found in other recent ITS amplicon sequencing studies of bulk soil fungal communities in the region, which may exceed 25,000 (Elias, unpublished data), which we similarly attribute to these reasons. We identified differences in other mycelial fungal attributes between OG and SL forest, including a significant increase in the relative abundances of AM and ErM fungi in SL forest. The vast majority of forest mycorrhizal reads comprised EcM fungi (Figure 3). The canopy of lowland rainforests of Borneo is typically dominated by species belonging to the Dipterocarpaceae, which are an EcM-associating family (Taylor and Alexander, 2005; Brearley, 2012). Although our finding may relate to presence of host tree species, high EcM relative to AM fungal relative abundance may also result from amplification bias associated with ITS primers (Hart et al., 2015). Although mean dipterocarp basal area (Table 2) and relative abundance of EcM fungi (Table 3; Figure 3) were lower in SL forest, we did not detect significant differences in either metric due to the high variability in both. We nevertheless consider the relative increase in abundances of AM and ErM fungi in SL forest (the only other mycorrhizal types identified; Figure 3) to indicate increasing evenness in mycorrhizal types and diminished dominance of EcM fungi. This pattern reflects differences in vegetation indicated by the increased proportion of non-dipterocarp pioneer species in SL forest, which typically form arbuscular mycorrhizas (Mcguire et al., 2008). This shift may have implications for nutrient and C cycling, as recent correlative studies have suggested EcM- rather than AM-dominated communities may promote soil C storage on a global scale (Averill et al., 2014) through their ability to access N and P (Liu et al., 2018) from organic sources, potentially competing with free-living soil decomposers requiring N for the breakdown of organic matter (Averill et al., 2014). An increase in relative abundances of AM fungi may, therefore, have implications for the release of C from the soil as N becomes more available for soil saprotrophs. However, our analysis revealed no relationship between EcM fungal relative abundance and soil C concentration. It is possible that patterns reported in other studies (e.g. Averill et al., 2014) reflect much longer timescales in dominance of a particular mycorrhizal type than seen in our study. Alterations in soil C-cycling processes in SL dipterocarp forest have been observed, although as litter decomposition has been found to be slowed in SL forest due to subtle changes in forest microclimatic conditions (Both et al., 2017), the overall consequences for biogeochemical cycling and ecosystem functions are unclear.

Overall fungal community structure within in-growth bags was significantly different between OP plantation and OG and SL forests, and the impact of OP as a driver of differences over land-use types for saprotrophic, mycorrhizal and pathogenic mycelial fungal communities were also observed (Figures 2A–D), which corroborates the effects of forest conversion on bulk soil fungal communities (Mcguire et al., 2015). Mycorrhizal richness was also strongly affected (Figure 4A), with much lower values in oil palm relative to forest plots. As we found no differences between OG and SL forest, overall mycelial fungal communities appear to follow a similar pattern to soil bacterial communities, which have been shown to be generally unaffected by logging but distinctly different in oil palm plantation (Lee-Cruz et al., 2013; Tripathi et al., 2016).

In contrast to expectations (H1), mycelial production did not significantly differ between OG and SL forest, but was significantly less in OP relative to both forest land-use types (Figure 4B). We also found that hyphal length did not correspond to shifts in mycelial community attributes, soil properties and vegetation community between SL and OG forest, which may reflect resilience of some soil functions to selective logging (e.g. enzyme activity; Mcguire et al., 2015). The resilience of both mycelial production and relative abundance of mycelial EcM fungi implies their extramatrical mycelium may be largely retained in SL forest, when individual hosts remain. Mycorrhizal mycelial networks are crucial for tree seedling establishment (Nara, 2006), and may facilitate interplant exchange of resources with trees sharing the same mycorrhizal partnerships (Gorzelak et al., 2015). Our findings indicate implications for the potential restoration of these degraded forests, as EcM mycelial networks vital for supporting dipterocarp recruitment and survival may remain even after individuals of this family have been removed. Although basal area was significantly highest in OP, we found significantly lower mean root NPP in OP relative to SL (Table 2). There was also a complete removal of dipterocarp species through forest conversion to agricultural plantation. Furthermore, although canopy NPP was not measured in the OP site, significantly lower soil organic layer depth (Table 2) may indicate reduced litter inputs – especially as dead palm fronds are typically removed and collected in localized areas as standard in OP. This suggests a combination of a lack of EcM fungi, due to replacement by AM-associating oil palm species (Phosri et al., 2010), diminished belowground allocation of C by trees and reduced organic matter inputs may be down-regulating mycelial production rates in this land-use type. This corresponds to a ten-fold increase in the availability of limiting nutrients (inorganic P; Table 2) in OP due to the use of fertilizer. As the OP plantation was established after the selective-logging of these forests, lower hyphal length values, mycorrhizal abundance and richness appears to result from magnitude of EcM host removal rather than time since disturbance. The significant reduction in the soil (mycorrhizal) mycelium through the extreme disturbance of forest conversion represents a substantial barrier to the restoration of these systems.

Across OG and SL forest, mycelial fungal community attributes significantly varied with vegetation and soil properties (H2), but did not correspond with mycelial production (H3). Soil microbial community assemblage, particularly fungi, have been shown to correlate with aboveground plant taxonomic and phylogenetic structure in addition to local soil properties (Barberán et al., 2015). Soil pH was a significant predictor of community composition for all fungal groups, especially mycorrhizal and pathogenic fungi, a pattern observed across large scales in temperate ecosystems (Dupont et al., 2016). Mycorrhizal fungal community structure was surprisingly not affected by soil texture, previously found to be a strong predictor of EcM communities in Bornean lowland dipterocarp rainforest (Essene et al., 2017). It is possible that the use of sand inside in-growth bags masked the influence of bulk soil texture. Tree and mycorrhizal fungal community structure were significantly correlated at the plot-level (Figure 5), likely driven by the significant relationship between mycorrhizal fungal community dissimilarity and total dipterocarp basal area (Table 5). Previous studies in Bornean lowland dipterocarp rainforest have shown EcM fungal communities to be related to dipterocarp species assemblage in addition to strong effects of soil properties, suggesting dipterocarp species distributions across different soil types to be mediated by assemblage of their mycorrhizal partners (Essene et al., 2017). However, we found no significant correlation between mycorrhizal fungal and dipterocarp community dissimilarities, indicating overall abundance rather than composition of dipterocarps is more important in structuring actively foraging mycorrhizal fungal communities in the present study.

Mycelial productivity was significantly positively related to soil inorganic P (Figure 6), although mycelial biomass has previously been shown to be greater in more P-deficient forests (Potila et al., 2009). However, in some cases the proliferation of mycorrhizal hyphae into nutrient ‘hotspots' has been observed both in laboratory experiments and natural systems (Ekblad et al., 2013), whilst the role of inorganic P in controlling rates of mycelial production has been demonstrated in forests using hyphal in-growth bags augmented with mineral P, and field manipulations using fertilizer applications (Hagerberg et al., 2003; Potila et al., 2009; Ekblad et al., 2013; Camenzind et al., 2016). In the current study, the tenuous nature of this relationship (being mainly driven by one sampling point in SL; Figure 6) limits interpretation of underlying mechanisms, and should be approached with caution. Further study is required to confirm this link, which may be achieved with alterations in sampling design discussed below.

Some limitations must be recognized with respect to the interpretation of the results in this study. Primarily, soil sampling was conducted 12-14 months before installation and 18-20 months before collection of hyphal in-growth bags. Any alterations in soil properties arising from climatic factors, vegetation dynamics (i.e. regeneration of SL forest) or OP management during this period will not have been accounted for. Unlike rainforests in other tropical regions, the forests of Borneo do not experience a defined wet or dry season (Richards, 1996; Dykes, 1997) reflected in relatively stable canopy litter production throughout the year (Kumagai and Porporato, 2012). Differences in time of year sampled would not be expected to have a pronounced effect (i.e. due to coupled precipitation-driven seasonal dynamics in litter inputs, soil conditions and soil microbial communities observed other tropical systems; Buscardo et al., 2018), although the potential influence of supra-annual climatic and phenological variation (Itioka and Yamauti, 2004) and irregular inter-annual dry spells (Walsh and Newbery, 1999; Kumagai et al., 2005) should be noted. Furthermore, soil data represent bulk soil conditions at the subplot scale (10 × 10 m), rather than those within bags to which fungi were directly exposed. The combined temporal and spatial decoupling of soil and fungal community-productivity data may underlie the general lack of associations found between soil characteristics and mycelial productivity (also problematic for confidence in relationships observed, i.e. inorganic P; Figure 6) or fungal community attributes (e.g. little overall variation in community dissimilarities explained by soil factors; Table 5). Future studies should aim to concurrently sample for soil physicochemical properties and soil fungal communities and mycelial productivity (during installation, incubation and harvesting of hyphal-ingrowth bags). Corresponding molecular analysis of bulk soil fungal communities for comparison with in-growth bags, in tandem with a temporally coupled approach, would help to confirm the reciprocal mechanistic relationships between soil properties and fungal community-productivity patterns.

The present study analyzed data from 27 composite samples from 9 plots distributed across three land-use types at the landscape scale (Supplementary Information 1). Although composite samples comprised material from 72 individual in-growth bags to control for local variation, an unexpectedly high level of variability in fungal community attributes and mycelial productivity was observed across land-uses (particularly within the two forest types; Figure 4B). This variability (Table 3) makes it difficult to ascertain effects (or lack thereof) of management, and future surveys should incorporate greater plot-level replication in light of these findings. Similarly, site effects cannot be excluded from differences detected in fungal community attributes and mycelial productivity between OP and forest samples, as replication was restricted to only one 1 ha plot in this land-use type. While consistency with previous studies in peninsular Malaysia and Borneo (e.g. major reductions in EcM relative abundance in OP relative to OG and SL dipterocarp rainforest; Kerfahi et al., 2014; Mcguire et al., 2015) suggests findings to be representative of the broad impacts of forest conversion to OP on actively-foraging fungal communities, caution should be exercised in their interpretation.

The in-growth bag method has been found to adequately represent mycorrhizal mycelial productivity, as demonstrated by δ¹³C (Wallander et al., 2001) and molecular (Kjøller, 2006) analyses of bags in EcM-dominated systems. However, our results indicate that this may not be the case in tropical forest, with the largest proportion of overall reads belonging to ASVs identified as saprotrophs. As such, hyphae measured cannot be directly attributed to one fungal group or mycorrhizal type in this study. Incorporation of biomarker analyses to partition biomass by fungal types, or the inclusion of additional measures to reduce saprotrophic mycelial ingress are needed improve understanding of mycorrhizal productivity alone (Wallander et al., 2013). We found a high level of variability in mycelial fungal community dissimilarities across forest plots. Our analysis focused on the topsoil, as this is where the vast majority of nutrient turnover occurs. Here, fine-scale environmental heterogeneity can drive greater spatio-temporal fungal community variation compared to underlying mineral soil where communities may differ (Bahram et al., 2015). Future studies considering effects of forest degradation on mycelial fungal communities and productivity at different soil depths may identify different patterns along the soil profile.

In conclusion, selective logging did not significantly shift mycelial fungal community structure or productivity from OG forest, suggesting mycelial fungal communities and function may be relatively resilient to forest degradation compared to fungi in bulk soil. Results indicate the extramatrical EcM mycelium to be largely retained in selectively logged forest, with positive implications for potential restoration of dipterocarp forest by providing existing mycorrhizal networks for tree seedling establishment. However, SL forest was associated with higher relative abundances of AM and ErM fungi relative to OG, which may have consequences for soil C cycling and storage in lowland rainforests in Borneo over longer timescales than those reported here. Mycelial production was not related to any fungal community attributes, and broadly independent of vegetation, soil and environmental characteristics apart from a weak positive relationship with soil inorganic P concentration. In contrast, land conversion to oil palm plantation was shown to have a significant effect on overall, saprotrophic, mycorrhizal, and pathogenic mycelial fungal community structures, reducing mycorrhizal relative abundance and diversity. Mycelial production was significantly lower in oil palm relative to forest, corresponding to soil inorganic P concentrations of an order of magnitude greater under this land-use type. The impact of forest conversion on fungal community attributes and mycelial productivity may have wide repercussions for nutrient dynamics, plant-soil interactions and restoration potential over large scales as oil palm plantation has rapidly become a major land-use across Southeast Asia.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher. The molecular datasets generated for this study can be found in the European Nucleotide Archive with study accession no. PRJEB357.

Author Contributions

SR, DE, DJ, SB, NPM, and NO contributed to the initial conceptual design of the study. SR and DE established study plots and collected soil samples and data in the field. TG conducted DNA sequencing analysis and bioinformatics to generate molecular dataset. RG and TG provided conceptual approach to statistical analysis of microbial communities. SR, DJ, DE, and NPM contributed to hyphal extraction method. DE, NPM, and NM conducted soil analysis and provided soil characteristics dataset. SB provided tree community and functional trait dataset. TR provide tree productivity dataset. SR conducted statistical analysis and wrote the first draft of the manuscript. DE, SB, TR, NPM, RG, NO, DJ, and TG contributed writing to the manuscript. All authors provided critical revision of intellectual content, and approved the manuscript for submission.

Funding

This study was funded through the UK Natural Environment Research Council (NERC) Human Modified Tropical Forests programme (grant number NE/K016377/1) and Sime Darby Foundation funding to the SAFE project. This publication is a contribution from the NERC Biodiversity and Land-use Impacts on Tropical Ecosystem Function (BALI) consortium. DJ receives partial support from the N8 Agrifood programme.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We express our gratitude to the Sabah Biodiversity Council (SaBC), Yayasan Sabah (Sabah Foundation), the Danum Valley and Maliau Basin Management Committees, Glen Reynolds and the South East Asia Rainforest Research Partnership (SEARRP), and Rob Ewers for support and granting permission to access field sites [SaBC Access License no. JKM/MBS.1000-2/2 JLD.3 (173)]. We are extremely grateful for support from Laura Kruitbos, Unding Jami, and all research assistants in Sabah, without whom this study would not have been possible.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2020.00064/full#supplementary-material

References