Samplings were carried out during the KH-17-4 cruise on the RV Hakuho-Maru (Japan Agency for Marine-Earth Science and Technology) at three stations in the North Pacific Ocean between August 13 and October 2, 2017. Stations 1, 2, and 3 were located in the California Current (CC) region, the eastern subtropical gyre (ESG) region, and the western subtropical gyre (WSG) region, respectively (Figure 1A). Vertical profiles of water temperature and salinity between 0 and 500 m were recorded with a conductivity, temperature, and depth (CTD) system (SBE-911 plus; Sea-Bird Electronics). Chlorophyll a (chl a) was extracted with N, N-dimethylformamide, and the concentration was analyzed using a Turner fluorometer (Turner Designs).

Euphausiids were collected at night by oblique tows of Matsuda-Oozeki-Hu Trawl (MOHT; Oozeki et al., 2007) with a 2 m² mouth opening and a 1.95 mm mesh size from approximately 500 m depth to the surface (Table 1). A flow meter was attached to the net to measure the water volumes filtered. Samples were split on board promptly after net landing, and one-fourth of aliquots of the specimens were placed into a −80°C deep freezer immediately after splitting and kept frozen before subsequent laboratory operations for feeding ecology analysis. Another one-eight or one-fourth of aliquots of samples were fixed with 5% buffered formaldehyde for community structure analysis. To characterize the environmental communities as potential prey of euphausiids, seawater was collected at the surface (0 m) using an acid-cleaned bucket and at 10 m, the subsurface chlorophyll maximum (SCM), and at 100 m using acid-cleaned Teflon-coated Niskin bottles attached to the CTD. A 5 L aliquot of the water sample at each depth was immediately filtered on a Nuclepore polycarbonate membrane filter (Whatman, NJ, United States) with a 3.0 μm pore size, and the filters were also strictly preserved at −80°C. However, the seawater samples from station 3 were probably not processed timely, resulting in a large proportion of bacteria in the environmental communities according to the metabarcoding results. Thus, the environmental community data at station 3 was excluded from this study.

For the isotopic analysis, suspended particulate organic matter (POM) at 0 m was used as data of primary producers. In addition to the euphausiid samples, we analyzed net-plankton (size ranged from 0.2 to 1.0 mm) and myctophid (size ranged from 15.2 to 85.5 mm) in each station. For a sampling of POM, aliquots of 3–5 L surface seawater (0 m) collected using acid-cleaned buckets were filtered onto pre-combusted GF/F filters via gentle vacuum after pre-filtration with a 200 μm nylon mesh. Net-plankton samples were collected from a 200 m depth from the surface using vertical tows of a NORPAC net (45 cm mouth diameter, 100 μm mesh size), and then immediately size-fractionated into small (0.2–0.5 mm) and large (0.5–1.0 mm) fractions. Myctophids were selected from the micronekton samples collected by the same tow using a MOHT net for the euphausiid samplings. All samples were frozen at −20°C until laboratory processes.

All adult euphausiids at each station were identified to the lowest possible taxonomic level according to Baker et al. (1990) under a stereomicroscope (Nikon SMZ-745T), and the abundance for each species was standardized as individuals (ind.) 1,000 m^–3. The Shannon index (H’; Shannon, 1948) and the Simpson’s diversity index (D; Simpson, 1949) were calculated to compare the euphausiid species diversity among the stations.

The feeding ecology of the adults of the dominant euphausiid species at each station was further analyzed. Frozen samples were naturally thawed on ice, and 20–30 euphausiid individuals of the target species were randomly selected during the thawing. In the dissection process, the cephalothorax section was firstly excised to avoid the potential net feeding influence from the ingested food in the stomach of a euphausiid. To minimize DNA degradation, the whole gut was quickly removed and immediately transferred into a 30 μL of precooled 5% Chelex buffer (BIO-RAD) in a 1.5 mL Eppendorf tube kept on ice. Eyes of one of the selected individuals of each species were used as negative controls and treated in the same method as the gut contents. The remains of the euphausiid tissues were preserved back to −80°C until the stable isotope analysis. All the dissection facilities were sterilized every time before touching a different part of one euphausiid individual and dissecting another individual to avoid cross-contamination.

