The study area, including the Bering Sea, Bering Strait, and Chukchi Sea, is characterized by a broad, shallow, and continental shelf region (De Robertis et al., 2017). This region is the only pathway connecting the Pacific Arctic region and the northern Pacific Ocean that provides nutrients and heat as water currents (Woodgate et al., 2005; Stabeno et al., 2018). The features of the water masses that flow into the Chukchi Sea vary markedly according to the season (Jinping et al., 2006; Norcross et al., 2010; Eisner et al., 2013; Gong and Pickart, 2015). Although Pacific winter water, including cold water and relatively salty water, such as newly ventilated Pacific winter water (NVPWW), dominates in the winter and early spring, Pacific summer water has a warmer and fresher water inflow into the Chukchi Sea through the Bering Strait in the summer (Woodgate et al., 2005; Pickart et al., 2019). Two water masses flow into the Chukchi Sea through the Bering Strait in summer: Alaskan coastal water (ACW) and Bering Summer Water (BSW) (Lin et al., 2019). ACW consists of the warmest and freshest water flowing into the Chukchi Sea along the eastern coastline of the Bering Strait (Pickart et al., 2019). BSW is relatively colder and saltier than the ACW flow into the Chukchi Sea through the central and western coastlines of the Bering Strait (Lin et al., 2019; Pickart et al., 2019). These water masses discharge around the Barrow Canyon, transporting different pathways: ACW flows into the eastern Chukchi Sea, whereas BSW flows through the western and central Chukchi Sea (Gong and Pickart, 2015). Furthermore, the Siberian Coastal Current (SCC) has freshwater flow into the western Chukchi Sea along the Siberian coastline (Weingartner et al., 1999; Spall, 2007). Although the seasonal cycle of salinity in the Chukchi Sea is controlled by Pacific-originated water such as BSW, the SCC also often influences the western Chukchi Sea (Weingartner et al., 1999; Spall, 2007). Atlantic water (AW) also flows into the Pacific Arctic region along with the shelf break of the Eurasian basin (Karcher and Oberhuber, 2002; Li et al., 2022). In the northern Chukchi Sea, the AW occasionally upwells along the shelf-break region (Bourke and Paquette, 1976; Gong and Pickart, 2015). Different water masses flow into the Chukchi Sea, with northward AW and Pacific Water from the south, and eastward-flowing SCC.

A multidisciplinary survey around the Pacific Arctic and sub-Arctic was carried out on the Korean icebreaking research vessel (IBRV) Araon from 26 July to 26 August 2020. Ichthyoplankton were collected using a ring net (mesh size, 505 μm; mouth area, 1.77m²) at 18 sampling stations (Figure 1). The net was obliquely hauled with a speed of 60 m min⁻¹, reaching a maximum depth of 344 m (Supplementary Table 1). In the shallow stations below 200-m depth, the net was towed from 5 m above the bottom to the surface. Tow speed and duration were about 1.5–2 knots and 20–30 min, respectively. The net depth was determined in association with the depth of sound scattering layers collected by a scientific echosounder (EK60). The collected samples were immediately subsampled using a Folsom plankton splitter, fixed in 4% filtered seawater formalin, and transported to the laboratory. Specimens were identified on the basis of morphological characteristics, such as habitus, preanus length, and melanophores, or ecological characteristics, such as habitats and spawning periods. Systematic accounts and morphology followed the methods by Matarese et al. (1989) and Okiyama (1989). The abundance was determined in individual numbers per cubic meter (ind./m³), based on the volume filtered by the net, which was measured using the revolution count of a flow meter attached to the center of the net mouth. Total length of specimens was measured with a caliper.

Vertical profiles of seawater temperature (Temp), salinity (Sal), dissolved oxygen (DO), and fluorescence (Flu) were obtained using a conductivity-temperature-depth sampler (Sea-Bird, SBE 911 plus) with an oxygen sensor at each sampling station. For the following analyses, the surface, average, and bottom values were used for Temp, Sal, DO, and Flu. Classifications between surface, average, and bottom values are represented by subscripts (e.g., Temp_surface, Temp_average, and Temp_bottom). The “surface” values were determined from the 5-m to 10-m depth, the “average” value from the surface to the maximum hauled depth at each station, and the “bottom” value from the maximum hauled depth to 10 m above it at each station. The water column of the present study area was divided into five water masses, AW, ACW, BSW, meltwater/runoff (MW), and winter water (WW), according to Gong and Pickart (2015) and Pickart et al. (2019).

