Small copepods could channel missing carbon through metazoan predation

Abstract Global ecosystem models are essential tools for predicting climate change impacts on marine systems. Modeled biogenic carbon fluxes in the ocean often match measured data poorly and part of this could be because small copepods (<2 mm) are modeled as unicellular feeders grazing on phytoplankton and microzooplankton. The most abundant copepods from a seasonal upwelling region of the Eastern North Atlantic were sorted, and a molecular method was applied to copepod gut contents to evaluate the extent of metazoan predation under two oceanographic conditions, a trophic pathway not accounted for in global models. Scaling up the results obtained herein, based on published field and laboratory estimates, suggests that small copepods could ingest 1.79–27.20 gigatons C/year globally. This ignored metazoan‐copepod link could increase current estimates of biogeochemical fluxes (remineralization, respiration, and the biological pump) and export to higher trophic levels by 15.6%–24.4%. It could also account for global discrepancies between measured daily ingestion and copepod metabolic demand/growth. The inclusion of metazoan predation into global models could provide a more realistic role of the copepods in the ocean and if these preliminary data hold true at larger sample sizes and scales, the implications would be substantial at the global scale.

rates by 30% compared with laboratory values to match the mesozooplankton observations (Hernández-León & Ikeda, 2005) or consider an "artificial" (i.e. largely below the optimal prey/predator ratio) carbon ingestion coming from particulated organic carbon and bacteria (Buitenhuis et al., 2006). For example, in the COBALT model used to predict the fate of anthropogenic C, microzooplankton ingestion by mesozooplankton such as copepods was increased to >97% of its standing stock to artificially match observed and modeled mesozooplankton production (Stock et al., 2014). Global models thus omit a significant carbon source that might be channeled through copepod predation but has not yet been accounted for.
It is well known that several large copepods species prey upon smaller zooplankton including copepods (Kleppel, 1993;Turner, 2004), but it is less known that small copepods (<2 mm in size) apparently also ingest other small metazoans. The ingestion of metazoans by small copepods has been confirmed by visual analysis of gut contents in small calanoid (Paffenhöfer & Lewis, 1990), cyclopoid (Turner, 2004), and poecilostomatoid copepods (You-Bong, Bong-Cheol, & Makoto, 1998). Recently, molecular methods have also confirmed the detection of metazoan prey in small copepods in laboratory (Durbin, Casas, Rynearson, & Smith, 2008) and field studies (e.g. Cleary, Durbin, Rynearson, & Bailey, 2015;Hu et al., 2014;Yi et al., 2017). This information about the diet of small copepods could be F I G U R E 1 The changing paradigm in global ecosystem models. (a) Global ecosystem models consider small copepods as unicellular feeders preying upon phytoplankton (dark green) and heterotrophic protists (light green). (b) Metazoan prey detected in copepods with DNA-based methods (Cop1-Cop 16, see Table 1 for species details) represented in a maximum-likelihood tree associated with their closest matches. Bootstrap values above 60% after 1,000 replications are shown.  important because they are typically modeled as unicellular feeders in global ecosystem models ( Figure 1a), feeding upon phytoplankton (Paffenhöfer, 1988), ciliates, and heterotrophic dinoflagellates (Calbet & Saiz, 2005;Saiz & Calbet, 2007, 2011, but ignoring the potential carnivorous input from metazoans. To make significant further progress in modeling mesozooplankton in global models, the missing carbon source associated with the food preference and grazing rate as a function of food quality in small copepods need to be quantified (Mitra, Castellani, et al., 2014;Rose et al., 2010). As this carnivorous predation in small copepods has been disregarded as an energy source in marine ecosystem models, here we: (a) test the extent of the metazoan pathway in the copepod community of a seasonal upwelling region under two oceanographic conditions; and (b) produce a preliminary estimate of how much carbon and nitrogen could be missing in the pelagic realm at a global scale if we ignore this trophic link.

| Zooplankton sampling
Zooplankton samples in this study are from long-term sampling in the CAIBEX project off the Ría de Vigo, NW Iberian Peninsula Samples collected at the surface within the Ría de Vigo on July 14th and September 30th were selected for the following reasons: (a) samples belonged to the coastal mesozooplankton community, which is the community with highest diversity and abundance of small copepods (Roura et al., 2013); and (b) these samples reflected two contrasting oceanographic conditions (high and low productivity), which enabled testing of whether metazoan predation in copepods was influenced by the productivity of the system. We defined The dominant copepod species were sorted from samples after the downwelling event (seven species) and after the upwelling event (10 species, Table 1). We selected one adult female from each of the dominant species with sizes ranging from 0.8 to 2.7 mm length, commonly found in the zooplankton communities of upwelling ecosystems worldwide (Peterson, 1998) were included as a positive control for the presence of metazoans in the diet since they are well-known carnivores (Paffenhöfer & Lewis, 1990;Turner, 2004).

| Molecular identification of metazoan prey in copepods
To avoid possible contaminants from the body surface, animals were washed with sterile distilled water, and appendages removed (mandibles, maxillules, maxillas, maxillipeds, and swimming legs). The digestive tract from each copepod species was dissected or suctioned with a pipette and the DNA extracted  Simon et al., 1994). This set of primers was previously used to identify 20 different metazoan prey in Octopus vulgaris, but failed to amplify copepod DNA (Roura, González, Redd, & Guerra, 2012).
A semi-nested PCR approach was employed to maximize the amount of prey DNA available. An initial PCR was carried out with the universal 16Sar and the primer 16Scrur to increase the proportion of prey DNA. The second PCR (semi-nested) was performed with primers 16Scruf and 16Scrur to preferentially amplify the prey DNA. PCR reactions were performed in a final volume of 25 µl, containing 2.5 µl 10X PCR reaction buffer, 0.5 µl dNTPs (Roche), 0.75 µl each primer (10 µM), and 0.025 U/µl Taq polymerase (Roche). The first PCR contained 50-100 ng of template DNA and the nested PCR contained 1 µl of a 1:10 dilution of the product of the first PCR.
Cycling conditions for the primers 16Sar-16Scrur consisted of an initial denaturation at 94ºC for 2 min followed by 39 cycles of denaturation at 94ºC for 30 s, annealing at 58ºC for 35 s, extension at 72ºC for 40 s, and a final step of 7 min at 72ºC. Cycling conditions for semi-nested PCR with primers 16Scruf-16Scrur consisted of an initial denaturation at 94ºC for 2 min followed by 33 cycles of denaturation at 94ºC for 30 s, annealing at 56ºC for 35 s, extension at 72ºC for 40 s, and a final step of 7 min at 72ºC. PCR amplifications were carried out in a TGradient thermocycler (Biometra), and visualized on a 1.5% agarose gel.
Semi-nested PCR products with the expected size (around 200 bp) were cloned using TOPO TA Cloning kit (Invitrogen) following the manufacturer's protocol. Five clones per copepod were sequenced using 200 ng of plasmid DNA and the universal primer T7. Sequences recovered from clone libraries were considered to be part of the same operational taxonomic unit (OTU) if there was <1% sequence divergence, allowing for intra-specific variation and Taq polymerase errors. OTUs were compared to sequences found in GenBank using the BLAST algorithm to identify metazoan prey.
A phylogenetic tree was constructed to assign unknown sequences to the highest taxonomic level and to verify the BLAST identifications ( Figure 1b). The tree contained all OTUs obtained from the copepods, together with the five closest matches of each OTU that were downloaded from GenBank. Sequences were aligned using MAFFT v5.7 with default settings (Katoh & Standley, 2013

| Estimating the metazoan-copepod link: how much carbon could be missing?
We first obtained an estimate of the weight specific ingestion rates (WSIR) of adult calanoid and cyclopoid copepods that were fed carnivorous and omnivorous diets. The dataset used included predation on heterotrophic dinoflagellates, ciliates, and metazoans ( Figure 2a,b) from field (n = 122) and laboratory (n = 37) studies (Saiz & Calbet, 2007), to obtain an average estimate of metazoan ingestion rates for small and large copepods. We assumed that 60% of the copepods at a given time are adults and older stages (Turner, 2004), and these are likely to feed on metazoans more than nauplii and younger stages. Using the mean weight of calanoids and cyclopoids as a proxy, we obtained an estimate of the ingestion rate for these two groups (Figure 2a,b). The greater variability in field studies, including experimental methods, food concentrations, and copepod assemblages present, results in greater variability in estimates for small copepods (5-10 µgC, Figure 2a). Since most laboratory studies are at near saturation conditions, the WSIR measured in the laboratory ( Figure 2b) were between 17 and four times higher than those obtained in the field for calanoids and cyclopoids, respectively.
We then multiplied the WSIR by the adult copepod and older copepodite standing stock in the first 100 m of the water column globally (Figure 2c). This value was estimated using the following assumptions: zooplankton standing stock biomass in the upper F I G U R E 2 Global estimates of the metazoan-copepod link. Weight specific ingestion rates (WSIR) of carnivorous and omnivorous adult calanoid (WSIR cal ) and cyclopoid (WSIR cyc ) copepods (Saiz & Calbet, 2007) were estimated using their average weight from field (a) and laboratory studies (b). Adult calanoid and cyclopoid copepod standing stocks in GtC were calculated using different calanoid-cyclopoid percentages obtained for coastal (NRS-IMOS database) and oceanic realms (COPEPOD database) (c). Estimated annual carbon ingestion through the metazoan-copepod link (d), which results from multiplying the WSIR obtained in a and b by the adult copepod standing stock shown in c 100 m = 0.31 gigatons of C (GtC) globally (Bogorov, Vinogradov, Voronina, Kanaeva, & Suetova, 1968;Moiseev, 1971). Biomass was corrected by increasing it by 1/3 since mesozooplankton biomass is underestimated with 200 μm nets (Gallienne & Robins, 2001). About 80% of the zooplankton biomass is copepods (Kiørboe, 1998) and ~60% of the copepods at a given time are adults and older copepodites (Turner, 2004). Given that calanoid and cyclopoid copepods displayed significantly different WSIR (Figure 2a,b), we then split the adult copepod biomass into calanoids and cyclopoids to better estimate their contribution to the total biomass. We used two datasets to obtain an approximation of the different percentages of calanoids and cyclopoids from coastal and oceanic environments ( poecilostomatoids. The small poecilostomatoid genera Corycaeus and Oncaea are well known to be carnivorous (Turner, 2004) and have been traditionally considered cyclopoids and have recently been placed once again in this order, we included them in our calculations as cyclopoids. Finally, the product of WSIR by the adult copepod biomass in the upper 100 m was multiplied by 365 days to estimate the amount of carbon ingested through the overlooked metazoan-copepod link annually (Figure 2d).

| Impact of the metazoan-copepod link in pelagic fluxes of C and N
To estimate the consequences that this overlooked input of carbon could have on the trophic web in the top 100 m (Figure 3), fluxes of C and N were calculated using the following data: standing stock of phytoplankton = 45 GtC (Falkowski, Barber, & Smetacek, 1998). Values of C and N coming from phytoplankton and ciliate grazing (Calbet & Saiz, 2005), but corrected considering that the zooplankton biomass Field-estimated contributions of the metazoan-copepod link are shown in blue and the consequent increase in copepod respiration, growth, indirect primary production, ingestion, fecal pellet production, biological pump, and export to higher trophic levels (HTL) is represented in bold. Copepod predation on phytoplankton and ciliates (Calbet & Saiz, 2005) were calculated under the standing stock assumptions of this study (see Section 2 for details). Copepod growth, remineralization, dissolved organic matter (DOM), and transfer to HTL estimated from the COBALT model (Stock et al., 2014) are shown in italics was underestimated by 1/3 (Gallienne & Robins, 2001) and 80% of that biomass is composed of copepods (Kiørboe, 1998). Key rates for copepods (Ikeda & Motoda, 1978): assimilation efficiency = 0.7, growth efficiency = 0.3, ingestion rate = 2.5 respiration, and growth rate = 0.75 respiration. Excreted NH 4 + and N released in debris and pellets ~20% and 40% of ingested N (Kiørboe, Mohlenberg, & Hamburger, 1985). Global mean C:N of zooplankton = 6.86 (Martiny et al., 2013), C:N of phytoplankton = 6.6 (Ho, 2003), and C:N of ciliates = 5.0 (Stoecker & Capuzzo, 1990). The contribution of copepods to higher trophic levels (fisheries, 8%), remineralization (8%), and dissolved organic matter (DOM, 3%-5% of copepod biomass) were obtained from the COBALT model (Stock et al., 2014). Finally, carbon sinking was 3% of egestion (Turner, 2015).

| Metazoan predation in copepods
Between one and four OTUs were detected within the digestive tracts of the copepod species investigated ( Table 1). Most of the OTUs (10 of 16) showed homologies >96% with sequences available on GenBank and were identified to species level. The taxonomic identity of the remaining OTUs (78%-93% homology) was assigned on the basis of their topographical position on the consensus tree using bootstrap methods (Figure 1b). In total, 16 OTUs belonging to three different phyla were detected: Crustacea (n = 13), Mollusca (n = 2) and Chordata (n = 1). A total of seven and 12 OTUs were detected in the copepods collected under downwelling (n = 7) and upwelling conditions (n = 10), respectively (

| Estimating the metazoan-copepod link: how much carbon could be missing?
Based on copepod diets that include predation on heterotrophic dinoflagellates, ciliates and metazoans and mean copepod mass of calanoids and cyclopoids, we estimated that the WSIR was 0.023 and 0.063 for field studies (n = 122) and 0.396 and 0.241 µgC prey/ µgC copepod per day for laboratory studies (n = 37), using the mean body mass of calanoids and cyclopoids, respectively (Figure 2a,b).
The estimated standing stock of calanoid/cyclopoid copepods in the first 100 m was 0.12/0.078 GtC for the coastal area and 0.181/0.013 GtC for the ocean (Figure 2c) based on the percentages obtained from the two large public datasets.
Assuming the WSIR obtained as an estimate of the metazoancopepod link, then multiplying the WSIR by the copepod standing stock and converting to an annual estimate (Figure 2d)

| Impact of the metazoan-copepod link in the pelagic fluxes of C and N
Estimates from field studies (1.79-2.80 GtC/year, Figure 2d) were used to examine the hypothetical impact of the metazoan-copepod link in terms of fluxes of C and N in the pelagic trophic food web (shown in blue in Figure 3). Within this preliminary global-scale estimate, the metazoan-copepod link would represent 4.0%-6.2% of the primary production. This pathway of carbon and nutrients could support an equivalent 15.6%-24.4% increment in copepod growth, pellet production, and carbon sinking to deeper layers through the biological pump, as well as more carbon exported to higher trophic levels. The amount of N excreted through ammonia or fecal pellets would increase the nutrients available to PP in the photic layer from 17 to 20.7. Another indirect consequence of this disregarded trophic link is the amount of PP that indirectly reaches copepods through the consumption of zooplankton ( Figure 3). Thus, assuming that their zooplankton prey mostly feed on phytoplankton through a single trophic level, then 18.9%-29.7% of PP would be assimilated by the zooplankton and ultimately transferred to copepods.

| D ISCUSS I ON
Copepods are the most abundant metazoan animals on Earth and are a key link in pelagic ecosystems, preying upon unicellular organisms, and serving as prey for higher trophic levels (Gattuso et al., 2015;Pauly et al., 2002). Our preliminary results support that this current paradigm of "herbivorous" copepods could be underestimating their true role in the marine ecosystem by missing a key energy source. The molecular technique applied here showed that all the copepods analyzed from a coastal upwelling systems under two contrasting oceanographic conditions (upwelling vs. downwelling) had metazoan prey within their guts. Using a meta-analysis of published copepod feeding rates, we estimated that such metazoan-copepod link at global scale could be 1.79-27.20 GtC/year. This source of carbon-currently missing in global ecosystem models-could have cascading effects in the pelagic environment, increasing current estimates of biogeochemical fluxes (remineralization, respiration, and biological pump among others) as well as increased export to higher trophic levels.
Despite the relatively low number of clones sequenced and the rapid digestion of prey in copepods (Kleppel, 1993), we were able to detect 1-4 OTUs in the digestive tract of all copepods analyzed, suggesting that ingesting metazoans is common for small copepods and Unfortunately, molecular methods cannot distinguish developmental stages of the prey ingested (eggs, nauplii, larvae, or adults) and this can have repercussions in terms of carbon ingested depending on the size of the prey. Considering that the mean size of the copepods analyzed in this study was close to 1.5 mm total length, it is likely that the majority of metazoan prey ingested would be early stages (e.g. nauplii, zoeae, veliger). Nonetheless, the physical limitations involved in swallowing prey larger than the mouth would involve sloppy feeding, which is the breakage and partial ingestion of the prey, releasing particulate, and dissolved carbon to the water column (Møller, 2005). Despite the partial ingestion of the metazoan prey, the amount of carbon incorporated into copepod biomass could range from one to three orders of magnitude higher than that of preying upon ciliates/flagellates and phytoplankton (Saiz, Griffell, Calbet, & Isari, 2014).
Despite the high sensitivity of molecular techniques that allow detection of small quantities of DNA, little is known of the origin of the DNA (King, Read, Traugott, & Symondson, 2008). We cannot be certain whether copepods directly prey upon metazoans or upon fecal pellets of other metazoans (secondary predation through scavenging). Field studies provide some support for predatory behaviour in some small copepods as crustacean mandibles, fish eggs, and fish larvae have been observed within their guts and fecal pellets (reviewed in Kleppel, 1993, Turner, 2004. Secondary predation is a potential source of error when tracking trophic links (Sheppard & Harwood, 2005) and may be significant in copepods owing to their feeding strategy. However, since cyclopoid copepods are known to feed on fecal pellets where DNA is highly degraded (Turner, 2004), and detecting secondary predation heavily decreases with increasing time since ingestion (Sheppard & Harwood, 2005), the likelihood of detecting prey DNA from a scat within a copepod would be expected to be low compared with that of a recently consumed prey item. It is thus more likely that the DNA detected within the digestive tract of the copepods originates from fresh metazoans (either predated or grazed) rather than from scavenging (fecal pellets, marine snow or molts).
Although it has been assumed that in productive environments the linear diatom-copepod-fish food chain predominates after an upwelling event (Kleppel, 1993), it is now evident that most primary production is consumed by microzooplankton and small metazoans , with diatoms contributing only up to 8% of the copepod diet globally (Mitra, Castellani, et al., 2014). We expected that metazoan predation by copepods would be higher under downwelling conditions, since less primary production was available and they would have to switch to carnivorous feeding (Kleppel, 1993). In a recent study, the relative frequency of metazoan prey-measured as the relative abundance of mesozooplankton or gelatinous plankton sequences-detected in the small Pseudocalanus species was higher under low chlorophyll concentrations (<2 μg/L, equivalent to the downwelling conditions of our study), than under medium (4-6 μg/L) or high chlorophyll concentrations (>20 μg/L, Cleary et al., 2015). However, in our preliminary study, there was no significant difference in the mean number of prey detected per copepod species under downwelling or upwelling conditions. This result might be explained due to the low number of OTUs analyzed per copepod, as well as the low number of copepods analyzed.
Our study suggests that trophic pathways in upwelling ecosystems involve complex, highly interactive webs (Figure 1c), instead of the conventional linear trophic chain. Moreover, the metazoancopepod link detected is bidirectional, not only affecting the mortality and export to higher trophic levels but also promoting PP by releasing nutrients through sloppy feeding and reducing grazing pressure (Mitra, Castellani, et al., 2014). The metazoan-copepod link ( Figure 2d) represents a 15.6%-24.4% increase in the ingestion of carbon by small copepods compared to the combined ingestion of phytoplankton and ciliates (8.67 and 2.82 GtC/year, respectively, Calbet & Saiz, 2005).
Feeding rates used in the calculations in Figure 2 included values for copepods consuming ciliates and flagellates-where copepods are most efficient-but also large metazoan prey that may represent just one capture per day. Despite marked differences in the ingestion rate of these organisms, we assumed a mean carnivorous ingestion rate for all copepods (small and large), as estimates would be unrealistically underestimated if we only considered ingestion rates obtained for large copepods feeding upon large metazoans. Optimal prey/predator ratios reported for suspension-feeding calanoid copepods (2%-6%) and cyclopoids (4%-10%) (Saiz et al., 2014) suggest that small copepods (<2 mm) might optimally capture prey ranging from <0.04 mm up to 0.2 mm, and even upon prey larger than themselves such as fish larvae (Calbet, Carlotti, & Gaudy, 2007).
Therefore, we have been conservative by including ingestion rates obtained for smaller copepods feeding upon unicellular organisms and small metazoans. This enabled us to obtain a mean feeding rate that includes the "carnivorous" repertoire of copepods ranging from small nauplii to large fish larvae.
Taking into account the difficulties of obtaining a mean feeding rate for copepods that display an broad variety of feeding behaviors (ambush predators vs. suspension feeders) and the capability to switch between them (Kiørboe, 1998) across a wide range of sizes (from <0.5 mm to >1 cm), the results obtained in this study suggest remarkable consistency between our field estimates (1.79-2.80 GtC/year, Figure 2d) and those obtained with the ecosystem model COBALT (Stock et al., 2014), where 2 GtC/year of small copepods (<2 mm) are ingested by the larger mesozooplankton fraction (>2 mm). Also notable is that the carbon ingestion estimated here for the top 100 m based on laboratory studies (24.16-27.20 GtC/year, The metazoan-copepod pathway appears to represent a significant input of carbon ( Figure 2d) that could account for the discrepancies usually found between quantified unicellular ingestion and metabolic demands of copepods, supporting production rates and survival when there is limited available food in terms of phytoplankton (Calbet, 2001;Hirst & Bunker, 2003;Mayor, Sommer, Cook, & Viant, 2015). It will also help to clarify discrepancies found in global models between the observed and modeled copepod biomass. If the copepod-metazoan link would be included in these models, 15.6%-24.4% more biomass could be available for copepod production ( Figure 3) and there would be no need to artificially increase the copepod grazing rates on bacteria (Buitenhuis et al., 2006) and phytoplankton (Hernández-León & Ikeda, 2005) or consider an "artificial" (i.e. largely below the optimal prey/predator ratio) carbon ingestion coming from microzooplankton (Stock et al., 2014).
The complex interaction between bottom-up and top-down effects currently hampers our capacity to predict how the inclusion of the metazoan-copepod link will affect climate change predictions. Nonetheless, our preliminary calculations show that a more realistic interpretation of copepod trophic interactions could uncover a carbon pathway that has either been overlooked or thought to be minor (Buitenhuis et al., 2006;Kiørboe, 1998;Ohman & Hirche, 2001), with cascading effects across the entire pelagic and benthic ecosystems increasing current estimates of biogeochemical fluxes (remineralization, respiration, and biological pump) and export to higher trophic levels by 15.6%-24.4% ( Figure 3). The metazoan-copepod link is roughly quantified here and requires further investigation. If further studies confirm its magnitude as suggested here, the metazoan-copepod link needs to be accounted for in global ecosystem to represent accurately the zooplankton component . This could ultimately enhance our capacity to predict the fate of increased anthropogenic emissions in the oceans likely impacting current predictions of climate change or fisheries yields.

ACK N OWLED G M ENTS
We are indebted to Mariana Cueto for her technical assistance, Isabel Teixeira for the phyto-and microzooplankton pictures (CESAM, University of Aveiro) and Alexandra Castro for some of the zooplankton pictures. We also thank the captain, crew, and technicians of R/V "Mytilus" (IIM, CSIC Vigo) for their assistance in collecting the zooplankton samples. This study was supported by the projects CAIBEX (Ministerio de Ciencia e Innovación CTM2007-66408-C02), LARECO (CTM2011-25929) and CALECO (CTM2015-69519-R).
The first author was supported by a "Fundación Pedro Barrié de

AUTH O R S' CO NTR I B UTI O N S
AR conceived and designed the experiment, performed the experiments, wrote the first draft of the manuscript; AR and AFG collected samples; AR and AJR analyzed the data; AR, JMS, AFG, and AJR contributed reagents/materials/analysis tools; and all authors contributed substantially to revisions.

DATA ACCE SS I B I LIT Y
Prey detected in copepods is available on GenBank under accession numbers FR851241-FR851253 and LN614707-LN614709.
Copepod data were obtained from the public databases NRS-IMOS and COPEPOD.