We used soil samples collected from 325 locations worldwide. The sampling locations and sample collection protocol are described by Vasar et al. (2022). The sampling protocol was designed, following Taberlet et al. (2012), to minimise the influence of local heterogeneity by incorporating soil from multiple cores. The DNA mixture extracted from replicate soil samples is thus expected to represent collective local biodiversity as closely as possible. Briefly, at each sampling location, a site that was little disturbed by human activities was identified. About 20 g of topsoil (1–5 cm) was collected from up to the 40 points within an approximately 50 × 50 m sampling area. For further analysis, the samples were pooled per site, dried, and homogenised. For molecular analysis, a 2 g subsample of soil was collected from the pooled sample.

DNA was extracted from 2 g of dried soil using the PowerMax Soil DNA Isolation Kit. The short fragment of the trnL (UAA) intron (the P6 loop, 10-143 bp) was amplified using primers trnL-g (5′-GGGCAATCCTGAGCCAA-3′) and trnL-h (5′-CCATTGAGTCTCTGCACCTATC-3′) (Taberlet et al., 2007) for identification of plants. The P6 loop of the trnL intron is short enough to allow amplification of degraded plant DNA from soil samples (Taberlet et al., 2007). Thermal cycling included an initial hot-start denaturation at 95°C for 15 min, 35 cycles of denaturation for 30 s at 95°C, annealing for 30 s at 55°C, elongation for 10 s at 72°C and storage at 4°C. Barcode selection and PCR product preparation follow Vasar et al., 2022. Libraries were sequenced on the Illumina MiSeq platform, using a 2 × 150 bp paired‐end sequencing approach at Asper Biogene (Tartu, Estonia). Raw reads from this Targeted Locus Study have been deposited in the NCBI SRA (BioProject PRJNA659159).

Estimates of mean annual temperature (MAT), mean annual precipitation (MAP) and precipitation seasonality (SeaPrec) at sample locations were taken from the CHELSA database (Karger et al., 2017). These climatic variables are expected to be important drivers of plant communities (Sabatini et al., 2022) and have large-scale dynamics making them relevant to plots of variable size (GBIF, Kreft and Jetz, 2007 and Cai et al., 2023). Only variables that have annual measurements were selected from the CHELSA database to remove the possible effect of seasonality.

Plant trnL data were analysed using the gDAT pipeline (Vasar et al., 2021). Combining paired-end reads was conducted with FLASh (v1.2.1071, Magoč and Salzberg, 2011) using the default parameters (at least 10bp overlap with 75% identity) and enabling the option for outie alignments, which allows overlaps at the beginning of reads, rather than exclusively at the end. As the insert sizes (when outie alignment occurs) are shorter than MiSeq read length, sequences result in palindrome reads where adapter read-through occurs (Bolger et al., 2014). Adapters were removed from both ends using pipeline built-in functionality. Reads were demultiplexed, allowing 1 mismatch for the barcode and primer for both pairs. Only pairs where both reads had an average quality score of >30 were retained (after the removal of barcode sequences). Reads < 20bp (excluding primer length) were removed as they contain very little information for defining OTUs, assigning taxonomy and detecting potential chimeric sequences. Orphan reads (paired end reads that did not meet the conditions for combination and cleaning) were removed from the analyses, leaving 2,050,981 cleaned sequences. The VSEARCH (v2.14.166, Rognes et al., 2016) chimera filtering algorithm was used to remove putative chimeric reads in de novo mode, yielding 2,042,297 chimera-free sequences. Reads were clustered with VSEARCH at 97% identity into 32,127 OTUs (excluding singletons). Representative sequences (OTU centroids) for each non-singleton OTU were taxonomically classified using a BLASTn (v2.11.0+, Camacho et al., 2009) search against the INSDC (International Nucleotide Sequence Database Collaboration, status October 2022) database using up to 300 best hits for each query sequence, which were used to identify the lowest common ancestor by the pipeline based on the following criteria: > 80% alignment length; identification for species level > 97% identity, for genera > 95% and for family > 90%; consensus achieved when 51% of hits shared the same taxonomy at each rank. A total of 29,271 OTUs were identified as plants; 43 OTUs were identified as Bacteria, Fungi or Metazoa, which, along with OTUs that could not be identified, were removed from the data set prior to analysis.

