Big or small, patchy all: Resolution of marine plankton patch structure at micro- to submesoscales for 36 taxa.

Despite the ecological importance of microscale (0.01–1 meter) and fine-scale (1 to hundreds of meters) plankton patchiness, the dimensions and taxonomic identity of patches in the ocean are nearly unknown. We used underwater imaging to identify the position, horizontal length scale, and density of taxa-specific patches of 32 million organisms representing 36 taxa (200 micrometers to 20 centimeters) in the continental and oceanic environments of a subtropical, western boundary current. Patches were the most frequent in shallow, continental waters. For multiple taxa, patch count varied parabolically with background density. Taxa-specific patch length and organism size exhibited negative size scaling relationships. Organism size explained 21 to 30% of the variance in patch length. The dominant length scale was phylogenetically random and <100 meters for 64% of taxa. The predominance of micro- and fine-scale patches among a diverse suite of plankton suggests social and coactive processes may contribute to patch formation.

INTRODUCTION

Spatial heterogeneity in plankton abundance occurs across a range of scales in aquatic ecosystems (–). While the etymological root of plankton (i.e., “wanderer” or “drifter”) implies a random distribution, multiple studies have inferred from trophic arguments (–) or direct observations (–) that plankton are heterogeneous in space. This spatial heterogeneity, where planktonic organisms exhibit higher standardized densities (e.g., individuals or biomass per cubic meter) relative to the background density, is described as patches, aggregations, or clusters in the literature. In the aquatic environment, patches can form and persist at a number of spatial scales, from the microscale (0.01 to 1 m) to the mesoscale (10 to hundreds of kilometers) in both the vertical () and horizontal dimensions (). Moreover, they can exist at multiple scales simultaneously, with smaller patches comprising larger aggregations (, ).

The wide range of spatial scales over which plankton patchiness can occur indicates that a variety of mechanisms underlie patch development and maintenance, with consequences for individuals, populations, and ecosystem dynamics (). Physical processes (e.g., fronts, internal waves, thermal stratification, eddies) often drive the formation of plankton patchiness (–), but biological processes are also important drivers—particularly at the smaller scales (e.g., <100 m) where behaviors such as diel vertical migration can play a key role (). Theoretical and laboratory studies have demonstrated that fine-scale (1 to hundreds of meters) patchiness is critical for individual fitness, including locating mates, finding adequate prey to support metabolic needs and growth, and avoiding predators (, ). At fine vertical scales (i.e., “thin layers”), in situ plankton aggregations are often composed of distinct taxonomic assemblages () whose members interact with one another (). This biological spatial heterogeneity has implications for modeling marine ecosystem dynamics as plankton assemblage structure influences community organization and stability (), nutrient cycling (), grazing rates (, ), larval fish survivorship (–), and trophic energy transfer efficiencies (, ). Moreover, early evidence suggests that including details of predator and prey distributions at the fine scale in mechanistic ecosystem models could improve fishery predictions, particularly in coastal environments ().

Despite the importance of plankton patchiness over a range of spatial scales, we have a relatively poor understanding of the mechanisms underlying the generation and maintenance of patches. Difficulty in simultaneously sampling a broad range of organisms has led to a poor understanding of whether and how patch characteristics (size and density) vary with taxonomic group or organism size. This difficulty arises despite the extensive work describing how plankton physiology, behavior, and ecology vary as a function of organismal size [(, ) and references therein]. Little is known about the relationship between individual plankter size and the spatial scale of their aggregations (i.e., taxa-specific patch size), and the mechanisms that may drive such relationships. It is sensible to imagine that such a relationship exists: If patch formation helps small organisms “appear” larger to predators (e.g., fish schooling) and there is a finite number of larger predators, then advantages of patch formation should be greater for smaller organisms. Alternatively, if patch formation was primarily for the purposes of identifying favorable feeding environments (), then it may be easier for larger organisms with greater sensory and behavioral abilities to aggregate. A recent modeling effort indicates that the collective behavior on the part of a diversity of organisms may increase ecosystem stability (). One difficulty of examining the size scaling relationship of plankton patch formation is that it requires the simultaneous measurement of organisms spanning multiple orders of magnitude and their patches. Resolving this relationship requires sampling not only where there is a range of sizes and taxonomic groups present but also where a single group does not overwhelm the abundance of others, namely, a dynamic, subtropical marine ecosystem.

The importance of fine-scale structure in planktonic systems has motivated the development of high-resolution acoustic and optical technologies that can simultaneously sample across a broad range of spatial scales (, ). Traditional plankton nets integrate organisms over the sampling distance and depth, producing a spatial mean rather than measures of fine-scale variability. Moreover, nets do not sample ecologically important, but delicate, gelatinous plankton (, ). Acoustic systems, while able to yield high-resolution data for hard-bodied zooplankton such as euphausiids, cannot resolve gelatinous plankton that have a low acoustic impedance contrast (i.e., low backscatter). Optical systems like the Video Plankton Recorder, Shadow Image Particle Profiling Evaluation Recorder, Lightframe On-sight Keyspecies Investigation, Zooplankton Visualization and Imaging System, and Underwater Vision Profiler have produced outstanding results with regard to sampling gelatinous and nongelatinous plankton (, ). However, these systems do not adequately quantify rare taxa such as larval fishes (which typically occur at densities of ca. 0.01 to 0.001 individuals per liter) or larger macroplankton (2 to 20 mm) because of either limited sampling volumes (i.e., ≈2 liters/s) or deployment mode (). Other optical sampling technologies such as laser optical counters have yielded information about submesoscale plankton patch structure but cannot resolve specific taxa (). For these reasons, a concurrent understanding of how fine-scale patchiness varies within and among taxonomic specific plankton groups that span the size mesozooplankton spectrum (0.2 to 20 mm) is lacking.