All multivariate analyses were carried out using the PRIMER v7 package (Clarke et al., 2014; Clarke and Gorley, 2015) and PERMANOVA+ for PRIMER (Anderson et al., 2008). For the environmental analyses, draftsman plots were used on the basis of environmental data to apply an appropriate transformation (Valesini et al., 2014). As a result, Sal_surface, Sal_average, Sal_bottom, DO_surface, DO_average, and DO_bottom were transformed by log (X+1) prior to analyses. All environmental values were normalized before the analyses. The spatial environmental status of the 18 sampling stations was described using CLUSTER based on Euclidean distance from the log-transformed and normalized environmental values. PERMANOVA was performed to confirm statistically significant (P < 0.05) differences between separated groups. Biological-environmental (BIOENV) analysis was applied to verify which environmental factors best explained the distinction between the ichthyoplankton communities.

For the biotic analyses, a square root transformation was applied to the data before the analyses. A resemblance matrix based on Bray–Curtis similarity was calculated to identify dissimilarities between sampling stations. Hierarchical cluster analyses (CLUSTER) were conducted on the basis of the Bray–Curtis similarities using group-average linking. Similarity profile (SIMPROF) permutation test was simultaneously conducted with CLUSTER to confirm whether the divided groups represent statistically significant (P < 0.05) differences. PERMANOVA was carried out to confirm whether ichthyoplankton communities had a significant (P < 0.05) distinction according to the divided groups. Similarity percentage (SIMPER) analyses were conducted to verify which species contributed to the differences between groups. Canonical analysis of principal coordinates (CAP) was carried out to visualize patterns in divided groups. The correlations between species and environmental values and stations were plotted using CAP.

The ranges of the environmental factors among the sampling stations are summarized in Table 1; Figure 2. The highest temperature was 10.88°C at the surface of station A11. The water temperature tended to increase southward. The salinity ranged from 25.14 to 34.80 PSU. The highest salinity was observed at station B43, whereas the lowest salinity was observed at station B86. The reason for the large salinity gap between close stations was thought to be that the salinity values varied greatly according to the depth in this area. Additionally, the surface and average salinity decreased northward, but the bottom salinity showed a relatively even distribution pattern except for station 49. Although DO values were missing at stations B01, B03, B05, and B14, DO values were distinguished according to location. That is, the DO in the northern stations (>70°N) was much higher than the DO in the southern stations (<70°N), except for station B14. DO showed relatively low values in B14. Fluorescence indicated relatively high values at the surface of southern stations (A01, A05, A06, A07, A11, A15, and A18).

When the environmental values were analyzed using hierarchical cluster analysis, the stations were divided into five groups at three distance levels that were statistically significant (SIMPROF, P < 0.05, n = 999) (Figure 3A). At the 5.86 distance level, the stations were divided into two groups. One group included B20, B43, B49, B58, B70, B81, and B86. These stations were in the northern part of the survey region (>71°N). On the other hand, another group included southern stations (<71°N), A01, A05, A06, A07, A11, A15, A18, B01, B03, B05, and B14. In the first group, including northern stations, B49 separated from other stations, showing quite a high distance level. In the second group, including the southern station, B14 separated from other stations, representing quite a high distance level. The other stations were subdivided into two subgroups. One group included B01, B03, and B05, whereas another group is composed of A01, A05, A06, A07, A11, A15, and A18. As a result, the five groups could be distinguished from each other (Figure 3A). In the PERMANOVA, the null hypothesis that there was no difference in the environmental values between the five groups was rejected at a significance level of 5% (pseudo-F = 29.426, P < 0.05) (Table 2). However, in the pairwise tests, there were significant differences at a level of 5% only among three groups: B, D, and E (Table 2). Consequently, the environmental values could be distinguished according to the geographical distribution pattern, namely, north and south. In the southern stations, stations were subdivided into Bering Strait and Bering Sea.

Comparing the composition between the five water masses, AW, WW, ACW, BSW, and MW, according to the abovementioned groups, group A consisted of only MW. Group B was composed of subequal amounts of WW and MW, whereas some ACW and BSW were confirmed (Table 3; Figure 4). In group C, AW and WW dominated. The Pacific-originated water masses, ACW and BSW, dominated in groups D and E. Comparing water masses with the groups, AW, WW, and MW represented high correlations with group B, namely, northern stations, whereas ACW and BSW showed a positive correlation with groups C, D, and E, the southern stations (Figures 3B, 4B). Additionally, ACW only was confirmed at the B03 and B14. In the nMDS plot (stress, 0.04), the stations were closely grouped according to the cluster analysis grouping (Figure 3B). Temperature, salinity, and fluorescence showed a relatively high correlation with groups D and E. DO values positively correlated with group B. Groups A and C indicated negative correlation with the temperature and salinity.

