Zooplankton Abundance Reflects Oxygen Concentration and Dissolved Organic Matter in a Seasonally Hypoxic Estuary

<jats:p>Ocean deoxygenation, warming, and acidification resulting from global change and increasing nutrient inputs are major threats to marine ecosystems. Despite this, spatial and temporal patterns of oxygen availability and their impacts on marine life are understudied compared to warming and acidification, particularly in coastal tropical ecosystems. We measured the abundance of major groups of zooplankton in the context of five covarying environmental parameters [temperature, salinity, dissolved oxygen (DO), dissolved organic matter (DOM), and chlorophyll concentration] in a tropical estuary for one year. Partial least squares demonstrated that environmental conditions explained 20% of the variation in the community and found that temperature, salinity, DO, and DOM were most important (Variable Importance in Projection [VIP] &gt; 0.8). A generalized linear model identified depth, DO, salinity, and chlorophyll as significant main effects, and temperature and DOM were also significant via two-way interactions (p &lt; 0.05). When examined separately, the abundance of each zooplankton group was explained by a slightly different combination of environmental factors, but in all cases DO had large, significant effects, and in most cases DOM or its interactions were also significant. These results demonstrate that the seasonal cycle of hypoxia in this system significantly impacts the abundance of major zooplankton groups and likely also recruitment of benthic fauna through impacts on meroplankton and benthic-pelagic food webs.</jats:p>

section_content: Sufficient oxygen is necessary for animal life, yet oxygen concentrations in oceanic and coastal waters are falling [1, 2], thus exposing animals to more frequent or more severe hypoxia. Chronic marine hypoxia is already evident in large expanses of the open ocean in oxygen minimum zones, which are growing and in which the already low levels of oxygen are declining [2, 3]. Deoxygenation is also evident in coastal marine ecosystems, where stratification, alone or in combination with localized nutrient inputs and heating, can generate hypoxic water masses [1, 4]. Dissolved oxygen concentration (hereafter "DO"), therefore, does not vary independently of other physical conditions and has been demonstrated to correlate positively with pH and negatively with temperature in many nearshore environments [5] [6] [7] [8]. Acidification is expected to have the most detrimental effects on marine life at high latitudes, while deoxygenation is expected to have more severe impacts in the tropics [5]. Heating may also have a greater impact in the tropics where many organisms live near their critical thermal maximum [9] [10] [11] [12]. Here, we examine how spatial and temporal variation in DO and other co-varying environmental conditions influence the abundance and composition of zooplankton in a tropical estuary. 

Zooplankton act as important links in pelagic food webs and benthic-pelagic coupling. Zooplankton play a vital role in facilitating nutrient cycling [13, 14], transferring carbon to the seafloor via their feces and bodily remains [15, 16], and limiting phytoplankton blooms via grazing [13]. Zooplankton serve as trophic links between phytoplankton and nektonic communities, serving as important food sources for the larger animals in the water column, which in turn provide economic benefits through fisheries. While much of the zooplanktonic community is made up of organisms that spend their entire lives in the water column (holoplankton), there can also be substantial numbers of meroplankton: the planktonic larval stages of primarily benthic animals, which serve to promote the dispersal of animals, especially those with sedentary adult stages [17]. Poor conditions for larval growth and survival can significantly impact larval settlement and recruitment, potentially driving variation in adult demographics and shaping supply-side limitations observed by benthic ecologists [18] [19] [20] [21]. It is therefore evident that environmental conditions that impact zooplankton abundance and diversity are likely to have a number of knock-on effects in both benthic and pelagic ecosystems. Yet, the impacts of hypoxia on zooplankton have been studied in detail for few systems, and information is particularly lacking in the tropics. 

It is fairly well-documented that hypoxia at depth results in some amount of vertical habitat compression of the zooplankton community, at least in temperate estuaries [19, 22], fjords [23], and the Gulf of Mexico [24]. In the Gulf of Mexico overall plankton biomass moved up in the water column by as much as 7 m when deep waters were hypoxic, and large zooplankton were absent from hypoxic water [24]. A similar result was reported in Chesapeake Bay, where copepods, normally abundant in deeper waters, were absent from the same depths when it was hypoxic [25]. By contrast, moderate hypoxia in a coastal fjord in Washington State was not found to impact the distribution of fish or zooplankton [26], while a more recent study in the same area found lower total abundance and biomass at hypoxic sites [23]. Meroplankton can also concentrate above the oxycline, possibly preventing settlement on the deeper benthos [27]. 

