UVP5 data ( Figure 1 ) were compiled from all oceans, covering a 10 year period (2008-2018). A detailed description of the operation of the UVP5 is given in Picheral et al. (2010). All particles large than ≈ 100 μm in Equivalent Spherical Diameter (ESD) were measured and counted, but only images of particles (zooplankton and aggregates) larger than ≈ 600 μm ESD were kept by the UVP5 for further processing because smaller objects contained too few pixels to be identifiable. Acquisition of metadata (geographic location, date, etc.) and processing of all 8.46 million images (95% being detritus) were carried out by the ZooProcess software which provided information on 42 morphological features associated with each object (area, major and minor axis, etc.). The results were imported into EcoTaxa (Picheral et al., 2017), an application which allows a taxonomic classification of images via supervised learning algorithms, followed by manual validation (Irisson et al., 2022). As 61% of the profiles have a maximum depth ≤500 m, only images of organisms between 0-500 m were kept and the overall estimates of biomass were restricted to this depth range. To ensure that profiles were representative, a filter was also applied to only keep profiles that covered at least 80% of the layer of interest.

Living organisms were separated from detritus (aggregates, fibers, fecal pellets) as well as artifacts (e.g. bubbles) and classified according to their taxonomic identity. Recognition and sorting of organisms can be a source of bias depending on the levels of perception and experience of the people who perform them. Several cognitive factors biases such as boredom, fatigue or a classification biased towards the most used groups have been presented by Culverhouse (2007) and Culverhouse et al. (2014). To reduce the risk of poor identification, a shared UVP5 taxonomic guide was used to homogenize image sorting into 119 taxonomic groups. The image data were thereafter grouped into 25 broader taxonomic groups ( Table S1 ), and a subset of the resulting dataset was checked for homogeneity of sorting within these groups. A minimum of 51 images and a maximum of 10% of all images were extracted from each group and were independently checked after the assembly of the final data set. The maximum error or uncertainty rate per taxon was 9.8% and a vast majority of taxa were under 2.5%. We checked the classification and if accuracy was <95%, we rechecked the categories to assure proper sorting. In addition, only fully validated profiles were used for this analysis. The resulting global data set consisted of 466,872 images from 3,549 stations. Under-sampled groups with less than 500 images in the dataset which could not be used for a global study were not included in the analysis.

We computed the organisms’ size spectrum to detect the size range within which the UVP5 can be used to properly quantify their distribution. The concentration of objects in the ocean is expected to decrease with size; when this is computed as a normalized size spectrum, the relationship is expected to be linear (Forest et al., 2012). A peak in the size spectrum at the lower size range generally reflects the minimum size of efficient detection by in situ imaging while high variability in the large size range reflects the poor ability to detect rare large objects (Stemmann and Boss, 2012). With that in mind, the spectrum was linear for the size range 1.02-50 mm and organisms outside this range were not included in the analysis since large mobile fauna (including large crustaceans) are likely to be undersampled and small zooplankton organisms close to the UVP5’s threshold of detection are difficult to identify. This size range selection ensures that the data used in this study was properly quantified by the UVP5.

To avoid errors due to incorrect ellipse fits (around appendages of organisms rather than their body, ellipse fitted to non-ellipsoidal organisms, etc.), we chose the spheroid method: it is based on the area ( Table 1 ), which is more consistently measured by the image analysis performed in ZooProcess.

