Dominant shrub species in temperate deserts exhibit specialized adaptations and divergent evolutionary strategies in response to varying and extreme environmental stresses. However, it remains unclear how shrub species balance growth and carbon storage to cope with abiotic combined stresses across extensive spatiotemporal gradients. Guided by the Soil-Plant-Atmosphere Continuum (SPAC) theory and combined with a 20-year monitoring of non-structural carbohydrates (NSC), we conducted extensive field surveys across a representative temperate desert area. Using ensemble learning on 60 integrated environmental variables, the region was automatically classified into four SPAC systems that reflect gradients of combined temperature, precipitation, radiation, soil properties, and other factors. Results revealed divergent trade-offs between growth and carbon storage of shrubs mediated by intensity and combination of stresses. Shrubs in the Qinghai-Tibet Plateau faced severe temperature-water stresses, with growth limited by carbon storage. In contrast, shrubs in the Ningxia-Shanxi region tended to promote growth in minimal water stress. NSC mobilization and internal transport capacity were key determinants of shrub resilience to extreme climate events. These findings suggest that long-term evolutionary processes have shaped flexible carbon allocation strategies along environmental gradients. Therefore, understanding these adaptive strategies is crucial for predicting vegetation dynamics and ecosystem resilience under future climate scenarios.

Study area

The field sampling regions are mainly situated in the desert region of central Asia, spanning 37 degrees of longitude and 20 degrees of latitude, stretching from 78° E, 32° N to 115° E, 52° N. This expansive area encompasses several typical temperate deserts, arranged from west to east: Gurbantunggut desert, Taklamakan desert, Kumtag desert, Qaidam desert, Badain Jaran desert, Tengger desert, Ulan Buh desert, Kubuqi desert, Mu Us sandy land and Hunshandake sandy land. The orientation of the sample zone ran from northwest to southeast, essentially perpendicular to the precipitation contour in the northwest desert region of China. Within this region, the average annual precipitation exhibited considerable variation, ranging from 40 mm to 394 mm. The average annual temperature spanned from 4.32 ℃ to 9.79 ℃.

Crucially, the remote and largely uninhabited nature of these deserts makes them an unparalleled natural laboratory. Having grown quietly for millennia across these diverse landscapes with minimal direct human interference, the shrub communities here exhibit true self-adaptive strategies honed by natural selection. They are therefore ideal subjects for investigating the fundamental principles of plant physiological responses to environmental variation.

Field sampling

In July 2022, during the growing season of most desert shrubs, field surveys and sample collections were conducted in the designated study area. Sampling sites were chosen based on their ability to accurately represent the local plant diversity. Additionally, areas with stable topography were prioritized, while those influenced by human activities or environmental disturbances were excluded. At each site, three large plots measuring 30 m × 30 m were established, with a minimum distance of 1 km between them. Within each large plot, smaller subplots of 5 m × 5 m were positioned diagonally to focus on shrub sampling. In these subplots, only mature, healthy shrubs with uniform growth were selected and marked for study, with their length, width, and height recorded. For each shrub, 2–3 vigorously growing branches were pruned, specifically from the southern side and mid-height of the individual. Leaves (or assimilative branches, collectively referred to as leaf in this study) were removed from the branches, which were then cut into segments approximately 3 cm in length. The collected leaves were cleaned of surface dust, sorted, labeled, and immediately placed in a refrigerated container at 0–4°C with ice packs for transportation to the laboratory. Due to the limited number of leaves on some branches, the study ultimately established 94 sampling sites, 282 subplots, and collected a total of 1,356 leaf samples and 1,476 branch samples.

Measurement of non-structural carbohydrate

To halt all enzymatic activity as quickly as possible, all fresh samples were processed on the same day of collection (within 12 hours) by being microwaved at 800 W for 90 second.**** After transporting the samples to the laboratory, all samples were placed in a constant-temperature oven at 60°C and dried until a stable weight was achieved. this is an established and standard protocol for ensuring consistency in large-scale ecological studies. The dried samples were then ground into fine powder using a ball mill (5100 Mixer Mill, Metuchen, NJ, USA). Given that the sum of total soluble sugars and starch typically accounts for over 90% of total NSC content, this sum was used as an estimate of total NSC content. The NSC content was measured using the anthrone-sulfuric acid method. In brief, 0.05 g of powdered sample was extracted with an 80% ethanol solution for 12 hours, and the soluble sugar content was determined using the supernatant obtained after two rounds of centrifugation. The remaining precipitate was boiled in distilled water, hydrolyzed thoroughly with hydrochloric acid, and the resulting solution was centrifuged again to measure starch content. The absorbance of soluble sugars and starch was measured at a wavelength of 620 nm using the anthrone-sulfuric acid method on a multifunctional microplate reader (HR 7000; Hamilton, Reno, NE, USA). Their concentrations were quantified based on a standard curve and expressed in mg/g.

Data acquisition

The ERA5 global reanalysis dataset provides comprehensive meteorological data, offering optimal climate-related variables crucial for modeling NSC content. SPEI values were obtained from the Standardized Precipitation-Evapotranspiration Index database. Spectral characteristics within the study area were derived from various bands of Landsat 8 OLI and Landsat 7 ETM satellite imagery. Soil properties were sourced from the SoilGrids dataset, while parameters derived from the Digital Elevation Model (DEM) were calculated using data from the Shuttle Radar Topography Mission (SRTM). Vapor Pressure Deficit (VPD) data were obtained from the TerraClimate dataset, which provides high-resolution (∼4 km) monthly climate and water balance data from 1958 to the present, integrating multiple observational and reanalysis products to support ecological and hydrological research. The Global OCO-2 Solar-Induced Chlorophyll Fluorescence (GOSIF) dataset, which exhibits moderate to high performance in open shrublands, serves as a robust proxy for ecosystem-scale gross primary productivity and was used in this study to approximate shrub photosynthetic assimilation capacity. Finally, the high-resolution images that were used to calculate the field coverage of desert shrubs were derived from satellite data from the Chinese high resolution earth observation system. To ensure a unified analytical framework, all datasets were resampled to a common resolution of 1000 meters using the nearest neighbor method.