Hypoxia can also alter community composition. Some groups are more sensitive to hypoxia than others, resulting in changes in community structure and trophic transfer. For example, fish in the Baltic Sea are restricted to the well-oxygenated shallower water, while crustacean zooplankton distributions peak in lower oxygen or hypoxic water, suggesting that crustaceans find a refuge from predation where oxygen levels are too low for fishes to perform well [28]. Hypoxic water is also often associated with increases in or dominance of gelatinous zooplankton [22, 23]. Significantly reduced abundance in hypoxic water has been documented for several groups, including larvaceans [29, 30], calanoid copepods [22, 25, 29], and bivalves [29, 31]. 

Despite the importance of holoplankton and meroplankton on marine ecosystem function, there is still a lot to learn about how hypoxia shapes plankton abundance and distribution in the coastal environment. One of the major challenges to this work is the difficulty of disentangling the effects of hypoxia from other co-varying environmental factors (e. g., temperature, pH). Low DO is often associated with lower pH and warmer temperatures, obscuring which of these potential stressors may be driving any observed patterns. It is also notoriously difficult to find robust patterns of association with plankton diversity and abundance, as zooplankton are often patchy across space and time [32, 33]. 

To try to disentangle these factors and understand which planktonic groups are most sensitive to extremes in physical conditions, as well as which of those conditions have the greatest influence on the zooplankton community, we sampled zooplankton and environmental conditions weekly for one year in a tropical estuary that experiences seasonal hypoxia. Our goal was to answer the following questions: (1) How does abundance of each major group of zooplankton vary across sites and between mid-and deep-water? (2) Which environmental factors, including hypoxia, are the most important in explaining variation in the plankton community as a whole across the bay? (3) Which environmental factors best explain the variance in abundance of each group of zooplankton? Because DO co-varies with many other environmental factors, and because count data is often non-normal and zero-inflated, our data were not well suited to common statistical analysis. We therefore take multiple analytical approaches and discuss the commonalities across their results.
section_name: Introduction

section_name: Methods

section_content: Almirante Bay (Figure 1 ) is an enclosed lagoon with estuarine characteristics on the western Caribbean coast of Panama [34] [35] [36]. Seasonal rainfall, combined with low winds and small tidal amplitudes, result in seasonal stratification of the water column, thermal inversions, and hypoxia [36] [37] [38]. Further onshore, oxygen is reduced at depth throughout the year, but hypoxia generally expands and shoals in early summer and then dissipates in December. The spatial and temporal extent of shoaling of the hypoxic bottom waters varies between years and, in extreme cases, reaches local coral reefs and results in significant die-offs of benthic fauna [39, 40].
section_name: Study Site

section_content: We sampled zooplankton weekly at 3 sampling sites (Figure 1 ), which were selected to capture the range of hypoxic conditions in the bay: STRI is nearest to the mouth of the bay and experiences light, episodic hypoxia, with DO seldom dropping below 2 mg L -1 even in deeper water; Cristobal is in the middle of the bay and experiences moderate hypoxia, with low DO in deeper water during part of the year; and Pastores is furthest onshore, experiencing extreme, persistent hypoxia with low DO throughout the water column seasonally and in deeper water almost year-round [41]. Plankton samples were collected at ~11 m (mid-water column) and ~21 m (near the bottom) every week for one year. 

Immediately prior to plankton collection at each site, hydrographic conditions were measured using a multiparameter water quality sonde (YSI EXO2), set to log data every 5 s. These measurements were: temperature, salinity, dissolved oxygen (DO), pH, chlorophyll, and dissolved organic matter (DOM). The sonde was calibrated biweekly following the manufacturers' instructions; however, the pH sensor suffered from substantial sensor drift, and therefore pH data were excluded from subsequent analyses. Nevertheless, pH is chemically linked to water temperature and oxygen content [1], and these parameters have been reported to co-vary in Almirante Bay [7]. 

To obtain measurements, the sonde was slowly lowered by hand through the water column in ~1 m increments and held at each depth for 5 s. On 8 dates, the sonde was unavailable on the sampling day, and we used sonde data taken at the same site at a similar time of day collected earlier that week. On 3 other dates, the sonde was not available during the week, so we used a handheld oxygen, salinity, and temperature meter (YSI Pro 2030). In these instances, water was collected from ~11 and ~21 m in a Van Dorn bottle immediately before plankton samples were collected, and measurements were made in the boat. On 1 other date, sonde data was corrupted, and values were interpolated from the weeks immediately before and after. 

