Decoding pelagic ciliate (Ciliophora) community divergences in size spectrum, biodiversity and driving factors globally spanning five temperature zones

The community structure of microzooplanktonic ciliates – encompassing size spectrum, biodiversity and biotic–abiotic interplay – is critical for unravelling their ecological role in marine ecosystems, yet it remains challenging to elucidate on a global scale. To address this knowledge gap, we conducted field observational studies across five temperature zones (North Frigid Zone, NFZ; Sub-Arctic Zone, SAZ; North Temperate Zone, NTZ; Torrid Zone, TZ; South Frigid Zone, SFZ). Our analysis demonstrates a sharp decline in ciliate abundance and biomass below the 100 m layer, with distinct vertical distribution patterns observed in each climate region. Moreover, although abundance of ciliate size spectra exhibited a decrease trend from small to large size spectra globally, there were steeper slope lines observed in both polar zones (NFZ and SFZ) compared to the other temperature zones. Latitudinally, ciliate abundance and tintinnid biodiversity exhibited an anti-phase relationship, where the TZ hosted peak biodiversity, while the polar seas showed the highest abundance. Furthermore, a multivariate biota–environment analysis indicated that temperature has a primary influence on ciliate community constitution in the global marine ecosystem, and the bottom-up control plays a key role in shaping assemblages. In conclusion, these results underscore the unprecedented divergences in ciliate trait structure among five temperature zones and can be taken as a guideline for assessing the potential effects of climate change on pelagic ciliates in future marine realms.

1 Introduction

The Earth is traditionally partitioned into five temperature zones based on established climate classifications: the North Frigid Zone (NFZ), North Temperate Zone (NTZ), Torrid Zone (TZ), South Temperate Zone (STZ) and South Frigid Zone (SFZ) (Köppen, 1936; Trewartha et al., 1967). Therein, each temperature zone possessed a unique ocean circulation pattern and concurrent specific plankton biome structures (Longhurst, 2007; Spalding et al., 2012). Although a myriad of prevailing research exists relevant to plankton biogeography and its interplay with environmental drivers highlighting its importance in disentangling marine ecosystems and biogeochemical cycles (e.g. Wang et al., 2020; Darnis et al., 2022; Segaran et al., 2023; Tagliabue et al., 2023), substantial global-scale studies have predominantly relied on modelling frameworks (Spalding et al., 2012; Blanchard et al., 2017; Anderson et al., 2021; Benedetti et al., 2021; Heneghan et al., 2023; Atkinson et al., 2024). To date, an explicit and comprehensive representation of plankton community trait structure using data-derived statistical analysis originating from field surveys remains unresolved.

A holistic paradigm of plankton biogeography across marine ecosystem is crucial for deciphering global ecological connectivity (Hillman et al., 2018) and predicting how ecosystems respond to stressors induced by climate change (Darnis et al., 2022). Over recent decades, anthropogenic CO₂ emissions have led to increased atmospheric concentrations and greater global radiative forcing (IPCC, 2023), triggering diverse ecological feedbacks worldwide, for instance poleward distribution shifts (Neukermans et al., 2018; Oziel et al., 2020; Benedetti et al., 2021), adjustments in phenology (Poloczanska et al., 2013; Atkinson et al., 2015; Chust et al., 2024) and reductions in mean body size (Daufresne et al., 2009; Verberk et al., 2021; Wang et al., 2023a, b). In this sense, extensive existing studies put emphasis on biotic community response to climate change in the polar and adjacent seas, owing to their higher susceptibility compared to tropical, subtropical and temperate seas (Serreze et al., 2009; Screen and Simmonds, 2010; IPCC, 2023; Noh et al., 2024). Unfortunately, informative research related to environmental affinity of plankton, particularly microzooplankton, is not sufficiently understood in the aforementioned five temperature zones.

In the realm of microzooplankton, pelagic ciliates stand out as the predominant biological entities, spanning in size from 10 to 200 µm, and hold significant sway over both biodiversity and abundance, particularly in the polar and adjacent seas (Taniguchi, 1984; Strom and Fredrickson, 2008; Lu and Weisse, 2022; Kohlbach et al., 2023; Wang et al., 2023b, 2024a, b). Taxonomically categorized within the phylum Ciliophora, class Spirotrichea, and subclasses Oligotrichia and Choreotrichia, pelagic ciliates, including aloricate ciliates and tintinnids, are ubiquitous single-cell protozoans found in various aquatic environments worldwide (Lynn, 2008). Furthermore, ciliates play an irreplaceable role in marine trophodynamics (carbon cycle and energy transfer) through prey–predator interactions, serving as both phytoplankton grazers and prey for metazoans (Stoecker et al., 1987; Dolan et al., 1999; Calbet and Saiz, 2005; Gómez, 2007; Weisse and Sonntag, 2016). Specifically, owing to their simple life cycle, fast reaction to environmental changes and strong adaptability, pelagic ciliates, particularly tintinnids, are widely recognized as ideal bioindicators for assessing various sea conditions (e.g. Kato and Taniguchi, 1993; Jiang et al., 2013; Wang et al., 2021; Yu et al., 2022).

