Global estimate of mesopelagic mesozooplankton biomass

Mesozooplankton are a pivotal component of marine trophic webs, ranging in length from 0.2 to 20 mm¹. They consist primarily of crustacean plankton (copepods, amphipods, euphausiids), meroplanktonic larvae, rhizaria, and smaller gelatinous zooplankton^2,3. Mesozooplankton play a crucial role in pelagic ecosystems by serving as consumers of lower trophic levels⁴ (e.g., ciliates, heterotrophic dinoflagellates, microzooplankton) and prey for higher trophic levels (e.g., fish and macrozooplankton), influencing energy flow, carbon sequestration, and nutrient cycling within pelagic ecosystems⁵.

The mesopelagic zone, also known as the twilight zone, extends from 200 to 1000 m below the surface and is characterized by low light levels, supporting a diverse community that includes mesozooplankton. Despite the importance of mesopelagic zooplankton in oceanic food webs and their roles in biogeochemical cycling, global biomass estimates are scarce and mostly spatially limited^6,7,8. Some ocean basin biomass estimates exist, but they are based on limited studies^9,10. Vereshchaka et al.⁹ discovered a strong relationship between surface chlorophyll concentrations and zooplankton biomass, aligning with theoretical expectations of connectivity between surface and deep-sea layers. The first available global ocean maps of zooplankton biomass are outdated and based on limited data, hand-drawn, and cover only the top 100 m of the epipelagic layer^11,12. Strömberg et al.¹³ addressed this issue by developing a model that related the flow of energy from primary production to zooplankton biomass, providing a map of global net zooplankton distribution in the epipelagic zone (0–200 m) and reported a mean global biomass of 5.52 ± 8.94 mgC m^− 3.

Recent estimates of epipelagic mesozooplankton standing stocks from global databases range from 0.19 to 0.49 PgC^3,14,15. These estimates are difficult to compare due to redundancies, differences in units, and variable sampling techniques. Furthermore, the range of global stock estimations of epipelagic mesozooplankton biomass varies significantly both spatially and with depth, complicating the incorporation of a global range into models focused on specific oceanic basins. Studies have also assessed mesozooplankton biomass in the mesopelagic zone^14,16,17. The global integrated biomass in the epipelagic and upper mesopelagic zone (0–500 m) varied between 0.14 and 0.40 PgC. Hernández-León et al.¹⁷ is the only study examining mesopelagic zone (200–1000 m) distribution, finding a strong link between dark ocean biomass and average epipelagic primary productivity. They estimated total oceanic mesozooplankton biomass at 1.4 Pg C, with ~ 0.5 and ~ 0.7 PgC located in epipelagic and mesopelagic layers, respectively¹⁷. This assumes uniform mesopelagic zone depth globally. However, Reygondeau et al.¹⁸ revealed that the vertical division of the zone is not constant over the global ocean but varies between ocean basins and latitudes.

Zooplankton biomass is influenced by temperature^16,19, primary productivity^9,17, and predatory pressure²⁰. Net Primary Production (NPP) is an important factor, but not all NPP is exported to the mesopelagic zone, with fractions ranging from < 10 to > 50% in central gyres and polar regions²¹. Another assessment method considers the standing stock of particulate organic matter (POM), which includes decomposing and sinking organic matter²². Particulate Organic Carbon (POC), a key POM component, represents the carbon fraction within the suspended organic material in the water column. POC is important because it can serve as a source of nutrition for many organisms and may provide a more direct and inclusive proxy for the zooplankton biomass.

This study aimed to determine the global mesopelagic mesozooplankton biomass (MMB) using NPP and POC standing stocks as predictors and compare these findings with previously published data. It also examined correlations between epipelagic and mesopelagic zooplankton biomass and explored the relationship between mesozooplankton diversity and biomass in the mesopelagic zone.

Models

We fitted an ordinary least squares (OLS) regression model to predict MMB with NPP data (formula: \(\:{\text{log}}_{10}Bio\sim{\text{log}}_{10}NPP\)) using the data from variable mesopelagic depth (Fig. 1A). The model explained a statistically significant and moderate proportion of the variance (R² = 0.24, F(1, 42) = 13.03, p< 0.001). The effect of logged NPP was statistically significant and positive (β = 2.08, 95% CI [0.92, 3.24], t(42) =3.61, p < 0.001). The same model applied to 200–1000 m depth data (Fig. 1B) explained a statistically significant and moderate proportion of the variance (R² = 0.18, F(1, 40) = 9.06, p = 0.005). The effect of logged NPP was statistically significant and positive (β = 1.62, 95% CI [0.53, 2.70], t(40) = 3.01, p = 0.005).

Relationships between NPP (1st column) or POC (2nd column) and mesozooplankton biomass for Variable Depth Models (1st row) and Classical Depth Models (2nd row). Solid blue lines represent regression lines, and light blue areas are 95% confidence interval. Note that both axes are on a logarithmic scale.

Full size image

