Macroscale patterns of oceanic zooplankton composition and size structure

Ocean plankton comprise organisms from viruses to fish larvae that are fundamental to ecosystem functioning and the provision of marine services such as fisheries and CO2 sequestration. The latter services are partly governed by variations in plankton community composition and the expression of traits such as body size at community-level. While community assembly has been thoroughly studied for the smaller end of the plankton size spectrum, the larger end comprises ectotherms that are often studied at the species, or group-level, rather than as communities. The body size of marine ectotherms decreases with temperature, but controls on community-level traits remain elusive, hindering the predictability of marine services provision. Here, we leverage Tara Oceans datasets to determine how zooplankton community composition and size structure varies with latitude, temperature and productivity-related covariates in the global surface ocean. Zooplankton abundance and median size decreased towards warmer and less productive environments, as a result of changes in copepod composition. However, some clades displayed the opposite relationships, which may be ascribed to alternative feeding strategies. Given that climate models predict increasingly warmed and stratified oceans, our findings suggest that zooplankton communities will shift towards smaller organisms which might weaken their contribution to the biological carbon pump.


Supplementary
: Maps and latitudinal patterns of the abundance (cubictransformed ind.m 3 ) of the zooplankton groups sampled with the three plankton nets (WP2, Bongo and Régent). The solid curves on the zonal side plots illustrate the prediction from the Generalized Additive Model (GAM) fitting transformed abundance against latitude. The explanatory power of the GAM (adjusted R 2 ), the number of samples used and the significance of the smooth term (p < 0.001 = ***, p < 0.01 = ***, p < 0.05 = *, p > 0.05 = ns) are reported on the plots. The grey ribbon illustrates the standard error of the prediction. Supplementary Document S3: Distribution of the relative contribution (%) of the main mesozooplankton groups to total zooplankton community abundance across the three plankton nets used.

Supplementary Document S4:
Comparing mesh size and sampling depth effects on the abundance distributions of the main zooplankton groups. Figure S4.1: Comparing the cubic-transformed abundance distributions of the main zooplankton groups studied between the three plankton mesh sizes (200µm, WP2; 300µm, Bongo; 680µm, Régent) and sampling depth of the net hauls (0-100m, WP2; 0-500m, Bongo and Régent). a) Total zooplankton, b) Copepoda, c) Rhizaria, d) Cnidaria, e) Tunicata, f) Eumalacostraca, g) Pteropoda, h) Chaetognatha, and i) Ostracoda + Cladocera. For each group and mesh size, non parametric variance analyses (Kruskal-Wallis tests) were performed to test if the nets displayed significant different zooplankton abundances. Conventional labels illustrate the level of significance of the test: n.s. = not significant; ** = p-value < 0.01; *** = p-value < 0.001. Figure S4.2: Pairwise correlations of cubic-transformed abundance distributions of the main zooplankton groups studied between the sampling depth of the net hauls (0-100m, WP2; 0-500m, Bongo). a) Total zooplankton, b) Copepoda, c) Rhizaria, d) Cnidaria, e) Tunicata, f) Eumalacostraca, g) Pteropoda, h) Chaetognatha, and i) Ostracoda + Cladocera. The results from the Spearman's rank correlation tests are shown on the plots.  Tables 1 and 2. Ranks range between 0 (least important covariate) and 1 (most important covariate). The F statistic of the smooth term associated to each covariate was extracted and was used to quantitatively rank the covariates according to their relative significance (the higher the F value, the more significant) within their corresponding GAM. For each GAM, ranks were normalized by dividing the F value of every smooth term by the maximum F value across al smooth terms. The complete statistics of the GAMs are given in Supplementary Table  S7.    Any sampling station whose spatial coordinates fell into one of the four EBUS was labelled as a "EBUS station". The corresponding abundance estimates of the main zooplankton groups were compared to those stemming from the "non EBUS stations" that displayed a similar latitudinal range (i.e. 29°N-46°N and -30°N-1°N) as the EBUS stations (see below). This way, high latitude sampling stations were discarded from this analysis so we can focus on those abundance patterns that mostly reflect effects due to the EBUS conditions.

California
Benguela NW Africa Peru Figure S12.2: Comparing the cubic-transformed abundance distributions of the main zooplankton groups studied across three plankton nets (WP2,200µm;Bongo,300µm;Régent,680µm) between the EBUS stations in blue (see Fig. SXX1) and the comparable non EBUS stations in red (i.e. stations located between 29°N -46°N and -30°N-1°N). a) Total zooplankton, b) Copepoda, c) Rhizaria, d) Cnidaria, e) Tunicata, f) Eumalacostraca, g) Pteropoda, h) Chaetognatha, and i) Ostracod + Cladocera. For each group and net type, non parametric variance analyses (Kruskal-Wallis tests) were performed to test if the EBUS stations displayed significant higher zooplankton abundances. The groups and nets for which a significant (p-value < 0.05) was found between EBUS and non EBUS stations are highlighted with a black star. Conventional labels illustrate the level of significance of the test: * = p-value < 0.05; ** = p-value < 0.01; *** = p-value < 0.001. Figure S12.3: Maps and latitudinal patterns of the abundance (cubic-transformed ind.m3) of the main zooplankton groups sampled with three plankton nets (WP2,200µm;Bongo,300µm;Régent,680µm), excluding those sampling stations that fall within an eastern boundary upwelling system (EBUS). The solid curves on the zonal side plots illustrate the prediction from the Generalized Additive Model (GAM) fitting transformed abundance against latitude. The explanatory power of the GAM (adjusted R2), the number of samples used and the significance of the smooth term (p < 0.001 = ***, p < 0.01 = ***, p < 0.05 = *, p > 0.05 = ns) are reported on the plots. The grey ribbon illustrates the standard error of the prediction. As in Supplementary Figure S2, the order of appearance of the groups' maps+plots panel is always as follows: WP2 data, Bongo data and then Régent data.
The effect of the EBUS stations on the cubic-transformed abundance latitudinal patterns modelled above can be examined by comparing these maps and zonal latitudinal plots to their counterparts shown in Supplementary Figure S2. Very few differences were found, which suggests that EBUS have a relative weak effect on the macroscale patterns of zooplankton abundances. However, we acknowledge that the sampling design of the Tara Oceans cruises does not allow to fully investigate this effect, as no stations sampled the northwestern African EBUS and relatively few stations were taken in the other three EBUS (Fig. S12.1). The most notable differences in latitudinal abundance patterns were found for the:  Copepoda (WP2): the modelled increase in abundances near the equator was way less marked as the relatively high abundances observed in the Peru and Benguela EBUS were discarded.  Rhizaria (WP2): the modelled decrease in abundances near 30°S was removed as the relatively low abundances observed in the Benguela EBUS were discarded.  Tunicata (WP2): the modelled increase in abundances near the equator was less marked as the relatively high abundances observed in the Peru EBUS were discarded.
All the other abundance latitudinal patterns modelled were left unchanged by the removal of the EBUS stations. Therefore, although the EBUS seem to have a positive influence of the abundances of most of the main zooplankton groups studied here ( Fig. S12.2), they do not strongly impact the modelled latitudinal gradients in zooplankton abundance.