Marine ciliates (phylum Ciliophora) are heterotrophic and mixotrophic protists that are of high biodiversity and playing important ecological roles in coastal and marine ecosystems (Song et al., 2003, 2009; Lynn, 2008). In pelagic food webs, marine ciliates usually act as microzooplankton that transfers carbon and energy to higher trophic levels (e.g., metazoans) by predating on bacteria, pico- and nano-sized protists (Stoecker and Capuzzo, 1990). Many ciliate species are also abundant in marine benthos, especially the water-sediment interface layer (Fenchel, 1969; Gao et al., 2017; Su et al., 2018), due to high adaptation to low oxygen contents, or even hypoxic environments (Park et al., 2008).

Studies of diversity and distribution of marine benthic ciliates were mainly based on microscopical identification and enumeration. The species richness and abundance of benthic ciliates were investigated using cell staining techniques for different sites along coastal line of the Black Sea (Azovsky and Mazei, 2005) and Yellow Sea (Meng et al., 2012; Zhou and Xu, 2016). While the application of these traditional methodologies has contributed greatly in ciliate diversity and ecology, a large part of ciliate diversity might have been missing. For example, fragile and nano-sized ciliates sometimes are abundant in marine habitats (Dong et al., 2014), but could be easily ignored due to loss of morphological integrity, cell disruption, and difficulty in telling from sediment particles (Caron and Countway, 2009). In the past decade, molecular approaches, especially high-throughput sequencing, have circumscribed these methodological difficulties and been widely applied in studying benthic protist diversity and biogeography (Doherty et al., 2010; Massana et al., 2015; Zhao and Xu, 2016). Recently, using molecular method based on high-throughput sequencing, Gong et al. (2015) have revealed that depth shapes the diversity and community composition and structure of microbial eukaryotes in the surficial sediments of the Bohai Sea and Yellow Sea of China, and ciliates appear to be important members in the benthic microeukaryotic communities. However, the variations in ciliate genetic diversity and community in these regions have not been analyzed, and it remains to be explored what factors determined the community assembly of coastal benthic ciliates across the Bohai Sea and Yellow Sea.

In our study, we focus on ciliates by extracting their sequences from the existing 454 sequencing dataset for benthic microeukaryotes in the Bohai Sea and Yellow Sea (Gong et al., 2015). We are particularly interested in the following questions: (1) What ciliate phyla or groups exist in these basins, and how the ciliate community was assembled and varied across space and time? (2) Is the biogeographic pattern of benthic ciliates different from the protists or microeukaryotes?

Previously, the benthic ciliate diversity and community composition and structure in the Yellow Sea (including NYS and SYS) were studied using morphology-based methods (MBM), e.g., Ludox density centrifugation and quantitative protargol stain (Meng et al., 2012; Zhou and Xu, 2016). The present study expands the sampling regions (including the Bohai Sea) and revisits the topic from a molecular perspective, revealing several consistent characteristics in comparison with the previous work. For example, the regional pattern of OTU-based ciliate species richness, much like these described in the MBM studies, exhibits a decreasing tendency from NYS to SYS (Fig. 3); the benthic ciliate richness was positively correlated with bottom water temperature, but negatively with bottom water salinity and water depth (Table 3). The depth-dependent pattern of alpha diversity of ciliates is much like that previously observed for microeukaryotes (Gong et al., 2015), in which limitation of sunlight and primary productivity may be also responsible for the reduction of ciliate species number in the benthos.

There are also notable discrepant results between these two types of investigations. Our study using MiSeq approach obtained a higher species-level richness (~480 OTUs) than the MBM studies (200–250 morphospecies), indicating the possibility that some benthic ciliate species may escape from the fixation and microscopical observations. Furthermore, Spirotrichean ciliates are known to have highly diverse morphological characters and ecological adaptations, and many new marine ciliate species estimated in recent years come from this class (Gao et al., 2016, 2017). Our molecular survey of the Bohai Sea and Yellow Sea showed that over 40% OTUs were unclassified within the class Spirotrichea (Fig. 2), indicating that the species discovery of spirotricheans (including the highly diverse choreotrichians) in marine habitats, especially from benthos, is far from being ended, and that potentially many new taxa or poorly known morphospecies are needing to be morphologically described and linked to molecular information. In fact, the hidden diversity of ciliates is not just restricted within a few classes, and there were about 16% ciliate OTUs unclassified at the class level (Fig. 2). BLASTing these unclassified OTUs against GenBank showed that 215 OTUs (about 45% of all OTUs detected in this study) shared a similarity not greater than 97% with those of the morphologically identified ciliate species (Fig. S2, Tables S2 and S3), highlighting the potential for discovering new ciliate taxa in the coastal sediments. However, high levels (up to 1% difference) of intragenomic polymorphisms of rRNA genes in ciliates were reported (Gong et al., 2013). Selecting a proper OTU cutoff (e.g., 97% or lower) could reduce this risk of species diversity overestimation. In the future, much effect has yet to be made, for example, to cultivate and isolate the morphospecies in marine benthos, and to obtain longer sequences for resolving their systematic ranks.

It is well known that the proportion of 18S rRNA gene sequences of a taxon does not necessarily correlated to their relative abundance in cells, because there are great copy number variations within and between protist species (Gong et al., 2013). Our recent single-cell analysis indicate that the sequence-based relative proportion of a taxon with protistan communities more reflect their relative biomass rather than cell abundance (Fu and Gong, 2017), providing a reasonable basis for bridging these cell-based and rDNA-based characteristics. In comparison with the morphological surveys of benthic ciliate communities in the Yellow Sea, our molecular investigation showed some marked differences in the ciliate community composition and structure. The MBM studies generally found the classes Prostomatea and Karyorelictea dominated in biomass (Meng et al., 2012; Zhou and Xu, 2016), which are contrast with the present study indicating the dominance of Spirotrichea sequences (biomass) (Fig. 2 and Table 2). Approximately a quarter of sequences in this study were affiliated with the two spirotrichean subclass Choreotrichia and Oligotrichia, which are typical planktonic groups (Santoferrara et al., 2016; Zhang et al., 2017), and were hardly observed in the MBM studies involving fixation and morphological examination (Meng et al., 2012). It is possible that these planktonic cells have been broken during physical mixing with sediment particles and extraction using Ludox density centrifugation, thus lost the integrated morphology for proper identification under microscopes, because these protists generally have a relatively thinner cortile membrane and more fragile than most other ciliates (Song et al., 2003, 2009). Previous studies have also demonstrated microscopic observations may miss a large number of the existent species (e.g., Tamura et al., 2011). Alternatively, these ciliates presented as cyst forms in the surface sediments (e.g., Chao et al., 2013), thus could not be readily identified and enumerated in previous morphology-based surveys. However, every coin has two sides. The rDNA-based surveys could not tell which community members are viable or dormant, thus have to be interpreted as a reflection of recent situation of the benthic community (Zhu et al., 2018). It is also likely that the eukaryotic universal primers used in this study are biasing against prostomatean ciliates, thus fail to amplify their 18S rDNA sequences. Another cause for the absence of prostomatean rDNA sequences may be that they have been PCR amplified but not properly classified due to the lack of enough reference sequences in the publically available databases. By building phylogenetic trees containing reference and sequences yielded in this study (data not shown), and BLASTing against the GenBank (Tables S2 and S3), we did not obtain any evidence for the notion that some prostomatean sequences might be lumped into the category of unclassified ciliates. This possibility can thus be ruled out. To some extents, the MBM may have an advantage over rDNA-based approaches in detecting the most viable members, which are assumed to be directly involved in the ecological and biogeochemical processes in the microbial food webs at the sampling time points.

Our analyses point to a distinct pattern of regional variations in ciliate community structure among the BHS, NYS and SYS (Table 1), which was not previously recognized. In the MBM studies, the benthic ciliate community structure was not significantly different between the NYS and SYS (Meng et al., 2012; Zhou and Xu, 2016), indicating that the inactive members may be quite differently assembled among these regions. This implies that classical taxonomic work has to expand the geographic range of sampling sites from coastlines to off-shore locations for species discovery. In addition, the present study provides new data on benthic ciliate diversity in the BHS, and shows there is much higher species richness of ciliates in the benthos of this basin than in the NYS and SYS, illustrating that the shallow basin could be a hotspot for benthic ciliate diversity.

It is interesting to note that the relative proportion of Litostomatea sequences was significantly higher in the deeper SYS than in the shallow NYS and BHS, a statistically supported characteristic for the community structure changes among the three basins (Table 2). The litostomateans are carnivorous ciliates primarily feeding on other protozoans even small-sized metazoans. The increasing biomass contribution of these carnivorous in deeper sites reflects a decrease of bacteriophages and herbivores, which are sustained by bacterial and microalgal food (microphytobenthos) dwelling in or depositing on the surface sediments. With the increasing water depth, the declining temperature and light will increasingly limit the growth of microphytobenthos and hence bacterial and algal predators, resulting in the depth-dependent pattern of litostomateans observed in this study. However, again this pattern disagrees with the MBM studies in the similar regions, in which the carnivorous biomass seemed to be decreasing in the SYS (Meng et al., 2012; Zhou and Xu, 2016), suggesting that the changes in relative abundance/biomass of a given taxon in the communities does not necessarily follow its biological abundance in the samples, as we previously discussed (Zhu et al., 2018).

The spatial scale of the present study is relatively large (spanning about 500 km), it is thus important to consider geographic distance in addressing the relative importance of spatial factors (horizontal and vertical distances that are related to biological dispersal) and environmental filtering in shaping the benthic ciliate communities in the BHS, NYS and SYS (Martiny et al., 2006; Gong et al., 2015). Our study shows that water depth is the most important factor in benthic ciliate community assembly, a result similar to the biogeographic study of microeukaryotes (Gong et al., 2015), but different from that of fungi, which showed distinct seasonality (hence environmental filtering is the most important driver) in the similar studied area (Wang et al., 2017). This demonstrates that different microeukaryotic groups may have different biogeographies. Life history, trophic style, growth kinetics, and ecological association may be key factors explaining these differences, which is an interesting issue warranted further investigations.