Zooplankton samples were collected in the Ría de Vigo (NW Spain) onboard RV “Mytilus” in 2012 and 2014. The timing of the sampling was based on previously identified periods of maximum paralarval abundance (Rocha et al., 1999; González et al., 2005) in 2012 and 2014: we carried out ten nocturnal surveys each year, four in summer (July), and six in early autumn (September and October). Additionally, diurnal surveys were conducted in summer and in autumn 2012 (one per season). Sampling surveys were conducted along four transects (Figure 1A). For each transect, a Multinet® Hydrobios Mammoth of 250 μm mesh size, fitted with two electronic flow meters, was lowered at 2.5 knots to the sea floor and lifted up gradually to the surface. We defined seven depth layers: from 105 to 85 m, Z7; 85 to 55 m, Z6; 55 to 35 m, Z5; 35 to 20 m, Z4; 20 to 10 m, Z3; 10 to 5 m, Z2; and 5 m to the surface, Z1; see Figure 1B). Within each layer, the Multinet® filtered up to 200 m³ of seawater (approximately from 5–10 min for each layer and hence in total between 20 and 70 min), and collected independent samples. The collected zooplankton was fixed onboard in 96% ethanol and frozen at −20°C until sorting. In the laboratory, all cephalopod paralarvae were separated and preserved individually in 70% ethanol and stored at −20°C.

Dorsal mantle length (DML) was measured to the nearest 0.05 mm on the dorsal side of all octopus and squid paralarvae using a Leica M205C stereomicroscope and Leica Application System image analysis software (Leica Microsystems, Germany). All octopus paralarvae (n = 492) were identified as O. vulgaris based on morphological characters following Sweeney et al. (1992). Due to the difficulty of identifying squid paralarvae using morphological characters, all loliginid paralarvae (n = 163) were identified genetically. Molecular identification relies on previous work with adult specimens identified morphological and subsequently genetically.

Briefly, DNA from the mantle of each loliginid paralarva was extracted with a QIAamp DNA Micro Kit (QIAGEN) following manufacturer's instructions, with the exception of two steps: Digestion at 56°C was done overnight and the final elution was done in two steps using 15 μl buffer AE in each elution. The barcoding region of the Cytochrome c Oxidase subunit I (COI) was amplified with the universal primers HCO2198 and LCO1490 (Folmer et al., 1994) and PCR products were sequenced by Sanger sequencing (Stab Vida, Portugal). Each sequence was compared to the following GenBank reference sequences using the BLAST algorithm (Altschul, 2014): Alloteuthis media, EU668085 (Anderson et al., 2008); A. subulata EU668098 (Anderson et al., 2008), and L. vulgaris, KF369142 (Lobo et al., 2013).

Loliginid paralarvae were identified as A. media (n = 93), A. subulata (n = 35), and L. vulgaris (n = 22) (Olmos-Pérez et al., unpublished data). For dietary analysis, a total of 64 O. vulgaris paralarvae (summer, n = 26; autumn, n = 38) and 32 A. media (summer, n = 16; autumn, n = 16) were selected. Alloteuthis media was selected for molecular diet analyses because it was the most abundant loliginid present. Samples were chosen to maximize information from different seasons, transects and depth.

Digestive glands of all 96 paralarvae were dissected out, cleaned with sterile distilled water and placed into DNA-free tubes (Suzuki et al., 2006). DNA was extracted with a QIAamp DNA Micro Kit (QIAGEN, Hilden, Germany) following the modifications in the elution step as stated before. DNA purity and concentration were controlled with NanoDrop 2000c UV-Vis Spectrophotometer (Thermo Fisher Scientific Inc., Massachusetts, USA).

