Introduction

Space determines the nature of and scale over which individuals meet and interact. The characteristics of the discrete spatial habitat which an organism occupies affects individual competitive success with a bottom-up effect on population-wide colonisation, speciation and extinction [1]. There is a dynamic link between spatial ecology and competitive success where transitive (where species A outcompetes B, which outcompetes C) communities with a strict competitive hierarchy become intransitive (A > B; B > C; C > A, like the game of rock–paper–scissors) when competing in a spatially more complex system [2], allowing individuals outcompeted under some scenarios to coexist with their competitors [3]. Dimensionality of habitat landscapes influences individual behaviour [4], and stochasticity of species interactions results in changes to the pool of community-produced metabolites [5], altering individual combative ability, community succession and structure between 2- and 3-dimensional landscapes [6]. Despite these findings, the mechanisms that influence stability and succession in the context of how communities occupy and exploit space are rarely adequately quantified as most ecological study systems are too complex and largely intractable. A model system is needed that allows such quantification, and understanding of how altered dynamics, coexistence and community-scale biodiversity in the context of space underpins changes in community function.

The functional diversity–area relationship, i.e., the correlation between increased habitat size and greater functional diversity, is one theory explaining how space mediates function in biodiverse communities [7]. However, the model does not account for effects of distributions and patch dynamics of species within habitats of varied area, yet these factors cause competitive communities to shift between hierarchical transitive and non-hierarchical intransitive relationship states [2, 8]. Non-hierarchical intransitivity is an established mechanism of coexistence [9, 10] and is thought to be associated with an intrinsically related functional-diversity mechanism [11, 12]. Wood decay fungi offer an ideal ecologically relevant model for the study of these processes as individual mycelia occupy columns of decay forming complex 3-dimensional communities in wood and, through their decomposition activities release carbon and nutrients [13]. They typically form a hierarchical community structure (tertiary/late stage colonists outcompete secondary colonists which in turn outcompete primary/earliest colonists) [14] and competitive interactions between species can be easily observed [15]. In addition to compounds that primarily function in the exploitation and decomposition of lignocellulose [16, 17], these fungi produce a plethora of potentially antagonistic compounds which function in changing territory, and differ in quantity and identity during interspecific competition [6, 18,19,20].

Here, we used a tractable system of wood decay fungi to quantify the impact of space on the mechanisms of coexistence and community composition, in the context of its occupation and exploitation. We compared the combative abilities of fungi in linear 2-dimensional systems and species richer 3-dimensional systems. The study system allowed more detailed analysis of community dynamics between 3-dimensional systems where fungi were dispersed in evenly spaced patches and 3-dimensional systems where the weakest member of the community occupied the same volume but as a larger adjacent patch size. Previously, our data revealed an emergent property where intransitivity promoted biodiversity in more spatially diverse 3-dimensional systems where territory was less fragmented [2]. Here, we assess the underlying mechanisms causing community stability and coexistence dynamics to change. We do this by measuring the metabolic response of an individual to changing coexistence dynamics across spatial scales. Our large-scale untargeted metabolomics and other chemical methods analysed a comprehensive network of intracellular, extracellular and gas-phase metabolic products produced during community interactions. We hypothesised that stochasticity of species would influence functional biochemical processes, and that changes to metabolites involved in pathways for resource utilisation and antagonism would alter coexistence dynamics and community composition. Spatially heterogeneous systems containing non-hierarchical communities promote biodiversity [2, 21], and here we deepen our understanding of this concept with the novel finding that space occupied alters metabolic function and coexistence, therefore, moderating the diversity–function relationship.

Methods

Experimental design and sampling