To visualize environmental conditions throughout the water column and across the year at each site, sonde data were averaged by 1 m depth bins (because not every cast had only one value per depth bin, and plotting software required identical axis values) and plotted using the R function levelplot from the latticeExtra library [42]. Any missing depths were linearly interpolated using same-day data from the nearest depth bins above and below. For the 4 dates without sonde measurements, data were interpolated using values at the same depth from the weeks immediately before and after. Plankton were sampled weekly from 10 January to 19 December 2017 in midwater (11. 17 ± 0. 06 m) and near the benthos (21. 4 ± 0. 14 m), always between 8:00 and 15:00. These 2 depths were chosen to differentiate between (1) reduced overall abundance due to seasonal variation or, for meroplankton, variation in adult reproduction and larval supply, and (2) reduced abundance from only the hypoxic portion of the water column due to avoidance or mortality. Surface waters were not sampled due to the periodic occurrence of a freshwater lens and salinity stratification resulting from the ~3 m of rain per year experienced in this region [37]. Plankton were collected under collection permits SE/A-79-16 and SE/APHBO-9-18, issued by the Panamanian Ministry of the Environment. 

To obtain plankton samples, we used a gasoline-powered water pump to collect 379 L of seawater via an 80 mm intake hose, with an additional sensor (either a Star-Oddie DST CTD or a DiveNav BlueBuddy) affixed just above the hose intake on 36 of the 50 weeks to measure exact sampling depth. Pump rates were approximately 20 to 50 L min -1. Water was filtered through a 125 µm mesh plankton net, and samples were immediately preserved in ethanol. In the laboratory, a total of 298 plankton samples were sorted using a stereomicroscope, and the major groups of zooplankton and meroplankton were enumerated. We focused on groups with robust skeletons, as sampling with these kinds of pumps generally damages or destroys soft-bodied groups (e. g., polychaete larvae) making it impossible to identify the material. In addition, strong swimmers (e. g., crab zoeae and fish larvae) may be able to avoid this type of pump intake. We found very few of these larvae, and they were therefore excluded from our analysis. With these limitations in mind, holoplankton groups that were abundant enough to analyze and in sufficiently good condition to identify to phylum or major taxonomic group were: copepods, larvaceans, chaetognaths, and uncoiled pteropods (almost exclusively casseids). Meroplankton were separated into the following groups: bivalve veligers, gastropod veligers, plutei (echinoid and ophiuroid), barnacle nauplii, and barnacle cyprids. As it is not possible to identify individuals from most of these groups to a lower taxonomic level, and as we wanted to obtain the broadest possible view of the community, we performed our analysis on the abundance of these major groups (i. e., on the total number of individuals per site, depth, and zooplankton group). There is no doubt that the patterns we recover are the result of compounded distributions and abundance patterns of different species and developmental stages. Therefore, our focus is on broadly describing the community as a whole, to the extent that it can be sampled with our methods.
section_name: Sampling

section_content: All sonde data from ±1 m of each plankton-sampled depth were averaged to yield conditions experienced by plankton at the time of sampling. In these instances, for the 4 dates without sonde data, handheld YSI data were used when available for oxygen, temperature, and salinity, and values for DOM and chlorophyll were interpolated as described above. 

To explore the covariation of environmental conditions at 11 and 21 m, we conducted a principal component analysis (PCA) on DO, temperature, salinity, chlorophyll, and DOM data, using the rda function in the vegan package in R [43], conducted on scaled data to account for different measurement units. Community data was similarly visualized using a PCA, conducted with the same rda function, using log(x + 1) transformed counts to meet the assumptions of linear relationships between variables. Results were plotted separately for each site. To explore links between the ordinations of the biological community and the environmental variables, we used the function envfit from the vegan package to project vectors of the environmental parameters [DO, temperature, salinity, chlorophyll, and DOM, as well as depth (mid or deep), and site] onto the community ordination plot. 

To explore which plankton groups were more commonly found in mid-or deepwater, which could affect the range of environmental conditions they were exposed to, we calculated the difference in group-specific abundance between mid-and deep-water samples for each site and week. To further assess whether observed decreases in plankton abundance in deep samples were the result of (a) behavioral avoidance of (or mortality due to) stressful conditions, or (b) decreased abundance throughout the water column because of, for example, reproductive cycles, we calculated the Kendall Tau correlation coefficient between log(x + 1) transformed abundances at 11 m and 21 m for each group. 

To better examine the drivers of community composition and abundance over space and time, we applied a generalized linear model (GLM) using the R package mvabund [44, 45], which allows for model analysis of community data. For predictor variables, we included DO, temperature, salinity, chlorophyll, and DOM, all two-way interactions between them, as well as depth (mid or deep). Model selection was carried out via iterative use of the mvabund manyglm function [44], with stepwise removal of non-significant terms until all model parameters were significant at a threshold of p = 0. 05. The same model selection process was carried out for each plankton group separately, using the glm. nb function from the MASS package in R [46], in order to explore whether environmental predictors varied between plankton groups. In all GLM models, we assumed a negative binomial distribution, to account for the quadratic mean-variance relationship common to count data. 