Total DNA of the gut contents and eyes were extracted according to the methods of Nagai et al. (2012). For the environmental samples, total DNA extraction was carried out with the MagAttract PowerWater DNA/RNA Kit (QIAGEN). Polymerase chain reaction (PCR) amplifications were conducted for all samples using KOD Plus version 2 (Toyobo), with the eukaryotic universal primers 1389F and 1510R (Amaral-Zettler et al., 2009) according to the methods of Hirai et al. (2017). Adaptors used to differentiate samples during sequencing on Illumina MiSeq and dual-index target sequences were added in the 2^nd and 3^rd rounds of PCR. The PCR products of every round were diluted 20 times with diethylpyrocarbonate (DEPC)-treated water and used as templates for the next round. The final PCR products were purified using the QIAquick PCR Purification Kit (QIAGEN). Concentrations of the final PCR products were measured using a Qubit 3.0 Fluorometer (Invitrogen), and 100 ng of each purified PCR product was pooled. A sequencing run on Illumina MiSeq was carried out using MiSeq Reagent Kit v2 (FASMAC Co), and raw data of 2 × 250 bp paired-end sequence reads were obtained.

All raw sequence reads were quality-filtered using Trimmomatic (Bolger et al., 2014). Short- and low-quality fragments were eliminated, sequencing adaptors were clipped, and all technical sequence reads were trimmed into appropriate length using the following settings: CROP: 135; MINLEN: 50; LEADING: 20; TRAILING: 20; and SLIDINGWINDOW: 30:30. Paired-end sequences were merged, and further quality-filtrating was performed using MOTHUR (Schloss et al., 2009) according to the methods of Hirai et al. (2017). Briefly, the reads with more than six homopolymers and uncertain bases were removed, and qualified sequences were aligned against SILVA 138.1 database (Quast et al., 2013). Aligned sequences were single-linkage pre-clustered and screened to remove potential chimeric sequences with a reference dataset in UCHIME (Edgar et al., 2011). In this study, the reference database V9_PR2 from the Tara Oceans project (de Vargas et al., 2015) was used to classify the remaining high-quality sequences into different taxonomic groups using a naive Bayesian classifier (Wang et al., 2007), with a threshold greater than 80%. Only the Eukaryotic sequences were preserved for further analysis, and sequences classified into “Malacostraca” and “Arthropoda_unclassified” were excluded as a possible read from the euphausiid host. Mammalian, fungi, and parasitic taxa, which were determined according to Cleary and Durbin (2016), were also removed. Finally, preserved sequences were clustered into operational taxonomic units (OTUs) at a 99% similarity threshold.

After quality controlling and taxonomic classification, to present the gut content composition of euphausiid individuals and the environmental community structure of different stations, the relative sequence read abundances of all taxonomic supergroups in Eukaryota were measured for each euphausiid and environmental sample. A Kruskal-Wallis test with Dunn post hoc test was conducted on the relative abundances of main supergroups in each sample within and across different euphausiid species using IBM SPSS statistics 23 (Zar, 1999). Moreover, to evaluate the gut content composition differences among different euphausiid species, the Bray-Curtis similarity index was first calculated on euphausiid gut content composition of square root-transformed data, and the homogeneity of dispersion was evaluated using permutational analysis of multivariate dispersions (PERMDISP). Differences among species were then tested using permutational analysis of variance (PERMANOVA) and a posterior pairwise test. PERMANOVA was conducted using Type III sum of squares and the unrestricted model, and number of permutations was 9999 for all permutation-based tests. To detect those prey supergroups responsible for the possible observed differences between euphausiid diets, a similarity percentages (SIMPER) analysis was also performed. All these analyses were conducted in PRIMER version 7 with the add-on package PERMANOVA+ (Anderson et al., 2008). Furthermore, to clarify main preys of all the euphausiid species with detailed taxonomy information, the percentage of each prey OTU in each sample was computed and then averaged among all euphausiid individuals at each station. For all three euphausiid species, detected prey OTUs comprising more than 1.5% of the total gut content reads and with occurrence frequency higher than 40% were deemed as dominant prey OTUs (Harms-Tuohy et al., 2016; Deagle et al., 2019). Every dominant prey OTU was identified to the lowest possible taxonomic level by BLASTN searches against the NCBI “nt” database.

Particulate organic matter (POM), net-plankton, and euphausiid samples for stable isotope analysis were oven-dried at 50°C overnight. Filter samples of POM were exposed to HCl fumes for 2 h to remove the carbonate and dried again at 50°C. All euphausiid individuals and 1–10 mg of the net-plankton samples in the two fractions were homogenized, and approximately 0.5–1.5 mg of each homogenized sample was measured (two euphausiid individuals were pooled when the sample weight was insufficient). For the myctophid samples, the dorsal muscles of 11–14 specimens from each station were dried and homogenized. Half of the sample was lipid-extracted using chloroform/methanol (2:1, v/v) for the measurement of δ¹³C (Folch et al., 1957).