We constructed pair-wise interactions of 2 cm3 Fagus sylvatica (beech) blocks that had been precolonised for 12 weeks by placing on agar (5 gL−1 malt extract, 15 gL−1 agar; Lab M, UK) cultures of field isolates (fruit body/wood) of Vuilleminia comedens (strain VcWvJH1; a primary coloniser), Trametes versicolor (strain TvCCJH1; a secondary coloniser) and Hypholoma fasciculare (strain Hf DD3; a late secondary coloniser) wood decay fungi (maintained in the Cardiff University fungus culture collection), which all co‑exist in nature, at 20 °C (as in [2]). Interactions were performed in all combinations including, self-pairings (n = 10). We also constructed 3 × 3 × 3 27-block cubes containing 9 blocks of each species. In total, 27-block cubes were arranged: (1) with all three species dispersed and no two blocks of the same species in contact (n = 10); (2) so that all 9 blocks of the weakest competitor, V. comedens, occupied an adjacent volume while the other two fungi were dispersed (n = 10); (3) entire 27-block assemblages containing each fungus alone (Fig. 1) [2]. Blocks in pair-wise interactions were arranged with cut vessel ends touching, and in 27-block cubes blocks were joined such that some cut vessel ends were touching but others were not, but all vessels were parallel (Fig. 1). Fungal species and specific strains were selected based on their successional order and expected combative hierarchy in the natural environment: H. fasciculare > T. versicolor > V. comedens [22,23,24] (Supplementary Table 1). Interacting combinations of wood blocks were incubated individually at 20 °C in 70 and 500 ml polypropylene pots for pair-wise and 27-block cubes, respectively, and were laid upon a layer of perlite (20 and 85 ml) containing 2 and 12 ml of water respectively which was maintained by weekly addition of water to retain the original water content, as detailed in [2]. After 1 and 28 days volatile organic compound (VOC) production was measured (described below), and we deconstructed n = 5 systems. Each individual block was split along the grain into quarters, and three of the quarters from each block were foil wrapped, flash frozen in liquid nitrogen and stored at −80 °C. These were subsequently analysed for enzyme activity and metabolomics (see below). For the remaining quarter, two chips (~2 mm3) were taken from an inside face and inoculated onto 2% malt agar, incubated at 20 °C, then emerging mycelia identified morphologically, to determine which species occupied the wood (Fig. 2). We removed the chip excised face from the quarter by splitting with a chisel, and determined final density as dry weight (80 °C for in excess of 72 h) per fresh volume (cm3) (blocks sampled at 28 days only), and the rate of decay estimated by comparison with density of blocks scarified (n = 10) at the end of the precolonisation period.

Fig. 1: Spatial distribution of species interactions.
figure 1

Interactions were constructed in pair-wise (a–d) and 3-dimensional 27-block (e–g) arrangements. a–d pair-wise interactions in all conceivable combinations, plus T. versicolor self-pairing; e dispersed cube (fungi were dispersed throughout the system and arranged so that no two blocks containing the same species had adjacent faces); f ‘wall’ distribution cubes (all fungi occupied the same total volume of wood but the adjacent territory occupied by V. comedens was larger; the other two competitors were arranged so that no two blocks containing the same species were adjacent); g single species, T. versicolor 27-block cube. Cut vessels in pair-wise interactions and rows within 27-block layers were touching so that the wood grain ran in the same parallel direction as denoted by arrows.

Full size image
Fig. 2: Sampling and deconstruction of wood blocks.
figure 2

Individual blocks were split along the grain into quarters, and three of the quarters from each block were foil wrapped, flash frozen in liquid nitrogen and stored at −80 °C. For the remaining quarter, two chips (~2 mm3) were taken from an inside face and inoculated onto 2% malt agar, then emerging mycelia identified morphologically.

Full size image

Overview of metabolite analysis

To quantify changes to metabolic function (antagonistic chemicals and compounds for habitat exploitation/wood decomposition) associated with different spatial dynamics in response to changing community composition, we extracted and measured: the profile of VOCs from the headspace of interactions, and the activity of 12 targeted enzymes, chosen because they are directly involved in interspecific competition [6, 19]. We also conducted ultra-high performance liquid chromatography-mass spectrometry (UHPLC-MS) metabolomics analysis. All analyses were conducted after 1 and 28 days interactions in blocks originally colonised by T. versicolor (n = 5), and in 27-block interactions ‘pseudo’-replicates representing different spatial locations within 3-dimensional cube arrays were analysed (all n = 5). We chose to target T. versicolor since the fungus has been well characterised in the past [16, 19, 25], and because it neither dominated systems (as did H. fasciculare) nor was it driven to near extinction (as was V. comedens). There were no significant differences (ANOVA: p > 0.05) in enzyme activity or small metabolite levels between pseudoreplicates of T. versicolor in 3-dimensional cubes, so we pooled activities of all pseudoreplicates for each system (therefore, n = 15).

