This study was conducted over a period of 39 months, from April 2021 to July 2024, across four representative coastal harbors in South Korea: Yangpo (East Sea), Tongyeong (Korea Strait), Bieung (Yellow Sea), and Jeju (offshore Jeju Island). These sites were selected to capture regional variations in oceanographic conditions and harbor characteristics (Fig. 1). Yangpo (East Sea) represents an open-coast harbor setting on the Korean east coast, where coastal hydrography is strongly influenced by the East Korea Warm Current (a branch of the Tsushima Warm Current) and associated alongshore variability (Pak et al. 2023; Lee et al. 2024). Tongyeong (Korea Strait) is located in a relatively sheltered, semi-enclosed southern-coast embayment with extensive harbor infrastructure and high vessel activity; hydrography in this region is strongly affected by Kuroshio-derived warm currents entering through the Korea Strait and seasonally varying stratification (Jiang et al. 2023; Kang et al. 2024). Bieung (Yellow Sea) represents a macrotidal, well-mixed west-coast environment where strong barotropic tides and tidal mixing can increase turbidity and drive short-term variability in nearshore conditions relevant to settlement (Lin et al. 2021; Du et al. 2022). Jeju (offshore Jeju Island) reflects a more oceanic island setting influenced by warm, saline Tsushima/Jeju Warm Current pathways and energetic tidal currents around the Jeju Strait, which can create pronounced spatial variability in temperature and salinity (Cha and Moon 2020; Kim et al. 2022; Kim and Choi 2024).

At each site, artificial settlement substrates were deployed to quantify percent cover of Watersipora subtorquata as an integrated measure of colony establishment and growth on the plates. Substrates consisted of square acrylic plates (300 mm × 300 mm × 5 mm), a material commonly used in harbor infrastructure. Plates were grouped in sets of ten and suspended vertically using polypropylene ropes, with adjacent plates spaced 30 cm apart, spanning approximately 1–4 m depth along the water column (uppermost to lowermost plates). Plates were redeployed in the same vertical configuration after each survey and were not intentionally reassigned among depths between sampling events. Because plates were deployed as standardized vertical arrays rather than as fixed-depth experimental treatments, percent cover values were interpreted as an integrated response across this representative 1–4 m near-surface depth range, and depth was therefore not included as an independent predictor in the GAM analyses. Three deployment points were established per site. Points were selected in areas with high vessel berthing and traffic, and were positioned closer to the seaward side of each harbor (near the harbor mouth) rather than nearshore inner sections, to better represent conditions experienced by floating and fixed harbor infrastructure. These plates remained submerged for a standardized exposure period that was consistent across all sites; only plates with a minimum exposure duration of 6 months were included in the analysis. Plates were deployed at least 6 months prior to the first survey (i.e., before April 2021), and the monitoring period (April 2021–July 2024) represents repeated observations of already established assemblages on these plates. Therefore, all observations used in the analyses corresponded to plates with a minimum immersion time of 6 months. Surveys were conducted in April, July, and October in 2021; in April, June, August, and October in 2022; in April, July, and October in 2023; and in April and July in 2024. Within each survey cycle, all four sites were sampled within approximately 2 weeks of each other (not necessarily on the same day) to ensure temporal comparability among sites (Table 1). During each field survey, plates were retrieved to the surface and photographed out of the water on shore using a digital camera (Olympus Tough TG-5, Japan), after which they were redeployed. These images were later analyzed using ImageJ software (Schneider et al. 2012) to determine the percent cover of W. subtorquata calculated as the percentage of the plate area covered by the species.

Percent cover ( % ) = ( Number of pixels in the area occupied by the target species Number of pixels in the total area of the attachment plate ) × 100

Colonies of Watersipora subtorquata observed on settlement plates were identified in the field based on colony morphology and diagnostic characters. Identification followed the taxonomic keys and illustrated guides for bryozoans in Korean waters (Seo 2005; Seo 2010) and was cross-checked against relevant taxonomic literature. To strengthen taxonomic reliability, identifications were additionally reviewed through consultation and collaboration with a Ph.D.-level taxonomic specialist in the field.