Stable nitrogen and carbon isotope ratios were measured using a DELTA V advantage mass spectrometer (Thermo Fisher Scientific Inc.) connected to an elemental analyzer. Amino acids with different isotopic values (alanine, histidine, and glycine), which were certified using the standard sample provided by the International Atomic Energy Agency, were used to obtain the calibration curves. Alanine (δ¹⁵N = −1.2‰, δ¹³C = −19.9‰) was measured once every 10 samples for correction of a shift of the baseline. As for the trophic position analysis, Minagawa and Wada (1984) reported an average stepwise increase in δ¹⁵N of 3.4 ± 1.1‰ per trophic level across a wide range of habitats, including marine habitats. Thus, the trophic levels of the analyzed samples were calculated as:

Differences in the trophic levels at the different sampling stations were analyzed by one-way analysis of variance (ANOVA) with Fisher’s least significant difference (LSD) post hoc test using IBM SPSS Statistics 23.

A lower temperature and salinity were observed at station 1 in the CC compared to the other two stations in the subtropical gyre (Figure 1B). At station 1, the temperature decreased gradually with increasing depth, from 17.6°C on the surface to 5.3°C at 500 m, and salinity increased slightly from 33.2 to 34.1. Stations 2 and 3 in the subtropical gyre shared more similar environmental conditions, especially in the epipelagic layer: the water temperature of station 3 was higher at the surface and below 200 m, but was nearly the same as station 2 at 50–200 m (average of 20.2°C in station 2 and 19.7°C in station 3); the salinity of these two stations were in common (average of 34.9) in the upper 200 m. A higher chl a concentration (0.13 μgL^–1) at the surface (0 m depth) was observed at station 1 than at station 2 (0.07 μgL^–1) and 3 (0.06 μgL^–1). Deeper SCM depths were observed at station 2 (117 m) and station 3 (115 m), compared with station 1 (70 m). The chl a concentration at the SCM depth of each region was the highest at station 3 (0.59 μgL^–1), followed by station 1 (0.37 μgL^–1) and station 2 (0.24 μgL^–1).

The total abundance of euphausiids was the highest at station 1 (152 ind. 1,000 m^–3) and the lowest at station 3 (48 ind. 1,000 m^–3 at station 2, 47 ind. 1,000 m^–3 at station 3), and 22 euphausiid species belonging to six genera were identified (Figure 2 and Table 2). E. pacifica was the most abundant species at station 1, contributing approximately 80.0% of the abundance (121 ind. 1,000 m^–3) of the euphausiid community. Other euphausiids at station 1 were mainly attributed to E. gibboides and Nematoscelis difficilis. E. brevis and E. hemigibba were the dominant euphausiid species at stations 2 and 3, respectively. However, in comparison with E. pacifica at station 1, dominance of these two subtropical species was considerably low, 24 ind. 1,000 m^–3 for E. brevis and 11 ind. 1,000 m^–3 for E. hemigibba, accounting for 50.0% and 22.8% abundance at stations 2 and 3, respectively. The euphausiid community compositions at stations 2 and 3 were similar at the genus level. Euphausiid species in the genera Stylocheiron and Thysanopoda were detected frequently at both stations 2 and 3 but not at station 1. Abundances of these genera were higher at station 2 than at 3. In contrast to the species abundance, H’ and D were higher in the subtropical gyre (stations 2 and 3) than at station 1 of CC area. The highest species diversity was observed at station 3 (Table 2).

