Trophic Structure of Neuston Across Tropical and Subtropical Oceanic Provinces Assessed With Stable Isotopes

Introduction

The neuston, one of the less described and known aquatic ecological groups, is paradoxically the closest to our sampling platforms as it inhabits the upper centimeters of the ocean. The term neuston was coined in 1972 (Hempel and Weikert, 1972) to define the pelagic organisms that occupy the vicinity of the surface layer, albeit often in a temporally-restricted and variable manner. Neuston occupies a delimited ecological niche and is generally grouped into three ecological categories: (a) euneuston: organisms with maximum abundance in the vicinity of the surface on which they reside day and night; (b) facultative neuston: organisms concentrating at the surface only during certain hours of the day, usually during darkness; and (c) pseudoneuston: organisms with maximum concentrations at deeper layers but reaching the surface layer at least during certain hours (Marshall and Burchardt, 2005). The neustonic community structure is conditioned by sunlight and an array of endogenous (organic matter, respiratory, photosynthetic, decompositional processes) and exogenous (atmospheric deposition, inorganic matter, winds, wave action, precipitation, UV radiation, oceanic currents, surface temperature) variables and processes affecting nutrient inputs and recycling (Marshall and Burchardt, 2005; Rawlinson et al., 2005; Rezai et al., 2019). Furthermore, the neuston provides a food source to the zooplankton migrating from deeper layers to the surface (Hempel and Weikert, 1972), as well as to seabirds roaming over the oceans (Cheng et al., 2010). For these reasons, the neustonic community is believed to play a critical role on the structure and function of marine food webs. Yet, research on neuston communities to date focused predominantly on geographically-limited regions of the ocean (Zaitsev, 1971; Hempel and Weikert, 1972; Holdway and Maddock, 1983; Ebberts and Wing, 1997; Rezai et al., 2019) or coastal areas (Brodeur, 1989; Le Fevre and Bourget, 1991; Padmavati and Goswami, 1996). Consequently, neuston complexity is still poorly understood as studies on the community structure and the taxonomical composition of organisms inhabiting this ecological niche remain few (Rezai et al., 2019), and global scale analyses are yet lacking.

The neustonic animals form a subset of the zooplankton community, which plays a pivotal role in the functioning of marine ecosystems. Zooplankton are partially responsible for the active energy flux between superficial and deep layers of the ocean (Turner, 2002; Jónasdóttir et al., 2015; Hernández-León et al., 2020). Zooplankton species composition, biomass, and secondary production influence a wide range of trophic levels in marine communities, as they constitute a link between primary production and secondary consumers (Litchman et al., 2013; Benedetti et al., 2016; de Oliveira Sodré and Bozelli, 2019). Copepods constitute the most abundant zooplankton taxon in terms of biomass and diversity worldwide (Kiørboe, 2011; Neumann-Leitão et al., 2018); therefore, changes in their community composition can thus impact the biogeochemical cycles (Bianchi and Mislan, 2016) and might be indicative of climate variability impacts on ecosystem functioning (Hooff and Peterson, 2006).

Historically, zooplankton assemblages research has focused mainly on taxonomic studies and those related to community structure (Pomerleau et al., 2015). However, recently, research has veered toward an alternative trait-based approach (Pomerleau et al., 2015; Benedetti et al., 2016; Campos et al., 2017), providing a perspective more focused on groups of species with analogous functional traits. This allows individuals to be classified into types characterized by the presence/absence of certain alleles of a gene, into size classes, ecological guilds, or functional groups (FGs; Tuomisto, 2010). Functional traits are phenotypes affecting organism fitness, growth, survival, and reproductive ability (Violle et al., 2007; de Oliveira Sodré and Bozelli, 2019). These are regulated by the expression of genes within species, and the expression of traits regulate, in turn, the species fitness under contrasting biotic and abiotic circumstances (Barton et al., 2013). Moreover, a specific functional trait can also develop from the interactions between other traits and environmental conditions (Kiørboe, 2011), leading to a given trait grouping being favored under certain conditions. Zooplankton traits can be classified in accordance to ecological functions – feeding, growth, reproduction, survival, and other characteristics such as morphology, physiology, behavior, or life history (Litchman et al., 2013; Hunt et al., 2015; Brun et al., 2016). Particularly, feeding strategies and trophic groups are relevant to ascertain feeding efficiency and associated predation risk (Brun et al., 2017). Additionally, they facilitate the understanding of ecosystem services associated with zooplankton, such as the distribution of fisheries or biogeochemical cycling (Prowe et al., 2019) while also allowing the positioning of zooplankton taxa in the food web (Benedetti et al., 2016, 2018).