For Rhizaria, biovolume (mm ³) to carbon (mgC) conversions were done using factors from the literature ( Figure S1 and Table S2 ). For other groups, the conversion from individual volume to individual wet weight assumed a density of 1 g cm^-3 (Kiørboe, 2013). Then the conversion from individual wet weight to individual biomass in carbon units (mgC) was calculated using taxon-specific linear conversion factors from McConville et al. (2016); when several conversion factors were available for a taxon, their median was used for each group. To take into account differences in density of some parts of the organisms, the Appendicularia group was actually split into Appendicularia_body and Appendicularia_house, whereby the “body” group contains images with only the animal and the “house” group contains the house and the animal. For the images labeled Appendicularia_house, we used the relationship of house diameter (major axis) to Appendicularia trunk length from Lombard and Kiørboe (2010). We then converted this body size equivalent into carbon weight using the corresponding relationship from Lombard et al. (2009). For the images labeled Appendicularia_body, we converted the biovolume of the organism into carbon weight using the corresponding relationship from Lombard et al. (2009). Two groups also have been created to separate the Collodaria into solitary Collodaria and colonial Collodaria. This choice was done based on the fact that solitary Collodaria are smaller than colonial ones and have a different vertical distribution (Faillettaz et al., 2016). For solitary collodarians with a dark central capsule (subgroup of solitary Collodaria) described in Biard et al. (2016), the estimation of carbon (0.189 mgC mm ^-3) by Mansour et al. (2021) was done on the capsule of the organisms. As Zooprocess measures the area of the whole organism, we determined the ratio a r e a   w h o l e   o r g a n i s m a r e a   c e n t r a l   c a p s u l e = 0.713 and applied this factor to avoid overestimation of carbon biomass for this group. For the rest of the collodarians, the estimation of Mansour et al. (2021) was directly applied.

In this work we used boosted regression trees (BRTs) to predict the biomass of different zooplankton groups as they show different advantages over other commonly used machine learning approaches for the nature of our dataset and intended application (Elith and Graham, 2009). This ensemble method uses regression trees, models that link a response (here biomass) to predictors (environmental variables) by successive dichotomous separations (Breiman et al., 1984; Hastie et al., 2001). Regression trees automatically select the relevant explanatory variables, can deal with categorical or continuous inputs, are not sensitive to the distribution of the continuous ones, can represent relations of arbitrary form and naturally include interactions among explanatory variables (Elith et al., 2006). With so-called surrogate splits, they can also deal with missing values in the explanatory variables. They are therefore very convenient to use, but their predictive power is often limited and they have difficulties to capture smooth relationships. Boosting is a way to overcome these drawbacks (Schapire, 2003). It is based on the fact that it is easier to find many rough rules of thumb than to find a single, highly accurate prediction rule (Schapire, 2003). BRTs combine many short regression trees in succession, each new tree being adjusted to consider the observations poorly predicted by the previous ones (Elith et al., 2006; Leathwick et al., 2006; Elith et al., 2008). This improves predictive performance and the smoothness of the prediction (Leathwick et al., 2006). In addition, only a random subset of the input data is used to fit each tree and this stochastic component reduces the variance of the final model ensemble (Friedman, 2002).

Boosted regression trees (BRTs) have an ability to handle a large number of variables and - other than Generalised Linear Models (GLMs, Nelder and Wedderburn (1972)) or Generalised Additive Models (GAMs, Hastie and Tibshirani (1986); De’ath (2007); Elith et al. (2008)) - do not seek to fit one single model portraying the relationship of the response variable (here biomass) and its predictors (environmental variables). Various recent studies (González Carman et al., 2019; Chen et al., 2020; Hu et al., 2021) have compared BRTs results to other modeling tools such as GAMs, GLMs, Random Forests (RFs), Maximum Entropy modeling (Phillips et al., 2006; Elith and Graham, 2009) or neural networks and have obtained better predictive performance with BRTs. Other studies (Zhang et al., 2018; Son et al., 2018) used complementary GAMs and BRTs to study the effects of explanatory variables. However, BRTs could be slower than RFs (Chen et al., 2020) and training parameters need to be chosen carefully to avoid overfitting (Leathwick et al., 2006; Elith and Graham, 2009). BRTs were chosen over RFs because of their capacity to reduce both the bias and the variance of model results (Hastie et al., 2001). BRTs are also less sensitive to the effect of extreme outliers and the inclusion of irrelevant predictors (Leathwick et al., 2006). This makes them suitable for plankton datasets, as sometimes very high plankton biomass values do occur during blooms (Brodeur et al., 2018; Pettitt-Wade et al., 2020). BRTs also have the ability to handle sharp discontinuities which is not the case of the GAMs (Elith et al., 2008). This is important when modeling taxa which can have a narrow habitat.