Because GLMs are known to be less robust against cross-correlated independent variables, we also applied Partial Least Squared analysis (PLS) in JMP 16. We chose a PLS because the environmental variables have strong covariance structure, and this approach is not as highly influenced by covariation between the explanatory variables as are approaches, such as multiple regression [47]. Prior to analysis, count data were log(x + 1) transformed, and residual plots were examined to ensure approximate normality. We report the coefficients of the centered and scaled data for variables with a Variable Importance in Projection (VIP) > 0. 8 (the default in JMP) and more conservatively with VIP > 1. 0 [48]. 

Finally, based on results from the above analyses, we used linear models to describe how log(x + 1) transformed abundances varied with log-transformed DO. To better identify threshold values in DO, we used the segmented package in R [49] to model log(x + 1) transformed abundance versus DO with 1 estimated breakpoint, and we report subsequent breakpoints for all models with R 2 > 0. 25. Additionally, based on results from breakpoint analysis of pluteus larvae, we re-ran the linear model for this group to instead describe log(x + 1) transformed abundance as a function of (DO) 2. All analyses were conducted in R v 4. 1. 1 [50] unless otherwise stated.
section_name: Analysis

section_name: Results

section_content: As expected from previous studies, the environmental data demonstrated a seasonal pattern in temperature, salinity and DO. Early in the year (January-May), the water was relatively cool, salinity and DO were high, and there was little stratification in the water column (Figure 2 ). This period coincided with the dry season in Panama, when the trade winds blow across the Isthmus. During the rest of the year, the country experiences monsoonal conditions, which were evident in reduced surface salinity, increased stratification, and decreased DO at depth (Figure 2 ). During September and October, there is a period locally known as the small summer, which is characterized by low rainfall, high water temperatures, and calm seas. This was the period of peak hypoxia. Coinciding with these warm, hypoxic waters were increased levels of chlorophyll and DOM (Figure 2 ). The onshore-to-offshore gradient in hypoxia is also evident in Figure 2, with Pastores experiencing the most severe hypoxia at depth and STRI experiencing the least. 

Figure 2. Environmental conditions at each site over the period of the project, as measured weekly by the YSI sonde: Dissolved oxygen (DO, mg L -1 ); temperature ( • C), salinity, chlorophyll (µg L -1 ), and dissolved organic matter (DOM, QSU). Value ranges for visualization were chosen to maximize visible contrast of conditions present at 11 and 21 m (highlighted with cyan rectangles), resulting in the truncation of some outlier values in salinity, chlorophyll, and DOM, mostly at the surface and benthos. 

Principal component analysis (PCA) of environmental conditions across 11 and 21 m demonstrated that the largest amount of variation, as explained by Principal Component 1 (PC1; 36. 3%), arranged points with respect to DO versus other variables, with higher DO occurring at lower values of PC1 and higher values of temperature, salinity, DOM, and chlorophyll occurring at higher values of PC1 (Figure 3 ). PC2 arranged the data with respect to temperature and salinity versus chlorophyll and DOM. Overall, the points were not well separated, but rather clustered around 0, and only ~60% of the variation in conditions was explained by the first two principal components.
section_name: Environmental Conditions

section_content: A summary of the abundances of each taxon at each site and depth are given in Figures 4 5 6 and Table S1. Copepods were the most abundant in the bay during the dry season when the water column was well mixed, and the overall abundance of the community was reduced during the periods of deoxygenation later in the year (Figures 2 and 4-6 ). Mapping community composition at the six site-depth combinations on the same principle component axes, generated by a global analysis of the plankton community, provides a more nuanced view of the differences between sites and depths (Figure 7 ). Overall, the communities at the three sites fell into the same PC space, showing that despite these subtle differences described below, the overall community compositions were fairly similar (Figure 7B-D ). Principle component 1 (PC1; 36. 4% of total variation) reflects abundance, with higher abundance of all taxa increasing to the left (Figure 7A ). PC2 (17. 9% of total variance) reflects a split between pteropods, larvaceans, and plutei in the negative direction and bivalves, gastropods, chaetognaths and barnacle larvae in the positive direction (Figure 7A ). When environmental data were mapped onto the community ordination plot, all variables (DO, temperature, salinity, chlorophyll, DOM, depth, and site) were indicated as significant predictors (p < 0. 01; Figure 7A ), though the relative explanatory power was greatest for DO (R 2 = 0. 52) and DOM (R 2 = 0. 37). 

