Below-ground Fungal Biodiversity Can be Monitored Using Self-Supervised Learning Satellite Features

Robin Young; Michael E. Van Nuland; E. Toby Kiers; Tomáš Větrovský; Petr Kohout; Petr Baldrian; Srinivasan Keshav

cs.LGApr 10, 2026

Below-ground Fungal Biodiversity Can be Monitored Using Self-Supervised Learning Satellite Features

Authors: Robin Young, Michael E. Van Nuland, E. Toby Kiers, Tomáš Větrovský, Petr Kohout, Petr Baldrian, Srinivasan Keshav

Organizations: Department of Computer Science and Technology, University of Cambridge, Cambridge, Cambridgeshire, UK. · Society for the Protection of Underground Networks, Dover, DE, US. · Amsterdam Institute for Life and Environment (A-LIFE), Section Ecology & Evolution, Vrije Universiteit Amsterdam, Amsterdam, The Netherlands. · Institute of Microbiology, Czech Academy of Sciences, Videnska 1083, Prague, 14200, Czech Republic. ·

Abstract

Mycorrhizal fungi are vital to terrestrial ecosystem functioning. Yet monitoring their biodiversity at landscape scales is often unfeasible due to time and cost constraints. Current predictions suggest that 90% of mycorrhizal diversity hotspots remain unprotected, opening questions of how to broadly and effectively map underground fungal communities. We show that self-supervised learning (SSL) applied to satellite imagery can predict below-ground ectomycorrhizal fungal richness across diverse environments. Our models explain over half the variance in species richness across ~12,000 field samples spanning Europe and Asia. SSL-derived features are the most informative tested predictor group, and outperform each of the established climate, soil, and land cover baselines. We achieve a 10,000-fold increase in spatial resolution over existing techniques, moving from 1km landscape averages to 10m habitat-scale observations. As satellite observations are dynamic rather than static, this enables temporal monitoring of below-ground biodiversity at landscape scales for the first time. We apply this approach across two UK National Park woodlands, where ancient forests are predicted to have high but declining ectomycorrhizal diversity, marking key areas for further field verification. These results establish SSL satellite features as a scalable tool for extending sparse field observations to continuous, high-resolution biodiversity maps for monitoring the invisible half of terrestrial ecosystems.

Figures & tables

Figure 1: Overview of the satellite-to-fungi biodiversity prediction framework. Ground-truth ectomycorrhizal fungal richness estimates, derived from metabarcoding of soil eDNA samples across Europe and Asia, are paired with co-located satellite time series from Sentinel-1 (radar) and Sentinel-2 (optical) imagery. A self-supervised learning foundation model (Tessera) processes the multi-temporal, multi-modal satellite data to produce high-dimensional embeddings that encode vegetation dynamics, phenology, and environmental conditions without requiring hand-crafted features. These learned representations, along with optional traditional environmental predictors (climate, soil, land cover), are used to train a gradient boosting model to predict fungal species richness. The output is a continuous, high-resolution prediction of below-ground fungal species richness at habitat scale (10m).

Figure 2: Predictive performance and feature importance across models and feature sets. (a) Out-of-sample

R^{2}

for each feature set, by model. Points are the mean over 50 runs, error bars are 95% confidence intervals. (b) Share of total feature importance carried by each feature group, as mean

\pm

1 standard deviation over 50 runs.

Figure 3: Visual comparison of model predictions for fungal richness in a mountainous landscape. (A) Satellite basemap provides the landscape context. (B) The high-resolution (10m) prediction from the SSL model with water or glacial areas masked. (C) The coarse-resolution (1km) prediction from the climate-only model.

Figure 4: Map showing geographic distribution of absolute prediction errors of satellite SSL only model, with point sizes and colors representing error magnitude