In addition, in regression trees, the loss function, used to determine which dichotomous split to perform, can be changed to be adapted to the distribution of residuals. Here we explored the classic mean squared error, which assumed a somewhat normal distribution of the residuals, as well as a Tweedie loss adapted to zero-inflated data (Zhou et al., 2019), and a Poissonian loss, which considered data as discrete counts, also including many zeros. To use the Poisson loss, the biomass was scaled so that the value of the 1% quantile was ≥ 1 and then rounded to the nearest integer; the inverse scaling was performed after prediction. This later approach proved to produce the best fits and more robust models in a few test taxa and all models were therefore fitted with Poisson loss. The models and statistics were computed using the xgboost package (Chen et al., 2021) in R version 4.1.2 (R Core Team, 2021).

Individual biomass values derived from UVP5 images and environmental data measured at various layers were both averaged over a depth range of interest and matched geographically, on the 1°×1° grid. Biomass values matched to the same 1° pixel, and therefore associated to exactly the same environmental data, were averaged.

We hypothesized that an association between biomass and environment investigated at a fine scale could be more efficiently learned by the model because is contains less noise, so we divided the data vertically between the epipelagic (0-200 m) and mesopelagic (200-500 m) zones and also tried a finer partition, into 100 m depth bins between 0 and 500 m. Evaluating separate models for each layer could allow to focus on finer subgroups within our quite coarse taxonomic units (some species being mostly present in one of the layers) and therefore define biomass-habitat relationships at a finer, more relevant biological level.

For the same reason, we also built models on subsets of data partitioned geographically. Indeed, polar copepods have a different thermal niche compared to tropical ones (Rombouts et al., 2009; McGinty et al., 2021). So, in addition to a model fitted on the global dataset (world), we trained models on data from the region between 40°S and 40°N (low latitude) and from the data collected outside of this latitudinal band (high latitude). Out of the 3,549 profiles composing the UVP5 dataset, 2,837 are located between 40°S and 40°N and 712 were done outside of this latitudinal band.

For each taxon in each spatial partition, the data was split to distribute 80% of it in a training and validation sets, on which the model was fitted and assessed, and 20% to a test set, on which predictive performance was evaluated. This split was stratified according to the deciles of biomass in the data, to ensure that both the learning and test sets contained low and high biomass points.

To choose model hyperparameters (i.e. parameters of the model adjustment algorithm) and to evaluate the variability in the prediction due to the constitution of the training set, each 80% portion set was resampled through five-fold cross validation repeated 20 times [i.e. 100 resamples; (Hastie et al., 2001)]. For each cross-validation fold, the model was actually trained on four folds and validated on the last one. The splits into the five folds were also stratified according to the deciles of biomass, for the same reason invoked above.

To extract as much information from the data, while avoiding overfitting, various combinations of hyperparameters were tested for each model (Elith et al., 2006). They included: 1) the learning rate per tree determining the contribution of each tree to the ensemble model (0.05, 0.08 and 0.1 were tested); 2) the maximum depth of a tree (2, 4 and 8 were tested); 3) the minimum number of elements per leaf (which also limits the depth of the trees; 1, 3 and 5 were used); 4) the number of trees used for the prediction (values up to 600 were tested). For each combination, the model was fitted to the training set and evaluated on the validation set of each of the 100 resamples; the loss was then averaged over the 100 resamples. The best set of hyperparameters is usually the one for which this average loss is minimal. The differences around that minimum are often small and not always meaningful; to be sure to avoid overfitting, we applied an early stopping criterion whereby the increase in the number of trees was stopped when the error did not decrease by more than 1% after adding 10 trees.

