Food Web Trophic Structure at Marine Ranch Sites off the East Coast of Korea

Understanding the trophic ecology of the giant Pacific octopus Enteroctopus dofleini is challenging in developing marine ranches and in reestablishing its regional stocks against the severe stress of fishing. We adopted carbon and nitrogen stable isotope techniques (termed δ13C and δ15N, respectively) to identify the trophic niche (i.e., pathways and positions) of this species systematically in the entire food webs of two marine ranches off the east coast of the Korean peninsula. While a slight spatial shift in the isotopic nestedness of faunal communities was observed, the δ13C and δ15N values of consumers were distinct and separate among functional groups at both ranches. The consumer δ13C values spanned a broad range between pelagic and benthic sources of organic matter, and their δ15N values recorded a stepwise trophic-level enrichment, indicating that suspension feeders and herbivore-deposit feeders served as baselines of pelagic- and benthic-based trophic pathways, respectively. The δ13C values of predators, including E. dofleini, were arrayed between the two primary consumer groups. Neither δ13C nor δ15N values showed any remarkable variations with increasing octopus weight. Dietary mixing-model calculations indicated that E. dofleini is a generalist predator relying on both benthic- and pelagic-affinity prey, similar to some teleost species that consume a diverse spectrum of prey. In contrast, other teleost groups showed prevalent trophic links with either pelagic- or benthic-based pathways. The trophic-level estimations revealed that E. dofleini occupies an intermediate position slightly below the teleosts. A lack of discrete trophic positions between E. dofleini and teleosts seemed to be indicative of the released teleost predation but instead reflects the imposed food competition. Overall, the results demonstrated that despite compositional changes in the taxa constituting individual trophic groups, E. dofleini occupied a very similar trophic niche in both ranching systems. Finally, as extracted from information based on octopus marine ranches launched on natural rocky bottoms, our isotopic evidence provides a greater understanding of the trophic ecology of this octopus species in nearshore natural habitats along the southwestern margin of its distribution range.


INTRODUCTION
Roughly 300 octopus species of the 845 cephalopod species known worldwide have been described to date (Hoving et al., 2014;Jereb et al., 2014). Octopuses are known for their intelligence and for being masters of camouflage (Mather et al., 2010). They serve as generalist feeders (Hanlon and Messenger, 1996;Mather et al., 2010) in various oceanic regions. Octopuses are also an important fisheries resources, with an annual global production of 305,851 to 373,846 metric ton (t) between 2005 and 2014 (8-12% of the total cephalopod catch; Sauer et al., 2019). Because octopuses are animals of ecological and economic importance, scientists have long sought to better understand their biology and ecology (Hanlon and Messenger, 1996;Scheel et al., 2014Scheel et al., , 2016Caldwell et al., 2015).
The giant Pacific octopus Enteroctopus dofleini Wülker 1910 (hereafter referred to as GPO) is considered to be the largest octopus species, and is distributed in the North Pacific ocean, in coastal areas spanning the United States, Russia, Japan, and Korean peninsula (Scheel, 2002;Cosgrove and McDaniel, 2009). Seasonal onshore-offshore migratory patterns are also found in the eastern and western Pacific (Hartwick et al., 1984;Sato and Yorita, 1999;Rigby and Sakurai, 2005). This species has a short lifespan (3-5 year) and fast growth, occupying benthic habitats at its subadult stage (Hartwick and Barriga, 1997). The GPO often use dens and select habitats with a preference for the presence of boulders and dense kelp cover in shallow coastal waters (Scheel, 2002;Scheel and Bisson, 2012). As observed by intermittent mobile tracking devices, they are stationary and in hiding >90% of the time, and larger individuals have greater movement ranges (Scheel and Bisson, 2012). The feeding ecology of E. dofleini could be a key component in understanding ecological processes associated with its unique life history, habitat selection and use, and movement patterns. E. dofleini is considered to be a generalist predator that feeds on both benthic and pelagic prey items such as decapods, bivalves, gastropods, fish, cephalopods, and, occasionally, even seagulls and sharks (Mottet, 1975;Hartwick et al., 1981;Scheel and Anderson, 2012). The diet of E. dofleini varies depending on prey abundance, with the proportions of major prey items differing with habitat type (Vincent et al., 1998). Although the prey size generally increases as the GPO grows (Hartwick et al., 1981), E. dofleini individuals take a wide range of prey sizes, reflecting foraging behavior with no preysize selectivity (in individuals weighing >2.5 kg, for example; Vincent et al., 1998). Therefore, it can be summarized that the dietary composition of E. dofleini might reflect various intrinsic and extrinsic factors in different habitats (Vincent et al., 1998;Greenwell et al., 2019).
Although E. dofleini can avoid predation by remaining in a protected den, camouflaging itself, or hiding in kelp beds (Hartwick et al., 1978), juveniles and adults are preyed on by a variety of predators, such as large fish, marine mammals, and other octopuses (Gillespie et al., 1998). GPO also exhibits interspecific competition with sea otters (Vincent et al., 1998). In addition, E. dofleini is fished commercially and is consumed as a source of protein by humans around the Pacific (Takeda, 1990;Gillespie et al., 1998;Lee et al., 2014;Sauer et al., 2019).
Despite important roles as major foraging predators with prey preferences and avoidances but no individual specialization (Marsteller, 2018) and as prey in benthic marine communities, and the possible impact of fishing on ecosystem functioning (Sauer et al., 2019), systematic studies on the GPO trophic niche (trophic pathways and positions) in coastal marine food webs are still lacking compared with squids (Hobson and Cherel, 2006;Merten et al., 2017;Murphy et al., 2020).
The identification of prey items in the diet of GPO has commonly employed analyses of prey contents in the gastric tract (Mottet, 1975;Smith, 2003;Ibáñez and Chong, 2008;Marsteller, 2018) and of hard-shell prey remains around their dens (i.e., midden piles; Hartwick et al., 1978;Scheel and Anderson, 2012). However, the validity of these approaches has been challenged (Michener and Kaufman, 2007). Stomach contents exhibit prey compositions that differ markedly from midden counts for shelled mollusks and soft-bodied organisms, and flaws in prey identification have been identified for both methods (Smith, 2003). Rapid digestion and external ingestion can make the stomach contents visually unidentifiable (e.g., Bo et al., 2019). Midden contents left around dens after the GPO processes prey are often highly fragmented, making accurate identification difficult (Boyle and Rodhouse, 2005). The GPO may move considerable distances on foraging paths from their dens, so the midden contents might not reflect their complete diets (Hartwick et al., 1984;Scheel and Bisson, 2012). More importantly, their diet compositions may seldom capture topological descriptions of food webs to characterize their trophic niche . Alternatively, carbon and nitrogen stable isotope analyses (referred to here as δ 13 C and δ 15 N values, respectively) have been applied successfully to the identification of diets, feeding habits, and trophic ecology of cephalopods (Cherel et al., 2009;Merten et al., 2017;Murphy et al., 2020). The dietary resources of a species within a marine food web can be resolved by the δ 13 C values of its tissues, as the latter contains information on the source (benthic or pelagic) of prey items (Kang et al., 2008;Park et al., 2020). Similary, their δ 15 N values provide an estimation of their position in the food web, increasing with the assimilation of resources at relatively higher trophic levels (TLs; Post, 2002;Caut et al., 2009). Therefore, stable isotope-based metrics can be used to characterize the trophic niche of an organism as well as the structure, function, and dynamics of food webs (Layman et al., 2007;Cucherousset and Villéger, 2015).
Artificial reefs have been commonly deployed to increase fishery productivity around the world (Seaman, 2019). The installation of artificial reef structures is being implemented to improve habitat quality and thereby support marine ranching for E. dofleini off the east coast of the Korean peninsula (see later Study Site). Despite its known relationships between resources and consumers, improving our understanding of its trophic niche within the marine food web will need a definition of ecological role and further management of the marine ranching systems. In this study, we aimed to delineate food web structures in association with the trophic ecology of E. dofleini in the nearshore marine ranching ecosystems off the east coast of the Korean peninsula. To do this, we examined δ 13 C and δ 15 N values of the GPO and other sympatric species collected from artificial reefs deployed to create nursery and spawning grounds to determine whether its trophic interaction varies among habitats where community compositions differ from one another. Given the documented foraging behavior of GPO on available resources (Greenwell et al., 2019), it is expected that δ 13 C and δ 15 N values of E. dofleini would reflect trophic pathways supporting its nutrition and trophic position within the marine food web. We tested this hypothesis by comparing the isotope composition of E. dofleini and different community components in two octopus marine ranches. Our ultimate goal was to identify the generality of its trophic niche and highlight its ecological role (as predator, prey, and trophic competitor) in the nearshore communities.