VOC extraction and data preprocessing

We collected VOCs from the headspace of interactions after 1 and 28 days (n = 3) by inserting pots individually and lidless into a multi-purpose roasting bag (46 × 56 cm; Lakeland, UK), which was sealed for 30 min to allow VOCs to equilibrate in the headspace. Then, 500 ml headspace gas was collected onto thermodesorption (TD) tubes (Tenax TA & Sulficarb, Markes International Ltd) using an EasyVOC manual pump (Markes International Ltd, UK).

VOCs were desorbed using a TD100 TD system (Markes International Ltd, UK) with the following settings: tube desorption 10 min at 280 °C, at a trap flow of 40 ml min−1; trap desorption and transfer 40 °C s−1 to 300 °C, with a split flow of 20 ml min−1 into gas chromatograph (GC; 7890A; Agilent Technologies Inc., USA). VOCs were separated over 60 m, 0.32 mm I.D., 0.5 μm Rx5ms (Restek, UK) with 2 ml min−1 helium as carrier gas under constant flow conditions using the following temperature programme: 35 °C for 5 min, 5 °C min−1 to 100 °C, hold 5 min. Mass spectra were recorded from m/z 30 to 350 on a time-of-flight mass spectrometer (BenchTOF-dx, Markes International Ltd, UK). C8–C20 alkane standard (0.5 μl, Supelco) was loaded onto a blank TD tube as a retention standard and quality control (QC).

GC-MS data checked with MSD ChemStation software (E.02.01.1177; Agilent Technologies, Inc.) and chromatograms were deconvoluted and integrated with AMDIS (NIST11) using a custom retention-indexed mass spectral library. MS spectra from deconvolution were searched against the NIST 2011 library (Software by Stein et al., version 2.0 g, 2011). VOCs scoring > 80% in forward and backward fit and a retention index (RI) match of +15 were included into the custom mass spectral library as putatively identified VOCs; VOCs scoring > 80% in forward and backward fit and no RI match were included as chemical class, e.g., alkane, alkanol and recurrent components that did not show either the required mass spectral fit or RI match were added as ‘unknown’. Peak list from integration with AMDIS were aligned using the pivot function in Excel in preparation for subsequent statistical analysis.

Enzyme assays

For enzyme assays, we freeze dried the T. versicolor frozen blocks for 48 h (Edwards Modulyo, UK), then ground them to sawdust using a spice and coffee grinder (Wahl James Martin, UK). In total, 0.5 g of sawdust was added to 5 ml of 50 mM sodium acetate buffer and shaken overnight at 4 °C. For pair-wise interactions, T. versicolor blocks from all interactions and one block from its respective self-pairing were used (n = 5). In 27-block cubes, each fungus occupied 3 different spatial positions in both of the mixed-species assemblages, and T. versicolor occupied 4 spatial positions within the assemblage which it fully occupied. For each spatial position for T. versicolor blocks within 27-block systems (excluding the central-cube position which T. versicolor occupied when it was the sole occupant, as this was not represented in the mixed-species assemblages), 5 replicates were used for assays (full details in Supplementary Table 2).

The activities of the following terminal hydrolases were measured using 4 methylumbelliferol (MUF)-based substrates: β-glucosidase (EC 3.2.1.21), α-glucosidase (EC 3.2.1.20), cellobiohydrolase (EC 3.2.1.91), β-xylosidase (EC 3.2.1.37), N-acetylglucosaminidase (EC 3.2.1.30), phosphodiesterase (EC 3.1.4.1), phosphomonoesterase (EC 3.1.3.2) and arylsulfatase (EC 3.1.6.1). Briefly, substrates (40 μl in dimethylsulfoxide) at final concentration of 500 mM were combined with three technical replicates of 200 µL of samples (diluted 1:10) in a 96 well plate. Background fluorescence was determined by combining 200 µL sample (diluted 1:10) with 40 μl MUF standards. The 96 well plates were incubated at 40 °C and fluorescence recorded at 5 and 125 min using a Tecan Infinite microplate reader (Tecan, Switzerland) with an excitation wavelength of 355 nm and an emission wavelength of 460 nm. Quantitative enzymatic activities were calculated after blank subtraction based on a standard curve of MUF. One unit of enzyme activity was defined as the amount of enzyme releasing 1 nmol of MUF min−1.