Totally, 79 adult euphausiid individuals (21 E. pacifica, 27 E. brevis, and 31 E. hemigibba) were selected for the feeding ecology study. The total length of euphausiids was 7.6–15 mm for E. pacifica (mean ± standard deviation (SD): 10.0 ± 2.1), 7.0–9.1 mm for E. brevis (mean ± SD: 8.2 ± 0.5), and 9.4–13.2 mm for E. hemigibba (mean ± SD: 11.4 ± 1.1). An overall number of 384,698 sequence reads from eight environmental samples in CC and ESG provinces and 262,758 reads from euphausiid gut samples were obtained. However, 20,533 reads from three negative-control (eye) samples still remained after removing the host and contaminant sequences. Furthermore, >98% of the remaining sequences of negative controls as well as 47.8% of the gut content sequences were classified into “unclassified Crustacea.” We, therefore, double-checked all the OTUs in the negative controls and found that ∼80% of the unclassified OTUs in negative controls also appeared in the gut contents, taking up 95–97% of the unclassified sequences in gut contents. The unclassified OTUs were eventually verified as euphausiid species and other arthropods or remained unclassified based on the BLASTN search against the “nt” database (Supplementary Figure 1). The unclassified sequences discovered in the gut contents were highly likely attributed to the host sequences; thus, all the unclassified sequences in the gut were removed manually. In addition, for each euphausiid species, gut samples that contained insufficient reads (<0.5% of the total gut content sequence reads) were excluded. At last, a total of 12 E. pacifica, 27 E. brevis, and 28 E. hemigibba were used, and sequence reads for final data analysis were 112 to 46,697 for euphausiid guts, and 32,381 to 61,099 for environmental samples.

Gut contents of all three euphausiid species were mainly dominated by the phytoplankton taxonomic groups including Dinophyta, Stramenopiles, and Archaeplastida (Figure 3). Within the phytoplanktonic supergroups, the Kruskal-Wallis test suggested that Dinophyta sequences had significantly higher (p < 0.05) relative abundances than other prey taxa in the guts of all three euphausiid species and were detected in every individual, accounting for 3.3–87.8% (average 45.8%) in E. pacifica, 1.6–97.6% (average 40.5%) in E. brevis, and 4.2–97.6% (average 56.3%) in E. hemigibba. Dinophyta was also abundant in the environmental communities of both CC and ESG regions, averagely representing 29.1% at CC (station 1) and 36.8% at ESG (station 2) in sequence proportions. Stramenopiles was the second most abundant in sequence proportions across the phytoplankton found in the guts, averaging ∼10% in all three species (8.7% in E. pacifica, 8.4% in E. brevis, and 11.7% in E. hemigibba), but comprised a smaller proportion in the water columns (average 4.2% in CC and 1.6% in ESG). Archaeplastida showed relatively small proportions of sequence reads in the water columns, with greater abundance in CC (average 1.8%) than ESG (average 0.2%) (Figures 3B,C), however, Archaeplastida displayed higher percentages in gut contents of euphausiids, making up on average 2.8% in E. brevis, 6.7% in E. hemigibba, and 8.4% in E. pacifica. Kruskal-Wallis test indicated that the relative abundance of Archaeplastida was significantly different between E. brevis and the other two species with p < 0.05, and percentage of Archaeplastida in E. pacifica and E. hemigibba were similar with each other. Moreover, Archaeplastida displayed a significant lower (p < 0.05) relative abundance than Dinophyta and Stramenopiles within E. brevis. Except of phytoplankton, Copepoda sequences represented a similar minor portion in all three species of euphausiids (average 5.7% in E. pacifica and E. brevis, 4.7% in E. hemigibba), but was the most abundant in the environmental community sequence reads. Sequence reads from Rhizaria were relatively abundant in ESG environmental communities (average 7.0% in ESG, 4.8% in CC), and significantly higher percentage (p < 0.05) of Rhizaria was discovered in guts of E. brevis (average 6.3%) than in guts of E. hemigibba (average 1.4%). Rhizaria was found to present a significantly lower percentage than other preys within E. hemigibba. Additionally, Hydrozoa sequences, with a trace amount in the water column (<0.5%), was strikingly discovered in several euphausiid individuals, resulting in a prominent proportion of the gut contents of all three euphausiid species (average 16.6% in E. pacifica, 25.6% in E. brevis, and 8.4% in E. hemigibba). PERMDISP showed no significant differences (p > 0.05) in dispersions among three species, and PERMANOVA indicated a significant difference (p < 0.05) in the gut content composition between E. brevis and E. hemigibba (Supplementary Table 1). The dissimilarity mainly attributed to Dinophyta, Stramenopiles, and Rhizaria with relative greater average dissimilarities and distributed evenly (larger dissimilarity/standard deviation ratios) among euphausiid samples. Hydrozoa also largely contributed to the dissimilarity between two species, but the contributions were inconsistent among all the samples (lower Diss/SD ratio; Supplementary Table 2).