Study Site
The annual catch of E. dofleini has declined from 5,500 to 3,700 t since the late 1990s on the shallow shelf of the East Sea of Korea and exploitation of its wild population has been quite intense for small individuals (less than 10-15 kg, 50% group-maturity weight, Lee et al., 2014). The Republic of Korea government has developed marine ranching to reestablish the stock by deployment of artificial reefs specified as a suitable habitat for breeding and nursery ground of GPO. The installation of these marine ranches are acknowledged as a pilot program to enhance the sustainable productivity of the octopus fisheries by habitat improvement, nursing, and resource management. Once site selection by evaluation of the habitation and migration characteristics, and the following ranch creation are completed, any type of exploitation (fishing, recreational diving, etc.) is prohibited within the ranches. Their relevant activities include deploying artificial reefs, constructing kelp beds, and releasing juveniles of GPO. Further information of food web structure (i.e., ecological and trophic interaction) in the marine ranch ecosystems is needed to develop effective management strategies to maintain the habitats and enhance the GPO stocks. This study was conducted at two marine ranching sites (Mukho coast-MH, approximately 37 • 32 N, 129 • 8 E; and Pohang coast-PH, 36 • 0 N, 129 • 35 E) in the shallow subtidal zone off the east coast of the Korean peninsula (Figure 1). The MH site is often covered by the North Korea Cold Current and the PH site by the Tsushima Warm Current. In total, 3,142 and 330 reef blocks covering 110 ha and 122 ha at the MH and PH sites, respectively, have been deployed at a water depth of 10-40 m to attract GPO since 2018. Artificial reefs consist of truncated square concrete pyramids with holes and round pots into which GPO can crawl. Given the expected heterogeneity of community composition occupying benthic and pelagic habitats, we chose these two study sites to test the generality of the trophic niche of this species. Only small patches of kelp occurred at the two marine ranches.

