Diversity and Seasonality Dynamics of Ciliate Communities in Four Estuaries of Shenzhen, China (South China Sea)

: Ciliates are fundamental components of microzooplankton, with important ecological roles. However, ciliate communities are particularly difﬁcult to monitor using conventional morphological approaches. New molecular tools, such as DNA metabarcoding, can facilitate the study of these communities. This study used high-throughput sequencing to examine the diversity and seasonal dynamics of ciliate communities in four estuarine ecosystems in the South China Sea from June 2019 to March 2020. The ampliﬁcation of the V4 region of 18S rDNA using ciliate-speciﬁc primers identiﬁed a total of 1645 amplicon sequence variants (ASVs), corresponding to 13 ciliate classes, 97 families, and 157 genera. The dominant species across all four sampling stations were spirotrichs (including choreotrichs, oligotrichs, and stichotrichs), oligohymenophorean scuticociliates, litosto-mateans Didinium , and prostomateans Cryptocaryon . Signiﬁcant differences in ciliate diversity and community composition in the four stations were mainly due to differences in rare, rather than abundant, ASVs. Analysis of the ciliate communities and seasonal patterns in their composition revealed that variations in habitat and environmental conditions have a greater effect than seasonal changes on community composition.


Introduction
Estuarine ecosystems form an important interface between terrestrial rivers and open marine ecosystems. They undergo highly dynamic changes and have variable environments because of alternating inputs of freshwater and seawater. Estuaries contain higher concentrations of inorganic and organic nutrients compared with most aquatic ecosystems [1]. Organisms that live in estuaries must adapt to the dynamic environment, including variations in water chemistry (such as salinity), tides, and river emissions. Ciliates are unicellular organisms that belong to the Alveolata supergroup; they have been found in freshwater, marine, soil, and extreme environments [2,3]. They dominate many microeukaryote communities in terms of their diversity and abundance [4,5]. In aquatic ecosystems, ciliates display enormous functional diversity as producers, consumers, and decomposers of food webs [6]. Moreover, their diversity and community dynamics contribute to the ecosystem stability and function [4,7]. Therefore, studies on the community composition and dynamics of estuarine ciliates will improve our understanding of how the distribution and function of ciliates contribute to estuarine ecosystems and inform future assessments of the environmental status of estuaries.
Data on estuarine ciliate communities are rather scant, but previous studies indicate remarkable diversity [8][9][10][11][12]. In general, choreotrichs, oligotrichs, scuticociliates, and tintinnids were found to be the most abundant assemblages [11][12][13][14]. Other studies focusing on particular groups of ciliates, such as tintinnids, have reported seasonal changes in species richness, with the highest level of richness in different seasons for different locations [11,15,16]. Most of these studies used traditional microscopy approaches to identify and count species. However, most ciliate species cannot be cultured and do not survive in the laboratory long enough to be accurately identified. Moreover, fixation can alter the morphological characteristics of ciliates (such as their shape, color, and size), with cells often mistaken for detrital materials [17][18][19]. In recent years, high-throughput sequencing has been used to study ciliate community structures in many aquatic environments, including freshwater, deep sea, and polar environments [20][21][22][23]. These studies revealed greater taxonomic richness, including non-culturable and rare taxa, compared with microscopic surveys [24]. Indeed, there have been relatively few studies on estuarine ciliates, despite recent molecular studies revealing that they account for most of the microeukaryote diversity and abundance in this environment [25,26]. Therefore, further molecular studies are needed to better evaluate the diversity and seasonal dynamics of ciliates in estuarine environments.
In this study, 48 seasonal samples of surface water from four estuaries in coastal zones of the northern area of the South China Sea in Shenzhen, southeastern China, were analyzed by 18S rDNA high-throughput sequencing to explore the diversity and seasonality of estuarine ciliate communities. The results improve our understanding of the diversity and dynamics of ciliate communities in estuarine ecosystems and provide a reference for the ecological assessment of estuary, marine, and coastal ecosystems using ciliates as indicator organisms.

