This study was conducted in the Eurasia steppe region of the Inner Mongolia Autonomous Region (IMAR) of northern China (Figure 2), which has a semi-arid and arid continental climate in the middle temperate zone. The study sites encompass three grassland types from east to west; meadow steppe (MS; 48.66°N, 119.13°E), typical steppe (TS; 44.14°N, 116.51°E), and desert steppe (DS; 41.77°N, 110.57–61°E). Over 1,125 km wide, the elevation ranges from 600 m in the east to 1,300 m in the west. Based on long-term meteorological observations (1980–2017), the mean annual temperature (MAT) ranged from −2.9°C to 6.5°C, with averages of −0.9°C, 3.1°C, and 4.8°C for MS, TS, and DS, respectively. The mean annual precipitation (MAP) ranged from 142.2 mm to 590.8 mm, with averages of 332.5 mm, 309.2 mm, and 279 mm for MS, TS, and DS, respectively. Most precipitation occurs from May to August. The soil is classified as chestnut and calcic brown (Bai et al., 2008). The MS is dominated by species Leymus chinensis, Carex duriuscula, and Stipa baicalensis. The TS is dominated by Stipa krylovii, Cleistogenes squarrosa, and Leymus chinensis. The DS is dominated by Stipa breviflora, Agropyron desertorum, and Cleistogenes songorica.

This study utilized multispectral satellite image data, species ground-truth survey data, and vertically projected RGB images of the quadrats. The multispectral satellite image data was obtained from the Sentinel-2A satellite, which was launched by the European Space Agency in 2015. Every 10 days, the satellite revisits the same viewing angles. A multispectral instrument (MSI) with 13 bands is mounted on the sentinel-2A satellite. Four bands are at 10 m spatial resolution, and six bands are at 20 m spatial resolution. In addition, there are three bands for atmospheric correction and cirrus cloud detection at 60 m spatial resolution. The reflectance bands at the finest resolution of 10 m were used in the study; blue (492.4 nm), green (559.8 nm), red (664.6 nm), and near-infrared (832.8 nm). Sentinel-2A Level 2A images were orthographic-corrected reflectance products at the bottom layer of the atmosphere and acquired from the European Space Agency Copernicus Open Access Hub.¹ The collected images were free of clouds and cirrus contamination over the study sites. The acquired multispectral images were dated 7 July for DS, 21 June for TS, and 11 July for MS, which aligns with the timing of the fieldwork in 2020. We extracted the pixel reflectance values from the four bands using the central coordination of quadrats in the R-package “raster” (Hijmans, 2022; R Core Team, 2022).

Grassland species were investigated during fieldwork in the mid-growing season in June and July 2020. Three 100 m*100 m plots were set up in each steppe type, oriented in a north–south direction, and quadrats (1 m*1 m) were demarcated along two diagonal lines within each plot at intervals of approximately 14 m (Figure 2). This interval of 14 m was selected to ensure that each quadrat aligns to one 10 m-pixel from the multispectral image. In total, there were 37 quadrats in MS, 49 quadrats in TS, and 47 quadrats in DS. Each quadrat was situated around a central point identified by Real Time Kinematic- Global Navigation Satellite System (RTK-GNSS) to form a north–south oriented square with an accuracy of 0.01 m. After the demarcation, we recorded the name and individual count of each species falling into each quadrat. The height of individual species was measured in succession three times. The aboveground biomass of each species was weighed after harvesting and drying.

Significant differences were observed in the four metrics of α-diversity along the precipitation gradient, which ranged from 332.5 mm for MS to 309.2 mm for TS to 279 mm for DS (p < 0.05; Figure 4). V-shaped patterns were observed in the three Hill numbers (i.e., species richness, expShannon, and invSimpson indices). There were large differences in species richness between steppes, and the TS was the least diverse (Figure 4A). When expShannon index was used, there was an insignificant difference between MS and DS (p > 0.05; Figure 4B), but when invSimpson index was used, the difference between MS and TS was insignificant (p > 0.05; Figure 4C). With the three Hill numbers, which evaluated rarity to a lesser degree than the Pielou index, MS, which had many rare species, became less diverse compared to TS and DS. DS had higher species richness in the expShannon and invSimpson indices than TS. The Pielou index showed the opposite pattern with the peak in TS (Figure 4D). Overall, the α-diversity metrics followed divergent patterns of variation over the region.