Many dietary studies recommend the use of restriction enzymes (Blankenship and Yayanos, 2005) or blocking primers (Vestheim and Jarman, 2008; Deagle et al., 2009; Leray et al., 2013a) to avoid amplifying predator DNA and maximize detection of prey DNA. However, the huge number of sequences currently obtained with NGS platforms allows the use of universal primers that facilitate the detection of unexpected prey (Boyer et al., 2013) without the necessity of using restriction enzymes or blocking probes (Piñol et al., 2014). Accordingly, we employed the universal pair of primers HCO2198 (Folmer et al., 1994) and mlCOIintF (Leray et al., 2013b), to amplify 315 base pairs (bp) of the barcoding region of the mitochondrial cytochrome c oxidase subunit I (mt-COI) gene (Table 1). Cycling conditions for the touch-down PCR with COI primers were: initial denaturation at 95°C 3 min, 10 initial cycles of denaturation at 95°C for 30 s, annealing for 30 s at 57 (−1°C per cycle), and extension at 72°C for 40 s, followed by 29 cycles of denaturation at 95°C for 30 s, annealing at 47°C for 30 s, extension at 72°C for 40 s, and a final elongation for 4 min at 72°C. PCR amplification was performed in a total volume of 25 μl:1 μl (10 μM) of each forward and reverse primers, 12.5 μl Thermo Scientific™ Phusion™ High-Fidelity PCR Master Mix with HF Buffer (Thermo Fisher Scientific Inc., Massachusetts, USA), 1 μl of DNA (20 ng/μl), and 9.5 μl H₂O.

Since decapods, krill and fishes had been previously detected in the digestive tract of O. vulgaris paralarvae (Passarella and Hopkins, 1991; Vecchione, 1991; Roura et al., 2012), we also employed two pairs of specific primers to amplify the 16S mitochondrial gene of malacostracan crustaceans (mt-16Sa, Table 1) and chordates (mt-16Sb; Table 1). Cycling conditions for both 16S pairs of primers were: initial denaturation at 94°C 15 min, 33 cycles of: denaturation at 94°C for 20 s, annealing at 48.7°C for 90 s, extension at 72°C for 45 s, and a final elongation step at 72°C for 2 min. PCR amplification was performed in a total volume of 25 μl:2 μl (10 μM) of each forward and reverse primers, 12.5 μl of Promega GoTaq® Green Master Mix (Promega Corporation, Wisconsin, USA), 1 μl of DNA (20 ng/μl) and 7.5 μl H₂O.

All primers were synthesized following New Zealand Genomics Ltd. (NZGL) recommendations, with an overhang on the 5′ ends to permit the ligation with Illumina multiplexing indices and sequencing adapters (Table 1). The optimum annealing temperature with the overhangs was determined with a gradient PCR. For each primer, 2 μl of PCR product were checked on 1.5% agarose gels. Those that presented a clear band of expected size were cleaned up with Agencourt AMPure beads following the manufacturer's protocol (Beckman Coulter Life Science Inc., México). Afterwards, PCR products were quantified using a Qubit™ 3.0 fluorometer (Thermo Fisher Scientific Inc., Massachusetts, USA). Purified PCR products of the same individual with concentration higher than 1.0 μg/ml were pooled together. Library preparation with 96 Nextera Index Primers, quantification, normalization, and pooling were performed by New Zealand Genomics Ltd. (NZGL) in their laboratories. The library was then sequenced with MiSeq Reagent Kit V3 in MiSeq sequencer (Illumina Inc., USA).

MiSeq reporter was used to separate and remove the adapters for the 96 samples. The software Fastq-Multx (Aronesty, 2011) was used to demultiplex amplicons according to the primer sequences. Due to the wildcard characters in the primer sequences, up to 7 base pair mismatches were permitted. Software SolexaQA++ 3.1.4 was used to ensure the reads were still paired (Cox et al., 2010). The paired end reads (read 1 and read 2) were merged using VSEARCH 1.9.5 (Rognes et al., 2016). Paired reads that did not meet the following quality filtering were discarded (Edgar and Flyvbjerg, 2015): (i) reads with quality score over 3, (ii) reads longer than 140 bp, (iii) reads with less than one expected error in the primer sequence or barcodes. Unique sequences were clustered using a 97% identity threshold and remaining singleton Molecular Taxonomic Units (MOTUs) were discarded. Chimeric sequences were then removed using UCHIME (Edgar et al., 2011). Using the final list of representative sequences, each MOTU was searched against the GenBank database using BLAST 3.2.31 (Camacho et al., 2009).