Once the best set of hyperparameters had been chosen, the relevance of the corresponding model was quantified by the Pearson correlation between the observed biomass data in the test set and the predicted biomass, where prediction is the average of the predictions of the 100 models fitted to the resamples. This metric captures the model’s ability to correctly represent general trends and patterns in the data set and is one way to compute the R². The significance of this correlation can also be tested and quantified with a p-value. These metrics can be readily compared across the various spatial partitions of the data because they represent the skill of the models on an independent data set, not the quality of the fit to the training data (like the way the R² is usually computed). To compare the worldwide and regional approaches fairly, it is important to focus on the same regional subset. To this effect, two additional R² were computed for the global model: on the test data located inside the 40°S-40°N latitudinal band and on those outside of it (world low latitude and world high latitude).

To identify which environmental variables drive the change of biomass in each specific model, the percentage of variance explained by each variable was calculated as the sum of the effects of the variable at each node of each tree where it was used. To describe the shape of the effect of each variable, univariate partial dependence plots were computed as the average ± standard deviation marginal effect of the variable in the 100 resamples. Practically, the variable of interest was set at a given value at all training points, the other variables were left at their original values, the average biomass predicted over all points was computed, for each resample; then the mean and standard deviation of those averages were computed across resamples. Finally, the variable was set to another value and so on. To describe the full range of each variable, the partial dependence was estimated at 10% quantile.

To obtain global maps of predicted biomass, the regression between UVP5 biomass data and environmental variables was applied to all points in the corresponding partition of the world, in depth and space. Because 100 models were fitted to the resamples of the training data, the standard deviation of biomass among the 100 predictions (σ_b ) can be computed in addition to the mean (m_b ), and the coefficient of variation (CV), defined as C V = σ b m b , then gives an indication of the uncertainty of the model predictions.

To get a robust estimate of global zooplankton biomass in the 1.02 mm to 50 mm size range, we chose to be conservative (i.e. ad minima): only the taxonomic groups in the global partition for which the correlation between predicted and observed biomass was significant were used. The surface area of each 1°×1° cell was computed using the following formula:

with the area A in m ², the south and north latitudinal limits of the cell in radians and ℛ, the earth radius (6,378.137 km). For each group used, the biomass was integrated over the relevant layer in each 1°×1° cell by the following calculation

where b ^ is the estimated biomass in mgC.m ^-3, A in m ² is defined above, l is the layer thickness in m and therefore b ^ t is the total biomass in mgC. Finally, the global ad minima zooplankton biomass estimate was computed by adding up the biomass for all selected groups and the 0-200 and 200-500 m depth layer.

Copepoda is one of the best predicted groups in the epipelagic (R² = 0.66), likely because it is the most abundant. The structuring environmental variables were different for the epi- ( Figures S8A, B ) and mesopelagic layers ( Figures S8C, D ): temperature (33%) and oxygen (19%) for the former and temperature (29%), bathymetry (19%) and chlorophyll a (15%) for the latter. The highest copepod biomass in the top 200m was found in high latitudes ( Figure 5A ), where water temperature is low and oxygen concentrations are relatively high. In the mesopelagic layer ( Figure 5B ), high copepod biomass was associated with shallow coastal and cold water masses. The patterns of distribution predicted by the global models were similar in both layers ( Figures 5A, B ), with the highest predicted biomass values in the Baffin Bay, Labrador Sea and Greenland Sea as well as at the Southern Ocean polar front region. The lowest predicted biomass was predicted at oceanic gyres and in the Arctic, north of 80°N. For both layers, the highest values of the coefficient of variation ( Figure 5C ) were found north of Canada and Greenland, as well as south of 60°S, especially for the epipelagic layer. These high values depict disagreement among the 100 models fitted to the data resamples and therefore inform on the uncertainty of the model in these zones. Caution is therefore advised regarding the interpretation of the very low values of biomass predicted in those regions. In the northern hemisphere, except for the Arctic ocean, the values of the coefficient of variation were rather low at locations where either low or high biomass values were predicted. In the southern hemisphere, model predictions varied relatively strongly at the level of the Antarctic polar front ( Figures 5C, D ).