Figure 5: Modeled baseline ectomycorrhizal (EcM) fungal richness. Maps show predicted baseline EcM fungal richness (2017–2019 mean) for woodland areas within the Lake District (top) and Cairngorms (bottom) National Parks. White outlines delineate designated Ancient Woodland boundaries. Distribution plots (right) show the predicted baseline richness across all forested pixels, separated by historical forest type. Grey bars indicate the central 99%, 75%, 50%, 25%, and 10% of pixels (light to dark) and white circles show the median.

Figure 6: Example spatio-temporal analysis of predicted EcM fungal richness trajectories. Maps display predicted EcM forest categories for woodland areas within the Lake District (top) and Cairngorms (bottom) National Parks. White outlines delineate designated Ancient Woodland boundaries. Forested pixels were classified into five categories based on their 2017–2019 baseline richness percentiles and 8-year regression trends in predicted richness: High-stable, High-decreasing, Low-increasing, Low-decreasing, and Intermediate/stable. Line plots track the 2017–2024 annual mean predicted richness (solid lines) and linear trends (dashed lines) for each category, with grey shading indicating values below the 25th and above the 75th baseline richness percentiles for a given park. Stacked bar charts show the proportional extent of each category, separated by historical forest type.

Setting	Value	Notes
Representation
Embedding dimension	128	per 10m pixel, per year
Encoder width $d_{model}$	512	per modality branch
Transformer blocks	8	per modality branch, feed-forward width 4,096
Modality fusion	concatenation	Sentinel-1 and Sentinel-2 encoder outputs
Projector hidden width	16,384	discarded after pretraining

Table S1: Configuration of the trained Tessera checkpoint used in this work. Pretraining settings are those of the released Tessera v1.0 model (Table 8 of Feng et al., 2026 ); the remaining rows describe our feature-extraction pipeline.

Category	Source	Resolution	Variables & processing
Climate	WorldClim 2.1 ( Fick and Hijmans, 2017 )	$\sim$ 1km	19 bioclimatic variables (BIO1–BIO19); annual mean, min, and max of monthly tavg , tmin , tmax , prec , srad , vapr , wind
Topography	WorldClim 2.1 DEM ( Fick and Hijmans, 2017 )	$\sim$ 1km	Elevation (ma.s.l.), slope (°), and aspect (°) derived via richdem ( Barnes, 2016 )
Soil properties	SoilGrids 2.0 ( Poggio et al., 2021 )	250m	pH, organic carbon, nitrogen, clay/sand/silt fractions, CEC, and bulk density at 0–5cm and 5–15cm depths; WRB soil class (one-hot encoded)
Land cover	ESA WorldCover 2021 ( Zanaga et al., 2022 )	10m	11 land cover classes (one-hot encoded)
Geography	—	—	Latitude and longitude (continuous)

Table S2: Summary of environmental predictor variables used as baseline features for comparison with SSL satellite features.

Component	PC1	PC2	PC3	PC4	PC5	PC6	PC7	PC8	PC9	PC10
Explained variance (%)	14.4	12.0	7.0	6.3	4.8	3.4	3.1	2.9	2.5	2.3
Cumulative (%)	14.4	26.4	33.4	39.7	44.4	47.8	50.9	53.8	56.3	58.6

Table S3: Variance explained by the leading principal components of the SSL embedding.

Variable	PC1	PC2	PC3	PC4	PC5	PC6	PC7	PC8	PC9	PC10
Topography
Elevation	0.16	0.38	0.03	$-$ 0.05	$-$ 0.50	0.09	0.30	$-$ 0.37	0.14	0.04
Temperature
Bio 1 (Ann. Mean Temp)	$-$ 0.30	$-$ 0.22	$-$ 0.76	0.14	0.00 †	0.05	0.16	0.07	0.00 †	0.06
Bio 2 (Mean Diurnal Range)	0.20	$-$ 0.05	$-$ 0.16	$-$ 0.05	$-$ 0.02 †	$-$ 0.02 †	0.39	$-$ 0.10	0.03	0.05
Bio 3 (Isothermality)	$-$ 0.03	0.09	$-$ 0.61	$-$ 0.03	$-$ 0.24	$-$ 0.04	0.27	$-$ 0.21	0.04	0.19