Laccase (phenoloxidase; EC 1.10.3.2) activity was determined by monitoring the oxidation of 2,2′-azino-bis(3-ethylbenzothiazoline-6-sulfonic acid) diammonium salt (ABTS) in citrate phosphate buffer (100 mM citrate, 200 mM phosphate, pH 5.0), by monitoring the formation of green colouration spectrophotometrically at 420 nm. Three technical replicates were performed for each sample.

Manganese peroxidase (MnP; EC 1.11.1.13) activity was determined by monitoring spectrophotometrically at 595 nm the purple colouration from oxidative coupling of 3-methyl-2-benzothialone-hydrazone hydrochloride (MBTH) and 3-(dimethyl amino)-benzoic acid (DMAB) in succinate-lactate buffer (100 mM, pH 4.5). Three technical replicates were performed for each sample. The results were corrected by activities of samples without manganese, and with ethylene diamine tetraacetate to chelate any Mn2+ present in the samples, allowing detection of Mn2+-independent peroxidases (versatile peroxidase). The results were also corrected by activities of samples in the absence of H2O2, allowing detection of oxidase (but not peroxidase) activity.

For each enzyme, n = 3 technical replicates were performed and enzyme activities were normalised to the protein content of each sample, which was determined using QubitTM fluorometric assays (Thermo Fisher Scientific Inc., UK).

Metabolomics analysis

As for enzyme assays, we selected blocks precolonised by T. versicolor (n = 5) upon which to perform UHPLC-MS (Supplementary Table 2). In total, 0.5 g of sawdust was added to 1666 μl of each of H2O, methanol and chloroform, vigorously vortexed and sonicated for 15 min (Elmasonic S150, Singen, Germany). The extracts were allowed to sit until the polar (containing H2O and methanol) and non-polar (containing mostly chloroform) layers separated, and we then removed 1500 μl of the upper layer containing the polar metabolite extracts. The extracts were centrifuged for 5 min at 17,000 × g (Biofuge, Thermo Fisher Scientific, MA, USA), and 200 μl of supernatant removed and dried in vacuo (Thermo Savant, NY, USA) for ca. 3 h. We then stored extracts at −80 °C until metabolomics analysis.

An intrastudy QC sample was prepared by pooling small aliquots of all samples, and the single extract was removed and dried down by centrifugation (20,000 × g for 10 min at 4 °C, Biofuge). UHPLC-MS-based metabolomics was performed (20 µl per sample) on a Thermo Dionex Ultimate 3000 RS system with a Thermo Scientific Q Exactive Orbitrap mass spectrometer. Samples (20 µl) were separated over a 100 × 2.1 mm, 1.9 µm particles, C18 column (Thermo Hypersil Gold) at a flow rate of 400 µl min−1 using a 14 min linear gradient programme from 0.1% formic acid in water to 0.1% formic acid in methanol. MS acquisition started at 0.1 min, with the flow up to 0.45 min directed towards waste. We acquired data in positive ion and profile mode from m/z 100–1000 at 70,000 resolution. Samples were analysed in a controlled randomised order, with the intrastudy QC sample repeatedly analysed equidistantly between the biological samples.

Metabolomics data processing

To process the UHPLC-MS data, the Thermo.raw data files in profile mode were converted into mzML format in centroid mode using MSConvert (Proteowizard 3.0.7665). Data were then aligned using an R (3.2.0) based XCMS and CAMERA script (both R packages [26, 27]) which resulted in a csv file intensity matrix (containing 9309 features, i.e., peaks in the mass spectra). This matrix was imported into MatLab and inserted into a direct infusion mass spectrometry SIMStitch workflow [28] where a blank filter of >2x sample over blank signal was applied, and a sample filter of peak presence of at least 50% of all samples [29]. The matrix was further processed by probabilistic quotient normalisation and subsequently missing values were imputed using K-nearest neighbour with k = 5. The imputed data matrix was used as an input to univariate statistics, including calculation of fold changes (FC). For multivariate statistics, a g-log transformation of the imputed data matrix was additionally applied, using an assessment of the technical variance across the repeated measurements of the intrastudy QC sample [29].

We putatively annotated UHPLC-MS features by inputting m/z values and their associated mean intensities into MI-Pack software version 2 beta [30], where metabolites were compared against the Kyoto Encyclopaedia of Genes and Genomes (KEGG) database. A set of highly significant metabolites (ANOVA: p < 0.001) were further searched against the KEGG database to determine possible roles within metabolic pathways.

Statistical analysis

To analyse the rate of decay and progression of interactions, we used R statistical software [31]. For each interaction, we assigned every species an individual score of combative ability, expressed as the percentage of the total system that it finally occupied. Briefly, each competitor scored between 0 and 2 for each block within a system (since two regions of every block were isolated; no outgrowth of a competitor from either isolation point scored 0, outgrowth from one isolation point scored 1, outgrowth from both isolation points scored 2). Scores for all blocks within a system were combined for each competitor individually, normalised to the number of replicates and converted to a percentage of the total system colonised. The data were analysed using a General Linear Model combined with Tukey post hoc tests, with individual block position (i.e., number of faces of that block involved in direct combat) and access to water (water was added to the perlite, as such, in 27-block interactions the layer laid on the perlite has greatest access to water) factored into the model. The rate of decay of wood in all interactions was compared using a one-way ANOVA followed by Tukey post hoc tests.

For enzyme activities, we used a one-way ANOVA with Tukey post hoc tests to compare mean activity (from five replicates), or Kruskal–Wallis tests followed by a Dunn’s test post hoc procedure when data were non-normally distributed, in R statistical software [31].

For GC-MS (VOCs) data, the entire data set was analysed by principal components analysis (PCA) to check for clustering of the QCs and, therefore, robustness of the data set (see Supplementary Fig. 1 for QC clustering), using MetaboAnalyst 3.0 [32]. We then removed QCs from the matrix, and orthogonal projection to latent structures-discriminant analysis (OPLS-DA) (Q2 = 0.93, R2 = 0.96), chosen for its cross-validation method which reduces false-positive results [33], was applied to the standardised binned data to determine the degree of separation between the four major sample groups: pair-wise samples 1 and 28 days after interaction set up, and 27-block samples 1 and 28 days after interaction set up. Next, we separated the data and applied OPLS-DA to pair-wise (Q2 = 0.57, R2 = 0.61) and 27-block (Q2 = 0.37, R2 = 0.55) sample groups separately. The modelled covariance and correlation were used to identify the features contributing most to the discriminant model separation, and one-way analysis of variance (ANOVA) with a 5% Benjamini–Hochberg false discovery rate (FDR) correction for multiple comparisons [34], and Tukey post hoc tests were applied to those features.

Lastly, for UHPLC-MS data we applied PCA to the g-log transformed data to explore the separation between control samples (T. versicolor growing alone), interaction samples and QCs. The median Relative Standard Deviation for the intrastudy QC samples was 11.15%, indicating that the MS data were of sufficiently high quality for further statistical analysis (see Supplementary Fig. 2 for QC clustering). After removing the QCs, ANOVA simultaneous components analysis (ASCA) was applied, and the model was permutation tested (5000 permutations) to determine the significance of factors (sample day, and block position within systems) [35, 36]. An additional ASCA model was tested to determine the significance of species distribution patterns within cubes, i.e., species being dispersed, V. comedens occupying a larger adjacent volume, or T. versicolor comprising the entire system, but did not include spatial location within assemblages as a factor. Pair-wise comparisons using ASCA were carried out for post hoc testing of significant effects, and the p values were adjusted for multiple testing using a 5% Bonferroni–Hochberg multiple testing correction [37]. Univariate ANOVA was applied to the whole normalised matrix with a 5% FDR correction [34] to test for significant metabolites. Finally, we determined FC between significant groups.

Network analytics

We investigated the synergy of metabolites produced by T. versicolor during all interactions by creating a co-occurrence Force Atlas2 [38] network analysis plot using Sci2 [39] and Gephi [40]. Data were filtered such that only significantly abundant features (ANOVA/ASCA: adjusted p < 0.05) relative to the baseline were included in the analysis, and abundancies of <10% the maximum were removed from the matrix, resulting in 908 retained variables. The weighted degree of nodes was calculated, and nodes were partitioned based on their weighting to facilitate removal of those which did not cluster into a discrete module. A final network was constructed from the refined data set, with edges weighted by count of occurrence and clusters coloured by weighted degree. The average abundance of all features within a cluster was calculated, and we used a one-way ANOVA with Tukey post hoc tests with a 5% FDR correction [34] to compare mean abundance between interacting systems after 28 days, in R statistical software [31].

Results

Hierarchy and coexistence dynamics

Typical transitive hierarchy of the focus decay species, T. versicolor, was established in paired block interactions (Fig. 1), i.e., H. fasciculare > T. versicolor, H. fasciculare > V. comedens, and T. versicolor > V. comedens (Fig. 3). This hierarchy reflected the general niche occupied in wood decomposition: least competitive V. comedens an early/primary decay species, T. versicolor a secondary colonising species, and H. fasciculare a later/tertiary decay species, the most competitive. Transitivity was also exhibited when the three fungi were dispersed throughout more spatially heterogeneous 3-dimensional systems (Fig. 1e): 15% of the original territory (defined as the relative proportion of a block occupied by a single fungus) of T. versicolor was captured by H. fasciculare, and 69% of the original territory of V. comedens was captured by its competitors. However, when spatial dynamics within the 3-dimensional cubes were changed such that V. comedens occupied a larger adjacent, but same total, volume as its competitors (Fig. 1f), V. comedens was displaced from just 6% of its original territory, and T. versicolor was displaced from 46% of its original territory. The three species coexisted without V. comedens being driven to near extinction (within a closed community network loop) as it was when territory was patchy, which is a characteristic of an intransitive relationship (Fig. 3).

Fig. 3: Interaction progression after 28 days in pair-wise and 3-dimensional interactions.
figure 3

Bars show the proportion of territory occupied by each species at the end of the experiment (mean of n = 5, SEM = ±12.2) in pair-wise interactions (T. versicolor (Tv) against V. comedens (Vc); T. versicolor against H. fasciculare (Hf); H. fasciculare against V. comedens), and 3-dimensional cubes in which all three species were evenly dispersed (Fig. 1e), and where V. comedens occupied a larger adjacent volume (Wall; Fig. 1f).

Full size image

Metabolite network function

In addition to our 12 targeted enzyme assays, our untargeted analyses yielded 67 VOCs, and 2825 LCMS signals of which 1597 were putatively annotated against existing compound libraries (Supplementary Table 3). Univariate ANOVA analysis (with FDR correction of p values) and ASCA multivariate analyses provided evidence of significant differences in the quantity of individual functional metabolites from blocks originally colonised by T. versicolor between community systems which showed divergent coexistence dynamics. Network analysis of the detectable metabolome, VOCs and enzymes throughout all interactions (average weighted degree = 3230), revealed 10 distinct clusters of synergistic metabolites (both antagonistic and lignocellulose decaying), plus an additional 11 independently functioning compounds with a high network degree. Namely, toluene (VOC: C66), an undefined oxidase enzyme, MnP, manganese independent peroxidase (peroxidase), arylsulfatase, phosphodiesterase, and five unidentified small metabolites (Fig. 4; cluster identity Table 1).

Fig. 4: Metabolic network of the full complement of significant compounds produced by T. versicolor throughout all interactions.
figure 4

Synergistic clusters of putatively annotated metabolites are clearly visualised and grouped based on their weighted degree (WD) (number of weighted edges, i.e., connections to other nodes, for a node. Average WD = 3230) and grey nodes labelled with compound names denote metabolites that did not form clusters (note that unlabelled grey nodes lack putative identifications). Clusters and their WD are detailed in Table 1. Edges are weighted by the count of occurrence of synergistic metabolites within samples, and node sizes and cluster colours represent weighted degree.

Full size image
Table 1 Details of a subset of putatively annotated metabolites clustered into 10 cluster sets plus unclustered metabolites, with the weighted degree (WD) of clusters given.
Full size table

By comparing the average abundance of individual clusters of metabolites produced in blocks precolonised by T. versicolor in pair-wise competitive systems, we found that production of seven clusters of antifungals was induced in significantly greater quantities when H. fasciculare was T. versicolor’s opponent, compared to three clusters when it was paired against the weaker V. comedens, relative to single species controls (Fig. 5). This particular response was consistent, and in 3-dimensional community interactions there was a stronger competition response shown in T. versicolor precolonised blocks when V. comedens was more combative when occupying a larger adjacent patch size (Clusters 6 and 9 were significantly more abundant (ANOVA: p < 0.05)), compared to when all fungi were dispersed in 3-dimensions. Cluster 6 contained the putatively annotated metabolite swainsonine, which functions in the biosynthesis pathway of piperidine- and pyridine-based antimicrobial alkaloids, and isolongifolene, a sesquiterpene with known antifungal properties, featured in Cluster 9 (Table 1). Although only putatively identified, the increased abundance of these combative/defensive compounds in T. versicolor-colonised blocks correlates with the increased combative ability of competitor V. comedens and longer coexistence of the three fungi within an intransitive relationship loop (Fig. 3).

Fig. 5: Significant differences in the average abundance of metabolic clusters (details of composition of clusters in Table 1) between interactions.
figure 5

Orange denotes cluster abundance is significantly higher (ANOVA: p < 0.05) in the interaction on the left of the colon; blue denotes significantly lower (ANOVA: p < 0.05) cluster abundance in the interaction on the left; and no colour denotes no significant difference (ANOVA: p > 0.05) between a pair of interactions. Tv, T. versicolor; Vc, V. comedens; Hf, H. fasciculare; Wall, 3-dimensional cube where V. comedens occupied a larger adjacent volume; Dispersed, 3-dimensional cube where all three species were dispersed; TvCube, 3-dimensional cube colonised entirely by T. versicolor (colour figure online).

Full size image

In addition to our findings in competitive systems, we found distinct metabolic differences between 2-dimensional and more spatially complex 3-dimensional controls where T. versicolor was the sole occupant of these resource habitats. For example, the unclustered enzyme arylsulfatase and Cluster 10, which comprised another six similarly functioning enzymes (Table 1), were significantly more highly abundant (ANOVA: p < 0.05) when T. versicolor solely occupied a 3-dimensional cube compared to when it was solely paired in 2-dimensions (Fig. 5).

The presence of other species affected metabolic function, as for example Clusters 7 and 8, amongst others, comprising compounds such as antifungal sesquiterpenes (Table 1), were significantly more abundant (ANOVA: p < 0.05) when T. versicolor was paired in an interaction compared to when it decayed wood alone. Similarly, when we compared whole T. versicolor 3-d cubes with 3-d community interactions, clusters 2, 3 and 4 which contain metabolites such as ankorine and fortimicin involved in antibiotic biosynthesis pathways, were significantly more abundant (ANOVA: p < 0.001) in the more species-rich systems. Within this experimental time frame the process of decomposition was not affected by spatial dynamics or species diversity (ANOVA: p > 0.05; Supplementary Fig. 3).

Discussion

Our results indicated that coexistence dynamics and metabolic function are directly affected by spatial occupation and patchiness of territory and can be translated into mechanistic functional processes. Furthermore, we provide evidence of a metabolic response to shifts in community structure as a result of altered connectivity. Landscape architecture changes combative mechanisms [2] and consequently may cause variation in expenditure of metabolic products between systems of varied levels of spatial heterogeneity. Community composition may be altered by a change in combative mechanisms as a result of landscape structural complexity [2], possibly due to stochastic effects within community assemblages [3] or different gaseous regimes throughout 3-dimensional structures and altered edge effects in 2-dimensions vs 3-dimensions. A study in which individuals of a community were paired against each other in artificial media [21] found negligible effects of species diversity alone on community function, but that a diverse community comprising weak competitors with high intransitivity exhibited a positive diversity–function relationship, i.e., the structure of a competitive network impacts community-level function. The more realistic complexity of the model system used in the present study revealed very different relationships: T. versicolor changed its mechanism of combat in systems where V. comedens occupied a larger adjacent volume, as clusters of metabolites functioning in the biosynthesis pathways of antagonistic antimicrobial compounds were produced in greater abundance. Presumably this emergent property was as a result of the greater combative strength of V. comedens (i.e., more antifungals were needed to attack the connected V. comedens), or as a result of a change of strategy by H. fasciculare which focused its attention on antagonising the now weakest competitor with whom it shared the most antagonistic fronts, T. versicolor, which responded to this change with increased antimicrobial compound production, or, the emergent property could have resulted from both simultaneously. Additionally, in this scenario H. fasciculare was able to capture some of the territory of T. versicolor which would result in a change to species-specific biochemical production. When V. comedens was dispersed it was less competitive than when connected, and connectivity of V. comedens altered the community dynamic such that the usually stronger competitor, T. versicolor, came under survival pressure. The change in mechanisms led to coexistence of the three species with more similar relative abundances and a closed community network loop (characteristic of intransitivity), compared to transitive community dynamics when the individuals were dispersed, where V. comedens was outcompeted to near extinction. It is worth noting that in natural dead wood species diversity would be greater than that presented here, and the presence of other fungal species as well as bacteria would influence interaction outcomes and, therefore, metabolism [41]. Some fungal–bacterial interactions are mutualistic [42], which could give individual fungal species a competitive advantage against their competitors, and further alter expected community hierarchies.

In the present study, production of seven substrate processing enzymes functioning in resource exploitation was boosted when spatial scale was increased from linear 2-dimensions to more heterogeneous 3-dimensions. The differences in functional metabolic processes across spatial scales, and between the more diverse community exhibiting intransitive characteristics and the transitive community is, therefore, reflective of a relationship between diversity and function that is regulated by space. The impact of the nature (space occupancy) of communities in directing community structure should be considered when looking at complex communities where outcomes are less predictable, and metabolic quantification to potentially inform predicted outcomes is, therefore, very useful.

That species diverse communities promote coexistence and positively impact community function is well known [21, 43,44,45]. Our results confirm this relationship between function and diversity, i.e., metabolic function changed significantly between systems with different numbers of species (single species assemblages vs multi-species assemblages), but we also show that spatial scale and distribution of species (patch dynamics) affect metabolic function as well. While these effects were not directly translated into ecosystem services (i.e., the rate of wood decay was not significantly affected), decay in the natural environment occurs over much larger time scales than used here [46]. So it might be predicted that given a greater length of time spatial heterogeneity and species diversity would have resulted in changes to the rate of decomposition, since substrate utilisation was affected over the short time scale measured in this study. The decomposition activities of wood decay fungi determine the rate of nutrient cycling in forest ecosystems which impacts forest function [47]. The relationship between spatial dynamics, species diversity and function highlighted in our study is, therefore, a key mechanism in the release of carbon from organic substrates into the carbon cycle, which drives global change [48]. Microbial communities in every global ecosystem carry out an array of functions as important as that of wood decayers [49,50,51], but the effects of spatial dynamics and species diversity on these functions have not previously been measured and quantified by experimental studies. The ecologically pertinent systems presented here are pioneering in their quantitative capture of 3-dimensional spatial dynamics into the experimental study of microbial, community and landscape structural ecology, which can be adapted, reconfigured and reimagined for the study of communities with a range of interaction types (e.g., neutral, mutualistic, facilitative). Our 3-dimensional experimental design, and the finding that spatial dynamics directly impact coexistence, diversity and function, are not only translatable to the understanding of diversity in existing microbiomes but may also provide key insights into extinctions and predictions of future ecological trends and community-level evolution.