Stable isotope analysis (SIA) has been widely used to explore the food web structure; to identify an organism’s trophic position; to quantify carbon, nitrogen, and energy fluxes; and to characterize trophic niches (Fry, 2006; Bouillon et al., 2011; Middelburg, 2014). Generally, stable isotopes undergo a predictable trophic enrichment between prey and consumer (Minagawa and Wada, 1984) and reflect the organism’s diet over a considerable period of time (Vander Zanden et al., 2015). Ratios of carbon isotopes (δ¹³C) generally have a low trophic enrichment and are commonly used to identify carbon sources (DeNiro and Epstein, 1978; Vander Zanden and Rasmussen, 2001; Post, 2002a). Nitrogen isotope ratios (δ¹⁵N) show progressive enrichment between prey and consumers (DeNiro and Epstein, 1978; Minagawa and Wada, 1984; Post, 2002b; McCutchan et al., 2003) and were thus employed to estimate trophic positions (Fry and Sherr, 1984; Minagawa and Wada, 1984; Peterson and Fry, 1987; Owens, 1988; Post, 2002b). Thus, by measuring the ratios of δ¹³C and δ¹⁵N, it is possible to infer the trophic structure of marine food webs (Fry, 2006).

Recently, SIA has been increasingly employed for the characterization of trophic niche (Layman et al., 2007a,b, 2012; Hunt et al., 2015). The trophic niche of a single species, community, or ecosystem is the aggregate of the interactions between its constituents and the ecosystem (Elton, 1927), hence representative of the characteristics of its habitat and trophic position (Leibold, 1995). The trophic niche can be inferred from the isotopic niche (e.g., the δ¹³C–δ¹⁵N bi-plot space) to uncover relevant aspects of trophic structure. For instance, niche-based quantitative metrics can be used to test ecological theory and trophic responses to anthropogenic impacts (Layman et al., 2007a; Schmidt et al., 2007), including the degree of overlap between distinct trophic niches (Jackson et al., 2011).

To our knowledge, the use of isotopic niche metrics on the study of the neustonic zooplankton community has never been attempted before. Here, we bridge this gap by exploring the variations in the size of the trophic niche of neustonic zooplankton across subtropical and tropical oceanic provinces. Our main hypothesis is that the trophic structure of the neustonic community is preserved in spite of differences in nutrient sources and productivity along oceanic provinces. More specifically, in this study we (1) characterize ratios of C and N for neuston across subtropical and tropical oceanic provinces, (2) determine trophic structure similarities across oceanic provinces, and (3) analyze the relationships between selected environmental variables and the trophic structure among oceanic provinces. We do so based on the samples collected along the Malaspina Circumnavigation Expedition, which circumnavigated the subtropical and tropical ocean in 2010–2011 (Duarte, 2015).

Materials and Methods

Sample Collection

The neuston samples were collected along the Malaspina 2010 Expedition, which circumnavigated the globe and was carried out between December 2010 and July 2011 across tropical, subtropical, and temperate regions of the Atlantic, Indian, and Pacific Oceans between 35° N and 40° S (Duarte, 2015). Sampling stations were distributed to characterize pelagic communities across regions of the open ocean in the northern and southern hemisphere (Duarte, 2015). In order to allow for the intercomparability of the observations and avoid adverse weather during sampling, the cruise was scheduled to visit most regions during their spring-summer months. The sampling locations were assigned to Longhurst Biogeochemical provinces (Longhurst, 2007). There, four biomes (Polar, Westerlies, Trades and Coastal) and a total of 56 provinces were identified, based on the characterization of primary production, mixed depth layer, nutrients availability, photic depth, algal biomass, Brunt-Väisälä frequency, and the Rossby radius of internal deformation. Specifically, the neuston samples reported here included 10 oceanic Longhurst provinces (Longhurst, 2007). Four provinces in the Atlantic [Northeast Atlantic subtropical gyral (NASE; n = 16); North Atlantic tropical gyral (NATR; n = 14); Western tropical Atlantic (WTRA; n = 15); South Atlantic gyral (SATL; n = 12)], four in the Pacific [South Pacific gyre (SPSG; n = 18); Pacific equatorial divergence (PEQD; n = 12); North Pacific tropical gyre (NPTG; n = 14); North Pacific equatorial counter current (PNEC; n = 15)], and two in the Indian Ocean [Indian South subtropical gyre (ISSG; n = 18); South subtropical convergence (SSTC; n = 16); Figure 1; Supplementary Table S1]. Between three and four stations were sampled within each province. Samples were collected twice a day at 12 pm and 4 am, by towing a neuston sampler with a mouth opening of 80 × 30 cm and a mesh size of 200 μm, at 2–3 knots for 10–15 min, which sampled the first 15 cm of the water at a distance of 5 m from the starboard of the vessel (González-Gordillo et al., 2012). The content of each sample was stored in 4% formaldehyde until analysis. Environmental data were acquired using CTD casts deployed at each station from the surface to 4,000 m or 100 m above the seafloor when this was shallower than 4,000 m.

Sample Processing, Taxonomical Identification, and Functional Grouping