Table S4: Pearson correlation between the leading ten principal components of the SSL embedding and every environmental variable. Correlations with

|r|\geq 0.30

(before rounding) are in bold. Correlations whose 95% bootstrap confidence interval (2000 bootstraps) includes zero are marked † . All others are significant at

\alpha=0.05

Component	Strongest correlates	$r$	95% CI	$p$
PC1	Bulk density (0–5cm)	$-0.381$	$[-0.400,-0.362]$	$<10^{-100}$
	Clay content (5–15cm)	$-0.378$	$[-0.397,-0.357]$	$<10^{-100}$
	Clay content (0–5cm)	$-0.373$	$[-0.393,-0.353]$	$<10^{-100}$
PC2	Bulk density (depth gradient)	$+0.444$	$[+0.428,+0.460]$	$<10^{-100}$
	Vapor Pr. (SD)	$-0.410$	$[-0.432,-0.390]$	$<10^{-100}$
	Elevation	$+0.377$	$[+0.357,+0.398]$	$<10^{-100}$

Table S5: Three strongest environmental correlates of each principal component, with 95% bootstrap confidence intervals.

Scheme	$R^{2}$	MAE	Spearman $\rho$	Median train–test separation
Random	0.507 $\pm$ 0.018	29.9	0.765	0m
Grouped, 25m	0.316 $\pm$ 0.047	40.0	0.612	2,187m
Block, 25m	0.296 $\pm$ 0.057	40.8	0.610	1,920m
Block, 50m	0.298 $\pm$ 0.064	41.1	0.614	2,035m

Table S6: Performance under random and spatially explicit cross-validation. LightGBM trained on all samples, five folds

\times

three seeds. The last column is the median great-circle distance from each test sample to its nearest training sample.

	2022	2023	2024
Within-year $R^{2}$	0.4905 $\pm$ 0.0243	0.4995 $\pm$ 0.0099	0.5000 $\pm$ 0.0172
Within-year MAE	33.58	33.59	33.36
Mean cosine similarity to 2024 embeddings	0.872	0.892	—

Table S7: Effect of embedding year. Satellite-only, LightGBM, five seeds, on European samples. Within-year trains and tests on the same year.

Window	Features	$R^{2}$	MAE	Spearman $\rho$	$\Delta$ MAE vs 3 $\times$ 3 (95% CI)
1 $\times$ 1	128	0.219 $\pm$ 0.057	40.78	0.652	$+2.87$ $[+1.43,+4.27]$
3 $\times$ 3	1,152	0.303 $\pm$ 0.034	37.91	0.680	—
5 $\times$ 5	3,200	0.320 $\pm$ 0.020	37.32	0.687	$-0.59$ $[-1.69,+0.53]$
7 $\times$ 7	6,272	0.340 $\pm$ 0.015	36.90	0.700	$-1.01$ $[-2.23,+0.20]$

Table S8: Patch-size ablation. All sizes are center-crops of the same 7

\times

7 patch, so the samples and folds are identical across rows and comparisons are paired. 784 locations across five seeds.

\Delta

MAE is against the 3

\times

3 window used throughout the paper.

Family	$n$	Definition
Sentinel-2 band statistics	50	Mean, median, 10th percentile, 90th percentile and standard deviation of each of the ten bands (B2, B3, B4, B5, B6, B7, B8, B8A, B11, B12)
Vegetation indices	24	NDVI, EVI, NDWI and NBR. For each, same five statistics and a seasonal amplitude term defined as the 90th minus the 10th percentile
NDVI phenology	4	Day of year of the annual maximum and of the annual minimum, the integrated annual NDVI, and the mean absolute successive difference between consecutive observations as a measure of trajectory roughness
Sentinel-1 backscatter	12	Mean, median and standard deviation of VV and VH, computed separately for ascending and descending orbits

