Data: Decadal-scale observations are key to detecting the stabilizing effects of plant diversity in natural ecosystems

Regional Data Our regional transect data includes 160 sampling sites across alpine grasslands on the Tibetan Plateau, China (Fig. 1a). Sampling occurred during peak growing seasons (July-August) in 2019, 2020, and 2022. Sites span latitudes 28°35’24” to 33°12’9" N and longitudes 79°48’43" to 96°44’24" E, with a minimum 15 km distance between neighboring sites. Elevations range from 3614 m to 5268 m, averaging around 4500 m. We used a 20-year dataset (2000-2020) from the China Meteorological Data Service Centre (CMDC; data.cma.cn) for mean annual precipitation (MAP) and mean annual temperature (MAT) at the sampling sites. MAP values ranged from 49 mm to 646 mm, and MAT values from -5.7 to 7.2℃. Evapotranspiration (ET) data, obtained from the Penman-Monteith-Leuning Version 2 dataset (CAS IGSNRR PML V2), showed mean annual ET values from 250 mm to 625 mm. We calculated the aridity index by dividing mean annual ET by MAP. We selected four key plant functional traits: plant height (HT), specific leaf area (SLA), leaf water content (LW), and leaf chlorophyll content (SPAD value). Plant height was measured in 10 quadrats per site per species. SLA (cm²/g) is the leaf area to dry mass ratio, and LW (%) is the water content to dry weight ratio. Leaf chlorophyll content (SPAD value) was measured using a SPAD-502 chlorophyll meter on ten healthy leaves per species per site. Functional diversity was quantified using functional richness (FRic), functional evenness (FEve), and functional divergence (FDiv), along with community-weighted mean (CWM) trait values. Species diversity was calculated using the Shannon-Wiener diversity index (H' = −Σ Pi ln Pi), where Pi is the proportion of each species within a quadrat. Faith’s phylogenetic diversity (PD) was calculated using the V.PhyloMaker2 package in R, based on a PhytoPhylo megaphylogeny. EVI data for each site, sourced from the MOD13Q1 product on Google Earth Engine (earthengine.google.com), has a 250 m × 250 m spatial resolution and a 16-day temporal resolution. The maximum EVI during the growing season (June-August) was used as a proxy for primary productivity. Temporal stability of primary productivity was assessed by calculating the ratio of mean annual maximum EVI value (μ) to the standard deviation (σ) of annual maximum EVI values (2000-2022). Coefficient of Variation (CV) of climatic variables， calculated as the interannual variability of precipitation and temperature over the same time periods used for assessing productivity stability. Global Data We utilized the most comprehensive and accurate predictive maps of taxonomic (species richness), phylogenetic (Faith’s PD), and functional diversity in terrestrial angiosperms. This database (https://doi.org/10.6084/m9.figshare.21583083) compiles distribution data for 332,293 angiosperm species, with species richness determined by counting species occurrences per geographic unit.To match the GPP stability resolution, plant diversity data were resampled from 25 × 25 km to 0.5° × 0.5°resolution. We also used the FLUXCOM GPP dataset (1958-2017), derived from FLUXCOM observations and upscaled through machine learning. The dataset, with a monthly temporal resolution and a spatial resolution of 0.5×0.5 degrees (approx. 55.5 km × 55.5 km). We computed the ratio of mean annual maximum GPP value (μ) to the standard deviation (σ) of annual maximum GPP values (1951-2017) to represent GPP temporal stability over 5 to 67 years. Abbreviations: AI: Aridity index FRic: Functional richness FEve: Functional evenness FDiv: Functional divergence CWM.SPAD: Community-weighted mean SPAD values CWM.SLA: Community-weighted mean specific leaf area CWM.Height: Community-weighted mean plant height CWM.Leafwater: Community-weighted mean leaf water content Shannon: Shannon-Wiener diversity index TD: Taxnomic diversity FD: Functional diversity PD: Phylogenetic diversity Stability: Productivity stability over years (5-20 or 67)