Samples were stained with Bengal Rose in order to facilitate the identification of the organisms to taxonomic level of a family which, whenever possible, was determined under a stereomicroscope (Olympus SZX16) using appropriate guides (Boltovskoy, 1999; Castellani and Edwards, 2017). The most abundant and representative taxa selected for SIA were copepods of the families Acartiidae, Calanidae, Corycaeidae, Oncaeidae, and Pontellidae, and Phylum Chaetognatha. The latter were included as a representative of top planktonic predators. Each of these taxa have specific trophic roles, according to the literature. For instance, copepod families were assigned to FGs according to Benedetti et al. (2016, 2018). Four FGs were considered and taxa were assigned as follows: FG1 = detritivores (Oncaeidae); FG2 = herbivores/omnivores (Acartiidae and Calanidae); FG3 = carnivores (Corycaeidae and Pontellidae); and FG4 = predators (Chaetognatha). All specimens selected for analysis were adults according to their morphological aspect under the microscope.

Stable Isotope Analysis

Analysis of the natural abundance of carbon and nitrogen isotopes were performed on previously dried (50°C, 48 h) neuston samples. Even when specimens were classified at genus or species level, we pooled individuals from the same sample station, including day and night samples, at family (Acartiidae, Calanidae, Corycaeidae, Oncaeidae, and Pontellidae) or phylum (Chaetognatha) level in order to obtain ca. 1 mg dry weight for isotopic determination and to achieve a minimum of three data points for each functional group and minimize the constraints of small sample size (Jackson et al., 2011). Dried samples were packed into tin capsules and measured in an elemental analyzer (Carlo Erba CHNSO 1108) coupled to an isotope-ratio mass spectrometer (Finnigan Matt Delta Plus). Isotopic analyses were performed by the Servicio de Análisis Instrumental of the Universidade da Coruña (Spain). Samples were not acidified in order to remove carbonates since the selected taxa were only slightly calcified, and it has been demonstrated that this procedure could impact nitrogen measurements (Mateo et al., 2008). No corrections were made for the possible effect of formaldehyde on the stable isotope composition as most studies specifically made on marine zooplankton samples point out to no significant effects in δ¹⁵N and generally a decrease of less than 2‰ in δ¹³C after several years of storage (e.g., Mullin et al., 1984; Bicknell et al., 2011; de Lecea et al., 2011). Taking into consideration that all samples have been preserved similarly and that the time between collection and SIA largely exceeded 1 year, we assumed that differences in stable isotope values reflect genuine variations in the provinces and selected FGs.

Values of natural abundance of stable isotopes were expressed as δ¹³C and δ¹⁵N (‰) relative to Vienna Pee Dee Belemnite and atmospheric nitrogen, respectively (Coplen, 2011). Certified isotope standards (USGS40 and L-alanine) were analyzed along with internal acetanilide and sample standards with standard deviation (SD) between certified and measured values <0.1‰. The precision [standard error (SE)] of replicate determinations of standards and samples was <0.05‰ for both isotopes (n = 4). As the C:N mass ratio of most samples exceeded 3.5 (thus suggesting a significant and variable lipid content), δ¹³C values were normalized using an empirical linear regression with the sample C:N value determined for aquatic organisms (Post et al., 2007). This procedure aimed at removing the effect of the low δ¹³C associated to lipids and was preferred due to low biomass constraints that hindered the use of lipid removal methodologies such as dichloromethane and an accelerated solvent extraction system (Bodin et al., 2009).

Environmental Variables

