1 Introduction

Seagrass meadows are valuable coastal and marine ecosystems on the planet, providing a range of critical environmental, economic, and social benefits (Cullen-Unsworth and Unsworth, 2018; Jayathilake and Costello, 2018). While covering approximately 0.1% of the Earth's seafloor, seagrass meadows support a wide range of biodiversities, stabilize sediment, filter water, provide coastal protection, produce more oxygen than rainforests, form the basis of the world's primary fishing grounds, and play a vital role in mitigating climate change and stabilizing the carbon cycle (Duarte et al., 2013; Unsworth et al., 2018). However, seagrass has been disappearing globally at a rate of 7% per year since 1990. Now totally more than 29% of seagrass beds have been lost, and approximately 14% of the species are at risk of extinction (Waycott et al., 2009; Duarte et al., 2013). In addition, seagrasses have received comparatively little consideration in scientific research and conservation, and seagrass meadows are among the world's least-known ecosystems (Unsworth et al., 2018).

Globally, the cumulative impacts of multiple stressors, such as coastal development, population growth, serious pollution, and climate change, play an important role in determining the distribution of the seagrass meadows at present and in the future (Adams et al., 2020). Several previous studies focused on seagrass meadow monitoring (Short et al., 2014; Telesca et al., 2015), genetic evolution (Olsen et al., 2016; Kendrick et al., 2017), species and distribution management, and restoration practices (Kenworthy et al., 2018), have provided scientific basics for the protection of seagrass population. However, increasing recognition of the conservation importance of seagrasses is calling for accurate estimates of global seagrass distribution. Thus, scientists and ocean and coastal managers have focused on mapping the distribution of seagrasses (Short et al., 2011) by using different methods, such as remote sensing (Gumusay et al., 2019), aerial photography, and underwater videography (Schultz, 2008). With the increasing model development and data accumulation in recent years, computerized modeling has become a significant research method for seagrass population management and conservation (Jayathilake and Costello, 2018; Staehr et al., 2019). Jayathilake and Costello (2018) modeled the global distribution of seagrass biome by combining over 43000 occurrence records and 13 environmental variables in the Max-Ent modeling software with 30 arcsec resolution and updated its predicted distribution area to 1646788 km². By rationalizing and updating a range of existing datasets of seagrass distribution around the globe, McKenzie et al. (2020) estimated with moderate-to-high confidence that the global seagrass area is 160387 km².

Obviously, the distribution of seagrasses is influenced by different parameters, such as physical variables (e.g., water depth, temperature, salinity, water clarity, wave height, and photosynthetically active radiation), chemical parameters (e.g., pH, phosphate, nitrate, and dissolved oxygen concentration), and biological factors (competition, predation, and genetic factors) (Unsworth et al., 2019). Different seagrass species usually have relatively specific distribution areas; for instance, the eelgrass Zostera marina L. is widely distributed in shallow waters of the northern hemisphere (Olsen et al., 2016), whereas Posidonia oceanica is an endemic seagrass to the Mediterranean Sea (Houngnandan et al., 2020). Given the successful application of species distribution models (SDMs) in the marine realm (Martínez et al., 2018), the variability of the offshore environment and the specificity of seagrass distribution make the model requirements more complicated (Assis et al., 2018). In particular, climate change is the most widespread threat to seagrass ecosystems, and the distribution pattern of seagrass and its relationship with marine environment changes are key problems that need to be solved (Chefaoui et al., 2018; Olsen et al., 2018; Unsworth et al., 2019).

Compared with studies in Australia, Europe, and North America, studies on seagrass in China are still at the starting stage (Larkum et al., 2018). In addition, the lack of basic information hinders national conservation and restoration programs for seagrasses (Zheng et al., 2013). Several researchers suggested that four families, 10 genera, and 22 species of seagrasses are distributed in Chinese coastal waters, with a predicted area of 87.65 km² (Zheng et al., 2013; Meng et al., 2019). In recent years, the researchers from China have participated in some international seagrass monitoring programs (e.g., SeagrassNet) and investigated seagrass population dynamics (Zhang et al., 2019), physiological ecology and evolutionary mechanism (Yu et al., 2018), and protection and restoration (Wang et al., 2019; Jiang et al., 2020). In 2015, a National Science and Technology Basic Work Program (2015FY110600) was funded for a national seagrass habitat survey in China, and many new recorded seagrass populations were found in different sea areas (Zhang et al., 2019). However, an evaluation of seagrass distribution across China has not been conducted, and the magnitude and direction of changes in seagrass remain unclear.