Water temperature and salinity were measured at each survey site using a multiparameter water quality meter (Pro Plus, YSI, USA). Measurements were taken simultaneously with plate retrieval to ensure temporal consistency between biological and environmental data. Water temperature and salinity were measured at approximately 2 m depth adjacent to the plate arrays during each survey. Because these measurements were collected at discrete survey times rather than continuously, they represent contemporaneous conditions during each sampling event. Site-specific summaries of observed temperature and salinity across surveys are provided in Suppl. material 1, and temporal changes in water temperature and salinity at each harbor are shown in Suppl. material 4. Because these measurements were collected only at the time of each survey, they do not represent integrated environmental histories prior to sampling. Accordingly, potential lagged effects of temperature and salinity on percent cover were not evaluated in this study.

Outliers in plate percent cover were identified and corrected using the Interquartile Range (IQR) method. For each sampling event at each site, plate-level percent cover values were pooled across plates (within the same location and survey cycle) to compute the IQR and identify outliers. Following Tukey’s rule, values above the upper fence (Q3 + 1.5 × IQR) were winsorized (i.e., capped rather than removed) by replacing them with the upper fence value (Tukey 1977). Values ≤ 0 were set to 0.01 to allow model fitting under a log link. To improve transparency of the raw observations and the extent of zero values and outlier handling, we provide a site-specific scatter plot of plate-level W. subtorquata percent cover over the monitoring period (Suppl. material 2). Observations flagged as outliers using Tukey’s IQR rule within each harbor × survey month (values > Q3 + 1.5 × IQR) are indicated in the figure. Raw zero values are shown as 0 in Suppl. material 2; however, for Gamma GAM fitting under a log link, zeros were replaced with 0.01 as described above. For statistical modeling, plate-level percent cover values were summarized to a single mean percent cover for each harbor and survey date after outlier treatment. Thus, although outlier detection and correction were performed on plate-level observations, the response variable used in the GAM analyses was the survey-event mean percent cover rather than individual plate values.

Generalized Additive Models (GAMs) were fitted separately for each harbor to evaluate the relationships between environmental conditions and W. subtorquata percent cover. For each harbor and survey date, plate-level percent cover values were first summarized to a single mean value, and the corresponding water temperature and salinity measured during that survey were used as predictors. Therefore, the statistical unit in the GAM analyses was the harbor-by-survey event rather than individual plates. No lagged environmental covariates were included; therefore, the GAMs were designed to evaluate associations between percent cover and temperature/salinity measured contemporaneously at each survey event. Analyses were conducted in R using the mgcv package (R version 4.4.1) (Zuur et al. 2009; Wood 2017). Because summarized percent cover values were positive and right-skewed, models were fitted with a Gamma error distribution and a log link. When necessary for model fitting under the log link, zero values were replaced with 0.01. For additive effects, the site-specific model was: percent cover ~ s(water temperature, k = 10) + s(salinity, k = 10), where smooth terms were represented using thin-plate regression splines. To assess the temperature–salinity interaction, we additionally fitted a tensor-product smooth model for each harbor: percent cover ~ te(water temperature, salinity, k = 3). Smoothness parameters were optimized automatically in mgcv. Model performance was evaluated using adjusted R-squared and deviance explained (%), and the statistical significance of smooth terms was assessed using F-statistics and associated p-values (α = 0.05). Because each harbor was represented by repeated survey events over time, the GAM results are interpreted as survey-event-level associations rather than as independent plate-level responses.

Results were visualized using the ggplot2 package in R (Wickham 2016). Smooth term plots were used to illustrate the individual effects of environmental factors, and contour plots were generated to display temperature–salinity interaction patterns across continuous gradients. The supplementary materials include site-specific environmental summaries (Suppl. material 1), raw plate-level percent cover observations by survey month, with zero values and IQR-identified outliers highlighted (Suppl. material 2), temporal variation in mean percent cover (Suppl. material 3), and temporal variation in water temperature and salinity (Suppl. material 4).