MOTUs with BLAST query coverage under 60% or BLAST identities lower than 74% were also deleted from the database. Potential contamination and predator MOTUs (i.e., A. media and O. vulgaris) were removed from the database. Potential prey MOTUs were assigned using the following criteria to taxonomical categories: MOTUs with identity higher than 97% were determined at species level, MOTUs between 93 and 97% were assigned to genus, and MOTUs with identity below 93% were assigned to family.

For each predator (A. media and O. vulgaris), we analyzed separately the MOTUs identified by different pairs of primers. Then, we calculated the proportion of reads for each MOTU in relation to the total number of reads (PR). We also calculated the frequency of occurrence for each MOTU (FM: percentage of number of samples tested positive for a given MOTU in relation to the total number of samples) and the frequency of the occurrence of each family (FF: percentage of number of samples tested positive for a given prey family in relation to the total number of samples).

Frequency of occurrence was calculated for higher taxonomic levels by combining information for all MOTUs falling within the relevant taxon (Table 2). Moreover, for those taxa that were detected with at least two pair of primers and to species level, we calculated the overall frequency of occurrence (percentage of samples which tested positive for a given taxon in relation to the total number of samples, Table 2).

Redundancy Analysis (RDA) was used to detect patterns in the diet of O. vulgaris and A. media and determine which explanatory variables influenced those patterns. We included the occurrences of each family (FF) detected in the analysis as response variables, considering the presence-absence of each family in each paralarva. All octopus paralarvae presented three suckers per arm and were therefore probably less than 10 days old (Villanueva, 1995; Garrido et al., 2016b). All squid paralarvae were less than 41 days old based on statolith ring measurements (Olmos-Pérez et al., unpublished data). However, since we did not have complete age data, size-at-age is very dependent on environmental (pre- and post-hatching) temperature, and food ingestion is likely more dependent on size that on age, we categorized DML into three different classes to facilitate detection of ontogenetic changes during paralarvae growing. Thus, O. vulgaris paralarvae were categorized as small (1.20–1.74 μm; n = 21), medium (1.75–1.98 μm; n = 21), or large (1.99–2.28 μm; n = 22), and A. media were categorized as small (1.42–1.99 μm, n = 9), medium (2.00–2.99 μm; n = 15), or large (3.00–6.02 μm; n = 8). Paralarval size class (i.e., small, medium, large), transects (i.e., transects T2, T3, T4, T5), seasons (i.e., summer, early autumn), and depth (z1, z2, z3, z4, and z5) were included as nominal explanatory variables. We used the correlation triplot (α = 0, species conditional triplot), and the correlation matrix for the response variables. A significance test was applied with 4,999 permutations.

The effects of the season, transects, depth, and DML on dietary diversity and the presence of particular prey families in the diet were analyzed with generalized additive modeling (GAM) using a Poisson distribution for diversity and a binomial distribution for the other response variables and logit link function (link=logit), with season and transect as fixed factors and DML and depth effects fitted as smoothers (setting the bases dimension using k = 4 to avoid overfitting). Only the most frequently predated families (i.e., those detected in at least 10% of the predators, FF > 10) were used in this analysis. Models were fitted using backwards selection. The goodness-of-fit of the models was assessed with the Akaike Information Criterion (AIC). When the difference in AIC between two models (i.e., with and without one explanatory variable) was less than 2, an F-test was employed to select the best model (in case of a significant F-value the more complex model was preferred). All statistical analyses were performed with Brodgar 2.7.4. (Highland Statistics Ltd., UK).

Finally, “discovery curves” were plotted to determinate if the number of samples was sufficient to determine the importance of the most frequently detected prey species, for each combination of predator species (O. vulgaris and A. media), and primers (COI and 16sa). For each prey species, predator and primers, sets of 0, 1, 2… n samples were drawn at random from the available n samples and the proportional occurrence of the prey type was calculated for each sample size. Ten replicates were used to generate means and confidence limits, which were then plotted against sample size.

