Patterns of summer ichthyoplankton distribution, including invasive species, in the Bering and Chukchi Seas

A multidisciplinary survey was carried out in the Paci ﬁ c Arctic and sub-Arctic regions of the North Paci ﬁ c Ocean on the Korean icebreaking research vessel Araon . During this survey, ichthyoplankton ﬁ shes in the Paci ﬁ c Arctic and sub-Arctic region ranged from the Bering Sea to the northern Chukchi Shelf in summer. The most dominant species was Gadus chalcogrammus , followed by Pleuronectes quadrituberculatus and Boreogadus saida . Gadus chalcogrammus and P. quadrituberculatus were particularly abundant near the Bering Sea and Bering Strait, whereas B. saida was dominant in the Chukchi Sea. Hierarchical cluster analysis revealed four distinct ichthyoplankton communities in Paci ﬁ c Arctic and sub-Arctic regions based on geographical regions. However, Eleginus gracilis , which was previously known to be seen between latitudes 66.5°N and 69.5°N, was found above 70°N, suggesting that its distribution extends further north. Furthermore, we noticed that Benthosema glaciale , which is usually found in the Atlantic sector of Arctic Ocean, was observed in the northern Chukchi Sea. In addition to these unusual species distributions, several species that are mainly observed in coastal areas are observed in the Chukchi Sea region. The observed in ﬂ ux of various uncommon ﬁ sh species into the Chukchi Sea can be attributed to multiple factors, including freshwater in ﬂ ow from the East Siberian Sea and the intrusion of warm Atlantic and Paci ﬁ c waters, which are strongly affected by global warming. Consequently, it is imperative to conduct rigorous monitoring of the Paci ﬁ c Arctic region, with a particular focus on the Chukchi Sea, to better understand the implications of global warming.


Introduction
The Pacific Arctic Ocean has been undergoing dramatic environmental change related to sea ice melting and increasing surface temperatures (Mueter and Litzow, 2008;Wang et al., 2018;Baker et al., 2020Baker et al., , 2023)).Heat advected into the Arctic region resulted in a promoted decrease in sea ice coverage in the Arctic Ocean (Polyakov et al., 2017;Richards et al., 2022;Wang and Danilov, 2022).As a result, the minimum sea ice boundary in the western Arctic Ocean has gradually retreated northward, and the mean sea ice cover has been consistently below the mean sea ice cover over the last decade (Wood et al., 2015).Although the Bering Sea was excepted from these trends, increased warm southerly winds and air temperature influenced on the sea ice retreat in the Bering Sea, recently (Stabeno and Bell, 2019;Baker et al., 2020).Furthermore, more extended sea ice free periods are predicted in the Chuck Sea and Bering Sea under the global warming (Wang et al., 2018;Baker et al., 2023).This environmental change could affect the entire Arctic marine ecosystem, such as the community, behavior and phenology, growth and condition, abundance and demography, and range and regime shift (Wassmann et al., 2011;Wisz et al., 2015;Levine et al., 2023).
Ichthyoplankton are essential in Arctic pelagic ecosystems for evaluating ecosystem functions, assessing stocks, and investigating recruitment processes (Kendall and Duker, 1998;Kendall, 2000).They are also important sentinels to understand the consequences of Arctic climate change, as they rely on planktonic prey and lead to phase shifts in ecosystem dynamics.For this reason, many studies have been conducted in different regions of the Arctic Ocean, focusing on aspects of Arctic ichthyoplankton ecology.Oceanographic, geographic, and atmospheric features influence ichthyoplankton ecology, such as sea ice conditions (Bouchard and Fortier, 2008), water masses (Atwood et al., 2010;Busby et al., 2014;David et al., 2016;Randall et al., 2019), bathymetry (Doyle et al., 2002), currents (Siddon et al., 2011), precipitation (Boeing and Duffy-Anderson, 2008), freshwater input (Bouchard and Fortier, 2011), wind (Rodriguez et al., 2011), and air temperature (Genner et al., 2010).Furthermore, shifts in the assemblage of species, including the movement of sub-Arctic species into the Arctic region have been noted in recent years (Baker and Hollowed, 2014;Vestfals et al., 2019;Huntington et al., 2020;Marsh et al., 2020;Baker, 2021;Wildes et al., 2022).
Despite the ecological significance of ichthyoplankton, our knowledge of Arctic ichthyoplankton is insufficient in the highlatitude Arctic Ocean because of the harsh sea ice conditions that limit access.Few studies have been conducted around the East Siberian Shelf in the Pacific Arctic region, whereas most studies on ichthyoplankton have focused on the Pacific sub-Arctic and Arctic waters of the Bering, Chukchi, and Beaufort Seas (Benoit et al., 2008;Eisner et al., 2013;Falardeau et al., 2014;Randall et al., 2019;Stevenson and Lauth, 2019;Cooper et al., 2023;Levine et al., 2023).
Our aim of this study was to describe the community structures of ichthyoplankton and examine the relationship with environmental conditions in the Bering Sea and Bering Strait, Chukchi Shelf, northern Chukchi Shelf, and northern East Siberian Shelf in austral summer in 2020.This result can provide further knowledge to understand the variation in the species assemblages of larvae and juvenile fishes, as well as the ecology of adult fish in the Pacific Arctic Ocean.