Percent cover of W. subtorquata in Yangpo Harbor demonstrated marked seasonal and interannual variability, reflecting strong environmental modulation over the study period. The annual average percent cover ranged from 0.86% in 2021, peaking at 2.48% in 2022. It then decreased to 0.90% in 2023 before it slightly increased again to 1.76% in 2024. Percent cover was particularly elevated during spring seasons. April 2022 had the highest percent cover (7.87%), followed by April 2021 (1.83%) and July 2024 (2.64%), indicating that higher percent cover was observed during several spring sampling events, particularly under the contemporaneous conditions recorded during those surveys (Suppl. material 3: fig. S3A).

Environmental monitoring recorded mean water temperatures of 18.0 ± 5.0 °C, with seasonal extremes ranging from 8.1 °C in April 2023 to 25.9 °C in July 2021. Salinity showed a moderate fluctuation, averaging 31.7 ± 1.7 PSU (range, 28.9 PSU in July 2024 to 34.3 PSU in April 2022) (Suppl. material 1).

GAM modeling identified a statistically significant negative linear effect of water temperature on percent cover (EDF = 0.9731, F = 0.5750, p = 0.0318), suggesting that percent cover decreased modestly with increasing temperatures (Table 2, Fig. 2A). However, salinity exhibited no significant influence on percent cover (EDF = 2.4825, F = 19.1478, p = 0.6151) despite its nonlinear relationship shape (Fig. 2B).

Additionally, the interaction between temperature and salinity was not significant (EDF = 2.1458, F = 0.7076, p = 0.1065), indicating independent environmental effects in this harbor (Table 3, Fig. 3A).

Percent cover in Tongyeong Harbor followed a bimodal seasonal trend, showing elevated cover during both early spring and mid-summer. The highest monthly percent cover occurred in April 2021 (2.08%) and July 2023 (1.39%). However, annual averages showed a general decline from 1.09% in 2021 to 0.21% in 2024, suggesting reduced percent cover over time (Suppl. material 3: fig. S3B).

Water temperatures varied from 14.3 °C to 23.7 °C, with a mean of 19.4 ± 3.4 °C. Salinity ranged from 28.4 PSU to 34.0 PSU, with an average of 31.5 ± 1.9 PSU (Suppl. material 1).

GAM analysis (Table 2) revealed a significant nonlinear effect of temperature (EDF = 7.7853, F = 300.5334, p = 0.0426) on percent cover, exhibiting a U-shaped percent cover response with higher percent cover observed at temperatures below 16 °C or above 21 °C (Fig. 2C). In contrast, salinity exhibited no significant effect (EDF = 0.0000, F = 0.0000, p = 0.4223), indicating that variation in salinity did not influence percent cover patterns at this site (Fig. 2D).

Additionally, no statistically significant interaction was found between temperature and salinity (EDF = 0.2435, F = 0.0394, p = 0.3142), suggesting that the effect of each variable was independent in Tongyeong (Table 3, Fig. 3B). These findings suggest that Tongyeong’s W. subtorquata population is very sensitive to specific temperature windows but insensitive to salinity variations, possibly reflecting its local adaptation or physiological tolerance.

Percent cover at Bieung Harbor was the lowest among the study sites. Monthly percent cover seldom exceeded 0.84%, often remaining under 0.40%. Despite low averages, localized seasonal peaks were observed in April 2022 (0.84%) and October 2022 (0.75%), suggesting episodic increases in cover under specific conditions (Suppl. material 3: fig. S3C).

Water temperatures in Bieung ranged from 15.5 °C to 27.0 °C, with an average of 21.4 ± 3.2 °C across the monitoring period. Salinity in Bieung ranged from 29.7 PSU to 34.5 PSU with a mean of 31.2 ± 1.6 PSU (Suppl. material 1), and showed a slightly wider fluctuation than in other sites.