A total of 19 OTUs were verified as primary prey of three euphausiid species, including phytoplankton (dinoflagellates, pelagophytes, and green algae) and zooplankton (copepods, hydrozoans, and rhizarians) (Table 3). Hydrozoan OTUs, generally identified as siphonophores, were one of the most dominant OTUs in all three euphausiid species, present in >90% of the euphausiid individuals. Among phytoplankton OTUs, dinoflagellates were the uppermost food items in all three euphausiid species (average of 4.6–14.2% sequence reads), with higher occurrence in E. hemigibba (71%) than E. pacifica (50–58%) and E. brevis (52%). Kareniaceae dinoflagellates were predominant dinoflagellate prey for E. brevis and E. hemigibba. Pelagophytes OTUs were also identified as prime prey for E. pacifica (2.0–4.1% average reads) of CC region and E. hemigibba (3.9% average reads) of WSG region. Pelagomonadales OTUs were common prey of both E. pacifica and E. hemigibba; however, Coccoid pelagophyte OTU was detected only in 58% of the E. pacifica individuals. Green algae OTUs were only listed as dominant in the E. pacifica with a high frequency (58–75%). As for zooplankton prey, Pontellidae copepods for E. brevis, Clausocalanidae copepods for E. hemigibba and E. pacifica were prime food of three species of euphausiids. These copepod OTUs played more important roles in the gut contents of E. brevis and E. hemigibba in the subtropical areas (2.1–2.4% sequence reads; 74–77% occurrence frequency) than E. pacifica (1.7% sequence reads; 67% occurrence frequency). Sphaerozoidae rhizarian was uniquely dominant prey for E. brevis in ESG.

The highest mean nitrogen isotope ratio was in E. pacifica (10.9 ± 0.8‰) and the lowest in E. brevis (3.6 ± 0.2‰), which was consistently 5–5.5‰ higher than that of POM (Figure 4A). The mean δ¹⁵N values of the euphausiid species were higher than those of the small or large net-plankton and were lower than those of the myctophid fish. This relationship was also verified between the mean δ¹³C values of POM, net-plankton, euphausiids, and myctophid fish in both ESG (station 2) and WSG (station 3) regions, however, the mean δ¹³C value of the E. pacifica was higher than that of the myctophid fish in CC region. The highest value of δ¹³C was observed in E. brevis (−20.3 ± 0.2‰), and the lowest was observed in E. pacifica (−21.9 ± 0.5‰). The mean trophic levels of E. pacifica, E. brevis, and E. hemigibba were 2.32, 2.48, and 2.57 (Figure 4B), respectively, indicating that the three euphausiid species were between primary and secondary consumers. The trophic levels of E. brevis and E. hemigibba were significantly higher than those of E. pacifica; however, no significant difference was found between E. brevis and E. hemigibba. A relatively wider trophic level spectrum was observed in E. hemigibba (2.07–2.77) in WSG province, compared with the other subtropical species E. brevis (2.37–2.61) in ESG province.

This study aimed to reveal the diets of dominant species of euphausiids in the CC and subtropical areas, and the community structures of euphausiids in three water areas were firstly investigated. As reported by Mauchline and Fisher (1969), the hydrographic characteristics of water masses were important factors governing the zonation pattern of euphausiids species in our study, and the euphausiid community in the CC region was distinct from that in the subtropical area. Furthermore, the sea surface temperature is known to be an important driver of diversity in pelagic taxa (Letessier et al., 2009, 2011; Tittensor et al., 2010), agreeing with the diversity patterns of euphausiids in our study. Station 1 (CC area) was in a transition zone with the merging of subarctic water, subtropical water, and equatorial water of the North Pacific, and E. pacifica were highly dominant. A preferred water temperature of E. pacifica ranged between 5 and 16°C, and their abundances were centered in the CC region (Taki, 2007; Lara-Lopez et al., 2012; Tao et al., 2015). Because the warm temperature in the subtropical region around 20–25°N is unfavorable for E. pacifica (Brinton, 1962), this species was not found in both WSG and ESG in our study. On the other hand, most of the predominant euphausiid species of station 2 and 3 were absent in CC province. Roden (1991) noted that the southern boundary of the transition zone (32°N) plays a role as a geographical barrier preventing subtropical species, such as E. brevis from northward extending their distribution. The subtropical species E. brevis and E. hemigibba usually co-occur in the oligotrophic subtropical gyre together with N. atlantica and Stylocheiron suhmi (Moore, 1951; Fager and McGowan, 1963; Wiebe and D’abramo, 1972; Lindley, 1977). However, there were environmental differences in water temperature and chl a concentration within the North Pacific subtropical gyre, possibly leading to different dominant species between ESG and WSG. Based on the results of the present study, more sampling stations in the subtropical regions should be explored for further investigation of euphausiid community structure in related to the hydrographic gradients, and more information about trophic interaction between dominant and minor euphausiid species would provide more insights. The community structures and dominant species, whose diets are further discussed in the following, observed in this study were thus typical in the North Pacific.