Amidst the swift advancement in sensor technology, scholars have developed various metrics rooted in the interplay among soil, vegetation, and spectral data, applying these metrics to assess soil and vegetation cover in specific regions. After a thorough literature review and integrating findings from previous studies, this study selected seven widely used spectral signature metrics as environmental variables to investigate NSC content in dominant desert shrubs: Normalized Difference Vegetation Index (NDVI), Soil adjusted vegetation index (SAVI), Carbonate index (CAI), Wetness index (WI), Brightness index (BI), Salinity index (SIA), Salinity index (SIB) [60]. The formula are as follows:

SAVI = (1 + L) * (NIR - R) / (NIR + R + L)

CAI = (SWIR1 - NIR) / (SWIR1 + NIR)

SI = (NIR * R) * 0.5

BI = sqrt(R^{2 + NIR}2)

SIA = G / NIR

SIB = (G - NIR) / (G + NIR)

where B, G, R, NIR, SWIR1, and SWIR2 correspond to reflectance in the blue, green, red, near infrared, shortwave infrared 1, and shortwave infrared 2 bands, respectively, derived from Landsat 7 ETM and Landsat 8 OLI satellite data.

Statistically analysis

To analyze NSC content in desert shrub communities, we normalized NSC levels across genera. First, we normalized data for each genus at different sampling sites, then calculated a weighted average at each site based on sample counts. The formula for community-level NSC content at each sampling point is:

NSC_community = Sum(NSC_i * Weight_i) / Sum(Weight_i)

where NSC_community is the community-level NSC content at the sampling site. Sum represents the summation from i=1 to n (where n is the total number of genera at the site). NSC_i is the standardized NSC content for genus i. N_i is the number of samples collected for genus i (which serves as the weighting factor).

The volume of an individual shrub was estimated based on the recorded dimensions of length, width, and height during sampling. The formula is as follows:

Volume = (pi / 6) * L * W * H

where Volume  is the estimated volume of the shrub, L, W and H represents the recorded length, width and height, respectivily.

SHapley Additive exPlanations (SHAP) values elucidate the decision-making process of a model, transforming it from an opaque "black box" into an interpretable system. Their primary advantage lies in detailing how each variable influences individual predictions, thereby enhancing the visualization and comprehensibility of the model's operations. For each prediction, the model assigns a SHAP value to each feature, quantifying its impact on the outcome. Compared to GAIN, SHAP provides a more consistent and theoretically grounded approach based on cooperative game theory, ensuring fair attribution of feature contributions across all possible feature combinations. This method not only increases model transparency but also enables precise explanations of how features drive predictions. The formula for the effect of SHAP value on the model is as follows:

y_hat_i = y_m + f(x_i1) + f(x_i2) + ... + f(x_ij)

where, y_m is the mean prediction and y_hat_i is the SHAP value for feature , x_ij with positive values increasing and negative values decreasing the prediction.

Moreover, we used the box-and-whisker plot method to identify the extreme values of 60 environmental variables over 20 years for SPAC system 1 through 4. Extreme values are defined as those that fall below the lower bound or above the upper bound of the interquartile range (IQR). The IQR method calculates the lower and upper bounds as follows:

Lower bound = Q1 - 1.5 * IQR

Upper bound = Q3 + 1.5 * IQR

where Q1  is the first quartile (25th percentile), Q3  is the third quartile (75th percentile), IQR is the interquartile range, defined as the difference between Q1 and Q3.

The network graphs effectively highlight the differences in environmental relationships across different areas and the varying associations among NSC content and environmental factors. We randomly sampled 2,000 pixels from each of the four sub-regions within the study area and extracted 60 environmental parameters, as well as LNSC, BNSC, and GOSIF values for these pixels. First, we conducted Pearson correlation analysis on these environmental, NSC, and SIF parameters using the "LinkET" package in R version 4.3.3. After calculating the correlation coefficients and significance values for all pairs of variables, only the significant results were retained. Next, we performed Principal Component Analysis (PCA) on each category of environmental variables using the "FactoMineR" package. We then conducted the same correlation analysis between the first principal component values and NSC content, as well as SIF values. Finally, we visualized the networks for different regions using Cytoscape software version 3.10.

Data conforming to a normal distribution were analyzed using one-way ANOVA and Tukey HSD method, while those not conforming were analyzed using the Dunn's test following a Kruskal–Wallis test. The dimensions and spatial resolution of environmental variables were aligned using the "terra" package. The "xgboost" package  was employed to develop XGBoost models, with grid search for the optimal hyperparameter combinations conducted via the "caret" package. SHAP values were calculated using the "shapviz" package. The MaxEnt model version 3.4.4 was used to predict the suitable zone of desert shrub NSC content sites. The U-net model was constructed using the “tensorflow” and “keras” libraries in Python 3.11, with image processing operations performed using the “opencv”, “pillow”, and “numpy” libraries. Spatial mapping and analysis were carried out using ArcGIS Pro version 3.3.

Changes after Feb 9, 2026: Renamed the "Shrub_data.xlsx" file and included additional data (Source_data_file.xlsx).