Of 8,274,658 raw reads, 5,734,163 were successfully demultiplexed and contained both read 1 and read 2. Then, 5,119,926 paired end reads were successfully merged, and 4,752,768 reads remained after quality filtering. A total of 4,124,464 reads was clustered into 1,155 MOTUs using a 97% threshold. Of the total reads, 3,189,247 corresponded to COI, 744,474 reads to the primers set 16Sa and 190,743 reads to the primer set 16Sb. Of these, 405 MOTUs (31,604 reads) did not match any GenBank sequence and 1,155 MOTUs, 750 (4,092,860 reads) had a match on GenBank database.

After meeting the thresholds of query coverage and identity, 465 MOTUs (131,843 sequences) were removed and 285 MOTUs (3,961,017 reads) were classified as: contamination (48 MOTUs, 41,422 reads), O. vulgaris (60 MOTUs, 2,150,967 reads), and A. media (55 MOTUs, 1,534,538 reads) and were removed prior to the analysis. In total, 122 MOTUs (234,090 reads) were then considered as potential prey.

Finally, 66 prey MOTUs (112,015 reads) were detected with the pair of primers COI, 53 prey MOTUs (122,029 reads) with pair of primers 16Sa and 3 prey MOTUs (46 reads) with pair of primers 16Sb (Supplementary Material 1). The total number of sequences in each category detected in both predators by different pair of primers was, the average and the range are presented in Supplementary Material 2.

Genetic analyses with primer COI revealed the presence of prey in the digestive glands of 56 O. vulgaris (56/64 = 87.5%) and 25 A. media (25/32 = 78.1%). Primers 16Sa revealed the presence of prey in the digestive tracts of 51 O. vulgaris (79.7%) and 10 A. media paralarvae (31.3%). Finally, 16Sb revealed the presence of prey in 2 O. vulgaris (3.12%) and 7 A. media paralarvae (21.87%). Primers COI amplified a wide spectrum of species detected in gut content of cephalopods paralarvae, belonging to a minimum of 7 phyla, 15 orders, and 32 families. Within these higher taxa, we were able to distinguish 28 genera and 27 species (Table 2). Primers 16Sa amplified mainly decapods, but also species belonging to other taxonomic groups. In total, they amplified taxa belonging to a minimum of 3 phyla, 6 orders, and 22 families. Within these, we were able to distinguish 21 genera and 24 species (Table 2). Primers 16Sb amplified in total 2 phyla, 3 orders, 3 families. Within these, we were able to distinguish 3 genera and 3 species (Table 2). In total, 21 families were exclusively identified by COI primers, 9 by 16Sa primers, and 2 by 16Sb. Thirteen families were detected with both primers COI and 16Sa (Table 2).

In total, considering all pair of primers together, prey were detected in 62 (96.9%) O. vulgaris and in 26 (81.3%) A. media paralarvae. The number of different prey taxa identified in individual O. vulgaris paralarvae range between 0 and 9 (mean ± standard error, 2.1 ± 0.267) with COI and between 0 and 7 (2.22 ± 0.232) with 16Sa. In individual A. media the number of prey taxa identified with COI primers were between 0 and 8 (2.09 ± 0.334) and between 0 and 11 (0.94 ± 0.378) with 16Sa.

The most abundant prey reads detected with primers COI (Figure 2A) matched with the crabs Goneplax rhomboides (order Decapoda), an unknown species of the family Portunidae (order Decapoda) and Pilumnus hirtellus (order Decapoda), as well as an unknown ophiuroid of the family Euryalidae (order Euryalida). With 16Sa primers (Figure 2B) the most abundant prey reads matched with the crabs Carcinus maenas (family Carcinidae) and P. hirtellus (family Pilumnidae) (Supplementary Material 1).

COI primers identified a total of 54 unique MOTUs belonging to 6 phyla, 14 orders, and 31 families. Within these, we were able to distinguish 28 genus and 20 species while 16Sa primers identified a total of 47 MOTUs belonging to 3 phyla, 5 orders, and 20 families. Within these, we were able to distinguish 18 genus and 19 species (Table 2).

Using COI primers, the most frequently detected MOTUs (FM) in O. vulgaris were the crab P. hirtellus (family Pilumnidae), the copepod Paracalanus sp. (family Paracalanidae) and the crab G. rhomboides (family Goneplacidae). With 16Sa, the most frequently detected MOTUs (FM) were the decapods C. maenas and P. hirtellus. The remaining MOTUs were detected in less than 10 octopus paralarvae (Supplementary Material 1).