Study Sites and Sampling
The four estuarine sampling stations (S1-S4) used in this study were located in three bays in Shenzhen, China ( Figure 1): S1 is the estuary of an unnamed urban river that forms a brackish open bay; S2 is an estuary of the Xixiang River, which is located in the semienclosed DaChan Bay; and S3 and S4 are located in Shenzhen Bay, also a semi-enclosed bay. S3 is an estuary of the Dasha River and S4 is a joint estuary of the Shenzhen and Xinzhou rivers. S1, S3, and S4 are shallow coastal waters within urban parks; in contrast, S2 is a deeper bay that is used as a fishing port, and thus might be polluted. The average annual temperature for all four estuaries is 25.6 • C. The highest temperatures occur in June-August and the coldest months are from December-February. Most rainfall occurs from May to September, when short rainstorms often occur.
In June 2019 (summer), September 2019 (autumn), December 2019 (winter), and March 2020 (spring), surface water samples were collected using sterilized plastic bottles and processed immediately in the laboratory. For each sample, 1 L of water was prefiltered with a 200 µm nylon mesh to remove large particles, and mesozooplankton and ciliate cells were collected by filtration onto 0.2 µm pore polycarbonate membranes (Millipore, Darmstadt, Germany). The membranes were immediately frozen and stored at −80 • C until DNA extraction. Three replicate samples were collected from each station, for a total of 48 samples. Environmental factors (water temperature, dissolved oxygen concentration, salinity, and pH) were measured at each station using a YSI water quality instrument (YSI, Yellow Springs, OH, USA). total of 48 samples. Environmental factors (water temperature, dissolved oxygen concentration, salinity, and pH) were measured at each station using a YSI water quality instrument (YSI, Yellow Springs, OH, USA).