All the three Euphausia species in this study presented an evident omnivorous feeding. Both phytoplankton (dinoflagellates, pelagophytes, and green algae) and zooplankton (copepods and hydrozoan) were frequently observed from the gut contents by 18S V9 metabarcoding, and omnivorous feeding was also confirmed by the stable isotope analysis. These results also indicated that the feeding pattern of euphausiids was largely influenced by the food resource availabilities of phytoplankton in the environment. An omnivore feeding has been well studied in E. pacifica. In addition to major prey items such as diatoms and ciliates, E. pacifica has been verified to feed on dinoflagellates (Gyrodinium sp. and Protoperidinium sp.) and small copepods (Pseudocalanus sp.) using microscopic observation, stable isotope analysis, and rearing experiments (Dilling et al., 1998; Nakagawa et al., 2001, 2004; Park et al., 2011; Du and Peterson, 2014; Okazaki et al., 2020). In addition, a similar trophic level of 2.5 for E. pacifica in the Oyashio Current was reported by Sogawa et al. (2017) using stable isotope analysis. However, feeding of subtropical euphausiid species have been poorly reported, and only a few studies found foraminifera and unclassified copepods in E. brevis (Mauchline, 1967). Thus, our study was the first to report detailed feeding habits of the dominant species of E. brevis and E. hemigibba in the North Pacific subtropical gyre.

Dinoflagellates are widely recognized as important prey items of pelagic taxa (McClatchie, 1988; Turner et al., 2001; Pethybridge et al., 2014), and as an essential energy supply of all developmental stages of E. pacifica using dinoflagellate-specific fatty acid biomarkers 18:4n-3 and 22:6n-3 (Kohlbach et al., 2019; Fisher et al., 2020). Although some species of Kareniaceae (e.g., Karlodinium sp.), which were detected as dominant OTUs in the present study, could deter grazing from predators by producing karlotoxins, zooplankton such as Acartia tonsa would potentially graze on dinoflagellate strains with lower toxin levels (Waggett et al., 2008). Furthermore, pelagophytes and green algae have not yet been reported in euphausiid diets possibly due to their small size, however, our results showed their importance as food sources of dominant euphausiids. Copepods were also detected as a major prey, even though proportions of copepods were much lower in gut contents of euphausiids than environmental communities in metabarcoding analysis. Cleary et al. (2018) similarly discovered copepod sequence reads in the stomach contents of E. superba, indicating a carnivory supplementary diet of this species. The importance of copepods was verified by Nakagawa et al. (2002) based on carbon content analysis, and E. superba was supposed to feed on Ctenocalanus copepods (Clausocalanidae) according to their overlapping vertical distribution (Marrari et al., 2011).

In the present study, all dominant species belonged to the genus Euphausia. According to comparative morphological studies, the thoracic legs of Euphausia species have a nearly uniform length and form a well-developed feeding basket for omnivorous filter-feeding (Mauchline, 1967), in comparison with other euphausiid species with one or more pairs of elongated thoracic legs (Berkes, 1975; Dalley and McClatchie, 1989; Suh and Choi, 1998). For example, N. difficilis and Nematobrachion boopis displayed higher δ¹⁵N values indicating more carnivorous feeding than Euphausia species in the same water area (Kinsey and Hopkins, 1994; Sogawa et al., 2017). Moreover, although the feeding basket structure of Thysanopoda spp. is similar to that of Euphausia spp., mesopelagic Thysanopoda species are more carnivorous as they keep away from the primary producers (Nemoto et al., 1977). Thus, the diet spectrum of the minor euphausiid species was bottlenecked by the inconvenient limbs for filter feeding and their mesopelagic habitat; Euphausia species, however, have their nutritional needs satisfied by an omnivorous dietary ranging from picophytoplankton to microzooplankton. Omnivorous feeding habit enables Euphausia species to adapt to various prey availabilities and become the dominant species of euphausiid communities.