The number of species and abundance is summarized in Table 4. In this study, 20 ichthyoplankton taxa were confirmed. Among the sorted specimens, the mean total length was 35.55 mm ranging from 6.11 mm (Hippoglossoides dubius) to 109.74 mm (Chauliodus macouni). The minimum abundance was indicated at station A07, whereas the maximum was represented at station B05. Gadus chalcogrammus was the most dominant taxon in the present study, followed by Pleuronectes quadrituberculatus and Boreogadus saida. Gadus chalcogrammus showed a relatively high abundance at station B05. Pleuronectes quadrituberculatus indicated high abundance at station B01. Boreogadus saida showed high abundance at station B43. Geographically, G. chalcogrammus and P. quadrituberculatus dominated near the Bering Strait, whereas B. saida dominated around the Chukchi Sea (Figure 5). These three species accounted for approximately 50% of the total abundance in the present study.

Through hierarchical cluster analysis, we divided the ichthyoplankton community into four groups (Figure 6). Group A included only station B14. Group B consisted of two stations (B58 and B70). Group C included five stations (B20, B43, B49, B81, and B86). Of the stations in group C, station 49 indicated relatively low similarity with the other stations. Group D was composed of ten stations (A01, A05, A06, A07, A11, A15, A18, B01, B03, and B05). The stations in group D showed a relatively low latitude distribution (< 70°N) (Figure 6B). The four groups were separated from each other at a similarity level of 0%. This result meant that the four groups had different species compositions. In the PERMANOVA, the null hypothesis that there were no significant differences in the ichthyoplankton community between the separated groups was rejected (pseudo-F = 3.6917, P < 0.05) (Table 5). However, pairwise comparisons represented no significant difference between the A & B, A & C, and A & D pairs of groups (P > 0.05), whereas the B & C, B & D, and C & D pairs of groups showed significant differences (P < 0.05) (Table 5). In the SIMPER analysis, the average similarity of group A was not determined because this group was composed of only one species (Eleginus gracilis) (Table 6). The average similarity of group B was 95.62%, showing that one species (Benthosema glaciale) contributed to the group similarity of 100%. Group C represented 58.11% of the average similarity with the 100% contribution of B. saida. This result meant that B. saida was the main species in group C. The average similarity of group D was 10.84%, and two species (Stenobrachius leucopsarus and G. chalcogrammus) contributed to the group similarity of the upper 70%. Considering the average similarity of group D, the community composition of group D was relatively more varied than that of the other groups. Comparing the group dissimilarity between the groups, the average dissimilarity between group A and the other groups was 100%, and E. gracilis commonly contributed to the group dissimilarity (Table 7). Additionally, E. gracilis was confirmed only in group A at station B14 (Table 4; Figure 5). The average dissimilarity of groups B and C was 100%, and two species (B. saida and B. glaciale) contributed to the group dissimilarity of the upper 70%. Among the two species, B. saida was confirmed only in group C, whereas B. glaciale was confirmed in only group B. The average dissimilarity of groups B and D was 100%, and five species (B. glaciale, S. leucopsarus, G. chalcogrammus, P. quadrituberculatus, and Sebastes reedi) contributed to the group dissimilarity of the upper 70%. Only B. glaciale was confirmed in group B, and the remaining species were confirmed in group D. Between groups C and D, the average dissimilarity was 100%. Seven species (B. saida, G. chalcogrammus, S. leucopsarus, P. quadrituberculatus, S. reedi, Gymnocanthus tricuspis, and Stomiiformes sp.) contributed to the group dissimilarity of the upper 70%. Among the seven species, B. saida and G. tricuspis were observed only in group C, whereas other species were confirmed only in group D.

In the CAP analysis, the first two squared canonical correlations clearly discriminated the four groups that were separated by CLUSTER analysis (Figure 7). The first squared canonical axis with a large correlation (δ² = 0.9949) separated ichthyoplankton communities in group B from other community groups. The second squared canonical axis with a large correlation (δ² = 0.9949) discriminated the ichthyoplankton communities of the groups, although groups A and D were not clearly distinguished from each other. Among the species that contributed to group similarity and dissimilarity, E. gracilis, G. chalcogrammus, P. quadrituberculatus, S. reedi, S. leucopsarus, and Stomiiformes sp. showed a positive correlation with groups A and D. Two species (B. saida and G. tricuspis) indicated a positive relationship with group C. Benthosema glaciale represented a high positive correlation with group B (Figure 7A). Comparing the relationship between ichthyoplankton groups and environmental factors, temperature, salinity, and fluorescence had a positive relationship with groups A and D, although Sal_bottom indicated a negative correlation with groups A and D. Sal_bottom and DO showed a positive relationship with group C. Group B did not indicate any strong relationship with the environmental values (Figure 7B). In the BIOENV analysis, the ichthyoplankton community was best described by a combination of Temp_surface, Temp_average. Sal_surface, Flu_average, and Flu_bottom (ρ = 0.498, P = 0.01) (Table 8).