Table S9: Composition of the spectral-index baseline, 90 features per pixel. Statistics are computed over the cloud-free observations of a single year, the same year from which the embeddings are extracted.

Feature set	$R^{2}$	MAE
Tessera embeddings	0.5094 $\pm$ 0.0154	30.06
Spectral indices and statistics	0.4866 $\pm$ 0.0230	30.89

Table S10: Learned embeddings against a spectral-index baseline computed from the same imagery at the same resolution. LightGBM, ten seeds.

	Transferred ( $R^{2}$ )		Transferred ( $\rho$ )
Feature set	eur $\rightarrow$ asia	asia $\rightarrow$ eur	eur $\rightarrow$ asia	asia $\rightarrow$ eur
Satellite	$-9.8\pm 2.3$	$-0.051\pm 0.045$	$+0.04\pm 0.05$	$+0.10\pm 0.07$
Climate	$-8.2\pm 1.9$	$-0.221\pm 0.038$	$+0.14\pm 0.02$	$-0.17\pm 0.11$
Soil	$-13.3\pm 3.4$	$-0.023\pm 0.029$	$-0.06\pm 0.06$	$+0.22\pm 0.03$
Land cover	$-14.5\pm 1.9$	$-0.280\pm 0.085$	$-0.01\pm 0.06$	$+0.17\pm 0.04$

Table S11: Cross-continental transfer for each feature set. all samples, LightGBM, three seeds.

Quartile of predicted richness	$n$	Coverage	Mean width	Mean absolute error
1 (lowest)	2,961	0.837	60.7	6.7
2	2,961	0.863	77.4	18.4
3	2,961	0.870	131.6	34.7
4 (highest)	2,961	0.842	202.4	71.5

Table S12: Interval behavior by quartile of predicted richness.

Figure S1: Diagnosing the impact of prediction outliers on model performance. This sensitivity analysis plots the model’s mean $R^{2}$ ( $y$ -axis) as a function of the percentage of the test set retained after progressively removing the highest-error points ( $x$ -axis). Data are from the best-performing model (LightGBM, all features) over 50 runs; the shaded area shows the standard deviation.

Diagnostic (top decile)	Lift
Highest predicted richness	1.94
Lowest sampling density within 1km	1.67
Greatest novelty in embedding space	1.65
Greatest distance to nearest training sample	1.60

Table S13: Share of total absolute error captured by flagging the top decile under each diagnostic, expressed as a multiple of the 10% expected by chance.

Full Models (4 Feature Sets)
Data Source(s)	Model	$R^{2}$	RMSE	MAE	ME (Bias)
Sat + Clim + Soil + WC	LightGBM	0.550 $\pm$ 0.017	46.22 $\pm$ 1.54	28.29 $\pm$ 0.70	-0.07 $\pm$ 0.97
	XGBoost	0.545 $\pm$ 0.017	46.47 $\pm$ 1.58	28.55 $\pm$ 0.76	-0.09 $\pm$ 0.97
	Random Forest	0.529 $\pm$ 0.018	47.27 $\pm$ 1.50	28.34 $\pm$ 0.77	1.07 $\pm$ 0.91
Single-Feature-Set Models (Ablations)
Satellite only	LightGBM	0.535 $\pm$ 0.018	47.01 $\pm$ 1.61	29.07 $\pm$ 0.78	0.03 $\pm$ 0.96

Table S14: Performance of models and feature set ablations for predicting fungal biodiversity. Metrics are mean

\pm

standard deviation over 50 runs varying the random seed of the data split. The best value in each column is in bold. Data sources: Sat = satellite SSL features; Clim = WorldClim; Soil = SoilGrids; WC = ESA WorldCover. Metrics:

R^{2}

= coefficient of determination; RMSE = root mean squared error; MAE = mean absolute error; ME = mean error (bias).

