Modelling the potential impacts of climate change on the distribution of ichthyoplankton in the Yangtze Estuary, China

Species distribution models (SDMs) are an effective tool to explore the potential distribution of terrestrial, freshwater and marine organisms; however, SDMs have been seldom used to model ichthyoplankton distributions, and thus, our understanding of how larval stages of fishes will respond to climate change is still limited. Here, we developed SDMs to explore potential impacts of climate change on habitat suitability of ichthyoplankton.

Ichthyoplankton, comprising fish eggs, larvae and juveniles, plays an essential role in the recruitment success, temporal and spatial variations of fish populations, as well as in marine food webs (Miller & Kendall, 2009;Richardson, 2008). Being characterized by weak swimming abilities, ichthyoplankton possess a high intrinsic vulnerability to ambient stimuli, such as predation, starvation and variations in physical conditions; therefore, these planktonic stages should be more sensitive to future climate change than their adult counterparts (Hunter, 1981;Miller & Kendall, 2009;Pankhurst & Munday, 2011). A number of experimental and field studies have demonstrated that climate change will affect various aspects of ichthyoplankton, from their metabolism and phenology to their overall abundance and general distribution (Asch, 2015;Edwards & Richardson, 2004;Hsieh, Kim, Watson, Di Lorenzo, & Sugihara, 2009;Hsieh, Reiss, Hewitt, & Sugihara, 2008;Pankhurst & Munday, 2011;Walsh et al., 2015). For instance, it has been demonstrated that an increase in water temperature leads to higher egg mortality, shorter incubation periods, increased metabolic rates and enhanced developmental rates in ichthyoplankton (reviewed in Pankhurst & Munday, 2011). Moreover, using larval fish data from a 50-year temporal series, Hsieh et al. (2009) found that climate change has led to substantial changes in their distribution and abundance in Southern California. Considering the ecological importance of ichthyoplankton in marine ecosystems and the severity of climate change impacts therein, it is of utmost importance to understand how future climate change will influence ichthyoplankton for the purpose of better protecting marine fisheries resources (Richardson, 2008;Robinson et al., 2011).
Species distribution models (SDMs) can evaluate habitat suitability for target species, by relating species distribution data to a set of explanatory predictors, such as climatic and topographic variables (Franklin, 2009;Guisan, Thuiller, & Zimmermann, 2017). For example, these models have been used to evaluate habitat suitability under present and future climate scenarios for a wide variety of organisms, including terrestrial (Dyderski, Paź, Frelich, & Jagodziński, 2018), subterranean (Mammola, Goodacre, & Isaia, 2018), freshwater (Capinha, Leung, & Anastácio, 2011) and marine species (Robinson et al., 2011). Future maps of habitat suitability based on these models can effectively be used to make inference-based decisions about conservation strategies, especially for designing future areas of conservation and, ultimately, for the purpose of better protecting and managing marine fisheries resources (Robinson et al., 2011). Thus far, the majority of SDMs applied to marine ecosystems have focused on the adult stage of marine organisms, while limited effort has been devoted to modelling the distribution of ichthyoplankton (but see Pattrick, Strydom, Harris, & Goschen, 2016;Sundblad, Härmä, Lappalainen, Urho, & Bergström, 2009;Vanhatalo, Veneranta, & Hudd, 2012). Moreover, in spite of the pervasive effect of climate change across all components of marine life (Poloczanska et al., 2013), SDMs have seldom, if ever, been employed to examine climate change impacts on ichthyoplankton distributions.
Here, we constructed SDMs for five ichthyoplankton species in the Yangtze Estuary, China (Figure 1), and investigated potential climate-induced changes in their habitat suitability. As the largest estuary in China, the Yangtze Estuary provides a variety of important ecosystem services, including spawning, nursery and foraging habitat for a wide range of fisheries species (Luo & Shen, 1994). This estuary represents a coherent biogeographic F I G U R E 1 Sampling location in the present study. Sampling stations 1-34 and station 40 are in the open sea; sampling stations 35-39 are in the river area that can be used as a semi-closed model system in which to explore the potential effects of climate change on the distribution of ichthyoplankton species. Indeed, the Yangtze Estuary has experienced an accelerated rate of warming exceeding the global average (Belkin, 2009). Additionally, Yang et al. (2016) (He et al., 2008). Hypoatherina valenciennei is distributed throughout the south-western Pacific and as far north as Japan (Ivantsoff & Kottelat, 1988); this species is a common and sometimes dominant species in the Yangtze Estuary (Zhang, Xian, & Liu, 2015, 2019. Being one of the target species of bottom trawling (Li, Shan, Jin, & Dai, 2011), L. polyactis is an economically important fish resource. In China, this species is mainly distributed in the Bohai Sea, the Yellow Sea and the East China Sea. Salanx ariakensis is an annual fish distributed in estuaries and adjacent coasts from the Yellow Sea to the East China Sea (Hua et al., 2009). Chelidonichthys spinosus is an offshore warm temperature fish which could be found in China coastline, representing a keystone element in the marine ecosystem food web (Li, Xu, Jiang, & Zhu, 2010).
Forty sampling stations (five stations within the river and thirty-five outside the river mouth) were established in the Yangtze Estuary and monitored between 1999Estuary and monitored between and 2013Estuary and monitored between (05/1999Estuary and monitored between , 05/2001Estuary and monitored between , 05/2004Estuary and monitored between , 05/2007Estuary and monitored between , 05/2009Estuary and monitored between , 05/2010Estuary and monitored between , 05/2011Estuary and monitored between , 05/2012Estuary and monitored between and 05/2013. The spatial distribution of the sampling stations is given in

| Environmental explanatory variables
Marine ecosystems are exposed to a wide range of environmental change stressors, including ocean warming, acidification and hypoxia (Gruber, 2011;Poloczanska et al., 2013). In addition, the distribution of ichthyoplankton can be influenced by multiple physical or chemical factors, particularly tides, dissolved oxygen and water temperature (Miller & Kendall, 2009;Pattrick et al., 2016;Sundblad et al., 2009;Vanhatalo et al., 2012). In the present study, a total of 13 marine predictor variables were initially considered: water depth, distance to shore, pH, annual mean and range of sea surface temperature, annual mean and range of sea surface salinity, annual mean and range of sea surface current velocity, annual mean and range of sea surface dissolved oxygen, and annual mean and range of sea surface primary productivity. All predictors were obtained from the Bio-ORACLE v2.0 dataset (http://www.bio-oracle.org; Assis et al., 2018), except data for water depth and distance to shore, which were derived from the Global Marine Environment Datasets (http:// gmed.auckl and.ac.nz; Basher, Bowden, & Costello, 2014). The spatial resolution of all predictors was 5 arcmin, corresponding to 9.2 km at the equator. We checked multicollinearity among these 13 predictor variables by calculating pairwise Pearson's correlation coefficients (r), using an |r| > 0.70 to cull collinear predictors (Dormann et al., 2013). According to the results of the collinearity analysis ( Figure   S1) and data availability under both present and future time periods, four predictors were retained to model habitat suitability for each ichthyoplankton species: distance to shore, annual mean sea surface temperature, annual mean sea surface salinity and annual mean sea surface current velocity.
Projections of annual mean sea surface temperature, salinity and current velocity for the future [time periods 2040-2050 (2050s) and 2090-2100 (2100s)] under four different emission scenarios-representative concentration pathways (RCPs; RCP26, RCP45, RCP60 and RCP85)-were derived from three atmosphere-ocean general circulation models (AOGCMs): CCSM4, HadGEM2-ES and MIROC5 (Assis et al., 2018). Average outputs of the three AOGCMs were used as future climate conditions. These were also retrieved from the Bio-ORACLE v2.0 dataset (Assis et al., 2018). We assume that distance to shore remains constant in future projections.

| Modelling procedure
An ensemble modelling approach was adopted in this study to reduce uncertainties associated with single modelling algorithms (Araújo & New, 2007;Qiao, Soberón, & Peterson, 2015 We fitted these ten algorithms with their default parameters using the package biomod2 (Thuiller, Georges, & Engler, 2014) in the R software environment (version 3.4.3; R Development Core Team, 2017). Binary data (i.e. presence-absence of species) are needed for several SDM algorithms; owing to the small number of true absence records, we used pseudo-absence points instead (Guisan et al., 2017;Thuiller et al., 2014). Pseudo-absence points were randomly sampled from the study region; for each species, the number of pseudo-absences equalled ten times the number of presences (Barbet-Massin, Jiguet, Albert, Beaumont et al., 2009). Predictive abilities of the ten algorithms were evaluated using the true skill statistic (TSS; Allouche, Tsoar, & Kadmon, 2006) and the area under the receiver operating characteristic curve (AUC; Swets, 1988) via a fivefold cross-validation technique (Guisan et al., 2017;Thuiller et al., 2014). As with Engler et al. (2011) and Araújo, Pearson, Thuiller, and Erhard (2005), algorithms with values of TSS ≤ 0.4 and AUC ≤ 0.7 were disregarded to minimize uncertainty resulting from modelling algorithms with poor predictive ability.
Relative importance of the four explanatory variables was determined by a randomization procedure (see Guisan et al., 2017;Thuiller et al., 2014). Habitat suitability for five ichthyoplankton species in the Yangtze Estuary, under present and future climates, was evaluated using all binary data. Projections of selected single modelling algorithms were ensemble by a committee averaging technique (Guisan et al., 2017;Thuiller et al., 2014). The projected habitat suitability values range from 0 to 1,000 with 0 representing the lowest occurrence probability (i.e., 0) and 1,000 representing the highest occurrence probability (i.e., 1). For each species, the projected continuous probability maps were converted into binary presence-absence maps by selecting a probability threshold maximizing the TSS value (Franklin, 2009;Guisan et al., 2017;Jiménez-Valverde & Lobo, 2007;Thuiller et al., 2014). In addition to biomod2 package, different R packages were used in our analyses for data manipulation and exploratory analyses, including corrplot (Wei & Simko, 2017), maptools (Bivand & Lewin-Koh, 2013), maps (Brownrigg, Minka, & Deckmyn, 2018) and raster (Hijmans, 2019).

| Model performances and the importance of explanatory variables
Predictive performances of the ten modelling algorithms varied depending on the ichthyoplankton species considered (Figure 2).

| Habitat suitability under present and future climate scenarios
The predicted habitat suitability for the five ichthyoplankton species was generally consistent with each species' known distribution in the Yangtze Estuary (Figure 3). The predicted climatically suitable areas for H. valenciennei, L. polyactis and S. ariakensis were relatively larger than the areas for C. mystus and C. spinosus ( Figure 3). The species C. mystus is predicted to have high habitat suitability close to the mouth of the Yangtze Estuary, while the other four species are expected to have high habitat suitability in the inner estuary (Figure 3). The ensemble of SDM projections suggested that future climate scenarios will have different impacts on habitat suitability for the five ichthyoplankton species: C. mystus was consistently projected to expand its range under future climates, while distributions of the other four species will likely contract in the future (Table 3). Given the identical trends in changes of habitat suitability under the four different emission scenarios (Table 3), projections under the RCP45 scenario-a midrange emission scenario-were displayed to show the potential climate change impact on the distributions of the ichthyoplankton species.
Projected changes in occurrence probability of the ichthyoplankton species suggest that, among the five species tested here, there will likely be two different responses to climate change in both 2050s (Figure 4) and 2100s ( Figure 5). Climates in the 2050s and 2100s will seemingly favour the expansion of C. mystus in the Yangtze Estuary, and habitat suitable for this species is predicted to still occur near the estuary mouth (Figure 4a, Figure 5a). In contrast, in addition to their range contractions, the other four species are predicted to shift their distributions northward in response to future climate change, meanwhile losing suitable habitats in the current core of their distribution (Figures 4b-e, 5b-e). Accordingly, owing to changes in the distribution patterns of ichthyoplankton caused by climate change, the biodiversity of ichthyoplankton, as represented by the five species considered, is projected to overall decline in the Yangtze Estuary in the future ( Figure 6).

| D ISCUSS I ON
In this contribution, we used ensemble SDMs to explore present and future habitat suitability for five ichthyoplankton species in the Yangtze Estuary, the largest estuary in China. Our study represents The ecological importance of ichthyoplankton to marine ecosystems is widely recognized, as well as their great vulnerability to climate change (Asch, 2015;Hsieh et al., 2009;Hunter, 1981;Miller & Kendall, 2009;Pankhurst & Munday, 2011;Poloczanska et al., 2013;Richardson, 2008). Thus, a few previous studies have proposed the potential utility of SDMs in investigating present and future habitat suitability for ichthyoplankton (Dambach & Rödder, 2011;Richardson, 2008). Despite this attention, only few SDM studies to date have focused on planktonic larvae, and these were exclusively devoted to studying the present-day habitat suitability based on a single SDM algorithm (Pattrick et al., 2016;Sundblad et (Araújo & New, 2007;Qiao et al., 2015;Thuiller et al., 2019). A multi-model ensemble approach, whereby the predictions of multiple modelling algorithms are synthesized, is typically used to minimize these uncertainties (Araújo & New, 2007;Guisan et al., 2017;Thuiller et al., 2019). In our case, the ensemble SDMs for five ichthyoplankton species exhibited predictive abilities superior to those of single algorithms, and the ensemble of SDM projections of habitat suitability under the present climate was highly consistent with the known distributions of each species.
Despite the potential significance, our modelling approach has at least two limitations. First, we should notice that in addition to environmental variables, other factors including species dispersal capacity are also important in regulating species distributions (Guisan et al., 2017). Ichthyoplankton dispersal is a complex process which can be influenced primarily by passive drift with currents and secondarily by active swimming (Miller & Kendall, 2009) (Hernandez, Graham, Master, & Albert, 2006;Wisz et al., 2008). To improve model accuracy, it would be useful to establish more sampling stations in the Yangtze Estuary.
The five species we studied are common species in the Yangtze Estuary and have great ecological and economic importance. These waters. Conversely, C. mystus is an estuarine migratory fish that commonly lives in shallow marine habitats but migrates to brackish estuarine waters in spring as sexually mature individuals (He et al., 2008). Previous studies have proposed that euryhaline species F I G U R E 4 Changes in occurrence probability of five ichthyoplankton species in 2050s under RCP45 scenario. Shades of blue indicate areas in which the probability of occurrence will decrease, and vice versa for red areas. Insets on the top right of each graph represent suitable habitats under present-day and future climates. Red areas are projected to be suitable in the future, green areas are projected to be suitable under both present-day and future climates, and blue areas represent present-day suitable habitat that will become unsuitable in the future are generally eurytopic and should be more resistant to environmental stresses than stenohaline marine organisms (Boesch, 1974;Wright, Kennedy, Roosenburg, Castagna, & Mihursky, 1983). Tang, Hu, and Yang (2007)  The molecular work would support our SDM projections that the response pattern of C. mystus to future climate change is range F I G U R E 5 Changes in occurrence probability of five ichthyoplankton species in 2100s under RCP45 scenario. Shades of blue indicate areas in which the probability of occurrence will decrease, and vice versa for red areas. Insets on the top right of each graph represent suitable habitats under present-day and future climates. Red areas are projected to be suitable in the future, green areas are projected to be suitable under both present-day and future climates, and blue areas represent present-day suitable habitat that will become unsuitable in the future F I G U R E 6 The sum of projected presences of the five ichthyoplankton species considered in this study under present conditions and future (RCP45) climate change scenarios. The continuous SDM predictions were converted into binary maps by selecting probability thresholds maximizing the True Skill Statistics (TSS) value (see Guisan et al., 2017;Thuiller et al., 2014). Map scale ranges from 0 (no species is predicted to be present) to 5 (all five species are predicted to be present). (a) present period, (b) RCP45 in 2050s (2040-2050), (c) RCP45 in 2100s (2090)(2091)(2092)(2093)(2094)(2095)(2096)(2097)(2098)(2099)(2100) expansion. Yet, it would be also useful to investigate physiological variations among different ichthyoplankton species in the Yangtze Estuary to further corroborate this projection.
Projected changes in habitat suitability induced by climate change will likely result in losses and northward shifts of the ichthyoplankton biodiversity in the Yangtze Estuary. This conclusion supports previous studies that found support for changes in species distribution and biodiversity as a result of climate change (Hsieh et al., 2009(Hsieh et al., , 2008Walsh et al., 2015). In addition to range shifts, changing climates can also affect species abundance (Barrett et al., 2018;Doney et al., 2012;Richardson, 2008). Declines in the abundance of ichthyoplankton have already been observed in the Yangtze Estuary (see Zhang et al., 2015Zhang et al., , 2016Zhang et al., , 2019. It is import- ant, yet especially challenging, to accurately model species abundance (Oppel et al., 2012;Pearce & Boyce, 2006). We strongly advise that further efforts be made to estimate the impacts of climate change on the abundance of planktonic larvae in this region.
It has been revealed that climate change can influence species in a wide variety of ways, such as behavioural changes, range shifts, changes in phenology and alterations in species interactions (Asch, 2015;Doney et al., 2012;Edwards & Richardson, 2004;Vergés et al., 2016;Walsh et al., 2015); hence, future studies should address  We thank Mr. Yan Qin for editing Figure 1. We express our sincere thanks to the editor and two anonymous reviewers for their valuable comments and suggestions, which helped improve this manuscript. SM is supported by Bando per l'Internazionalizzazione della Ricerca-Anno 2018 (Compagnia di San Paolo). The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript.

CO N FLI C T O F I NTE R E S T
The authors declare there are no competing interests.

B I OS K E TCH
Zhixin Zhang is interested in studying impacts of climate change on the distribution of a breadth of species, from commercially important fisheries species, along with species of conservation concern, to invasive species.

Stefano
Mammola is an ecologist whose scientific activity focuses primarily on subterranean biology and ecological modelling.
Hui Zhang is a marine biologist interested in marine fisheries diversity and ecology. He relies on multiple methods, such as field survey and molecular markers, to explore the dynamics of marine fish diversity.
Weiwei Xian is an estuarine ecologist, interested in understanding the long-term environmental changes in these peculiar and highly dynamic ecosystems.
Author contributions: H. Z. and W. X. provided data for this study; Z. Z. and H. Z. conceived the idea; Z. Z. performed data analyses and wrote the first draft of the manuscript with constructive comment and suggestions from S. M., H. Z. and W. X.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section.