Field Sampling and Laboratory Processing
Samplings for organic matter sources and consumers were conducted simultaneously at the MH site during March-April 2019 and at the PH site during April-May 2019. On each sampling occasion, water samples for suspended particulate organic matter (SPOM) were collected at each site using a 20-L van Dorn water sampler from 1 m below the water surface and prefiltered, onboard, through a 180-µm Nitex screen, to remove zooplankton and any large particles. Water sample processing continued immediately (within 1 h of collection) after transport to the wet laboratory, located on the coast. About 10 L of seawater for measurements of stable C and N isotopes in the SPOM were filtered through precombusted (450 • C for 2 h) Whatman GF/F filters (diameter, 25 mm; pore size 0.7 µm). After filtration, filters for the C isotope ratios were acid-fumed overnight in an HClcontaining desiccator to remove inorganic carbonates but filters for the N isotope ratios were not acidified. Then all the filters were folded, wrapped in aluminum foil, and then kept in a vacuum desiccator until analysis.
For isotopic measurements, samples for epilithon were collected by SCUBA divers. The divers collected 5-6 small-sized rocks randomly and brought to the diving boat. Then epilithon were sampled on board by scraping the rock surface gently using a toothbrush. Materials recovered and washed with ultrapure water were prefiltered through a 100-µm mesh net to remove large particles, and transported to the laboratory. Immediately after this, the water samples containing epilithon were filtered through precombusted (450 • C for 2 h) Whatman GF/F filters and treated following the procedure used for SPOM samples. Macroalgae were sampled by SCUBA divers and washed with filtered seawater on board. After being transported to the laboratory, macroalgae were rinsed with Milli-Q water, sorted to identify them, and air dried. Stable isotope values previously reported for epilithon and macroalgae at the MH site were used in the present study (Kang et al., 2008;Park et al., 2015).
Macrobenthic invertebrates for isotopic measurements were collected by SCUBA divers by scraping the rocky substratum using a stainless steel knife. Specimens of microbenthic invertebrates were placed in iceboxes filled with filtered seawater and immediately transported to the laboratory. They were sorted live under a dissecting microscope and kept alive overnight in filtered seawater from the sampling site to evacuate their gut contents. Only live and intact organisms were dissected to collect tissues.
Nektonic and benthonic organisms (including fish, GPO, crabs, and shrimps) for isotope measurements were captured using a gill net (50 m long × 2.5 m high, 9 cm mesh) and pots (60 cm long × 30 cm diameter, 1 cm mesh net) with frozen mackerel as bait. These fishing devices were anchored on the rocky bottom for 12 h (overnight). A few GPOs were also captured by the divers. All the GPO and nektonic samples collected were retained in iceboxes filled with ice during the transportation to the laboratory. In the laboratory, invertebrates (including GPO) were dissected to collect substantial quantities of tissues for isotope analysis. For fish, muscle tissues from the dorsal region were collected after the removal of their viscera by dissection and defatted in a solution of methanol, chloroform, and water (2:1:0.8; Bligh and Dyer, 1959) to avoid a bias in the δ 13 C value by the effect of the contents of isotopically lighter ( 13 C-depleted) lipids compared with other biological macromolecules (i.e., carbohydrates, proteins, and nucleic acids, Focken and Becker, 1998). Invertebrate samples were not defatted because of the low contents of lipids. Given the negligible amount of carbonates, animal tissue samples were not acid treated. They were rinsed with Milli-Q water and kept frozen in the laboratory until subsequent treatment.
Epilithon, macroalgae, invertebrate (including octopuses), and fish samples were lyophilized and pulverized to a fine powder with a ball mill (Retsch MM200 Mixer Mill, Hyland Scientific, Stanwood, WA, United States) and then stored in a vacuum desiccator until isotope analysis.

Stable Isotope Analysis
Aliquots of the powdered samples (0.5-1.0 mg for animals; 1.5-2.0 mg for algae) were placed into tin combustion cups (8 × 5 mm, Elemental Microanalysis, Devon, United Kingdom) and filters containing SPOM samples were wrapped in tin disks (Thermo Scientific, Bremen, Germany). The sealed samples were introduced into an elemental analyzer (Eurovector 3000 Series, Milan, Italy) connected to a continuous-flow isotoperatio mass spectrometer (CF-IRMS; Isoprime, GV Instruments Ltd., Manchester, United Kingdom) to determine their δ 13 C and δ 15 N values. The stable isotope data were expressed as the δ notation: i.e., deviations from conventional standard reference materials (Vienna Pee Dee Belemnite for carbon and atmospheric N 2 for nitrogen) following the formula: δX ( ) = [(R sample /R standard )−1] × 10 3 , where X is 13 C or 15 N and R is 13 C/ 12 C or 15 N/ 14 N. International standards for sucrose [ANU C 12 H 22 O 11 ; National Institute of Standards and Technology (NIST), Gaithersburg, MD, United States] and ammonium sulfate ([NH 4 ] 2 SO 4 ; NIST) were used to calibrate the internal reference material (urea) and the instrument. One measurement of this laboratory reference material was repeated after every five samples to ensure analytical precision. The average difference among samples was <0.15 for δ 13 C and <0.2 for δ 15 N based on replicate measurements.