A decade’s worth of work using the novel In Situ Ichthyoplankton Imaging System (ISIIS) has partially addressed this knowledge gap. The ISIIS’s 13 cm by 13 cm field of view and 50-cm depth of field allows undisturbed imaging of an array of plankton taxa, including delicate gelatinous zooplankton. ISIIS is capable of imaging particles 200 μm to ~20 cm in size, at a rate of 140 liters/s and per image sampling volume of ≈6 liters. The resulting images have a pixel resolution of 66 μm. (). Initial research using this system demonstrated habitat partitioning across the shelf-slope front by gelatinous plankton () and pelagic tunicates and larval fishes (). Subsequent studies revealed the importance of phytoplankton thin layers, internal waves, hypoxia, and lateral salinity gradients on submesoscale (≈3 to 70 km) patchiness of larval fishes, copepods, chaetognaths, shrimps, and gelatinous plankton (, , ). Through the resolution of submesoscale spatial patterns across a suite of trophically linked, taxonomic groups (including rare and delicate ones), this collective work has set the stage to answer questions of patchiness in a broad range of plankton in contrasting marine environments.

An important contrast in the world’s oceans is continental versus oceanic habitats. On average, continental marine habitats tend to be more physically dynamic, have higher nutrient concentrations and turbidity, and support greater levels of biological productivity than adjacent, oceanic waters (). The abundance, biomass, taxonomic, and size composition of plankton guilds (e.g., phytoplankton and zooplankton) reflect these environmental differences. Broadly, plankton in continental habitats are more abundant, have larger body sizes, and assemblages have lower diversity from reduced evenness (i.e., a few taxa dominate) in comparison to oceanic waters (–); notable exceptions include coastal water advected offshore, frontal zones, and cold-core eddies where nutrient concentrations tend to be elevated. These biophysical contrasts strengthen with increasing latitude as seasonality effects intensify (), yielding stronger habitat partitioning between continental and oceanic plankton assemblages.

To quantify plankton patchiness, we deployed the ISIIS in 2014 and 2015 in both continental and oceanic regions of the Straits of Florida (SOF). Flow through the SOF is dominated by the Florida Current that enters the southwestern SOF as part of the Loop Current from the Gulf of Mexico (GOM), with maximal flow along the western boundary. Dynamically, this flow pattern creates mesoscale cross-strait and coarse-scale vertical gradients in the flow field () onto which temperature, salinity, and the biological environment map onto. Specifically, the isopycnals tilt upward on the western edge, creating shallow, dense, cool, salty water and an upwelling of nutrients (). Consequently, the eastern (oceanic) and western (continental) areas of the SOF are distinct environmentally, and this difference is reflected in the regional distribution, abundance, and growth of the planktonic community (, ). Towing the ISIIS at 5 knots enabled the sampling of large distances rather quickly to limit aliasing. We towed ISIIS at discrete depths (15 m, shallow; 30 m, mid-depth; and 50 m, deep in 2014; and 15 and 30 m in 2015) along an average transect length of 10.8 km (± 2.6 km; 1 SD) in continental (n = 21 transects) and oceanic (n = 18 transects) areas. ISIIS imagery was automatically classified using a spatially sparse convolutional neural network trained to sort vignettes (). These analyses were performed to identify (i) the frequency of patchiness among plankton taxa, (ii) whether patch characteristics are related predictably to organism size and how patch size varies among plankton taxa, (iii) the dominant length scale of patches, and (iv) whether relationships exist between taxa-specific patch occurrence, patch size, and abundance.

RESULTS

More than 32.2 million vignettes (i.e., individual particle images) were segmented from 104,125,200 TIFF (tagged imaged file format) images (357 terabytes of imagery data). Vignettes were sorted into 125 classes at an error rate of <5%. Classified vignettes were assigned to one of 40 taxonomic groups (table S1). Aggregations of living plankton groups (n = 36, excluding artifacts, detritus, fecal pellets, and unknown) were identified using a modified distance to next encounter (DNE) method (, ) because each transect’s vignettes (i.e., individual organisms) were sequential observations in horizontal space. A classified individual organism was assigned to an aggregation if the maximum distance between it and its nearest neighbor belonging to the same taxon was less than or equal to a threshold value, Dmax (table S1). Each taxon’s Dmax was found by comparing a theoretical, random distribution of taxon-specific DNE values to all observed taxon-specific DNE values (details in Methods). Aggregations were designated as patches if the density of individuals per cubic meter in the aggregation was greater than or equal to the mean + 95% confidence interval (CI) of SOF background densities (table S1). Organism size was expressed as equivalent spherical diameter (ESD; in millimeters).

Patch frequency

Patchiness was highly variable among plankton taxa (Fig. 1). There was no evidence of a phylogenetic pattern in patch frequency across a diverse set of planktonic taxa ranging from cyanobacteria to larval fishes. Patterns of overall horizontal plankton patchiness varied more between continental and oceanic waters than among the three fixed depths (15, 30, and 50 m). In both years, plankton patches (with all taxa combined) were significantly more frequent in continental versus oceanic waters (fig. S1; χ22014 = 12.59, df = 2, P = 0.002; χ22015 = 84.50, df = 1, P < 0.001). This pattern was especially evident in shallow waters in 2014 (fig. S1). Patches of plankton were also more frequent in oceanic, mid-depth waters than expected by chance in 2015 (P < 0.0001). Parcels of water with plankton randomly distributed (i.e., where DNE between taxonomic-specific individuals was greater than Dmax) were more frequent in oceanic, mid-depth waters in 2014 (χ22014 = 10.32, df = 2, P = 0.006) than expected by chance. However, random plankton parcels were distributed evenly across the SOF and with water depth in 2015 (fig. S1).

