Contribution of aquatic invertebrates to global nutrient supplies
Global capture and aquaculture statistics
To estimate the contribution of invertebrates to global animal capture fisheries and aquaculture production, we used 2019 reconstructed marine fisheries landings data from the Sea Around Us website25 and reported inland fisheries and aquaculture production from the Food and Agriculture Organization26. Reported inland fisheries, which are probably an underestimate of true inland fisheries production28, were analysed separately from marine capture fisheries throughout the Article, owing to the absence of globally standardized reconstruction datasets for freshwater systems. We used 2019 data in the main text because it is the latest year available with catch reconstructions. However, as a sensitivity analysis, invertebrate contributions were also estimated for other years (2014–2019). We found our invertebrate contribution results were consistent across years (Extended Data Fig. 3). The Sea Around Us uses officially reported landings from international and national fisheries statistics authorities as a baseline. Reconstructions are then performed by adding estimated unreported catches (such as landed illegal catches) using various literature sources. FishStatJ inland fisheries and aquaculture production26 data of territories and land areas were reported at the country level. Aquaculture and inland fisheries production species, reported using the Aquatic Sciences and Fisheries Information System (ASFIS) taxonomic reference system in FishStatJ, were converted to scientific species and species groups names (for example, rainbow trout to Oncorhynchus mykiss)60. Aquatic invertebrates in production data included mollusks (such as clams, mussels, oysters, scallops, cockles, snails, abalone, whelks, conchs, octopus, squid and cuttlefish), crustaceans (such as shrimp, prawns, crabs, lobsters and crayfish), sea cucumbers, sea urchins, sea worms, sponges and jellyfish (Supplementary Table 4). Marine capture fisheries and aquaculture production estimates shown throughout this Article (Fig. 1a and Extended Data Fig. 2) are global estimates derived from a total of 208,728 and 2,507 observations, encompassing 2,304 and 492 unique species (not only invertebrates) and 282 and 206 unique countries, respectively. Similarly, invertebrate-specific estimates for marine capture fisheries and aquaculture production shown throughout the Article (Fig. 1b and Extended Data Fig. 4) are global estimates derived from a total of 28,358 and 624 observations, encompassing 552 and 150 unique species and 252 and 124 unique countries, respectively.
Assigning nutrient composition data to aquatic foods
Nutrient concentration estimates per 100 g and edible proportions for each individual entry in global data (for example, fish and invertebrates) were assigned using raw muscle tissue samples within the AFCD12. We selected 30 nutrients that are important for public health13: 13 minerals (calcium, chromium, copper, iodine, iron, magnesium, manganese, molybdenum. phosphorous, potassium, selenium, sodium and zinc), 11 vitamins (A, C, D, E, B1 (thiamin), B2 (riboflavin), B3 (niacin), B5 (pantothenic acid), B6, B9 (folate) and B12 (cobalamin)), 5 versions of essential fatty acids (total monounsaturated fatty acids (MUFAs), total omega 3 fatty acids, total omega 6 fatty acids, DHA and EPA, and ALA), and protein. We included total protein and not specific essential amino acids because the dataset that we are using is currently being updated for those nutrients (for example, validating units and conversions). Moreover, we did not include other nutrients, such as cobalt and arachidonic acid, because they lacked established RNI reference values. Mean nutrient composition estimates and edible proportions were assigned hierarchically on the basis of the closest taxonomic resolution (Supplementary Fig. 2). In other words, for any given nutrient, if a species had species-specific nutrient composition observations in AFCD, we assigned that value. However, if no observed nutrient concentrations were available at the species level, we assigned the mean of the next taxonomic level (such as genus). To obtain nutrient supplies or yields, nutrient concentrations per 100 g were multiplied by equivalent units in live weight. Note that we used live weights for our main analyses (that is, extrapolating the nutrient content from muscle tissue to live weight volumes), but also performed a sensitivity analysis using edible weight (for example, multiplying live weights by edible proportions; Extended Data Fig. 2). We grouped nutrient supplies and total catch or production volumes separately for invertebrates and fish, and calculated the percentage of nutrient supplies coming specifically from invertebrates. We separately did this for different sectors, taxa, years and countries or EEZs. Note that, similar to finfish5, country- or EEZ-specific nutrient yields are not strongly correlated with the average nutrient concentration of their catches (Supplementary Figs. 7 and 8), indicating, for example, that a country or EEZ with high nutrient yields does not necessarily target the species with most nutrient concentrations.
Estimating public health relevance of aquatic foods
Nutrient supplies were converted to the number of yearly requirements met by dividing nutrient-specific nutrient supplies from aquatic invertebrates by their estimated RNI, averaged across available demographic groups (gender and age; Supplementary Table 3). When available, we used the recommended dietary allowance61. However, when such data were not available, we used adequate intake61 estimates or values reported in the literature: 433 mg per day for DHA and EPA (that is, average between different studies reviewed)62 and 48.8 g per day for MUFAs (that is, 20% of total energy recommended intake assuming that a gram of MUFA represents 9 kcal)63,64. Note that for total omega 6 fatty acids, we used recommended intakes reported for linolenic acids; and for total omega 3 fatty acids, we used the sum of ALA and DHA and EPA recommended intakes (that is, 1,633 mg per day). Yearly nutrient supplies were divided by yearly requirements, ensuring that units were the same and assuming 365 days (for example, yearly requirements in tonnes = daily requirements in tonnes × 365). Our approach estimates the number of nutrient requirements met exclusively from the nutrients available from aquatic invertebrates. While we acknowledge that people eat other foods and do not meet their nutritional adequacy exclusively from aquatic invertebrates, this enabled us to compare the importance of aquatic invertebrates across nutrients with a standardized unit relevant for public health.
Predictive model of invertebrate nutrient concentrations
To estimate the variability and potential drivers of nutrient content in aquatic invertebrates, we first updated and validated species-specific invertebrate nutrient concentration samples within AFCD. A total of 13,888 samples from 465 invertebrate species were compiled, updated and validated (for example, reviewing source studies or food composition tables for accuracy, and adding extra information (for example, sample preparation, relative weight)). Specifically, food part (that is, the body part that was sampled for nutrient concentration) and food processing (that is, mechanical or chemical processes that transform the animal to the form before consumption) were validated (for example, examining categorizations for accuracy and disaggregating processing from sample preparation; Supplementary Table 1). Sample sizes and species varied by nutrient (Supplementary Table 2).
We next merged species-specific nutrient data with ecological and environmental trait information available from SeaLifeBase27. Traits were assigned hierarchically, using species-specific data when available or the mean or most common category of the closest taxonomic group (for example, genus). Taxonomic level assignment of traits is shown in Supplementary Fig. 9. We included available traits related to energetic demand, thermal regime, habitat and environment that are likely to influence the nutrient composition of aquatic invertebrates (Supplementary Table 1).
Similar to work done on ray-finned fishes5, we developed a series of Bayesian hierarchical models to predict the nutrient concentration of aquatic invertebrate species on the basis of their taxonomy and ecological and environmental traits, while controlling for other factors thought to impact the nutrient concentration of samples (for example, food part or processing form; Supplementary Table 1). Note that owing to limited species-specific invertebrate sample sizes (Supplementary Table 2), from the nutrients highlighted, molybdenum, vitamin C and D were not included in this modelling section. Given the structure of our data and the importance of taxonomic identity in explaining the variability in nutrient content65, taxonomy was included in our model hierarchically (that is, nested random effects of samples within genus, within families, within orders, within class and within phyla). All other factors (Supplementary Table 1) were included in the model as fixed effects, with (1) maximum depth log-transformed to normalize the spread of a highly skewed distribution (to increase model efficiency); (2) all continuous variables standardized by subtracting the mean and dividing by two times the standard deviation (s.d.); and (3) maximum length and length at maturity standardized within the taxonomic level of class because for invertebrates, different length types are recorded depending on the taxonomic group. Dividing by two times the s.d. facilitates a more direct comparison of regression coefficients between continuous and binary predictors within the same model66. Some variables were slightly correlated (all Pearsons’s correlation coefficients < 0.6). For example, length of maturity was positively correlated to maximum length (Pearson’s correlation coefficient = 0.43) and trophic level was positively correlated with maximum depth (Pearson’s correlation coefficient = 0.57), indicating that (1) invertebrates that reach larger maximum lengths tend to mature at larger lengths; and (2) invertebrates that have higher trophic level tend to reach greater depths. While these correlations influence how marginal posteriors were interpreted, they did not cause problems to model fits.
Nutrients were modelled independently, owing to different sample sizes (Supplementary Table 2). For each nutrient, we tested three alternative model structures through cross validation67: a null model that includes only the intercept; a hierarchical model that includes only the intercepts reflecting the taxonomic nested structure of our data; and a full model that includes taxonomy hierarchically and ecological and environmental traits as covariates. First, we performed leave-out-one cross validation using pareto smoothed importance sampling (PSIS-LOO). However, as this process yielded some high Pareto k diagnostic values, we also performed k-fold cross validation67. To be able to perform k-fold cross validation, for some nutrients, covariate groups with limited sample sizes (for example, exoskeleton for body parts) were merged into another category (for example, whole/mix). Our k-fold cross validation results reinforced those obtained from leave-out-one cross validation: for a minority of nutrients examined (for example, vitamin B12), the inclusion of the taxonomic hierarchy alone predicted the out-of-sample data as well as the full model that includes trait data (Supplementary Tables 5 and 6). However, across nutrients (that is, sum of leave-out-one or k-fold cross validation information criteria assuming response variables are conditional independent), the full model was significantly preferred in terms of predictive accuracy (Supplementary Tables 5 and 6), supporting the inclusion of both taxonomic hierarchy and ecological and environmental traits in predicting nutrient concentration for aquatic invertebrates. Thus, for each nutrient, the best linear model structure was:
$$\begin{array}{l}{\mu }={\beta }_{0,\mathrm{GEN}}+{\beta }_{1}\times {\mathrm{ENVTEMP}}_{\mathrm{subtropical}}+{\beta }_{2}\times {\mathrm{ENVTEMP}}_{\mathrm{temperate}}+{\beta }_{3}\\ \,\,\times {\mathrm{ENVTEMP}}_{\mathrm{cold}}+{\beta }_{4}\times \mathrm{TL}+{\beta }_{5}\times \mathrm{DEPTH}+{\beta }_{6}\times \mathrm{LMAX}+{\beta }_{7}\times \mathrm{LM}\\ \,\,+{\beta }_{8}\times {\rm{K}}+{\beta }_{9}\times {\mathrm{ENV}}_{\mathrm{freshwater}}+{\beta }_{10}\times {\mathrm{ENV}}_{\mathrm{mixed}}+{\beta }_{11}\times {\mathrm{DEMERSPELAG}}_{\mathrm{benthopelagic}}\\ \,\,+{\beta }_{12}\times {\mathrm{DEMERSPELAG}}_{\mathrm{pelagic}}+{\beta }_{13}\times {\mathrm{DEMERSPELAG}}_{\mathrm{demersal}}\\ \,\,+{\beta }_{14}\times {\mathrm{DEMERSPELAG}}_{\mathrm{reef}}+{\beta }_{15}\times {\mathrm{DEMERSPELAG}}_{\mathrm{sessile}}\\ \,\,+{\gamma }_{1}\times {\mathrm{PART}}_{\mathrm{gills}}+{\gamma }_{2}\times {\mathrm{PART}}_{\mathrm{skin}}+{\gamma }_{3}\times {\mathrm{PART}}_{\mathrm{whole}}+{\gamma }_{4}\times {\mathrm{PART}}_{\mathrm{viscera}}\\ \,\,+{\gamma }_{5}\times {\mathrm{PART}}_{\mathrm{reptissue}}+{\gamma }_{6}\times {\mathrm{PART}}_{\mathrm{exoskeleton}}+{\gamma }_{7}\times {\mathrm{PROC}}_{\mathrm{frozen}}+{\gamma }_{8}\times {\mathrm{PROC}}_{\mathrm{dried}}\\ \,\,+{\gamma }_{9}\times {\mathrm{PROC}}_{\mathrm{canned}}+{\gamma }_{10}\times {\mathrm{PROC}}_{\mathrm{boiled}}+{\gamma }_{11}\times {\mathrm{PROC}}_{\mathrm{baked}}+{\gamma }_{12}\\ \,\,\times {\mathrm{PROC}}_{\mathrm{cooked}}+{\gamma }_{13}\times {\mathrm{PROC}}_{\mathrm{unknown}}\\ \,{\beta }_{0,\mathrm{GEN}} \sim N({\beta }_{0,\mathrm{FAM}},{\sigma }_{\mathrm{GEN}})\\ \,\,{\beta }_{0,\mathrm{FAM}} \sim N({\beta }_{0,\mathrm{ORD}},{\sigma }_{\mathrm{FAM}})\\ \,\,{\beta }_{0,\mathrm{ORD}} \sim N({\beta }_{0,\mathrm{CLASS}},{\sigma }_{\mathrm{ORD}})\\ \,\,{\beta }_{0,\mathrm{CLASS}} \sim N({\beta }_{0,\mathrm{PHY}},{\sigma }_{\mathrm{CLASS}})\\ \,\,{\beta }_{0,\mathrm{PHY}} \sim N({\beta }_{0},{\sigma }_{\mathrm{PHY}})\end{array}$$
where β0,…. represents the estimated intercepts at different taxonomic hierarchical levels (GEN, genus; FAM, family; ORD, order; CLASS, class; and PHY, phylum) for average continuous and most common (that is, tropical, benthic and marine) covariate categories that are separate categorical predictors, β0, is the global estimated intercept, β1–15 represents the estimated parameters for different covariates: thermal regime (ENVTEMP, subtropical, temperate or cold), trophic level (TL), maximum depth (DEPTH), maximum length (LMAX), length at maturity (LM), von Bertalanffy growth parameter (K), environment (ENV, freshwater or mixed) and preferred habitat (DEMERSPELAG, benthopelagic, pelagic, reef-associated, sessile or demersal), γ1–13 represent the estimated parameters for different nuisance variables: food part (PART, gills, skin, whole or multiple parts, viscera, reproductive tissue or exoskeleton) and food processing (PROC, frozen, dried, canned, boiled or steamed, baked, grilled or smoked, cooked other or unknown food preparation); and σ… represent the estimated s.d. for the hierarchical taxonomic intercepts.
For nutrients with a limited number of zeroes (<2%; Supplementary Table 2) we only used positive data and used a Student’s t family distribution on the variable’s natural logarithm:
$$\log ({N}_{i}) \sim {\rm{S}}{\rm{t}}{\rm{u}}{\rm{d}}{\rm{e}}{\rm{n}}{\rm{t}}-t(\nu ,\mu ,\tau ),$$
whereas for nutrients with a higher number of zeroes (Supplementary Table 2), we used a hurdle log-normal data likelihood distribution:
$$\begin{array}{l}{\rm{if}}{N}_{i}=0,\,{N}_{i} \sim {\rm{bernouilli\_logit}}\,(\delta )\\ {\rm{if}}{N}_{i} > 0,\,{N}_{i} \sim {\rm{lognormal}}\,(\mu ,\,\sigma )\end{array}$$
where Ni is the nutrient concentration of sample i, µ is the mean nutrient concentration informed by the linear model structure above (in log scale), ν and τ are the degrees of freedom and scale parameters for the Student’ t distribution, δ is the estimated probability of observing a zero for a given nutrient if it had >2% of zeroes, and σ is the s.d. of the variable’s natural logarithm.
Models were run in Rstan68 through the brms package69 and using the following priors:
$$\begin{array}{c}{\beta }_{0} \sim N(0,10)\\ \beta .. \sim N(0,2)\\ \gamma .. \sim N(0,2)\\ \nu \sim \mathrm{gamma}(2,0.1)\\ \delta \sim \mathrm{beta}(1,1)\\ \sigma .. \sim {N}^{+}(0,1)\\ \tau \sim {N}^{+}(0,1)\end{array}$$
Four chains were run for each scenario using 10,000 iterations (5,000 warmup and a thinning of 5), leaving 4,000 samples in the posterior distribution of each parameter. Convergence was monitored by running four chains from different starting points, examining posterior chains and distribution for stability, checking that the potential scale reduction factor (also termed R hat) was close to 1 (below 1.01) and examining the effective sample sizes (>400). Some parameters (such as the s.d. of hierarchical structures) were informed by limited sample sizes; we therefore used a large number of iterations to ensure parameters informed by limited sample sizes exceeded the recommended effective sample sizes. Bayesian learning was examined by inspecting posteriors vs. prior distributions and by calculating posterior contraction values70. We examined posterior predictive distributions to check for model fit. Moreover, we checked model residual patterns against variables not included in our model, owing to high missingness (Supplementary Figs. 5 and 10). All final models reported converged, fit the data relatively well, had high posterior contraction for key parameters and had relatively well performing probability integral transformations (PIT) plots (Supplementary Figs. 11–14).
Furthermore, to test whether (1) models of less complexity (for example, with less fixed effects) had better predictive accuracy than the full model; and (2) there were a minimal set of variables which could provide similar predictions to the full model, we performed projection predictive variable selection assuming positive data followed a log-normal distribution71. This revealed that (1) subsets of the full model did not provide better predictive accuracy than the full model (that is, the expected log-predictive density of the full model was always better or the same as simpler models); (2) the ranking of covariates being selected as most important to predict nutrient concentrations for invertebrates differed across nutrients; and (3) for some nutrients excluding some covariates could provide similar predictive performance to the full model, but the covariates to exclude varied by nutrient (Supplementary Fig. 1). Combined, these analyses supported the use of our full model to predict concentrations in invertebrates across nutrients.
Predicting nutrient composition of invertebrates in SeaLifeBase
Finally, we used model parameters estimated in the above section to predict nutrient composition of all macroinvertebrates available in SeaLifeBase27. While SeaLifeBase probably underestimates all invertebrate species (for example, freshwater invertebrate species), it is the most comprehensive life-history trait database available for non-fish species. The initial database contained a total of 71,591 species. From these, in total, 55,864 species were identified as macroinvertebrates on the basis of taxonomic groupings: 281 were confirmed as edible based on literature72,73,74,75,76,77 and 50,807 were identified as potentially edible (excluding 834 non-edible species because of known toxicity and 4,223 species from the phyla Bryozoa, Platyhelminthes, Porifera and Priapulida, which we assumed were not edible, owing to the absence of documented human consumption and/or biological constraints, such as chemical defences, toxicity or parasitic lifestyles, that limit their suitability for consumption (for example, antipredatory metabolites in bryozoans and sponges or tetrodotoxin in some flatworms)). This designation does not imply cultural acceptance, safe consumption or practical harvestability across contexts. On the basis of available literature25,27,72,73, some potentially edible macroinvertebrate species (1,083) were also categorized on the basis of their fisheries importance into the following groups: industrial, highly commercial, commercial, minor commercial, subsistence fisheries and bycatch. We provide nutrient composition estimates for all invertebrate species; but, in our analyses, we filtered the database to include only species classified as potentially edible macroinvertebrates (Fig. 3). Note that the database also included ascidians (phylum Chordata), which are not invertebrates. However, they were included in model predictions given (1) similarities in some ecological and environmental traits with invertebrates, and (2) that some are edible78.
Some species had missing information. However, to perform predictions for all available species we took some assumptions. For example, if a specific group (such as Cnidaria) did not have a group-specific estimated intercept (that is, we did not have that group for that nutrient when building the Bayesian hierarchical models), we used the upper-level intercept (for example, global intercept) in combination with ecological and trait information. Similarly, if a species was missing some ecological or environmental trait information, we assumed the baseline category (for example, mean for continuous variables and most common separate categories: benthic, marine and tropical). To represent the data visually, we divided predicted nutrient concentrations by per capita daily RNIs (see above), restricting values to one when they were above such value (that is, if 100 g of muscle tissue provided over 100% of RNIs).
Ultimately, we used matched predictive nutrient concentrations in terms of RNIs of all potentially edible macroinvertebrate species to species-specific marine occurrence data35 to show the geographical public health potential of macroinvertebrates (that is, what nutrient concentrations from invertebrates may be available in different geographical contexts).
Case studies
We developed six case studies representing a diverse array of geographies, target species, spatial scales and management contexts to highlight the multiple social and ecological benefits of aquatic invertebrates, including their environmental, economic, cultural and equity related benefits (Fig. 4). These case studies highlight the indirect contribution of aquatic invertebrates to nutrition by linking ecological functions and production systems to social outcomes such as income generation, gender equality and climate resilience.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