A negative relationship was observed between species richness and Pielou index (r = −0.37, p < 0.001; Figure 5A) over the three steppes. In particular, there was a strong negative relationship between species richness and the Pielou index in TS (r = −0.54, p < 0.001; Figure 5C). In contrast, species richness and the Pielou index had insignificant relationships in MS (r = 0.09, p = 0.585; Figure 5B) and DS (r = 0.09, p = 0.526; Figure 5D). The aggregate data representing the full region showed that the expShannon index correlated most to species richness (r = 0.56, p < 0.05; Figure 5A), and the invSimpson index correlated most to the Pielou index (r = 0.68, p < 0.05; Figure 5A). Specifically, the expShannon and invSimpson indices correlated to the Pielou index more strongly than to species richness in MS (r = 0.94 and 0.94, p < 0.05; Figure 5B) and DS (r = 0.76 and 0.85, p < 0.05; Figure 5D). In contrast, the expShannon index correlated to species richness more than to the Pielou index in TS (r = 0.61, p < 0.05; Figure 5C). This implies that the relationship of α-diversity indices at a local scale could strongly influence the overall relationship across the region.

Correlations between α-diversity indices and the many variables varied in their strength, with r ranging from 0.01 to 0.69 (Figure 6). H_mean, H_std, SA_rN, and PR significantly explained the variation in species richness, with r-values of −0.69, −0.63, −0.64, and 0.61, respectively. H_mean, H_std, PR, NDVI, and EVI significantly correlated to the expShannon index, with r-values of −0.45, −0.45, 0.49, −0.46, and − 0.46, respectively. H_mean, H_std, PR, and EVI significantly explained variation in the invSimpson index, with r-values of −0.24, −0.26, 0.35, and − 0.30, respectively. The Pielou index was significantly explained by H_mean, DIVI, LPI, and SIDI, with r-values of 0.32, 0.44, −0.41, and 0.41, respectively. Vertical structure and SVH variables explained species richness variation among α-diversity metrics well, while community pattern variables characterized the Pielou index better.

Responses of vertical structure and community pattern variables to species richness and the Pielou index were asynchronous. H_mean and H_std correlated strongly and negatively to species richness (r = −0.69 and  − 0.63, p < 0.001), but they correlated slightly positively to the Pielou index (r = 0.32 and 0.25, p < 0.01). Meanwhile, PR and LPI correlated moderately and positively to species richness (r = 0.61 and 0.35, p < 0.001), but they correlated negatively to the Pielou index, with r-values of −0.11 and − 0.41, respectively. The scenario also appeared in SVH variables. Overall, SA_rN, NDVI, and EVI correlated strongly and negatively to species richness (−0.64 < r < −0.46, p < 0.001) but weakly and positively to the Pielou index (0.06 < r <  0.20, 0.02 < p < 0.50).

The slopes of the two-variable regressions also indicated the asymmetrical responses of explanatory variables (i.e., vertical structure, community pattern, and SVH) to variation in α-diversity (Figure 7). The fitted surfaces of species richness were negatively related to those of the Pielou index (Figures 7A,D). In contrast, slopes of the fitted surfaces of the expShannon and invSimpson indices are gentle compared to the sharp and opposite slopes for species richness and evenness (Figures 7B,C). Variation in species richness (0.44 < R² < 0.51; Figure 7A) and the expShannon index (0.20 < R² < 0.24; Figure 7B) were better explained than the invSimpson (0.05 < R² < 0.15; Figure 7C) and Pielou indices (0.088 < R² < 0.14; Figure 7D). H_mean negatively related to species richness (slope = −0.69, p < 0.001) and positively affected the Pielou index (slope = 0.53, p < 0.01). H_std was more influential on the Pielou index (slope = −0.23, p = 0.25). H_mean and H_std were important and complementary in describing α-diversity. For the community pattern variables, SHDI was better at describing the Pielou index (slope = 0.38, p < 0.001). In addition, NDVI outperformed other SVH variables in characterizing the variation in Pielou index (slope = 0.76, p < 0.001) and Hill numbers, with the slope of −1.59, −0.71, and − 0.18 for species richness (p < 0.001), expShannon index (p < 0.001) and invSimpson index (p = 0.37), respectively.