GAM results (Table 2) revealed that both water temperature and salinity had statistically significant effects on percent cover in Bieung Harbor. Water temperature showed a significant nonlinear effect (EDF = 3.5402, F = 4.2444, p = 0.0453), indicating that percent cover varied nonlinearly across the observed temperature range (Fig. 2E). Salinity also had a statistically significant positive effect on percent cover (EDF = 0.8088, F = 0.6952, p = 0.0251), suggesting that higher salinity levels in this semi-enclosed estuarine setting were associated with higher percent cover of W. subtorquata (Fig. 2F).

Additionally, no significant interaction between temperature and salinity was observed (EDF = 0.3654, F = 0.0946, p = 0.1635) (Table 3, Fig. 3C), indicating that both temperature and salinity contributed to variation in percent cover in Bieung, although their interaction was not significant. Although temperature and salinity showed significant additive associations with percent cover in Bieung, the explanatory power of models based only on these variables remained limited, suggesting that additional unmeasured factors may also contribute to temporal variation in percent cover at this site.

Jeju Harbor exhibited the greatest variability in percent cover and environmental conditions among the four study sites. Percent cover increased over the study period, with the highest monthly value reaching 1.45% in July 2024. Annual averages of percent cover increased from 0.21% in 2021 to 0.78% in 2024. However, because temperature and salinity were measured only at discrete survey times, this temporal increase should not be interpreted as direct evidence of ocean-warming effects (Suppl. material 3: fig. S3D).

Environmental conditions were notably dynamic, with water temperature fluctuating between 15.8 °C and 27.7 °C (mean: 21.1 ± 4.1 °C) and salinity ranging from 25.6 PSU to 36.6 PSU (mean: 31.2 ± 2.9 PSU). These broad ranges are characteristic of offshore sites exposed to variable oceanographic processes (Suppl. material 1).

GAM results (Tables 2, 3) showed that temperature had a statistically significant positive effect on percent cover (EDF = 0.7255, F = 0.5383, p = 0.0230), indicating that higher percent cover was associated with higher contemporaneous temperatures under the salinity conditions observed at this site (Fig. 2G–H). Although salinity alone was not significant (EDF = 0.5923, F = 0.2719, p = 0.0625), its interaction with temperature was significant (EDF = 1.5439, F = 0.2334, p = 0.0141), indicating a synergistic effect between salinity and temperature under specific environmental thresholds (Fig. 3D).

Percent cover was maximized when temperatures exceeded 25 °C and salinity remained within 25–27.5 PSU, revealing a narrow but optimal environmental window associated with higher percent cover of W. subtorquata in Jeju. These findings suggest that percent cover was highest under a relatively narrow range of contemporaneous temperature–salinity conditions in Jeju, although future expansion cannot be inferred directly from these snapshot measurements alone.

The influence of water temperature and salinity on percent cover varied by location. Results are summarized in Table 4. While temperature showed significant effects across all four harbors, salinity played a meaningful role only in Bieung Harbor. Notably, a statistically significant interaction between the two variables was only observed in Jeju Harbor.

Percent cover of Watersipora subtorquata differed among the four harbors (Friedman test across 12 survey cycles: χ² = 9.70, df = 3, p = 0.0210), although values were generally low at all sites. This variation appears to reflect differences in local temperature and salinity. Temperature showed a significant association with percent cover across sites. Percent cover slightly decreased with increasing temperature in Yangpo. In Tongyeong, there was a more complex response, with higher percent cover at both low and high temperatures. In Bieung, percent cover also varied significantly with temperature, showing a nonlinear response across the observed range. In Jeju, the relationship was statistically significant, although the trend was relatively stable. These associations were broadly consistent with previously reported temperature ranges under which W. subtorquata has been observed to perform well (Boyd et al. 2002). However, because temperature and salinity were measured only at the time of each survey, our data do not allow direct inference about growth or recruitment processes, and the observed relationships should be interpreted as associations with contemporaneous conditions rather than as evidence of integrated seasonal drivers. These associations should therefore not be interpreted as evidence of causal environmental drivers of percent cover variability. Salinity effects were less consistent among harbors. At most sites, no strong association was found. However, in Bieung, percent cover increased with salinity, suggesting that local salinity variability may help explain cover dynamics at this site. At the same time, these results should be interpreted cautiously, because temperature and salinity alone are unlikely to capture all of the drivers of percent cover variability. In this more enclosed setting, additional factors such as local hydrodynamics, turbidity, propagule supply, and substrate-related differences may also contribute to variation in percent cover. Similar patterns have been noted for other fouling organisms in estuarine systems (Dunstan and Johnson 2004). Interpretation of these site-specific patterns should also consider differences in local harbor settings that can influence propagule delivery and colony development. For example, open-coast versus semi-enclosed geomorphology can alter water residence time and the retention of larvae near plate surfaces, and differences in harbor infrastructure and vessel traffic may modify substrate availability and introduction pressure among sites.