Looking at the data by site, in Pastores and Cristobal, low DO is associated with reduced abundance (i. e., samples from weeks and depths with lower DO fall further to the right on PC1; Figure 7C, D ). At STRI point and Cristobal, the deep samples tend to have a higher score on PC2, indicating increased relative abundance of most meroplankton and chaetognaths compared to shallow and more oxygenated samples, which have higher relative abundance of pteropods, larvaceans, and plutei (Figure 7A-C ). The relationship between abundances at 11 m and 21 m (Figure 8 ) can indicate if changes in abundance are related to shifts in depth distributions or in more global shifts in abundance, such as those resulting from reproductive phenologies. A positive slope for this relationship would indicate that abundances at the two depths were positively correlated, a flat slope would indicate no relationship, and a negative slope would indicate a shift in depth distribution. Kendall Tau correlation coefficients for each taxon per site are reported in Table 1. All of the 14 statistically significant pairwise comparisons had positive correlation coefficients, and only two of the 13 non-significant relationships had a negative correlation, providing strong evidence that abundances at the two depths were either un-related or tended to experience parallel changes (Table 1 ). Distributions of the differences in abundance between 11 and 21 m for each taxon per site are shown in Figure 8, which demonstrate that pteropods, larvaceans, gastropods, and plutei were usually more abundant in midwater, while copepods were generally more abundant in deeper water. Notably, the depth at which bivalves and barnacle cyprids were most abundant differed between sites (Figure 8 ).
section_name: Patterns of Abundance across Sites and Depths

section_content: To understand how environmental factors shaped the entire plankton community, we applied both PLS and GLM models to the dataset. PLS analysis counld only include continuous variables (and therefore could not block for site or depth). The model explained only 19% of the overall variance in log(x + 1) transformed plankton abundances. Temperature, salinity, DO, and DOM were considered important (VIP > 0. 8), but only DO and DOM crossed the VIP > 1. 0 threshold, indicating that they were the most important factors shaping the variation in the total plankton community across the seascape of Almirante Bay. The GLM analysis, which blocked for depth, generated a model that included significant effects of all five environmental variables either as main effects or through two-way interactions with other variables (Table 2 ).
section_name: Total Community Composition and Environmental Factors

section_content: Univariate PLS and GLM analyses focused on the importance of the different environmental variables on the abundance of each specific plankton group. The PLS analysis run across all sites and depths generated significant fits for all plankton groups other than barnacle nauplii (Figure 9 ). For the eight other groups, DO yielded a VIP > 1. 0 for all groups and was positively related to abundance (Figure 9 ). DOM also had a VIP > 1. 0 for all groups except for pteropods and chaetognaths, which each had a VIP > 0. 8 (Figure 9 ). For all eight groups with significant fits, DOM was negatively related to abundance (Figure 9 ). Salinity was strongly negatively related to abundance of pteropods, and temperature was strongly positively associated with abundance of pteropods and negatively related to abundance of copepods and gastropods (Figure 9 ). Because abundances and community composition, and therefore the realized range of environmental conditions, varied with depth (see above), we repeated the PLS analysis for all of the mid and deep samples separately. When analyses were separated in this way, the environmental conditions did not significantly explain patterns of abundance of copepods, chaetognaths, bivalves, or cyprids at 11 m, but they explained between 5% and 37% of the variation in the other taxa at 11 m (Figure 9 ). At 21 m, 13-44% of the variation was explained for every group, and the details of the relationships were broadly similar to those in the total datasets: DO had a positive coefficient and a VIP > 1. 0 in every analysis with a significant fit, except for nauplii at 11 m, in which the PLS only explained 5% of variation. DOM had a VIP > 1. 0 in all taxa at 21 m, but the direction varied with a negative coefficient for gastropods, plutei, bivalves, nauplii, and copepods and a positive coefficient for the other three groups (Figure 9 ). 