Fig. 1.

Position and length of plankton patches for 36 taxa.

The position and length of patch (colored bars with black outline) and aggregated parcel (opaque colored bars) for each taxonomic group along a single example set of 10-km ISIIS transects in the shallow and mid-depth waters of continental and oceanic regions of the SOF. The bar color gradient indicates the log10-transformed density of individuals (individuals per cubic meter). Multiple groups had patches or aggregated parcels that extended the entirety of the 10-km transects (e.g., diatom).

Patch size

Taxonomic-specific patch length was inversely related to organism size (ESD) in both continental and oceanic areas (Fig. 2A), meaning the smallest organisms had the largest patches. The size scaling relationships were best fitted with power law functions (Fig. 2A). The continental (kcontinental = −3.27) and oceanic (koceanic = −4.08) slopes were statistically similar. Thus, the relative change in log10 of plankton patch length versus log10 of planktonic organism ESD in the SOF can be described using weighted mean slope of kSOF = −3.64. This inverse relationship was somewhat reflected in the taxonomic-specific dominant length scale (i.e., the patch length mode), although it varied considerably across taxonomic groups, ranging from 0.6 to 244 m (Fig. 2B). Like plankton patch frequency, there was no phylogenetic pattern in the dominant length scale of patches. Eleven taxonomic groups regularly had patches exceeding 10 km in size—our standard transect length (Fig. 2B). These groups were among those with the highest mean abundances (table S4). As a result, the dominant length scale of these 11 groups exceeded 10 km, and the patch size mode was not calculated because the true length of these long patches was beyond our sampling range.

Fig. 2.

Plankton patch length scaled negatively with organism size.

(A) In the SOF continental (green cross) and oceanic (blue circle) regions, mean patch length (m) increased as planktonic organism size (ESD, mm) decreased. Organism size explained 21 and 30% of variability in mean patch length in continental and oceanic waters, respectively. Shaded regions represent the 95% CI. (B) Dominant horizontal patch length scale. The dominant patch length scale for plankton groups for which we could accurately estimate the largest mode (value to right of each bar). Parenthetical values are each plankton group’s mean ESD (mm).

The length scale of patchiness varied in both horizontal and vertical dimensions and between years (Fig. 3 and table S2). Horizontal contrasts were stronger than vertical ones (Fig. 3), with differences accentuated in 2014. Patches were larger in continental (mean ± 1 SD, 2014: 188.5 ± 879.2 m; 2015: 153.9 ± 715.0 m) versus oceanic waters (mean ± 1 SD, 2014: 102.7 ± 640.0 m; 2015: 56.1 ± 313.9 m) when considering all plankton patches [analysis of variance (ANOVA) linear mixed effects, F1,13600 = 169.9, P < 0.0001, Akaike information criterion = 26605.3). This cross-strait contrast was significant for taxa such as larval fish, physonect siphonophores, polychaetes, pteropods, decapod shrimps, and other crustaceans.

Fig. 3.

Distribution of patch lengths.

Distribution of patch lengths of representative planktonic taxa in the SOF continental and oceanic waters. Blue shaded boxplots are for 2014, and red shaded boxplots are for 2015. Colors darken with increasing depth (note absence of deep data in 2015).

Multiple groups exhibited considerable variance in patch length among water depths—patch length sometimes differed by an order of magnitude (Fig. 3 and fig. S2). Patch length declined with increasing water depth, regardless of habitat, for larval fish, physonect siphonophores, decapod shrimps, sea star larvae, and narcomedusae (table S3). Only patches of radiolarian and urchin larvae increased in size with water depth. Vertical gradients in patch size tended to be steepest in continental waters in 2014 compared with 2015, although between-year comparisons were constrained by the absence of “deep” transects in 2015. Overall plankton patch size decreased with increasing water depth in continental waters, but not in oceanic waters (table S3). These vertical gradients were likely driven by significant variability identified in continental waters for taxa such as larval fish, physonect siphonophores, polychaetes, decapod shrimps, and other crustaceans (table S3).

Patch size differed between 2014 and 2015 for 12 groups. Among them, larval fish, salps, physonect siphonophores, polychaetes, other crustaceans, pelagic tunicates, urchin larvae, sea star larvae, tornaria larvae, and radiolarian patches were larger in 2014 than in 2015 (Fig. 3 and fig. S2). The mean and variance in patch size for other groups such as cydippid ctenophores, cestid ctenophores, and other shrimps were relatively consistent in space and time (Fig. 3 and fig. S2). In contrast, patch size for salps, doliolids, and narcomedusae was highly variable across time and space.

Plankton patch and background densities

Our patch identification method inherently meant patches would be more concentrated than background abundances for a given planktonic taxon. However, this method only sets a minimum threshold; patches were far more densely packed in many instances (Fig. 4). Taxon-specific plankton patch concentrations were significantly higher than background densities (Fig. 4 and table S4). This pattern was highly consistent across taxonomic groups, with stronger differences for taxa such as euphausiids and cydippid ctenophores (Fig. 4). Multiple plankton patches had densities in the top 75th quantile. In some cases, these intense concentrations of plankton extended over an order of magnitude (Fig. 4).

Fig. 4.

