Shandong is bordered by the Bohai Sea and the Yellow Sea ( Figure 1 ). It boasts a coastline spanning 3,345 km and shallow sea areas within a depth of 20 m, covering 29,731 km². This area accounts for 37% of the total area of the Bohai Sea and the Yellow Sea (Peoples Government of Shandong Province, 2023), making it highly suitable for marine aquaculture and has the most production of Pacific oysters in the world. The rivers in Shandong particularly the Yellow River, transport significant amounts of freshwater and sediment to the Bohai Sea, affecting the salinity, nutrient, and sediment concentration near the estuary (Zheng et al., 2021).

The offshore areas of Shandong are primarily influenced by the north Shandong coastal current, which originated from the Bohai Sea and flows eastward to the Yellow Sea eastward. This current runs parallel to the northern coastline of the Shandong peninsula, turning southward along the eastern coast before eventually moving southwestward along the southern coast of the Shandong peninsula. The north Shandong coastal current correlates well with wind speed (Zheng et al., 2021), becoming more pronounced in winter when northerly winds prevail (Zhang et al., 2018), causing sea level rise in LZ Bay (Li et al., 2015), and transporting large amounts of freshwater and sediment from the Yellow River eastward along the northern coast (Yang et al., 2011; Wang et al., 2020), resulting in low salinity and high turbidity in the coastal waters along the north coast of the Shandong Peninsula (Yang et al., 2011).

The climate of Shandong is controlled by the East Asian monsoon (Song et al., 2021). Winter winds are predominantly blow from the north while summer winds are predominantly southerly. This monsoonal climate is characterized by distinct seasons, higher summer temperatures with concentrated precipitation (60%–70%), with the temperature gradients increasing from the southeast coast to the northwest inland and precipitation patterns showing the opposite trend (Peoples Government of Shandong Province, 2023).

The environment, climate, meteorology, and socioeconomic data utilized in this study, along with their sources and resolutions, are listed in Table 1 . To assess the spatiotemporal distribution of suitable areas for Pacific oyster marine aquaculture along the Shandong offshore, we acquired and processed monthly average data of factors from May 2011 to December 2022 and correlated these with climate events. The climate data used for future climate impact analysis, including its source, model, scenarios, time series, and resolution, are detailed in Table 1 . Among them, 12 factors, including Chl-a, TSS, SST, DO, pH, SO, WW, Res, VO, WS, Tas, and Pr, have time series data from May 2011 to December 2022. Bathymetry, bioclimatic, distance to city, pier, and WWTPs lack time series data. Mean data will be used for both time series and climate scenarios.