The Global Biodiversity Information Facility (GBIF) collates a variety of information about organisms and their geographical locations. We took information about the distributions of vascular plant species in GBIF from Tamme et al. (2021) including 156,200,298 georeferenced vascular plant occurrences for 289,295 species (available for download: https://doi.org/10.15468/dl.4nqoev). The dataset was cleaned by the authors to group the species with synonym entries, reducing occurrences to 266,074 species, which were normalised into hexagonal cells using dggridR package in R (Bousquin, 2021), producing 65,612 equally sized grids cells (cell area of approximately 7,774 km2) worldwide.

We also took empirical estimates of plant alpha diversity worldwide from the information provided by Kreft and Jetz (2007), which presented geographic patterns of global vascular plant diversity at the species level. It used an exhaustive dataset of 1032 regions worldwide to generate expert-opinion-based continental to global maps of plant species richness. The dataset provides a cokriging estimate (cokrig), which was used as a richness estimate for our study. Finally, we considered the recently updated empirical estimates of plant alpha diversity worldwide produced by the Cai et al. (2022) (hereafter Cai) who used machine learning to improve the global models. We compared Cai with Kreft & Jetz dataset to assess the improvements.

To normalise eDNA sequence count data, we implemented the variance stabilising transformation (in R package DESeq2 v1.28.175), as suggested by McMurdie and Holmes (2014). The method uses fitted dispersion-mean relationships to normalise data with respect to sample size (sequencing depth of individual samples) and variance.

For each eDNA sampling location, the intersecting cell from the GBIF, Kreft & Jetz dand Cai et al., 2022 datasets was identified. Multiple data subsets were then generated with the following criteria: eDNA reads - raw read counts from eDNA sequences; eDNA c97% - raw reads clustered at 97% similarity threshold and singletons removed; eDNA VST - clustered reads at 97% similarity, variance stabilised transformed; eDNA family - INSDC BLASTn hits extracted at family level; eDNA species - INSDC BLASTn hits extracted at species level; Kreft & Jetz – Kreft & Jetz cokrig value; GBIF family - GBIF dataset grouped at family level; GBIF species - GBIF dataset grouped at species level; Cai – Cai et al., 2022 machine learning predicted richness using ensemble prediction. The data subsets were then correlated using Pearson’s correlation and plotted (using chart.Correlation from R package PerformanceAnalytics). The significance of coefficient estimates was estimated using Dutilleul et al.’s (1993) correction for spatial autocorrelation. The correction was implemented using a version of the modified.ttest function from the R package SpatialPack (Vallejos et al., 2020) that was adapted to incorporate great circle distances, i.e., distances over the surface of a sphere such as the Earth; the modified function is available from Davison et al. (2022).

Climatic drivers of eDNA and GBIF communities were identified using distance-based redundancy analysis (dbRDA; vegan package in R, Oksanen et al., 2013). Variation in sample distance matrices (Bray-Curtis distance following variance stabilising transformation) was modelled against all measured climatic variables. The significance of the effects was measured using permutation (n=999). Sampling locations were assigned to biogeographic realms following Olson et al. (2001).

We predicted the global distribution of the four different diversity measures (eDNA, GBIF, Kreft & Jetz and Cai) using generalised additive models (GAMs) and the spline-over-the-sphere algorithm in the R package mgcv, with the method ‘sos.smooth’ and the default arguments except k = 90 (Wood, 2003). This model can predict smooth variation in diversity values over the globe without producing edges. To make the spatial patterns comparable, we used for GBIF, Kreft & Jetz and Cai data only the subset of cells that intersected the eDNA sampling locations.

We estimated the availability of plant taxon reference sequences in the INSDC by searching “trnL[Gene Name] OR tRNA[Gene Name] OR Leu[Gene Name]”. We recorded 138,286 sequences, including 64,513 plant species and 529 plant families (accessed on October 2022). For each sampling location, we estimated the fraction of recorded GBIF families and species represented in INSDC; and the fractions of eDNA reads and OTUs getting a match against families and species in INSDC. We also calculated the overlap of family- and species-level taxa recorded in eDNA and GBIF. We then used GAM models, as described above, to assess geographical variation in the match between plant occurrence data (eDNA and GBIF) and reference sequence data.

Richness estimated using eDNA was weakly-to-moderately positively correlated with the corresponding estimates derived from GBIF (family: r = 0.20, p < 0.001; species: r = 0.11, p = 0.045), Kreft & Jetz (eDNA family: r = 0.30, p < 0.001; eDNA species: r = 0.25, p < 0.001) and Cai (eDNA family: r = 0.26, p < 0.001, eDNA species: r = 0.17, p < 0.01) at both the species and family level, with family-level correlations systematically stronger (Figure 1, Figure S1). The eDNA metrics were highly correlated with each other (r = 0.50-0.89, all p < 0.001), as were the GBIF metrics (r = 0.86, p < 0.001), and Cai and Kreft & Jetz (r = 0.74, p < 0.001). eDNA had a stronger correlation with Kreft & Jetz than with GBIF and Cai. Different taxonomic levels of GBIF were positively correlated with Kreft & Jetz (GBIF family: r = 0.55, p < 0.001; GBIF species level: r = 0.48, p < 0.001). Compared with Kreft & Jetz, Cai was more strongly correlated with GBIF family (r = 0.57, p < 0.001) but more weakly correlated with GBIF species (r = 0.44, p < 0.001). Pearson’s correlations with significance values corrected for spatial autocorrelation are shown in Table S1.

Generalised additive models identified similar global hotspots and coldspots of richness for all four datasets (Figure 2). dbRDA analysis revealed similar relationships between eDNA- or GBIF-derived plant family composition patterns and climatic factors, with colder and warmer regions generally clustering separately (Figure 3).

The annotated eDNA reads comprised 412 families and 4550 species, while 496 families and 106,490 species were recorded in the intersecting cells of the GBIF dataset. Altogether 304 families were shared between sets, while 108 and 192 were unique for eDNA and GBIF, respectively. In turn, 2663 species were shared between sets, while 1887 and 103,827 species were uniquely recorded in eDNA and GBIF, respectively.

The success of taxonomic assignment varied between different geographic regions. The highest rate of species-level assignment was 82.3%, with 82 out of 99 OTUs identified in a sample; the lowest rate was 5.6%, with only 2 OTUs out of 36 OTUs identified in a sample. In general the proportion of OTUs being assigned to species level was 38.7%. The proportion of reads and OTUs being successfully matched against species represented in INSDC was notably high in the peri-arctic region and lower elsewhere (Figure 4); the family-level pattern was less pronounced but still characterised by successful assignment at northern latitudes (Figure 4). Among GBIF occurrence records, there was similar geographic variation in the proportion of taxa (families and species) represented among trnL sequence data in INSDC (Figure S2). At the species level, representation in INSDC was highest at northern latitudes and decreased towards the equator, remaining low in the Southern Hemisphere. Among GBIF families, representation was highest in Europe, North America, and Australia. The proportion of eDNA reads and OTUs representing species and families that were also present in the corresponding GBIF data also exhibited substantial geographic variation: family-level matching was highest in Europe and North America with the strongest discrepancy in parts of Asia; species-level matching was highest at high northern latitudes (Figure S3).