Study area background
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 shelfbreak 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.

Data sampling and processing
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 mm; mouth area, 1.77m 2 ) at 18 sampling stations (Figure 1).The net was obliquely hauled with a speed of 60 m min −1 , 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 3 ), 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).

Statistical analyses
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.Biologicalenvironmental (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.

Spatial variation in physical and biological properties
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 The spatial pattern of the environmental values among the sampling stations.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.

Spatial pattern of the ichthyoplankton community
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).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 relationship between the ichthyoplankton community and environmental variables
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 (d 2 = 0.9949) separated ichthyoplankton   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 (r = 0.498, P = 0.01) (Table 8).

Physical properties of the Pacific Arctic and sub-Arctic regions in summer
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   (Jinping et al., 2006;Hu et al., 2019).This result well agreed with the future prediction that Arctic Ocean Comparison of the dominant species abundance (ind./m 3 ) at each survey station.

B A FIGURE 6
Dendrogram (A) indicating separated groups and map (B) representing divided groups according to the CLUSTER analysis.
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).

Ichthyoplankton community structure
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   (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).

Conclusions
In the Pacific Arctic, mainly in the Chukchi Sea, B. saida was the most dominant species, which is a typical Arctic fish (Eisner et al., 2013;Marsh et al., 2020;Baker et al., 2022;Axler et al., 2023;Cooper et al., 2023;Levine et al., 2023).In the Pacific sub-Arctic region, including the Bering Sea and Bering Strait, G. chalcogrammus and P. quadrituberculatus were the dominant and subdominant species, respectively (Eisner et al., 2013;Vestfals et al., 2019;Axler et al., 2023).Therefore, it is thought that the Bering Strait represented a more similar community structure to the Bering Sea than the Chukchi Sea in the summer season.This conclusion is supported by the biotic analyses of the close clustering of the stations around the Bering Strait with the Bering Sea.Additionally, the species composition was also well matched to typical species in the Pacific sub-Arctic region.Taken together, the typical ichthyoplankton communities in each region were confirmed in this survey.However, unusual community structures were also observed at several stations.At station B14, E. gracilis was observed.Considering a previous study, E. gracilis is more distributed northward than in previous reports because this species is mainly distributed between 66.5°N and 69.5°N (Vestfals et al., 2019;Baker, 2021;Levine et al., 2023).According to the previous studies, they are moving northward as the water temperature increases due to global warming (De Robertis et al., 2017;Baker et al., 2020;Axler et al., 2023;Levine et al., 2023).Furthermore, B. glaciale, typically  distributed in the Atlantic Arctic region, was observed in the present study.Zhang et al. (2022) argued that this species is likely to invade the Pacific Arctic region through the expansion of warm AW because cold water restricts the expansion of B. glaciale to the North Pacific Arctic region (Zhang et al., 2022).This result is associated with "Atlantification, " which is triggered by the warm surface temperature and melting sea ice because sea ice melting leads to AW upwelling caused by the increased wind effect on the northern Chukchi Sea (Lin et al., 2016;Corlett and Pickart, 2017;Jung et al., 2021;Wang and Danilov, 2022).In addition to these unusual community compositions, several species, I. spatula and L. fabricii, moved to the Chukchi Sea region from the East Siberia Sea in this study as the SCC, including freshwater, expanded to the western Chukchi Sea.Considering that Siberian river discharge from permafrost thaw is increasing, water inflow from these rivers could often occur (Nummelin et al., 2016;Wang et al., 2021;Kanamori et al., 2023).All the abovementioned uncommon compositions are related to global warming, such as increased water temperature and sea ice melting (Richards et al., 2022;Wang and Danilov, 2022).That is, various fishes are transported into the Chukchi Sea as the water temperature increases and sea ice retreats.Species transportation from the Atlantic and East Siberian regions in this study is consistent with the findings of Wisz et al. (2015), who argued that transfers from the Atlantic to the Pacific Ocean would be primarily facilitated through the Northwest Passage.This phenomenon has been confirmed in various taxa (Beaugrand et al., 2002;Perry et al., 2005;Reid et al., 2007;Wisz et al., 2015;Kim et al., 2020, Kim et al., 2022).However, Gadus macrocephalus, another dominant species in the Pacific sub-Arctic regions, has not been observed in this survey although this species is also known as extending their distribution northward (Stevenson and Lauth, 2019;Baker, 2021;Cooper et al., 2023).Considering there have been various studies reported these days that record a dramatic shift in its geographical distribution, shift of the spatial distributions of the ichthyoplankton in the Pacific region (Mueter et al., 2017;Vestfals et al., 2019;Baker, 2021;Drinkwater et al., 2021).Further monitoring of the ichthyoplankton community in this region is required because of their commercial values (Stevenson and Lauth, 2019;Marsh et al., 2020;Baker, 2021).As the Chukchi Sea water masses are influenced by a variety of factors, including freshwater flow from the East Siberian Sea and warm water intrusion from the Pacific and Atlantic Oceans (Figure 8) (Grebmeier et al., 2006;Norcross et al., 2010), comprehensive monitoring of interbasin exchanges in the Chukchi Sea is critical in the context of global warming.