Eumalacostraca contains mostly vignettes of euphausiids, amphipods and decapods. They were predicted globally with an R² of 0.31 for the epi- and 0.1 for the mesopelagic layer, both with significant p-values (p-value < 0.05; Figure 3 ). In the epipelagic, high biomass of these organisms was associated with high concentrations of phosphate (22%) and low concentrations of silicate (17%) ( Figures S9A, B ). In the mesopelagic layer, the distribution of this group was associated with low concentrations of silicate (16%), bathymetry (15%) and high chlorophyll α (15%) ( Figures S9C, D ). In terms of spatial distribution, high biomass is predicted in eastern boundary currents, especially in the Peruvian and Californian upwelling systems. Low biomass is predicted in high latitudes and in the oceanic gyres, especially in the North Atlantic. Similar patterns were predicted in the mesopelagic layer, but with lower biomass values. The model uncertainties are highest in the zones of low biomass (high latitudes and oceanic gyres).

Solitary Collodaria were predicted globally with an R² of 0.1 for the epi- and 0.09 for the mesopelagic layer, both with significant p-values (p-value < 0.05; Figure 3 ). In the epipelagic, the distribution of solitary Collodaria were mainly associated with low salinity (21%, between 35 and 37) and bathymetry (14%) ( Figures S10A, B ). In the mesopelagic, high abundances of this group were associated with distance to shore (18%) and high chlorophyll a (17%)( Figures S10C, D ). In this layer, 65% of the biomass was predicted at less than 1,000 km from the coast. Solitary collodaria were mainly located between 50°N and 50°S, in a rather diffuse manner ( Figure 7 ) with maximum biomass predicted at the equator. In the intertropical region, the highest biomass was found in the epipelagic zones of productive areas such as the upwelling regions off the western coast of Africa (Cape Verde and Angola) and of the eastern boundary of the Pacific Ocean (Peru and California). The model also predicted high biomass in the Mediterranean Sea. The importance of the environmental variable “distance to coast” in the learning process created unusual patterns in the prediction map such as a hexagonal shape in the Pacific Ocean. North of 50°N and south of 50°S, environments that are typically characterized by water masses with low salinity (1^st most structuring variable in the epipelagic) and high nitrate (4^th variable), the predicted biomass was rather low especially in the epipelagic layer.

For this group, the worldwide epipelagic model was statistically significant (p-value < 0.05; Figure 8 ) with an R² of 0.27, but the mesopelagic model was not (p-value > 0.05; Figure 3 ). Therefore, only the 0-200 m layer is displayed ( Figure 8 ). In this layer, Phaeodaria was one of the best predicted groups ( Figure 3 ) especially in the upper 200m. The predicted epipelagic distribution of Phaeodaria is associated with low values of salinity (38%) followed by bathymetry (11%), surface chlorophyll a (10%), oxygen and temperature (8% each) ( Figures S11A, B ). This is visualised on the map of global prediction ( Figure 8A ) on which high biomass was mainly predicted in the Californian upwelling (characterized by low salinity, cold and coastal waters), with lower biomass north of the upwelling up to the Gulf of Alaska. High biomass values were also predicted in the Bay of Bengal and Adaman Sea. The coefficient of variation in zones of high biomass is very low, providing strong confidence in this pattern. The lowest predicted biomass for this group are found in oceanic gyres and high latitudes of the northern hemisphere.

