2.1 Study area
We surveyed 30 continental shelf islands (Fig. 1 and Table S1) in southeastern Xiangshan County, Zhejiang Province, eastern China (28°55′N ~ 29°25′N, 121°45′E ~ 122°10′E). The distinct geographical and climatic conditions among 30 islands contribute to notable plant diversity variations. Recent research have indicated that the biodiversity of this region is shaped by both natural elements and human activities (Hu et al. 2011; Yu et al. 2019). The average area of surveyed islands is 0.1812 km², ranging from the largest at 1.3371 km² to the smallest at 0.014 km². The terrain is predominantly hilly, composed mainly of laterite and yellow soil. Vegetation types commonly encompass evergreen coniferous and broadleaf forests, alongside shrub-herb communities. The region is characterized by a subtropical monsoon climate with distinct seasons. The frost-free period reaches 248 days. The average annual temperature is 18.3°C, while rainfall exceeds 1574.8 mm, and annual sunshine ranges from 1482 to 1656 hours.
2.2 Plant survey
The study of island biogeography generally takes the species richness of the whole island as the research object (Triantis et al. 2012; Schrader et al. 2020). This study is also based on the analysis of the species richness of different life-forms of plants on the entire island, therefore, we need to investigate the plant species on the entire island. According to studies on species-area relationships in subtropical regions, the minimum sample volume required to describe community structure and diversity in forest ecosystems is 20 m × 20 m (Song et al. 2013). In addition, in order to fully investigate the plant species on the islands, we set up more quadrats (Fig. 2a) on larger islands by referring to the method of Xu et al. (2023). Based on the above considerations, we set three plots in the larger islands (area ≥ 1 km2). One exposed on the windward slope of the island, while another on the leeward slope, and the third is at the highest elevation position. In contrast, two plots were positioned at the summit and mid-slope, respectively, in middle-sized islands (0.5 ≤ area < 1 km²). For small islands (area < 0.5 km²), only a single plot was placed at the summit. Within each plot, five 1 m² sub-plots were designated for herb surveys at four corners and the center. This sampling design yielded 36 woody plant plots across the 30 islands, with a total of 180 herbaceous subplots. In view of the fact that the species of the entire island cannot be fully investigated by the quadrat, and the difference in the number of quadrats between different islands may also lead to changes in the number of island species. This study also adopts the transect survey method (Fig. 2b) by referring to Yu et al. (2012) as a supplementary means, the specific method was to set transect lines around the island and along the middle of the island, and record the plant species that appeared within 5 meters around the transect lines. Repeated species are no longer recorded until no new species appear. At the same time, in inaccessible areas such as the edge of the island and cliffs, UAV aerial photography (Fig. 2e-2g) is used to assist the investigation. Data obtained through the above three methods can effectively survey species richness data across the island, thereby reducing the impact of quadrat number differences on survey results. Most of the plant species in the survey were directly identified and recorded through field observation. For species that could not be accurately identified in the wild, specimens were collected and sent to the Botanical Geography Laboratory at Ningbo University for further identification. The identification of species and the classification of life-forms were referred to the Atlas of Plants in Ningbo (Li et al. 2021) and the Flora Reipublicae Popularis Sinicae (http://www.iplant.cn/) and alien invasive plants were referenced from the Checklist of the Chinese alien invasive plants (Yan et al. 2014).
2.3 Island attributes, human activities and climatic factors
We sourced island latitude, longitude, area, and elevation from the Ningbo Institute of Oceanography (Table S1). Based on these data, we computed the remoteness and shape indexes to measure island accessibility and shape complexity, respectively. Area size and remoteness are the drivers of island biogeography (MacArthur and Wilson, 1963, 1967), but the shape complexity is also considered to be an important factor affecting island biodiversity because islands with a high value have more habitat within them and can accommodate more species (Xie et al. 2023). Our shape complexity is replaced by shape index. The minimum value of shape index is 1, which indicates that the island shape is a perfect circle. After that, as the island shape becomes more complex and irregular, its shape index values continue to escalates, indicating higher diversity of island habitats (Hoffmeister et al. 2005). The shape index is determined using Formula 1. C denotes the island’s perimeter, while A signifies the area. Similar with the previous studies, area, remoteness, elevation and shape indexes are summarized as island attributes that affecting diversity change (Hu et al. 2011).
$$SI= \frac{C}{\left[2*{\left(\pi *A\right)}^{0.5}\right]}$$
1
As suggested by previous research (Chen et al. 2021; Xu et al. 2023), we employed geographic distance (DM /km), a minimum distance from the island’s center to the nearest mainland coastline (Li et al. 2019), to quantify island remoteness using the “sf” package in R version 4.3.0. Additionally, to accurately reflect island remoteness, according to the method of Weigelt and Kreft (2013), we also calculated the stepping-stone distance (SDM / km) that accounting for both the direct mainland distance (DM) and the proximity to the nearest large island. After preliminary analysis, we found that SDM has high collinearity with DM. Additionally, these two measures of remoteness have similar impacts (Fig. S1). Thus, we only used DM as a representative of remoteness. Importantly, we observed there no significant correlation between remoteness (DM or SDM) and island area (r = -0.196, p = 0.3; r = -0.295, p = 0.113), corroborating the independence of island area effects and remoteness effects.
We use land use intensity and habitat heterogeneity to reflect human activities on islands (Xu et al. 2023; Liu et al. 2023). Land use intensity (LUI) quantifies human impact on the land surface (Foley et al. 2005). As suggested by Xu et al. (2023), the LUI is calculated by the percentage of buildings, farmland, and orchards that encompass the total area of the islands, which was assessed using the ArcGIS 10.4 based on the PlanetScope Scene remote sensing imagery with 3 m spatial resolution. The Shannon index in (2) was used to calculate the habitat heterogeneity (H) of each island (Formula 2). Where m is the number of land use types on an island, pi is the proportion of land use type i in the island area. A large H-value suggests a high measure of habitat heterogeneity on the island.
$$H=- \sum _{i=1}^{m}{p}_{i}\text{ln}{p}_{i}$$
2
Similar with numerous previous studies, we obtained the bioclimatic variables (Bio1–Bio19) in 30 m resolution from WorldClim v2.0 (http://www.worldclim.org) to explore their link to species diversity (Kreft et al. 2008; Elith and Leathwick 2009; Rijsdijk et al. 2014). These variables are annual mean temperature (Bio1), mean diurnal range (Bio2), isothermality (Bio3), temperature seasonality (Bio4), maximum temperature of warmest month (Bio5), minimum temperature of coldest month (Bio6), temperature annual range (Bio7), mean temperature of wettest quarter (Bio8), mean temperature of driest quarter (Bio9), mean temperature of warmest quarter (Bio10), mean temperature of coldest quarter (Bio11), annual precipitation (Bio12), precipitation of wettest month (Bio13), precipitation of driest month (Bio14), precipitation seasonality (Bio15), precipitation of wettest quarter (Bio16), precipitation of driest quarter (Bio17), precipitation of warmest quarter (Bio18), and precipitation of coldest quarter (Bio19). All the above variables selected for this study are shown in Table S2.
2.4 Statistical analysis
According to the previous studies, a power function model with logarithmic transformation was constructed to delineate the species-area curve for species richness (total, tree, shrub and herb) relative to island area (Formula 3) (Triantis et al. 2003, 2012; Losos and Ricklefs 2009). Where A denotes the island area, S is species richness, z and logc are the fitting parameters for the slope and intercept, respectively. z denotes the rate of change in species richness with area, and is often used as a measure of the sensitivity of species and communities to fragmented habitats (Yu et al. 2012).
$$\text{log}S= \text{log}c+z*\text{log}A$$
3
Since the species-area curve were only applicable to island areas, we also used linear regression to analyze the relationship of varied richness (total, tree, shrub and herb) with remoteness, habitat heterogeneity and LUI. Standardized major axis regression was used to determine the differences in the relationships between richness and both island area and remoteness among the three plant life forms.
After that, according to Sanchez (2013) and Wang et al. (2016), the Partial Least Squares Path Modeling (PLS-PM) was used to elucidate the influencing mechanism of species richness. To test which variables had a significant effect on species richness, we examined the correlation between island attributes (area, remoteness, elevation and shape index), human activities (LUI and H), and climatic factors (Bio1-Bio19) with species richness using Pearson's correlation coefficients (p < 0.05, Fig. S2), and excluded uncorrelated variables. Subsequently, we conducted a collinearity diagnosis using the Variance Inflation Factor (VIF), only retaining those with VIF values below 10 (Table S3). After these two calculations, the seven factors including island attributes (area and remoteness), human activities (LUI and H) and climatic factors (Bio4, Bio7 and Bio18) were screened out to construct the PLS-PM for exploring the driving mechanism of biodiversity.
Finally, we employed the Akaike Information Criterion (AIC) to select the optimal multiple linear regression model among the remaining factors, which was used in constructing the pathways for PLS-PM (Table S4). In the modeling, the same variables were categorized and then performed 5,000 iterations to refine the path coefficients and their significance. The total effect of each influencing factor was used to measure its influence on richness (Sanchez. 2013; Wang et al. 2016). PLS-PM modelling fit was evaluated using the Goodness of Fit index (GoF). The GoF value ≥ 0.36 indicates better model alignment (Henseler and Sarstedt 2013). All data analyses were completed in R version 4.3.0. Linear regression, the standardized major axis regression and PLS-PM were conducted using the “ggplot2”, “smart” and “plspm” packages, respectively.