Pilumnidae was the most frequently detected family in O. vulgaris with primers COI and 16Sa (FF = 47 and 67%, respectively) (Figure 3). With primers COI, other families detected in more than 10% of O. vulgaris paralarvae were Paracalanidae, Goneplacidae, Paguridae, Portunidae, Euryalidae, Sididae, and Polybiidae (Figure 3A). Primers 16Sa detected the families Carcinidae, Goneplacidae, Inachidae, Paguridae, and Polybiidae in more than the 10% of octopus paralarvae (Figure 3B).

Analysis with the primer 16Sb revealed that Rosacea flaccida (Order Siphonophorae) was present in two O. vulgaris paralarvae (four reads).

The most abundant prey reads with COI primers (Figure 2A) matched with the copepod Paracalanus sp. (order Calanoida), followed by the hydrozoan Obelia geniculata (order Leptothecata), and an unknown ophiuroid of the family Euryalidae (class Ophiuroidea). With 16Sa (Figure 2B), the most abundant prey reads matched with the crabs C. maenas (family Carcinidae), Pisidia longicornis (family Porcellanidae), and P. hirtellus (family Pilumnidae) (Supplementary Material 1).

In A. media, COI primers identified a total of 29 unique MOTUs, belonging to 5 phyla, 10 orders, and 17 families. Within these, we were able to distinguish 16 genera and 11 species, while 16Sa primers identified a total of 18 unique MOTUs, belonging to 2 phyla, 5 orders, and 13 families. Within these, we were able to distinguish 12 genera and 14 species (Supplementary Material 1).

The MOTUs most frequently detected (FM) in A. media were the hydroid O. geniculata (family Campanulariidae), the copepod Paracalanus sp. (family Paracalanidae), and the siphonophore Muggiaea sp. (family Diphyidae). The remaining 31 MOTUs were detected in less than 15% of squids (Table 2). 16Sa primers revealed that the most frequent detected MOTUs were P. longicornis and P. hirtellus. Seventeen MOTUs were detected in <10% of the squids (Supplementary Material 1).

The most frequently detected families detected with COI primers were: Campanulariidae (order Leptothecata), Paracalanidae, Clausocalanidae (order Calanoida), Diphyidae, Euphausiidae, and Euryalidae (Figure 4A). The most frequently detected families detected with 16Sa primers were Pilumnidae and Carcinidae (Figure 4B).

In A. media, 16Sb revealed the presence of the hydrozoan O. geniculata (FM = 12.5%, 32 reads), the siphonophore R. flaccida (FM = 9.3%, 7 reads), and the salp Thalia democratica (FM = 3%, three reads) (Supplementary Material 1).

RDA analysis showed that season and transect, but not depth or individual size, significantly affected the prey families detected in the diet of O. vulgaris (Table 3). The sum of all canonical eigenvalues was 0.192 and the first two axes accounted for 51.53% of the fitted variation (i.e., the 9.91% of the total variation in the family data; Figure 5). In A. media, transect, but not season, size, or depth, significantly affected the families detected in their diet (Table 3). The sum of all canonical eigenvalues was 0.309, and the first two axes accounted for 48.73% of the fitted variation (i.e., the 15.05% of the variation in the family data; Figure 5).

The GAM analysis for the most frequently occurring families (detected in at least 10% of paralarvae) in O. vulgaris revealed that the copepod family Paracalanidae (FF = 39%) was more frequently predated in autumn than in summer (p < 0.001; Table 4). Predation on the crab family Pilumnidae (FF = 47%) differed across transects (p = 0.028; Table 4), being more frequent in T5 and T4 than the rest of transects (Figure 6). The ophiuroid family Euryalidae (FF = 14%) was more frequently predated in autumn than in summer (p = 0.010, Table 4) and was more frequent in smaller individuals (p = 0.020; Figure 7). The cladoceran family Sididae (FF = 13%) was more frequently found in small individuals (p = 0.030, Figure 7) and only detected in autumn (Table 4). The decapod families Goneplacidae (FF = 23%), Portunidae (FF = 14%), Paguridae (FF = 14%), and Polybiidae (FF = 11%) did not differ between seasons, among sizes or transects (p > 0.05 in all cases). The number of families in the diet of O. vulgaris differed between seasons (p < 0.001; Table 3) and among transects, (p = 0.023; Table 3). Thus, a wider range of families was predated in autumn than in summer, and in T3 than in T2 or T5 (Table 3). DML or depth did not affect the number of prey families detected in O. vulgaris (p > 0.05 in both cases).