For the environmental analyses, sampling stations could be divided into five groups. Geographically, groups A, B, and C were in the Chukchi Sea, whereas groups D and E were in the Bering Sea and Bering Strait, respectively. In other words, the former groups were included in the Pacific Arctic region, whereas the latter groups were in the Pacific sub-Arctic region. Group B included most of the Pacific Arctic region stations. Therefore, group B was considered to represent the summer physical properties of the Pacific Arctic region. However, group A was more closely grouped with group B than other groups, whereas group C was closer with groups D and E. Given that Pacific-originated water, ACW and BSW, have been confirmed in B14, the southern Chukchi Sea was influenced by the Pacific Ocean in this period. On the other hand, groups D and E were regarded as the typical Pacific sub-Arctic by including the stations located on the Bering Sea and Bering Strait. This clustering result corresponded well with the recent studies indicating that the Pacific Arctic and sub-Arctic biogeographic boundary now appears to the north of the Bering Strait (Baker et al., 2020, Baker et al., 2022). However, groups D and E can be subdivided into two groups according to the geographical region. Group D was composed of the Bering Sea station, whereas group E consisted of Bering Strait stations. In the nMDS analysis, the Pacific Arctic region stations, group B showed a positive correlation with the DO values. The Pacific sub-Arctic region stations, groups D and E, indicated a positive relationship with the water temperature, salinity, and fluorescence values. These results were derived from the effect of the warm surface water flow into the Chukchi Sea from the Pacific Ocean in summer (Jinping et al., 2006; Hu et al., 2019). This result well agreed with the future prediction that Arctic Ocean temperatures will increase derived from higher heat fluxes through the sub-arctic region (Drinkwater et al., 2021). The higher fluorescence in the Pacific sub-Arctic region than in the Pacific Arctic region is believed to be caused by nutritional supplementation from the Pacific Ocean (Clement Kinney et al., 2022). Comparing the water masses between groups, group B was represented by the Pacific Arctic region composed of various water masses. That is, AW, WW, and MW were mixed, although the proportion differed according to the stations. Groups D and E, located in the Pacific sub-Arctic region, consisted of BSW, although B03 and B14 also included ACW. Of the Pacific Arctic groups, although the water mass composition of group C was similar to that of the other Pacific Arctic groups, Pacific-originated water masses, ACW and BSW, were also observed in group C. Therefore, it is thought that Pacific-originated summer water masses expanded to group A in this period. On the other hand, another Pacific Arctic group, group A, was composed only of the MW. Considering that meltwater could inflow from the SCC in the summer season (Weingartner et al., 1999; Semiletov et al., 2005; Spall, 2007) and freshened water is extending further east (Clement Kinney et al., 2022), group A could be affected by freshwater likely SCC containing relatively low salinity. Taken together, in the summer season, the survey region could be distinguished into the Pacific Arctic and sub-Arctic regions although climate changes such as warming and decreasing sea ice are changing water properties (Hirawake et al., 2021). The temperature increased southward as the warm surface water expanded from the Pacific Ocean (Jinping et al., 2006; Gong and Pickart, 2015; Pickart et al., 2019). The salinity decreased northward on the surface, although the average and bottom salinity were similar to each other. Considering that sea ice melt provides freshwater, the northwardly low surface salinity seems to be related to sea ice melting around the Arctic region (Gong and Pickart, 2015). DO was much higher in the Pacific Arctic region than in the Pacific sub-Arctic region. Given that winter water normally contains much oxygen supplied from the atmosphere, the high DO value in the Pacific Arctic region is related to the high winter water composition (Codispoti and Richards, 1971; Nishino et al., 2016). Considering water masses, although WW or MW dominated in the Pacific Arctic region, Pacific-originated summer water, such as ACW and BSW, flowed into the Pacific Arctic region in the southern Chukchi Sea (Pickart et al., 2019). The freshened MW that originated from the SCC flowed into the western Chukchi Sea from the East Siberian Sea (Semiletov et al., 2005; Spall, 2007).