By integrating three categorized variables using SEM, the explanations were improved by 5, 1, and 14% for species richness, the expShannon, and Pielou indices, respectively (Figure 8). As with the negative relationship between richness and Pielou index (Figure 5), the SHDI variable showed a strongly opposite response to species richness compared to the Pielou index (r_p = −0.20 versus 0.33), but a diminished relationship was observed in both the expShannon and invSimpson indices (r_p = 0.14 and 0.21) in the category of community pattern. As for the vertical structure category, H_mean and H_std achieved the highest path coefficient with species richness and evenness (r_p = −0.97 and 0.72). Meanwhile, similar situations with asymmetrical and disproportional relationships among these four metrics were observed (−0.97 < r_p < 0.72). In contrast, NDVI in the SVH category showed moderate performance in explaining variation in the four metrics of α-diversity. As for the interplay of three categorized variables per se, we found that NDVI was strongly dependent on vertical structure (r_p = 0.87). Variation in R², which ranged from 0.13 to 0.56, indicated that the responses of the explanatory variables were inconsistent and mirrored the complexity of the α-diversity description. Overall, the variation in species richness and Pielou index (R² = 0.56 and 0.28) were modeled better by SEM than the expShannon and invSimpson indices were (R² = 0.25 and 0.13).

Precipitation has been demonstrated to be a dominant environmental factor influencing the species diversity of grassland ecosystems by promoting species richness, heterogeneity, evenness, functional divergence, and biomass (Kang et al., 2007; Reynaert et al., 2021). However, species are idiosyncratic and highly depend on their local plant community in their response to precipitation variation across various heterogeneous patches and niches (Tsafack et al., 2021). This exposes the intrinsic complication in predicting how grassland communities respond to climatic change (Barnett and Facey, 2016). Up from the prior meteorological observations of 2010–2019, the precipitation in 2020 increased by 30.96% in TS and 22.18% in DS (Figure 9A). As a result, the proportion of annual species increased by 27.27 and 28.57% (Figure 9B), respectively. This strong and disproportional response of annual species to precipitation fluctuation between TS and DS bloated species richness and the expShannon index to abnormally high levels. This scenario was especially promising in DS, resulting in higher species richness, expShannon, and invSimpson in DS compared to TS. The same process also occurred in the abundance of dominant species for TS and DS (Figure 9B). The concurrent increases in annual species and abundance of dominant species led to evenness (i.e., Pielou index) declining, emerging the negative relationship between species richness and evenness in TS. The differing sensitivity of species diversity to precipitation fluctuation in DS and TS and the adverse variations in species richness and evenness highlight the complexity and variability of α-diversity in response to climate change (Zhang et al., 2014).

Consistently, MacDonald et al. (2017) also noted asymmetric variation in α-diversity in response to annual precipitation and the growing season temperature, and this resulted in an increased abundance of rare species as well as a disproportionally increased abundance of common species. These cases demonstrate the prevalence of the effect of environmental change on the pattern of α-diversity (Kang et al., 2007; Wilsey and Stirling, 2007; Smith et al., 2022). In contrast to precipitation increases in TS and DS, annual species were constrained in MS due to declining precipitation (4.80% decrease), which aligned with the significantly reduced aboveground biomass (AGB) (>Supplementary Table S1). These above findings speak to the first objective, where the differing sensitivity of α-diversity metrics and interactive behavior between richness and evenness reflect intricate and asynchronous responses of α-diversity to precipitation variation across various heterogeneous patches and niches. In addition, herbivores tend to prefer perennial and tall plants to annual and short plants (Díaz et al., 2007), which accordingly leads to a decrease in the abundance of perennial forbs but is insignificant for the perennial bunchgrasses (Liang et al., 2018). When precipitation is in short supply, increased consumption of net primary production by grazing herbivores aggravates the community more strongly (Liang et al., 2018). However, due to the lack of investigation on grazing practices over the region, we are limited in our ability to examine the integrated effects of grazing on variation in species diversity.