Univariate GLM models, which included two-way interaction terms as well as depth category, also predicted only ~16-41% of variation in plankton abundance (Table 3 ). Notably, DO was the only predictor to remain in all abundance models: DOM as both a main effect and interaction was dropped from the model of copepod abundance but remained in the models for the other taxa; temperature and salinity remained as important predictors in all models except for barnacle nauplii, chlorophyll remained in all models except for bivalve and plutei, and depth remained as a significant predictor for only pteropods, larvaceans, copepods, plutei, and barnacle nauplii (Table 3 ). Out of the interaction terms, the interaction between DO and DOM remained in the most models (seven of nine; pteropods, chaetognaths, larvaceans, bivalves, gastropods, and both barnacle nauplii and cyprids), followed by the interactions between DO and temperature (six of nine; pteropods, chaetognaths, copepods, gastropods, plutei, and barnacle cyprids), temperature and salinity (five of nine; larvaceans, copepods, bivalves, gastropods, and plutei), and temperature and DOM (five of nine; larvaceans, bivalves, gastropods, plutei, and cyprids; Table 3 ). In other words, DO, DOM, and temperature-as main effects or interactions-were most consistent in their explanatory power for abundance of individual zooplankton groups. On the other hand, all models retained at least four of the six main effects, indicating that abundance was best explained by the combination of environmental factors, rather than by one or two alone. Examination of the relationships between abundance and DO (Figure 10 ) demonstrated that, for several taxa, abundance decreased gently with decreases in DO between 3-6 mg L -1, and that an abrupt decrease in abundance occurred at or below 2 mg L -1. This relationship was most evident at the deeper sampling sites, as the mid-water sites seldom experienced very low DO (Figures 2 and 10 ). Piecewise linear regression analyses (Figure S1 ) of abundances in deep samples identified breakpoints between 1. 6 and 4. 2 mg L -1, with copepods exhibiting the lowest breakpoint and plutei exhibiting the highest (Figure 10 and S1 ). All other group breakpoints clustered between 2 and 3 mg L -1 (Figure 10 and S1 ). 

Because oxygen availability typically varies with temperature, we also explored relationships between abundance and percent saturation of oxygen, but regression model fits were equal to or worse than the respective model using mg L -1 DO (Figure S2 ). Similar relationships were explored between abundance and DOM, as this was identified as a significant predictor in PCA and PLS analyses (Figures 7 and 9 ), but resulting regressions were noisy, and linear model fits were poor (Figure S3 ), indicating that much of the impact of DOM likely plays out through interactions with other variables. This is supported by the fact that DOM as a main effect was not significant in four of the nine univariate GLMs (Table 3 ). Temperature alone explained <15% of variation in abundance for seven of the nine groups at 21 m, and for all nine groups at 11 m (Figure S4 ). PC1 of the environmental conditions (Figure 3 ) explained a greater amount of variation in abundance, but models using DO alone performed better for all groups at 21 m and for six of the nine groups at 11 m (Figures S5 and 10 ). ) and order-of-magnitude concentration of the nine groups of zooplankton. Trendlines lines show the best fit for the abundances at each depth (red for 11 m and blue for 21 m), assuming concentration ~log(DO) for all groups except plutei at 21 m, which assumes concentration ~(DO) 2. Vertical dotted lines denote 'breakpoints' (±SE; grey rectangles) for taxa at 21 m, as identified in piecewise linear regression (Figure S1 ) and shown only for models with R 2 > 0. 25. Point shape denotes site.
section_name: Environmental Factors Shaping the Abundance of Each Group

section_content: Environmental conditions in coastal environments are highly variable, even in the tropics. Although the wide range of observed values for each environmental factor should provide sufficient variation to support rigorous statistical assessment of its role in shaping organismal diversity and abundance, the non-independence of the factors makes this a challenge. Hypoxia is physically linked to temperature, as warmer water holds less oxygen, and biologically linked to acidification, as community respiration consumes DO and releases CO 2 [51]. Hypoxia is also causally linked to eutrophication, which can drive DO down by fertilizing microbial respiration [51] [52] [53] [54] [55]. Low salinity water, which can lead to stratification, can also promote hypoxia [4]. These relationships make it virtually impossible to design a sampling scheme that can easily disentangle the effects of these different factors. 

Thus, as expected, the environmental variables we measured demonstrated a complex pattern of correlations. Temperature, pH, and DO have been shown previously to co-vary in Almirante Bay [7], with pH and DO exhibiting a positive correlation and each exhibiting a negative correlation with temperature. It is unfortunate that we were unable to use the pH data for this analysis, but it is thought that pH should be less important in the tropics than in temperature systems [5], and experimental work with benthic taxa supports this conclusion over the range of extreme values currently observed in the bay [7]. Our present analysis also demonstrates that temperature and salinity co-varied, as did chlorophyll and DOM (Figure 3 ). Our approach, especially the PLS, was designed to disentangle these correlations as much as possible, to distinguish which of these covarying factors were most important in shaping plankton abundances. Amongst the five environmental factors examined here, DO and DOM appeared to be the main drivers of differences in plankton community composition, as they had the largest explanatory power in the PCA, had the highest VIP scores in the PLS, and were the most highly significant in the GLM. In addition to the results of the community analysis, these two factors were the most consistently important in describing the abundance of each taxon separately. Below, we discuss each of these in turn.
section_name: Discussion