The group Acantharea was predicted with low total biomass ( Figure 4 ). This group was well predicted in the world model fitted with the epi- (R² = 0.26) and mesopelagic (R² = 0.63) layers ( Figure 9 ). In the epipelagic layer, nitrate (18%), salinity (15%) and phosphate (12%) were the main driving variables ( Figures S12A, B ). In the mesopelagic layer, the link between biomass and environment ( Figure 9B ) was defined by the influence of several variables: silicate (19%), phosphate (12%) followed by chlorophyll a (12%) ( Figures S12C, D ). The highest epipelagic biomass ( Figure 9A ) was predicted in the intertropical range, in productive areas such as the upwellings off the West coast of Africa (Cape Verde, Angola) and America (Peru and California). These high biomass patches are associated with a salinity around 35 as the 2^nd most structuring variable, as well as with high nitrate and phosphate concentrations (respectively 1^st and 3^rd). Intermediate biomass values were predicted mostly between 50°N and 50°S in a diffuse way, except in the oceanic gyres where the predicted biomass was lowest. The largest uncertainty was present in the Southern and Artic Oceans, Bering Sea and Gulf of Alaska where low biomass values were predicted ( Figure 9C ). In the mesopelagic layer, biomass was predicted to be 16.7-times lower overall ( Figure 9B ), with highest values found in the Gulf of Alaska and the Bering Sea. Intermediate biomass values were predicted for the upwelling regions and the Southern Ocean. In this layer, the high biomass estimates correspond with low coefficient of variation values ( Figure 9D ).

Copepoda biomass was predicted to be highest in high latitudes in both epi- and mesopelagic layers of the global models. The lowest values were predicted at the gyres and an increase of biomass was observed centered at the equator. In the global models, temperature always appeared within the top three environmental factors explaining the distribution of copepods (except for 0-100 m model where it appeared 4^th), which is in agreement with previous work suggesting that surface temperature and thermal tolerance of marine ectotherms, including copepods, are important constraints for their distribution and abundance (Beaugrand et al., 2009; Sunday et al., 2012). We also predict significant Copepoda biomass centered at 50°S in the Southern Ocean, at the location of the strongest horizontal gradient of temperature within the epipelagic layer. This geographic pattern is in agreement with earlier observations of high Copepoda occurrence along the Polar front (Pinkerton et al., 2020). Hence, despite a low number of UVP5 profiles in this latitudinal band, the model is able to retrieve this fundamental pattern. Higher values of the coefficient of variation ( Figure 5C ) are found in the Arctic Ocean, as well as south of 60°S. More data from these regions could help to further reduce the uncertainty of our models.

The distribution of the predicted Eumalacostraca biomass showed high values in coastal areas mainly on the eastern boundary currents of the Atlantic and Pacific Oceans and low values at high latitudes and at the locations of the oceanic gyres. Due to the low image resolution, a finer taxonomic resolution than Eumalacostraca (mostly euphausiids, decapods and amphipods) is not possible for UVP5 vignettes. Euphausiids are well known for their ability to escape standard oceanographic plankton nets (Brinton, 1967; Wiebe et al., 1982; Sameoto et al., 1993) and even low noise gliders (Guihen et al. 2022). This behavior might also be dependent on the species and stage development while the UVP5 mostly detects small Eumalacostraca (≤ 50 mm) for which taxonomic identification is not possible. Nevertheless, as Euphausiids are the second most abundant crustacean taxon after copepods (Castellanos et al., 2009), they may compose a large fraction of the biomass in this group. They are described as widely distributed in high numbers in the world ocean between 0-300 m with the exception of the eastern Canadian Arctic and the Arctic Ocean (Castellanos et al., 2009). This is consistent with our predictions of higher biomass in the epipelagic zone (0.058 PgC) compared to the mesopelagic (0.049 PgC), and low values predicted for the Arctic Ocean. The high Eumalacostraca biomass predicted in the North Atlantic also consistent with other observations that reported high abundances of krill in this region (Edwards et al., 2021). Euphausia superba and Euphausia mucronata have been respectively described as keystone species of the Antarctic and the Humboldt Current System (Antezana, 2010). The comparatively low values of biomass predicted in the Antarctic in the epipelagic layer ( Figure 6A ) might be too low, as Euphausia superba is known to show a patchy distribution (Siegel, 2005; Siegel, 2016). Since we only have very few samples from the Antarctic Ocean, we probably under-sampled this region and specifically krill. The high coefficient of variation in this region seems to reflect this problem. Overall, our observations and models likely underestimate the abundances of Euphausiids and of Eumalacostraca, due to their escape behaviors, the comparatively small sampling volume of the system and the low sample size in the Southern Ocean.