Data Analysis
Prior to performing the statistical analyses of data, the assumptions of normality and homogeneity of variance were assessed using the Shapiro-Wilk procedure and Levene's test, respectively. We tested for significant differences in δ 13 C and δ 15 N values between pelagic and benthic producers as sources of organic matter to identify the base of trophic pathways. Because the δ 13 C or δ 15 N values of producers did not meet the assumption of homogeneity of variance in both sites, a permutational multivariate analysis of variance (PERMANOVA, based on Euclidean distance) was employed to test the significance of variabilities between sites and sources. When a significance was detected, the two isotope values were then analyzed separately using a Mann-Whitney nonparametric U test or Student's t test. After measuring the body weight (as a proxy of size) and stable isotope values of GPO individuals, we employed a regression analysis to examine the relationships between weight vs. δ 13 C and δ 15 N values to analyze their size-based dynamics in food sources and TLs.
To identify patterns of trophic pathways, comparisons of the δ 13 C and δ 15 N values of consumers were made to test the levels of isotopic overlap between communities in the two sites and the isotopic niches of trophic groups within a site. First, the isotopic niches between consumer groups of both sites were compared by their overlap (position and area) in the isotopic space. For this comparison, isotopic nestedness (INes) was calculated after scaling to 0-1 range of stable isotope values (Cucherousset and Villéger, 2015). The calculation of INes is based on a convex hull area that is a measure of the isotopic richness of each sampling site (Layman et al., 2007). The INes value represents the ratio between the niche volume overlapped and the minimal volume filled by a group on the two-dimensional 0-1-scaled space (δ 13 C and δ 15 N). Isotopic overlap metrics are measured using the convex hull volume (represented by the areas of the two tested groups) and the volume of isotopic space they shared. The index varies from 0 (no overlap between the trophic niches of the tested groups) to 1 (a complete overlap or the group with the lowest volume is located within the larger one). Second, a hierarchical cluster analysis (Bray-Curtis similarity, average grouping method) was carried out on the mean δ 13 C and δ 15 N values of each species to determine trophic groups of consumers and depict the overall configuration of the trophic web structure of these artificial reef communities (Davenport and Bax, 2002;Kang et al., 2008). A further graphical representation of trophic groups was produced using a δ 13 C-δ 15 N bi-plot. The trophic groups classified (grazer, detritivore, herbivore, predator, suspension feeder, deposit feeder, omnivore, and scavenger) were identified by previously published guilds (Ruppert et al., 2004;Hong et al., 2006; http://www.fishbase.org/search.php; http://eol. org; http://www.marinespecies.org). Because the average isotope values of trophic groups met the assumptions of equality of covariance by Box's test, a multivariate analysis of variance (MANOVA) was then run on two main factors (sites and trophic groups) to test significant differences among trophic groups. After the MANOVA was performed, a one-way analysis of variance (ANOVA) followed by a Tukey honestly significant different (HSD) multiple comparison test was conducted for each dependent variable (δ 13 C and δ 15 N) to test for differences among trophic groups. For all the tests, significance was set at P = 0.05. PERMANOVA was performed using PRIMER version 6 combined with PERMANOVA + PRIMER add-on (Anderson et al., 2008), INes was estimated using the R code published by Cucherousset and Villéger (2015), and the other analyses were performed using IBM SPSS Statistics software (version 21.0; IBM Corp., Armonk, NY, United States).
We estimated the relative contribution of potential prey to the nutrition of secondary and higher-level consumers (i.e., predators) using a Bayesian isotopic mixing model (available as free Stable Isotope Analysis software in R: versions SIAR 4.2 and R 4.0.2; www.r-project.org). Benthic-vs. pelagic-affinity primary consumer groups were identified as prey items by the abovementioned cluster analysis to identify major trophic pathways. The SIAR mixing model allows the inclusion of sources of uncertainty (Parnell et al., 2010). The variability in the isotope values (mean and standard deviation, SD) of potential prey species, isotope values of each individual of target species, and trophic enrichment factors (TEFs, isotopic shifts between predators and their prey) should be incorporated into the model. We employed TEF values of 1.3 ± 0.3 for δ 13 C and 3.3 ± 0.26 for δ 15 N as the average fractionation factors (McCutchan et al., 2003). Given the uncertainties of the parameters and the within-individual variability, the probability distribution of prey contributions to each target predator group was expressed as 50, 75, and 95% confidence intervals.
Finally, the TLs of consumers were estimated using the following formula (Vander Zanden and Rasmussen, 1999;Post, 2002): where TL i and δ 15 N i represent the trophic level and the mean δ 15 N value of species i, respectively; δ 15 N PC represents the mean δ 15 N values of primary consumers obtained at each site in the present study [i.e., herbivores Acmaea pallida = 5.3 (MH site) and Strongylocentrotus nudus = 7.8 (PH site); suspension feeders Mytilus galloprovincialis = 6.2 (MH site); and Balanus sp. = 7.7 (PH site), respectively]; λ PC is their TL (=2 in this study); and TEF is the trophic enrichment factor (=3.3 ± 0.26 ).

Stable Isotope Ratios of Primary Sources of Organic Matter
In the present study, SPOM, epilithon, and macroalgae were considered the most possible sources of organic matter. SPOM was sampled as a representative of pelagic source at both sites. Benthic producers (i.e., epilithon and macroalgae) were collected only at the PH site and their isotopic values at the MH site were compiled from the previous studies (Kang et al., 2008;Park et al., 2015). Mean δ 13 C and δ 15 N values of SPOM, epilithon, and seven macroalgal taxa were presented in Table 1. The mean δ 13 C values of primary producers displayed a broad range between −22.2 (SPOM) and −14.8 (Grateloupia elliptica), measured at the MH site. Their mean δ 15 N values fell within a relatively narrow range between 2.8 (Sargassum horneri, MH site) and 7.2 (epilithon, PH site). The isotope composition of producers differed significantly between sampling zones (as defined by pelagic and benthic zones; PERMANOVA, Pseudo-F = 19.945, P < 0.001) but not between sites (Pseudo-F = 1.964, P = 0.160, Table 2). The isotopic difference between pelagic and benthic producers within each of the sampling sites was generated by their δ 13 C values (pairwise Mann-Whitney U test or Student's t-test, P < 0.025), with overall mean values of −22.0 ± 0.2 and −20.9 ± 0.3 (pelagic) vs. −16.8 ± 1.8 and −16.5 ± 1.4 (benthic) at the MH and PH sites, respectively. No significant differences in mean δ 15 N values were found at both sites (P > 0.5).

Isotopic Compositions of Consumers and INes Between Communities
We collected 21 and 22 faunal taxa including E. dofleini for stable isotope analysis at the MH and PH sites, respectively (Tables 3, 4). The mean δ 13 C values of consumers exhibited relatively broad ranges of −21.9 ± 0.3 (Styela clava) to −15.8 ± 0.1 (Hemicentrotus pulcherrimus) and −20.0 ± 0.4 (Balanus sp.) to −15.3 (Platycephalus indicus) at the MH and PH sites, respectively, spanning the overall δ 13 C range of primary producers with the exception of two gastropods (−13.3 ± 0.6 , Chlorostoma turbinata and −13.6 ± 0.3 , Omphalius pfeifferi pfeifferi). The mean δ 15 N values of consumer taxa ranged from 5.3 ± 0.4 (A. pallida) to 11.1 (Sebastes schlegelii) and 7.5 ± 1.0 (Strongylocentrotus nudus) to 12.6 (Platycephalus indicus) at the MH and PH sites, respectively. The δ 15 N values of predators generally increased by 2 to 6 compared with those of primary consumers (i.e., grazers, herbivores, suspension feeders, and deposit feeders) at both sites. Based on the δ 13 Cδ 15 N space configuration of consumer isotope values (Figure 2), the INes value between communities of both sites was calculated as 0.46, highlighting a considerable overlap but a slight spatial shift in isotopic niche.
GPO individuals collected in this study fell within a narrow range of 0.24-3.92 kg in total weight: small compared with 0.34-20.93 kg in the north of the study area (Lee et al., 2014). No significant relationships between isotope values and total weight (i.e., proxy of size) of GPO were found (r 2 = -0.091and 0.118; ANOVA, F 1,11 = 0.000 and 2.599, P = 0.988 and 0.135 for δ 13 C and δ 15 N values, respectively, Figure 3). Their δ 13 C values were slightly less enriched at the MH site (mean −18.0 ± 0.5 ) than that (−16.4 ± 0.3 ) at the PH site (Student's t-test, t 11 = −5.782, P < 0.001; Tables 3, 4). In contrast, no significant difference in their mean δ 15 N values was observed between sites (t 11 = −0.344, P = 0.738), with 9.7 ± 1.0 and 9.9 ± 0.5 at the MH and PH sites, respectively.

Identification of Trophic Groups
The cluster analysis based on mean δ 13 C and δ 15 N values divided consumers into three main groups, A, B, and C, and two subgroups within group C (CI and CII; Figure 4). Some species were not assigned to any of the three groups (e.g., S. schlegelii, Rapana venosa, Kelletia lischkei at the MH site; Chlorostoma turbinate and Omphalius pfeifferi pfeiffer at the PH site). The trophic groups clustered in general agreement with their feeding guilds, postulated from the literature, rather than according to their taxonomic relationship. Groupings of consumers exhibited similar configurations between sites on the δ-space but it should be noted that groups CI and CII shifted in their position between MH and PH sites (Figures 5A,B). Group A included suspension feeders such as mussels, barnacles, and tunicates. The sea star Asterias amurensis and the longfin yellowtail Seriola rivoliana were also included in this group.  located close to this subgroup. The remaining fish species were assigned to subgroup CII. At the PH site, E. dofleini was also included in subgroup CI, together with mostly fish species. The predatory crustacean Ovalipes punctatus, gastropod K. lischkei, fish species S. pachycephalus, P. indicus, and Hexagrammos otakii were assigned to subgroup CII. Multivariate analysis of variance revealed that the isotopic signatures of consumers differed significantly among trophic groups identified by cluster analysis (Pillai's Trace: F 1,32 = 90.578, P < 0.001) as well as between sites (Pillai's Trace: F 3,32 = 67.013, P < 0.001), generating a significant interaction (site × trophic group) effect (Pillai's Trace: F 3,32 = 8.031, P < 0.001; Table 5). One-way ANOVA revealed that this was due to differences in both δ 13 C and δ 15 N values (F 7,32 = 69.83 and 38.77, P < 0.001 for both). A subsequent Tukey HSD test revealed a significant difference in mean δ 13 C values between cluster groups A and B constituting primary consumers at both sites (P < 0.05). Their δ 13 C values fell within the ranges for SPOM and benthic algae, respectively (we designate these as pelagicaffinity and benthic-affinity groups hereafter). Group A had slightly higher δ 15 N values than those of group B at each site and both groups had consistently higher δ 15 N values at the MH than at the PH sites. The Tukey HSD test indicated that the mean δ 13 C values of cluster groups CI and CII differed from those of primary consumer groups A and B at the MH site (P < 0.05), intervening between them. The mean δ 13 C values of both groups tended to overlap with benthic-affinity group B at the PH site (P > 0.05). Mean δ 15 N values of group C generally increased by 2 to 6 compared with those of primary consumer groups (i.e., suspension feeders and herbivores) at individual sites.

Mixing Model and Trophic Level Estimations
Given the similarity in the δ 13 C values of cluster groups A and B to pelagic and benthic sources of organic matter, respectively, a SIAR mixing-model calculation was performed using the δ 13 C and δ 15 N values of the pelagic-and benthic-affinity prey as the isotopic baselines of trophic pathways. Our mixing-model estimations revealed that benthic-affinity prey (group B) made a higher median contribution (54 and 60%; range of 51-57% and 56-63%; and for MH and PH sites, respectively) to the nutrition of consumer group CI than that of pelagic-affinity prey (group A; median contribution 46 and 40%; range of 43-49% and 37-44%; for MH and PH sites, respectively; Figure 6). The dietary contribution of pelagic-and benthic-affinity prey to consumer group CII was estimated at 84% (82-90%) vs. 16% (10-18%), respectively, at the MH site. In contrast, the relative contribution of both prey items to group CII was estimated to be 30% (27-35%) vs. 70% (65-73%), respectively, at the PH site.
The calculation of trophic levels indicated that consumer TLs ranged from 2.0 to 3.5 at both sites (Tables 3, 4). Fish were at the highest positions of 3.3-3.5 and 2.6-3.5 characterized by predators at the MH and PH sites. Predatory invertebrates (gastropods and crustaceans) occupied the intermediate trophic levels of 2.7-3.1 between primary consumers (suspension feeders, herbivores, and deposit feeders) and fish species. The only exception was Asterias amurensis, a predatory species whose TL (2.1) was similar to primary consumers. TLs of primary consumers were 2.0-2.5 and 2.0-2.8 at the MH and PH sites, respectively. E. dofleini overlapped their trophic level with other predatory invertebrates at the MH site but slightly higher than that of primary consumers at the PH site. Finally, the overall combination of the mixing-model and TL estimations characterize the food web architecture in these marine ranching ecosystems (Figure 7).

DISCUSSION
The annual catch of E. dofleini has declined from 5,500 to 3,700 t since the late 1990s on the shallow shelf of the East Sea of Korea and exploitation of its wild population has been quite intense for small individuals (less than 10-15 kg, 50% group-maturity weight, Lee et al., 2014). The Republic of Korea government has developed marine ranching to reestablish the stock by deployment of artificial reefs specified as a suitable habitat for breeding and nursery ground of octopuses. This study investigated food web structures in association with the trophic niche of E. dofleini at the artificially-created habitats of marine ranching areas for the GPO in the southwestern margin of its distributional range in waters off the Korean peninsula. There was a considerable overlap in isotopic niche of overall communities between sites. The stable isotope analysis of E. dofleini and its related floral and faunal community revealed that the GPO has a trophic niche as a generalist feeder  relying on both benthic-and pelagic-affinity prey in the two nearshore marine-ranching habitats. Although based on limited sampling sizes at one time of year, our results gave preliminary insight into carbon flow pathways through food web and trophic interactions between E. dofleini and its sympatric species in the marine ranching ecosystem, as conceptualized by food-web models (Figure 7).

Food Web Structure of Marine Ranching Ecosystem
Despite a considerable overlap in the isotopic niches of overall communities between sites as indicated by INes (0.46), species at the PH site seemed to have higher scaled isotopic values than those at the MH site (Figure 2). This result indicates a spatial heterogeneity in the contribution of particular food sources possessing higher isotope values than those of other producers to dietary sources of consumers. Two main groups constituting basal food sources included SPOM and benthic algae in both food webs. Their δ 15 N values did not differ but the δ 13 C values were lower in SPOM than in benthic algal sources. Benthic algal species displayed a wide δ 13 C range.
The shift in isotopic composition of primary consumers probably reflect their dependence on and the availability of the particular macroalgal species and SPOM within their habitats (Fredriksen, 2003;Kang et al., 2008). Of benthic producers,   Mean δ 13 C and δ 15 N values (±1 SD) of individual cluster groups at two sites. Bolds are significant at P < 0.001. Means followed by the same superscript letter in δ 13 C and δ 15 N values are not significantly different (P < 0.01, Tukey HSD test).
epilithon had higher δ 13 C values at the PH (mean −15.6 ) than at the MH (−18.8 ) site. δ 13 C values of SPOM (as a proxy of phytoplankton) was also slightly higher at the PH (−20.9 ) than at the MH (−22.0 ). The current-and windinduced upwelling creates a prevailing physical setting along the PH site of the Ulleung Basin and facilitates nitrate-fueled primary productivity (Hyun et al., 2009;Ji et al., 2019). In these nitratereplete upwelling conditions, the phytoplankton community is dominated by diatoms (Hyun et al., 2009), which generally exhibit high δ 13 C values compared with other phytoplankton groups (Nadon and Himmelman, 2006;Kang et al., 2009). Slight increases in the δ 15 N values of SPOM, epillithon, and macroalgae were found between sites, probably reflecting the assimilation of 15 N-enriched nitrates from upwelled water at the PH site. Spatial variations in the δ 13 C and δ 15 N values of both producers and consumers caused by varying nutrient availability (oligotrophic vs. eutrophic) conditions has been observed along the eastern and southern coasts of the Korean peninsula (e.g., Kang et al., 2007).
FIGURE 6 | Relative contributions (50, 75, and 95% confidence intervals as indicated by black, dark gray, and light gray colors, respectively) of benthic-vs. pelagic-affinity primary consumer groups to the nutrition of higher-level consumers at the MH (A) and PH (B) sites, which were estimated using the Bayesian isotopic mixing model. Prey groups A and B represent suspension feeder and herbivore-deposit feeder groups, respectively, as identified by the hierarchical cluster analysis. Groups CI and CII indicate the subgroups clustered for predatory taxa. It may be inferred that an increased nutrient availability might affect isotopic values of producers and consumers (Kang et al., 2008(Kang et al., , 2009). The broad δ 13 C ranges of consumers at both sites-spanning the values of primary sources of organic matter-suggest a relatively complex food web in these shallow marine ranch ecosystems, in which different origins of organic matter, from pelagic to benthic producers, were used by faunal members. Differentiation of the δ 13 C values of macrozoobenthic species was consistent with their feeding guilds reported in the literature (Ruppert et al., 2004; http://eol.org; http://www.marinespecies. org), revealing that consumers of different feeding guilds exploit different nutritional sources. More negative δ 13 C values for suspension feeders (cluster group A) than those for herbivores and detritivores (cluster group B) indicate that while the former feeds exclusively on SPOM, the latter group is dependent primarily upon benthic algae-derived organic matter. In addition, high δ 13 C values in detritivores may reflect the microbial degradation and colonization of sedimentary organic matter (Dunton et al., 2012). The use of different food sources according to feeding zone (pelagic vs. benthic) of primary consumers has been well recognized in coastal ecosystems (Deegan and Garritt, 1997;Kang et al., 2008). Taxa belonging to both suspension feeder and herbivore-deposit feeder groups, which are located at the base of the trophic web as indicated by their low δ 15 N values, differed considerably between sites. Compared to δ 15 N values of suspension feeders, the δ 15 N value of A. amurensis, which is known to prey primarily on bivalves (Ross et al., 2002), at the MH site probably indicate its omnivorous feeding strategy as a deposit feeder. The inclusion of the longfin yellowtail S. rivoliana, known as an opportunistic piscivorous predator preying mainly upon juvenile pelagic fish (Bareiros et al., 2003; http://www.fishbase.se), into suspension feeder group A thanks to its low δ 15 N value was notable. Minagawa and Wada (1986) found extremely low δ 15 N values (−2.1 to +1.0 ) of the N 2 fixing blue-green algae in the offshore oligotrophic area of the East China Sea compared to 2 to 4 higher δ 15 N values of mixed phytoplankton in the northern coastal sea, which provided a systematic enrichment of δ 15 N values in mixed plankton, zooplankton, and fish along the food chain. The low δ 15 N values of plankton and fish in the southern offshore region may explain a transient migration of the longfin yellowtail from southern offshore to the nearshore reef of the PH site. This fish species is known to be widely distributed in tropical and subtropical waters of the Indo-Pacidic and Atlantic Oceans (Kim et al., 1997), and thus its recent appearance probably suggests warming of the southern waters of Korea. Our result pointed out that the δ 13 C and δ 15 N values of the two cluster groups may serve as isotopic baselines of the pelagic and benthic trophic pathways.
Feeding habits of predatory consumers might also be important in determining their consumption of pelagic-and benthic-affinity prey items. Food web analysis based on the contribution of different prey groups to the predatory groups in this study indicated the consumption of different mixes of the two prey groups (Figure 7). Predator group CI, to which the GPO belongs, was composed mostly of fish species (P. indicus, S. pachycephalus, at site MH; Stephanolepis cirrhifer, Mitrella bicincta, Ditrema temminckii, Pagrus major, and Latridopsis ciliaris at site PH) that utilize pelagic-and benthic-affinity prey at nearly equal proportions at the MH site and slightly higher benthic-than pelagic-affinity prey at the PH site. This result seems to be supported by their diverse spectrum of prey items, including amphipods, crabs, shrimps, polychaetes, echinoderms, small fish species, gastropods, bivalves, and cephalopods, as reported (Hur et al., 2006; http://www.fishbase.se; https://www. nifs.go.kr/frcenter/species). Two predatory (or scavenging) whelk species, Rapana venosa and K. lischkei, with isotopic niches similar to this group at the MH site, are known to prey on a variety of mollusks including clams, mussels, and oysters (Giberto et al., 2011;Vaux et al., 2017).
At the MH site, another predator group CII consisted of Stichaeopsis epallax, Girella punctata, Stichaeopsis sp., Alcichthys elongates, Sebastes ijimae, Furcina osimae, and C. myriaster. Their δ 13 C values were slightly more negative than those of group CI. As indiated by the SIAR mixing-model calculation, their exclusive dependence (82-90%) on pelagic-affinity prey reflected the preference for such prey including zooplankton, shrimps, and small fish species (Choi et al., 2008;https://www. nifs.go.kr/frcenter/species). In contrast, at the PH site, predator group CII had slightly less negative δ 13 C values compared with group CI, exhibiting an increased dietary dependence (65-73%) on benthic-affinity prey. Of the taxa belonging to this group, three species (S. pachycephalus, Platycephalus indicus, and Hexagrammos otakii) were reported to feed on the benthic invertebrates crabs, shrimps, and polychaetes, and small fish species 1 . The crab O. punctatus forages on benthic mollusks with prey species preference (Du Preez, 1984). The Buccinidae whelk K. lischkei feeds on clams, carrion, and detritus as a carnivore and scavenger (Vaux et al., 2017).
Based on the δ 15 N values of herbivores and suspension feeders as benthic and pelagic baselines, respectively, and assuming a trophic fractionation of 3.3 ± 0.26 (McCutchan et al., 2003), our TL estimation denoted that the marine ranching ecosystem is characterized by a trophic length of 3.5, as observed previously in the adjacent restored macroalgal bed (Kang et al., 2008) and natural rocky habitat . The gap in δ 15 N values between primary (herbivores) and tertiary (predatory fish) consumers recorded on average 2.6 to 4.6 (MH site) and 1.8 to 4.6 (PH sie). Primary consumers (suspension feeders, deposit feeders, and herbivores) had on average 1.5 to 3.0 (MH site) and 2.5 to 4.1 (PH site) higher δ 15 N values than those of primary producers. The trophic shift in δ 15 N values between secondary and tertiary consumers (e.g., predatory fish species) was relatively small (<2 ). The low δ 15 N values of fish species indicate the consumption of a wide spectrum of prey items from primary and secondary consumer groups (Deegan and Garritt, 1997;Kang et al., 2008).

Trophic Niche of E. dofleini in the Marine Ranching Food Web
We investigated the trophic niche of E. dofleini in two marine ranching ecosystems off the east coast of the Korean peninsula. The GPO stock is small in the study area compared with that further north (above 37.5 • N; Sohn et al., 2014;Yoon et al., 2016) and its occurrence is limited to winter-spring with small individuals in the study area, compared with all the year round in the further north of the study area with large individuals (Lee et al., 2014). This approximately establishes the southwestern margin of its distributional range. Although based on limited sampling sizes at one time of year due to the limited seasonal occurrence of E. dofleini, our results gave preliminary insight into carbon flow pathways through food web and trophic interactions between this GPO and its sympatric species in these marine ranching ecosystems.
The δ 13 C values of E. dofleini did not show any trend with increasing body size (Figure 3), probably reflecting the obscured relationship between sizes of prey and GPO due to the limited size range of GPO available (<4 kg in wet weight compared with <9 kg, Hartwick et al., 1981;Vincent et al., 1998). Given that individuals less than 6 kg in total weight in the GPO population off the east coast of Korea are immature (Lee et al., 2014), the GPO collected in the present study may be characterized as juvenile and immature individuals (Rigby, 2004). Its δ 13 C values were rather close to benthic-affinity prey, displaying a slight spatial variation. This result indicates an opportunistic feeding habit, which primarily exploits the most abundant prey type available (Vincent et al., 1998;Scheel and Anderson, 2012;Greenwell et al., 2019). The δ 13 C values of E. dofleini were similar to some teleosts, indicating an overlap in the primary carbon sources utilized between them (Kang et al., 2008). Indeed, the SIAR mixing-model calculation supports this explanation that they utilize both pelagic-and benthic-affinity prey. In contrast, the shifts in δ 13 C values between E. dofleini and other teleost groups further distinguished the pattern of resource exploitation between the GPO, which exploits a broad spectrum of pelagicand benthic-affinity prey items, and the teleosts, which possess either pelagic-or benthic-oriented feeding strategies (Gillespie et al., 1998;Scheel and Anderson, 2012), as shown by our conceptual food-web models (Figure 7).
No increase in δ 15 N values corresponding with the growth of the GPO was observed at either of the marine ranching sites. Based on the respective δ 15 N baseline values at each site, their mean δ 15 N values suggest that they occupy identical trophic positions between sites. The δ 15 N values of predatory fish species and GPO fell within a narrow range (9.7 to 11.1 and 9.7 to 12.6 at the MH and PH sites, respectively) with the exception of S. rivoliana (9.0 ). The observed narrow δ 15 N ranges suggests the absence of true tertiary consumers, creating a high complexity in feeding relationships by which the GPO might play a role as prey and/or predator in the study areas. The δ 15 N values of E. dofleini corresponded to the lower limit of the range for teleosts, probably reflecting exclusive ingestion of primary consumers (i.e., benthic and bentho-pelagic invertebrates) with little consumption of teleosts based on other factors such as strong avoidance and predation risk (Vincent et al., 1998). As a result, their foraging behavior seemed to be associated with general prey availability and accessibility rather than preference of prey species, as indicated by their narrow δ 13 C range (Scheel and Anderson, 2012;Greenwell et al., 2019). On the other hand, the inter-individual variability in both δ 13 C and δ 15 N values and trophic fractionation of E. dofleini depending on starvation and/or food limitation can be expected (Doi et al., 2017;Gorokhova, 2018). However, it is not possible to interpret this effect on the isotopic variations in GPO s from the present result.
The TL estimation indicates that E. dofleini individuals occupy intermediate trophic positions (3.2 and 2.7 at the MH and PH sites, respectively), slightly below the teleost value of 3.5 TL in the studied marine ranching food web. This indicates that the GPO exploits primary consumers (TL = 2) comprising pelagicand benthic-affinity prey as a predatory consumer. The slight shift in TL values between sites appeared to reflect a spatially varying consumption of predatory invertebrates. A similar mid-trophic position has been reported for Octopus tetricus consuming large volumes of abalone in a marine ranching food web (Greenwell et al., 2019). Earlier studies reported that GPOs are preyed upon by large fish, marine mammals (such as sea otters, harbor seals, sea lions, and killer whales), sharks, and even other octopuses (Haley, 1978;Hart, 1980;Gillespie et al., 1998). They are sometimes exposed to predation risk, losing their arms during mating (Huffard and Godfrey-Smith, 2010) and dead octopuses are scavenged by crabs and sea stars (Gillespie et al., 1998). However, as discussed earlier by the close proximity of their δ 15 N values, a lack of discrete trophic positions between E. dofleini and fish species appeared to be indicative of reduced predation by predatory fish species and instead the imposed competition for food between them.
The conceptualized food web architecture in these marine ranching ecosystems (Figure 7) demonstrates that autochthonous pelagic and benthic production supports E. dofleini and fish production through trophic mediation of primary consumers. Despite compositional changes in the taxa constituting individual trophic groups between the study sites, the presence of pelagic-and benthic-affinity prey (i.e., suspension feeders and grazer-deposit-feeder groups) in these systems provides a critical habitat feature to underpin feeding strategy of the GPO. Our isotopic analysis highlights that while GPO and some teleosts had a consistent trophic link with both pelagic-and benthic-based pathways between sites, other teleost groups showed prevalent trophic links with either pelagic-or benthic-based pathways, reflecting their feeding habits.
The installation of these marine ranches are acknowledged as a pilot program to enhance the sustainable productivity of the octopus fisheries by habitat improvement, nursing, and resource management. Once site selection by evaluation of the habitation and migration characteristics, and the following ranch creation are completed, any type of exploitation (fishing, recreational diving, etc.) is prohibited within the ranches. Their relevant activities include deploying artificial reefs, constructing kelp beds, and releasing juveniles of octopus. Further information of food web structure (i.e., ecological and trophic interaction) in the marine ranch ecosystems is needed to develop effective management strategies to maintain the habitats and enhance the octopus stocks.

CONCLUSION
This study investigated overall food web structure in association with the trophic niche of E. dofleini in the two nearshore marine ranches that were established for the GPO. Our findings may provide a great understanding of trophic ecology of the GPO in the natural rocky-bottom food webs. Our isotopic evidence revealed that the GPO had a trophic niche as a generalist feeder, relying on both benthic-and pelagic-affinity prey as natural diets at both marine ranches. Considering the opportunistic foraging behavior of E. dofleini (Vincent et al., 1998;Scheel and Anderson, 2012), the observed pattern of resource utilization indicates that nearshore natural rocky bottoms of the marine ranching areas retain faunal diversity and abundance (with different functional roles) to support both pelagic-and benthicbased food webs. This study offers an insight into their important role as predators of benthic and pelagic fauna and their possible role as a trophic competitor with upper trophic-level organisms in the functioning of the nearshore marine ecosystem in the southwestern margin of its range. Finally, the artificial reefs of marine ranches might provide shelter for GPO, and thereby the combination of shelter and dietary faunal composition appeared to provide an effective mechanistic role for the operation of marine ranching for the GPO.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
The animal study was reviewed and approved by Ministry of Oceans and Fisheries, South Korea.

AUTHOR CONTRIBUTIONS
HK and C-KK designed the experiments and identified overarching research goals. HK, C-KK, and D-LC wrote the manuscript. Y-JL, DK, D-HK, and J-HK conducted field observations and sample preparation and processing. C-KK and D-LC contributed reagents/materials/analytical tools. CK and DK performed stable isotope-ratio mass-spectrometer measurement. HK, Y-JL, and CK analyzed data and performed statistics. All authors reviewed the manuscript.

FUNDING
This work was supported by the GIST Research Institute (GRI) grant funded by the GIST in 2021.