section_content: The negative impacts of coastal hypoxia and ocean deoxygenation are well documented across an array of taxa and habitats [19, 24, 40, 56], and Almirante Bay is well-known to suffer from seasonal hypoxia. We therefore expected to detect an impact of hypoxia on the plankton community. What was surprising to discover was that a relatively large portion of the variation in abundance of some groups at 21 m could be explained by DO alone (>30% in larvaceans, copepods, bivalves, and gastropods). It was also surprising that large reductions in abundance of some groups, particularly larvaceans and gastropods, can occur in water as shallow as 11 m (Figure 10 ). Groups that were less abundant and therefore less well sampled demonstrated a poorer fit between abundance and DO, but nevertheless, a marked and rapid reduction in abundance of most groups is evident at conditions below 2 mg L -1. 

Long-term monitoring data show that DO concentrations fall below this benchmark 16% of the time in Almirante Bay at 21 m depth (Adelson et al. submitted). This 2 mg L -1 definition of hypoxia is based on organismal responses in numerous experimental studies that demonstrate significant impacts on most organisms at this oxygen level [57]. However, there is also ample evidence that sub-lethal effects can occur at higher concentrations for many macroscopic marine animals [58, 59]. Our data suggest this is also the case for planktonic animals, as the abundance of several taxa decreases with decreasing DO over the range of 2-4 mg L -1 (Figure 10 and Figure S1 ), which is equal to ~30-60% oxygen saturation at the temperatures and salinities observed in the bay. 

It is difficult to compare previous studies of the relationship between zooplankton and oxygen concentrations in coastal systems with our results. Previous studies have generally tested for multivariate associations between entire communities and various environmental variables, rather than explicitly examining at what level of oxygen availability decreases are evident. For example, Keister et al. (2020) found that DO did occasionally drop below 2 mg L -1 and that hypoxic sites were associated with lower full water column abundance of zooplankton [23]. Likewise, Roman et al. (2012) binned samples into normoxic (>2 mg O 2 L -1 ) and hypoxic (<2 mg O 2 L -1 ) categories and found that overall median zooplankton volume was 2. 5 times higher in the normoxic than hypoxic samples [24]. Binning significantly reduces our ability to understand what DO levels impact organismal abundance, which is vital for understanding how DO impacts ecosystem function and species distributions. In order to obtain this information, a univariate correlative approach, such as the one taken here, can be applied. Alternately, an experimental approach can be taken, which documents the physiological tolerances of the organisms and relates that to abundances in the plankton. This has approach has been usefully applied to understanding occurrences of one coastal copepod [60], but overall few invertebrates and fishes have been studied in this way [61, 62].
section_name: Impacts of Hypoxia on Plankton and the Almirante Bay Ecosystem

section_content: DO and DOM have a complex relationship in the coastal zone. It is broadly understood that excessive nutrient inputs can lead to phytoplankton blooms, resulting in increased particulate organic matter, the decomposition of which draws down oxygen [58]. Such inputs can come from natural terrestrial sources, anthropogenic sources, resuspension of nutrients in marine sediments, and even directly from the atmosphere [58, 63]. The resulting low DO can alter the redox potential of the sediment, which increases the dissolution of NH 4 and PO 3 into the water, creating a positive feedback loop [64]. Another positive feedback loop between DO and DOM is mediated by benthic animals, whose burrows promote nitrification and denitrification, thereby drawing down nitrogen concentration of nearby water [56, 64]. As the bottom becomes more hypoxic, the number of these infaunal animals, and this mechanism for the removal of nitrogen from the water, decreases [64, 65]. Our results suggest that the reduction in meroplankton at low DO could exacerbate this situation through reduced recruitment of important infaunal organisms, such as bivalves, possibly leading to supply-side limitation in population size or diversity as well as limiting recolonization after local extinction events. As larvae are more sensitive to stressors than adults [66] [67] [68], it is possible that DO levels that are not directly harmful to adults could reduce benthic infaunal diversity through limiting recruitment. 

DOM can also interact in complex ways with zooplankton. For example, high levels of terrigenous DOM can reduce light levels as well as fertilize primary producers, which can have a bottom-up effect on zooplankton grazers [69] and references therein. Mesocosm experiments demonstrated that high terrigenous-dissolved organic matter (tDOM) treatments resulted in an increase in crustacean zooplankton biomass and C:N ratio [69]. The authors suggest that this increase likely resulted from an increase in phytoplankton biomass despite the light reduction caused by the tDOM [69]. It seems likely that natural nutrient concentrations, as well as the addition of DOM, could result either in overall positive or negative impacts. Shifts in the relative abundances of different microbes is also likely to impact specific zooplankters differently, as these taxa employ different feeding mechanisms, selecting different size classes of particles or even utilizing DOM directly [70] [71] [72] [73]. Our results suggest a general negative relationship between zooplankton abundance and DOM for most groups sampled in Almirante Bay, but this relationship is likely mediated heavily through interactions with other environmental variables (Table 3 ).
section_name: DOM and Zooplankton in Hypoxic Environments

