Dataset DOI: 10.5061/dryad.h9w0vt4wz

Description of the data and file structure

This README describes the supporting datasets, analytical outputs, R codes, supplementary materials, and methodological information associated with the manuscript:

“Climate change will exacerbate the mismatch between plant diversity hotspots and protected areas in the Hindu Kush Himalaya.”

Prepared for Submission to: Dr. Yad

1. Overview

This dataset supports the study evaluating spatial and temporal patterns of vascular plant diversity, phylogenetic diversity, and endemism across the Hindu Kush Himalaya.

The dataset includes:

• 464,407 vascular plant distribution records;

• 22,058 accepted vascular plant species;

• a completed analytical phylogeny containing 22,058 species;

• current, Last Glacial Maximum, and future climate data;

• spatial taxonomic and phylogenetic diversity metrics;

• Random Forest and Generalized Linear Model projections;

• multidimensional biodiversity hotspot analyses;

• current-period CANAPE analysis;

• protected-area gap assessment;

• species-level complementarity analysis;

• exploratory Partial Least Squares Path Modeling (PLS-PM); and

• ecological vulnerability assessment.

  External data sources are listed below with links/DOIs. Users should obtain raw data directly from the original providers under their terms.

  Source	Use in this study	Link / DOI	Included in Dryad?
  POWO	Taxonomic standardization	https://powo.science.kew.org/	No raw data; standardized names only
  GBIF	Occurrence records	https://www.gbif.org/ ; download DOIs are listed in Supporting Information Table S2	No raw records; derived grid presence only
  WorldClim 2.1	Current and future bioclimatic variables	https://www.worldclim.org/data/worldclim21.html	No raw layers; derived grid summaries
  WorldClim CMIP6 / MIROC6	Future projections under SSP2-4.5	https://www.worldclim.org/data/cmip6/cmip6_clim30s.html	No raw layers; derived projections
  PaleoClim / NCAR-CCSM4	Last Glacial Maximum bioclimatic variables	https://www.paleoclim.org/	No raw layers; derived grid summaries
  WDPA	Protected-area gap analysis	https://www.protectedplanet.net/en/thematic-areas/wdpa?tab=WDPA	No raw boundaries; derived protected according to study area map.
  PREDICTS	Biodiversity Intactness Index	PREDICTS database (Hudson et al., 2017)	No raw data; derived grid values
  NASA LP DAAC MODIS	NDVI, EVI, and vegetation quality layers	NASA LP DAAC at the USGS EROS center (MODI3Q1), NASA LP DAAC at the USGS EROS center (MODI3A2), NASA LP DAAC	No raw layers; derived grid values
  CGIAR-CSI	Aridity index	Global Aridity Index and Potential Evapotranspiration database version 3.1 (Zomer et al., 2022)	No raw layers; derived grid values
  Human footprint	Anthropogenic pressure	Gassert et al. 2023 (Global 100 m resolution dataset)	No raw layers; derived grid values
  Population density	Anthropogenic pressure	Global 1-km gridded population dataset (L. Liu et al., 2024)	No raw layers; derived grid values
  Habitat heterogeneity	Habitat heterogeneity indices	Metrics of spectral habitat heterogeneity. (Tuanmu & Jetz, 2015)	No raw layers; derived grid values
  Land use / land cover	Global Human Modification	Land Use/Land Cover classification (Kennedy et al., 2020)	No raw layers; derived grid values
  R packages	Analysis and visualization	CoordinateCleaner, V.PhyloMaker, picante, Biodiverse, caret, glmulti, iml, VennDiagram, ggupset, dplyr, tidyr, lmtest, broom	Code only

2. Core Workflow Summary

The analysis was conducted through the following sequential workflow.

2.1 Species-occurrence compilation

Plant distribution records were compiled from:

• national and regional floras;
• taxonomic checklists;
• herbarium datasets;
• published occurrence databases; and
• Global Biodiversity Information Facility (GBIF) downloads.

The complete list of national and regional floras and all GBIF download DOIs is provided in the Supporting Information.

2.2 Taxonomic standardization

Species names and synonyms were standardized against Plants of the World Online, facilitated by the Royal Botanic Gardens, Kew.

Invalid coordinates, zero coordinates, non-numeric coordinates, and marine records were removed. Verified occurrence records were converted into presence records within 50 × 50 km grid cells across the Hindu Kush Himalaya.

2.3 Phylogenetic-tree preparation

The source phylogeny contained 22,152 species and was derived from a previously published dated phylogeny of the Pan-Himalayan flora.

After matching the occurrence dataset with the source phylogeny:

• the pruned backbone tree contained 19,944 species;
• 2,114 missing species were incorporated using V.PhyloMaker; and
• the final analytical phylogeny contained 22,058 species.

Ten replicate completed trees were generated to evaluate sensitivity to random phylogenetic grafting.

2.4 Spatial phylogenetic analyses

Grid-level taxonomic and phylogenetic diversity metrics were calculated using Biodiverse and R.

The calculated indices included:

• species richness;
• phylogenetic diversity;
• phylogenetic endemism;
• weighted endemism;
• relative phylogenetic diversity;
• relative phylogenetic endemism;
• standardized effect size of phylogenetic diversity;
• net relatedness index; and
• nearest taxon index.

2.5 Environmental variables and temporal projections

Current, Last Glacial Maximum, and future climate variables were prepared for diversity projections.

Random Forest-based variable selection and Generalized Linear Models were used to project four principal diversity indices:

• species richness;
• phylogenetic diversity;
• phylogenetic endemism; and
• weighted endemism.

Random Forest models were fitted using 500 trees. Predictor selection was performed using Random Forest-based recursive feature elimination with ten-fold cross-validation.

Model performance was assessed using:

• coefficient of determination (R²);
• root mean square error (RMSE); and
• mean absolute error (MAE).

Future projections were generated after refitting the selected models using the complete current dataset.

2.6 Hotspot identification

Multidimensional biodiversity hotspots were identified based on the spatial overlap of:

• species richness;
• phylogenetic diversity;
• phylogenetic endemism; and
• weighted endemism.

Two hotspot thresholds were assessed:

• top 5 %, representing the highest-priority hotspot areas; and
• top 10 %, representing broader hotspot regions.

Temporal changes in hotspots were evaluated through:

• spatial loss;
• gain;
• retention;
• turnover;
• directional shifts; and
• elevational shifts.

2.7 CANAPE analysis

CANAPE was conducted only for the current period to identify centers of:

• neo-endemism;
• paleo-endemism; and
• mixed-endemism.

These analyses characterized current evolutionary patterns but were not incorporated into the Level I or Level II conservation-gap framework.

2.8 Conservation-gap analysis

Current and future multidimensional biodiversity hotspots were overlaid with the existing World Database on Protected Areas.

Conservation gaps were defined as hotspot areas occurring outside the current protected-area network.

2.9 Complementarity analysis

Species-level complementarity analysis was conducted to evaluate representation of plant species within existing protected areas and identified conservation gaps.

The analysis:

• identified species with less than 10 % of occupied grid cells represented within protected areas;
• quantified their representation within Level I and Level II conservation gaps;
• identified narrow-endemic species occupying five or fewer grid cells; and
• evaluated evolutionary distinctiveness among poorly represented species.

2.10 Driver and vulnerability analyses

Exploratory Partial Least Squares Path Modeling (PLS-PM) was used to examine relationships among:

• evolutionary factors;
• vegetation;
• climate;
• topography;
• anthropogenic pressure;
• species richness; and
• phylogenetic endemism.

An ecological vulnerability index (EcoVI) was calculated using pressure, resilience, and sensitivity indicators.

The abbreviation EVI refers only to the remotely sensed enhanced vegetation index.

3. Repository structure

The complete supplementary materials are provided in:

Supplementary_Information_and_R_codes_R2.zip

This archive contains two folders:

Folder 1: R codes

The R scripts used for data processing, analyses, and visualization are provided as follows:

Code 1_VphyloMaker
R code for phylogenetic tree preparation and incorporation of missing species using V.PhyloMaker.

Code 2_SESPDNTIandNRI
R code for calculating phylogenetic structure metrics, including SES.PD, NTI, and NRI.

Code 3_Prediction-LGM-RF
R code for Random Forest modelling, recursive feature elimination, GLM modelling, and Last Glacial Maximum and future diversity projections.

Code 4_Upset plot and Venndiagrams
R code for generating UpSet plots and Venn diagrams.

Code 5_Plspm model
R code for Partial Least Squares Path Modeling (PLS-PM) analysis.

Code 6_EVI
R code for calculating the ecological vulnerability index (EcoVI).

Code 7_Complementarity analysis
R code for species-level complementarity analysis and conservation-gap assessment.

Folder 2: Supporting Information files 1–24

This folder contains supplementary datasets, tables, figures, and methodological information associated with the manuscript.

Supporting Information 1
Species list and locality records.

Supporting Information 2
List of fossil calibrations and records.

Supporting Information 3
Tactual_HKH.tre — initial phylogenetic tree for the Hindu Kush Himalaya dataset.

Supporting Information 4
List of species incorporated using V.PhyloMaker.

Supporting Information 5
Full completed phylogenetic tree.

Supporting Information 6
Variables used for the 2100 projection.

Supporting Information 7
Current plant diversity indices.

Supporting Information 8
Predictions from the Last Glacial Maximum model.

Supporting Information 9
Predictions for 2021–2040.

Supporting Information 10
Predictions for 2041–2060.

Supporting Information 11
Predictions for 2061–2080.

Supporting Information 12
Predictions for 2081–2100.

Supporting Information 13
Projected plant diversity indices for the Last Glacial Maximum and future periods.

Supporting Information 14
Gain and loss data among different periods.

Supporting Information 15
Multidimensional hotspot classifications for all periods using top 5 % and top 10 % thresholds.

Supporting Information 16
CANAPE categories.

Supporting Information 17
Phylogenetic Jaccard groups and clustering results.

Supporting Information 18
Protected and non-protected area grids.

Supporting Information 19
Variables used in correlation analysis for PLS-PM.

Supporting Information 20
Variables used in ecological vulnerability index (EcoVI) analysis.

Supporting Information 21
Complementarity analysis summary.

Supporting Information 22
Endemism summary.

Supporting Information 23
Poorly represented species.

Supporting Information 24
Poorly represented species with endemism information.

4. Variable descriptions

Variables occurring repeatedly across different temporal periods follow the same naming conventions.

Diversity variables

Temporal abbreviations

Environmental variables

Environmental predictors retain their original names from the source datasets and analytical files.

Examples:

Projection variables

Equivalent naming conventions apply for 2060, 2080, and 2100.

Conservation and complementarity variables

Variable definitions, transformations, and analytical procedures are further described in the associated manuscript Methods and Supporting Information.