Distribution taxon-specific patch and background densities.

The distribution of taxon-specific density (individuals per cubic meter) within patches and over transects (i.e., background) in the SOF. Transect-level densities were calculated using the volume of the entire transect sampled. The dashed red line is the taxonomic Dmax value (see table S1 for values for all taxa). Note that the y axis scale is logarithmic (base 10). A complete list of t test statistics is available in table S5.

Overall “patch intensity” was another metric considered. This measures the degree of spatial heterogeneity experienced by an individual. It was estimated as the ratio of patch mean density to background mean density. Spatial heterogeneity (i.e., “patchiness”) becomes more acute as patch intensity exceeds a value of one, which would indicate the density of individuals in patches and the background are equal. Patch intensity values varied considerably within and among taxa in the SOF continental and oceanic waters. Among the 29 taxa that had patches in both areas, 90% had higher patch intensities in oceanic waters (table S4). The three exceptions to this pattern were tornaria larvae, urchin larvae, and calanoid copepods that had greater patch intensity in continental waters. The difference in patch intensities between oceanic and continental waters was notable for multiple taxa, many of which were gelatinous. These taxa included siphonophores (4017× oceanic versus 609× continental), echinoderm larvae (689× versus 285×), pelagic tunicates (203× versus 32×), salps (98× versus 9.7×), and narcomedusae (19× versus 1.8×). Other groups also had greater patch intensities in oceanic waters, albeit with more muted contrasts. Cydippid ctenophore, calycophoran siphonophore, and other copepods patch intensities were similar in oceanic and continental (table S4). A few taxa only had patches identified in continental waters—the greatest patch intensity among those was Copilia spp. copepods at 836×.

The total number of patches of a particular plankton taxon increased monotonically with mean density (individuals per cubic meter) on a given transect to a threshold point and then subsequently declined. This pattern held across multiple taxonomic groups (Fig. 5), although for some groups, we captured only the increasing portion of this negative parabolic curve (e.g., cydippid and cestid ctenophores) and, for other taxa, the decreasing portion of the curve (e.g., larval fishes and pelagic tunicates). We hypothesize that the patch count generated a parabolic curve due to fewer (and longer) patches at higher densities (Fig. 5). Groups not shown on Figure 5 that also exhibited a significant, nonlinear relationship between patch count and background density included other shrimps, tornaria, other echinoderms, urchin larvae, other jellyfish, and doliolids (table S6).

Fig. 5.

Taxon-specific relationships between patch count and background density.

Taxon-specific relationships between the total number of patches and background density (indiv. m−3) for each transect sampled in the SOF in 2014 and 2015. These were most often parabolic and best fit with a second- or third-degree polynomial model. Circle size represents each taxa’s mean patch size (m) on a given transect. Mean patch size was scaled for each taxon, so the relative difference in point size is valid within individual plots. Circle color represents the mean density of organisms in patches.

DISCUSSION

A fundamental step in testing mechanistic hypotheses about factors affecting plankter fitness, behavior, community structure, and trophodynamics is determining whether a given phenomenon occurs across a wide range of scales or whether it is limited in ambit [sensu Levin ()]. Only after answering this question of scale can causal mechanisms generating spatial or temporal patterns be resolved. The patchy structure of plankton in the horizontal dimension, first postulated by Haeckel (), demonstrated by Hardy (), and synthesized by Hutchinson (), has been extensively measured in situ at scales exceeding 1 km for a diverse suite of taxa (). At these scales, potential mechanisms generating and maintaining marine plankton patches include reproductive (i.e., egg laying in the same area), vectorial (e.g., nutrient and salinity gradients), and stochastic (i.e., large, random events), sensu Hutchinson (). At spatial scales <1 km (micro-, fine-, and coarse scale), social (e.g., signaling and forage facilitation) and coactive (e.g., predation and niche partitioning) processes are expected to dominate.

Earlier efforts to quantify fine-scale plankton patchiness in the horizontal dimension often either fit data to theoretical models derived from a particular spatial distribution or analyzed its spatial variance (–). While these approaches succeeded in revealing the existence of spatial heterogeneity and patchiness in plankton, they overlooked the structure, taxonomic diversity, and precise spatial position of individual patches ()—parameters needed for identifying potential mechanisms driving aggregations. Only Wiebe () described fine-scale mean patch structure, and that effort focused on a few crustacean and pteropod species in a single habitat. We bridged that gap here through unobtrusive in situ sampling and automated classification of 32 million planktonic organisms. This approach allowed us to map the exact spatial position of each plankter in horizontal space simultaneously, identify whether they were in a patch, and estimate the structure of that patch. The dominant patch length scale we could conclusively measure that ranged from 0.6 to 244 m fell within the microscale (0.01 to 1 m) and fine scale (1 to 100 m). For each taxon, we quantified patch length, frequency of occurrence, and organism concentration at different depths in the same water masses—both those in proximity to the continental shelf and those more oceanic. Submesoscale currents, while potentially influential to ecological interactions that can generate spatial heterogeneity (), were assumed to be of secondary importance relative to dominant flow of the Florida Current. Taxa with defined patches ranged in size and rarity from small, highly abundant acantharia protists and Oithona spp. copepods to rarer larval fishes and large cestid ctenophores.

Size scaling of patch length

Patch lengths varied by several orders of magnitude within and among taxa. Notably, we found a significant and negative relationship between mean patch length and organism length, with a steeper negative slope (−4.08) in the oceanic compared with the continental environment (−3.27). Assuming a power of 3 relationship between carbon biomass and ESD [a reasonable assumption even considering the differences between crustacean and gelatinous zooplankton ()], this implies a negative sublinear to negative linear (−0.27 to −1.08) scaling between patch length and body mass. We also found that mean patch length was not only smaller in oceanic environments but also decreased faster relative to organism size than in coastal environments.