The GAM analysis of the most frequent families (detected in at least 10% of paralarvae) in A. media revealed that the occurrences of the families Campanulariidae (FF = 47%), Paracalanidae (FF = 31%), Clausocalanidae (FF = 19%), Diphyidae (FF = 16%), and Euryalidae, (FF = 13%) in the diet were not affected by size, depth, transect, or season (p > 0.005 in all cases). However, family Euphausiidae (FF = 13%) was only detected in autumn. The number of families in the diet did not differ significantly with paralarval size, depth, transect, or season (p > 0.05 in all variables).

The discovery curves for O. vulgaris (Figure 8) showed stabilization of the proportional occurrence estimates, when at least 45 of 64 paralarvae were sampled. The discovery curves for A. media (Figure 9), did not show any stabilization for the whole number of samples analyzed (n = 32).

Octopus vulgaris paralarvae mainly preyed on decapod species, that generally comprise <5% of the total zooplankton abundance in the Ría de Vigo (Roura et al., 2013; Buttay et al., 2015). Among decapods, the species most frequently detected in O. vulgaris were the crabs C. maenas, P. hirtellus, and G. rhomboides (families Carcinidae, Pilumnidae, and Goneplacidae, respectively), that are also the most abundant decapod species in the Iberian Peninsula coast (Paula, 1987; Fusté and Gili, 1991; Queiroga, 1996). Family Pilumnidae was less frequent in more oceanic transects (T5), probably because they migrate from estuarine zones to offshore waters during their larval development and there is higher concentration in more inshore waters. Moreover, species of this family were more frequently detected in paralarvae captured at depths between 5 and 10 m, probably because they migrate to the upper water layers at night (Dos Santos et al., 2008).

The second most frequently detected group in O. vulgaris gut contents were the Calanoid copepods, a group not detected in previous studies (Roura et al., 2012). In particular, Paracalanus sp. was the main copepod identified in O. vulgaris gut. In Galician zooplankton communities, Calanoid copepods in general represent more than 60% of total zooplankton abundance (Blanco-Bercial et al., 2006; Roura et al., 2013; Buttay et al., 2015). Zooplankton community studies in this area have also shown that high abundances of Paracalanus species are linked to low salinity values (Blanco-Bercial et al., 2006). In our study, this prey was more frequently detected in O vulgaris paralarvae captured in autumn. The upwelling conditions during this season (i.e., cold and low salinity waters), could have promoted high abundances of this species increasing their availability in the environment and thus facilitating the predation.

Brittle stars (family Euryalidae) and cladocerans (family Sididae) were both frequently detected in small O. vulgaris paralarvae, perhaps because they are an easier target than fast moving copepods and decapods. The cladoceran identified with COI primers was Penilia avirostris. This species has been highlighted as an indicator of warm waters, and high abundances have occasionally been described in the Ría de Vigo associated with an increase in water temperature (Figueiras et al., 2011). The sea surface warming trend observed in Galician coastal waters during recent years (Gómez-Gesteira et al., 2008) could be favoring the presence of this cladoceran species. Another cladoceran that is very abundant in the Ría de Vigo was detected by primers 16Sa, namely Podon intermedius (Roura et al., 2013; Buttay et al., 2015). It was identified also in small and medium individuals. The detection of abundant cladoceran species in octopus guts could suggest opportunistic predation on cladocerans, specifically by smaller paralarvae.

Only one fish species was identified with COI primers in a single O. vulgaris paralarvae, and no fish DNA was amplified with 16Sb primers that were specifically designed to amplify fish DNA (Deagle et al., 2009). This result suggests low predation on fish, perhaps because the high mobility of fish larvae makes it difficult for the paralarvae to capture them.