The present study aimed to provide a nationwide seagrass distribution map and to explore the dynamic changes of the seagrass population under the condition of climate change in China. We applied an SDM approach that combines information from different marine environmental conditions and seagrass occurrence data. This study represents the first attempt to map the environmental suitability for seagrasses in China under current conditions and future climate scenarios, which allows a fine-scale definition of priority areas in seagrass conservation measures for now and in the future.

2 Methods 2.1 Species Occurrence Data

As mentioned in the Introduction, 22 species of seagrasses from the families Cymodoceaceae, Hydrocharitaceae, Ruppiaceae, and Zosteraceae are distributed in Chinese coastal waters. Hence, this study focused on species belonging to these families (Table 1).

Seagrass occurrence data were from the Global Biodiversity Information Facility (GBIF, 2017) and Ocean Biogeographic Information System (Wang et al., 2021). Taxonomic names were reconciled with the World Register of Marine Species (Horton et al., 2020) and the Atlas of Seagrass (Green and Short, 2003). In case the insufficient occurrence records in China will reduce the simulation accuracy, the global seagrass occurrence data were downloaded for this study. In addition, we focused the current analysis during 2011 – 2020 to match the time scale of marine environmental datasets. Occurrence data were accessed as species-specific tables of latitude/longitude and imported into ArcGIS version 10.4 (ESRI, Redlands, California) with a WGS84 coordinate system (Bittner et al., 2020). After data preprocessing for excluding records with coordinate uncertainty, records falling on land, and duplicate records, the remaining 58961 species occurrence data were used in the subsequent modeling (Fig.1).

2.2 Marine Environmental Data

Twenty-seven abiotic variables related to the distribution of seagrasses were chosen as modeling parameters from the Global Marine Environment Datasets (GMED), which are the publicly available climatic, biological, and geophysical environmental layers featuring the present, past, and future environmental conditions (Basher et al., 2018; Table 2). All the layers were annual averages calculated over decades with a five arc-min spatial resolution (approximately 9.2 km at the equator) initially (Basher et al., 2018; Table 2). All interpolated layers were delimited to 100 m-depth thresholds for seagrass occurrence along Chinese coastal waters based on GEBCO_2019 Grid (http://www.gebco.net/) and then re-interpolated to 30 arc seconds (Jayathilake and Costello, 2018) to ensure model accuracy.

The past (last glacial maxima, 22 Myr) and future (the year 2100) abiotic layers were also from GMED, both of which contained four parameters (Basher et al., 2018) and were preprocessed as described above.

2.3 Modeling

The MaxEnt species distribution modeling software (version 3.4.1) was used to generate the seagrass SDM (Phillips et al., 2006, 2017). During data exploration, explanatory variables were assessed for correlation. When two variables were correlated at Pearson's r > 0.7, the variable which was the most proximal in determining the distribution of seagrasses was selected (Ferrari et al., 2018; Table 3). To obtain alternate estimates of which variables were the most important for different families and all seagrass species, we performed a jackknife test to improve predictive performance. In general, variables with a permutation importance of zero were deleted, and the remaining ones were more likely to be directly relevant to the species being modeled (Phillips et al., 2017).

In the current study, MaxEnt models were generated using 10 cross-validated replicate runs with the following parameters: random test percentage = 25, regularization multiplier = 1, max number of background points = 10000, maximum iterations = 1000, and convergence threshold = 10⁻⁵. MaxEnt models for each family and all seagrass species were generated separately to identify distribution differences between families. Post image processing for the map was carried out using ArcMap version 10.4, and areas with a probability value > 0.5 were taken as the potential distribution of seagrasses (Jayathilake and Costello, 2018).

Projections for the past (last glacial maxima, 22 Myr) and future (the year 2100) were applied to make reliable predictions of distribution under climate change. The potential distribution maps were cropped to Chinese coastal waters, and the geographic extent was 18˚ – 42˚N, 108˚ – 124˚E. Finally, centroids of the past, present and future projections were determined. The distance and direction between centroids were used to estimate how the distribution range will shift in respect to the present one. The receiver operating characteristic and area under the receiver operating characteristic curve (AUC) were assessed to test the performance of the model (Phillips et al., 2006).

3 Results 3.1 Environmental Variables

Correlation analysis of the environmental variables showed that seven variables, including particulate inorganic carbon, particulate organic carbon, primary productivity, slope, total suspended matter, wave height, and euphotic layer depth, had significant correlations ($\left| r \right|$ > 0.7), which were deleted in the subsequent modeling process (Table 3). Euphotic layer depth was significantly positively correlated with wave height, total suspended matter, and mean sea surface temperature (r = 0.99, 0.96, and 0.91, respectively) while negatively correlated with maximum sea surface temperature (r = −0.95). Similarly, significant positive correlations were found between wave height and total suspended matter (r = −0.99), surface current and slope (r = −0.86), mean sea surface temperature and maximum sea surface temperature (r = 0.87), particulate inorganic carbon and particulate organic carbon (r = 0.80), salinity (r = 0.94), silicate (r = 0.92), nitrate and slope (r = 0.93), and silicate and primary productivity (r = 0.90) (Table 3). In addition, significant negative correlations were observed between maximum sea surface temperature and total suspended matter, wave height, and euphotic layer depth (r = −0.95, −0.96, and −0.95, respectively; Table 3).

The relative importance of the environmental variables showed obvious divergences between each family and all seagrasses (shown in Fig.2). For Cymodoceaceae, water depth (42.4%), minimum sea surface temperature (19.5%), and nitrate (17.2%) had the highest contribution in seagrass distribution modelings. For Hydrocharitaceae, the most important variables were water depth (34.1%), calcite (15.8%), and dissolved oxygen (13.7%). Nevertheless, mean sea surface temperature (31.4%), land distance (16.4%), and phosphate (14.5%) were much more important for Ruppiaceae distribution modeling. For Zosteraceae and all seagrasses, land distance (39.3%), minimum sea surface temperature (25.4%), maximum sea surface temperature (18.3%), water depth (46.7%), mean sea surface temperature (22.3%), and land distance (19.1%) were the three most critical variables, respectively (Fig.2). Furthermore, variables with a permutation importance of 0 were deleted to perform specific distribution modeling for each family. Hence, the MaxEnt model showed that the most effective single variable of predicting the distribution of seagrasses had obvious inter-family differences in Chinese coastal waters.

3.2 Predicted Distribution and Areas

The MaxEnt model for each family and all seagrasses demonstrated high predictive strength, indicated by AUC = 0.961 ± 0.036 (mean ± standard deviation; Table 4).

The MaxEnt model for family Cymodoceaceae indicated a probability of distribution from 0.500 to 0.584, with a predicted area of 20.54 km² along the southern Hainan and western Taiwan coast (Table 4 and Fig.3A). In 2100, the distribution areas of Cymodoceaceae will decrease slightly to 16.26 km² (Table 4 and Fig.3C), and no suitable area in China was found during the last glacial maximum (Table 4 and Fig.3B).

The family Hydrocharitaceae showed good agreement with Cymodoceaceae but had a larger and relatively stable distribution area (about 33.38 km²) in China at present and in the future (Table 4, Figs.3D and 3F).

The family Ruppiaceae was mainly distributed in the coastal areas along the East China Sea, and the predicted area was 74.47 km² at present, which might decrease by about 10% in 2100 (Table 4, Figs.3G and 3I).

The family Zosteraceae was predicted to have the largest potential distribution in China coastal waters. It was more than the sum of the three other families (412.59 km²; Table 4), and mainly distributed in the Bohai Sea, the Yellow Sea, and the East China Sea (Fig.3J). The predicted distribution showed a trend of migration to the north by 2100, and the area would reduce to 372.36 km² (Fig.3L). In addition, Zosteraceae was the only family predicted to be distributed in China during the last glacial maximum, with the area reaching 10.27 km² (Table 4 and Fig.3K).

When considering all the seagrass species, the predicted areas would increase to 790.09 km² and reach 923.62 km² by 2100 (Table 4, Figs.3M and 3O).

Notably, the maps derived from the four families showed an under-prediction of seagrass area, especially in many coastal waters of the Bohai Sea, the Yellow Sea, and the South China Sea (Figs.3A – 3D). Importantly, the predicted distribution areas for nearly all families showed a decreasing tendency by 2100 (Figs.3K – 3N). Nevertheless, the map derived from seagrasses had a larger predicted distribution area, with an increasing tendency by 2100 (Fig.3E and 3O).

3.3 Centroid Migration of Suitable Habitat for Seagrass

Except for Hydrocharitaceae, all the other families tend to migrate northwest in the future. For individual families, Cymodoceaceae migrated nearly 200 km, which was the one that might move the longest distance. When modeled using all seagrass species, we obtained a relatively consistent result. The centroid of suitable habitat would shift slightly to the northwest (Fig.4).

4 Discussion

Researchers and managers have been calling for increased frequency and accuracy in the mapping of seagrass distribution to benefit seagrass conservation for over a decade (Green and Short, 2003; Bittner et al., 2020). As mentioned above, Chinese researchers have exerted remarkable efforts on seagrass ecosystems conservation and restoration in recent years (Xu et al., 2021). A high-accuracy nationwide seagrass distribution map will promote the management and protection of seagrass in China. Basing on marine environmental datasets and seagrass occurrence data, we attempted to predict the potential distribution of known seagrass species in Chinese coastal waters and reveal the past and future trends based on the MaxEnt model. However, we used the global datasets as a proxy because of the inadequacy of Chinese marine environmental data and seagrass monitoring data.

For the present, the predicted distribution area was far beyond the distribution range of seagrasses known in China (Zheng et al., 2013). The basic monitoring is not well developed, and previous studies may have greatly underestimated the distribution area of seagrasses in China. In addition, many new seagrass meadows have been discovered in several provinces in recent years (Zhang et al., 2019). In terms of the potential distribution area, our model was consistent with the known distribution of seagrasses in China, such as the specific distribution of Cymodoceaceae in Hainan and Taiwan and the distribution of Zosteraceae along the northern coast of China (Zheng et al., 2013). However, some hotspots of seagrass were missed in our model, such as the Miaodao Islands in the Shandong Peninsula for Zosteraceae (Wang et al., 2019). In the present study, all the models showed high AUC values, which indicated good prediction accuracies. However, similar to the results of Jayathilake and Costello (2018), the individual family models could under-predict the distribution of seagrass owing to the variations in species environmental niches, and the model using all species records provided the most spatially accurate map. For instance, distribution models can be strongly influenced by records. Additional data directly from Chinese coastal areas will greatly improve the accuracy of SDMs.

Physical variables (e.g., temperature, salinity, depth, wave height, and photosynthetically active radiation) and chemical parameters (e.g., pH and nutrients such as phosphate, nitrate, and dissolved oxygen concentration) control the distribution of seagrasses (Waycott et al., 2009; Jayathilake and Costello, 2018). Our analysis of variable contributions showed that the physical ones (e.g., depth, land distance, and sea surface temperature) had the greatest contribution in predicting seagrass distributions, regardless of families or seagrasses, which was consistent with the results of other studies (Jayathilake and Costello, 2018; Bittner et al., 2020). However, these results were obviously diverse in different spatial scales for different seagrass species when considering a single variable. For instance, benthic light availability represented one of the most influential predictors in the Texas Gulf Coast (Bittner et al., 2020), whereas maximum sea surface temperature had the greatest contribution in predicting global species distributions (Jayathilake and Costello, 2018). Importantly, our model clearly showed that seagrass is likely to appear in shallow areas (Staehr et al., 2019), where light conditions are within the optimal range for this species (Krause-Jensen et al., 2021).

In addition to physical variables, nutrient concentrations ranked among the key influential predictors for habitat suitability in our study area because nutrient availability controls seagrass growth, distribution, and metabolism (Papaki et al., 2020; Leiva-Dueñas et al., 2021). Chlorophyll-a concentration also had an important influence on the distribution of seagrasses, especially for Cymodoceaceae. Chlorophyll-a concentration is the key parameter determining seagrass presence/absence in certain areas. However, some issues cannot be ignored. For instance, sediment conditions are critical for determining the possible coverage of seagrass (Staehr et al., 2019; Stankovic et al., 2021). Sediment properties were not included in the dataset used in this study, which may have influenced the model results to some extent.

Accurate and effective evaluation of suitable habitats for species is essential for the identification and determination of priority protected areas (Rangel et al., 2018). Modeling can be employed to estimate the spatial distribution of potential habitat areas, where conditions in the past or in the future can be estimated. In recent years, SDMs have become an indispensable tool to describe the complex ecological relationships between species and their environment and to predict a distribution over multiple spatial and temporal scales. Meanwhile, SDMs allow us to identify the potential distribution range of seagrasses under changing conditions and to explore the impact of global climate change. MaxEnt (Phillips et al., 2006) uses presence-only data and outperforms other algorithms when applied to small datasets (de la Hoz et al., 2019). However, few studies focused on the potential distribution prediction of marine organisms in China. A good survey of the current status of marine biological distribution, integration of various biological and environmental datasets, and in-depth research on species distribution models will lay a theoretical and scientific foundation for coastal biological resources and ecological and environmental protection.

On the one hand, a correlation analysis was conducted to remove highly correlated variables, thereby improving modeling accuracy and avoiding common mistakes. In addition, variables with a permutation importance of 0 were deleted in the final models. For SDMs, the variable choice is important because inadequate variables reduce model accuracy and affect subsequent conclusions. For different families and seagrasses, the environmental variables used in the models were dissimilar and ultimately affected the simulation results. On the other hand, the environmental variables were re-interpolated to obtain higher spatial resolution. The availability of user-friendly, high-resolution environmental datasets is important for distribution modeling in aquatic environments (Assis et al., 2018). On the basis of the study by Jayathilake and Costello (2018), the spatial resolution of layers used in our study was 30 arcsec, which is reasonable for the representation of seagrass distributions at regional scales, especially when calculating the distribution areas.

The temporal resolution of environmental layers is another issue to be considered, especially when paying attention to dynamic patterns of distributions. In this study, the dataset used was averaging environmental variables over large time periods (Basher et al., 2018), which may mask the underlying dynamic patterns and produce a less realistic model (Fernandez et al., 2017). Several studies have suggested that the use of contemporaneous environmental data, such as daily or weekly fields, is preferable to fitting and projecting models on coarse-scale climatological fields, particularly in highly dynamic domains (Forney et al., 2012). Biologically relevant time scales can vary from thousands of years to minutes, indicating the need to elucidate how different temporal scales might affect SDMs in the marine realm (Fernandez et al., 2017).

Climate change is a growing concern for seagrass ecosystems (Short et al., 2016; Edwards, 2021). Obvious evidences and previous results showed that dramatic alterations driven by global climate change have occurred in oceans, and such changes are the principal cause of seagrass ecosystem degradation (Wernberg et al., 2016; Perry et al., 2019). Overall, the potential distribution of seagrass varies in different regions in the future. For instance, Chefaoui et al. (2018) indicated a dramatic loss of seagrass habitat under projected climate change in the Mediterranean Sea. By contrast, our results showed an increasing trend for seagrass in Chinese coastal waters by 2100. Although limited by the quality of environmental datasets, the present study still provides a scientific basis for the dynamic changes in seagrass distribution under global climate change. Interestingly, the suitable habitat for almost all seagrass species will shift northwest in the future despite the difference in travel distances. This finding might be evidence that seagrass meadows are moving to the polar or deep sea in response to global warming.

The distribution limit of seagrass changed tightly in China under climate change, especially for Zosteraceae. In specific, the southern limit of Zosteraceae might move 300 km to the north in the future. Similarly, Valle et al. (2014) indicated that increasing seawater temperature would trigger a poleward shift of Zostera noltii habitat suitability. Meanwhile, Franssen et al. (2011) proposed transcriptomic resilience to global warming in the seagrass Zostera marina. This phenomenon did not appear when other families or all seagrasses were used in modeling. Many studies have indicated that ocean warming will continue to drive latitudinal abundance shifts in marine species (Hastings et al., 2020; Trisos et al., 2020; Nguyen et al., 2021). Further work is warranted to elucidate the impacts of global climate change on seagrass distribution.

Awareness of the presence of seagrass meadows is the first step for protection, restoration, and sustainable management (Staehr et al., 2019; Unsworth et al., 2019). Combined with in situ investigation, this model provides a tool to identify the potential regional occurrence of seagrass at a national scale and highlight the areas where restoration efforts are likely to be successful (McKenzie et al., 2020). This modeling technique is also affected by autocorrelation because species presence sampling is inherently biased (Townhill et al., 2018). Consequently, identifying where seagrasses do not occur is also very important because many areas might have not been observed, e.g., vast areas of South East Asia (McKenzie et al., 2020). Finally, our model will be strengthened when high-quality data layers (both for marine environment datasets and seagrass occurrence data in local areas) are available and will benefit from any advances in ecological theory, statistical methods, and increases in computational power.

5 Conclusions

This study is the first to explore the potential distribution of seagrass meadows in China by using species distribution models and examine their shift trends under global climate change. Results showed that the potential distribution of seagrass is much larger than currently known distributions, suggesting that seagrass monitoring and protection still have a long way to go. Equally noteworthy is that the suitable habitat for seagrasses will shift northwest in the future, which will provide fresh insights for developing management policies and conservation strategies for seagrasses according to the climate changes. Aside from high modeling accuracy, the lack of high-quality seagrass occurrence data and marine environment data is still a key constraint for accurately predicting seagrass distribution in China.

Acknowledgements

This work was supported by the National Key R & D Program of China (No. 2019YFC1408405-02), the National Natural Science Foundation of China (No. 6207070555), and the Youth Foundation of the Shandong Academy of Sciences (No. 2019QN0026).