We fitted an OLS model to predict MMB with POC (formula: \(\:{\text{log}}_{10}Bio\sim{\text{log}}_{10}POC\)) using the data from variable mesopelagic depth (Fig. 1C). The model explained a statistically significant and substantial proportion of the variance (R² = 0.30, F(1, 42) = 18.21, p < 0.001). The effect of logged POC was statistically significant and positive (β = 2.71, 95% CI [1.43, 3.99], t(42) = 4.27, p < 0.001). A similar formula for 200–1000 m depth data (Fig. 1D), explained a statistically significant and moderate proportion of the variance (R² = 0.24, F(1, 40) = 12.42, p = 0.001), slightly lower than the variable depth model. The effect of logged POC was statistically significant and positive (β = 2.18, 95% CI [0.93, 3.43], t(40) =3.52, p = 0.001). Models’ coefficients and diagnostics are presented in Tables S1-2.

Projection into the mesopelagic zone

Overall, all approaches showed similar MMB distributional patterns (Fig. 2) with higher biomass values observed in the northern hemisphere, along the west coasts of continents, and in the equatorial and 50°S bands compared to other regions. Low biomass was found in ocean gyre centers in both hemispheres. NPP-derived biomass decreased at higher latitudes (Fig. 2A,E) while POC-derived values increased toward the higher latitudes in the southern hemisphere and stayed constant in the northern hemisphere (Fig. 2C,G). The latitudinal gradient for NPP-derived biomass showed a greater range with a sharp decrease towards higher latitudes (Fig. 2B,F). Both POC and NPP-based values showed low biomass at 30°S and 30°N with a slight increase in the equatorial area. Biomass increased from 30° to higher latitudes.

NPP- (1st column) and POC- (2nd column) based estimates of mesopelagic mesozooplankton biomass in the global ocean using variable depth models (1st row) and classical depth models (2nd row). Latitudinal biomass gradients for each model are shown in plots (B), (D), (F) and (H). Solid black line in latitudinal graphs are the mean biomass value per latitude. All maps are in Eckert IV global equal area projection.

Full size image

Global MMB in the mesopelagic zone was estimated at 0.29 ± 0.06 PgC using NPP-derived estimates and three-fold higher values of 0.91 ± 0.08 PgC using the POC-based estimates (Table 1). Using the classical definition of the mesopelagic zone (200-1,000 m), NPP-based values were lower than POC-based values (0.20 ± 0.06 PgC and 0.45 ± 0.08 PgC, respectively). Biomass estimates for 40°S-40°N were 68–70% of global estimates for NPP values and 55% for POC-derived estimates.

Biomass and diversity

We used biomass estimates derived from the variable depth POC model to explore the relationship between biomass and diversity. Average biomass values per mesopelagic province are shown in Fig. S1.

Log-transformed biomass and species richness followed a logarithmic trend: rapid species richness increase at lower biomass values, followed by slower changes as biomass increased (Fig. 3). Lower biomass and species richness were associated with areas of low POC.

Global relationship between mesopelagic zooplankton biomass (mgC m^-2) and potential species richness colored by the POC value for each cell. Note the log₁₀ transformation on the \(\:x\) and \(\:y\) axes and colour scale. The black line is fitted gam line (formula = y~s(x, bs = “cs”)).

Full size image

High biomass and species richness were observed in the Arctic, North Pacific, coastal and equatorial areas, and Southern Ocean (Fig. S3A). Low biomass and low species richness were found in the North, South Pacific, and Indian Ocean basins. However, the center of the North Atlantic basin had low biomass but high species richness. Low biomass and high species richness were also recorded in the Mediterranean Sea and western boundaries of the subtropical gyres. Areas of high biomass, but low species richness were in subtropical regions (both sides from the equatorial band), around the gyre edges, and some areas in the North Pacific and the 50°S band.

The latitudinal pattern of median MMB resembled that of species richness (Fig. S3B) but with a larger magnitude of change. The lowest species richness and biomass values were found around 23.5°N and 23.5°S, gradually increasing towards the equator. In the northern hemisphere, biomass increased through the temperate region and stabilized at high values in the polar region. Species richness showed a similar increase in the northern temperate region, but the increase was not gradual, with some rises and declines around 45–50° N. In the southern hemisphere, both biomass and species richness increased throughout the temperate region, with a small peak around 30°S. Biomass continued increasing towards the southern pole, while species richness showed a moderate increase.

A skewed hump-back relationship was observed between log-transformed MMB and the total-range size rarity (Fig. S3A). Biomass remained low until total rarity reached ~ 0.3 (three in ten species are rare) and then declined. Species richness was highly and positively correlated to total rarity. The average rarity for each cell was calculated to remove the effect of species richness on the total rarity (Fig. S3B). Biomass was positively associated with average rarity, increasing as rarity increased. Biomass variability increased with higher average rarity scores. A peak in biomass was observed at average rarity scores of 0.3, associated with the highest species richness values. Biomass was generally lower at high average rarity scores. This pattern was evident in almost all ocean basins except the Arctic Ocean, where only high average rarity cells were found (Fig. S3B).