Regarding squids, in A. media different prey species and different frequencies of occurrence were detected compared to O. vulgaris: Cnidarians were detected in A. media paralarvae of all sizes. Cnidarians are not very abundant in zooplankton community in Galicia (Buttay et al., 2015). Thus, these results could suggest selective predation on cnidarians, as also observed in turtles and sunfish (Dodge et al., 2011; Sousa et al., 2016). In contrast, cnidarians were only detected in three O. vulgaris. Their rare presence might be explained as a secondary predation effect (Sheppard et al., 2005) because high resolution of NGS, can detect small DNA amount present in the digestive tract of a prey captured by the paralarvae. It is also possible that hydroids are predated by O. vulgaris because they are easy to capture for slow recently hatched paralarvae (<10 days old, Garrido et al., 2016b). Moreover, squid paralarvae ingested up to ten copepod species, while only four were detected in octopus. This difference between A. media and O. vulgaris might be related with their hunting skills, which are developed during initial life stages (Villanueva et al., 1997). Alloteuthis media also preyed on decapods, and species of this group were mainly detected with the primer pair 16Sa. Thus, this could imply that the amount of DNA present was low and it was only possible to amplify decapod DNA with the specific pair of primers.

Other prey detected in both cephalopod species such as amphipods, cladocerans, euphausiids, and fishes, had been previously detected in the paralarval digestive system, but with a lower taxonomical resolution (Passarella and Hopkins, 1991; Vecchione, 1991; Venter et al., 1999; Vidal and Haimovici, 1999; Roura et al., 2012). Additionally, the gut contents of paralarvae of both species included molluscs, echinoderms, chaetognaths, and a nemertean that had never been previously identified in cephalopod paralarvae. Finally, DNA of chaetognaths and nemerteans was detected in a small number of paralarvae gut contents, and thus could reflect opportunist predation or alternatively, their DNA might be present in an organism ingested by the paralarvae, as an effect of secondary predation as explained above.

Diet diversity for O. vulgaris was influenced by the season and distance to shore. Numerous studies have shown that zooplankton communities in Galicia change according to oceanographic and meteorological conditions (Bode et al., 2009; Roura et al., 2013; Buttay et al., 2015). Thus, diet variability observed in O. vulgaris paralarvae might be related to zooplankton changes in prey availability in the zooplankton community. In contrast, no relationship could be established between the diet of A. media and the environmental explanatory variables or individual size. This may be related to small number of samples analyzed: discovery curves in A. media, showed very wide C.I. and no stabilization of the proportional occurrence estimates for the whole number of samples analyzed (n = 32). In contrast, O. vulgaris discovery curves, showed narrower C.I. and a stabilization of the proportional occurrence estimates, when at least 45 paralarvae are sampled. These results suggest that the number of A. media paralarvae analyzed was insufficient for a comprehensive dietary analysis of this species. In contrast, results suggest that the number of paralarvae of O. vulgaris analyzed in this study could be enough for this dietary analysis.

Our results showed that O. vulgaris prey on a wide variety of decapod species, but also frequently prey on other taxonomic groups, including mollusks, ophiuroids, amphipods, cladocerans, copepods, chaetognaths, or cnidarians. However, the low number of samples analyzed in previous research could have prevented the identification of rarely detected prey, that would likely only be identified when increasing the number of paralarvae analyzed. Moreover, the employment of several primers targeting different genes, could have favored the detection of additional species with broader taxonomic range that previous studies.

Overall, our results showed the usefulness of the NGS approach with several primers targeting different genes to dietary analysis of wild cephalopod paralarvae. Results have shown that they feed on a wide diversity of prey, mainly decapods, copepods, and cladocerans, but also other taxa that have not been previously identified in wild cephalopod paralarvae such as mollusks, echinoderms, chaetognaths, salps, cnidarians, and a nemertean. This study provides essential data to elaborate more suitable diets for captive cephalopod paralarvae, with the aim of increasing their survival for economically sustainable farming. Further studies are needed, including use of a wider variety of prey, mainly copepods from the genus Paracalanus, Cladocerans, and different decapod species, to test the effect on the digestive gland performance, growth and survival of recently hatched paralarvae.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.