DNA Extraction, Polymerase Chain Reaction (PCR) Amplification, and High-Throughput Sequencing
Environmental DNA was extracted from all membranes using the PowerSoil DNA Isolation Kit (Qiagen, Germantown, MD, USA) according to the manufacturer's protocol. The DNA quality was checked by 1% agarose gel electrophoresis, and the DNA concentration and purity were determined using a NanoDrop 2000 UV-Vis spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). The hypervariable region V4 of the ciliate 18S rRNA was amplified using ciliate-specific primers for the 18S rDNA and V4 region [27]. PCR products were resolved by 2% agarose gel electrophoresis, purified using an AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) according to manufacturer's protocol, and quantified using a QuantiFluor ST Fluorometer (Promega, Fitchburg, MA, USA). Purified amplicons were pooled in equimolar concentrations and paired-end sequenced on an Illumina MiSeq PE300 platform using standard protocols. Raw sequence data are available in the Genome Sequence Archive (GSA, https://bigd.big.ac.cn/gsa/ (accessed on 2 December 2020)) of the China National Center for Bioinformation (CNCB) under BioProject accession number PRJCA003919.

Sequence Data Processing
Data analysis was performed using the QIIME2 pipeline [28], with additional analyses in R. Raw reads were denoised using the DADA2 plugin with default parameters and settings [29] and then clustered into representative amplicon sequence variants (ASVs). Reads presenting as a single copy were removed. The final dataset contained 1,849,216 sequences, representing 1645 ASVs for further analysis. Taxonomic assignment was performed against the SILVA 138 ribosomal RNA gene database using a confidence threshold of 70% [30].

DNA Extraction, Polymerase Chain Reaction (PCR) Amplification, and High-Throughput Sequencing
Environmental DNA was extracted from all membranes using the PowerSoil DNA Isolation Kit (Qiagen, Germantown, MD, USA) according to the manufacturer's protocol. The DNA quality was checked by 1% agarose gel electrophoresis, and the DNA concentration and purity were determined using a NanoDrop 2000 UV-Vis spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). The hypervariable region V4 of the ciliate 18S rRNA was amplified using ciliate-specific primers for the 18S rDNA and V4 region [27]. PCR products were resolved by 2% agarose gel electrophoresis, purified using an AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) according to manufacturer's protocol, and quantified using a QuantiFluor ST Fluorometer (Promega, Fitchburg, MA, USA). Purified amplicons were pooled in equimolar concentrations and paired-end sequenced on an Illumina MiSeq PE300 platform using standard protocols. Raw sequence data are available in the Genome Sequence Archive (GSA, https://bigd.big.ac.cn/gsa/ (accessed on 2 December 2020)) of the China National Center for Bioinformation (CNCB) under BioProject accession number PRJCA003919.

Sequence Data Processing
Data analysis was performed using the QIIME2 pipeline [28], with additional analyses in R. Raw reads were denoised using the DADA2 plugin with default parameters and settings [29] and then clustered into representative amplicon sequence variants (ASVs). Reads presenting as a single copy were removed. The final dataset contained 1,849,216 sequences, representing 1645 ASVs for further analysis. Taxonomic assignment was performed against the SILVA 138 ribosomal RNA gene database using a confidence threshold of 70% [30].

Statistical Analysis
Statistical analyses were conducted in R-3.3.1 (https://cran.r-project.org/ (accessed on 21 June 2016)). Diversity indices (richness, Simpson index, shannoneven) were calculated in R using the vegan and OTU table packages. The beta diversity of 48 samples from four stations was analyzed using non-metric multidimensional scaling analysis (NMDS) based on binary Euclidean distance matrices of ASV level annotation in R. To explore the correlation between the environmental factors (water temperature, dissolved oxygen, salinity, and pH) and the ciliate communities, canonical correspondence analysis (CCA) was performed in R. The correlations between richness, evenness, and community and environmental parameters were calculated using the Pearson correlation coefficient and considered significant at r > ±0.75 and p < 0.05. Multiple samples were compared using the nonparametric Kruskal-Wallis test, followed by Dunn's multiple comparisons test.

Environmental Parameters of Sampling Stations
During the sampling period (June 2019-March 2020), environmental parameters (water temperature, dissolved oxygen concentration, salinity, and pH) varied across all four estuarine ecosystems ( Figure S1). The water temperature ranged from 20.3 to 29.9 • C; as expected, it was higher in summer (June) and autumn (September) than in winter (December) and spring (March) (mean temperature, 29.4 vs. 21.8 • C). The dissolved oxygen concentration ranged from 3.3 to 6.6 mg/L across all samples, with most lower than 5.0 mg/L. The salinity varied more significantly across seasons, ranging from 2.9 to 19.5‰; low salinity was recorded in four sampling stations in summer and two in autumn (mean, 3.4 vs. 14.3 mg/L), probably because of increased rainfall in summer. At all stations, the pH of samples was relatively stable (7.6-8.3), with no apparent variations.

Overview of Sequencing Data
The diversity and community of ciliates in 48 water samples from all four seasons in four estuarine ecosystems was based on 18S rDNA V4 region amplicon sequencing. A total of 1,849,216 high-quality ciliate V4 sequences were clustered into 1645 ASVs. The number of sequences per sample varied between 2391 (S4_DEC_1) and 58,172 (S1_DEC_3), with an average of 38,525. The total number of ASVs obtained and used in further downstream analyses ranged between 23 (S4_MAR_3) and 306 (S4_SEP_3), with an average of 97.

Distribution of Ciliate ASVs among Sampling Stations
Overall, the four sampling stations differed in ciliate richness, diversity, and evenness ( Figure 3). The richness was higher in S1 than in S3 and S4, and lowest in S2. The same pattern was observed for diversity. The evenness had a similar pattern: on average, it was slightly higher in S2 than in S1 and S3, and lowest in S4. The ciliate communities also differed among the sampling stations ( Figure 4 and Figure S2). Only 82 of the total 1645 ASVs (4.98%) were common to samples from all four stations. Spirotrichea and Oligohymenophorea accounted for more than three-quarters of the 82 ASVs, followed by Litostomatea, Prostomatea, Heterotrichea, Nassophorea, Colpodea, and Phyllopharyngea. However, the latter group was not always abundant or simultaneously present in all samples. The most abundant species were from the genera Halteria, Parastrombidinopsis, Pelagostrobilidium, Spirotontonia, Strombidium, Tintinnidium, Tintinnopsis, an unclassified Choreotrichia species within Spirotrichea, Cryptocaryon within Prostomatea, Pseudovorticella, two unclassified species within Oligohymenophorea, and one Zosterodasys species within Nassophorea (>1% of the total reads). Furthermore, 122 ASVs were found in samples from three stations and 310 ASVs in samples from two stations. In addition, 1131 station-specific ASVs (68.75% of all ASVs; 346 for S1, 192 for S2, 327 for S3, and 266 for S4) comprised 39-49% of the total ASVs for each station (Figure 4). The community composition of station-specific taxa was different for each station. Spirotrichs (including choreotrichs, oligotrichs, and stichotrichs) and prostomateans (Cryptocaryon) were the most common dominant groups in the specific community for each station. According to our assigned genera, however, the most abundant genera differed among stations: Tintinnidium and Pseudouroleptus in S1; Helicostomella, Favella, and Tunicothrix in S2; Strombidium, Strobilidium, and Urospinula in S3; and Tunicothrix and Pseudouroleptus in S4. Hypotrichs were dominant in S1-S3 and absent in S4. In S2, scuticociliates, peritrichs, and peniculines of the class Oligohymenophorea and Odontochlamys of the class Phyllopharyngea were also abundant; they were represented by the genera Paramecium, Telotrochidium, Zoothamnium, Protocyclidium, Metanophrys, Cinetochilum, and Cyrtophoria. The specific abundant taxa in S2 were peniculines Stokesia, haptorians Teuthophrys, and nassophoreans Zosterodasys. Zosterodasys was also a dominant S4-specific genus, along with prostomateans Placus and Hemiophrys, and some unclassified haptorians.

Seasonal Variations in Taxonomic Groups
We pooled the annual data into 16 seasonal samples: four seasons for each of the four sampling stations. In all, the taxonomic richness was highest in autumn and summer, and lowest in spring. The ciliate diversity had a different seasonal pattern in each station ( Figure 5): in S1, it was highest in summer and lowest in spring; in S2, it was highest in winter and lowest in autumn; and in S3 and S4, it was highest in autumn and lowest in winter (S3) and spring (S4). The succession of the 13 ciliate classes in the four stations over the sampling period was determined as a proportion of reads and ASVs ( Figure S2). Of the 13 classes, nine were found in samples from all stations, and only five were detected in all 16 seasonal samples: Litostomatea, Oligohymenophorea, Phyllopharyngea, Prostomatea, and Spirotrichea. Spirotrichea and Oligohymenophorea were, as a rule, the groups with the most richness throughout the year. However, there were two exceptions: the classes Litostomatea and Prostomatea had the second highest richness in the autumn sample from S4, autumn and winter samples from S3, and the spring sample from S1. The classes with lower richness across all four stations were Nassophorea, Colpodea, Heterotrichea, and Armophorea. Moreover, these classes were absent in two of the 16 seasonal samples: Nassophorea (S2_SEP and S3_DEC), Colpodea (S1_DEC and S4_DEC), and Heterotrichea (S3_MAR and S4_MAR). Armophorea was absent in the autumn samples from stations S1 and S2, and the winter samples from stations S2-S4. The other classes, Plagiopylea, Protocruziea, Karyorelictea, and Licnophoriea, were defined as rare taxa; Plagiopylea and Protocruziea were identified only in the winter samples from stations S1 and S2 and the autumn sample from station S3, Karyorelictea in only the winter samples from stations S1 and S2, and Licnophoriea in only the winter samples from station S1 (Figure 6).

Seasonal Variations in Taxonomic Groups
We pooled the annual data into 16 seasonal samples: four seasons for each of the four sampling stations. In all, the taxonomic richness was highest in autumn and summer, and lowest in spring. The ciliate diversity had a different seasonal pattern in each station (Figure 5): in S1, it was highest in summer and lowest in spring; in S2, it was highest in winter and lowest in autumn; and in S3 and S4, it was highest in autumn and lowest in winter (S3) and spring (S4). The succession of the 13 ciliate classes in the four stations over the sampling period was determined as a proportion of reads and ASVs ( Figure S2). Of the 13 classes, nine were found in samples from all stations, and only five were detected in all 16 seasonal samples: Litostomatea, Oligohymenophorea, Phyllopharyngea, Prostomatea, and Spirotrichea. Spirotrichea and Oligohymenophorea were, as a rule, the groups with the most richness throughout the year. However, there were two exceptions: the classes Litostomatea and Prostomatea had the second highest richness in the autumn sample

Seasonal Variations in Taxonomic Groups
We pooled the annual data into 16 seasonal samples: four seasons for each of the four sampling stations. In all, the taxonomic richness was highest in autumn and summer, and lowest in spring. The ciliate diversity had a different seasonal pattern in each station (Figure 5): in S1, it was highest in summer and lowest in spring; in S2, it was highest in winter and lowest in autumn; and in S3 and S4, it was highest in autumn and lowest in winter (S3) and spring (S4). The succession of the 13 ciliate classes in the four stations over the sampling period was determined as a proportion of reads and ASVs ( Figure S2). Of the 13 classes, nine were found in samples from all stations, and only five were detected in all 16 seasonal samples: Litostomatea, Oligohymenophorea, Phyllopharyngea, Prostomatea, and Spirotrichea. Spirotrichea and Oligohymenophorea were, as a rule, the groups with the most richness throughout the year. However, there were two exceptions: the classes Litostomatea and Prostomatea had the second highest richness in the autumn sample The relative read abundance varied throughout the seasons of the annual survey. Spirotrichea was the most abundant group throughout the year, except in three samples, where Oligohymenophorea (S2_JUN and S2_SEP) and Nassophorea (S4_SEP) were the most abundant. The temporal dynamics of the 24 most dominate ASVs were assessed using the nonparametric Kruskal-Wallis test; these ASVs were not always abundant in a particular season or present in a similar pattern in all four stations ( Figure S3). At the lower taxonomic levels, Parastrombidinopsis, Spirotontonia, and Halteria were the dominant genera in the cooler seasons (spring and winter), whereas Tintinnidium, Tintinnopsis, Didinium, and Zosterodasys were more abundant in the warmer seasons (summer and autumn). and S2, and the winter samples from stations S2-S4. The other classes, Plagiopylea, Protocruziea, Karyorelictea, and Licnophoriea, were defined as rare taxa; Plagiopylea and Protocruziea were identified only in the winter samples from stations S1 and S2 and the autumn sample from station S3, Karyorelictea in only the winter samples from stations S1 and S2, and Licnophoriea in only the winter samples from station S1 (Figure 6).  The relative read abundance varied throughout the seasons of the annual survey. Spirotrichea was the most abundant group throughout the year, except in three samples, where Oligohymenophorea (S2_JUN and S2_SEP) and Nassophorea (S4_SEP) were the most abundant. The temporal dynamics of the 24 most dominate ASVs were assessed using the nonparametric Kruskal-Wallis test; these ASVs were not always abundant in a particular season or present in a similar pattern in all four stations ( Figure S3). At the lower taxonomic levels, Parastrombidinopsis, Spirotontonia, and Halteria were the dominant genera in the cooler seasons (spring and winter), whereas Tintinnidium, Tintinnopsis, Didinium, and Zosterodasys were more abundant in the warmer seasons (summer and autumn). and S2, and the winter samples from stations S2-S4. The other classes, Plagiopylea, Protocruziea, Karyorelictea, and Licnophoriea, were defined as rare taxa; Plagiopylea and Protocruziea were identified only in the winter samples from stations S1 and S2 and the autumn sample from station S3, Karyorelictea in only the winter samples from stations S1 and S2, and Licnophoriea in only the winter samples from station S1 (Figure 6).  The relative read abundance varied throughout the seasons of the annual survey. Spirotrichea was the most abundant group throughout the year, except in three samples, where Oligohymenophorea (S2_JUN and S2_SEP) and Nassophorea (S4_SEP) were the most abundant. The temporal dynamics of the 24 most dominate ASVs were assessed using the nonparametric Kruskal-Wallis test; these ASVs were not always abundant in a particular season or present in a similar pattern in all four stations ( Figure S3). At the lower taxonomic levels, Parastrombidinopsis, Spirotontonia, and Halteria were the dominant genera in the cooler seasons (spring and winter), whereas Tintinnidium, Tintinnopsis, Didinium, and Zosterodasys were more abundant in the warmer seasons (summer and autumn).

Ciliate Diversity in Relation to Environmental Factors
Correlation analysis was performed to investigate whether environmental variables were related to ciliate diversity and community. Ciliate richness was slightly related to salinity; evenness was slightly related to water temperature, salinity, and pH; and ciliate community composition correlated with pH. Statistical analysis revealed that dissolved oxygen concentration had no significant relationship with these variables (Table 1, Figure S4). We next evaluated the association between environmental variables and the dominant ASVs (>1% of the total reads per ASV) by correlation analysis ( Figure S5). Twenty-four dominant ASVs were clustered into four clades (I-IV). Three ASVs in clade I had a weak inverse correlation with pH. Clade II included six ASVs; most of these were inversely correlated with water temperature and directly correlated with salinity and dissolved oxygen concentration. The ASV from Parastrombidinopsis sp.1 was significantly correlated with salinity and dissolved oxygen concentration, and the ASVs from Spirotontonia turbinate, Askenasia sp., and Pelagostrobilidium sp.1 had a strong inverse correlation with water temperature. ASVs in clade III had a weak correlation with water temperature, dissolved oxygen concentration, and salinity, a direct correlation with pH, especially Parastrombidinopsis sp.2. The other 12 ASVs were clustered in clade IV; most of these had a direct correlation with water temperature and an inverse correlation with dissolved oxygen concentration and salinity (significant for the ASV from Tintinnidium balechi).

Discussion
Determining the diversity and community composition of eukaryotic ciliate microbes throughout the year is critical to understand the structure and function of estuarine ecosystems. High-throughput sequencing is efficient to reveal deeper knowledge of hidden ciliate diversity [17,[32][33][34]. This is the first long-term study of ciliate communities in four estuaries in Shenzhen by metabarcoding. A total of 1645 ASVs were identified. Some ASVs matched several reference sequences from the same species because of polymorphism in the 18S rDNA sequence; therefore, the number of ASVs probably overestimates the species diversity. Furthermore, most ciliates are difficult to culture and/or identify through microscopy in vivo, and molecular references are not available for these ciliates. Thus, in this study ciliates lacking molecular references were designated as unclassified taxa (1.64% of the ASVs).
Planktonic ciliates in estuaries are not well studied, and knowledge of the diversity and composition of estuarine communities is limited. Our results, based on a meta-analysis of four estuaries, revealed a high diversity of ciliates in this ecosystem. In an annual survey with a total of 48 samples, we identified 1645 ASVs belonging to 13 ciliate classes, 97 families, 157 genera, and 207 defined species; 1212 ASVs were assigned to unclassified or uncultured species. Species in the classes Spirotrichea, Oligohymenophorea, Prostomatea, Litostomatea, and Phyllopharyngea were identified in all four stations throughout the sampling period, representing about 96.86% of the total ciliate abundance. Spirotrichs (including choreotrichs, oligotrichs, and stichotrichs) and oligohymenophorean scuticociliates were the dominant ciliate taxa in estuarine environments in the present survey. Furthermore, sessile peritrichs of the genus Pseudovorticella and haptorian Didinium, nassophorean Zosterodasys, and prostomatean Askenasia genera were also abundant. These taxa were previously reported to be the dominant ciliated microorganisms in estuaries and coastal marine environments [4,13,[35][36][37]. The ciliates are capable of tolerating extreme changes of salinity, and some of them can sustain direct transfer from marine to fresh water [38][39][40]. Therefore, ciliate communities in estuarian environments probably include freshwater, brackish, and marine species, and they may have potential salinity tolerance.
The four estuarine sampling stations were located in three bays of the northern South China Sea. Ciliate diversity and community composition were significantly different at all four stations (Figure 4), even though stations S3 and S4 are both located in Shenzhen Bay. However, differences in urban runoff are likely to be an important factor in the estuarine water environment and diversity of ciliate microorganisms. Only 82 ASVs were identified in all four stations; of these, 68% were present in higher abundance (>0.05% of the total reads per ASV; Figure S6A). These ASVs correspond to widespread taxa with a strong ability to adapt to various environments that are common in estuaries and coastal waters [13,25,26,35,41,42]. In all, 7% of the total ASVs were identified in two of the sampling stations and 19% in three, corresponding to increasing numbers of ASVs with a lower abundance ( Figure S6B,C). Of the total ASVs, 68.75% were specific to a single station.
These tended to have a lower abundance; the relative abundance of 94% of these ASVs was <0.02% ( Figure S6D). Differences in ciliate communities among sampling stations seemed to be due to these rare ASVs rather than to more abundant ASVs. Indeed, rare taxa were found to be important contributors to diversity and have important ecological roles [43]. These rare taxa are hypothesized to be functionally redundant and sensitive to environmental perturbations or changes, and can increase in abundance to maintain the function of the ecosystem [44].
Our data showed that the diversity and community composition of ciliates from four estuarine stations underwent strong changes in different seasons ( Figure 5, Figure 6, and Figure S7). However, distinct seasonal variations in the ciliate community among the four stations indicated that environmental factors such as water temperature, salinity, and pH correlate with changes in ciliate diversity and community composition. Previous studies suggested that estuarine ciliate diversity is influenced by temperature and salinity [12,25,45,46]. In coastal areas, nutrients and weather changes (e.g., rainfall and wind patterns) also influence estuarine ciliate diversity and communities [21,47]. Furthermore, bacteria, planktonic algae, and mesozooplankton also have considerable effects on the structure of ciliate communities [43,48,49]. These various habitats of estuarine ecosystems respond differently to environmental conditions, thus further influencing ciliate diversity and community composition. Further studies are necessary in order to determine the specific correlation between ciliate communities and other abiotic and biotic parameters in estuarine ecosystems.

Conclusions
This study analyzed the seasonal dynamics of ciliate diversity and community structures in four estuaries in similar geographical areas. The results show that estuarine ecosystems contain a high diversity of ciliates, including members of most ciliate classes. Spirotrichs (including choreotrichs, oligotrichs, and stichotrichs), oligohymenophorean scuticociliates, haptorian Didinium, and prostomatean Cryptocaryon were the dominant taxa in the four sampling stations. The rarer taxa were important members of estuarine ciliate communities. Differences and seasonal changes in the composition of ciliate communities were found across the estuarine ecosystems. Sampling position variations had more influence on ciliate community composition than seasonal changes. The diversity and community composition of ciliates correlated with temperature, salinity, and pH. Overall, the present study improves our understanding of the diversity and seasonal dynamics of ciliate communities in estuarine ecosystems and provides a reference for the ecological assessment of estuarine, marine, and coastal ecosystems using ciliates as indicator organisms.