The community pattern can be conceptualized as a horizontal space that individual plant emergence, aggregation, and the physical arrangement relevant to resource facilitation and/or competition (Wesuls et al., 2013). Based on the empirical regression and SEM analysis (Figures 7, 8), we found that the community pattern significantly correlated with species richness, invSimpson, and Pielou indices. In contrast, vertical structure partly reflected the variation in species diversity because of the height heterogeneity existing among species (Torresani et al., 2020). The underperformance of vertical structure can be inferred to the grazing impact of large herbivores, which leads to a homogeneous plant height and then constrain the effort on differentiating α-diversity variation (Zhu et al., 2012). SVH has been widely used to detect plant biodiversity in a broad range of ecosystems, primarily by linking spectral properties to species and habitat distribution (Wang and Gamon, 2019). Currently, spectral variation is utilized to identify the plant functional group related to plant light absorption, and scattering is used to detect plant spectral properties (Cavender-Bares et al., 2017). In addition, the grassland community composition (e.g., species richness, evenness) can be sensitive to spectral variation, indicating that species richness and evenness together yield a more robust SVH relative to either diversity metric alone (Wang et al., 2018). Here, we found that SVH variables only explained 21 and 7% of the variation in the expShannon and invSimpson indices (Figures 7B,C). This weak explanation was attributed to the strong collinearity with properties of community pattern and vertical structure, which led to the insensitivity of spectral recognition. Accordingly, we questioned the value of exclusively using SVH to model α-diversity metrics. Additional factors for determining SVH are considerably constrained by atmospheric influence (Myneni and Asrar, 1994), limitations of spatial and spectral resolution (Bannari et al., 2002), the influence of topography (Zarnetske et al., 2019), soil background and the nature of optical saturation (Aneece et al., 2017). Altogether, the many limitations of SVH signify the caution of using SVH to explain α-diversity variation. This speaks to the last two objectives, where overall asynchronous variations in α-diversity metrics due to precipitation variation result in the asymmetric performance of the three categorized explanatory variables, thus demonstrating the limitations of using SVH to model variation in α-diversity metrics.

The common mission to depict the number of species and their relative abundances across different scales and biomes has led to the creation of multiple metrics of species diversity (Pielou, 1966; Hill, 1973). Among α-diversity metrics, the Hill numbers were meant to characterize α-diversity by focusing on the abundance of rare and common species (Chao et al., 2014; Roswell et al., 2021). A robust explanation of the expShannon and invSimpson indices is shown in Table 3, as a linear function of species richness and evenness, with uniform positive regression coefficients achieved over the region. MacDonald et al. (2017) suggested that coupling species richness with evenness can broaden the knowledge of the trends of fluctuating species diversity. However, we emphasize that the negative correlation between species richness and evenness led to unpredictable and inconsistent interpretations of α-diversity variation in the resulting models. Furthermore, richness in TS strongly affected the variation in the expShannon and invSimpson indices, while evenness in MS and DS contributed more to the variation in the expShannon and invSimpson indices (Table 3). This means that drawing out the distinctions between α-diversity metrics can be paramount for improving the performance of the integrated SEM scheme for depicting α-diversity variation, regardless of various heterogeneous patches and niches.

Climate change and anthropogenic activity have led to precipitous biodiversity loss and rampant variation in biodiversity distribution and composition (Ryabov et al., 2022), affecting the underlying mechanisms of biodiversity variation (Wilsey et al., 2005). Local species loss averages up to 14% (Newbold et al., 2015), while globally, over 20% of vascular plant species face severe threat of extinction (Schweiger et al., 2018). These losses have raised concerns about the consequential changes in ecosystem functioning and services, including nutrient cycling, climate regulation, and food production (Qiu and Cardinale, 2020). The accelerating declines in wild species and the complexity of the involved dynamics have sparked a concerted effort to answer the question about the distribution of biodiversity and the identification of high-priority areas for conservation to maintain essential ecosystem functioning and services (Wang and Gamon, 2019).

Airborne remote sensing techniques are widely utilized by scientific communities to model biodiversity variation (Coops et al., 2021). A solid method is needed for monitoring biodiversity variation over temporal and spatial scales. However, the existing theories that have been proposed based on spectral variation are still limited in addressing these issues. Ultra-resolution data portrayed by high-resolution airborne satellites (e.g., worldview 3) or unmanned aerial vehicles (UAVs) offer more information, reflecting subtle properties in plants (Wang and Gamon, 2019). State-of-the-art machine learning algorithms, such as support vector machine classification (SVM) and decision tree models, have been widely used to explore this data, and desirable results have been achieved (Delalieux et al., 2012; Lopatin et al., 2017). Meanwhile, high-density point clouds, scanned using the emerging technology of Light Detection and Ranging (LiDAR), are used for phenotypic property inversion through three-dimensional reconstructions (Sankey et al., 2021). Successful applications have emerged in forest ecosystems with species diversity detection by height heterogeneity (Marselis et al., 2020; Torresani et al., 2020), but applications in grazing ecosystems are still nascent. To understand the responding and interactive mechanism of these biodiversity metrics to climate change, it will be necessary to expand investigations into a heterogeneous grassland ecosystem (Gholizadeh et al., 2022). Further, grazing practices imposed by humans are another factor that should be assessed. The theoretical SEM proposed in this study is expected to combine emerging methods and techniques in a way that sheds light on biodiversity variation in various grassland ecosystems, coexisting with the above two forcing scenarios.