The Qinghai-Tibet Plateau (QTP) is ecologically significant due to its unique biodiversity and vulnerability to climate change and rapid urbanization. Among its emerging urban habitats, wastelands with relatively low anthropogenic disturbances offer key refuges and stepping-stones for spontaneous plants. However, the patterns and drivers of spontaneous plant diversity in these habitats remain poorly understood. To address this, we surveyed spontaneous plant communities across 17 cities on the QTP. We found that native species dominated urban wasteland flora (85.8%), whereas invasive species constituted over half (65.8%) of the non-native species. The results showed that land-use legacy effects and environmental filtering to shape community assembly. GLMM analyses further reveal that the climatic background constitutes a key factor shaping community diversity variations across the QTP, exerting the most significant influence on species richness, particularly through precipitation and wind speed. Urbanization and habitat quality factors jointly shape the diversity structural characteristics of communities. However, non-native and invasive species exhibit heightened sensitivity to local habitat quality. Furthermore, differences between communities across all groups were due to species turnover, though the driving factors differed between groups. Native species exhibited stronger overall ecological adaptability, whereas differences between communities of non-native and invasive species were primarily driven by human disturbance and habitat conditions. These findings underscore the ecological value of urban wastelands on the QTP for biodiversity conservation and invasion management under rapid urbanization and climate change.

Environmental variables

So, environmental variables, including climatic factors, urbanization metrics, and habitat quality, were collected to assess the drivers of spontaneous plant species diversity in urban wastelands. We consider two primary factors from classic climate characteristics: temperature and moisture. Additionally, we account for the unique conditions of high altitudes. Therefore, our climatic indices include: Mean annual precipitation (MAP), mean annual temperature (MAT), total annual solar radiation (TAR), and mean annual wind speed (MAW). Urbanization intensity was represented by seven metrics, including gross domestic product (GDP), GDP per capital (RGDP), and population density (PD), human activity intensity (HAI), urban population (UP), impervious surface ratio within a 500 m buffer surrounding each wasteland (Sealed.wl, it reflects the land-use patterns surrounding the wasteland and the impact of edge effects, thereby providing an indication of the availability of gene pools for species dispersal and escape from the site), and night light intensity (NLI). As for habitat quality, we used seven indicators derived from field surveys and spatial analysis, including factors that directly influence the plant diversity of wasteland: Total woody plants cover in each wasteland (TWC) is defined as the sum of the cover of all woody plant species within a given wasteland (non-mutually exclusive canopy overlap), calculated as the proportion of the area occupied by woody plants (primarily including shrubs, trees, and woody vines) relative to the total wasteland area, which reflects the overall dominance of woody vegetation and, to a certain extent, the baseline natural conditions of the area. It also reflects vegetation structure and habitat complexity, associated with ecological stability and biodiversity; Total bareness in each wasteland (TB), which serves as a proxy for disturbance intensity and habitat degradation. By analyzing land-use changes at each abandoned site using historical versions of Google Maps (1984-2025), it is possible to estimate the time of abandonment (AT). The study of wasteland covered a continuous gradient ranging from the early stages of abandonment (less than 2 years) to the late stages of abandonment (exceeding 20 years). This indicator can serve as an indicator of the successional stage and ecological recovery potential; Wasteland area (W.area) reflects the size of the habitat island; Soil organic carbon content in the 0 – 30 cm layer (Ct) reflects the soil quality of wasteland. It also includes factors that indirectly affect biodiversity: Perimeter-area ratio of wasteland (PA, PA = Wasteland perimeter / Wasteland area, the area and perimeter of the wasteland were calculated using ArcGIS 10.7) and landscape shape index (LSI), which indicates shape regularity and geometric complexity of wastelands.

Data analysis

Species richness directly quantifies the pool of native versus invasive taxa, addressing our core question of how climate, urbanization, and habitat structure shape wasteland communities. The Shannon–Wiener index weights rare species more heavily, capturing early-stage invasions and the persistence of uncommon native specialists under environmental stress. The Simpson index emphasizes dominant species and reveals whether communities are controlled by a few aggressive invaders or maintain balanced competitive hierarchies. Pielou's evenness isolates equitability independent of richness, distinguishing functionally diverse assemblages from those dominated by few species—a key contrast when evaluating invasion impact and native resilience. Together, these four complementary metrics capture the full spectrum of community structure—from species accumulation to compositional balance and dominance—ensuring we do not rely on a single, potentially misleading summary. The following four diversity indices were selected as community diversity indices.

Before model fitting, predictors were standardized by using the scale function in the base package to mitigate scale disparities. Then, a generalized linear mixed model (GLMM) was performed by using independent variables (climatic factors, urbanization metrics, and habitat quality) and nested random intercepts (CityID/WastelandID) with each dependent variable (four plant diversity indexes of total, native, non-native, and invasive species groups). Since richness is a count variable, the Poisson distribution was used by glmer function of lme4 package; the Gaussian distribution was used for other continuous data by using lmer function of lme4 package, and due to zero inflation in the non-native and invasive groups, the zero-inflation model was selected by using glmmTMB package with glmmTMB function (Bates et al., 2015; Brooks et al., 2017). After constructing the model, the variance inflation factor (VIF) of the model was further examined by using vif function in the car package, and variables with VIF greater than 5 were removed from the model one by one until all variables had VIF values less than 5 (James et al., 2014; Scott, 1997; Vittinghoff et al., 2005).

We employed full subset model selection using the dredge function from the MuMIn R package, and identified the best model based on the corrected Akaike Information Criterion (AICc). When multiple models fell within ΔAICc < 2 of the best model, model averaging was performed using the model.avg function to account for model uncertainty (Anderson and Burnham, 2004). To quantify the individual explanatory power of each fixed effect, we decomposed the marginal R² ( ) using the r.squaredGLMM function in the MuMIn R package, allowing us to estimate the relative contribution of each predictor to the overall model performance (Gao et al., 2023; García-Palacios et al., 2018). All the above calculations were performed in R 4.4.2 software (R Core Team, 2023).

To investigate the driving mechanisms behind community composition differences, we first decomposed overall β-diversity into species turnover (βsim) and nesting components (βsne) using the Sørensen β-diversity index (βsor) (Baselga, 2010). β-diversity calculations were performed using the R package betapart, with the beta.pair function generating three index matrices—βsor, βsim, and βsne—for species of different origins. Subsequently, to assess the explanatory power of environmental variables on spatial variation in β-diversity, we employed Redundancy Analysis (RDA) based on distance matrices (Legendre and Gallagher, 2001). RDA analysis was implemented using the vegan package. The rda function calculated constrained axes, and model significance was assessed via permutation tests (n = 999) using anova.cca. This further quantified the independent contributions of each environmental factor to community variation.