Recent escalation in global warming has imposed a cascade of impacts on aquatic ecosystems, presenting a formidable challenge to inherent holopelagic species that modify their relevant adaptive strategies (Stabeno et al., 2012; Yasumiishi et al., 2020; Carvalho et al., 2021; Atkinson et al., 2024). Accordingly, a prevailing viewpoint on phytoplankton, the cornerstone of the marine pelagic food web, is a major decline in both biomass and size spectra in the NTZ, TZ and STZ (Li et al., 2009; Lotze et al., 2019; Tittensor et al., 2021), leading to subsequent declines in higher trophic levels, termed “trophic amplification” (Kwiatkowski et al., 2019; du Pontavice et al., 2021). As a grazer of pelagic phytoplankton, response of microzooplanktonic ciliates to ocean warming in the polar and adjacent seas is substantial (Li et al., 2022; Wang et al., 2022a, 2023a, b, 2024a), yet comparative assessments amid their trait structure (e.g. size spectra, biodiversity and biotic–abiotic interplay) remain unexplored to date.

Consequently, elucidating microzooplanktonic ciliate size spectra, species diversity and biotic–abiotic interplay at a global scale is critical for projecting future marine ecosystem dynamics, particularly given their unresolved role in plankton response to climate changes. Here, we propose a hypothesis that hydrographic variability is likely responsible for the observed divergence in global ciliate trait structures. By optimizing field observational data and available methods, this study aims (1) to decode adaptive strategies of microzooplanktonic ciliate to heterogeneous hydrographic conditions across temperature zones and (2) to evaluate their potential response dynamics to accelerating climate change. Given the current foreseeable rapid climate change, this study will offer a benchmark for facilitating the phenological and bioclimatic progression of microzooplankton shifts in future global marine ecosystem realms.

2 Materials and methods

2.1 Study area and field sampling

Based on their latitudinal locations, field samplings of microzooplanktonic ciliates were conducted in five temperature zones (Trewartha et al., 1967): (1) the North Frigid Zone (NFZ), encompassing the Arctic Ocean, during July to August 2019 and 2023 aboard the R/V Xiang Yang Hong 1 and R/V Xue Long 2, respectively; (2) the Sub-Arctic Zone (SAZ), located in the Bering Sea, in July to August 2019 aboard the R/V Xiang Yang Hong 1; (3) the North Temperate Zone (NTZ), situated in the North Pacific, in September 2019 aboard the R/V Dong Fang Hong 3; (4) the Torrid Zone (TZ), which includes the tropical western Pacific in December 2016 and August 2017 aboard the R/V Kexue and the Indian Ocean in March 2021 aboard the R/V Xiang Yang Hong 6; and (5) the South Frigid Zone (SFZ), covering the Southern Ocean, from December 2020 to March 2021 aboard the R/V Xue Long 2 (Fig. 1). A total of 1117 samples (175 stations along 19 transects) were sampled.

Figure 1 Survey stations and transects (Tr.) in the tropical, temperate and polar seas. NFZ, North Frigid Zone; SAZ, Sub-Arctic Zone; NTZ, North Temperate Zone; TZ, Torrid Zone; SFZ, South Frigid Zone.

Seawater samples were collected with a rosette sampler carrying 24 Niskin bottles (each 12 L). All microzooplanktonic ciliate samples were collected at seven standardized depths (surface (2 m), 25, 50, 75, 100, 150 and 200 m) at each designated station, with the exception of SAZ stations, where bathymetry limited sampling to depths <200 m. Furthermore, each sample was fixed with Lugol's acid (1 % final concentration) and preserved in darkness at 4 °C until laboratory analysis.

2.2 Sample analysis

Laboratory processing involved concentrating each sample to approximately 200 mL through siphon-assisted supernatant removal following 60 h sedimentation. After two rounds of the siphon process, finally a 25 mL highly concentrated sample was obtained and then settled in a Utermöhl counting chamber (Utermöhl, 1958). Quantitative analysis was performed using an Olympus IX71 inverted microscope (100× or 400× magnification) to enumerate total ciliate abundance (including aloricate ciliates and tintinnids), measure morphometric parameters (body size) and document species richness across all five temperature zones by Chaofeng Wang. To ensure accuracy, cellular dimensions (e.g. length, width, shape) of aloricate ciliate or each tintinnid species were measured for at least 10 individuals if possible.

