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. · Laboratory of Microbial Ecology and Biogeography, Institute of Microbiology, Czech Academy of Sciences, Videnska 1083, Prague, 14200, Czech Republic.
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 R2 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 ± 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 dmodel
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 )
∼ 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 )
∼ 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∣≥0.30 (before rounding) are in bold. Correlations whose 95% bootstrap confidence interval (2000 bootstraps) includes zero are marked † . All others are significant at α=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
R2
MAE
Spearman ρ
Median train–test separation
Random
0.507 ± 0.018
29.9
0.765
0m
Grouped, 25m
0.316 ± 0.047
40.0
0.612
2,187m
Block, 25m
0.296 ± 0.057
40.8
0.610
1,920m
Block, 50m
0.298 ± 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 × 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 R2
0.4905 ± 0.0243
0.4995 ± 0.0099
0.5000 ± 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
R2
MAE
Spearman ρ
Δ MAE vs 3 × 3 (95% CI)
1 × 1
128
0.219 ± 0.057
40.78
0.652
+2.87[+1.43,+4.27]
3 × 3
1,152
0.303 ± 0.034
37.91
0.680
—
5 × 5
3,200
0.320 ± 0.020
37.32
0.687
−0.59[−1.69,+0.53]
7 × 7
6,272
0.340 ± 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 × 7 patch, so the samples and folds are identical across rows and comparisons are paired. 784 locations across five seeds. Δ MAE is against the 3 × 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
R2
MAE
Tessera embeddings
0.5094 ± 0.0154
30.06
Spectral indices and statistics
0.4866 ± 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 ( R2 )
Transferred ( ρ )
Feature set
eur → asia
asia → eur
eur → asia
asia → eur
Satellite
−9.8±2.3
−0.051±0.045
+0.04±0.05
+0.10±0.07
Climate
−8.2±1.9
−0.221±0.038
+0.14±0.02
−0.17±0.11
Soil
−13.3±3.4
−0.023±0.029
−0.06±0.06
+0.22±0.03
Land cover
−14.5±1.9
−0.280±0.085
−0.01±0.06
+0.17±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 R2 ( 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.
Data Source(s)
Model
R2
RMSE
MAE
ME (Bias)
Full Models (4 Feature Sets)
Sat + Clim + Soil + WC
LightGBM
0.550 ± 0.017
46.22 ± 1.54
28.29 ± 0.70
-0.07 ± 0.97
XGBoost
0.545 ± 0.017
46.47 ± 1.58
28.55 ± 0.76
-0.09 ± 0.97
Random Forest
0.529 ± 0.018
47.27 ± 1.50
28.34 ± 0.77
1.07 ± 0.91
Single-Feature-Set Models (Ablations)
Satellite only
LightGBM
0.535 ± 0.018
47.01 ± 1.61
29.07 ± 0.78
0.03 ± 0.96
Table S14: Performance of models and feature set ablations for predicting fungal biodiversity. Metrics are mean ± 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: R2 = 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×10−7
XGBoost
+ 0.0153
[+0.0082,+0.0223]
9.0×10−5
Random Forest
+ 0.0185
[+0.0106,+0.0265]
4.0×10−5
Soil only
LightGBM
+ 0.0204
[+0.0137,+0.0270]
2.0×10−7
XGBoost
+ 0.0232
[+0.0163,+0.0300]
4.2×10−8
Random Forest
+ 0.0157
[+0.0080,+0.0233]
1.7×10−4
Table S15: Paired comparison of the satellite-only configuration against each single environmental baseline. Differences are satellite minus baseline, in R2 , over 50 shared seeds.
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.
Jan Pauls, Karsten Schrödter, Sven Ligensa +7
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
Accurate estimation of above-ground biomass (AGB) from satellite imagery is essential for the large-scale monitoring of carbon stocks, yet it remains a challenging regression task at global scale. Geospatial Foundation Models (GFMs) have recently emerged as a promising machine learning paradigm to derive general-purpose representations from Earth observation (EO) data, but their utility for quantitative regression tasks like biomass estimation remains largely unexplored, as most benchmarks focus on classification and segmentation. Here, we benchmark 13 GFMs for above-ground biomass estimation, specifically GEDI L4A biomass values. We compare them against a fully supervised state-of-the-art (SOTA) baseline trained on A Global-scale Biomass Dataset (AGBD), a machine learning-ready benchmark spanning diverse biomes and geographies. Most of the GFMs evaluated substantially underperform with respect to the supervised baseline, whether frozen or fine-tuned; only TESSERA and AlphaEarth Foundations (AEF) features with certain prediction heads yield a better performance. Our SOTA model trained on AEF embeddings and selected raw EO features achieves the best overall result (RMSE of 50.79 Mg/ha vs. 53.73 Mg/ha on AGBD features). This feature combination also generalizes better than AGBD features: it reduces RMSE by up to 25% in selected regions excluded from training, and, trained on single-year data alone, it outperforms the model trained on multi-year AGBD features. With labels from only 11 regions around the world, the AEF-based model also performs on par with the globally calibrated ESA CCI biomass product on independent reference plots. These findings indicate that GFM features can be highly informative for biomass estimation, when the underlying model is trained on rich multi-modal data.
Ghjulia Sialelli, Linus Scheibenreif, Jan Dirk Wegner +1
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
Remote sensing is increasingly relied upon to deliver actionable science for forest and wildfire risk management across large landscapes. Wall-to-wall, annually updated maps are a persistent need for effective forest management. Many planning systems and data collections combine disparate data sources with different purposes, vintages, and prediction quality, which leads to confounding behavior in operational planning systems. We introduce the VibrantForests framework, developed and applied to map forest attributes and provide a coherent foundation for effective forest and wildfire planning. VibrantForests includes a satellite-based forest structure model trained on lidar-derived samples and applied across the contiguous United States to concurrently generate estimates of canopy cover, canopy height, aboveground live tree biomass, basal area, and quadratic mean diameter at 10-meter resolution. We demonstrate predictive capability spanning the full spectrum of forest conditions ranging from sparse-canopy/low-biomass to dense-canopy/high-biomass. Results show that our model extends the range at which saturation is commonly encountered in comparable passive-sensor models, and reduces regression-to-mean behavior that commonly produces overestimation of forest attributes in small/sparse conditions and underestimation in large/dense conditions. The VibrantForests framework addresses a key limitation in large-area forest and wildfire planning by delivering coherent wall-to-wall estimates of management-relevant attributes at annual cadence and 10m resolution.
Luke J. Zachmann, David D. Diaz, Vincent A. Landau +8
Vibrant Planet Public Benefit Corporation, Truckee, CA, USA