The photic‐aphotic divide is a strong ecological and evolutionary force determining the distribution of ciliates (Alveolata, Ciliophora) in the ocean

The bulk of knowledge on marine ciliates is from shallow and/or sunlit waters. We studied ciliate diversity and distribution across epi‐ and mesopelagic oceanic waters, using DNA metabarcoding and phylogeny‐based metrics. We analyzed sequences of the 18S rRNA gene (V4 region) from 369 samples collected at 12 depths (0–1000 m) at the Bermuda Atlantic Time‐series Study site of the Sargasso Sea (North Atlantic) monthly for 3 years. The comprehensive depth and temporal resolutions analyzed led to three main findings. First, there was a gradual but significant decrease in alpha‐diversity (based on Faith's phylogenetic diversity index) from surface to 1000‐m waters. Second, multivariate analyses of beta‐diversity (based on UniFrac distances) indicate that ciliate assemblages change significantly from photic to aphotic waters, with a switch from Oligotrichea to Oligohymenophorea prevalence. Third, phylogenetic placement of sequence variants and clade‐level correlations (EPA‐ng and GAPPA algorithms) show Oligotrichea, Litostomatea, Prostomatea, and Phyllopharyngea as anti‐correlated with depth, while Oligohymenophorea (especially Apostomatia) have a direct relationship with depth. Two enigmatic environmental clades include either prevalent variants widely distributed in aphotic layers (the Oligohymenophorea OLIGO5) or subclades differentially distributed in photic versus aphotic waters (the Discotrichidae NASSO1). These results settle contradictory relationships between ciliate alpha‐diversity and depth reported before, suggest functional changes in ciliate assemblages from photic to aphotic waters (with the prevalence of algivory and mixotrophy vs. omnivory and parasitism, respectively), and indicate that contemporary taxon distributions in the vertical profile have been strongly influenced by evolutionary processes. Integration of DNA sequences with organismal data (microscopy, functional experiments) and development of databases that link these sources of information remain as major tasks to better understand ciliate diversity, ecological roles, and evolution in the ocean.


I N T RODUC T ION
CILIOPHORA is one of the most diverse and ubiquitous phyla of unicellular eukaryotes. This clade originated at least one billion years ago, and some of the lineages that today prevail in the ocean appear to have radiated during the last 570 million years (Dunthorn et al., 2015;Fernandes & Schrago, 2019). Ciliates are morphologically diverse and present various lifestyles in marine waters, including free-living planktonic forms as well as parasites and commensals of invertebrates and fishes (Lynn, 2017). Ciliates form an important part of zooand mixoplankton in the nano-and micro-size fractions, with cell dimensions generally ranging between 10 and 100 μm, and nutrition strategies varying from exclusive phagotrophy of smaller protists and prokaryotes to combinations of phagotrophy and phototrophy via photosynthetic endosymbionts or kleptochloroplasts (McManus & Santoferrara, 2012;Weisse, 2017).
Planktonic ciliates act as trophic links between primary producers and larger consumers, contribute to photosynthesis, and participate in nutrient recycling (Calbet & Saiz, 2005;Pierce & Turner, 1992), but differences in taxon-specific functions remain vastly understudied. Despite the morphological and functional diversity of marine ciliates, their taxonomic composition and ecological roles are usually generalized from what is known about plankton in shallow and/or sunlit waters (Pierce & Turner, 1992). In these waters, ciliate assemblages are usually prevailed by algivorous and mixoplanktonic species of Oligotrichea, which are included in only one out of the 11 acknowledged major groups or "classes" of Ciliophora (Adl et al., 2019). Much less is known about ciliate taxa and their potential activities (bacterivory, detritivory, parasitism) in deeper waters, where sharp changes in major environmental factors (e.g. light, temperature, oxygen, food resources) are known to affect the diversity and composition of broader protist communities (Blanco-Bercial et al., 2022;Ollison et al., 2021;Schnetzer et al., 2011;Sogawa et al., 2022). This causes a biased view of ciliate contributions to the functioning of marine ecosystems, where these and other protists greatly influence productivity and biogeochemistry (Worden et al., 2015).
Microscope analyses of samples collected across vertical profiles have shown changes in ciliate assemblages with depth. Reports have shown decreases of one to three orders of magnitude from above to below the chlorophyll maximum in total ciliate abundance (e.g. Gomez, 2007;Pitta et al., 2000;Santoferrara et al., 2011;Tanaka and Rassoulzadegan, 2002;Wickham et al., 2011) and changes in species distribution based on the morphologically distinct tintinnids (Alder, 1999;Dolan et al., 2019;Wang et al., 2022). While light microscopy remains key in cell quantification, these analyses do not allow for accurate taxonomic classification of most (less conspicuous) ciliates, unless complex staining methods are used. These technical limitations, along with vast marine areas that remain logistically challenging to explore (Agatha, 2011), limit our knowledge on the diversity and distribution of marine ciliates in the deep ocean.
DNA metabarcoding offers the potential to quantify diversity at high taxonomic resolution by comparison of query 18S rRNA gene sequences to reference databases. Many marine metabarcoding surveys have explored protist communities in recent years (Santoferrara et al., 2020), but the staggering diversity of these microorganisms renders limited detail on ciliates, which usually end up analyzed in bulk. Other studies have focused only on Oligotrichea members Santoferrara et al., 2016Santoferrara et al., , 2018. The few available metabarcoding studies focusing on the vertical distribution of the phylum Ciliophora in both photic and aphotic waters have revealed interesting patterns about variations in assemblage composition (beta-diversity) with depth. First, a higher heterogeneity in the distribution of ciliate taxa across 2000-m vertical profiles than along 1300 km has been reported in the western Pacific Ocean (Zhao et al., 2017). Second, depth has been found to explain assemblage composition variations better than either geographical distance or a combination of standard environmental variables (e.g. temperature, dissolved oxygen, chlorophyll a, bacterial abundance) in the South China Sea (Sun et al., 2019). Third, ciliate assemblages show distinct taxonomic compositions across epi-, meso-, and bathypelagic waters of the Pacific, Atlantic, and Indian oceans, with the typical Oligotrichea prevalence shifting to increased contributions of Oligohymenophorea and enigmatic environmental clades (e.g. OLIGO5, NASSO1 sensu Boscaro et al., 2018) as depth increases (Canals et al., 2020). While patterns of beta-diversity with depth are thus evident, the vertical trends of alpha-diversity (e.g. richness) remain unclear given contradictory data in the literature (Canals et al., 2020;Sun et al., 2019;Zhao et al., 2017;see Discussion). Also, contemporary patterns of ciliate diversity and taxonomic composition in the vertical profile have not been reported in a phylogenetic context, even though depth can explain both community assembly and the evolutionary histories of other protists in the ocean (Junger et al., 2022;Weiner et al., 2012).
Here, we study the diversity, taxonomic composition, and phylogenetic partitioning of marine ciliates in photic and aphotic waters by using DNA metabarcoding at the most detailed vertical (12 depths across 1000 m) and temporal (monthly over 3 years) scales so far used for this goal. Taking advantage of a protist dataset obtained at one of the longest oceanic time-series, the Bermuda Atlantic Time-series Study (BATS) in the Sargasso Sea (Blanco-Bercial et al., 2022), and using phylogeny-based approaches, we tested the following hypotheses: (1) ciliate alpha-diversity decreases significantly with depth; (2) beta-diversity shifts significantly from photic to aphotic waters; (3) ciliate lineage distributions are partitioned in the vertical profile, possibly reflecting different trophic strategies and evolutionary histories across depth.

M AT ER I A L S A N D M ET HOD S Environmental data and DNA metabarcoding
The field and laboratory procedures for environmental characterization and DNA metabarcoding have been detailed by Blanco-Bercial et al. (2022). In brief, vertical profiles of water temperature, salinity, oxygen concentration, and chlorophyll fluorescence were measured with a CTD from the surface down to a 1000-m depth at the BATS station (31°40′N, 64°10′W) monthly from February 2016 to December 2018. On each date, water samples were obtained at 12 depths (0,40,80,120,160,200,250,300,500,600,800, and 1000 m) with Niskin bottles mounted on a rosette. For metabarcoding, 4 L of seawater were immediately filtered on 0.22 μm Sterivex filters (MERCK). Details on DNA extraction, PCR amplification of the hypervariable V4 region of the 18S rRNA gene (using primers by Stoeck et al., 2010), library preparation, and Illumina MiSeq sequencing using 2x250 PE V2 chemistry have been published (Blanco-Bercial et al., 2022). Raw reads for a total of 369 samples are available in the Sequence Read Archive (PRJNA769790).
About 170,000 reads distributed among 8637 ciliate variants were extracted for further analysis in this paper. The bioinformatic tools and code used in this study are available at https://github.com/LuSan tof/Cilio phora -Sarga sso-metab arcod ing.git. Identification was done with a naïve Bayes classifier (Bokulich et al., 2018) against the PR2 database v4.12.0 (Guillou et al., 2013) trimmed to contain only the V4 region. Taxonomic assignment was complemented with BLAST+ v2.12.0 (Camacho et al., 2009) and phylogenetic trees (built in the same way as described below in "Other phylogenetic analyses"), using the EukRef-Ciliophora curated dataset (Boscaro et al., 2018) and a curated alignment for Oligotrichea (Ganser et al., 2022). The complementation of three reference datasets was needed due to uneven curation of ciliate sequences in the PR2 version available at the time of analysis; all three references have now been reconciled and are available in PR2 v5.0.0 (Santoferrara and Ganser, in Vaulot, 2023). The three taxonomic assignment methods generally agreed, and discrepancies were solved by retaining the most conservative identification (i.e. up to the broadest taxonomic rank). Taxonomic classification follows the rank-less scheme of Adl et al. (2019), although we mention common "classes," "subclasses," and "orders" for practicality. Uncharacterized environmental clades (known only from the retrieval, curation, and phylogenetic analysis of environmental 18S rRNA gene sequences available in GenBank) are from Boscaro et al. (2018).
Ciliate variants were postclustered with LULU (Frøslev et al., 2017). This method groups variants based on similarity (BLAST results with 97% coverage and 97% identity as a cut-off) and co-occurrence (95% of samples or more), under the premise that rarer co-occurring variants represent intragenomic variability or technical errors derived from more abundant variants. Most "daughter" variants had less than 10 reads and their taxonomic affiliations were redundant with "parent" variants (i.e. no relevant distributional or taxonomic information was lost due to postclustering). Postclustered variants with less than five reads were eliminated, resulting in a final curated dataset of 2053 variants (almost 162,500 reads). The curated variants, taxonomy, and read abundances are provided in Table S1.

Diversity analyses
Multivariate and statistical analyses were done in QIIME2 v2019.7 (Bolyen et al., 2019) and RStudio using the vegan package (Oksanen et al., 2019). As described below, certain analyses included continuous or categorical environmental variables. Continuous variables (sample depth, temperature, salinity, oxygen concentration, and CTD fluorescence) were standardized to scale each variable to a zero mean and unit variance. Categorical variables include layer (11 layers: L0 to L2 and L3 to L10 corresponding to photic and aphotic waters, respectively) and season (mixed, spring, stratified, and fall), as assigned by Blanco-Bercial et al. (2022). The environmental metadata is provided in Table S2.
Ciliate sequence data were rarefied to a sequencing depth of 300 reads per sample (220 samples). Alphadiversity was computed as the Faith's phylogenetic diversity index (Faith, 1992), which considers variant richness and their phylogenetic relatedness (the sum of the lengths of variant branches on a tree built with MAFFT-aligned variants and FastTree under default parameters as implemented in QIIME2; Bolyen et al., 2019). Statistical differences in the Faith's index among layers and seasons were tested with the Kruskal-Wallis analysis of variance. Beta-diversity was evaluated with principal coordinates analyses computed with phylogenetically-aware qualitative (unweighted UniFrac) and semi-quantitative (weighted UniFrac) distance matrices (Lozupone & Knight, 2005). Because both distance matrices yielded very similar results, only results based on weighted UniFrac distance are shown. The significance of sample groups was evaluated with analysis of similarities (ANOSIM) using 999 permutations (Clarke, 1993). Pairwise correlations between selected variables were tested with Pearson's or Spearman's coefficients (R) with Bonferroni correction for multiple comparisons. The correlation between weighted UniFrac dissimilarity (sequence-based) and Euclidean distance (standardized continuous environmental variables) matrices was tested with the Mantel test using Spearman's coefficient and 999 permutations (Legendre & Legendre, 2012). The subsets of environmental variables that best correlated with sequencebased beta-diversity were identified with BIO-ENV (Clarke & Ainsworth, 1993).

Phylogenetic placement and correlations with depth
Phylogenetic placement and postplacement tools were used to test for correlation between the relative read abundance of major ciliate lineages (all classes, plus Oligohymenophorea subclasses and Oligotrichea orders) and depth. For phylogenetic placement, the Ciliophora reference alignment, reference tree, and tree model provided by Rajter and Dunthorn (2021) were used. Our variants and the reference sequences were combined and re-aligned with MAFFT v7 (FFT-NS-i algorithm; Katoh et al., 2019) and then analyzed with EPA-ng v0.3.8 (Barbera et al., 2019). The phylogenetic placement was then analyzed with GAPPA v0.7.1 (Czech et al., 2020). This tool was used to calculate and visualize correlations (Spearman's R) between each edge of the reference tree and depth (Czech et al., 2020;Czech & Stamatakis, 2019). Placements were visualized and edited in iTOL v6.6 (Letunic & Bork, 2021).

Other phylogenetic analyses
To study Discotrichidae sequences in more detail, a dataset was generated by combining our short V4 variants and longer or full 18S rRNA gene sequences. The dataset included: (1) Discotrichidae variants with 10 or more reads in our data, (2) all the Discotrichidae sequences available in PR2 (retrieved Nov 19, 2022 from https://app.pr2-datab ase.org/pr2-datab ase/), and (3) all the Nassophorea and select representative sequences from other Ciliophora classes extracted from a reference dataset (Rajter & Dunthorn, 2021). Sequences were aligned with MAFFT v7, using the G-INS-1 method (Katoh et al., 2019). Two versions were generated: trimmed with Gblocks v0.91b under default parameters (Castresana, 2000) and untrimmed. Gene trees were inferred with IQ-Tree (Trifinopoulos et al., 2016), using the substitution model GTR + F + I + G4 selected by Modelfinder under the Akaike information criterion (Kalyaanamoorthy et al., 2017). Branch support was calculated with 1000 iterations of the ultrafast bootstrap approximation (Minh et al., 2013). The tree generated with the untrimmed alignment better matched reported inferences (Rajter & Dunthorn, 2021), and it was thus kept for further analyses. The tree was visualized and edited in iTOL v6.6 (Letunic & Bork, 2021).

Alpha-and beta-diversity
The DNA metabarcoding data collected at the BIOS station monthly over 3 years yielded a total of 2053 ciliate variants. Variant richness decreased progressively from layer 0 to layer 10 (surface to 1000-m depth), and this decrease was globally significant based on a Kruskal-Wallis test ( Figure 1A). Variant richness also presented a significant, inverse correlation with sample depth ( Figure 1B). In turn, depth was negatively and significantly (p-value <0.00001) correlated with temperature (Pearson's R = −0.94), salinity (−0.94), oxygen (−0.84), and fluorescence (−0.46). There was no significant trend for variant richness with seasons ( Figure S1).
Beta-diversity analysis indicated a gradual change in the variant assemblage composition from the surface down to 1000 m, with a significant difference between samples obtained in the photic versus aphotic zones (ANOSIM; Figure 2A). This ordination of samples presented a significant correlation with the Euclidean distance matrix based on five continuous environmental variables (global Spearman's R = 0.44; p = 0.001). Of these variables, depth was identified by BIO-ENV as the single, most important factor influencing the observed beta-diversity patterns (Spearman's R = 0.49; p < 0.001).
In agreement with variant-level assemblage variations across layers (Figure 2A), taxonomic composition showed a shift in dominance from the class Spirotrichea (mostly the subclass Oligotrichea, including the orders Oligotrichida and Choreotrichida) in photic waters to the class Oligohymenophorea (mostly the subclass Apostomatia and the environmental clade OLIGO5) in aphotic waters ( Figure 2B). Other groups (other Spirotrichea, classes Litostomatea and Phyllopharyngea, the incertae sedis family Discotrichidae, and other Ciliophora) presented similar proportions in both photic and aphotic waters. The differences in variant distribution and taxonomic composition were even more obvious when the analyses were based on sequences grouped by layer ( Figure S2). In contrast, there were no clear trends in beta-diversity with seasons ( Figure S3).

Vertical distribution of prevalent variants
The progressive changes in the ciliate assemblage from surface down to 1000 m were also evident when considering the prevalent variants in each layer (Figure 3). Most prevailing variants in the upper three layers (L0 to L2) belonged to Oligotrichida, especially to the family Strombidiidae, followed by the family Tontoniidae and the environmental clades SPIOL1 and SPIOL3. Next, the Choreotrichida were represented by the aloricate genus Leegaardiella and representatives of the suborder Strobilidiina, and tintinnid variants identified to two genera and three species. Also present in the upper three layers were the incertae sedis genus Askenasia and variants within the environmental clades PHYLL4 (Phyllopharyngea) and OLIGO5 (Oligohymenophorea).
Layer 3 and below (L3 to L10) were prevailed mostly by Oligohymenophorea variants assigned to the subclass Apostomatia and the environmental clade OLIGO5. Among the prevalent variants in aphotic waters, there were some affiliated with the environmental clade NASSO1 (incertae sedis family Discotrichidae) and one assigned to the incertae sedis genus Protocruzia in layer L8. Oligotrichida had prevalent variants in all layers, except in layer L8.
Several variants prevailed in multiple layers, but always within either photic or aphotic waters.
Two variants (Var6548_Apostomatia and Var8745_ OLIGO5) were among the prevalent ones in all the aphotic layers examined. Layer 3 shared variants with both upper and lower layers. The pie chart (Figure 3) shows that all the top-15 variants could be identified to classes, except one only identified to the "superclass" CONTHREEP. About 40% of the variants could be identified at least to family, and only 6% could be identified to species.

Phylogenetic placement and taxa correlations with depth
The phylogenetic placement of variants was consistent with our taxonomic assignments and showed a strong correlation of several clades with depth ( Figure 4). Clades highly anti-correlated with depth (Spearman's R < −0.5) belonged to Oligotrichea (both Choreotrichida and Oligotrichida), followed by moderate anti-correlations (R ~ −0.5) in other Spirotrichea (some taxa in the subclass Hypotrichia), Litostomatea, Prostomatea, Phyllopharyngea, and the incertae sedis genus Cyclotrichium. Clades highly correlated with depth (R > 0.5) belonged to Oligohymenophorea, especially the subclass Apostomatia and some Astomatia (a subclass that does not form a monophyletic group in the reference tree and shows an opposite trend with depth for some branches); the subclasses Scuticociliatia, Peritrichia, and Peniculia show moderate correlations with depth (R ~ 0.5). Very weak or no correlation with depth was observed for some Hypotrichia, the incertae sedis genus Protocruzia, and the classes Colpodea, Nassophorea, and Heterotrichea.
The incertae sedis family Discotrichidae (represented almost fully by variants in the environmental clade NASSO1; Boscaro et al., 2018) showed a moderate correlation with depth based on GAPPA (R ~ 0.5; Figure 4), but a more detailed analysis shows that at least some of its subclades are differentially distributed across depths ( Figure 5). The combination of all the Discotrichidae sequences available in the PR2 database, Discotrichidae variants with 10 or more reads in our data, and representative reference sequences from other ciliate groups yielded a tree with most relationships showing >50% bootstrap support ( Figure 5) and that agree with current phylogenetic inferences (Boscaro et al., 2018;Rajter & Dunthorn, 2021). The Discotrichidae formed an independent clade, in agreement with previous reports (Fan et al., 2014). Most Discotrichidae sequences correspond to the environmental clade NASSO1 (Boscaro et al., 2018), which we recovered with 99% bootstrap support ( Figure 5). Our variants were distributed among most NASSO1 subclades, intertwined with sequences from photic and aphotic waters of the Atlantic, Pacific, Indian, and Arctic oceans. However, two of the subclades were almost exclusively restricted to either photic or aphotic waters. The photic subclade (75% bootstrap support) included 36 PR2 sequences from the Atlantic, Indian, and Pacific oceans (collected between 1 and 150-m depths) and 33 of our variants (with 91% of reads in photic waters). The aphotic subclade (100% bootstrap support) included 5 PR2 sequences from the Atlantic and Pacific oceans (collected between 500-and 2500-m depths) and 29 of our variants (with 93% of reads in aphotic waters). A few of our Discotrichidae variants found exclusively in aphotic waters formed independent branches within or outside NASSO1 (variants 10120 or 10154-5, respectively).

DI SC US SION
We studied ciliate diversity across epi-and mesopelagic oceanic waters in the Sargasso Sea. Our study site is located in the limits of the North Atlantic subtropical gyre, where nutrient availability, primary production, and microbial community composition are driven by a local seasonal cycle characterized by a stratified summertime photic zone and winter mixing that extends into the upper layers of the aphotic zone (Church et al., 2013). Seasonality also drives epipelagic changes in the diversity and structure of protist communities, which are prevailed by phyto-and mixoplanktonic taxa (especially within Dinophyceae, Pelagophyceae, Ochrophyta, and other Stramenopiles; Blanco-Bercial et al., 2022). Ciliates, however, showed no clear seasonal patterns in alpha-or beta-diversity in the water column or epipelagic layers ( Figures S2 and S3). Lack of seasonal patterns may be due to many factors (e.g. trophic differences as compared F I G U R E 4 Phylogenetic placement of variants and clade-specific correlations with depth (color-coded in the scale). Each branch represents a reference sequence (Rajter & Dunthorn, 2021). Colors are based on the number and relative read abundance of variants placed into each clade and their correlation with depth. Gray branches are reference sequences with no match in our variant dataset. Only taxa represented in our dataset are labeled, including classes, subclasses, orders, or incertae sedis (i.s.) as needed. Some clades with no variants assigned were collapsed.
to phytoplankton) and match reports in diverse areas of the Pacific for ciliates (Sun et al., 2019) and other protists (Kim et al., 2014;Ollison et al., 2021). In contrast, we found sharp changes in ciliate diversity and taxonomic composition with depth, the focus of the rest of the discussion.
Our analysis of V4 18S rRNA gene sequences from the surface down to 1000 m in the Sargasso Sea showed the presence of eight out of 11 ciliate classes (Adl et al., 2019). The prevalent classes in our data (Spirotrichea, Litostomatea, Oligohymenophorea, Prostomatea, Phyllopharyngea) are known to inhabit marine plankton, while groups we rarely detected (Colpodea, non-Discotrichidae Nassophorea, Heterotrichea) have been less commonly reported in marine waters (Dunthorn et al., 2014;Lynn, 2017). We also detected the incertae sedis Protocruzia (known for marine water; Jiang et al., 2017) and Discotrichidae (discussed in the last paragraph of the Discussion). Of the classes we did not detect, Armophorea and Plagiopylea inhabit anoxic environments and Karyorelictea are only common in estuarine and marine benthos (Lynn, 2017), although they are sometimes recovered in plankton metabarcoding surveys (Canals et al., 2020).

F I G U R E 5
Phylogenetic tree of Discotrichidae (dotted box) and other Ciliophora representatives. The environmental clade NASSO1 (yellow background) includes several subclades (collapsed), including two recovered in either photic (orange) or aphotic (blue) waters. The distribution of sequences from the PR2 database and variants from this study is shown.

Significant decrease in alpha-diversity with depth
Previous DNA metabarcoding studies presented contradictory conclusions about relationships between ciliate alpha-diversity and depth. The number of Oligotrichea operational taxonomic units (OTUs) was found to be higher at 850 m than in the photic zone in one station of the North Atlantic . Ciliophora OTU richness was found to peak in the Deep Chlorophyll Maximum (DCM) and 200-m layers, thus presenting a non-linear relationship with depth from the surface to 2000-m waters of the western Pacific Ocean (Zhao et al., 2017). Ciliophora OTU-and morphospecies-based richness values were reported as unchanged in epi-versus mesopelagic waters of the South China Sea (Sun et al., 2019). A sharp decrease in Ciliophora OTU richness and the Shannon index was observed in surface plus DCM versus meso-plus bathypelagic waters sampled during the Malaspina expedition (Canals et al., 2020).
Using the most detailed vertical (11 layers) and temporal (monthly over 3 years) scales so far used to investigate Ciliophora diversity across depth in the ocean, we found an inverse relationship between variant richness (Faith index) and depth (Figure 1). The gradual decrease in variant richness across depth layers ( Figure 1A) was slightly disrupted by moderate increases in certain layers (L1, L6, L8, L9), but surprisingly it was not affected by the DCM (layer 2; Blanco-Bercial et al., 2022). In fact, there was an overall linear, negative correlation between richness and depth ( Figure 1B). While inconsistent results among previous and current data may be due to differences in methodologies (e.g. the use of different alpha-diversity indices), local oceanography, and sampling scales, the comprehensive depth and temporal resolutions analyzed here support that ciliate variant richness decreases with depth. This pattern is likely explained by abiotic and biotic factors that correlate with depth; for example, in our study, depth was negatively related to temperature, salinity, oxygen, and fluorescence. This is in general agreement with changes in productivity and overall environmental conditions known to affect protist alpha-diversity in the vertical profile, especially from photic to aphotic waters (Blanco-Bercial et al., 2022;Schnetzer et al., 2011;Sogawa et al., 2022).

Significant differences in assemblage structure with depth
Our results support significant changes in Ciliophora assemblage structure with depth, in agreement with previous metabarcoding studies (Canals et al., 2020;Sun et al., 2019;Zhao et al., 2017). This is supported by several lines of evidence in our data: (1) the gradual change in variant assemblages from surface to the 1000-m depth (Figures 2A and S2A); (2) the significant shift in assemblage composition from photic to aphotic waters, which matched a switch in prevalence from Oligotrichea to Oligohymenophorea (Figures 2B and S2B); and (3) the changes in variant prevalence across layers (Figure 3). Our analyses identified depth as the single, most important factor influencing the observed beta-diversity patterns, likely because important parameters that influence ciliate ecological niches change with depth.
The Oligotrichea also include mixoplankton species, of which we detected Pseudotontonia simplicidens (prevalent in the uppermost layers 0 and 1, above the DCM; Figure 3). Cells of P. simplicidens show chlorophyll a fluorescence and contain plastid 23S rRNA gene sequences, which suggest acquired photosynthetic capabilities (Johnson & Beaudoin, 2019;Pitta et al., 2000). Mixotrophy is a favored strategy in the upper, nutrientdepleted waters in the Sargasso Sea (Blanco-Bercial et al., 2022); combination of phagotrophy and photosynthesis would allow mixoplakton to escape from the nutrient limitation conditions (Edwards, 2019). We found additional prevalent variants within the oligotrichid families Tontoniidae and Strombidiidae, which harbor most ciliate mixoplankton species (Mitra et al., 2023). However, mixotrophic capabilities are species-specific, and the limited linkage of known mixoplankton morphospecies with their 18S rRNA gene sequences prevents us from confirming additional identifications in our dataset. Another common mixoplankton ciliate, the incertae sedis Mesodinium rubrum, is not captured by primers used in V4 metabarcoding surveys (Stoeck et al., 2010), and thus we cannot establish its contribution in our samples.
Oligohymenophorea is one of the most diverse ciliate classes, and common representatives in plankton belong to the subclass Scuticociliatia (Pierce & Turner, 1992). However, this group was modestly represented across depths in our data, and instead, we found the prevalence of the subclass Apostomatia and the environmental clade OLIGO5 in aphotic waters (layers 3 to 10; Figures 2  and 3). Apostomatia are usually animal symbionts and were represented in our data by the Astomatophorida family Opalinopsidae, the Apostomatida families Pseudocolliniidae and Foettingeriidae (with the later including two Vampyrophrya variants), and other unidentified Apostomatia. Astomatophorida include known parasites of cephalopods, such as the Opalinopsidae genus Chromidina (Souidenne et al., 2016). Apostomatida include known parasites (some of them lethal) of copepods and krill, such as the genera Vampyrophyra (Hartley-Grimes & Bradbury, 2007), Collinia (Capriulo & Small, 1985;Gómez-Gutiérrez et al., 2006), and Pseudocollinia (Lynn et al., 2014). Some Apostomatida are known to spend most of their life cycle encysted on host crustacean exoskeletons (Bradbury & Trager, 1967). This suggests that the high proportion of Apostomatia found in our aphotic zone samples may correspond to active or inactive parasites that we collected as dispersal stages or as part of broken animal tissues. Our findings also suggest potentially important roles of these ciliates in recycling of carbon stored in animal biomass.
The other Oligohymenophorea prevailing in aphotic waters, OLIGO5, form a large clade of environmental sequences identified by the global meta-analysis of ciliate GenBank sequences done by Boscaro et al. (2018). No known morphospecies are included in this clade, which shows a weak phylogenetic relationship to divergent representatives of Scuticociliatia (Boscaro et al., 2018). Since the OLIGO5 clade was made evident by Boscaro et al. (2018), it arose as an abundant and widespread clade, being reported for example in Subtropical and Subantarctic mesopelagic waters of the Pacific (Gutiérrez-Rodríguez et al., 2022), mesopelagic waters of the Scotia Sea (Duret et al., 2020), coastal and open waters of the Gulf of Mexico (Snyder et al., 2021) and diverse stations of the Tara and Malaspina expeditions (Canals et al., 2020). It has been speculated that OLIGO5 may include small bacterivores, histophagous taxa, or parasitic forms, but this remains largely unknown (Snyder et al., 2021).
There were two unexpected findings in some aphotic layers. First, layer 8 corresponded to the oxygen minimum zone (OMZ, oxygen <160 μmol/kg, 600-850 m) and presented distinct communities when considering other protist groups (Blanco-Bercial et al., 2022). In contrast, ciliate assemblages were comparable to adjacent layers and included some of the variants found in all or at least some other aphotic layers (Figure 3). The only exceptions among dominant variants were the presence of Protocruzia, a taxon also reported in other OMZs around the world (Canals et al., 2020), and the absence of Oligotrichida (Figure 3). Second, layers 9-10 had three Oligotrichida among the dominant variants (one of them also found to prevail in layer 2), even though Oligotrichida had decreased contributions in the upper aphotic layers (Figure 3). One possibility is that these sequences represent inactive settling forms, in the form of detritus or dormant cysts (known for Oligotrichida in coastal waters; Kim et al., 2008). Alternatively, these sequences may represent active, vegetative cells. The latter is supported, for example, by morphology-based (quantitative protargol staining) observations of a single, unidentified Strombidium sp. among the top-3 most abundant ciliates both in both photic and aphotic waters across a 1000-m-depth profile (Sun et al., 2019). Although Oligotrichida are known mostly from coastal and sunlit waters, most of the diversity in this taxon is unknown from the morphological and functional points of view (Santoferrara et al., 2017), and it is thus possibly more flexible than currently believed. Rather than algivorous, many Oligotrichida are known as bacterivorous (Chen et al., 2020) and may actually be omnivores that can switch among food resources based on their availability (Weisse, 2017).

Phylogenetic partitioning of ciliate taxa across depths
Different from previous metabarcoding studies on Ciliophora vertical profiles, we added a phylogenetic dimension to our study by using phylogenetic-aware indices of alpha-and beta-diversity. Our analyses then consider both independent units (sequence variants) and their evolutionary relationships. Combined with variant postclustering (Frøslev et al., 2017), phylogenetically-aware indices attenuate biases caused by closely-related variants (many of which can represent errors or intraspecific variation; Santoferrara, 2019) by reducing inflation in alpha-diversity values and overdispersion in beta-diversity ordinations (Faith, 1992;Lozupone & Knight, 2005). We also added a phylogenetic component to the study of vertical distributions by placing variants in a reference tree and correlating these placements with depth (Barbera et al., 2019;Czech et al., 2020).
The contemporary differences in taxonomic composition observed across depth (Figures 2 and 3) likely have a deep evolutionary history (Figure 4). The two distinct clades that showed opposite prevalence in photic versus aphotic waters, Oligotrichea and Apostomatia, respectively, also show opposite correlations with depth ( Figure 4). This suggests that the photic and aphotic divide is a strong ecological and evolutionary force determining ciliate distribution. Also, the temporal and spatial persistence of the vertical divide (this and previous studies; Canals et al., 2020) suggests that distinct ciliate clades are differentially adapted to the conditions found in each zone. It could then be hypothesized that important niche parameters are availability of light and algal prey for Oligotrichea, or animal hosts for Apostomatia, but this remains largely unknown. For other planktonic protists, niche partitioning by depth has been shown as important in both contemporary taxa distributions and over evolutionary scales. For example, for radiolarians and foraminiferans, niche differences across depth generate persistent distribution patterns (Biard & Ohman, 2020;Boltovskoy, 1999) and are suggested as a cause for speciation (Ishitani et al., 2014;Weiner et al., 2012). Additional studies are needed to fully understand the ecoevolutionary processes that lead to distinct ciliate distributions in the vertical profile.
Beyond Oligotrichida and Oligohymenophorea, less prevalent ciliate groups had variable degrees of correlation with depth, including moderate negative (e.g. for Litostomatea, Prostomatea, Phyllopharyngea) to no (e.g. Colpodea, Nassophorea, Heterotrichea) correlations (Figure 4). The incertae sedis family Discotrichidae, however, presented a more complex relationship with depth ( Figure 5). Discotrichidae used to be placed in Nassophorea, but this class is now known as non-monophyletic (Fan et al., 2014;Zhang et al., 2022). Only two Discotrichidae morphospecies, which were isolated from marine coastal sediments (Lopezoterenia paratorpens and Discotricha papillifera; Fan et al., 2014), have been sequenced for the 18S rRNA gene. However, Boscaro et al. (2018) retrieved 284 unidentified environmental sequences from GenBank that form a wellsupported, monophyletic clade with L. paratorpens and D. papillifera. These sequences are now incorporated into the PR2 database as Discotrichidae and most of them (277 sequences) are included in the environmental clade NASSO1 (Boscaro et al., 2018). Combining the NASSO1 sequences in PR2 and our variant dataset, two subclades emerge that are present almost exclusively in photic or aphotic waters ( Figure 5). Thus, these NASSO1 subclades may partition their distribution by depth, although we know nothing about NASSO1 trophic preferences. The few described Discotrichidae species are not planktonic, but interstitial in sediments, where they consume algae via a characteristic cytopharyngeal basket (Fan et al., 2014;Foissner, 1997). Since NASSO1 was individualized by Boscaro et al. (2018) and then incorporated into the PR2 database, it has been reported with relatively high frequencies in coastal and open waters, many times in combination with OLIGO5 (Canals et al., 2020;Duret et al., 2020;Gutiérrez-Rodríguez et al., 2022;Snyder et al., 2021). NASSO1 ubiquity is also confirmed by a search in metaPR2 v2.0 (Vaulot et al., 2022), which indicates that this clade is present exclusively in marine environments, with 92 V4-based ASVs from 1523 samples (40% of all the marine samples in the database) collected in photic and aphotic waters of the five oceans and the Mediterranean Sea.
In addition to OLIGO5 and NASSO1, we also detected other environmental clades identified by Boscaro et al. (2018): SPIOL1 and SPIOL3 (Spirotrichea, subclass Oligotrichida) and PHYLL4 (Phyllopharyngea, order Cyrtophoria). These enigmatic clades, plus the facts that ciliate environmental sequences can rarely be identified into families or lower taxonomic ranks (e.g. only 40% of our variants; Figure 3) and usually present low similarity to known morphospecies (only 90%-93% similarity in average; Canals et al., 2020), highlight the insufficient knowledge we have about ciliates in the ocean, especially in aphotic waters. Integration of DNA sequences with organismal data (microscopy, functional experiments) and development of databases that link these sources of information remain as major tasks to better understand ciliate diversity, evolutionary relationships, and ecological roles in the ocean.