Baseline	Model	Mean difference	95% CI	Wilcoxon $p$
Climate only	LightGBM	$+$ 0.0174	$[+0.0110,+0.0238]$	$8.2\times 10^{-7}$
	XGBoost	$+$ 0.0153	$[+0.0082,+0.0223]$	$9.0\times 10^{-5}$
	Random Forest	$+$ 0.0185	$[+0.0106,+0.0265]$	$4.0\times 10^{-5}$
Soil only	LightGBM	$+$ 0.0204	$[+0.0137,+0.0270]$	$2.0\times 10^{-7}$
	XGBoost	$+$ 0.0232	$[+0.0163,+0.0300]$	$4.2\times 10^{-8}$
	Random Forest	$+$ 0.0157	$[+0.0080,+0.0233]$	$1.7\times 10^{-4}$

Table S15: Paired comparison of the satellite-only configuration against each single environmental baseline. Differences are satellite minus baseline, in

R^{2}

, over 50 shared seeds.

Biodiversity Monitoring

See all 23 papers

ECHOSAT: Estimating Canopy Height Over Space And Time
Feb 24, 2026Jan Pauls, Karsten Schrödter, Sven Ligensa +7Remote Sensing Remote Sensing Change Detection
Department of Information Systems, University of Münster, Germany · Laboratoire des Sciences du Climat et de l’Environnement (LSCE), France · Zuse Institute Berlin (ZIB), Germany +2
Multi-Modal, Multi-Temporal Geospatial Foundation Models Improve Global Above-Ground Biomass Estimation
Aug 5, 2026Ghjulia Sialelli, Linus Scheibenreif, Jan Dirk Wegner +1Geospatial Foundation Models Forest Aboveground Biomass Estimation
Photogrammetry and Remote Sensing, ETH Zurich, Zurich, 8049, Switzerland · ETH AI Center, Zurich, 8092, Switzerland · EcoVision Lab, Department of Mathematical Modeling and Machine Learning, University of Zurich, Zurich, 8057, Switzerland
Integrating national forest inventory, airborne lidar, and satellite imagery for wall-to-wall mapping of forest structure with computer vision
Jun 18, 2026Luke J. Zachmann, David D. Diaz, Vincent A. Landau +8Remote Sensing Forest Aboveground Biomass Estimation
Vibrant Planet Public Benefit Corporation, Truckee, CA, USA

Feb 24, 2026J/K move · Enter open · S save

ECHOSAT: Estimating Canopy Height Over Space And Time

Jan Pauls, Karsten Schrödter, Sven Ligensa, Martin Schwartz+6

Department of Information Systems, University of Münster, Germany · Laboratoire des Sciences du Climat et de l’Environnement (LSCE), France · Zuse Institute Berlin (ZIB), Germany+2

Remote Sensing Remote Sensing Change Detection

Open paper AskPDF arXiv

Forest monitoring is critical for climate change mitigation. However, existing global tree height maps provide only static snapshots and do not capture temporal forest dynamics, which are essential for accurate carbon accounting. We introduce ECHOSAT, a global and temporally consistent tree height map at 10 m resolution spanning multiple years. To this end, we resort to multi-sensor satellite data to train a specialized vision transformer model, which performs pixel-level temporal regression. A self-supervised growth loss regularizes the predictions to follow growth curves that are in line with natural tree development, including gradual height increases over time, but also abrupt declines due to forest loss events such as fires. Our experimental evaluation shows that our model improves state-of-the-art accuracies in the context of single-year predictions. We also provide the first global-scale height map that accurately quantifies tree growth and disturbances over time. We expect ECHOSAT to advance global efforts in carbon monitoring and disturbance assessment. The maps can be accessed at https://github.com/ai4forest/echosat.

Abstract

Biodiversity Monitoring

Explore similar work

ECHOSAT: Estimating Canopy Height Over Space And Time

Multi-Modal, Multi-Temporal Geospatial Foundation Models Improve Global Above-Ground Biomass Estimation

Integrating national forest inventory, airborne lidar, and satellite imagery for wall-to-wall mapping of forest structure with computer vision