In this study, the ichthyoplankton community was distinguished according to geographical region. That is, the ichthyoplankton community was divided into two major groups, groups C and D. Group C consisted of the Pacific Arctic stations. The species that most contributed to the group similarity in the SIMPER analysis was B. saida, which is a typical dominant species in the Pacific Arctic region (Eisner et al., 2013; De Robertis et al., 2017; Vestfals et al., 2019; Axler et al., 2023). On the other hand, group D was composed of Pacific sub-Arctic stations. The species that contributed to the group similarity were G. chalcogrammus and S. leucopsarus, which commonly dominate in the Pacific sub-Arctic of the northern Pacific region (Pearcy et al., 1979; Beamish et al., 1999; Moku et al., 2000; Eisner et al., 2013; Majewski et al., 2017; Axler et al., 2023). Consequently, the two groups were thought to be made up of the typical ichthyoplankton community structure in each region. Groups A and B were distinguishable from groups C and D, respectively. Eleginus gracilis was the only species in group A. Considering a previous study, this species dominates between 66.5°N and 69.5°N (De Robertis et al., 2017). However, their distribution was confirmed above 70°N in 2019 (Levine et al., 2023). In this study, the distribution of E. gracilis was also confirmed above 70°N. These results agreed well with De Robertis et al. (2017)’s argument that the distribution of E. gracilis is farther north than in earlier surveys (Vestfals et al., 2019; Marsh et al., 2020; Baker, 2021). Given that E. gracilis favors relatively high water temperatures, between 7°C and 9°C, and relatively low salinity, the inflow of relatively high temperature and low salinity water, such as the ACW, could influence the ichthyoplankton community of the Chukchi Sea (Gong and Pickart, 2015; De Robertis et al., 2017; Sigler et al., 2017; Lin et al., 2019; Pickart et al., 2019; Vestfals et al., 2019). This result agrees well with the present study that warm water and the ACW expanded to station B14 (Figures 2, 4B). Benthosema glaciale was the only species in group B. Considering that this species is mainly distributed in the Atlantic Arctic region, the observation of this species is unusual (Zhang et al., 2022). However, Zhang et al. (2022) confirmed the distribution of B. glaciale in the Chukchi Sea. These findings also correspond well with the present study. Considering that AW can flow into the Chukchi Sea through shelf-break upwelling along the continental slope (Lin et al., 2016; Corlett and Pickart, 2017; Jung et al., 2021; Wang and Danilov, 2022) and Atlantic species have been detected in the northern Chukchi Sea (Wildes et al., 2022), it is possible that this species was transported into the Chukchi Sea with shelf-break upwelling. The discovery of this species corresponds well with the route of the shelf-break jet (Corlett and Pickart, 2017; Zhang et al., 2022). Therefore, the invasion of Atlantic species is possible in the Chukchi Sea with AW upwelling. Additionally, although station 49 was classified within group C in the analysis of the ichthyoplankton community, it was separated from the other group C stations in the environmental analysis, indicating a relatively low similarity. Furthermore, the water mass composition of station 49 differed from that of the other Pacific Arctic region stations by consisting of only MW. In comparing the species composition of station 49, the three species, namely, Gymnelus hemifasciatus, Icelus spatula, and Liparis fabricii, were observed only at station 49 (Table 4). Of the three species, I. spatula and L. fabricii are mainly distributed along the coastline, showing that they favor fresh water (Tokranov and Orlov, 2005; Zorina and Chernova, 2022). Furthermore, although the two species show a wide distribution ranging from the Northern Pacific Ocean to the Atlantic Arctic region, considering that the spawning period of Icelus spatula occurs from late August in the Arctic region and after September in the sub-Arctic region, the larval specimens in this study likely originated from the eastern Siberian Sea (Tokranov and Orlov, 2005; Zorina and Chernova, 2022). Liparis fabricii is more abundant around the Atlantic Arctic region (Able, 1990; Walkusz et al., 2016; Smirnova et al., 2022). Currently, the record of this species has occasionally been reported in the East Siberian Sea and Chukchi Sea since Able (1990), whereas the abundance is low (Mecklenburg et al., 2007; Norcross et al., 2010; Smirnova et al., 2022). Considering that the SCC has freshwater inflow to the Chukchi Sea from the east, the two species at station 49 are thought to have invaded from the East Siberia Sea, whereas the MW expanded to the Chukchi Sea (Weingartner et al., 1999; Semiletov et al., 2005; Tokranov and Orlov, 2005; Spall, 2007).