Data and code from: Shrub species identity modulates plant-arthropod diversity coupling across a dryland aridity gradient
Data files
Sep 18, 2026 version files 98.72 KB
-
Dryad_resubmission.zip
81.95 KB
-
README.md
16.77 KB
Abstract
Shrubs are key ecosystem engineers in drylands, yet their effects on biodiversity range from strongly positive to negative depending on species identity and climatic stress. A critical knowledge gap remains in understanding how shrub species identity modulates plant-arthropod biodiversity and their coupling across aridity gradients.
To address this gap, we examined the diversity of herbaceous plants and ground-dwelling arthropods, along with soil properties, in microhabitats formed by native Artemisia ordosica, planted Caragana korshinskii, and adjacent open spaces. This study spanned three consecutive years along a pronounced aridity gradient (MAP 400 – 184 mm) in northern China.
Our findings revealed that C. korshinskii shrubs exerted a neutral effect on herbaceous diversity even under high aridity, whereas A. ordosica shrubs inhibited herbaceous diversity, with suppression intensifying as aridity increased. Contrastingly, neither shrub species affected arthropod diversity across the full aridity gradient. Furthermore, a consistently close and strong plant-arthropod diversity coupling was observed across the three microhabitats. Along the aridity gradient, the C. korshinskii shrubs maintained a stable and high coupling, whereas the A. ordosica shrubs maintained a low coupling, comparable to open spaces. Where, the coupling reduced along aridity, exhibiting a threshold (~ 270 mm) below which the coupling weakened.
Synthesis. Together, these results suggest that shrub species identity is a decisive regulator of plant and arthropod diversity and their coupling across aridity gradients. C. korshinskii maintains a stable and high-coupling strategy, while A. ordosica exerts a low-coupling strategy along aridification. These findings provide a mechanistic foundation for incorporating shrub identity into climate-adaptive restoration strategies in dryland ecosystems.
Plant and Arthropod Diversity, Coupling Index, and Structural Equation Modeling in a Dryland Shrubland System
Dataset DOI: 10.5061/dryad.2280gb67x
Overview
This dataset contains plant species richness, arthropod taxa richness, soil properties, and a plant-arthropod coupling index collected across three habitat types (open grassland, Artemisia shrub, and Caragana shrub) at three sites along a precipitation gradient in the drylands of northern China. Data were collected over three years (2019-2021).
The accompanying R scripts reproduce all statistical analyses and figures presented in the associated manuscript, including generalized linear mixed-effects models (GLMMs), zero-inflated beta regression, and piecewise structural equation modeling (pSEM).
Files in This Submission
Archive
Dryad_resubmission.zip
- Compressed archive containing all data files and R scripts listed below
- Extract this archive to access the individual files; all scripts assume the data files are in the same working directory
Data Files
Note: Variables in tabular data files are described in the section "Variable Descriptions" below.
HRrichness.xlsx
- Plant species richness data at the plot level
- 144 rows x 6 columns
- Used by scripts:
1code-plant_diversity_among_habitat.R,2code-plant_diversity_along_aridity.R
FRrichness.csv
- Arthropod taxa richness data at the plot level
- 144 rows x 6 columns
- Used by scripts:
3code-arthropod_diversity_among_habitat.R,4code-arthropod_diversity_along_aridity.R
coupling2.csv
- Plant-arthropod coupling index values at the plot level
- 144 rows x 6 columns
- Derived from
master_standardized_data.xlsx - Used by scripts:
5code-coupling_among_habitat.R,6code-coupling_along_aridity.R
sem.csv
- Full dataset for piecewise structural equation modeling (pSEM), including standardized richness values, coupling index, shrub type, and soil physicochemical properties
- 144 rows x 19 columns
- Derived from
master_standardized_data.xlsx - Used by script:
7code-pSEM.R
master_standardized_data.xlsx
- Master dataset containing all variables, including the standardized richness values (HR-standardized, FR-standardized) and additional columns (treatment, microhabitat) not used directly by the analysis scripts
- 144 rows x 21 columns
- Source file from which
coupling2.csvandsem.csvwere derived; included for transparency and data provenance
Analysis Scripts
1code-plant_diversity_among_habitat.R
- Negative binomial GLMM comparing plant species richness among three habitat types (OP, AT, CA)
- Produces diagnostic plots, boxplots, predicted bar plots with significance letters, trend plots, and result tables
2code-plant_diversity_along_aridity.R
- Poisson GLMM with habitat x precipitation interaction for plant species richness along the precipitation gradient
- Produces interaction plots, simple slope tests per habitat, zero-inflation diagnostics, and result tables
3code-arthropod_diversity_among_habitat.R
- Poisson GLMM with bootstrap cluster-robust standard errors comparing arthropod taxa richness among three habitat types
- Produces diagnostic plots, boxplots, predicted bar plots, and result tables
4code-arthropod_diversity_along_aridity.R
- Poisson GLMM with habitat x precipitation interaction for arthropod taxa richness along the precipitation gradient
- Produces interaction plots, simple slope tests per habitat, and result tables
5code-coupling_among_habitat.R
- Zero-inflated beta GLMM comparing the plant-arthropod coupling index among three habitat types
- Produces boxplot with jittered points colored by coupling type (zero vs. positive), predicted marginal means with significance letters
6code-coupling_along_aridity.R
- Zero-inflated beta GLMM with habitat x precipitation interaction for the coupling index along the precipitation gradient
- Produces interaction plots, simple slope tests per habitat, and result tables
7code-pSEM.R
- Piecewise structural equation model (pSEM) examining direct and indirect pathways from precipitation and shrub type through soil properties and plant/arthropod richness to the coupling index
- Runs three separate models: overall (OP vs. all shrubs), OP vs. Caragana, and OP vs. Artemisia
- Outputs a coefficient table (
sem_coefficients.csv)
Variable Descriptions
HRrichness.xlsx and FRrichness.csv
Both files share the same column structure. HRrichness.xlsx contains plant data; FRrichness.csv contains arthropod data.
| Column Name | Description | Unit / Values |
|---|---|---|
year |
Year of sampling | 2019, 2020, 2021 |
location |
Study site code (three sites along precipitation gradient) | YL (wet, 400 mm), YC (medium, 292 mm), ST (dry, 186 mm) |
precipitation |
Mean annual precipitation at the site | mm (400, 292, 186) |
habitat |
Habitat type | OP = open grassland; AT = Artemisia shrub; CA = Caragana shrub |
replicate |
Plot replicate number within each habitat-site-year combination | Integer (1-6 for OP; 1-5 for AT and CA) |
richness |
Species (plant) or taxa (arthropod) richness observed in the plot | Count (integer) |
coupling2.csv
| Column Name | Description | Unit / Values |
|---|---|---|
year |
Year of sampling | 2019, 2020, 2021 |
location |
Study site code | YL, YC, ST |
precipitation |
Mean annual precipitation at the site | mm |
habitat |
Habitat type | OP, AT, CA |
replicate |
Plot replicate number | Integer |
coupling |
Plant-arthropod coupling index, calculated as the ratio of standardized plant richness to standardized arthropod richness (or vice versa), bounded between 0 and 1; a value of 1 indicates perfect coupling, 0 indicates complete decoupling | Continuous [0, 1] |
sem.csv
| Column Name | Description | Unit / Values |
|---|---|---|
year |
Year of sampling | 2019, 2020, 2021 |
location |
Study site code | YL, YC, ST |
precipitation |
Mean annual precipitation at the site | mm |
habitat |
Habitat type | OP, AT, CA |
replicate |
Plot replicate number | Integer |
shrub |
Shrub type indicator (contrasts for pSEM) | 0 = OP (open grassland); 1 = Caragana shrub; -1 = Artemisia shrub |
plant.richness_standardized |
Plant species richness standardized by habitat and site (HR-standardized in master file) | Continuous |
arthropod.richness_standardized |
Arthropod taxa richness standardized by habitat and site (FR-standardized in master file) | Continuous |
coupling |
Plant-arthropod coupling index | Continuous [0, 1] |
pH |
Soil pH | pH units |
EC |
Soil electrical conductivity | uS/cm |
ST |
Soil temperature | degrees C |
SM |
Soil moisture | % (volumetric) |
clay |
Soil clay fraction | % |
silt |
Soil silt fraction | % |
sand |
Soil sand fraction | % |
TN |
Soil total nitrogen | g/kg |
OC |
Soil organic carbon | g/kg |
CN |
Soil carbon-to-nitrogen ratio (OC/TN) | Ratio |
master_standardized_data.xlsx
Contains all columns from sem.csv plus the following additional columns:
| Column Name | Description | Unit / Values |
|---|---|---|
treatment |
Unique treatment identifier combining year, location, habitat, and replicate | Text (e.g., "2019_YL_CA_1") |
HR-standardized |
Plant species richness standardized by habitat and site (same as plant.richness_standardized in sem.csv) |
Continuous |
FR-standardized |
Arthropod taxa richness standardized by habitat and site (same as arthropod.richness_standardized in sem.csv) |
Continuous |
microhabitat |
Microhabitat category | "Open" or "Shrub" |
Study Design
- Sites: Three sites along a mean annual precipitation gradient in the drylands of northern China:
- YL (Yulin): ~400 mm mean annual precipitation
- YC (Yanchi): ~292 mm mean annual precipitation
- ST (Shapotou): ~186 mm mean annual precipitation
- Habitats: At each site, three habitat types were sampled:
- OP: Open grassland (no shrub cover)
- AT: Artemisia ordosica shrub patches
- CA: Caragana korshinskii shrub patches
- Replication: 5-6 replicate plots per habitat per site per year
- Years: 2019, 2020, 2021
- Total observations: 144 plots
How to Reproduce the Analyses
Requirements
- R (version 4.0 or higher recommended)
- Required R packages:
lme4,glmmTMB(GLMM fitting)ggplot2,patchwork(visualization)dplyr,tidyverse(data manipulation)emmeans,multcomp(post-hoc comparisons and significance letters)DHARMa(model diagnostics)car(Type II ANOVA)boot(bootstrap)readxl(reading Excel files)piecewiseSEM(structural equation modeling)
Install all packages with:
install.packages(c("lme4", "glmmTMB", "ggplot2", "dplyr", "emmeans",
"DHARMa", "multcomp", "patchwork", "car", "boot",
"readxl", "piecewiseSEM", "tidyverse"))
Running the Scripts
-
Extract
Dryad_resubmission.zipand place all data files (.xlsxand.csv) and R scripts (.R) in the same working directory. -
Set the working directory in R to that folder:
setwd("/path/to/your/directory") -
Run each script in numerical order (scripts are independent but numbered for logical flow):
source("1code-plant_diversity_among_habitat.R") source("2code-plant_diversity_along_aridity.R") source("3code-arthropod_diversity_among_habitat.R") source("4code-arthropod_diversity_along_aridity.R") source("5code-coupling_among_habitat.R") source("6code-coupling_along_aridity.R") source("7code-pSEM.R")
Output Files
Each script generates:
- PNG and/or PDF figures in the working directory
- CSV result tables (estimated marginal means, pairwise comparisons, incidence rate ratios, significance letters, etc.)
- Console output with model summaries, diagnostics, and final results
The pSEM script (7code-pSEM.R) additionally outputs sem_coefficients.csv containing standardized path coefficients for all three models.
Notes
- All data files included in this submission are sufficient to reproduce every analysis; no external data sources are required.
- The
coupling2.csvandsem.csvfiles are derived subsets ofmaster_standardized_data.xlsx, with column names adjusted to match the expectations of the analysis scripts. - Habitat factor levels are consistently ordered as OP, AT, CA across all scripts for consistent visualization and post-hoc comparisons.
- The coupling index may contain exact values of 0 (complete decoupling) and 1 (perfect coupling); scripts using beta regression adjust these boundary values slightly (to 1e-6 or 0.9999) as required by the beta distribution.