Although all the three species are omnivorous, and similar gut content compositions were captured for E. pacifica and subtropical euphausiid species, E. pacifica in CC displayed a different feeding habit from E. brevis and E. hemigibba reflecting by lower trophic levels verified by stable isotope analysis. In addition, DNA metabarcoding detected more autotrophic preys (green algae) with greater abundance and occurrence frequency in E. pacifica by comparison with subtropical E. brevis. All the three dominant species conduct diel vertical migration: surface migrant E. pacifica in the northwestern Pacific were most abundant at 5–25 m for feedings during nighttime (Nakagawa et al., 2003; Endo and Yamano, 2006; Sogawa et al., 2016), and E. brevis were most abundant at upper 50 m in the Mediterranean Sea and Bermuda (Wiebe and D’abramo, 1972; Youngbluth, 1975). Despite a similar vertical distribution at night between E. pacifica and E. brevis, higher prey availability of pelagophyte (Stramenopiles) was observed in CC (0 m: 4.8%; 10 m: 4.7%) compared to that in ESG (0 m: 1.5%; 10 m: 1.6%), in addition to more green algae (Archaeplastida) prey (CC: ∼2.5%; ESG: ∼0.2%) (Figures 3B,C). Thus, E. pacifica was allowed access to more autotrophic prey. Moreover, E. pacifica would display a high growth rate when feeding on phytoplankton, such as diatoms, rather than copepods (Lasker, 1966; Ohman, 1984). Autotrophic pelagophyte and green algae are supposed to be the optimal prey of E. pacifica in the CC area, and E. pacifica would depend more on these prey that probably provides them with the optimum ratio of energy gain per energy expenditure (Schoener, 1971).

A lower chl a concentration in oligotrophic environments would lead to an increasing heterotroph (e.g., dinoflagellates, zooplankton) contribution to the whole diet of euphausiids, such as E. lucens, E. spinifera, E. superba, and E. pacifica. (Stuart and Pillar, 1990; Gibbons et al., 1991; Perissinotto et al., 2000; Nakagawa et al., 2002; Henschke et al., 2015). Non-autotrophic prey (copepods and mixotrophic/heterotrophic dinoflagellates) were important foods of E. hemigibba and E. brevis due to a potential poor autotrophic prey availability in the environment. In addition, the lower chl a level in ESG than that in CC (Figure 1B) triggered the lesser uptake of green algae of E. brevis in the present study. Furthermore, the gut content composition of two subtropical euphausiid species was significantly different (Supplementary Table 1): metabarcoding verified a relatively higher contribution of Dinophyta and Stramenopiles in E. hemigibba than E. brevis, while E. brevis consumed significantly more Rhizaria (Supplementary Table 2 and Figure 3). These differences might be attributed to a higher chl a concentration at the subsurface layer in WSG than ESG. Because E. hemigibba were found more frequently in deeper layers (100–150 m) in northeastern Japan and the Mediterranean Sea (Wiebe and D’abramo, 1972; Taki, 2008), they might utilize phytoplankton at the subsurface layer effectively in WSG. On the other hand, trophic level by stable isotope analysis was higher in E. hemigibba than E. brevis, with high variations and no significant difference. This discrepancy in metabarcoding and stable isotope analysis implies that feeding patterns of E. hemigibba and E. brevis change temporally according to food availability in the environments, which also underlines the need of the combining of DNA metabarcoding and stable isotope analysis in feeding habit investigations.

The high proportion and occurrence of hydrozoan sequences detected in the gut contents of the euphausiids were unexpected since hydrozoa (jellyfishes) have been known as predators of copepods, euphausiids, and fish larvae (D’Ambra et al., 2014; Fleming et al., 2015). Likewise, hydrozoan sequences have been recently detected in the gut contents of myctophid fishes of the Southern Ocean, eel larvae of the Northern Atlantic Ocean, and calanoid copepods especially in a low-phytoplankton environment (Yi et al., 2017; Ayala et al., 2018; Hirai et al., 2018; Chow et al., 2019; Yeh et al., 2020), suggesting a previously overlooked importance of gelatinous plankton in pelagic ecosystems. However, hydrozoan as natural prey of euphausiids is contradicted to the intermediate trophic positions of euphausiids uncovered by stable isotope analysis and the minute proportion of Hydrozoa sequences detected by 18S metabarcoding in the environment in our research. One possibility is that euphausiids feed on the early life stages of Hydrozoa or marine snow aggregates combining hydrozoa particles (Dilling et al., 1998; Turner, 2015; Yi et al., 2017; Ayala et al., 2018). Thus, marine snow, as another food source of euphausiids, is suggested to add to the sample collection and be analyzed in the same way as environmental samples in future studies. On the other hand, euphausiids are likely to take hydrozoa atypically in the cod-end during sampling. The feeding pattern difference between two subtropical euphausiid species were unevenly contributed by Hydrozoa among euphausiid individuals indicated by SIMPER, which actually supported this theory. Nicol (1984) noted that Meganyctiphanes euphausiids collected by towed nets within 10 min at sea surface would contain more copepods in their “food baskets” when compared with those collected by hand-held dip nets because potential food items become highly concentrated in towed nets; Judkins and Fleminger (1972) also found greater food diversity in sergestids caught by nets in comparison of those obtained from fish stomachs. Given that the sampling duration was relatively long in the present study, although the mouthpart and stomach of euphausiids were excluded during the dissection process, cod-end feeding of euphausiids might have happened and impacted the gut content analysis. Gut content analysis of euphausiids collected by shorter tows from the surface swarms and additional stable carbon and nitrogen isotope analyses dealing with hydrozoans inhabit in each study site are assured to improve the accuracy of euphausiid diet and trophic linkage study.