The water quality factors considered in this study included chlorophyll-a (Chl-a), total suspended sediment (TSS), sea surface temperature (SST), dissolved oxygen concentration in seawater (DO), seawater pH reported on the total scale (pH), seawater salinity (SO). Chl-a and TSS data were obtained from the Geostationary Ocean Color Imager (GOCI) onboard the Korean geostationary orbit satellite. Their data were downloaded from the Korea Ocean Satellite Center (http://kosc.kiost.ac.kr/index.nm) with a resolution of 500 m for Level 1 data. The GOCI Data Processing System 2.0 (GDPS 2.0) software’s Batch Process tool was used to calculate Chl-a and TSS bands. ENVI 5.3 software was used for geometric correction based on the official GLT files. After GOCI concluded its observation mission in March 2021, Chl-a and TSS data from April 2021 to December 2022 were obtained from the second Korean geostationary orbit satellite (GOCI-II). This study used Level 2 data with an initial resolution of 250 m. Monthly average data were computed using SeaDAS 8.0 software by selecting cloud-free data each month. SST data were sourced from the Moderate Resolution Imaging Spectroradiometer (MODIS) aboard NASA’s AQUA sun-synchronous polar-orbiting satellite. Data were obtained from the Ocean Color website (https://oceancolor.gsfc.nasa.gov/) with a resolution of 1 km for Level 2 data. Daily data were filtered to minimize cloud cover impact, atmospheric data was corrected using the OCSSW tool in SeaDAS 8.0, and monthly averages were calculated using ArcGIS 10.7 with interpolation and other processing performed in SeaDAS. SO, DO, and pH were obtained from the E.U. Copernicus Marine Service Information (CMEMS) Marine Data Store (MDS).

Hydrological factors include bathymetry, water velocity (VO), and sea surface wind wave significant height (WW). Bathymetry data were sourced from the General Bathymetric Chart of the Oceans (GEBCO, https://www.gebco.net/), offering a global resolution of approximately 15 arc seconds (500m). VO and WW were obtained from the CMEMS MDS. VO was calculated from the square root of the sum of the squares of the eastward and northward sea water velocities.

Meteorological data include wind speed at 10 m above the surface (WS), total precipitation (Pr), and 2-meter dewpoint temperature (Tas). They were sourced from the European Centre for Medium-Range Weather Forecasts’ ERA5 reanalysis monthly average data, these data were interpolated for missing values using MATLAB.

Climate data, the climate scenario data in 2010–2040 for SST and Tas are based on the projected temperature increases associated with the SSP1–2.6 and SSP5–8.5 mentioned in the sixth assessment reports of the Intergovernmental Panel on Climate Change (IPCC). The climate scenario data in 2040–2070, 2070–2100 for SST, Tas, and the future predictions for VO are sourced from the Biogeographic Oceanic and Regional Seas Environmental Predictive Model (Bio-ORACLE). Nineteen bioclimatic parameters (Karger et al., 2017) were derived from Climatologies at High resolution for the Earth’s Land Surface Areas (CHELSA), which offers high-resolution data (30 arc seconds, ~1 kilometer). We selected the GFDL-ESM4 model, the SSP1–2.6 and SSP5–8.5 scenarios in four periods: 1981–2010, 2011–2040, 2041–2070, and 2071–2100, these two scenarios can produce a more pronounced contrast to the suitability. The Oceanic Niño Index (ONI) and the Multivariate ENSO Index (MEI) (Wolter & Timlin, 2011) were used for climate events. ONI is available from NOAA’s Climate Prediction Center (https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php). The MEI can be obtained from NOAA’s Physical Sciences Laboratory (https://psl.noaa.gov/enso/mei/ ).

Socioeconomic data included the location of cities, piers, and wastewater treatment plants (WWTPs). Cities data were taken from global multi-temporal urban boundary data (Gong et al., 2020), which includes all global cities and surrounding settlements over 1 square kilometer, effectively capturing the contours of urban-rural edge areas. Pier data were extracted from coastal Points of Interest (POI) available in relevant map services. WWTP data came from the global database of HydroWASTE (Ehalt Macedo et al., 2022). Socioeconomic data on suitability distribution were obtained using the Euclidean Distance tool in ArcGIS 10.7.

All data were resampled to 500 m resolution using ArcGIS 10.7. Due to varying data sources, to ensure temporal continuity, 1000 random samples were generated within a 20 km range of the Shandong offshore using the ArcGIS10.7 Create Random Points tool to ensure temporal continuity. The relationship between datasets was quantified through linear fitting to ensure usability (Park et al., 2021) and consistency. Environmental (except Bathymetry) and meteorological data from May 2010 to December 2020 were statistically analyzed for the five Pacific oyster marine aquaculture areas and all offshore areas within 20 km of Shandong Peninsula to understand further the spatiotemporal variation characteristics of environmental and meteorological factors in Shandong.

Figure 2 illustrates the schematic diagram of the suitability evaluation model for Pacific oyster marine aquaculture, which primarily includes four categories of factors: water quality (Chl-a, SST, TSS, SO, pH, DO), hydrology (Bathymetry, VO, WW), climate and meteorology (WS, Tas, Pr, Bioclimatic), and socioeconomic factors (distance to city, pier, WWTPs). Each factor was assigned a score ranging from 1 to 8, with 1 representing the least suitable and 8 representing the most appropriate. The model employed the Analytic Hierarchy Process (AHP), a multi-criteria decision-making method proposed by American operations researcher Saaty in the early 1970s (Saaty, 1977), to determine the weights of each factor and each sub-model. The model utilizes ArcGIS 10.7 for Reclassify, Spatial Analyst, and Model Builder to enable batch processing of time series data, generating spatial distribution maps of suitability for Pacific oyster marine aquaculture. The ArcGIS raster calculator also calculated quarterly averages for all data from May 2011 to December 2022. A comparison was made between the highly suitable areas and the marine aquaculture areas identified through high-resolution satellite imagery (Wang et al., 2018; Liu et al., 2020b) to validate the accuracy of the model.

This study analyzed the impact of climate change on Pacific oyster marine aquaculture suitability from two aspects: the effects of climate events and the future climate change on suitability. Firstly, when investigating the influence of climate events such as El Niño on the suitability of Pacific oyster marine aquaculture, we compared the suitability of average years with El Niño (2015) and La Niña (2022) and utilized time series data from May 2011 to December 2022. We conducted Pearson correlation analysis and box plot analysis between time series data and the High Suitability Percentage (HSP), ONI, and MEI indices, with significance testing based on a significance level of 0.01 and 0.05 for t-tests. The causes of temporal and spatial variations in suitability were analyzed by combining existing research with correlation analysis results.

Mean data were used to investigate the impact of future climate change on suitability. Except for the predicted values and databases (SST, SO, Tas, Bioclimatic) shown in Table 1 , all other data were based on the mean values from May 2011 to December 2022, including mean data, Bathymetry, Bioclimatic, distance to city, pier, WWTPs, as history and future periods date. Bioclimatic factors data were used with the MaxEnt model to predict suitable historical and future suitable zones. Pacific oyster distribution data were sourced from the Global Biodiversity Information Facility (GBIF) (GBIF, 2023). Due to insufficient data in the offshore area of Shandong, distribution data from the northwest Pacific region (117°E-146°E, 30°N-45°N) were used for suitable zone prediction, supplemented with in-situ oyster reef distribution area and wild Pacific oyster population sampling point data ( Supplementary Table 1 ) from literature (Gu et al., 2005; Fang et al., 2007; Quan et al., 2012; Wang et al., 2014; Ran et al., 2018; Zhong et al., 2019; Li et al., 2020; Song et al., 2021; Zhang et al., 2021; Quan et al., 2022; Hong et al., 2023). ENMTools software was used to calculate spatial autocorrelation and exclude highly correlated (>0.8) features with a lower contribution, combining historical data for the MaxEnt model suitable zone prediction. The accuracy of the MaxEnt model was evaluated using the Area Under Curve (AUC) under the Receiver Operating Characteristic (ROC) Curve. The species’ habitable probability obtained was graded to determine suitability distribution.

Figure 3A illustrated the changes in the monthly average and Figure 3B illustrated the quarterly spatial distribution of time series factors across different regions from May 2011 to December 2022. Chl-a and TSS exhibited similar spatial distribution patterns, with concentration gradually decreasing from the coast toward the open sea. In winter, the concentration of Chl-a was highest, while the lowest concentration was observed in summer. LZ had the highest Chl-a concentration, and CD had the lowest, yet all areas have Chl-a concentration above the 20 km average. In contrast, all areas had lower TSS concentrations than the 20 km average, except for JM, which had higher TSS concentrations in summer.

SST was highest in summer and lowest in winter, with minimal spatial variation and almost no inter-regional differences. LZ Bay experienced higher SST in spring, while the area near RC exhibited lower temperatures in summer. SO increased spatially from the coast toward the open sea, with LZ Bay exhibiting the lowest yearly SO but noticeable seasonal fluctuations, being higher in summer and lower in winter. The spatial and temporal pH variation was minimal, with slight differences among regions, and the entire Shandong offshore had a weakly alkaline pH. DO was higher in winter and spring but lower in summer and autumn, with LZ Bay having the highest winter DO concentration.

WW increased from the coast toward the open sea, with the lowest wave heights observed in summer. The CD region exhibited the most significant variation in WW, reflecting its location in the open sea. VO was more robust in winter and spring, with spatial differences in VO, which can identify the coastal currents in the Shandong offshore, and LZ had the lowest water flow across all years.

Pr and Tas followed similar seasonal patterns, with more rainfall in the south than in the north during spring and autumn, though regional differences were not pronounced, and CD had the least rainfall. WS showed minor spatial variation, with higher speed in summer and winter.

Table 2 outlines the final evaluation model, with grading indicators determined based on existing research and the actual conditions of the Shandong offshore. This study set a high Chl-a concentration threshold detrimental to marine aquaculture at more than 5 mg/m³ (Terauchi et al., 2014). However, the Chl-a concentration in the Shandong offshore did not exceed 3.17 mg/m³. Therefore, in the grading indicators, the higher the Chl-a concentration, the higher the score. The grading indicators for TSS were based on the experimental levels: 0, 0.10, 0.25, and 0.5 g/L (Suedel et al., 2015a), with lower TSS concentrations receiving higher scores. SST in the Shandong offshore exhibited seasonal variations, with high summer temperatures causing mass mortalities of marine aquaculture oyster (Yang et al., 2021). The grading indicators for SST were based on the Arrhenius temperature of oyster (Van Der Veer et al., 2006). The grading indicators for SO were based on Wang et al (2023), with 25‰-35‰ representing the highest suitability, decreasing linearly outside this optimal range. The grading indicators for DO were based on the concentration needed for oyster reef restoration (Wang et al., 2023). The seawater pH in Shandong ranges from 7.9–8.2, showing minor seasonal differences. The grading indicators were set that higher pH values indicate higher suitability.

As the depth increased from the coast toward the open sea in the Shandong offshore but did not overall exceed 20m, and VO did not exceed 20 cm/s (except in RC), the grading indicators were set that greater depth and water flow indicated higher suitability. The resuspension of seabed sediments caused by wind waves could also determine the best suitable area. For the resuspension effect, which had a critical depth (h < 10WW), we used WW-h/10 as a grading indicator. If the effects was harmful, it could not cause resuspension; if it was positive, it could cause resuspension and was considered a positive factor, providing more nutrients. However, waves could negatively impact rafts, so a lower WW was considered better for suitability (Ogle et al., 1977; Goseberg et al., 2017).

WS was graded from less suitable to more suitable in an arithmetic sequence, as the average monthly WS in the Shandong offshore did not exceed 5 m/s, making its potential impact on marine aquaculture rafts relatively small. Tas had the same grading indicators as SST, based on the Arrhenius formula (Van Der Veer et al., 2006). Pr was graded from less suitable to more suitable in a linear relationship. After correlation analysis ( Supplementary Figure 1 ) and contribution contrast ( Supplementary Table 2 ), the bioclimatic factors finally selected include bio2 (mean diurnal range), bio5 (max temperature of the warmest month), bio8 (mean temperature of the wettest quarter), bio10 (mean temperature of the warmest quarter), bio14 (precipitation of the driest month), and bio16 (precipitation of the wettest quarter). The training and test set values for both historical and future periods were above 0.9 ( Supplementary Figures 2, 6, 10, 14, 18, 22, 26 ), indicating good prediction results and high reliability. Response curves ( Supplementary Figures 3, 7, 11, 15, 19, 23, 27 ), Variable contributions ( Supplementary Tables 3, 4, 5, 6, 7, 8, 9 ) and jackknife test of variable importance ( Supplementary Figures 4, 8, 12, 16, 20, 24, 28 ) of historical and future periods were in the Supplementary Material. The bioclimatic factors graded the habitability probability ( Supplementary Figures 5, 9, 13, 17, 21, 25, 29 ) obtained from the MaxEnt model.

In practical production, factors related to oyster harvesting, such as storage, transportation, waste disposal, and maintaining a certain distance to piers equipped with specific equipment (marine aquaculture equipment and vessels), are necessary to ensure profitability; otherwise, costs become significantly higher (Liu et al., 2014). This is particularly challenging for individual farmers and small enterprises (Barillé et al., 2020). Proximity to cities facilitates oyster transportation, processing, and sales. Our model incorporates socioeconomic factors that could quantify the specific impact of distance, such as distance to piers and cities. It also considers the negative impact of WWTPs, quantifying their impact on marine aquaculture activities using distance metrics. In socioeconomic factors, proximity to cities and piers indicated higher suitability, while proximity to WWTPs indicated lower suitability scores.

The weight value for each factor was obtained through a literature review and experts’ opinions (Radiarta et al., 2008, Radiarta et al, 2011; Saitoh et al., 2011; Cho et al., 2012; Liu et al., 2013, 2014; Aura et al., 2016; Snyder et al., 2017; Barillé et al., 2020; Liu et al., 2020a; Jiang et al., 2022). Environmental factors such as SST, Chl-a, TSS, SO, DO, and Bathymetry were commonly considered. Socioeconomic factors, including distance to the city, piers, land-based facilities, and constraints like harbors, townships, industrial areas, and river mouths, were also considered in the evaluation model.

Compared to socioeconomic factors, environmental factors hold a higher weight. Among the environmental factors, water quality, such as SST, Chl-a, and TSS, carries a significant weight of 40%. Chl-a concentration reflects the biomass and productivity of phytoplankton in seawater and can indicate food availability (Xing et al., 2017). Furthermore, it has been found to have a significant positive correlation with oyster growth (Mizuta et al., 2012). Oysters tend to increase their feeding activity to a maximum level and stabilize, increasing Chl-a concentration within a specific range (Tenore and Dunstan, 1973). TSS can negatively impact oyster pumping and clearance rates (Loosanoff and Tommers, 1948; Suedel et al., 2015). SST is a crucial parameter influencing biological processes, and high summer temperatures can lead to mass mortality in oyster marine aquaculture (Yang et al., 2021). While Pacific oyster are euryhaline species, SO does not significantly affect their growth (Nell and Holliday, 1988). However, increased freshwater influx can result in rapid, short-term salinity reduction, limiting oyster growth (Swam et al., 2022). In the Shandong offshore area, the maximum SO does not exceed 30‰, and the LZ Bay exhibits the lowest SO due to extensive river inputs, necessitating consideration of SO variations. DO in coastal ecosystems is experiencing more significant fluctuations than other environmental variables, and seasonal hypoxia events are rising (Diaz, 2001). Considering the decrease in DO with increasing marine aquaculture density, this factor needs consideration, albeit with a lower weight (Brito et al., 2023). pH reduction can impede the growth of early life stages of oyster (Ko et al., 2014). The pH value of seawater in Shandong ranges from 7.9 to 8.2, with minimal seasonal variation, warranting a lower weight.

Hydrological factors are also crucial in influencing the growth and marine aquaculture of Pacific oyster, and they hold a weight of 30% in this model, relatively lower than water quality factors. Suspended particles influenced by WS and WW, such as Chl-a and TSS, exhibit good temporal consistency with the mass concentration and diffusion intensity. Wind and waves play a significant role in the spatial and temporal distribution and diffusion of suspended particles, including Chl-a and TSS (Liu and Wang, 2019). Hence, the weight assigned to Res is relatively higher. Deeper waters contribute to higher VO (10–20 cm/s), resulting in faster water renewal and more favorable nutrient conditions. These conditions affect the feeding physiology of oyster, leading to higher growth rates (Lee et al., 2017). Offshore wind and wave conditions, persistent wave action, and strong ocean currents can impact marine aquaculture facilities (Ogle et al., 1977).

Meteorological factors such as WS, Pr, and Tas are less commonly considered in suitability site selection studies for oyster marine aquaculture. These climate and meteorological factors do not directly impact marine aquaculture and carry a lower weight of 20%. However, they exhibit strong correlations with environmental factors. For example, Tas is related to SST, and the Res is caused by wind and waves, leading to changes in TSS concentration. Pr affects river flow rates, influencing factors such as SO and river nutrient inputs, which can impact Pacific oyster marine aquaculture.

Furthermore, meteorological factors directly reflect variations in these factors that can characterize climate changes and extreme weather events. Due to the lack of relevant data, this study considers fewer socioeconomic factors, weighing only 10%. Among the socioeconomic factors, city areas and piers with established infrastructure often carry a higher weight.

The quarterly variations in the spatial distribution of suitability for each sub-model in the Shandong offshore are shown in Figure 4 . The results indicated that water quality and hydrology suitability in the Shandong offshore are suitable, with scores of 6 or above within the range. Water quality suitability decreases in the summer, while hydrological conditions are best in winter and spring. Compared to other areas, the water quality and hydrological conditions in LZ were relatively poorer. Climate suitability in Shandong is relatively low, exhibiting higher suitability only near the coast, which may be attributed to the limited availability of distribution data for wild Pacific oyster in the Shandong offshore.

Averaging 140 suitability distribution maps from May 2011 to December 2022, we obtained the final spatial suitability distribution ( Figure 5A ). The suitability scores were categorized into three levels: high suitability (scores of 6–8), medium suitability (scores of 4, 5), and low suitability (scores of 1–3). The HSP was 48% in 20 km, with all areas except CD (19%) having higher HSP than the overall 20 km level, including LZ (100%), RC (92%), RS (72%), and JM (100%). High scores in LZ all being score of 6.

The suitability distribution in RC Sanggou Bay ( Figure 5B ) was compared with high-resolution satellite imagery and the offshore marine aquaculture database obtained by Liu et al (Liu et al., 2020b) ( Figure 5D ). This database is based on Landsat 8 remote sensing images and object-oriented NDWI and edge feature extraction (Wang et al., 2018), along with manual interpretation methods, providing the spatial distribution of marine aquaculture within a 100 km range of China’s offshore areas. The red areas in Figure 5C represented actual offshore marine aquaculture areas. In RC, the high suitability areas in spring, summer, and autumn matched the locations of rafts in this database, proving the accuracy of the evaluation model.

The monthly average change in the HSP ( Figure 6B ) showed that LZ (85.1%—100%), RC (87.16%—100%), and JM (80.39%—100%) maintained high suitability throughout the year (above 80%). The CD was more suitable in spring (March-May) but was lower than 20 km in other seasons. RS had lower suitability in summer (June-August), consistent with actual conditions, with minimal seasonal variation and HSP above 80% in other seasons. CD and RS exhibited more considerable monthly variation in suitability. Overall, the offshore areas within 20 km of Shandong had higher suitability in spring and autumn, with the annual HSP fluctuating around 50%. Spatially, the northern part of Shandong had higher suitability ( Figure 6A ).

The time series of HSP from May 2011 to December 2022 ( Figure 7A ) indicated that summer (June-August) was the most susceptible to climate events. Comparing the spatial distribution of suitability in July of average years with El Niño and La Niña years ( Figure 7B ), the results showed a general decrease in suitability during El Niño years and an increase during La Niña years. LZ and CD were less affected by climate events. In contrast, RC, RS, and JM were more susceptible to climate events in summer, with a noticeable decrease in suitability (scores of 7) during El Niño years, and high HSP (scores 6–8) area expanded within the 20 km range during La Niña periods. The increase in suitability during La Niña could be related to positive anomalies of Chl-a forced by winds (Herrera-Cervantes et al., 2020).

Figure 8 showed that under the SSP1–2.6 scenario, overall suitability in the Shandong offshore gradually increases in all future periods compared to the historical period. Under the SSP5–8.5 scenario, suitability initially increases and then decreases. However, the HSP (scores of 6–8) remains higher than the historical average, with a sharp decline in the proportion of the score of 7.

Our research indicated that future climate change impacts would increase suitability in the northern offshore of Shandong. Importantly, high suitability had remained stable at 100% in LZ, primarily with a score of 6. Furthermore, a suitability score of 7 increased under the SSP1–2.6 scenario and maintained the highest in the near term under the SSP5–8.5 scenario, later being surpassed by RC and RS in 2041–2070. This stability and high suitability would provide some reassurance in the face of potential climate change impacts.

The suitability evaluation model for Pacific oyster marine aquaculture develop in this study is applicable for assessing operations in the Shandong offshore area and can be adapted to suit the specific conditions of other regions as well. However, it is essential to note that there is still scope for further development of this model. Due to data limitations, this study considered a limited number of socioeconomic factors. When conducting suitability assessments, it is essential to take into account areas where marine aquaculture is prohibited or conflicts with existing activities, such as protected areas, net fishing, touristic traffic, commercial traffic, etc (Brigolin et al., 2017; Barillé et al., 2020; Porporato et al., 2020). Additionally, it is necessary to consider the infrastructure that can be integrated with marine aquaculture activities, such as offshore wind farms (Buck et al., 2008; Benassai et al., 2014; Di Tullio et al., 2018) and oil drilling platforms (Ogle et al., 1977), which has been piloted and studied worldwide and may serve as potential areas for Pacific oyster marine aquaculture.

The MaxEnt model used in this study has demonstrated superior predictive accuracy compared to other species’ geographic distribution models (Phillips et al., 2006; Wang et al., 2007). It has shown promising results even dealing with limited species distribution data (Hernandez et al., 2006), addressing the challenge of predicting species distributions with small sample sizes in marine environments. The application of MaxEnt in marine studies has witnessed rapid development (Hu et al., 2021; Hughes and King, 2024; Yang et al., 2024; Yati et al., 2024). In this study, due to the scarcity of wild Pacific oyster distribution data in the Shandong offshore area, we utilized available data ( Supplementary Table 1 ) from the entire northwest Pacific region (117°E-146°E, 30°N-45°N) for prediction. However, this approach may lead to an expanded study scope and imbalanced species distribution data, potentially affecting the accuracy of the predictions (Soley-Guardia et al., 2024). Obtaining more species distribution data in future research will enhance the accuracy of MaxEnt predictions. Additionally, utilizing surrogate species with easily accessible data, can be employed to predict suitable habitats for species with insufficient survey data (Barata et al., 2024).

Pearson correlation analysis of the time series data for each factor ( Figure 9 ) explored their interactions and impact on suitability, categorizing the Pearson correlation coefficients into five levels: very strong (0.8–1.0), strong (0.6–0.8), moderate (0.4–0.6), weak (0.2–0.4), and very weak (0–0.2), and analyzed the correlation between HSP, ONI, and MEI with time series factor.

According to Pearson’s correlation analysis, HSP is moderately positively correlated with Chl-a, TSS, and pH, weakly negatively correlated with WW, Res, WS, and SO, weakly positively correlated with Tas, and very weakly correlated with SST. The result indicated that suitability was mainly affected by Chl-a, TSS, and pH.

Although the analysis showed no direct correlation between HSP and climate indices, research indicated that climate events could impact suitability through factors like SST, wind speed (Liu et al., 2014), Pr, and Chl-a (Liu et al., 2014; Liu et al., 2020a). In this study, the MEI index was weakly negatively correlated with SO and Pr, only weakly negatively correlated with Pr in LZ. The ONI index was weakly negatively correlated with SO and showed a weaker negative correlation with Pr. The MEI index had a stronger correlation with each factor than the ONI index. It could serve as a better standard for determining whether suitability was affected by climate events.

Based on the above correlation analysis, the ENSO could cause more extreme changes in Pr. ENSO typically matures in winter and effects on subsequent summer precipitation in eastern China (Liu, 2021; Liu et al., 2024c). The difference between the suitability of Shandong marine aquaculture in El Niño and La Niña years ( Figure 6B ) reflects the asymmetry of their effects on precipitation (Liu et al., 2024c). Compared with La Niña, El Niño had a more noticeable impact on precipitation in Shandong (Guo et al., 2017).

Pr was weakly or more associated with all water quality factors except SO. The study showed that the increase in precipitation would increase the dissolved inorganic nitrogen and inorganic phosphorus of semi-enclosed bay seawater through the way of atmospheric settlement and runoff input, thereby causing an increase in Chl-a (Han et al., 2023; Xiao et al., 2024) and TSS (Meng et al., 2022; Ma et al., 2024), the Pr negatively associated with Chl-a and TSS concentration, which may result from more summer precipitation diluting Chl-a and TSS in seawater.

Additionally, the increase in freshwater flow caused by precipitation may lead to a rapid decrease in seawater salinity in the short term (Swam et al., 2022); precipitation changes may further lead to other factors affecting the suitability. For example, the correlation between Pr and water quality factors in LZ was small, and the suitability was mainly affected by the Chl-a concentration. The Chl-a concentration remained high and changed less in a year, which was less affected by ENSO and had high stability.

According to the correlation analysis, Figure 10 showed a boxplot analysis of time series factors and the HSP The graph visually showed the distribution of outliers of time series factors; among them, SST, DO, Res, and pH, Tas had no outliers. SO, WW, VO, and WS had fewer outliers and differences across regions, and ENSO was less likely to cause outlier changes by larger spatiotemporal scales. Chl-a, TSS, and Pr had more outliers in all regions, and they strongly correlate with the HSP and the climate index. This suggested that the appearance of outliers may be due to ENSO, Chl-a had many low values, TSS had many high values, Pr had many high values, and HSP declining more, which was consistent with the above correlation.

El Niño was classified based on the Oceanic Niño Index (ONI) values: weak (+0.5°C to +0.9°C), moderate (+1.0°C to +1.8°C), and high (greater than +1.8°C), represented by shades of red from light to deep. La Niña classifications were similarly based on ONI values: weak (-0.5°C to -0.9°C), moderate (-1.0°C to -1.8°C), and high (less than -1.8°C), represented by shades of blue from light to deep in Figure 11 .

Combining the correlation analysis and the boxplot analysis, we calculated the difference of Chl-a, TSS, Pr, and HSP between the monthly average from May 2011 to December 2022 to contrast the occurrence of outliers with the occurrence ENSO. We observed that the abnormal increase in Pr generally occurred in the summer of the second year following the La Niña event (Wang et al., 2000) in 2011, 2012, 2017, 2020, and 2022. However, Pr did not increase significantly during El Niño events. During the development period of El Niño, Pr decreased in East Asia, whereas Pr increased during the decline period of El Niño (Cao et al., 2017).

Chl-a consistently showed low values in the summer, It had declined in the summer following La Niña events in 2011, 2012, 2018, and 2021, consistent with Pr being affected by climate events. In addition, El Niño also had an impact on Chl-a. During the strong El Niño year of 2015, RS and JM showed low values, with an overall decline in December 2018. The reduction in Chl-a concentration during that period may be due to the resuspension caused by wind and waves (Liu and Wang, 2019) rather than precipitation; During the strong El Niño of 2015, East Asian summer winds weakened (Liren et al., 1997), result in weak resuspension and thus low Chl-a concentration, these two factor have weakly positively correlation.

TSS exhibited outliers during the autumn and winter seasons of the ENSO events, for example, an unusual increase in LZ in 2011 and in different areas in 2016, 2020, and 2022. The mature phase of ENSO typically occurs during the northern winter and is accompanied by the weaker-than-average winter winds along the East Asian coast (Liren et al., 1997). Consequently, the weakened resuspension did not result in an increased TSS concentration, and the Pr did not decrease significantly during this time. Although there was no correlation between the SST and the climate index in this study, SST in the Shandong offshore was strongly negatively associated with TSS and had a moderate negative correlation with Chl-a. The results indicated that the changes of SST will also affect the concentration of Chl-a and TSS, but the specific mechanism still requires further study. We found a sizeable abnormal decrease in Chl-a concentration in the summer of 2014, which led to the reduction of HSP, this phenomenon may be due to data processing and still needs further research.

Due to the favorable water quality and hydrological conditions, the Shandong offshore maintained a high suitability in approximately 50% of the region throughout the year. LZ, RC, and JM exhibit exceptionally high suitability (over 80%), making them ideal for developing Pacific oyster marine aquaculture. It was recommended to increase the aquaculture scale in high-suitability areas (LZ, RC, and JM) while reducing it in low-suitability areas (CD, RS) to avoid issues such as the declining sustainability of existing marine aquaculture environments (Forrest et al., 2009; Wang et al., 2018; Gao et al., 2020; Brito et al., 2023; Huang et al., 2023). Simultaneously, it is crucial to closely monitoring factors related to suitability, such as Chl-a, TSS, and pH. Furthermore, based on the results of a suitability evaluation and monthly variations, it was possible to relocate aquaculture operations across multiple marine regions. For instance, Pacific oyster marine aquaculture could be conducted in CD during spring, while RS can be avoided during summer.

The existing RC, RS, and JM marine aquaculture areas are susceptible to ENSO. Therefore, during the summer following El Niño and La Niña events, it is necessary to monitor abnormal weather conditions, especially abnormal rainfall, and develop appropriate response measures, such as reducing the density and scale of Pacific oyster marine aquaculture. Moreover, in JM and RC, which have high suitability but are susceptible to climate events, Pacific oyster could be integrated into an Integrated Multi-Trophic Aquaculture (IMTA) system (Shpigel and Blaylock, 1991; Jiang et al., 2013). This approach could enhance the carrying capacity and stability of the marine aquaculture ecosystem while also adapting to climate change (Ahmed and Glaser, 2016) and ensuring system sustainability (Khanjani et al., 2022).

Future climate change may increase the suitability in open sea, providing opportunities for the Pacific oyster marine aquaculture. Specifically, the spring in CD exhibits higher suitability, with hydrological conditions similar to those in winter ( Figure 4B ), which may contribute to higher growth rates (Lee et al., 2021). In the future, marine aquaculture in open sea areas with increased suitability could help mitigate the impacts of climate change on bivalve aquaculture systems (Walker et al., 2006; Froehlich et al., 2018). One of the challenges open sea areas faces is their susceptibility to meteorological conditions, which can impact marine aquaculture infrastructure. Pacific oyster’s marine aquaculture can be integrated with future offshore structures and platforms, such as offshore wind farms (Buck et al., 2008; Benassai et al., 2014; Di Tullio et al., 2018) and oil drilling platforms (Ogle et al., 1977). These developments have created new opportunities for Pacific oyster marine aquaculture. Some research is underway to explore integrating physical and biological conditions for the site selection of these co-located systems.

Furthermore, from 2070 to 2100 under the SSP5–8.5 scenario, except for LZ, the suitability for aquaculture significantly decreased, indicating that suitability did not exhibit a continuous upward trend and necessitates ongoing assessment. LZ maintained a stable suitability throughout the ENSO and future climate change scenarios, making it suitable for long-term marine aquaculture of Pacific oyster.