The study area for the research was the region along the Western Pacific coast, between 95°E-145°E longitude and -10°N-45°N latitude. The study range for the three lineages of B. pectinirostris, AE1, AE2, and AES, was restricted to 123°E-145°E, 15°N-45°N; 95°E-123°E, 15°N-45°N and 95°E-145°E; -10°N-15°N, respectively, based on molecular phylogeography results ( Figure 1 ) (Chen, 2010; Chen et al., 2014a; Chen et al., 2014c). The occurrence records of B. pectinirostris were obtained from four sources, including the Ocean Biogeographic Information System (OBIS, http://www.obis.org), the Global Biodiversity Information Facility (GBIF, http://www.gbif.org), the iNaturalist website (http://www.inaturalist.org), and the records from the literatures ( Supplementary Table 1 ). Then the occurrence records were spatially thinned by using the R package spThin to reduce sampling bias and set at a 9.2 km radius, consistent with the resolution of the environmental layers (Aiello-Lammens et al., 2015). A total of 128 occurrence records were retained, including 21, 92, and 15 occurrence records for the AE1, AE2, and AES lineages of B. pectinirostris, respectively.

Given the availability of data and the correlation between species and environmental variables, we initially selected 10 environmental variables that could potentially impact the distribution of B. pectinirostris ( Supplementary Table 2 ). Notably, environmental variables such as water temperature, current velocity, salinity, dissolved oxygen, and ocean depth are key variables that restrict the distribution of marine fish species (Sperling et al., 2016; Bosch et al., 2017; Goldsmit et al., 2018). B. pectinirostris is herbivorous species that primarily feeds on benthic diatoms, which are usually represented by chlorophyll in marine ecology (Yang et al., 2003; Huot et al., 2007). Additionally, as B. pectinirostris mainly inhabits tidal flats, distance to shore was incorporated as a predictor variable, offering an inferential measure of the tidal flat expanses.

With a spatial resolution of 5×5 arc minutes (i.e., 9.2×9.2 km), all environmental variables ( Supplementary Table 2 ) were collected from the online database: Bio-ORACLE v2.1 (Assis et al., 2018) and the Global Marine Environment Datasets (Basher et al., 2018). The average values (monthly averages during the entire period) of environmental variables during 2000-2014 were taken to represent the present environmental conditions. The pairwise Pearson’s correlation coefficients between the 10 predictor variables were calculated in order to lessen the effect of collinearity on the evaluation of predictor variable importance ( Supplementary Figure 1 ). One variable was removed from each pair with a Pearson’s correlation coefficient greater than |0.7| (Dormann et al., 2013). The Dshore (distance to shore), Tm (mean water temperature), Sal (mean salinity), CV (mean current velocity), DO (mean dissolved oxygen), and Chl (mean chlorophyll concentration) were the final six environmental variables kept for subsequent analysis.

In addition, Bio-ORACLE provided future environmental projections using three atmospheric-ocean general circulation models (AOGCMs) under four representative concentration pathway (RCP) emission scenarios (RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5) (Assis et al., 2018). To decrease uncertainty, the average predictions of the three AOGCMs were used to represent future climatic conditions (Assis et al., 2018). RCP 2.6–RCP 8.5 represent the different scenarios of greenhouse gas emissions from low to high in the future (Moss et al., 2010). Considering the availability of future environmental data, only five environmental variables (Tm, Sal, CV, Chl, and Dshore) were used in this study for species distribution projection. And, in this study, two scenarios, RCP 2.6 and RCP 8.5, were selected to project the future distribution of B. pectinirostris in the different periods of the 2050s (2040–2050) and 2100s (2090–2100). In the future, the increase in greenhouse gas emissions will lead to a larger increase in water temperature as well as more significant decreases in salinity and chlorophyll concentration (Sun et al., 2022). Besides, this study assumed that the distance to shore would remain constant in the future (Zhang et al., 2021).

Firstly, the R package FactoMineR was used to perform a principal component analysis (PCA) on six environmental predictive variables (Dshore, Tm, Sal, CV, Chl, and DO) to reduce the dimensionality of the niche space (Le et al., 2008). The first four principal components were retained, as they accounted for over 80% of the total variance ( Supplementary Figure 2 ). Subsequently, we used the R package hypervolume to compute a four-dimensional hypervolume based on the retained principal components (Blonder et al., 2021). The size of the niche space occupied by a lineage can be indicated by the volume of the hypervolume. Finally, the hypervolumes of different lineages were overlapped using the R package BAT to assess niche differentiation among lineages (Cardoso et al., 2021). The total niche differentiation (βTotal) ranges from 0 to 1 and measures the degree of niche differentiation between two lineages (Carvalho and Cardoso, 2020; Mammola and Cardoso, 2020). Higher βTotal values indicated greater differentiation. And it consists of niche shift, which represents the replacement of space between the hypervolumes occupied by the populations, and niche contraction or expansion, which refers to the net differences in the spatial extent encompassed by each niche hypervolume (Carvalho and Cardoso, 2020; Mammola and Cardoso, 2020).

This study used 10 modeling techniques in the Biomod 2 R package (Thuiller et al., 2021). Considering the scarcity of true absence records, this study simulated 2000 pseudo-absence points in the contrasting environment conditions where the available occurrence records were not present (Barbet-Massin et al., 2012). Five environmental predictor variables (Dshore, Tm, Sal, CV, and Chl) were used in this study for modeling predictions. Since the AE1 and AES lineages of B. pectinirostris have only 21 and 15 occurrence records, this may lead to overfitting of the model. To overcome this limitation, this study applied the ensemble of small models (ESMs) strategy using the R package ecospat, which is a modeling strategy specifically designed for species with limited occurrence records (Lomba et al., 2010; Breiner et al., 2015; Di Cola et al., 2017). According to this strategy, this study selected two of the five environmental predictor variables at a time to construct many bivariate models, and then the predictions of these bivariate models were weighted and averaged according to the model performance, thus obtained the predictions of the ensemble of small models (ESM) (Breiner et al., 2015, 2018). Finally, the predictions of the ESM based on different algorithms were assembled into the final ensemble of small model (ESMs) predictions based on model performance (Breiner et al., 2015, 2018). A fivefold cross-validation strategy was used to assess the performance of the models (Guisan et al., 2017). The performance of each bivariate model and ESM was evaluated using the area under the curve (AUC) values (Lomba et al., 2010). We then used Somers’ D [D = 2 × (AUC − 0.5)] as weights to integrate the ESM into the final ensemble model (Lomba et al., 2010). In this process, any models with a Somers’ D value less than 0 (i.e., the AUC value below 0.5) will be excluded from the final ensemble model. In addition, the performance of the final ensemble model was evaluated in this study using the true skill statistic (TSS) and the Boyce index (Boyce) (Allouche and Kadmon, 2006; Hirzel et al., 2006). Ensemble models were constructed at both the species-level (species model) and the lineage-level (AE1 model, AE2 model, and AES model).

The relative importance of each predictor variable for modeling the distribution of B. pectinirostris was assessed according to the methods of Di Cola et al. (2017). The potential distribution of B. pectinirostris under different climate scenarios (RCP 2.6 and RCP 8.5) for the present and future (2050s and 2100s) was predicted using ensemble models constructed at the species and lineage levels, respectively. Continuous habitat suitability maps were generated based on the direct output of the ensemble model. In order to facilitate the interpretation of habitat suitability, this study converted the continuous habitat suitability predictions into binary maps by maximizing the probability threshold of the true skill statistic (Liu et al., 2013).

The multivariate environmental similarity surface (MESS) analysis was used to analyze the climate change in the study area in the future compared with the current time period and the key environmental factors driving climate change (Elith et al., 2010). First, we extracted the environmental data (Dshore, Tm, Sal, CV, and Chl) for the species distribution points (true habitat) in the current period as reference points. And we then used the R package ENMeval to calculate the environmental similarity (S) between all locations in the study area and the reference points for both the current period and the future period (taking the RCP8.5 scenario in the 2100s as an example) (Kass et al., 2021). The larger the S value, the higher the environmental similarity between the region and the reference points, indicating that the region may be a suitable habit for the species; when the S value is negative, it means that at least one environmental factor in the region is beyond the range of the reference set, indicating that the environment of the region and the reference point is more different, and it may be unsuitable for the survival of the species (Elith et al., 2010). Then the R package ENMeval was used to extract the key climate variable that contributes to the differences between each location within the study area during the future period and the reference set (Elith et al., 2010; Kass et al., 2021). This variable was referred to as the most dissimilar variable (MoD), which is the key environmental factor driving future climate change at a specific site (Elith et al., 2010). In addition to an overall analysis of climate change within the ranges of species, this study analyzed climate change within the respective ranges of the three lineages.

The four-dimensional hypervolumes for AE1, AE2, and AES were 693.34, 2201.73, and 627.18, respectively. Correlating with the size of their distribution ranges, the AE2 lineages with the largest distribution ranges had the largest hypervolume ecological niches. Comparison of the hypervolume results (βTotal) showed that there was marked differentiation in niches between the three lineages, and the differentiation in niches within the AE lineage was smaller than that between AE and AES ( Table 1 ). The niche differentiation between AE2 and AE1 or AES was primarily due to niche contraction/expansion (>82.3%), while the contribution of niche shift was minor (<17.7%) ( Table 1 ). The niche differentiation between AE1 and AES was mainly due to niche shift (>93.1%), while niche contraction/expansion only contributed less than 7.0% ( Table 1 ). The overall niche comparison among the three lineages indicated that niche differentiation was mainly influenced by PC1, which was mainly explained by DO, Tm, and Sal ( Figure 2 and Supplementary Figure 2 ). Pairwise niche comparisons showed that niche differentiation between AE1 and AE2 was mainly influenced by PC2, which was mainly explained by Dshore, CV, and Chl; niche differentiation between AE1 and AES was mainly influenced by PC1, which was mainly explained by Tm, Sal, and DO; niche differentiation between AE2 and AES was mainly influenced by PC1, which was mainly explained by DO, Tm, and Sal ( Figure 2 and Supplementary Figure 2 ).

The results of the performance evaluation of the ensemble of small models (ESM) constructed based on different modeling techniques were shown in Figure 3 . The results showed that the performance of these models exceeded the thresholds set in this study, and all these models will be used to construct weighted ensemble models (ESMs). The predictive performance of the four final ensemble models (ESMs) was shown in Table 2 , which indicated that all performance metrics exceeded 0.8, demonstrating the good predictive performance of the final ensemble models. The results of the lineage-level models showed that different lineages had different responses to environmental variables. The distribution of AE1 and AES lineages was greatly influenced by Dshore and Tm, while Dshore and Chl were the most important environmental variables affecting the distribution of the AE2 lineage ( Figure 4 ). In addition, the importance of Dshore exceeded the average importance of the five variables (i.e., 0.2) across all models.

In this study, we quantified and compared the niches of three lineages of B. pectinirostris, and the results showed that the niche space of the AE2 lineage was more expansive than that of the AE1 and the AES lineages. The broader intertidal zone within the AE2 lineage’s distribution region, compared to the AE1 and AES lineages (Murray et al., 2019), coupled with the extensive latitudinal gradient of the dispersal range, allowed it to encompass a wider range of climatic spectrums (Zhang and Yan, 2014). Furthermore, the three lineages of B. pectinirostris distributed in different seas had already produced local adaptations and had marked niche differentiation. And the niche differentiation between these lineages was caused by different reasons ( Table 1 ). The niche differentiation between the AE2 lineage and the other two lineages was mainly due to the expansion/contraction of the niche. Based on the previous description, this study speculated that the niche space of the AE2 lineage may include most of the niche space of the other two lineages. The AE1 lineage and the AES lineage were spatially distant from each other, and the environmental conditions of the regions where they were located differed considerably (Zhang and Yan, 2014). These different lineages had likely undergone evolutionary fine-tuning to their respective local conditions, thus settling in non-overlapping niche spaces. This divergence, predominantly attributed to niche shift dynamics, could potentially elucidate the niche differentiation observed between them.

The niche differentiation between the two lineages within the East Asian lineage was less pronounced than that between the East Asian lineage and the Southeast Asian lineage ( Table 1 ). Furthermore, the differences between them may also be reflected in morphology and genetics (Savolainen et al., 2013). Previous molecular biology studies had shown that the East Asian and Southeast Asian lineages of B. pectinirostris had a relative distant genetic relationship, and the Southeast Asian B. pectinirostris may be a cryptic species (Chen et al., 2014c). In addition, Chen (2010) mentioned that the adult fish of B. pectinirostris in East Asia were smaller in size than those in Southeast Asia. The results of the present study also supported the conclusion that the East Asian lineage and the Southeast Asian lineage of B. pectinirostris were distantly related to some extent. It can be seen that species distribution models not only play an important role in predicting potential habitats and preventing biological invasions but can also help discover new species (Tyberghein et al., 2012).

Species-level SDMs were often constructed based on the assumption of “niche conservation,” disregarding local adaptation and intraspecific variation during the natural selection process (Hällfors et al., 2016; Chefaoui et al., 2018; Zhao et al., 2020). However, an increasing number of studies have indicated that incorporating intraspecific variation and local adaptation into SDMs may improve the accuracy of model predictions (Collart et al., 2021; Zhang et al., 2021; Sun et al., 2022). Previous molecular biology studies had demonstrated genetic differentiation among the three lineages of B. pectinirostris (Chen, 2010; Chen et al., 2014a; Chen et al., 2014c), while this study has shown niche differentiation among them. This highlighted the necessity of constructing lineage-level SDMs to explain the local adaptation of different lineages.

The results of the research indicated that both species-level and lineage-level models predicted a negative effect of climate change on the AE2 lineage, with the species-level model predicting more optimistic results than the AE2 model. Some researchers have argued that species-level models may underestimate the suitable habitat range for a species as they do not take into consideration the phenomenon of local adaptation (Oney et al., 2013; Ikeda et al., 2017; Zhang et al., 2021). Nevertheless, the present study suggested that the species-level model might have overestimated the potential distribution of the AE2 lineage due to the fact that it assumed the lineage could endure all the climatic conditions that the species as a whole can tolerate. This meant that genetic variation within the species did not always mitigate the effects of climate change on its habitat suitability (Cacciapaglia and van Woesik, 2018; Hu et al., 2021; Li et al., 2022). This situation seemed to be more apparent for the AE1 lineage, as the lineage-level model predicted a contracting trend in the suitable habit range of the AE1 lineage, while the species-level model predicted an expanding trend. For the AES lineage, both the species-level and AES models projected an obvious expansion of the habitat range. Compared with the species-level model, the AES model projected a larger expansion of the habitat range, which might be attributed to the positive impact of local adaptation on the response of the AES lineage to climate change (Zhang et al., 2021; Sun et al., 2022).

Although local adaptation and intraspecific variation had different impacts on different lineages, we still tended to believe that incorporating these factors into the model could improve its accuracy (Li et al., 2022). For example, in the present study, the species-level model showed that the distribution of AE2 lineage could extend northward to the waters near Laizhou Bay and Bohai Bay in the current environmental conditions, while the lineage-level model predicted that the distribution of AE2 lineage could only extend to the waters near the southern part of the Shandong Peninsula. However, by consulting literature related to the distribution of B. pectinirostris (Wu and Zhong, 2008; Zhao, 2022), it was discovered that, possibly due to the water temperature limitations under current environmental conditions, the northern boundary of the AE2 lineage was in the waters near the southern part of the Shandong Peninsula. Thus, it could be inferred that constructing lineage-level SDMs might help to reduce prediction uncertainty when there is clear genetic or niche evidence of intraspecific phylogenetic heterogeneity (Zhang et al., 2021).

Despite the crucial role of SDMs in assessing the impact of climate change on the habitat suitability of B. pectinirostris, the present models still have some limitations (Rathore and Sharma, 2023). Although occurrence records of B. pectinirostris were collected from various sources, there were relatively few records of the AES and AE1 lineages in relevant databases and literature, which resulted in limited data for the construction of lineage-level models. When constructing species distribution models, insufficient data may cause an incomplete capture of all the niche space of a species, resulting in overfitting or underfitting of the model, which can affect the accuracy and reliability of predictions (Heikkinen et al., 2006; Elith and Leathwick, 2009). To overcome this limitation, the ensemble of small models (ESMs) strategy was used in this study to minimize the impact of limited data availability for the AE1 and AES lineages on model construction (Lomba et al., 2010; Breiner et al., 2015). Of course, if there were more occurrence data available, it could potentially lead to more precise predictions (Elith and Leathwick, 2009).

The three distinct lineages of B. pectinirostris exhibited geographic isolation across expansive oceanic expanses, resulting in their restricted genetic exchange (Chen, 2010; Chen et al., 2014a; Chen et al., 2014c). Long-term local adaptation may lead to niche differentiation among different lineages, and thus they respond differently to climate change (Sun et al., 2022). In terms of the degree of influence of environmental variables on the distribution of B. pectinirostris, both the results of the SDMs and the MESS analysis revealed that distance to shore was the most significant factor. This phenomenon elucidated that B. pectinirostris had exclusive confinement to the intertidal vicinity contiguous to coastal margins (Murdy, 1989; Takegaki, 2008; Polgar and Crosa, 2009; Polgar and Bartolino, 2010). In addition, water temperature and chlorophyll concentration were also identified as important factors influencing the distribution of B. pectinirostris.

The results of this study indicated that the continuous rise in water temperature might lead to habitat loss for the AE2 lineage in the waters near the Gulf of Beibu and the northern South China Sea because the base temperature in these waters was already relatively higher (Zhang and Yan, 2014). Furthermore, the decline in chlorophyll concentration in the waters near the southern part of the Shandong Peninsula might result in the loss of habitat for the AE2 lineage in this area (Hashioka and Yamanaka, 2007). For the AE1 lineage, the results of the two different level models predicted were different. The species-level model predicted a potential expansion of its habitat, while the lineage model predicted a potential loss of some of its habitat due to increased water temperatures. This study suggested that this may be attributed to the species-level model being influenced by the ecological niches of the other two lineages in its projection of the distribution of the AE1 lineage, thus overestimating the tolerance of the AE1 lineage to high temperatures. Therefore, the projections of the lineage-level model constructed solely based on the ecological niche of the AE1 lineage may be more accurate (Peñalver-Alcázar et al., 2021; Li et al., 2022). Notably, the AES lineage manifested a positive response to climate change. Building upon the results of this study, it was speculated that the AES lineage may have fostered heightened resilience to higher water temperatures due to its long-term adaptation to the high-temperature environment of Southeast Asia.

B. pectinirostris has important ecological and economic value, and its unique amphibious habits have also attracted considerable scholarly attention (Wu and Zhong, 2008; Aiello et al., 2023; Bian et al., 2023). During the 1970s and 1980s, the dramatic decline in resource abundance of the Japanese population of B. pectinirostris was primarily caused by habitat loss due to land reclamation and environmental pollution, coupled with the effects of overfishing (Takegaki, 2008). In the Ministry of the Environment Red List 2020 of Japan, its conservation status was listed as Endangered (EN) (http://www.env.go.jp). In China, the relatively wide distribution and abundance of resources of B. pectinirostris have resulted in less attention being paid to its conservation. Moreover, previous research on the conservation of B. pectinirostris had focused on the investigation of the genetic diversity of natural germplasm resources and the development and utilization of artificial propagation techniques, while relatively few studies had been conducted on the impacts of future climate change on the habitat of the mudskipper (Takegaki, 2008; Chen et al., 2014a; Chen et al., 2014b). Therefore, the results of this study will provide a new perspective for the conservation of B. pectinirostris. The results of this study suggested that the AES lineage might have greater resilience to climate change and might experience habitat expansion, whereas the AE1 and AE2 lineages might face habitat loss in certain regions. Considering the vulnerability of the AE1 and AE2 lineages to climate change, conserving these two lineages should be given a higher priority. In addition, the nearshore habitats of B. pectinirostris were also strongly affected by anthropogenic influences, which would exacerbate habitat loss and fragmentation and threaten its survival (Yoo et al., 2013; Murray et al., 2019; Chen et al., 2022). In conclusion, we recommended continuous monitoring of the distribution and population density of B. pectinirostris, as well as the development of more effective management measures, such as limiting fishing quotas and establishing marine protected areas, to ensure the long-term sustainability of this species.