DNA metabarcoding has been successfully applied to gut contents or feces of a series of organisms ranging from copepods and fish larvae to marine mammals to unveil their dietary diversity, dietary composition, and the contribution of different prey items to their overall diet (Hardy et al., 2017; Zamora-Terol et al., 2020; Kume et al., 2021). The most accepted advantage of DNA metabarcoding over traditional visual analysis is the lesser dependency on the digestion state of prey items. A stomach content analysis of catfish by DNA metabarcoding recovered 91.6% of prey items to species level compared to morphological identification (9.4% to species level, 12.1% to family level), and little difference was shown between the major composition of lightly digested piscine prey items and moderately digested ones (Aguilar et al., 2017). The firstly reported pelagophyte and green algae with evident larger contribution and higher occurrence frequencies in euphausiid diet of our study, once again highlighted the high efficiency and sensitivity of DNA metabarcoding. However, even though potential prey could be easily explored, the natural dietary composition might be biased due to the high sensitivity. Similar to cod-end feeding impacts, cross-contamination derived from body surface during dissection might be left in the gut content sequences (De Barba et al., 2014; Chow et al., 2019); fungi, symbiotic, and parasitic organisms originating from food or digestive tract are difficult to be verified as a true prey due to insufficient parasite knowledge and the limited registered sequences in 18S databases (Hirai et al., 2017; Kodama et al., 2017; Cleary et al., 2018). Furthermore, although V9 hypervariable region of 18S rDNA could reveal the extant diversity of eukaryotes, unclassified OTUs and prey OTUs with lower identify (<85%) during classification ought to be cross-checked by multiple genetic markers and supplementary databases (Da Silva et al., 2019; Choi and Park, 2020). Additionally, trace signals of secondary-prey are likely to remain because of the usage of universal primers, however, these items were generally in low occurrence, thus could be excluded (Hardy et al., 2017; Clarke et al., 2020).

A final challenge of not only feeding ecology study but all the research regarding DNA metabarcoding is how to interpret the sequence reads of different taxonomic groups in relation to their practical biomass (Pompanon et al., 2012; Bucklin et al., 2016). Despite sample size, PCR biases, amplicon efficiency, and different copy number of selected barcodes in different prey items are known to cause different degrees of recovery biases (Deagle et al., 2013). An additional weighing on samples based on gut content mass or fullness, importing blocking primers to inhibit host-sequence amplification, and excluding empty guts like the present study are supposed to contribute to the prey item quantification. By far, the stable isotope and fatty acid analysis are more committed methods for food source quantification, whereas DNA metabarcoding is more suitable to pilot diet investigation of organisms without prior knowledge (Stewart et al., 2017; Kohlbach et al., 2019). Combining the approach of DNA metabarcoding-based gut content analysis and stable isotope analysis is supposed to be the most appropriate way to mirror a convinced feeding ecology of studied consumers. Given all the inadequacies, in the future, feeding ecology study of euphausiids via DNA metabarcoding, we strongly recommend: (1) using disposable dissection tools on different individual specimens to avoid contaminations; (2) employing empty guts and other organs of euphausiids as negative controls to minimize the noises of fungi, parasites, and unclassified sequences; (3) importing additional marine snow samples to provide a more complete food resource availability; (4) using multiple metabarcoding genetic markers and additional reference databases to advance the integration of dietary data; (5) considering the frequency of occurrence as a criterion and relative read abundance for true-prey assessment instead of judging based on a presence/absence data alone.