Environmental attributes of each of the sampled stations were characterized by a number of variables collected in situ or satellite-derived. (Supplementary Table S2). These variables were already employed in a previous analysis of Malaspina 2010 cruise, where methodological details can be found (Mompeán et al., 2013, 2016b; Fernández-Castro et al., 2015). In brief, the stratification of the water column was represented by the depth of the mixing layer (MLD, m), the mean squared Brunt–Väisälä frequency (N², s⁻²), and the depth of the chlorophyll maximum (DCM, m), all estimated from vertical profiles of a CTD equipped with a fluorescence sensor (Fernández-Castro et al., 2015). Similarly, nutrient inputs from deep layers were estimated by diffusivity due to turbulence (K_T, m² s⁻¹), determined from vertical casts of a microstructure turbulence profiler (Fernández-Castro et al., 2015). Phytoplankton biomass was represented by surface and photic-zone integrated chlorophyll-a (Chla_s and Chla_i, mg m⁻³ and mg m⁻², respectively), determined from acetonic extracts of phytoplankton (Estrada et al., 2016). Microplankton (40–200 μm) characteristics were indicated by an abundance of the nitrogen-fixer Trichodesmium (Tricho, cells ml⁻¹), carbon biomass (C_40–200, mg C m⁻³), and natural abundance of nitrogen isotopes (δ¹⁵N_40–200, ‰) measured in samples collected by vertical tows of a plankton net between the surface and 200 m depth (Mompeán et al., 2013, 2016b). Satellite derived variables were considered to represent conditions prevailing over larger spatial and temporal scales than those considered during the specific sampling of each station and included annual averages of primary production (PP) for 2010 (mg C m⁻² d⁻¹) and mean monthly atmospheric dust deposition (MDU, g m⁻² month⁻¹). The former represented regional productivity and was derived from the data provided by the Ocean Productivity website (http://www.science.oregonstate.edu/ocean.productivity/index.php) in a grid of 0.17° × 0.17° including each station position. Dust deposition, a proxy for atmospheric inputs of key nutrients as Fe or P, was estimated from Aqua-MODIS Aerosol Optical Depth at 550 nm and Aerosol Small Mode Fraction data provided by the Giovanni online data system (NASA Goddard Earth Sciences) from a grid of 1° × 1° near each station (Mompeán et al., 2016b). In depth information on sampling and analytical methodology employed throughout the Malaspina 2010 cruise can be found in Moreno-Ostos (2012).

Statistical Analysis

Differences in δ¹³C and δ¹⁵N ratios of Oncaeidae among oceanic provinces were tested with ANOVA and Bonferroni post hoc test on Statistica 12 (StatSoft, Inc., Tulsa, OK, USA). Analysis of the differences in stable isotope values among oceanic provinces was assessed using a permutational multivariate analysis of variance (PERMANOVA) and pairwise tests on PRIMER 6.0 (Clarke and Gorley, 2006) and the add-on package PERMANOVA+ (Anderson et al., 2008; Supplementary Table S3). For subsequent analysis, a representative primary consumer was used to normalize stable isotope values according to Clark and Fritz (1997) and Stasko et al. (2018) to account for spatial heterogeneity in carbon and nitrogen ratios. To this end, consumer δ¹³C values were normalized relative to a pelagic baseline (Δ¹³C_pel) as: ${\Delta }^{13}Cpel=\left(\frac{1000+\delta w}{1000+\delta c}-1\right)\times {10}^3$ , were δ_C is the consumer δ¹³C value and δ_W is the province-specific mean δ¹³C value of the copepod family Oncaeidae, considered as the reference baseline. Consumer δ¹⁵N values were normalized by subtracting the province-specific mean δ¹⁵N value of the copepod family Oncaeidae from the consumer δ¹⁵N, represented thereafter as δ¹⁵N_n. The use of Oncaeidae as baseline is not intended to reflect the base of the food web but rather to set a homogeneous starting point along the δ¹³C and δ¹⁵N continuum. All statistical analyses were from this point on performed on Oncaeidae-normalized values.

Analysis of the isotopic niche size were made in the isotopic space defined by Δ¹³C_pel–δ¹⁵N_n of each province/functional group combination with the package SIBER (Jackson et al., 2011). Bayesian estimates of the standard ellipse area (SEAb) were calculated for each combination to characterize the full variability in foraging habits and used resources (Layman et al., 2007a, 2012). These estimates accounted for uncertainty due to the number of samples. A minimum of three samples was required for these estimations. SEAb statistics were computed from 10⁴ simulations per province/functional group arrangement, although only the maximum likelihood ellipses were used for graphical representations. Pairwise comparisons of SEAb (Supplementary Table S4) were performed by calculating the proportion of ellipse size that differed between two given combinations being interpreted as a direct proxy for the probability that one combination is different from the other (Jackson et al., 2011). Because of the limitations imposed by the available sample biomass for isotopic analyses, data of some FGs could not be obtained for all provinces.

Principal component analysis (PCA) was used to investigate the relationships between environmental variables and the niche size (as a proxy for trophic structure) of the FGs among oceanic provinces using PRIMER 6.0 (Clarke and Gorley, 2006). First, normalized environmental variables measured on each sampling site were employed to compute the eigenvalues, eigenvectors, and principal component (PC) coordinates. Then, Pearson’s correlation coefficient between SEAb average and 95% credible intervals (CI) range for each functional group/province combination and the environmental PC coordinates were computed using the corr() function on R 3.6.3 (R Core Team, 2020). Finally, the correlation coefficients between the niche metrics and the PCA axes were plotted on the variable space of the PCA first two principal component coordinates axis.

Results

δ¹³C and δ¹⁵N Ratios

The mean δ¹³C for Oncaeidae varied significantly among provinces (Figure 2A). Bonferroni post hoc test revealed three distinct groups (a, b, and c) of decreasing mean values (−21.2, −21.8, and −22.3‰, respectively). The highest values were found in NATR and SPSG provinces, while the lowest values occurred in NASE. Similarly, significant differences in mean values of δ¹⁵N between provinces were also found (Figure 2B). All δ¹⁵N values were positive, with maximum values in SSTC province, followed by those in PNEC and NPTG, and minimum values (<2‰) in NATR. Bonferroni post hoc tests also revealed three groups (d, e, and f) of decreasing mean δ¹⁵N (10.1, 5.4, and 3.5 ‰, respectively).

Isotopic Niche Size of Functional Groups

The analysis of the isotopic niches revealed differences between provinces and FGs in all oceans. The maximum likelihood ellipses were distinct from each other in all provinces and FGs but showed large overlaps between them in some provinces/functional group arrangements.

Atlantic Provinces

There was a clear separation of the ellipses of top predators (FG4) and detritivores (FG1) in the provinces where these groups were analyzed (Figure 3). In contrast, ellipses for herbivores/omnivores (FG2) and carnivores (FG3) showed large overlaps. The variability of SEAb (Bayesian standard ellipse area) was highest for FG3 and lowest for FG2 (Table 1).

Pairwise comparisons of the isotopic niche areas using SEAb revealed that only FG3 varied significantly from small values in NASE, intermediate in NATR and high values in SATL, while no significant differences resulted when comparing other FG among Atlantic provinces (Supplementary Table S4).

Indian Provinces

Distinct maximum likelihood ellipses were present in both province/FGs arrays; nonetheless, significant overlaps are discernible throughout the isospace (Figure 4; Table 1). SEAb was consistently larger for all FGs in the ISSG province and the largest trophic variability was found for FG2. Moreover, SEAb were the smallest for both FG1 and FG3 in SSTC.

There were no significant differences in SEAb of equivalent FG among Indian provinces (Supplementary Table S4).

Pacific Provinces

The isospace area of the Pacific provinces follows the trend described for Atlantic and Indian provinces, with perceptible differences in the distinct provinces/FGs arrangements, but the large overlap was found in many cases, especially in PNEC and SPSG (Figure 5). The largest values of SEAb were found for FG4 in SPSG and the minimum values for FG2 in PEQD (Table 1). NPTG and PEQD presented narrower SEAb compared to SPSG and PNEC. Pairwise comparisons of SEAb indicated that the values of FG1 were significantly smaller for PNEC compared to those for SPSG and that those for FG2 in NPTG and PEQD were smaller than those in PNEC and SPSG, while those for FG3 in PNEC and SPSG were larger than their equivalent areas in NPTG and PEQD. Finally, SEAb of FG4 in PNEC and SPSG were larger than their comparable group in NPTG (Supplementary Table S4).

Environmental Variables and Trophic Structure

The PCA on the environmental variables revealed distinctive conditions among provinces, as stations are clearly clustered within each province (Figure 6). The first two principal components of the PCA accounted for 66.4% of the total variance (Table 2). Surface chlorophyll had the largest positive loadings with the first component (PC1), followed by the microplanktonic carbon biomass and turbulence diffusivity (Figure 7). Conversely, PC1 was negatively correlated with the depth of chlorophyll maximum. The second component (PC2) was positively correlated with primary production, δ¹⁵N_40–200, and the mixed layer depth, while showing negative correlations with atmospheric dust deposition, Trichodesmium abundance, and photic-zone integrated Chl-a. The remaining environmental variables have a minor relevance in explaining the reported differences (Figure 7).

The size of the isotopic niche (indicated by median SEAb values and 95% CI) of FG4 was positively correlated with PC1, while those of FG3 showed negative correlations with this component. In turn, SEAb of FG1 and FG2 were mainly negatively correlated with PC2 (Figure 7).

Discussion

In the current study, we provide the first comparison of the trophic structure of the neustonic community in tropical and subtropical provinces of the Atlantic, Indian, and Pacific Oceans (Longhurst, 2007; Duarte, 2015). The analysis of stable carbon and nitrogen isotopes of neuston FGs allowed for the characterization of carbon and nitrogen ratios, the quantification of the size of the isotopic niche as a proxy of trophic structure, and its relationships with environmental variables across distinct oceanic provinces. Notwithstanding the limitations imposed by the analytical requirement and, thus, the small number of samples employed, as well as the samples preservation method, the results of this study allow us to establish an approximation of the trophic interactions in the neuston and their role in pelagic food webs at a global scale.

δ¹³C and δ¹⁵N Ratios

The differences found in carbon and nitrogen appear to reflect variations in the environmental regimes of the distinct provinces. Provided that different phytoplankton taxa favor the use of distinct carbon sources for primary production, it is plausible that it could lead to variations in δ¹³C that can be passed upon the food web. For instance, heavier δ¹³C values have been associated with higher phytoplankton fractionation of inorganic carbon, and subsequent sinking of organic matter relatively depleted in ¹³C from the surface waters led to high values in the upper ocean at subtropical and tropical waters (Gruber et al., 1999; Schmittner et al., 2013). Moreover, uptake of atmospheric CO₂ via air-sea exchange and consequent fractionation in the conversion process to HCO⁻ ₃ can result in high δ¹³C-dissolved inorganic carbon (Zhang et al., 1995). Such processes could explain the enrichment in ¹³C (i.e., higher δ¹³C) found for the isotopic baseline (Oncaeidae) in the SPSG and NATR. Additionally, the enrichment in δ¹³C measured in SPSG could also be related to the presence of diatoms, as reported in the region during the same cruise (Estrada et al., 2016). Although they were not the dominant phytoplankton taxa (Estrada et al., 2016), they may have disproportionately contributed to the metazoan food web, while the dominant picocyanobacteria would contribute to the microbial food web. While cyanobacteria uptake of inorganic carbon occurs via direct HCO⁻ ₃ transport, diatoms are also capable of transporting CO₂ derived from the catalyzed dehydration of HCO⁻ ₃ (Tortell and Morel, 2002). This implies a larger isotopic fractionation (i.e., less negative values of δ¹³C) in diatoms when compared with cyanobacteria or other microalgae (Fry and Wainright, 1991). Another hypothesis for the δ¹³C enrichment measured in NATR province is related to the abundance of Sargassum spp. on the westernmost stations (Gouvêa et al., 2020). This macroalgae has a δ¹³C ranging from −16‰ to −18‰, and it is conceivable that it could have influenced the δ¹³C baseline in NATR (Cabanillas-Teran et al., 2019). However, differences in δ¹³C could be the result of many factors associated with changes in phytoplankton CO₂ fixation, which is known to vary in relation to temperature, concentration of aqueous CO₂, phytoplankton composition, and the availability of dissolved inorganic carbon and nutrients (Wong and Sackett, 1978; Descolas-Gros and Fontugne, 1990; Francois et al., 1993; Burkhardt et al., 1999; Popp et al., 1999).

Both continental organic matter and atmospheric CO₂ are characterized by having lower values of δ¹³C than their equivalents in oceanic (Perry et al., 1999) or upwelled waters (Gruber et al., 1999). Our results agree with the report of higher δ¹³C values in zooplankton from oligotrophic regions of the subtropical North Atlantic during the same cruise (Mompeán et al., 2013). In contrast, more negative δ¹³C values of the Oncaeidae baseline, such as those measured in NPTG and NASE, can be explained by the dominance of dinoflagellates and coccolithophores in these provinces (Estrada et al., 2016) and by the influence of upwelling in NASE (Gruber et al., 1999; Mompeán et al., 2013).

Ratios of nitrogen also varied between provinces, according to the differences found in δ¹⁵N values of Oncaeidae. Low δ¹⁵N values could be due to diazotrophy since the atmospheric nitrogen presents δ¹⁵N = 0‰, and it can be traced through the food web (McClelland et al., 2003; Mompeán et al., 2016a; Bode and Hernández-León, 2018). This process is mediated by specialized prokaryotes – diazotrophs – which thrive in usually nitrate-poor, warm and stratified waters such as the ones of the subtropical and tropical gyres (Falkowski, 1997; Zehr et al., 2003; Capone et al., 2005; Luo et al., 2012) and introduce into the ocean bioavailable nitrogen depleted in δ¹⁵N (Somes et al., 2010). Additionally, it is also known that isotopic fractionation occurs during phytoplankton uptake, and this may cause low δ¹⁵N values (<5‰) in plankton when dissolved nitrogen concentrations are high, as reported at the initial phases of blooms (Waser et al., 2000). In turn, high concentrations of ¹⁵N-enriched nitrate of marine origin are similarly reflected in the δ¹⁵N values of zooplankton (Owens, 1988; Montoya et al., 2002), as found in the SSTC samples. Our results strongly suggest that the low δ¹⁵N values measured for Oncaeidae reflect the role of N₂ fixation by diazotrophs in supplying N in some provinces (mainly in NATR but also in NASE, WTRA, SATL, SPSG, and PEQD), while the high δ¹⁵N values measured in other provinces indicate inputs of nitrogen derived from denitrification processes (e.g., PNEC, NPTG, and SSTC). These assumptions are supported by the consistency with estimations of the abundance of the N-fixer Trichodesmium (Estrada et al., 2016; Mompeán et al., 2016b) and N fixation rates and NO₃ ⁻ diffusive fluxes obtained for the Malaspina cruise (Fernández-Castro et al., 2015), as well as recent estimates of marine nitrogen fixation and denitrification at the scale of the global ocean (Knapp et al., 2016; Bonnet et al., 2017; Gruber, 2019; Wang et al., 2019).

Moreover, the δ¹⁵N values reported here for Oncaeidae are consistent with those of microplankton collected simultaneously during the Malaspina cruise (δ¹⁵N_40–200) in the upper 200 m of the water column, thus supporting our assumption of the role of Oncaeidae as primary consumers and critical nodes in the food web. Furthermore, the general coherence between the distribution pattern of both Oncaeidae and microplankton δ¹⁵N across provinces (e.g., Mompeán et al., 2016b; Supplementary Figure S1) suggests similar turnover times of the stable isotopes (and hence growth rates) in epipelagic microplankton and, in our case, mesoneuston (>200 μm).

Spatial Differences in Isotopic Niche

The small layer of surface water the neustonic community inhabits would imply a scarce availability of food sources, particularly in more oligotrophic regions, such as most of provinces in this study. This suggests a strong competition for resources and a large overlap in the isotopic niche for the different FGs. Our results show a general overlap between the niches of the four FGs regardless of the province. More specifically, herbivores/omnivores and carnivores had large overlap between them and with other groups in most provinces, while chaetognaths and detritivores had small overlaps.

The large overlap in the isotopic niches in all provinces supports the hypothesis of a common trophic structure of the neustonic community, notwithstanding regional and local differences in nutrient sources and productivity. Such structure would be facilitated by the large plasticity of the isotopic (trophic) niche of the main FGs. Whenever present, chaetognaths occupied the highest trophic position as expected from their obligate carnivorous feeding (Pearre and Pearre, 1982; Kehayias et al., 1996; Baier and Purcell, 1997). Conversely, detritivores generally occupied the lowest end of the δ¹⁵N_n spectrum. Despite the low trophic fractionation generally associated to δ¹³C (e.g., Post, 2002b), the Oncaeidae baseline generally shows the lowest δ¹³C values, with a consistent increase in δ¹³C from FG1 to FG4 in most provinces. This is, in turn, reflected by a decrease in Δ¹³C_pel from FG1 to FG4 due to our baseline normalization. The intermediate FGs (FG2 and FG3) were particularly variable in δ¹⁵N_n, irrespective of the pre-assigned trophic category as herbivore/omnivore or carnivore, thus limiting the applicability of these classifications in the neuston (Benedetti et al., 2016, 2018). This is illustrated in the Atlantic provinces where both productive and oligotrophic regions are present. Here, the more productive regions (NASE and WTRA) show a defined hierarchical trophic structure, contrasting with the oligotrophic provinces of NATR and SATL.

Our results show the importance of omnivory in pelagic food webs, particularly in oligotrophic regions of the oceans, such as the majority of the provinces sampled in this study. This plasticity in feeding preferences increases trophic connectivity and ultimately the ecosystem resilience and stability (Madigan et al., 2012). Additionally, organisms of intermediate trophic positions generally have varied diets, hence expanding their options to feed on organisms from various trophic positions (Svanbäck et al., 2015). Therefore, a large degree of omnivory combined with opportunistic feeding would significantly increase the size of the isotopic niche, as observed in the studied FGs of neuston across the global ocean. This strategy may facilitate their survival in the oligotrophic environment of the sub-gyres and tropical gyres of the ocean and is further supported by the large niche overlap observed as a general pattern in our results.

However, the degree of niche overlap between FGs on each province was variable. For example, provinces such as NASE and SSTC presented a smaller degree of overlap compared to provinces such as WTRA and SPSG. The degree of niche overlap appears to be related with the food availability on each province, since the former provinces are in the proximity of coastal areas and under the influence of nutrient-rich upwelling waters, while the latter provinces are in oligotrophic regions of the ocean. However, large niche overlap is not a direct indication of competition for the same sources of food, as organisms may be preying on different resources that share similar isotopic signals (Layman et al., 2012). For example, Bode et al. (2015) found an important contribution of phytoplankton for some copepod species, while others showed relevant influence of bacteria, probably derived from preying upon ciliates and bacterivorous plankton. Furthermore, a study encompassing the same provinces in the Atlantic Ocean as our study observed that several of the most abundant mesopelagic fish presented comparable trophic positions but differed in their diets, the latter related to local prey availability (Olivar et al., 2019).

Environmental Variables and Trophic Structure

The results of the PCA on the environmental variables confirmed the expected singularity of samples within the biogeochemical divisions as proposed by Longhurst (2007). Furthermore, our results showed that the structure of FGs was indeed influenced by the environmental variables considered, as changes in some variables were more relevant to the size of the isotopic niches of certain FGs. The negative correlation of the size of the niche of FG1, and to a lesser extent FG2, with PP is symptomatic of trophic specialization, meaning that as PP increases the size of the niche of herbivores and detritivores decreases, since these copepods would be able to select the preferred feeding sources. Alternatively, the niche width would increase in provinces with relatively large amounts of dust deposition and Trichodesmium, an indication of diazotrophic zones (Fernández-Castro et al., 2015; Mompeán et al., 2016a), characterized by nutrient-depleted waters and lower abundance of prey (Carlotti et al., 2018). Transfer of diazotrophic nitrogen to zooplankton (including neuston) would occur not only by direct grazing of the cyanobacteria but also by indirectly preying upon organisms of the microbial loop (Sommer et al., 2006; Wannicke et al., 2013). The need for including more diverse prey would thus explain the increase in the niche size.

The size of the isotopic niche of carnivores was positively correlated with DCM depth and negatively with microplankton biomass and surface Chl-a. It has been reported that the optimum prey size of copepods ranges between 10 and 100 μm (Hansen et al., 1994; Li et al., 2007); therefore, increases in abundance of microplankton would allow for a prey selection better suited to their dietary preferences, thus reducing their isotopic niche size. In turn, the positive correlation of the niche size of chaetognaths with Trichodesmium, stratification, and dust deposition suggests an increase in the variety of their prey in diazotrophic provinces, as suggested for herbivores and omnivores.

Our results suggest a relationship between the size of the isotopic niche and their overlap with other FGs. In this way, the FGs of provinces closer to coastal areas or influenced by upwelling (e.g., NASE, SSTC, and NPTG) exhibit relatively small niches, suggesting that the abundance and diversity of available prey allows for a more specialized feeding according to the theoretical functional group of each organism. This is supported by the high importance of diffusive nitrate inputs described for Malaspina stations in these provinces, particularly in SSTC (Fernández-Castro et al., 2015). In contrast, highly oligotrophic provinces (e.g., SATL, ISSG, and SPSG) present generally broader isotopic niches.

The occurrence of opportunistic feeding on available prey would trigger variations in the isotopic signal of both prey and consumers throughout the food web, therefore resulting in larger isotopic niches in provinces characterized by low primary productivity and low plankton biomass. Similar observations have been reported for both the Atlantic Ocean and the Eastern South Pacific, where zooplankton communities of oligotrophic regions depicted larger niches compared to those of more productive regions (Bode and Hernández-León, 2018; González et al., 2019). An exception to this is the PNEC province, where two of the three stations sampled in this study were under the influence of the Costa Rica Dome, an enhanced biological productivity region, where phytoplankton and zooplankton biomass are larger than in neighboring tropical waters (Fiedler, 2002). Nonetheless, while being a highly productive system its phytoplankton community is dominated by the picophytoplankton Synechococcus (Taylor et al., 2015), and it has been reported that zooplankton are still able to match their metabolic needs by protistivory and detritivory (Stukel et al., 2018), thus increasing the variability of individual samples and consequently the estimated niche size. Comparable results have also been found in the Indian Ocean where the niche sizes of an oligotrophic seamount were smaller than those of a more productive one (Annasawmy et al., 2020). Furthermore, PNEC is affected by distinct oceanographic processes (upwelling, minimum oxygen zone), the confluence of different water masses (Lavin and Marinone, 2003; Rau et al., 2003; Fiedler and Talley, 2006), and nitrification-denitrification processes (Gruber, 2019).

Conclusion

Values of C and N isotopic ratios for the neuston are dependent on the environmental characteristics of the oceanographic regimes of each province. Uptake of atmospheric carbon were found to be associated with high δ¹³C values, whereas upwelling influenced regions or the presence of dinoflagellates and coccolithophores appear to be related with low δ¹³C. Additionally, high δ¹⁵N is a consequence of nitrogen uptake derived from denitrification processes and NO₃ ⁻ diffusive fluxes, while nitrogen fixation by diazotrophs resulted in low δ¹⁵N.

The hypothesis of a common trophic structure of the neustonic community is supported by the general niche overlap between FGs regardless of the province, although chaetognaths and detritivores generally occupied the highest and lowest trophic positions, respectively. The importance of omnivory in oligotrophic regions (the majority of provinces sampled in this study) is consolidated with our results since omnivores and carnivores exhibit varying trophic positions regardless of their nominal trophic category. Environmental conditions were depicted to influence the trophic structure of the FGs, as the size and overlap of isotopic niches were related to changes in some of the variables. Typically, FGs present smaller niche size and overlap as favorable conditions increase whereas disadvantageous conditions promote the opposite. The use of isotopic niche metrics allowed the comparison of the trophic structure of the neustonic community across the tropical and subtropical ocean, which is paramount for understanding neuston interactions and its role in pelagic food web. Further studies exploring additional FGs and/or taxonomical diversity, evaluating temporal variability and impacts on the trophic structure of the community within each province, as well as estimates of neustonic biomass, are required to better understand and improve the knowledge of the trophic structure of the neustonic community in the pelagic ocean.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

CD designed, coordinated, and led the Malaspina Expedition. All authors designed the research reported. RA, AB, and JG-G performed research. RA and AB analyzed data. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by project Malaspina-2010 (CSD2008-00077) funded by program CONSOLIDER-INGENIO 2010 (Ministerio de Ciencia e Innovación, Spain), by grant IN607A 2018/2 of the Axencia Galega de Innovación (GAIN, Xunta de Galicia, Spain). Thanks are also due to FCT/MCTES for the financial support to CESAM (UIDP/50017/2020+UIDB/50017/2020), through national funds. RA was supported by a Ph.D. fellowship funded by FCT (PD/BD/113483/2015).

Conflict of Interest

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.

Acknowledgments

We are grateful to all participants in the Malaspina expedition for their collaboration in plankton sampling and to A.F. Lamas for lab assistance with sample preparation. We also acknowledge the reviewers for their constructive comments and inputs which helped to improve this manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2020.606088/full#supplementary-material

References