The underlying drivers of plankton patch formation can be quite variable and can range from finding prey, avoiding predators, and locating mates (). We had two a priori hypotheses: that patch length would scale positively with size if identifying favorable feeding environments were the dominant mechanism. There are numerous studies that suggest that animal home ranges and thus foraging ranges scale positively with size (, , ). Alternatively, patch length would scale negatively with size if aggregation were driven by predator avoidance (e.g., plankton form patches to achieve safety in numbers or to appear larger to a potential predator). We found a significant negative size scaling relationship for patch length. While this does not give conclusive evidence for one driver over the other in determining patch size, our results imply that mechanisms that include increasing density or effective size to protect against predation are more likely than prey searching for the formation of plankton patches.

Taxonomic variations in patch size

Previous size estimates of fine-scale patches are available for only five taxa (e.g., larval fishes, pteropods, euphausiids, shrimps, and gelatinous planktons) in our diverse assemblage of 36 groups. In subtropical, northern GOM (nGOM) waters, patches of larval fishes ranged from 37 to 100 m (). In that study, the size of larval fish patches scaled with background abundance, i.e., larger patches occurred in areas with more larval fish. In the SOF, larval fish patch size also tended to increase with background densities. Median patch lengths ranged from 78 to 333 m in continental waters, where larval fish abundances were greater overall. In contrast, the median length of larval fish patches in oceanic waters (where background densities were lower) was shorter, ranging from 42 to 85 m. Pteropod and euphausiid patch length scales in the present study were in range of Wiebe’s () daytime patch radii for those taxa (pteropod: 14.4 m; euphausiid: 13.6 to15.0 m). The dominant patch length scales of SOF shrimp (12.3 to 18.9 m) and gelatinous plankton (8.7 to 244.1 m) were also coherent with Greer’s et al. () estimates from the nGOM (shrimp: 12 to 32 m; gelatinous plankton: 5 to 116 m). However, chaetognath patches in the SOF were 2000× larger than those in the nGOM (). One possible reason for this difference could be the Dmax thresholds (table S1); ours was 23.5 m for chaetognaths, while that of Greer et al. () was <1.5 m.

We could not accurately estimate the patch size for several taxonomic groups of plankton, including protists, diatoms, and calanoid, and cyclopoid Oithona spp. copepods, because their patches spanned the entirety of our 10-km transects. While fine-scale patchiness in phytoplankton in horizontal space is possible (), SOF diatom and protist patches exceeded 10 km in almost all cases. The scale of physical processes in this region (e.g., continental upwelling, passage of mesoscale eddies, and fast-moving western boundary current) may spread patches over a larger horizontal distance. Nonetheless, the length scales for phytoplankton and copepod aggregations were much larger than our study’s fine-scale focus. Longer sampling transects (50 to 100 km) are needed to elucidate the true patch length for copepods and phytoplankton groups and to determine whether the horizontal length scales of these trophically linked taxa are concordant in the SOF. Another consideration potentially affecting our estimates of patch length is if we happened to sample patches along a pycnocline. By sampling at a constant water depth, we would have entered and exited this continuous patch as its vertical position varied with the depth of the pycnocline. In those cases, our estimates of patch length in the horizontal dimension for that transect would be underestimated and patch count overestimated.

The density of organisms in patches (i.e., individuals per cubic meter) was significantly greater than the average background density for the same taxa. While this is not unexpected because of our definition of “patch,” in all cases, mean patch density was higher than the defined level, and in many cases, the patch density was more than 10× greater than the background density. These results are quantitative evidence that the mean background density of plankton measured by sampling gear that integrates over the water column (i.e., nets) does not represent the actual biological environment an individual plankter experiences. However, the distributions of plankters, their prey, and predators are extremely heterogeneous in both horizontal and vertical space. These densely packed horizontal plankton patches explain, in part, the success of critical ecological processes such as individual growth and reproduction otherwise not sustainable in a low productivity environment with homogenously distributed individuals ().

To the best of our knowledge, only two previous studies have estimated patch intensity for plankton in horizontal dimension. In the Svalbard polar front, patch intensities of calanoid copepods, Oithona spp. copepods, hydromedusae, chaetognaths, appendicularians, and meroplankton ranged from 3.2× to 16.7× (). SOF patch intensity estimates for these same groups tended to be lower (e.g., calanoid copepods, 2×). Only SOF relative patch intensities of Oithona spp. (2×) and hydromedusae (3×) were comparable. Because physical drivers, swimming ability, and locating prey play a role in patch formation and maintenance (), it is logical that the Svalbard front had higher intensity plankton patches compared with the subtropical SOF. Except for large copepods, Svalbard zooplanktons were concentrated at the frontal physical discontinuity (). They posit that the relatively enhanced swimming ability of large copepods [swimming velocity scales with organism mass ()] compared with other smaller taxa allowed the copepods to form biologically cued aggregations. While the SOF is physically dynamic, in this study, we did not sample distinct submesoscale features such as fronts or eddies that may drive patch formation. Wiebe () also calculated this metric for patches of a tropical pteropod (3.3×), chaetognath (3.1×), a poecilostomatoid copepod (2.6×), and three euphausiid species (3.1× to 5.1×). Relative patch intensity of pteropod and chaetognath in the SOF were marginally lower than his estimates, but the patch intensities for SOF Copilia spp. (835×) poecilostomatoid copepods and euphausiids (17×) were markedly higher than Wiebe’s values.