Additionally, the body sizes of both aloricate ciliates and tintinnids were categorized into 10 µm increments (10–20 µm, 20–30 µm, etc.) based on body length (Wang et al., 2020) and further classified into small (10–20 µm), medium (20–50 µm) and large (>50 µm) size fractions following Yang et al. (2019). Moreover, we did not distinguish the presence/absence of tintinnid lorica during the sample-counting process. Regarding species richness, tintinnid identification was assigned to the closest species, as described in Zhang et al. (2012). Furthermore, we select the average value (15, 25, 35, 45 µm etc.) of each size fraction of both loricate ciliate and tintinnid as the counting criterion for ciliate size spectra (Wang et al., 2024a). In addition, the slope or slope line means tendency of evaluating the decreasing trend from small to large size spectrum. Simultaneously, environmental factors of sampling depth (a pressure sensor to detect hydrostatic pressure, converted to depth via the formula $\text{depth}=\text{pressure}/[\mathit{\rho }\times g]$, where ρ is water density and g is gravitational acceleration) (van Haren et al., 2021), temperature (a thermistor, SBE-3 Plus; resolution is 0.0001 °C), salinity (derived from measured electrical conductivity (SBE-4C sensor) and temperature data, computed using the practical salinity scale algorithm) and chlorophyll a in vivo fluorescence (chl a, a fluorometer [SeaPoint] excites chlorophyll pigments with blue light and measures emitted red light intensity as a proxy for chl a concentration) were recorded by a multi-sensor profiler (CTD – SeaBird SBE 911, https://www.seabird.com/product.detail-cms.block.jsa?id=60761421595, last access: 5 June 2025) during each cruise.

2.3 Data processing

Ciliate volumes were estimated according to their appropriate geometric shapes (cone, ball, cylinder). Carbon biomass of each tintinnid was calculated by the following equation (Verity and Lagdon, 1984):

$$C=V_i\times \mathrm{0.053}+\mathrm{444.5},$$

where C (10⁻⁶ µg C) is the carbon biomass of individual tintinnid, and V_i (µm³) is the lorica volume. Additionally, a conversion factor ($\mathrm{0.19}\times {\mathrm{10}}^{-\mathrm{6}}$ µg C µm⁻³) was used for calculating aloricate ciliate carbon biomass (Putt and Stoecker, 1989). Concerning size spectra biomass, ciliate biomass was calculated based its specific organism volume and conversion equation and then categorized into each size spectrum as in Wang et al. (2024a). Furthermore, in order to better unravel tintinnid biodiversity spanning five temperature zones, the Margalef index (d_Ma) (Margalef, 1958) (1) and Shannon index ($H_{\mathrm{2}}^{\prime }$) (Shannon, 1948) (2) were conducted by the following equations:

$$\begin{array}{}\text{(1)}& d_{\text{Ma}}={\displaystyle \frac{S-\mathrm{1}}{\ln N}},\end{array}$$

where S is the number of species, and N is the total number of tintinnid individuals in the sample.

$$\begin{array}{}\text{(2)}& H_{\mathrm{2}}^{\prime }=-{\sum }_{i-\mathrm{1}}^S{P}_i{\log }_{\mathrm{2}}{P}_i,\end{array}$$

where S is the number of species, and N is the total abundance of tintinnid individuals in the sample. P_i ($N_i/N$) is the relative abundance of i species in a whole community.

Biogeographically, classification of tintinnid genera (cosmopolitan, warm water, boreal, austral and neritic) was based on Pierce and Turner (1993) and Dolan and Pierce (2013). Among them, tintinnid genera were further classified into oceanic (cosmopolitan, warm water, boreal and austral) and neritic types. Moreover, the average value of each parameter is represented as mean ± SD in the following text. Finally, seasonality is important to modulate protozoan communities, but this phenomenon was obvious in both temperate and polar seas. Regarding tropic seas in both the Pacific and Indian oceans, the community structure, including vertical distribution pattern, abundance and biomass values and species composition, was almost same (e.g. Sohrin et al., 2010; Li et al., 2018; Wang et al., 2019a, 2020, 2022c).

Hereinafter, sampling map was visualized by ODV (Ocean Data View, Version 4.7), and ciliate distributional data of size–diversity and temperature–diversity relationships were analysed using Surfer (Version 13.0), Grapher (Version 12.0) and OriginPro 2021 (Version 9.6). Moreover, in order to reduce deviation in the biotic–abiotic relationship in different temperature zones that may be mainly caused by the difference in the selection of sampling areas, rather than the fundamental differences between temperature zones, the internal correlations among each temperature zone at specific sampling depth (0, 50, 100 and 200 m) were compared in the following text. Meanwhile, the biota–environment analysis was performed based on Spearman's correlation between log-transformed abiotic parameters and square root-transformed abundance data (t test) using both PRIMER (Version 5.0) and OriginPro 2021 (Version 9.6). Additionally, the slope of the size spectrum (a straight line fitted through the size spectrum on a log–log plot) (Blanchard et al., 2017) was carried out to quantize its interplay with ciliate abundance at a discrete depth of aforementioned global seas (95 % confidence). In the following, based on the slope condition, we used the decreasing rate (Δ_D) or increasing rate (Δ_I) according to ciliate abundance or species richness and environmental variables to quantize their interplay in the global seas.

Figure 2 Variations in environmental variables and ciliate abundance and biomass at discrete depths (2, 50, 100 and 200 m) in each temperature zone.

3 Results

3.1 Hydrography and ciliate abundance and biomass

Each environmental parameter (temperature, salinity and chl a) displayed distinct spatiotemporal variations globally (Figs. 2 and S1–S3 in the Supplement). Horizontally, at the surface and 50 and 100 m layers, both temperature and salinity peaked in the Torrid Zone (TZ), contrasting with chl a, which exhibited its lowest value in the same region (Figs. 2, S1 and S2). At 200 m depth, temperature peaked in the TZ, and chl a peaked in the North Frigid Zone (NFZ), contrasting with salinity patterns, which displayed high values in both the TZ and NFZ (Figs. 2 and S1). Vertically, both temperature and chl a declined in the NFZ and Sub-Arctic Zone (SAZ) (surface-peak pattern), while salinity increased from the surface to 200 m layers across all regions (Figs. S1–S3). Moreover, temperature displayed a low–high–low structure at inner stations of the South Frigid Zone (SFZ), and chl a peaked at subsurface layers in both the North Temperate Zone (NTZ) and TZ (Fig. S1).

Pelagic ciliate abundance ranged from 22–9142 $\mathrm{ind}.{\mathrm{L}}^{-\mathrm{1}}$ in the NFZ, 182–9242 $\mathrm{ind}.{\mathrm{L}}^{-\mathrm{1}}$ in the SAZ, 65–886 $\mathrm{ind}.{\mathrm{L}}^{-\mathrm{1}}$ in the NTZ, 25–436 $\mathrm{ind}.{\mathrm{L}}^{-\mathrm{1}}$ in the TZ and 44–5866 $\mathrm{ind}.{\mathrm{L}}^{-\mathrm{1}}$ in the SFZ, whereas their biomass ranged from 0.0–39.3 µg C L⁻¹, 0.3–24.0 µg C L⁻¹, 0.1–1.1 µg C L⁻¹, 0.0–1.1 µg C L⁻¹ and 0.0–26.1 µg C L⁻¹ in aforementioned regions, respectively (Figs. 2 and S1–S3; Table S1 in the Supplement). Horizontally, both high abundance (≥2000 $\mathrm{ind}.{\mathrm{L}}^{-\mathrm{1}}$) and biomass (≥5.0 µg C L⁻¹) of ciliates were observed in surface layers of the NFZ, SAZ and SFZ, coinciding with high chl a levels. At 50, 100 and 200 m layers, the SAZ and TZ had the highest and lowest abundance, respectively (Figs. 2 and S1). Vertically, both ciliate abundance and biomass exhibited a surface-peak pattern in the NFZ, SAZ and SFZ, whereas in the NTZ and TZ, this pattern transitioned to subsurface-peak and bimodal-peak distributions, respectively (Figs. S1 and S2).

Meanwhile, aloricate ciliates dominated the ciliate community, accounting for ≥90 % of total abundance at each depth in the NFZ, NTZ, TZ and SFZ. However, in the SAZ, tintinnid played a more significant role in the ciliate community, with an average relative abundance at most sampling depths exceeding 10 % (Fig. S4). In terms of aloricate ciliates in the horizontal direction, small (10–20 µm) and medium (20–50 µm) size fractions in the SAZ exhibited the highest average abundance at the surface and 50, 100 and 200 m layers, whilst the largest (>50 µm) size fraction had the highest average abundance at the surface and 50 and 100 m layers in the SFZ (Fig. S5). Additionally, except for the NTZ, the abundance and relative abundance of the medium size fraction were highest in the other four regions at both the surface and 50 m layers. At 200 m depth, the small size fraction predominated among the aloricate ciliates (Figs. S5). Vertically, the relative abundance of the large size fraction (>50 µm) exhibited a decreasing trend, whereas the small size fraction displayed an increasing trend across the five temperature zones (Figs. S5).

Figure 3 Variations in body-size spectra of ciliate normalized abundance at discrete depth (2, 50, 100 and 200 m) in each temperature zone.

3.2 Notable variations in pelagic ciliate size spectrum composition

The abundance and biomass of pelagic ciliate size spectra displayed significant variations across global seas (95 % confidence) (Figs. 3 and 4). Generally, the slopes of the normalized abundance size spectra varied from −2.13 to −0.87 (average −1.60 ± 0.33), and relevant biomass values varied from −0.99 to −0.08 (average −0.53 ± 0.25), with the former slope line much steeper than the latter (Fig. 3). Therein, ciliate abundance decreased from small (15 µm) to large size spectra (>100 µm), with the slope line of the normalized abundance size spectra in both the NFZ (−2.13 to −1.93, average −2.01 ± 0.09) and SFZ (−2.01 to −1.63, average −1.80 ± 0.17) being steeper than in the other three regions at each depth (Fig. 3). Additionally, a secondary peak in abundance, featuring large size spectra (>100 µm), was observed at the surface layers of the NFZ, SAZ and SFZ (Fig. 3).

Figure 4 Variations in body-size spectra of ciliate normalized biomass at discrete depth (2, 50, 100 and 200 m) in each temperature zone.

In contrast, the distribution characteristics of ciliate biomass within size spectra did not align with the abundance trend (Fig. 4). Notably, the 65 µm size spectrum exhibited the highest values at both surface and 50 m layers of the NFZ, followed by the SFZ (55 µm) and SAZ (55 µm), with the TZ (35 µm) and NTZ (25 µm) showing lower values (Fig. 4). Moreover, the slope lines of the normalized biomass size spectra in the SFZ (−0.99 to −0.77, average −0.86 ± 0.10) were steeper than those in the SAZ (−0.74 to −0.43, average −0.62 ± 0.13), NTZ (−0.63 to −0.44, average −0.53 ± 0.09), TZ (−0.74 to −0.25, average −0.47 ± 0.22) and NFZ (−0.37 to −0.08, average −0.21 ± 0.12) (Fig. 4). Interestingly, the highest biomass of ciliate size spectra at the surface and 50 and 100 m layers of the TZ corresponded to the 35 µm size spectrum, while at the 200 m layer, the 15 µm size spectrum became dominant (Fig. 4).

3.3 Dynamics in tintinnid species richness and diversity indices

Tintinnid assemblages exhibited significant spatial heterogeneity in both species richness and diversity metrics (Margalef index d_Ma and Shannon index $H_{\mathrm{2}}^{\prime }$ are quantitative measures of species richness in ecological communities) across five temperature zones (Figs. 5 and S6). Horizontally, species richness, Margalef index (d_Ma) and Shannon index ($H_{\mathrm{2}}^{\prime }$) were notably high at discrete layers in both the NTZ and TZ, followed by the SAZ, NFZ and SFZ (Figs. 5a and S6). To enable cross-regional comparison, we excluded neritic genera (restricted to SAZ and NFZ) from species richness calculations, revealing higher species richness in the SFZ versus NFZ (Fig. 5a). Vertically, elevated values of tintinnid species richness, d_Ma and $H_{\mathrm{2}}^{\prime }$ were primarily observed in the upper 50 m waters of the NFZ, SAZ and SFZ, while these values peaked at 75 and 100 m in the NTZ and TZ, respectively (95 % confidence) (Fig. 5b). Notably, we observed an inverse relationship between ciliate abundance and tintinnid species richness across five temperature zones (Fig. S7), suggesting potential competitive exclusion or niche partitioning dynamics.

Figure 5 Variations of tintinnid species richness, Margalef index (d_Ma) and Shannon index ($H_{\mathrm{2}}^{\prime }$) in latitudinal (a) and vertical (b) direction of all regions.

3.4 Biotic–abiotic interplay and its variations

Ciliate abundance and tintinnid species richness exhibited varying correlations with environmental parameters across the five temperature zones (Figs. 6 and S8–S10). In terms of the biotic–abiotic interplay trend, our results revealed that only the NFZ and SAZ exhibited an increasing trend (Δ_I≥0.03) in abundance–temperature correlation at both surface and 50 m layers compared to other three temperate zones (Fig. S9). Concerning all sampling layers, only the SFZ, differing from the trends observed in the other four temperature zones, displayed a decrease in ciliate abundance with increasing temperature (${\mathrm{\Delta }}_{\mathrm{D}}=-\mathrm{0.26}$, R²=0.06) (Fig. S10). Moreover, only the TZ and SFZ exhibited an increase (Δ_I≥0.29) and a decrease (${\mathrm{\Delta }}_{\mathrm{D}}\le -\mathrm{0.01}$) trend at each sampling layer in abundance–salinity correlation, respectively (Fig. 6b). Furthermore, only SFZ showed an increase (Δ_I≥0.02) trend at each sampling layer in abundance–chl a correlation (Fig. S8), which was aligned with trends in other four temperature zones at all sampling layers (Δ_I≥0.06) (Fig. S10). Regarding species richness–temperature correlation, the highest increasing trend occurred at 50 m of the NFZ (Δ_I=0.26, R²=0.44), while the highest decreasing trend was found at 100 m of the SAZ (${\mathrm{\Delta }}_{\mathrm{D}}=-\mathrm{0.28}$, R²=0.09) (Fig. S9). As for all sampling layers, only the NFZ and TZ exhibited an increasing trend in species richness–temperature correlations, with the former (Δ_I=0.15, R²=0.26) being higher than the latter (Δ_I=0.06, R²=0.23) (Fig. S10). Moreover, concerning biotic–salinity correlations, only the SAZ exhibited an increasing (Δ_I≥0.06) trend at each sampling layer (Fig. S9). In addition, only the polar seas exhibited an increasing trend (Δ_I≥0.01) in species richness–chl a correlation at each sampling layer (Fig. S9).

Figure 6 Variations in principal component analysis (PCA) (a) and Spearman's rank correlation (b) between environmental parameters (depth; temperature, T; salinity, S; chl a) and ciliate (tintinnid, Tin; aloricate ciliate, AC; total ciliate, TC; tintinnid species richness, SR) in five regions. The x axis is the first PCA axis, and the y axis is the second PCA axis. Environmental variables and ciliates are indicated by red lines and black lines, respectively. Grey dots are sampling points. ${}^{\ast \ast }$: p<0.01, ^∗: p<0.05, t test.

To further quantize the physical–biological interplay in five temperature zones, we conducted both principal component analysis (PCA) and Spearman's rank correlation using abundance of aloricate ciliate, tintinnid and total ciliate, as well as tintinnid species richness to test abiotic influence (Fig. 6). The PCA revealed that two principal components effectively differentiated the environmental conditions among five temperature zones. These components accounted for a substantial proportion of the biotic variation in the NFZ (62.85 %), SAZ (67.83 %), NTZ (64.75 %), TZ (72.68 %), SFZ (63.84 %) and all regions (61.42 %) (Fig. 6a). Akin to PCA, Spearman's rank correlation reflected that abundance of aloricate ciliate, tintinnid and total ciliate in all five temperature zones displayed a strong significant negative and positive correlation with depth (p<0.01) and chl a (p<0.01), respectively (Fig. 6b). Furthermore, both aloricate ciliate and tintinnid featured a significant positive correlation with temperature in the SAZ, NTZ and TZ (p<0.05). However, in the SFZ, the relationship between aloricate ciliate and temperature shifted to a significant negative correlation (p<0.05) (Fig. 6b). Except that, tintinnid species richness exhibited a strong significant negative correlation with salinity in both the NFZ and SFZ (p<0.01), which was inconsistent with that in the NTZ, where it changed into a strong significant positive correlation (p<0.01) (Fig. 6b).

4 Discussion

In a nutshell, this study presents a first holistic epitome of microzooplanktonic ciliate community divergences and corresponding biotic–abiotic interplay among five temperature zones (NFZ, SAZ, NTZ, TZ, SFZ) spanning the global scale, revealing significant divergence in trait-based assemblages driven by temperature zone-specific physicochemical conditions. Simultaneously, it is noteworthy that our data-driven multivariate analyses demonstrated pronounced heterogeneity in ciliate trait structures (including vertical distribution patterns, latitudinal dynamics, size spectrum and biodiversity metrics) among five temperature zones (Figs. 2–4). Among these, abiotic parameters, particularly temperature, likely played a significant role in driving these variations, as hypothesized (Chapin et al., 1997; Anderson et al., 2021; Tanioka et al., 2022; Jiao et al., 2024). However, the current dataset remains geographically constrained, particularly lacking representation from Atlantic Ocean ecosystems where ciliate communities may exhibit distinct adaptive strategies. Hence, future research should prioritize comparative studies in Atlantic systems to test the global applicability of these findings. Additionally, more emphasis should be put on uncovering the trophy mode of pelagic ciliates in marine ecosystems.

4.1 Significant divergences in functional trait of ciliate size spectrum

Plankton size spectrum, which represents the distribution of individuals within a community or ecosystem by numerical abundance or biomass across size classes typically displayed on log axes, plays a crucial role in modulating various microbial processes, such as the carbon cycle driven by prey–predator interactions (García-Comas et al., 2016; Andersen, 2019; Trombetta et al., 2020; Serra-Pompei et al., 2022; Antoni et al., 2024; Atkinson et al., 2024). Simultaneously, the size spectrum provides insights into the ecological functions within marine food webs (Vandromme et al., 2012). In this sense, although empirical evidence has elucidated both the functional traits of plankton size spectra and valuable concurrent models, the majority of integrative analyses have primarily focused on biomass density within the size spectrum rather than on the abundance distribution across different trophic levels (Sprules et al., 2016; Blanchard et al., 2017; Atkinson et al., 2024; Stukel et al., 2024). Currently, research on specific zooplankton assemblages, such as microzooplanktonic ciliates (Wang et al., 2024b), is rarely studied on a global scale. Similar to Stukel et al. (2024), our study revealed that the slopes of abundance size spectra in both the NFZ and SFZ were steeper in polar seas than the other three regions latitudinally (Fig. 3). Furthermore, the consistently steeper slopes at the surface compared to the 200 m layer across all regions (Fig. 3) suggest the following: (1) a depth-dependent shift in pelagic ciliate community size structure and (2) greater accessibility of prey for meso-/macro-zooplankton in surface waters compared to the 200 m layer, thereby influencing carbon flux efficiency to higher trophic levels (Stukel et al., 2024).

In addition, Stukel et al. (2024) depicted that the slopes of the normalized biomass size spectra varied from −1.6 to −1.2 (median slope was −1.4), spanning over 5 orders of magnitude from phytoplankton to macrozooplankton in plankton communities in the tropical and subtropical seas. In contrast, our findings revealed the median slope was about −0.53 for the biomass size spectrum (no clear straight line on a log–log plot) across all discrete depths of the global seas (Fig. 4). We deem that the finer-scale monospecific trophic group, spanning 1 order of magnitude (10–200 µm, microzooplankton), might be too small to accurately calculate the slopes of the normalized biomass size spectra (Sheldon et al., 1972). Conversely, it is noteworthy that the slopes of the abundance size spectrum exhibited an inverse relationship between abundance and body size (Fig. 3), resembling the pyramid-of-numbers concept (Elton, 1927; Trebilco et al., 2013; Blanchard et al., 2017). Hence, we posit that the slope of the abundance size spectrum may be more informative than its biomass counterpart in covering 1 order of magnitude within the plankton community.

4.2 Tintinnid biodiversity dynamics and its underlying formation mechanisms

By virtue of its critical role in regulating ecosystem processes and resource utilization efficiency, plankton species diversity plays a crucial role in marine ecosystem functioning and biogeochemical cycling (Chapin et al., 1997). Similarly, a higher functionally similar species diversity enhances stability in resistance and resilience aspects of marine ecosystem processes (Ibarbalz et al., 2019; Benedetti et al., 2021; Chust et al., 2024). Consistent with both observational and modelling studies, tintinnid biodiversity was highest in the tropical and subtropical seas and was lowest in the polar seas (Fig. 5) (e.g. Sherr et al., 1997; Dolan et al., 2014, 2016; Righetti et al., 2019; Benedetti et al., 2021; Wang et al., 2020, 2024b; Li et al., 2016, 2018, 2022). Two explanations may account for this phenomenon. On the one hand, the intrinsic mechanism is the endosymbiosis (Kutschera and Niklas, 2005). After a long-term genetic DNA exchange and evolution process driven by closely prey–predation interaction (Chen et al., 2012), more diversified phytoplankton in tropical zone (Tian et al., 2024) is probably responsible for subsequent higher tintinnid biodiversity compared to polar zones through endosymbiosis mechanism (Margulis and Sagan, 2002; Clark et al., 2023).

On the other hand, physical barriers constitute a fundamental extrinsic mechanism governing plankton biogeography (Amargant-Arumí et al., 2024; Antoni et al., 2024; Chust et al., 2024). Generally, large-scale hydrographic features, particularly oceanic gyres and distinct water masses, create biogeographic discontinuities that disrupt ecological connectivity despite physical ocean connectivity (Yang et al., 2020). These mesoscale structures establish unique ecoregions with characteristic environmental sensitivities (Longhurst, 2007), as evidenced by pronounced tintinnid community differentiation across the North Pacific Gyre, Subarctic Gyre and Beaufort Gyre systems (Wang et al., 2020). Therein, our results revealed that tintinnid biodiversity was highest in the tropical (West Pacific and Indian Ocean) and temperate (North Pacific) seas, followed by the Sub-Arctic (Bering Sea) and polar seas (Arctic Ocean and Southern Ocean around Antarctic) (Fig. 5) consistent with Wang et al. (2020), proving that plankton biogeography was deeply affected by oceanic gyres. Ultimately, elucidating biodiversity patterns across diverse temperature zones provides critical insights into microzooplankton adaptive affinity potential under climate change scenarios, particularly regarding niche conservation versus ecological plasticity in response to shifting oceanographic boundaries.

4.3 Physicochemical factors determine the habitat of microzooplankton

Hydrography habitat conditions formed by large gyres (horizontal) or water masses (vertical) are critical factors in reshuffling sophisticated species composition of the microbial food web (Lennartz et al., 2024). Conventionally, temperature can impact plankton biodiversity through regulating intrinsic temperature-dependent metabolic processes, which further determine which kind of species can live in a specific temperature environment (Archibald et al., 2022; Lukić et al., 2022; Weisse, 2024). Coincidentally, the statistically positive correlation observed between tintinnid species richness and temperature (Fig. 6) fully supports the abovementioned ecological process. In this perspective, we conclude that temperature determines organism mortality by affecting their thermal affinity within biogeochemical cycles (Knies et al., 2009; Stuart-Smith et al., 2015; Archibald et al., 2022; Chust et al., 2024) through an indirect effect (Weisse and Sonntag, 2016; Weisse, 2024). Similarly, through modulating osmotic pressure, salinity plays a crucial role in shaping the species composition of the microbial food web (Pedrós-Alió et al., 2000; Zang et al., 2024) and in hindering the dispersal of Pacific species into the Arctic Ocean (Wang et al., 2019b, 2022b). Our study, along with others, indicates that ciliate inhabiting higher salinity environments in both the TZ and NTZ (Fig. S8) compared to polar regions might be a reflection of their higher osmotic pressure affinity.

Furthermore, the chl a roughly represents the phytoplankton at a specific sampling layer, which further influences marine ecosystem stability through both quantitative (abundance) and qualitative (nutrient composition) pathways via the fundamental prey–predator interplay (Šolić et al., 2010; Våge and Thingstad, 2015; Holm et al., 2022). As direct micro-grazers of phytoplankton, both the abundance and species richness of ciliates exhibit a significant positive correlation with chl a (Figs. 6 and S8–S10), aligning with the aforementioned viewpoint regarding the ecological role of chl a (Li et al., 2024). As outlined above, coupling with our results about multivariate analyses revealed strong hydrographic–ciliate relationships (Fig. 6), while the observed trait plasticity in ciliate communities (Yu et al., 2022) further supports the predominance of bottom-up control mechanisms (resource availability, prey quality) (Lu and Weisse, 2022; Wang et al., 2023c, 2024c) over top-down regulation (predation pressure from microcrustaceans) (Power, 1992; Calbet et al., 2001; Worm and Myers, 2003) in structuring global pelagic ciliate communities.

4.4 Prediction for microzooplanktonic ciliate community to future global warming

Global warming, primarily stemming from anthropogenic CO₂ emissions, has caused enduring and irreversible impacts on marine ecosystems globally, impelling a suite of threats to biodiversity and marine ecosystem, such as phenology evolution and adaptation, species poleward dispersal and body-size miniaturization (Daufresne et al., 2009; Poloczanska et al., 2013; Atkinson et al., 2015; Hastings et al., 2020; Møller and Nielsen, 2020; Yasumiishi et al., 2020; Wang and Wu, 2022; Qian et al., 2023; Wang et al., 2024a). To date, contemporary biogeographic observations have revealed marked increases in planktonic abundance and biodiversity across polar and subpolar seas (Ershova et al., 2015; Wassmann et al., 2015; Hunt et al., 2016; Kim et al., 2020; Lewis et al., 2020; Mueter et al., 2021; Wang et al., 2022a, 2023b), reflecting rapid thermal niche expansion under current warming regimes. Nevertheless, it should be mentioned that future global warming is expected to induce species extirpations by both compelling species beyond their thermal limits (Benedetti et al., 2021) and disrupting optimal survival habitats (Wang et al., 2024a).

Unfortunately, surface-dwelling ciliates (Kršinić, 1982; Wang et al., 2019a, 2023b, 2024a) are particularly vulnerable to recent, more frequent extreme temperature events, especially in tropical seas. Benedetti et al. (2021) projected a median speed of approximately 35 km per decade for the poleward shift of species dispersal under a high CO₂ emission scenario by the end of this century. In this perspective, our study provides a fundamental benchmark for understanding the adaptive strategies (extirpation, dispersal or adaptation) of ciliate to rapid warming processes in global seas. Meanwhile, unlike “winner” pioneer species possessing strong adaptation abilities (Casoli et al., 2020; Boutin et al., 2023), native species characterized by lower adaptive ability, such as the Arctic endemic tintinnid species Ptychocylis urnula, may either migrate passively to new environments (Wang et al., 2022a, 2023b, 2024a) or collapsed by a combination of warming and competition (Chust et al., 2024). Moreover, combined with our results that only the NFZ and SAZ exhibited an increasing trend (Δ_I≥0.03) in abundance–temperature correlation at surface layers compared with other three zones (Fig. S9), we predict that the pelagic surface–dweller ciliates in both the sub-Arctic and Arctic seas will benefit from the future global warming. Furthermore, the dynamics of the pelagic ciliate community in future trophic food webs and biogeochemical flux in the global marine ecosystem will heavily rely on how indigenous and/or intrusive species adjust to a warmer ocean state amidst multiple ecosystem stressors.

5 Conclusions

Our results provide comprehensive disparities in microzooplanktonic ciliate trait structure focused on size spectrum, biodiversity and biotic–abiotic interplay based on 1117 water samples from 175 stations across five temperature zones from the North Pole to the Southern Ocean (Antarctic). Concerning ciliate size spectrum, the slope of the normalized abundance value displayed an inverse relationship between ciliate abundance and body size, resembling a pyramid norm, while the biomass-size spectrum showed relatively smoother slopes. Additionally, tintinnid biodiversity was highest in tropical and subtropical seas and lowest in polar seas, likely influenced by endosymbiosis (intrinsic mechanism) and physical barriers (extrinsic mechanism). Furthermore, the interplay between biotic and abiotic factors manifested in temperature exerting a primary influence on the ciliate community structure. Under the current foreseeable rapid global warming process, we conjecture that bottom-up control (resource limitation) will play a more primary role through an indirect way in the global marine ecosystem.