FIGURE 1
FIGURE 1 Map of the survey area showing the sampling stations in the Pacific Arctic and sub-Arctic regions from July to August 2020.
FIGURE 3Dendrogram (A) and nonparametric multidimensional scaling (nMDS) plot (B) of the environmental values.
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 FIGURE 4 Potential temperature-salinity scatter plot (T-S diagram) for the sampling stations (A) and relative composition (%) of water masses (B).ACW, Alaskan coastal water; AW, Atlantic water; BSW, Bering Sea water; MW, meltwater/runoff; NVWW, newly ventilated Pacific winter water; RWW, remnant Pacific winter water.

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

TABLE 1
Summary of the environmental values collected during the survey in the study area.

TABLE 2
Results of the PERMANOVA test based on the environmental values according to groups and between (pairwise tests) groups from each group.
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.

TABLE 3
Relative composition (%) of water masses of the divided groups.
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

TABLE 4
Species composition and abundance (ind./m 3 ) of the ichthyoplankton community in the Pacific Arctic and sub-Arctic regions.

TABLE 5
Results of the PERMANOVA test based on the ichthyoplankton abundance data according to groups and between (pairwise tests) groups.

TABLE 6
Result of the SIMPER analysis showing species that contribute to the group similarity of the divided groups.

TABLE 7
Result of the SIMPER analysis, representing species that contribute to the group dissimilarity between the divided groups.

TABLE 8
Summary of the results of the BIOENV analysis showing the best-matched combination between environmental values and the ichthyoplankton community..498Temp surface , Temp average , Sal surface , Flu average , Flu bottom 0.01 rBest combination of environmental factors P r, Spearman correlation coefficient; P, statistical significance level.Significance differences (P < 0.05) were bolded.