Global models in epi- and mesopelagic layers predicted a widespread distribution of solitary Collodarians between 50°N and 50°S, from oligotrophic to eutrophic zones. Their distribution can be explained by the selective advantage of their mixotrophy, since all collodarian species live in symbiosis with photosynthetic microalgae (Suzuki and Not, 2015; Biard et al., 2016). Consistently with the models’ prediction of solitary Collodaria as the third most important group in terms of global biomass in 0-200 m, it has been shown by Biard et al. (2016) that Collodaria contribute most to the biomass of the Rhizaria between 0-100 m.

The distribution of Phaeodaria shows a latitudinal pattern with three peaks in biomass, at 50°N (with high biomass values at the level of the subarctic gyres), at 5°N and at 60°S. These three peaks were not observed by Biard et al. (2016). The highest values being predicted in the subarctic gyre are consistent with Steinberg et al. (2008) who estimated their mean biomass there as 5.5% (range 2.7–13%) of the metazoan biomass sampled using a MOCNESS (Wiebe et al., 1985). The distribution of this group in the epipelagic (high biomass in coastal regions especially around the Californian upwelling and low biomass in the gyres conditions) could be related to food availability which might not be abundant enough in the open ocean. In the models’ output, this group only accounted for to ~ 1.2% of the global biomass in the epipelagic. This is consistent with previous work describing these organisms as being distributed in water below 150-200 m (Stemmann et al., 2008; Nakamura and Suzuki, 2015; Boltovskoy et al., 2017; Biard and Ohman, 2020). The high (R² = 0.50) and low latitude (R² = 0.39) models for the mesopelagic layer reveal similar patterns as the ones shown for the epipelagic layer in Figure 8 . This pattern of high biomass predicted in the North Pacific can be put in perspective with a previous study (Ikenoue et al., 2019) which highlighted Phaeodaria in the Western North Pacific as one of the major carriers of carbon in the twilight zone (200-1000 m (Buesseler and Boyd, 2009)), with an organic carbon standing stock reaching its highest value at depths between 200-500 m. A maximum in abundance of Phaeodaria was observed in the lower epipelagic or mesopelagic zone in the Sea of Japan by Nakamura et al. (2013) as well as in the Antarctic beneath the sea ice with similar abundances as the North Atlantic and Pacific (Morley and Stepien, 1984). In the regional mesopelagic predictions, the mean biomass in the Sea of Japan is not particularly high, but it reached higher values in the Southern Ocean.

Here, we present results on large Acantharea only, but it should be kept in mind that most species are smaller than 600 μm (Biard et al., 2016). Most Acantharea species are associated with symbiotic algae (Michaels, 1991) which could explain the rapid observed biomass decline with depth. Indeed, the biomass predicted is 16.7-times lower in the mesopelagic (1.36 10^-5 PgC) compared to the epipelagic layer (2.27 10^-4 PgC). These mixotrophs are present throughout the world oceans (Suzuki and Not, 2015) and commonly distributed in intertropical latitudes (Bottazzi and Andreoli, 1982) mostly in the surface with an abundance rapidly declining below 20-50 m depth (Michaels, 1988). The model confirmed this biomass diminution in the epi- and mesopelagic layers ( Figure 9 ). We also observed latitudinal patterns with the highest biomass in intertropical areas consistent with these previous studies. The highest biomass of Acantharea predicted by the mesopelagic global model in the Gulf of Alaska coincides with a large number of organisms imaged by the UVP5. This is surprising knowing the above described distribution patterns. More observations from this region are required to clarify whether this was a temporally limited occurrence or whether it represents a region of permanent abundance maxima. The predicted biomass in Antarctic waters in this depth layer is also surprising. Acantharea are marine planktonic unicellular eukaryotes in the Rhizaria group and produce a mineral skeleton made of strontium sulfate (Michaels, 1991; Decelle and Not, 2015). The surprisingly high abundance at high latitudes might be important for studies done on the strontium biogeochemical cycle (Bernstein et al., 1987; Decelle et al., 2013).