Comparison of the number of patches of a given plankton taxon on a transect to the mean density of that taxon over the entire transect yielded an interesting and consistent pattern. For multiple taxa, including decapod shrimp, pteropods, and physonect siphonphores, the relationship between frequency of patches on a given transect and the mean transect density was nonlinear. While initially increasing with mean transect density, in most instances, patch counts declined after reaching an inflection or vertex point (Fig. 5). This taxon-specific inflection point combined with increasing patch lengths after the inflection point indicate that there is a mean abundance at which smaller patches begin to merge into larger (but fewer) patches. As the patch length and corresponding patch volume increased, the densities of individuals within patches tended to decline.

Four taxonomic groups (i.e., cydippid ctenophores, cestid ctenophores, pelagic tunicates, and larval fishes; not all shown in Fig. 5) did not follow this parabolic pattern and instead exhibited significant linear relationships. For cydippid ctenophores, cestid ctenophores, and pelagic tunicates, the patch count rose consistently with background density. In contrast, the number of larval fish patches declined with increasing background densities. We considered three interpretations of these exceptions. First, mechanisms driving these patterns could differ among taxa, and the linear relationships may be expressions of that variability. This interpretation seems unlikely given the consistency of the pattern among the 14 other taxonomically diverse groups. A second hypothesis is that we only captured a portion of a response curve due to a sampling limitation (e.g., not sampling a large enough range in background densities). A third possibility is that the response of those four taxa was constrained by an unknown mechanism. The second and third hypotheses are more consistent with the nonlinear pattern exhibited by the other groups. For example, similar to physonect siphonophores and pteropods, larval fish mean patch size became larger as the number of patches declined (Fig. 5). The dominant length scales of patches for all taxa with significant relationships (including the divergent four) are at the fine scale and on the same order of magnitude (Fig. 2B). A reproductive driver is unlikely given marked variability in reproductive and life history strategies among the 18 taxonomic groups for which we identified a pattern. Meroplankton (e.g., larval fish) and holoplankton (e.g., cydippid ctenophore) are represented, as are taxa that have alternating sexual and asexual reproductive stages (e.g., salps). Many species within these groups also consume primary production. The fine scale at which patchiness was observed may be due to coactive or social processes (e.g., signaling, predation; competition, response to prey) as hypothesized by Hutchinson ().

Continental versus oceanic patch structure

Our results demonstrate that while patch metrics varied greatly among plankton taxa at fine- to coarse scales, plankton patches were more frequent and larger in continental than in oceanic waters. A greater number of patches in an environment where plankton are more abundant is somewhat intuitive, as more individuals are available to aggregate, although such relationships are not necessarily the case as substantially higher abundances could potentially contribute to a higher uniformity. Of the 36 plankton groups analyzed, 91% were more abundant in continental waters relative to oceanic areas.

The continental versus oceanic contrast in patch frequency and size suggests that the drivers of plankton patchiness, such as physical gradients and discontinuities, and the interplay between those conditions and taxon-specific individual organism behaviors are stronger in more physically dynamic waters. Maximal flow along the western boundary of the Florida Current generates horizontal and vertical structure in the flow field; specifically, the isopycnals tilt upward on the western edge, creating an upwelling of cooler, salty water and nutrients (). This physical structure and associated gradients in vertical current shear and nutrient concentrations stimulates phytoplankton production and, in turn, the aggregation of zooplankton consumers optimizing foraging opportunities (, , ). Enhanced phytoplankton productivity via upwelled nutrients would attract a multitude of zooplankton grazers sampled in our study, including calanoid and cyclopoid Oithona spp. copepods, shrimp, and euphausiids.

Localized vertical shear can also aggregate large gelatinous predators such as scyphozoan medusa when individuals swim into the current to increase prey encounter rate (). In contrast, smaller gelatinous predators such as ctenophores actively avoid hydrodynamic mixing (). This divergent behavioral response may explain why patches of physonect siphonophores (10- to 20-cm total length) were generally longer and more variable in continental waters than cydippid ctenophores [0.5- to 1-cm diameter; ()].

As anthropogenic climate change continues to alter the physical nature of continental and oceanic environments in ways directly affecting planktonic systems, including temperature (), acidification (), and upper water-column stratification (), the ecological consequences of these physical changes are numerous. The upper ocean ecosystem has experienced reduced net primary productivity, mesozooplankton biomass, and fish abundance (, ). Many temperate, subtropical, and tropical marine ecosystems have exhibited range and phenological shifts, as well as assemblages characterized by more smaller, warm-water species (), with the further expansion of the subtropics projected under climate change (). One implication of the size scaling of patch length that we observed is that under climate change, patch length would be expected to decrease as organism size decreases (). This change would have broad implications for overall energy expenditure and foraging efficiency for higher trophic-level predators such as fish, seabirds, and marine mammals.

Conclusion

Patch frequency, density, and size in the horizontal dimension varied greatly among the 36 plankton taxa sampled in contrasting subtropical continental and oceanic waters. The dominant length scale of patches was at the microscale (0.1 to 1 m) and fine scale (1 to 100 m), with most length scales between 10 and 30 m. We found a significant negative size scaling relationship for patch length in both continental and oceanic waters, suggesting that increasing aggregation size to protect against predation was more likely than prey searching for patch formation. Patches tended to be more frequent and larger in continental than in oceanic waters for multiple taxa such as larval fishes, physonect siphonophores, and pteropods. For some other taxa (e.g., euphausiids and cydippid ctenophores) differences were not significant between water masses. The density of individuals within patches was 1.2× to 4017× greater than the background density of the same taxon. Patches of Copilia spp. copepods, euphausiids, salps, and cestid ctenophores were especially concentrated. These results clearly demonstrate that the mean background density of a plankton group does not represent the actual biological environment an individual plankter experiences. Last, the relationship between total patch count and mean background density is consistently nonlinear for multiple plankton taxa, with many exhibiting inflection points, after which patch counts declined with increased background densities. Our empirical quantification of in situ patch structure for individual taxonomic groups at the micro and fine scale sets the stage to ask future, advanced ecological questions about coactive drivers, including taxonomic groups coexisting in the same patch space, predator-prey interactions, competition, and niche partitioning. Further, this deeper understanding of plankton patch structure will help refine complex biogeochemical models to predict ecological and population responses to changing oceanographic conditions.

MATERIALS AND METHODS

Experimental design

Plankton were sampled during 28 May to 5 June 2014 and 18 to 26 June 2015 in continental and oceanic waters of the SOF using an ISIIS (, ). Because the ISIIS is towed at 5 knots, a large area can be sampled relatively quickly, limiting aliasing. A drogue was deployed at sunrise each day before sampling, and ISIIS tows were conducted within sight of the drogue to ensure the same water mass was sampled during each tow. We towed ISIIS behind the R/V Walton Smith at discrete depths (15 m, shallow; 30 m, mid-depth; and 50 m, deep in 2014; 15 and 30 m in 2015) along an average transect length of 10.8 km (±2.6 km, 1 SD) in continental (n = 21 transects) and oceanic (n = 18 transects) areas. ISIIS vignettes were segmented and classified using an automated spatially sparse convolutional neural network (sCNN) trained to sort vignettes ().

Image processing

Training library

We used the same sCNN classifier training library as in (). Training set development was an iterative process where experts sorted plankton vignettes into classes, trained the sCNN, tested the classifier on the same set of vignettes, and then manually resorted vignettes or added classes based on the performance of the classifier (e.g., creating three classes of appendicularians based on tail shape). Most classes had hundreds or thousands of vignettes, with a few containing 20 to 100 (National Centers for Environmental Information, doi: 10.7289/v5d21vjd).

Classes were added or modified after each training-testing iteration if confusion matrix metrics indicated consistent classifier errors. The training library was finalized after the classifier performance reached an error rate <5% for 400 epochs. The final library had 61,445 vignettes in 125 classes. Classes were assigned to broader taxonomic groups and confusion matrix metrics estimated ().

Class predictions

The sCNN classifier generated for each plankton vignette and class combination a probability value ranging from 0 (least likely) to 1 (most likely). Because these values represented the likelihood a vignette belonged to that class, vignettes were assigned to the class with the maximum probability. To maintain an overall precision of 90% for a taxonomic group, group- and year-specific probability filtering thresholds were applied (, ). Filtering thresholds were identified by (i) manually classifying vignettes (separate from the training library and confusion matrix sets) drawn randomly from each sampling year (2014: 71,969, 2015: 66,405), (ii) iteratively setting each classification score within the group’s vignettes as a threshold, and (iii) calculating the precision for vignettes with scores greater than the threshold. A Loess curve was fitted through these precision score threshold data. The classification score corresponding to a 90% precision was set as the group’s filtering threshold (). Vignettes with scores less than the filtering threshold value were reclassified as unknown. A count of one was then applied to each classified vignette. True zeros for each group (i.e., frames with no or all unknown plankton vignettes) were included.

Confusion matrix

The overall ISIIS confusion matrix was generated from 143,418 vignettes randomly sampled from data collected in 2014 and 2015 in the SOF (2014: 48,633 vignettes, 2015: 94,785 vignettes. Sampling years were combined to provide the largest number of images to evaluate the sCNN classifier’s performance. Because of the random sampling approach, representation of taxa in the confusion matrix sample was proportional to their natural frequency of occurrence. Each vignette was manually identified by one human expert and then reviewed by a second to confirm the classification. The confusion matrix and associated metrics of classifier performance (e.g., precision) were generated at the group level using the validated vignettes after applying filtering thresholds. The average precision was 89.0% across all groups; mean recall (number of correct positive results/number of positive results that should have been returned) was 50%. The F1 metric, a metric of the classifier’s accuracy, was 58.8%. Classifier performance varied across taxonomic groups (). Precision for larval fish, a rare but key group, was high (88.2%); however, recall (41.2%) was low, yielding an F1 score of 56.2%. Diatoms, appendicularians, tornaria larvae, and Oithona spp. copepods all had F1 values >90% (). The correction factor for each group (precision/recall) was calculated for each group and applied to scale the number of individuals identified in each three-frame bin.

Distance to next encounter

Individual observations were binned for each plankton group by finding the total number of organisms in three, sequential frames representing 24.6 liters of water. This binning step was necessary to perform the DNE analysis as the variance with the data at the original single-frame frequency was too high.

The distance between each pair of sampling points was calculated using points within 10 m of target depths (15, 30, and 50 m). Despite this range, variation from targets was minimal (mean ± 1 SD, 2014: shallow, 14.9 ± 2.0 m; mid, 30.1 ± 1.8 m; and deep, 50.4 ± 4.3 m; 2015: shallow, 15.1 ± 0.7 m and mid, 29.9 ± 2.2 m). Splitting the data by group, sampling points ordered by time were binned every 60 s, and the distance (m) traveled between each bin was obtained using the Law of Cosines, which allowed us to account for the horizontal (x) and vertical change (i.e., depth, z) in vehicle position (Eq. 1). The mean distance traveled between sampling points in each bin was then calculated. Mean distances were added cumulatively along the entire transect to find each sampling point’s distance from start (m). This is a modification of existing DNE methods (, )

The maximum distance allowed between individuals to be assigned to an aggregation was set by a threshold value (Dmax). The Dmax was found for each group by comparing a theoretical, random distribution of DNE values to observed DNE values. The group-specific theoretical DNE distribution was estimated by generating 1000 random distributions (histogram bin size = 0.5 m) for each group-transect combination assuming a uniform distribution using presence-only observations and setting the range as the DNE minimum and maximum. A grand theoretical distribution was calculated by finding the mean and lower 95% CI of counts for each bin. This grand distribution was merged by bin with mean real DNE counts (rounded to the nearest integer) from all transects. The smallest bin with mean real DNE counts equal to zero and significantly lower than the theoretical lower 95% CI count was denoted as the group’s Dmax (table S1).

Patch definition

Plankton aggregations were identified for each taxonomic group in each three-frame binned sampling point (hereafter, sampling point) using the modified DNE analysis. Transect sections where individuals were highly aggregated (“patches”) were identified as those where the DNE of sampling points’ less than or equal to Dmax, the group’s within-patch distance threshold, and where a minimum density (individuals per cubic meter) threshold was exceeded (table S1). This threshold was defined as the 95% CI above a group’s grand mean density (i.e., across the entire study; table S1). Breaks between patches were identified as locations where the distance between a pair of sequential points (in patches) was greater than or equal to Dmax. A distinction is made between sections of sampling points where organisms were present but with DNEs ≥Dmax (“random” parcels) and those with no organisms at all (“zero” parcels).

Statistical analyses

All analyses were performed in “R” () (R Core 2020).

Patch frequency

We tested contrasts in the weighted mean patch counts of “living” plankton groups between continental and oceanic areas using the “chisq.test” function from the “stats” package. A weighted patch count was used because the sampling effort in these two areas was slightly unbalanced during 2014 and 2015. The weighting factor was the quotient of the total distance assuming equal effort and the actual total distance sampled. For 2014, the equal effort distance was 124.7 km, and actual distances were 125.2 km (continental) and 124.3 km (oceanic), yielding weighting factors of 0.997 and 1.003. For 2015, the equal effort distance was 85.4 km, and actual distances were 105.6 km (continental) and 65.2 km (oceanic), yielding weighting factors of 0.809 and 1.309. The total patch counts for each area and year were multiplied by their respective weighting factor. Patch frequency was considered significantly different than random at α = 0.05. We used the “chisq.posthoc.test” function with the “Bonferroni” method for P value adjustment for post hoc, pairwise comparisons among classes (i.e., area and depth).

Patch length

We tested for relationships between mean patch length and organism size using the mean ESD of the taxonomic groups. ESD was determined by measuring the area of the largest particle in each vignette using the Python 2.7 OpenCV library and then converting to diameter assuming a perfect circle. Vignettes were obtained during the same 2014 sampling period and region as the present study but from a separate 0- to 50-m undulation transect. Despite the ESD data coming from a separate transect, the imaging data were processed in the same way as in the present study. In addition, the number of vignettes used for ESD measurements totaled 24 million, giving us high confidence that the mean group-level ESD measurements are representative of those of the present study.

We tested for depth and area contrasts in patch length for individual taxonomic groups (e.g., larval fishes) using linear mixed effects models (“lme” function from the “nlme” package). Models were fit by maximizing the restricted log likelihood. Untransformed patch length data were checked for normality using Shapiro’s test (“shapiro.test”) from the stats package, and heteroscedasticity by examining the model residuals. In all cases, patch length data needed to be log10 transformed to improve normality. The fixed main factors were depth and area, and year was set as a random factor (random = ~1|year). lme models were also initially fit with an interaction term for depth and area. We used an ANOVA function for linear mixed effect models (“anova.lme”) from the nlme package to evaluate the significance of the fixed and random terms to patch length. The interaction term was removed if it were nonsignificant (P > 0.05), and the model was refit using only the fixed main factors of depth and area and the random year factor. If the overall model were significant (P < 0.05), then a Tukey post hoc test was performed to evaluate the effect of area or depth using the general linear hypothesis “glht” function from the “multcomp” package. The Bonferroni method was used to adjust the post hoc test values. The effect of the random factor year was tested by comparing models with and without the year factor using the “anova” function from the stats package. Models without the random year factor were fit using the generalized least squares “gls” function from the nlme package. Ten taxa had insufficient sample sizes to perform statistics: siphonophore other (n = 4), calanoid copepod (n = 3), diatom (n = 5), protist (n = 4), trichodesmium (n = 3), appendicularian (n = 5), chaetognath (n = 4), Oithona spp. copepods (n = 2), lobate ctenophore, and acantharian protist (n = 4). Copilia spp. copepod and echinoderm larvae patch data were too unbalanced.

Patch density and patch count versus background densities

Patch intensity was measured as the ratio of patch mean density to background mean density in a domain (e.g., a transect). We compared patch versus background densities (individuals per cubic meter) for each taxonomic group using a t test (“t.test” function from the stats package of the mean patch and transect density for each transect).

The relationship between total patch count along a transect and background density for each taxonomic group was evaluated using the “lm” function from the stats package. Linear, quadratic, cubic, and quartic polynomial models were fit to the data and compared using the “anova” function from the “stats” package and the models’ R2 values.