A significant interaction between temperature and salinity was only observed in Jeju. In Jeju, the predicted temperature–salinity response surface indicated higher percent cover under warmer temperatures (e.g., > 25.0 °C) combined with relatively lower salinity (approximately 25.0–27.5 PSU) within the observed range (Fig. 3D), consistent with a significant interaction between the two variables. However, this interaction should be interpreted cautiously because environmental variables were measured as snapshot conditions and do not reflect continuous environmental exposure. The relatively open coastal structure and stable conditions around Jeju might have contributed to this pattern. These findings are in line with studies showing that certain invasive species respond strongly to specific combinations of environmental factors (Ricciardi and MacIsaac 2000). In contrast, other sites showed no significant interaction effects possibly due to reduced environmental variability or stronger physical isolation from open water. In Yangpo and Tongyeong, both of which are open-coast harbors, interaction signals may have been harder to detect due to the limited sampling frequency and the relatively narrow observed range of one or both variables during the survey windows. Differences in connectivity and nearshore larval delivery among harbors may also affect whether interaction signals are detectable, particularly if establishment is supply-limited or episodic at some sites. A further consideration is potential cryptic diversity within Watersipora. Mackie et al. (2012) reported latitudinal segregation among Watersipora lineages (“new sp.” and additional clades) and suggested that lineage occurrence can be correlated with temperature. Because our identifications were based on morphology, we cannot fully exclude the possibility that percent cover responses at some sites reflect contributions from a different lineage. Future work integrating morphometric and DNA-based identification would help resolve lineage-specific cover–environment relationships.

These results highlight the need for site-based monitoring strategies. Jeju appears particularly suitable for W. subtorquata establishment under certain conditions, and ongoing warming may be consistent with increased risk, although this cannot be directly inferred from the present data. Bieung also warrants attention given its sensitivity to salinity variability, as changes in estuarine inputs or circulation could alter establishment patterns. We did not use satellite-derived products because coastal harbor environments are strongly influenced by local circulation and freshwater inputs, and satellite pixels may not reliably represent nearshore conditions at the spatial scale of the deployment sites. A limitation of this study is that temperature and salinity were measured intermittently during field surveys rather than continuously, and thus may not fully capture short-term variability. In addition, because environmental covariates were measured only at the time of sampling, our analyses do not account for lagged effects or prior environmental exposure that may also influence percent cover. In addition, the present models considered only temperature and salinity, and therefore did not explicitly incorporate other potentially important drivers such as hydrodynamic conditions, turbidity, propagule pressure, or substrate-related variation. Another limitation is that plates spanned an approximately 1–4 m depth range and were analyzed as an integrated vertical array, so potential fine-scale depth effects within this near-surface window were not evaluated explicitly. A further limitation is that each harbor was represented by a relatively small number of repeated survey events over time and thus potential temporal autocorrelation among observations was not explicitly modeled in the present analyses. Future studies using high-frequency temperature–salinity loggers, denser temporal sampling, and plate-specific depth records would better resolve percent cover dynamics, lag effects, within-array depth variation, and temporal autocorrelation.