section_content: One clear result of our univariate tests is that each zooplankton group shows a different response to the suite of environmental factors. This is evident in both the importance and direction of coefficients for each factor in the PLS, and in the pattern of significance in the GLMs. Most notably, the PLS demonstrated a diversity of responses to temperature. Gastropod abundance demonstrated a clear negative response to increasing temperature (both overall and at 11 m), while pteropods, bivalves, and both stages of barnacle larvae demonstrated an important positive response to temperature at 21 m, and the remaining groups gave little evidence of any important role of temperature in shaping abundance. It should also be noted that the impact of temperature and DO in shaping bivalve larval abundance has already been analyzed in detail [31]. As for other environmental factors, cyprids and chaetognaths demonstrated a positive response to chlorophyll overall, and at 21 m, nauplii and plutei demonstrated a positive response at 11 m, gastropods demonstrated a negative response at 11 m, and no other group appeared to be impacted by this factor. Lastly, pteropods demonstrated a surprising, strong negative response to increased salinity. This diversity in responses explains the low explanatory power for the fully multivariate models, which demonstrate that most variables are important or significant while only low percentages of variation in the overall community are explained. Similarly, the univariate GLMs indicated the important role played by two-way interactions between the different environmental variables (e. g., temperature, DOM). It is highly likely that three-way and higher-level interactions are similarly affecting plankton communities in this and similar systems, though they are difficult to conceptualize. 

The diverse responses among taxa, and the fact that abundances of most taxa are shaped by more than one factor, has important implications for future community responses to environmental change. Modeling approaches aimed to predict changes in organismal distributions generally use a single explanatory factor, usually temperature, e. g., [74, 75]. Recently, the use of the metabolic index, a metric reflecting the oxygen and temperature dependence of ectothermic metabolism [76], has been used to effectively model abundances of anchovies in the California current system [62], as well as to explain the geographic distributions of a variety of marine ectotherms and the habitat occupancy of tropical brittle stars [8, 76]. Our results suggest that in the coastal zone, inclusion of DOM might significantly improve predictive modeling of the distributions of plankton.
section_name: Taxon Specific Responses

section_content: Current environmental conditions in Almirante Bay, particularly dissolved oxygen, are already sufficiently extreme that impacts can be observed on all of the holo-and mero-plankton groups studied here. However, much of the impact plays out through interactions between environmental conditions, likely due to complex biogeochemical feedback loops and the inherently patchy dynamics of plankton, such that no statistical model explained >52% of variation in community or group-specific abundance. Dissolved oxygen concentration had the most explanatory power for these taxa. All models agreed that dissolved organic matter, alone or through interactions with other parameters, played an unexpectedly large role in shaping zooplankton communities alongside dissolved oxygen, suggesting that more attention should be paid to these relationships in future studies Supplementary Materials: The following supporting information can be downloaded at https: //www. mdpi. com/article/10. 3390/jmse10030427/s1: Spanish translation of abstract; Table S1 : Summary statistics for plankton count data;
section_name: Conclusions

section_name: Acknowledgments:

section_content: We thank the staff, students, and volunteers who provided support at the Bocas del Toro Research Station, and especially Isis Ochoa and Tania Romero who helped count the copepods. We also thank Robert Thacker for his suggestions of statistical analyses, Phil Yund and four anonymous reviewers for their comments on the initial manuscript, and Maycol Madrid for the Spanish translation of the abstract. This is contribution 100 from the Smithsonian's MarineGEO and Tennenbaum Marine Observatories Network.
section_name: Acknowledgments:

section_name: Institutional Review Board Statement: Not applicable.

section_content: Funding: This research was funded by a Fulbright Fellowship awarded to J. B. W. and a Smithsonian Tropical Research Institute internship awarded to L. V.

section_content: Informed Consent Statement: Not applicable.
section_name: Institutional Review Board Statement: Not applicable.

section_content: Data Availability Statement: Data are available at www. mdpi. com/article/10. 3390/jmse10030427/ s1. Code for data analysis and visualization are available at https://github. com/jbweinstock/Bocaslarvae, accessed on 1 March 2022.

section_name: Conflicts of Interest:

section_content: The authors declare no conflict of interest.
section_name: Conflicts of Interest: