2.1 Soil sampling and laboratory analyses

We used 577 soil samples from the Biomes of Australian Soil Environments (BASE) project (Bissett et al., 2016). In that project, sampling was undertaken from in that supports diverse plant communities across Australia. The sampling was carried out during the growing season when hydrothermal conditions are most conducive to typical plant growth. In the higher rainfall forested regions of the continent, the soil samples were collected mostly in spring and summer from September to February. In the shrublands and grasslands of the semi-arid and arid interior, soil samples were collected in spring from September to November. In the transitional zone between the south-east coast and the more arid interior, soil samples were collected in mainly autumn from March to May. Samples came from two soil depths (0–0.1 and 0.2–0.3 m), covering five typical Australian ecosystem types comprising cropland, forest, grassland, shrubland and woodland (Fig. 1a). Woodlands in Australia represent ecosystems which contain widely spaced trees, the crowns of which do not touch. Woodlands consist of areas with fewer and more scattered trees than forests. In temperate Australia, woodlands are mainly dominated by Eucalyptus species. Temperate woodlands occur predominantly in regions with a mean annual rainfall of between 250 and 800 mm, forming a transitional zone between the higher rainfall forested margins of the continent and the shrub and grasslands of the arid interior. Each sample was partitioned into subsamples for DNA sequencing (see below) and air-dried and crushed to a particle size of ≤2 mm for physicochemical analyses. The soil properties analysed were total organic carbon and soil nutrients (e.g. ammonium, nitrate, phosphorus, potassium), pH, exchangeable cations (aluminium, sodium, magnesium, calcium) and texture (sand, silt, clay). The methods are described in Bissett et al. (2016). Subsamples of the ≤ 2 mm portions were used for the spectroscopic analysis (see below).

Figure 1(a) Sampling sites and the range of ecosystem types across Australia. (b) The mean relative abundances of dominant fungal phyla and unclassified “Others” taxa in five ecosystem types. Individual abundance of each phylum and their cumulative abundance were shown in the graph.

2.2 Derivation of fungal abundance and diversity

The methods for DNA extraction and sequencing are detailed in Bissett et al. (2016). Briefly, the soil DNA was extracted in triplicate following methods used in the Earth Microbiome Project (http://www.Earthmicrobiome.Org/emp-standard-protocols/dna-extraction-protocol/, last access: 22 March 2022). Sequencing occurred with an Illumina MiSEQ, which is described in the BASE protocols (https://ccgapps.Com.Au/bpa-metadata/base/information, last access: 22 March 2022). Summarising, amplicons targeting the fungal ITS region were prepared and sequenced for each sample. The ITS amplicons were sequenced using 300 bp paired end sequencing. ITS1 regions were extracted using ITSx (Bengtsson-Palme et al., 2013). Sequences comprising full and partial ITS1 regions were passed to the operational taxonomic units (OTUs) selection and assigning workflow (Bissett et al., 2016), which followed guidelines described in the BASE protocols (https://ccgapps.com.au/bpa-metadata/base/information, last access: 22 March 2022). These are based on the most current version of UNITE database (version 8.2, updated 15 January 2020) for molecular identification of fungi (Nilsson et al., 2018). We used the final sample-by-OTU data matrix and annotated taxonomy file for the analyses of fungal diversity and composition.

To eliminate bias in the diversity comparison caused by unbalanced sequencing, samples were resampled at the same sequencing depth using functions of the RAM library in R software (R Core Team, 2014). The BASE dataset sought to produce as many sequences as resources allow with a minimum sequencing number of 10 000 per sample. Here, 11 000 sequences (the median number of sequences in the samples) were used as the resampling depth, because the majority of samples only had this amount of sequences, but also because the rarefaction curves started to flatten out for all 577 samples at this sequencing depth. This suggested that the sequencing number was sufficient (Fig. S2 in the Supplement). To quantify community diversity, we then calculated the abundance-based coverage estimator (ACE) index (Lozupone and Knight, 2008) from the resampled sample-by-OTU matrix. The relative abundance of fungal phyla were then determined using the ratio of the number of sequences classified at each phylum to the total number of sequences of each sample.

2.3 Soil visible–near-infrared spectroscopy

We measured the diffuse reflectance spectra of all air-dried ≤ 2 mm soil samples with the Labspec^® vis–NIR spectrometer (Malvern Panalytical, Boulder, Colorado, USA) following the protocols described in Viscarra Rossel et al. (2016). The spectral range of the spectrometer is 350 to 2500 nm. Due to the low signal-to-noise ratio at the start and end of each spectrum, for our analysis, we kept only spectra in the range between 380 and 2450 nm. As the spectra are highly collinear, to reduce redundancy in the data, we re-sampled them to a resolution of 10 nm. The measurements were performed with the instrument's high intensity contact probe (Malvern Panalytical, Boulder, Colorado, USA), and a Spectralon^® white reference panel was used for calibration once every 10 measurements.

For the modelling and interpretation, we first transformed the reflectance (R) spectra to apparent absorbance, using $A={\log }_{\mathrm{10}}(\mathrm{1}/R)$, and then used the Savitzky–Golay method with a window of size 7, a quadratic polynomial and first derivative method (Savitzky and Golay, 1964), in order to remove baseline effects and to improve the signal-to-noise ratio. To visualise the spectra, we further fitted each reflectance (R) spectrum with a convex hull and computed the deviations from the hull (Clark and Roush, 1984). These continuum removed (CR) spectra help to visualise the characteristic absorptions more clearly than the Savitzky–Golay first derivatives (SG1Der) absorbance spectra.

2.4 Modelling soil fungal abundance and diversity

We developed spectroscopic models, environmental models and spectrotransfer functions for estimating soil fungal abundance and diversity (see below). The spectroscopic models used only the vis–NIR spectra, the environmental models used only the publicly available soil and environmental data that represent the soil forming factors soil, climate, vegetation, terrain and parent material (Jenny, 1994) and the spectrotransfer functions used the vis–NIR spectra together with soil and environmental data.

We assembled a set of readily available soil and environmental maps that represented climate, terrain, vegetation and parent material. To relate the these covariates to the fungal data, we extracted values from these maps using the geographic coordinates of the sample set. The soil property data came from Australia-wide fine spatial resolution (90 × 90 m) digital soil maps of total organic carbon, total nitrogen, total phosphorus, bulk density, effective cation exchange capacity, available water capacity, pH and soil texture (sand, silt and clay; Viscarra Rossel et al., 2015), as well as maps of the clay minerals kaolinite, illite and smectite (Viscarra Rossel, 2011). To represent climate, we used data on mean annual temperature (MAT), mean annual precipitation (MAP), solar radiation and evapotranspiration (Xu and Hutchinson, 2011) and the Prescott index (PI; Prescott, 1950), which is calculated as the ratio of precipitation to evapotranspiration. To capture functional landscape characteristics, we used a digital elevation model (DEM) from the 3 arcsec Shuttle Radar Topography Mission (SRTM) and derived terrain attributes (Gallant et al., 2011). To represent vegetation, we used data on net primary productivity (NPP; Haverd et al., 2013) and on the fraction of photosynthetically active radiation intercepted by the sunlit canopy of the evergreen (Fpar-e) and woody (Fpar-r) vegetation (Donohue et al., 2009). To represent parent material, we used gamma radiometrics, which comprises data on potassium, uranium and thorium (Minty et al., 2009). Supplement Table S2 lists these data and their main characteristics.

The spectra and the covariates were centred and scaled before the modelling of fungal abundance and diversity. The algorithms that we tested were PLSR (Wold et al., 2001), Gaussian process regression (GPR; Rasmussen and Williams, 2005, support vector machines (SVM; Suykens et al., 2002), random forest (RF; Breiman, 2001), Cubist (Quinlan, 1992), extreme gradient boosting (XGBoost; Friedman, 2001) and optimised one-dimensional convolutional neural networks (1D-CNNs; Shen and Viscarra Rossel, 2021). The algorithms and their implementation are described in the Supplement linked to this article.

The predictability of the spectroscopic models and the spectrotransfer functions were assessed using 10-fold cross-validation. We evaluated the estimates using the Nash–Sutcliffe model efficiency, otherwise known as the coefficient of determination (R²), which represent the fraction of the explained variance based on the 1:1 line of estimated versus measured values (Janssen and Heuberger, 1995). The R² was computed as 1-RSS/TSS, where RSS is the residual sum of squares and TSS is the total sum of squares. The root mean squared error (RMSE) measures inaccuracy, the standard deviation of the error (SDE) measures imprecision and the mean error (ME) measures bias (Viscarra Rossel and McBratney, 1998). Inaccuracy (RMSE) embraces both the bias (ME) and the imprecision (SDE; Viscarra Rossel and McBratney, 1998). Their relationship is given by ${\mathrm{RMSE}}^{\mathrm{2}}={\mathrm{ME}}^{\mathrm{2}}+{\mathrm{SDE}}^{\mathrm{2}}$.

To interpret the models, we calculated their variable importance as follows. For the PLSR, GPR, SVM, Cubist, RF and XGBoost models, variable importance was calculated using the varImp function in the caret library (Kuhn et al., 2008) of the software R. To calculate the variable importance of the CNN models, we used permutation variable importance. In our case, we ran 1000 permutations and measured the decrease in RMSE after a predictor was permuted (randomly rearranged). The permutation breaks the relationship between the predictor and the response variables, and a reduction in RMSE indicates how much the model depends on the particular predictor. An advantage of this approach is that it can be applied on any estimator and does not require retraining the model (Breiman, 2001; Fisher et al., 2019). In order to compare the importance between different fungal phyla and diversity, we scaled the importance values between 0 and 1. In the results, we only report the variable importance of the model that performed best.

3.1 Modelling

With the different algorithms, the spectroscopic models (i.e. with only the vis–NIR spectra) could explain 9 %–45 % of the variation in fungal phyla relative abundance and diversity. Spectroscopic models of the Glomeromycota were the least successful, with R² values ranging from 0.09 using SVM to 0.30 using 1D-CNN, while those of the Mortierellomycota produced the largest R² values, ranging from 0.32 using XGBoost to 0.45 using 1D-CNN (Fig. 3). The models of diversity had R² values ranging from 0.14 with PLSR to 0.35 using 1D-CNN.

The models derived with the readily available soil and environment data could explain 14 %–60 % of the variation in fungal phyla relative abundance and diversity with the different algorithms. These environmental models generally performed better than the spectroscopic models, with an average 10 % additional variance explained.

Combining the vis–NIR spectra and soil and environmental data further improved the models and their explanatory power. The spectrotransfer functions (i.e. with the combined set of vis–NIR spectra and other soil and environmental data) performed, on average, 20 % better than the spectroscopic models and 10 % better than environmental models. Depending on the algorithm used, they could explain between 17 % and 73 % of the variation in fungal phyla relative abundance and diversity (Fig. 3). The spectrotransfer functions of Glomeromycota produced R² values, ranging from 0.17 using PLSR to 0.48 using 1D-CNN. The spectrotransfer functions of the Mortierellomycota and Mucoromycota produced the largest R² values ranging from 0.51 to 0.73 (Fig. 3).

Generally, PLSR and GPR were the least successful methods, while SVM, RF, Cubist and XGBoost were similarly successful at estimating fungal phyla relative abundance and diversity (Fig. 3). The 1D-CNN spectrotransfer functions were 13 %–31 % more successful compared to other machine learning methods as they could explain between 45 % and 73 % of the variation in fungal relative abundance and diversity (Fig. 3).

Figure 3Coefficient of determination (R²) for the vis–NIR spectroscopic models, soil and environmental models and the spectrotransfer functions that used a combined set of the vis–NIR and readily available soil and environmental covariates used to estimate soil fungal phyla abundance and diversity (n=577). The different statistical and machine learning methods were PLSR, GPR, SVM, RF, Cubist, XGBoost and optimised 1D-CNNs.

Download

3.2 1D-CNNs spectrotransfer functions

The final architecture and optimised hyperparameters of the 1D-CNNs are given in Supplement Table S3. As deep learning models are dataset dependent, the optimisation returned a different architecture for each response variable. Overall, the 1D-CNNs used simple architecture with less than four convolutional layers (Supplement Table S3). Scatter plots of the measured versus estimated values of relative abundance and diversity using 1D-CNNs spectrotransfer functions and their validation statistics are shown in Fig. 4. Estimates of the relative abundance of Mortierellomycota and Mucoromycota produced R² values ≥0.60, while estimates of Ascomycota and Basidiomycota produced values of $\mathrm{0.5}\le R^{\mathrm{2}}<\mathrm{0.6}$. Estimates of Glomeromycota and ACE produced values of $\mathrm{0.4}\le R^{\mathrm{2}}<\mathrm{0.5}$. The estimates were relatively unbiased (small ME), although generally small values were overestimated and large values were underestimated (Fig. 4). Imprecision contributed to the majority of the RMSE. The imprecision of our estimates was a result of the absence of repeated sampling and the high adaptability of soil fungi to the wide range of environments.

Figure 4Performance of the CNN spectrotransfer functions for estimation of the relative abundance of dominant fungal phyla and diversity index. The spectrotransfer functions used vis–NIR spectra with other publicly available data on soil environmental variables. The plots show measured vs. estimated values using a 10-fold cross validation. The grey points represent no overlap with any other points, and the black points represent at least two points that overlap.

Download

The important variables in the 1D-CNNs spectrotransfer functions of phyla relative abundance and diversity were vis–NIR wavelengths representing organic matter, iron oxide and clay minerals (Fig. 5).

Figure 5Important predictors of relative abundance of fungal phyla and diversity index measured by the variable importance of the 1D-CNNs spectrotransfer functions (n=577) derived with publicly accessible data that represent soil (S), climate (C), vegetation (V), terrain (T), parent material (PM) and vis–NIR spectra. The dots in red, orange and blue indicated the most, medium and least important levels. The importance value for the majority of wavelengths were low and close to zero, thus these wavelengths were not shown to make the figure clearer.

Download

The identified wavelengths mostly coincided with absorptions that are related to carbon functional groups found in organic matter, including C-H, N-H and C-O, with a smaller number of wavelengths coinciding with those that are related to clay minerals and iron oxides (Table 2). The organic functional groups, C-H alkyl and methyls, N-H of amines and C-O of carbohydrates, which might indicate the presence of relatively labile forms of carbon, were important in the models of fungal phyla but not of ACE diversity. The C=O of amides and carboxylic acids, which represent stable forms of carbon were not as important in modelling (Fig. 5). Other wavelengths that represent iron oxides and clay minerals were also important in the models, indicating the different ecological niches and physiological characteristics (Table 2).

Table 2Absorption wavelength assignment (in nanometres) for the most important vis–NIR wavelengths in the 1D-CNN models. The assignment of vis–NIR absorptions are based on Viscarra Rossel and Behrens (2010) and Stenberg et al. (2010).

Download Print Version | Download XLSX

Other soil properties, such as total organic carbon and pH were important variables in the spectrotransfer functions of Ascomycota and Basidiomycota, and fungal diversity. Total organic carbon and total nitrogen were important in the spectrotransfer functions of Mortierellomycota and Mucoromycota and bulk density was important in the spectrotransfer functions of Glomeromycota, Ascomycota and ACE diversity (Fig. 5). As well as soil properties, climatic factors such as PI and PET, and vegetation, represented by Fpar-e and NPP were also important in the modelling of fungal phyla relative abundance and community diversity. The variables that we used to represent terrain and the parent material exerted less influence in the models (Fig. 5).