Link between biomass in epi- and mesopelagic zones

The global ratio of epipelagic to mesopelagic biomass varied greatly due to different estimates of zooplankton biomass in the epipelagic layer (Table 2). The smallest estimate, 0.14 PgC, was from a digitized Bogorov’s map, covering only the top 100 m layer. Other epipelagic estimates were higher (0.19–0.49 PgC) and were reported from the 0–200 m depth bin. Mesopelagic biomass estimates also varied (0.45–0.91 PgC), influenced by different mesopelagic depth classifications. Using the classical mesopelagic zone definition, two estimates were 0.66 PgC from Hernández-León et al.¹⁷ and 0.45 PgC from this study (POC-based). We also included a POC-based estimate using variable depth for better representation.

Using Hernández-León et al.¹⁷ estimates, ratios ranged between 0.32 and 0.74 (mean 0.53 ± 0.18, median 0.48). Excluding Bogorov et al.¹¹ (top 100 m only), classical depth POC-model ratios ranged from 0.47 to 1.01 (mean 0.76 ± 0.24, median 0.71). Variable depth estimates had the lowest ratios, 0.21–0.53 (mean 0.38 ± 0.13, median 0.35).

Epipelagic to mesopelagic ratios varied by ___location. Using Strömberg et al. epipelagic estimates (Fig. 4A), the mean ratio was 0.89 ± 0.53 (range 0-4.66, median 0.85). High ratios (~ 2.2) were observed between 20° and 30° in both hemispheres, especially in ocean basin centers (Fig. 4B). Polar and equatorial areas had lower ratios (< 1). The lowest ratios were in upwelling regions along Africa’s west coast and near-coastal areas of the North Atlantic and North Pacific (Fig. 4B). Based on Bogorov et al.¹¹ estimates (Fig. 4D), ratios showed no clear pattern (Fig. 4E) due to the binned nature of original values (Fig. S2, with an average ratio of 0.98 ± 1.16 (median 0.67), with ranges 4 times higher than calculated based on epipelagic biomass estimates from Strömberg et al. (2009).

(A) Epipelagic net zooplankton biomass (in mgC m^-3) from Stromberg et al. (2009). Note the log-transformed colour scale. (B) Ratio between epipelagic and mesopelagic biomass using Stromberg et al. (2009) and epipelagic zooplankton and mesopelagic biomass from the variable-depth POC model. Note the square root of the scale for the colour bar. (C) 2D map depicting the relationship between mesopelagic mesozooplankton biomass derived from POC estimates and Stromberg et al. (2009). Quantile breaks in the legend were created at the 25th, 50th, and 75th percentiles. (D) Epipelagic net zooplankton biomass (mg m-3 wet weight) from Bogorov et al. (1968). (E) Ratio between epipelagic and mesopelagic biomass using Bogorov et al. epipelagic zooplankton and mesopelagic biomass from the variable-depth POC model. Note the square root of the scale for the colour bar. (F) 2D map depicting the relationship between the mesopelagic mesozooplankton biomass derived from POC estimates Bogorov et al. estimated the total zooplankton biomass in the top 100 m of the water column. Quantile breaks in the legend were created at the 25th, 50th, and 75th percentiles for mesopelagic biomass and < 50, 51–200, 201–500, > 500 mg·m⁻³ wet weight intervals for Bogorov’s et al. biomass estimates. All maps are in Eckert IV global equal area projection.

Full size image

In addition to computing the ratio between epipelagic and mesopelagic biomasses, we also utilized the maps with a 2D legend to simultaneously visualize both variables, effectively highlighting their interdependence and spatial relationships (Fig. 4C, F). Both epi- and mesopelagic biomass were low in ocean gyre centers and high at high latitudes, equatorial areas, and west coasts of continents. Low surface biomass and high mesopelagic biomass were found around gyre regions and the 50°S and 50°N bands.

Mesopelagic mesozooplankton biomass

The study used OLS regression to examine the relationship between MMB and environmental conditions using POC or NPP data. Despite model simplicity, the models could explain 18-30% of the variance in the biomass estimates depending on the method used. Simple OLS regression performed better (higher R² values) than more complex models (boosted regression trees) built on a more complex set of environmental variables¹⁶. NPP-based models showed lower R² values (18-24%) than POC-based models (24-30%) but lower than those reported for a similar model by Hernández-León et al.¹⁷. The lower R² values of NPP-based models could be due to climatological data not accurately representing 2011 conditions when biomass data were collected. The Malaspina data was gathered during a single year (2011), whereas climatological data represent average conditions over a decade, meaning year-specific variations may not align with the long-term mean.

When modeling presence/absence of mesopelagic mesozooplankton, NPP was one of the least important factors²³. However, this work demonstrated that NPP is a critical factor for modelling mesozooplankton biomass and can explain up to 24% of MMB variation. POC was a better predictor, explaining up to 30% of the variance (Table 1). Both variable and classical depth models produced similar MMB distribution patterns (Fig. 2A, E, C, G), with enhanced values in the northern hemisphere, along west coasts, in the equatorial region, and at the 50°S band (Antarctic Polar Front). Low biomass values were found in the center of ocean gyres. These patterns are consistent with study of Drago et al.¹⁶.

A latitudinal biomass distribution between the 30°N and 30°S showed lowest values at 30°N and 30°S, gradually increasing towards the equator (Fig. 2). This pattern matches the latitudinal distribution of Copepoda biomass reported by Drago et al.¹⁶. Crustaceans, comprising ~ 74% of total mesopelagic biomass¹⁶, show similar spatial patterns as reported in this study. NPP-derived biomass corresponded well with copepod biomass distribution estimated using boosted regression trees, while POC-derived biomass closely matched spatial patterns reported by in situ UVP5 samples¹⁶.

The main difference between NPP- and POC-derived projections was at higher latitudes, where NPP-derived values decreased, and POC-derived values increased or remained stable. NPP-based models showed a biomass peak around 40°-50° latitudes, while Drago et al.¹⁶ recorded it at 60° latitudes (BRT model). The discrepancy could be due to lower NPP estimates in polar regions (Fig. S5A), while POC values remained relatively high (Fig. S5B). The choice between NPP and POC and the discrepancy between estimates could be related to the omnivory proportion in mesozooplankton, particularly in the midwater realm^24,25. Omnivorous and carnivorous mesozooplankton outweighed herbivorous ones in the mesopelagic zone²⁶, where sinking marine snow and detritus are the main organic matter delivery modes²⁷. Thus, the POC-based biomass estimates would be, in our view, more realistic. This discrepancy may also stem from the greater importance of coprophagy and detritus consumption over herbivory in the mesopelagic zone, making POC a more accurate predictor of mesozooplankton biomass on a global scale.

Spatial coherence of biomass hotspots showed NPP-based variable-depth and classical depth model estimates were 68% and 56% lower than POC-based counterparts, respectively. This highlights the importance of considering the dynamic depth of the mesopelagic zone in biomass calculations because, in many instances, the real depth range of the mesopelagic zone was larger than inferred using the classically accepted boundaries that underestimated the total volume of the mesopelagic by 40% (18).

Vereshchaka et al.⁹ estimated the standing stock of the zooplankton community in the Atlantic Ocean at 70 Mt (wet weight), with the mesopelagic zone contributing between 13 and 16% or 9.1–11.2 Mt. They found a strong relationship between surface chlorophyll concentrations and zooplankton biomass⁹. Several recent global mesozooplankton standing stock assessments exist (Table 3). Moriarty & O’Brien³ estimated global epipelagic biomass at 0.19 PgC, with highest concentrations in the northern hemisphere and a slight decrease from polar to temperate regions. Buitenhuis et al.¹⁴ provided similar estimates, while Hatton et al.¹⁵ reported 41 Gt (wet weight) or roughly 0.49 PgC. None provided uncertainty information, making comparisons difficult.

Studies also investigated mesozooplankton biomass in the mesopelagic zone (Table 3). Drago et al.¹⁶ estimated the global biomass distribution of 19 zooplankton taxa at roughly 3500 stations using in situ imaging and machine learning. They found the highest biomass in polar regions and the equator, and lowest in oceanic gyres. Global integrated biomass in the epipelagic and upper mesopelagic zone (0–500 m) was 0.403 PgC, with Copepoda (35.7%) and Eumalacostraca (26.6%) being most abundant, and the upper mesopelagic zone (200–500 m) accounting for 0.173 PgC¹⁶. Buitenhuis et al.¹⁴ reported biomass values of 0.33–0.59 µgC L^− 1, approximating 0.14–0.25 PgC (Table 3). Hernández-León et al.¹⁷ estimated global oceanic mesozooplankton biomass at 1.4 Pg C, with 47% (0.66 PgC) in the mesopelagic and 34% (0.48 PgC) in the epipelagic layers (Table 3).

Our study found global mesopelagic mesozooplankton biomass varied between 0.20 and 0.91 PgC, depending on the method (Table 1). The classical depth NPP model estimated 0.20 PgC, similar to Drago et al.¹⁶ at 0.173 PgC. Drago et al. assessment was confined to the upper mesopelagic layer implying most mesozooplankton biomass is in the upper mesopelagic layer^28,29. All estimates, except the variable depth POC model, showed lower estimates than reported by Hernández-León et al.¹⁷. For the 200–1000 m depth range, they reported 0.21–0.46 PgC higher than POC- and NPP-derived estimates (Tables 1, 3). Lower biomass estimates can be explained by the exclusion in the analysis of highly productive coastal areas. Higher estimates in our study using the variable depth POC model highlights the larger volume of the mesopelagic zone usually underestimated by the classical division. Variability in estimated global MMB by Hernández-León et al.¹⁷. precluded statistical comparison with our estimates.

Latitudinally, MMB within the 40°S–40°N belt accounted for 68–70% of the global biomass estimates using NPP values, slightly higher than 55% from POC-derived estimates (Table 1). NPP-based models likely underestimate mesopelagic biomass at higher latitudes. POC-based models showed higher average biomass at higher latitudes (Fig. 2F, H), consistent with Fernandez de Puelles et al.²⁹, who found more biomass in high latitudes than in the tropics. Polar primary productivity seasonality produces intense periods several months a year³⁰, with the pelagic community relying on microbial and secondary production most of the year³¹. In such cases, POC-based models appear more realistic for determining MMB, particularly at higher latitudes. Additionally, overwintering (diapausing) zooplankton significantly contribute to seasonal mesopelagic biomass^32,33,34.

Biomass and diversity

This study showed a positive relationship between zooplankton biomass and species richness in marine environments (Fig. 3). Log-transformed species richness increased sharply with log-transformed biomass until reaching 1 gC m^− 2, then plateaued at ~ 300 species per cell. POC-based estimates showed high biomass regions tended to have high species richness (Fig. S3A), especially in the Arctic, North Pacific, coastal, equatorial areas, and Southern Ocean. Low biomass regions were often in North/South Pacific and Indian Ocean basins. Low biomass and high species richness values were found in the North Atlantic Basin center, Mediterranean Sea, and western boundaries of subtropical gyres. High biomass but low species richness were in subtropical regions around gyre edges and some North Pacific and southern hemisphere areas. While high biomass can indicate high species diversity, other factors like temperature, nutrient availability, and ocean circulation patterns also influence MM species distribution^17,23,35. Future studies should explore these factors to better understand the biomass-diversity relationship in mesopelagic zone. Species richness estimates in this study were based on species distribution models and may not reflect actual cell richness due to model accuracy dependence on species occurrence data³⁶. Rare species or areas with low sampling efforts may underestimate species richness values.

Latitudinally, median MMB mirrored species richness patterns (Fig. S3B). Both increased towards the equator, with lowest values around 23.5°N and 23.5°S. In the northern hemisphere, biomass increased throughout the temperate region, stabilizing at high values in the Arctic. Species richness showed a similar increase in the northern temperate region, with fluctuations around 45–50°N. In the southern hemisphere, biomass and species richness increased throughout the temperate region, peaking around 30°S. Biomass continued increasing towards the southern pole, while species richness showed moderate polar increases. These observations align with Drago et al.¹⁶, who observed high values north of 55°N and south of 55°S, and enhanced values around the equator (15°N-15°S). The main difference was observed in the Southern Ocean where Drago et al.¹⁶ reported high biomass between the Subantarctic Front and the southern Antarctic Circumpolar Current limit, with lower values north and south of this band. This discrepancy is explained by the choice of environmental parameters, particularly the use of different prey availability proxies: NPP (Fig. 2A and C) versus POC-based calculations (Fig. 2E and G).

High mean values for both biomass and species richness were found at higher latitudes compared to the equatorial region (Fig. S3B). The cold, dynamic, seasonal environment, and long periods of darkness in polar regions result in slower speciation, favoring fewer but highly specialized species. The enhanced average range-size rarity supports this suggestion (Fig. S4). Despite polar systems’ seasonal productivity, POC values reflect longer-term organic matter accumulation, less influenced by seasonal NPP fluctuations. Polar regions export about twice as much of their production compared to tropical regions³⁷, supporting higher biomass at high latitudes.

High biomass levels in polar regions are model projections, as biomass values were collected at lower latitudes. Environmental variables used to build the models were comparable globally, assuming POC or NPP relationships are consistent across latitudes. However, this may not be true as primary productivity modes vary drastically. We discuss only mean climatology values of POC and NPP, not reflecting annual or depth variability. Hence, obtained biomass data are snapshots of mean climatological conditions and may differ seasonally or annually.

In some cases, the high biomass of a dominant species can reduce diversity by outcompeting others. Similarly, low biomass does not always mean low diversity, as some species thrive in low-biomass environments. For instance marine phytoplankton shows a unimodal productivity-diversity relationship, with diversity peaking at intermediate productivity levels³⁸. A negative correlation between biomass and species diversity of marine zooplankton is found in regional studies, e.g., in the Indian Ocean³⁹.

The hump-back relationship between log-transformed total range-size rarity and biomass indicates highest biomass at intermediate biodiversity levels (Fig. S4A). A skewed humpback relationship showed biomass peaking at ~ 3/10 of total range-size rarity (three rare species out of 10). Species richness and total range-size rarity were generally low at low biomass levels due to limited resources⁴⁰. Low biomass values were associated with ocean gyres (Fig. 2A), where water tends to be warm and relatively stable with little mixing between the surface and deeper waters. This limits phytoplankton growth and enhances recycling, resulting in overall low food availability and keeping the MMB low. High biomass values were found in areas with few rare species and intermediate species richness (Fig. S4A) and potentially underlie the way different species interact with each other and with their environment. There is a balance (sensu tradeoff) between competition and resource availability at the intermediate species richness levels suggesting that organisms can exploit available resources most efficiently^41,42,43. Redundancy at the intermediate richness levels still sufficient to buffer environmental fluctuations leading to a stable ecosystem⁴⁴. High species richness intensifies competition, limiting growth and reproduction, decreasing overall biomass^45,46.

A pattern was observed between the mesozooplankton biomass, average rarity, and species richness (Fig. S4B). The average rarity of 3/10 seems to be an important threshold associated with a peak in biomass and species richness. The biomass decreased rapidly as the average rarity decreased from 3/10. This pattern was observed on both the global and ocean-basin scales (Fig. S4B). Peak species richness is established in areas of low average rarity (one rare species out of ten), found primarily along continental coastlines and 50–60°N in the North Pacific and North Atlantic.

In summary, we conclude that MMB was linked with POC values. In areas where POC is abundant, the mesozooplankton biomass tends to be high, whereas in areas where POC is scarce, the biomass levels are low. Regions with low biomass are often located in the center of ocean gyres and are characterized by low species rarity and/or low species richness. In contrast, areas with high biomass can only be found in locations where both the species richness and rarity are at a certain level. Specifically, high biomass is associated with high species richness and is equal to or lower than three rare species in every ten species. Therefore, to achieve a high mesozooplankton biomass, both high species richness and mid-level rarity are necessary conditions.

Link between biomass in epi- and mesopelagic zones

Using the classical depth POC model, we found that the mean epi-/mesopelagic biomass ratio was 0.76, suggesting that 76% of the epipelagic biomass was located within the mesopelagic zone. This was similar to the ratio (0.72) reported by Hernández-León et al.¹⁷. However, when the variable depth of the mesopelagic zone considered, the mean ratio decreased to 0.38, indicating that only 38% of the epipelagic biomass was found in the mesopelagic zone. It appears that the zooplankton biomass in the mesopelagic ocean may be much larger than previously assumed. It is important to note that these ratios are highly dependent on the epipelagic and mesopelagic biomass estimates. We found that even the well-studied epipelagic zones may yield highly variable biomass estimates. In the mesopelagic zone, only one study attempted to estimate the biomass, and even that was based on the classical definition of the mesopelagic zone. This variability in the estimates can also be explained by the lack of global coverage of zooplankton biomass in the upper layers of the ocean. We also acknowledge that the ratio can additionally be influenced by the conversion between units. For example, we used a single factor of 0.12 to convert wet weight to carbon mass (Table 3). However, the mesopelagic community composition is not spatially and temporally uniform and can be composed of variable taxa in different regions, and taxa-specific carbon content can vary between 0.003 and 0.15⁴⁷.

The link between the surface and mesopelagic biomass cannot be easily determined, as very few studies have provided a continuous global estimate of zooplankton biomass^11,13. Using Strömberg et al.¹³ data, we found the global epipelagic/mesopelagic biomass ratio was primarily low, with over 60% of cells having ratios < 1, the majority below 0.6. The mean cell ratio was 0.85, indicating higher mesopelagic biomass. The remaining cells showed epipelagic biomass exceeded mesopelagic biomass. These patterns are closely related to the mesopelagic-bathypelagic boundary depth. Areas with high epipelagic/mesopelagic ratios coincided with the deepest mesopelagic layers^18,48, found in low-productivity ocean basin centers (Fig. S5). In contrast, low ratios were found north of 30° in both hemispheres and in some equatorial upwelling areas. These regions, characterized by strong seasonality and water mixing, support high seasonal zooplankton biomass in the epipelagic zone and high POC values that support larger mesopelagic biomass. Additionally, several temperate and subpolar large calanoids, after accumulating carbon in lipids during a few summer months, undergo seasonal vertical migration, reallocating biomass between epipelagic and mesopelagic realms and supporting the carnivorous midwater food web⁴⁹.

Diel vertical migration (DVM) between epipelagic and upper mesopelagic layers must be considered⁵⁰. A portion of estimated MMB represents diel migrants⁵¹, but the global biomass proportion is unknown. Zooplankton DVM increases global export flux by 14%^32,52,53. To validate this estimate, we also compiled the day/night epipelagic mesozooplankton biomass ratios from various regions of the ocean^{54,55,56,57,58,59,60}. Day/night epipelagic mesozooplankton biomass ratios from various ocean regions varied between 0.04 and 0.39, with a mean of 0.20 (median 0.17). Applying a 17% median ratio, migrant biomass was 0.15 PgC (variable depth model) or 0.08 PgC (classical depth).

Limitations of the study

Estimating mesopelagic mesozooplankton biomass (MMB) is inherently challenging due to methodological and ecological limitations. One of the primary constraints is linked to net-based sampling methods, which do not fully capture the diversity of mesopelagic organisms. Certain groups, such as gelatinous zooplankton, are particularly underrepresented in net catches despite their ecological significance, affecting biomass estimates⁶¹. Additionally, the vertical distribution of mesopelagic organisms is not static; depth boundaries fluctuate seasonally, adding uncertainty in biomass assessments⁶². Zooplankton abundance declines with increasing depth across ocean basins, impacting the accuracy of net-based biomass extrapolations⁶³.

This study derived MMB estimates from biomass data collected during the Malaspina Expedition, which globally primarily surveyed tropical and subtropical regions^17,64. While this dataset provides valuable large-scale insights, its limited spatial coverage may lead to an underestimation of biomass at higher latitudes, where mesopelagic communities poorly documented. Moreover, the expedition’s reliance on net sampling introduces inherent biases, particularly in detecting fragile and rare deep-dwelling taxa, leading to an incomplete assessment of mesopelagic biomass.

To standardize biomass estimates, we defined the mesopelagic zone as the 200–1000 m depth range. However, this static representation oversimplifies the mesopelagic ecosystem boundaries. The latter vary seasonally⁶², and biomass distribution within this zone is not uniform, generally declining with depth due to decreasing food availability and light penetration⁶³. Consequently, our approach may not fully capture the vertical heterogeneity of mesopelagic communities, contributing to uncertainty in our estimates. While textbook framework allows comparability across datasets, it underscores the need for improved methodologies that account for spatial and temporal variability in mesopelagic realm. Future efforts should integrate dynamic sampling techniques, such as acoustic and optical methods, to refine MMB assessments and improve our understanding of mesopelagic ecosystem structure and function^17,61.

Methods

All statistical analyses, data manipulations, and visualizations were performed using R Statistical Software⁶⁶.

Biomass data

We used data from the 2010–2011 Malaspina cruise, as recorded by Hernández-León et al.^17,64. Data were collected using a consistent sampling method across the entire sampling area to minimize variations in biomass estimates. The study covers a vast area ranging from 40°S to 30°N. The Malaspina cruise’s zooplankton biomass (in mgC m^− 2) was downloaded from PANGAEA⁶⁴. To allow for easier comparison of different areas on a global scale, biomass data were gridded on Map of Life (MOL; https://mol.org/) equal-area grid with a per cell area of 3091 km². In addition, we used the mesopelagic depth boundaries estimated by Reygondeau et al.¹⁸ and projected onto the same equal-area grid.

The original data contained three depth intervals (0-200, 200–1000 and 1000–2000 m) where the 200–1000 m interval was attributed to the mesopelagic realm. However, Reygondeau et al.¹⁸ demonstrated that mesopelagic depth varies globally, with the epi-mesopelagic boundary ranges between 17 and 150 m (mean 42 m) and meso-bathypelagic boundary ranging from 775 to 4197 m (mean 1358 m). Thus, the actual mesopelagic zone depth can overlap with all three sampling intervals from Hernández-León et al.¹⁷. To determine which of the sampling depth intervals were located within the mesopelagic zone, according to Reygondeau et al.¹⁸, we calculated the overlap (in meters) between the sampling depth (0-200, 200–1000 or 1000–2000 m) and the depth bin of the mesopelagic zone.

Environmental data

NPP was downloaded from 1080 by 2160 Monthly HDF files from MODIS R2018 Data using the Standard Vertically Generalized Production Model⁶⁷ for July 2002-December 2021 (http://orca.science.oregonstate.edu/npp_products.php). The seasonal climatology data for POC were downloaded from NASA Earth Data (https://oceancolor.gsfc.nasa.gov/l3/order/) Standard AQUA-MODIS POC 9 km (mapped) for the period 2002–2022. Then the climatological means for NPP and POC were computed for each cell. The data were then projected onto an equal-area grid, similar to the gridding of biomass data. The distribution of environmental variables is shown in Fig. S6.

Models

Two types of models were constructed. A linear model (estimated using OLS) was fitted to predict the biomass using NPP or POC data (formula: \(\text{log}\left(\text{B}\text{i}\text{o}\text{m}\text{a}\text{s}\text{s}\right)\sim\text{log}\left(\text{N}\text{P}\text{P}\right)\) or \(\:\text{log}\left(\text{B}\text{i}\text{o}\text{m}\text{a}\text{s}\text{s}\right)\sim\text{log}\left(\text{P}\text{O}\text{C}\right)\:\).

The mesopelagic zone depth can overlap with all three sampling intervals reported by Hernández-León et al.¹⁷, we only used biomass estimates from depth intervals, where 80% of the interval overlapped with the dynamic mesopelagic zone depth. The total number of intervals employed in this study was eight from 0 to 200 m, 42 from 200 to 1000 m, and nine from 1000 to 2000 m. Since each cell’s mesopelagic zone depth varied from the depth bin found in the data from Hernández-León et al., a corrected biomass estimate was employed to take this into account. Thus, the original areal biomass units (mg C m^− 2) were first converted to volumetric units (mg C m^− 3) using the corresponding sampling depth bins (0-200, 200–1000 or 1000–2000 m). Units were then converted back to areal estimates using the variable depth bin of the mesopelagic zone (mg C m^− 2 mesopelagic). Since these models use data from the variable depth of the mesopelagic zone based on latitude and longitude to model biomass, they will be referred to as ‘Variable Depth Models’ hereafter.

For comparison, we ran the same models but used the original biomass estimates for 200–1000 m depth from Hernández-León et al.¹⁷ to represent global MM biomass using the classic definition of the mesopelagic zone. These models will be referred to as ‘Classical Depth Models.’ Assumptions of each linear model were checked with diagnostic plots provided in Fig. S7. The 95% confidence intervals (CIs) and p-values were computed using a Wald t-distribution approximation⁶⁸.

Projection into mesopelagic zone

We utilized the annual climatological data for POC and NPP to make predictions of mesozooplankton biomass on a global scale using the linear models we created. To ensure that the predictions corresponded only to the mesopelagic zone, coordinates were filtered using the mesopelagic zone boundaries established by Reygondeau et al.¹⁸. This allowed us to estimate zooplankton biomass within the mesopelagic zone based on a new set of environmental conditions. Global biomass estimates were expressed in PgC, with the standard error reported.

The global estimate of the mesozooplankton biomass (\(\:\widehat{p})\:\)was calculated as the sum of the predicted biomass values in each cell:

where \(\:{x}_{i}\) is an array of \(\:n\:\)new values of the explanatory variable, \(\:{\beta\:}_{o}\) and \(\:{\beta\:}_{1}\) are coefficient estimates, and \(\:{y}_{i}\) is the predicted value for \(\:{x}_{i}\) (POC or NPP) while Σ is their variance-covariance matrix. \(\:{S}_{x}\) is a simplified notation for \(\:\sum\:_{i=1}^{n}{x}_{i}.\:\)The standard error of this prediction was calculated as a square root of \(\:\widehat{p}\) using the following formula:

Since Hernández-León et al.¹⁷ data only cover the area between 40°S and 40°N, we also calculated a global MM estimate for that particular band. However, biomass estimates for the entire global region can be considered adequate, as the range of environmental variables used to build the models is a representation of global ranges (Fig. S6).

Biomass and diversity

To explore the relationship between MMB and diversity, we matched the predicted biomass estimates from the variable depth POC model with the potential species richness computed using species distribution models^23,35. We used 2D maps of MMB and species richness to explore the global patterns of species diversity and biomass along the corresponding latitudinal gradients.

The total and average range-size rarities per cell were calculated using the species richness data. The total range-size rarity, also referred to as weighted endemism⁶⁹, was calculated as the proportion of the distribution of a species found in a cell summed across all species⁷⁰. We then explored the relationship between the estimated biomass and total range-size rarity. The total range-size rarity is correlated with species richness; therefore, the average range-size rarity was calculated as a more appropriate measure of rarity⁷¹. Average range-size rarity was determined by dividing the total range-size rarity by the number of species in each cell. We then examined the relationship between global and per-basin biomass and the average range-size rarity.

Link between biomass in epi- and mesopelagic zones

We used a global estimate of mesopelagic biomass reported from various sources (Table 3), three estimates of mesopelagic biomass (one reported by Hernández-León et al.¹⁷, and two estimates from a POC-based model from this study) to compute the global ratio between the epipelagic and mesopelagic biomass.

Two maps were utilized to study the relationship between MMB and total zooplankton biomass in the epipelagic layer. The first was the global map of zooplankton distribution by Bogorov et al.¹¹, which was digitized as the original data were unavailable (Fig. 4D). The second map (Fig. 4A) was constructed based on data provided by Störmberg et al.¹³. Data were available through the NOAA portal: https://www.st.nmfs.noaa.gov/plankton/biomass/.

Originally, biomass estimates in Bogorov et al.. map were binned into six intervals: <25, 25–50, 51–100, 101–200, 201–500 and > 500 mg m^− 3 (Fig. S2), so we converted biomass bins to numeric estimates by taking a mid-value of the bin as an estimate of zooplankton biomass in that cell (for bin > 500 mg m^− 3 we used a 500–1000 mg m^− 3 range and a mean value of 750 mg m^− 3). Then, the values were converted to mgC using a conversion factor of 0.12⁷² and were further multiplied by a depth bin of 100 m to convert them to areal estimates. For Strömberg et al.¹³, units were already given in carbon weight; therefore, only conversion to areal units was performed by multiplying values by a 200 depth bin. The values were then gridded on an equal-area grid. Both maps were matched with the mesozooplankton biomass estimated using the variable-depth POC model. The ratio between epipelagic and mesopelagic biomass was calculated for each cell. In addition, the global biomass of epipelagic was computed for each map.