Diversity, Composition, and Activities of Nano- and Pico-Eukaryotes in the Northern South China Sea With Influences of Kuroshio Intrusion

Nano- and pico-eukaryotes play important roles in the diversity and functions of marine ecosystems. Warm, saline, and nutrient-depleted water that originates in the Kuroshio Current seasonally intrudes into the northern South China Sea (NSCS) from autumn to spring. To clarify the mechanisms in shaping the community structure of nano- and pico-eukaryotes as well as impacts of the Kuroshio intrusion on the NSCS ecosystem, genomic DNA and RNA were co-extracted from samples collected at two depths from nine stations, and then the V9 region of 18S rDNA and rRNA was sequenced with high-throughput sequencing. Our results showed that Dinophyceae was the most diverse and abundant nanoeukaryotic group during the study period revealed by both DNA and RNA surveys. In contrast, the relative read abundance of MAST, Pelagophyceae, and Dinophyceae in the size fraction of picoeukaryotes might be largely underestimated by the DNA survey. The RNA survey was the more reliable method to investigate the eukaryotic community structure. Environmental filtering played an important role in shaping the community structure, and the sampling depth became the governing factor of the beta diversity under the environmental setting of stratification during the study period. The spatial variations in the diversity of nanoeukaryotes were subject to the dispersal limitation under the size rule. The effects of the Kuroshio intrusion on the nanoeukaryotic community structure might also be explained by the dispersal limitation. Overall, neutral processes are critical in shaping the community structure of nanoeukaryotes. The relative metabolic activities of nanoeukaryotes were relatively stable in accordance with the high similarity of community structure between sampling sites. The responses of the relative metabolic activities of picoeukaryotes to environmental factors displayed two distinct patterns: positive correlations with salinity and nutrients and negative with temperature for Dinophyceae, MAST, and Pelagophyceae, while reversed patterns for Mamiellophyceae and Radiolaria. Our findings improve the understanding of the nano- and pico-eukaryotic communities in the NSCS and the mechanisms of their assembly.


INTRODUCTION
Nano-and pico-eukaryotes (size fractions of 3-20 and 0.2-3 µm, respectively) are present worldwide, abundantly, and persistently in various marine ecosystems (Massana et al., 2015). Much attention has been paid to their astonishing diversity (Caron et al., 2012;Piwosz et al., 2015;de Vargas et al., 2015;Caron et al., 2016) and more recently to their trophic modes and interactions (Piwosz et al., 2018;Adl et al., 2019;Sassenhagen et al., 2019). The photosynthetic groups of nano-and pico-eukaryotes are important contributors to primary production (Grob et al., 2007;Jardillier et al., 2010;Hartmann et al., 2012), which link trophic levels in the microbial food web , and mixotrophic groups are important bacterial grazers and have the potential to dominate the primary production and bacterivory in marine ecosystems (Unrein et al., 2014). Understanding the diversity, composition, and functions of these groups is therefore important to unveil critical biogeochemical cycles in the sea.
The main goal of microbial ecology is to reveal the mechanisms shaping and maintaining community assembly (Zhou and Ning, 2017). The framework described by Sloan et al. (2006) proposed four main processes for regulating the distribution patterns of microbial communities, including: selection, ecological drift, dispersal, and speciation. These four processes can be further grouped into two categories: deterministic processes (or selective processes) and neutral processes (or stochastic processes) (Sloan et al., 2006;Logares et al., 2013). The deterministic process refers to the governing of the microbial community structure by environmental conditions (e.g., temperature, salinity, and nutrient concentration) and biotic interactions (e.g., predation, competition, and mutualism) (Stegen et al., 2012;Yu et al., 2015;Zhou and Ning, 2017). The neutral process refers to the governing of the microbial community structure by the stochastic process of birth, death, extinction, colonization, and speciation (Hubbell, 2011). The neutral community model (NCM) is an adaptation of the neutral theory, which is very useful in quantifying the importance of neutral process in diverse environments (Logares et al., 2013;Chen et al., 2017). Normalized stochasticity ratio (NST), proposed by Ning et al. (2019), is also powerful in evaluating the community assembly processes . The relative importance of deterministic and neutral processes varies among regions, taxa, and ecosystems (Bahram et al., 2016;Wu et al., 2017b).
Eukaryotic diversity in the ocean has been extensively studied. Most studies targeted one specific size fraction, while a few studies explored diversity across size, usually in four size classes: picoplankton (0.2-3 µm), nanoplankton (3-20 µm), microplankton (20-200 µm), and mesoplankton (200-2000 µm) (de Vargas et al., 2015;Massana et al., 2015;Ramond et al., 2019). These studies disclosed cryptic diversity across oceans with many parasitic and phagotrophic species unseen by traditional methods. The deeply sequenced samples help reacquaint the trophic modes of eukaryotes in the ocean (Ramond et al., 2019). In spite of the dominance of heterotrophs in every size fraction, the ratio between heterotrophs and autotrophs increases with the nominal size of fractions (de Vargas et al., 2015), and cell size constrains fundamental requirements and physiological rates (e.g., growth rate) of plankton as well as biotic interactions (e.g., grazing and competition) between them ( Barton et al., 2013). In contrast with picoplankton studied by most research, nanoplankton have higher requirements for nutrients and food, in turn, nanoplankton contribute more to carbon export than picoplankton (Ward and Follows, 2016;Leblanc et al., 2018). Cell size is also an important factor for dispersal limitation to shape community structure. Dispersal limitation was reported to increase with cell size (de Vargas et al., 2015).
High-throughput sequencing (HTS) has revealed the astonishing diversity of protist (microbial eukaryotes) and greatly improved our understanding of their community composition in environmental samples (Edgcomb et al., 2011;Parris et al., 2014;Massana et al., 2015;de Vargas et al., 2015;Piwosz et al., 2018;Liu et al., 2019). Moreover, protistologists have generally accepted that protist assemblages inferred from 18S rRNA gene sequencing (hereafter DNA sequencing) represented the species present. Compared to DNA, RNA is unstable in extracellular conditions, microbial community structures inferred from rRNA sequencing (hereafter RNA sequencing) can exclude sequences from dissolved extracellular pool, dormant cells and dead materials . Thus, the RNA-based approach was used to identify the active species or groups in environmental samples (Massana et al., 2015;Hu et al., 2016;Wang et al., 2019). The RNA: DNA ratio (calculated by relative sequence abundance) has been used as an index for the relative metabolic activity in protist (Hu et al., 2016;Wu and Liu, 2018). However, sequences information derived from the rRNA indicates the possible ribosomal activity (i.e., protein synthesis), which is not a direct indicator of cell activity (Blazewicz et al., 2013). In addition, gene copy numbers, life strategies and histories are very different among species. The RNA: DNA ratio should be used with caution (Blazewicz et al., 2013). Comparison of RNA: DNA ratios usually used in the same species or group among different samples, but not for the comparison between different species or groups (Hu et al., 2016). A previous study showed that the RNA: DNA ratio could be a potential indicator of the upcoming Phaeocystis globosa bloom in the eastern English Chanel (Rachik et al., 2018). Xu et al. (2017) observed the varied metabolic activities of protist community at different sampling depths in the South China Sea. Wang et al. (2019) found the varied metabolic activities of picoeukaryotes at spatial scale in the Pacific Ocean.
The South China Sea is the largest marginal sea of the North Pacific (from the equator to 23 • N, from 99 to 121 • E), and its surface circulation is driven by monsoons (Liu et al., 2002). The North South China Sea (NSCS) exchanges water with the North Pacific through only the Luzon Strait channel. Kuroshio, the western boundary current of the North Pacific, is characterized by warm and saline waters and intrudes seasonally into the NSCS by passing through the Luzon Strait (Nan et al., 2015). Intrusion also generates mesoscale eddies in the NSCS . The mix of Kuroshio and NSCS waters dilutes nutrients (Du et al., 2013), but has been found to stimulate the activities of ammonia-oxidizing bacteria (Xu et al., 2018;Lu et al., 2019) as well as primary productivity . The diversity and composition of picoeukaryotes in the NSCS have been previously explored by high-performance liquid chromatography (Xiao et al., 2018), fluorescence in situ hybridization (Wu et al., 2017c), and DNA sequencing (Wu et al., 2014). However, little is known about the diversity and composition of nano-eukaryotes in this area. Moreover, neutral processes have been overlooked in the past, but these processes could be important for Kuroshiowater-influenced areas because ocean currents are always related to the dispersal and colonization of species. This study examines the diversity and community composition of nano-and picoeukaryotes in the NSCS using size-fractionated filtration and DNA/RNA high-throughput sequencing. The NSCS was sampled in late spring, when a significant proportion of the water mass is of Kuroshio origin, in order to evaluate the relative importance of environmental filtering and neutral processes in shaping communities. In particular, we assess how diversity and community composition of nano-and pico-eukaryotes are affected by deterministic processes by relating metabolic activity and environmental gradients.

Cruise, Sampling, and Environmental Parameters
The research cruise was conducted in the NSCS from 114 to 119 • E and 18 to 22.5 • N from May 15th to June 5th, 2016 (Supplementary Figure 1). A total of 36 samples (18 picoeukaryotes and 18 nanoeukaryotes) were collected from the surface layer (Sur, 5 m) and deep chlorophyll maximum (DCM) layers at nine stations (Supplementary Figure 1). For each sample, 10-20 L of seawater was collected using Niskin bottles (12-L PVC) mounted on a conductivity, temperature, and depth (CTD) rosette (Sea-Bird Electronics, United States). The seawater was pre-filtered through a 20 µm Bolting Cloth to remove large organisms. Then the pre-filtered seawater was sequentially filtered through 3-and 0.2 µm pore size polycarbonate membranes (142 mm diameter, Millipore, United States) using a peristaltic pump. To reduce the loss of RNA samples, the filtration rate in our experiment was controlled at 0.5-1 L/min. The filters (nanoeukaryotic samples: 3-20 µm, picoeukaryotic samples: 0.2-3 µm) were stored in 5 mL cryotubes containing 4.6 mL RNAlater (Thermo Fisher, Lithuania). The tubes were immediately frozen in liquid nitrogen. In this study, samples were named by the combination of sampling station and depth, for example, D2Sur means the sample was collected at the surface layer at Station D2. Seawater temperature and salinity at each sampling site were recorded using the SBE-911 CTD. Water samples (50 mL) were also collected at the sampling sites and preserved at −20 • C for nutrient analysis. For chlorophyll a (Chl a) measurement, 1 L of seawater from each sampling depth were filtered onto 25 mm GF/F glass fiber filters (Whatman, United States) under 200 mm Hg pressure. The filters were stored in liquid nitrogen on board. In the laboratory, the filters were submersed in 90% acetone in the dark at −20 • C for 24 h. The concentrations of Chl a extracted into the acetone were measured with a Trilogy fluorometer (Turner Design, United States) (Welschmeyer, 1994). A geostrophic current is the result of the horizontal pressure gradient and the Coriolis force. The direction of the geostrophic current is parallel to the contour of the sea surface height. To generate a synoptic view of how Kuroshio intrusion physically influences the SCS, sea level anomalies and geostrophic current vectors for the study area were downloaded from the US NOAA National Environmental Satellite, Data, and Information Service. 1

Kuroshio Fraction at the HTS Sampling Stations
An isopycnal mixing model described by Du et al. (2013) was used to quantify the impact of the Kuroshio intrusion (KI) at our sampling stations. This model assumes that the different water masses are dominated by isopycnal mixing, and diapycnal mixing is comparatively negligible (Du et al., 2013). In this model, two end-members (representing the Kuroshio water and SCS water, respectively, in this study) are required for index estimation. The P4 station (122.959 • E, 20.003 • N) was used to represent the Kuroshio water mass end-member, and the South East Asian time-series (SEATS) station (116 • E, 18 • N) was selected to represent the proper South China Sea (SCS) water mass end-member ( Figure 1A). Notably, the fractional contribution of the two end-members varies depending on the end-member selection. The proportion of the Kuroshio water (R K ) indicated the relative contribution from P4 instead of the "typical Kuroshio" water. Our HTS samples were all collected at the surface (5 m) and DCM layers (less than 100 m), and the average R K values of the upper 100 m (R Kave−100 , 1 m interval) were used to represent the KI fraction. Hence, an index was used to indicate the degree of KI influence (KI index) at each sampling station (Supplementary Table 1). According to the assumption of this model, Xu et al. (2018) who participated in the same cruise have calculated and published the KI index at all sampling stations. In this study, KI index ≥0.35 was considered to indicate more KI (MKI) affected stations, and the D2, D6, A3, A7, and K4 stations belonged to this category. KI index <0.35, was considered to indicate less KI (LKI), and 2A1, A10, SEATS, and K8 stations belonged to this category ( Figure 1A).

DNA and RNA Co-extraction, PCR Amplification, and Illumina High-Throughput Sequencing
For each sample, total DNA and RNA were extracted simultaneously using the AllPrep DNA/RNA Min Kit (Qiagen, #80204, Germany). The RNA extracts were purified using the gDNA wipeout buffer (Qiagen, #205311) to remove possible genomic DNA contamination, and then the purified RNA was reverse transcribed into cDNA using the QuantiTect Reverse Transcription Kit (Qiagen, #205311, Germany). The V9 hypervariable region of the eukaryotic 18S rRNA gene was amplified using the 1380F (5-CCCTGCCHTTTGTACACAC-3) and 1510R (5-CCTTCYGCAGGTTCACCTAC-3) primers (Amaral-Zettler et al., 2009). For each polymerase chain reaction (PCR) reaction, 25 µL mixture contained 12.5 µL 2 × Taq PCR mix (Takara, China), ≤1 µg DNA or cDNA template, 0.1 µM forward primer, and 0.1 µM reverse primer. The PCR program consisted of a 5-min initial denaturation step at 95 • C, followed by 34 cycles of denaturation at 94 • C for 1 min, annealing at 57 • C for 45 s, and elongation at 72 • C for 1 min, and a final extension at 72 • C for 10 min. Each PCR was performed in duplicates for each DNA and cDNA sample. The duplicate PCR products were pooled and purified using the QIAquick Gel Extraction Kit (Qiagen, #28704, Germany). Amplicons were sent with dry ice to the Novogene sequencing company (Beijing, China) where sequencing library, constructed using the TruSeq DNA PCR-Free Sample Preparation Kit, was prepared and sequenced on a single run using the Illumina HiSeq 2500 platform. High-throughput sequencing data from this study were deposited in the GenBank Sequence Read Archive database under the accession number PRJNA684778.

Sequencing Processing and Statistical Analysis
The raw tags were processed and checked using the Quantitative Insights Into Microbial Ecology (v.1.70) pipeline (Caporaso et al., 2010). Chimeras were identified and discarded after checking with the chimera search module (USEARCH, version 4.2) (Edgar et al., 2011) using the Silva 123 release (Quast et al., 2013) as a reference. The remaining filtered sequences were grouped into operational taxonomic units (OTUs) with 0.95 similarities using UPARSE (v. 7.0.100; Edgar, 2013). The most abundant sequence of each OTU was chosen as the representative sequence, and then taxonomically annotated against the V9_PR2 * database (de Vargas et al., 2015) using BLAST (E value = 10 −6 and minimum percent query coverage = 0.9) (Altschul et al., 1990). The OTU taxonomy was determined from the top hit based on the BLAST score. OTUs occurring in a single sample, OTUs with only one read (singleton), and OTUs assigned to Metazoan, Bacteria, Archaea, and organelles were also discarded from the downstream analyses. To normalize the sampling effort, the nano-and pico-eukaryotic OTU tables were rarefied to 29,799 and 28,698 reads per sample, respectively (the lowest sequence reads numbers in these tables).
All statistical analyses in our study were performed using the R software (version 3.3.2) (R Core Team, 2014). To investigate whether the community structure of nano-and pico-eukaryotes in the NSCS was affected by the Kuroshio intrusion, principal coordinate analysis (PCoA) was conducted based on the Bray-Curtis distance matrix. Variation-partitioning analysis (VPA, Borcard et al., 1992) was performed to disentangle the relative effects of spatial and environmental factors on variation in nano-and pico-eukaryotic DNA/RNA-based communities using a 3-way PERMANOVA (McArdle and Anderson, 2001). One of the environmental factors was the data matrix consisting of water temperature, NO 2 − + NO 3 − , PO 4 3− , SiO 3 2− , and Chl a, and the other was the KI index. It might represent some unmeasured environmental variables that were related to the fraction of Kuroshio water or salinity. The spatial factor was computed as the principal coordinates of neighborhood matrix (PCNM) by the pcnm function of the "vegan" package in R, the second PCNM was used for the VPA analysis in this study. The neutral community model (NCM) was used to determine the potential contribution of neutral processes to the nano-and pico-eukaryotic communities by predicting the relationship between OTUs occurrence and their relative read abundance (Sloan et al., 2006). The normalized stochasticity test (NST) was used to further confirm the importance of neutral process during nano-and pico-eukaryotic community assembly. The NST was conducted using the function of the "NST" R package, which is based on the Ruzicka index (NSTruzicka). This index ranges from 0 to 100%, indicating an increase from more deterministic to more stochastic assembly (Ning et al., 2019). Principal component analysis (PCA) was performed using the "vegan" package in R to access the distribution pattern of eukaryotic assemblages among sampling sites. OTUs only occur in both DNA and RNA sequencing results and appeared in all samples were selected, and the RNA: DNA ratios of the selected OTUs were calculated based on their relative read abundance in the RNA and DNA results (Wu and Liu, 2018).

Environmental Setting
Stations D2 and D6 were close to the Luzon Strait and characterized by higher salinity (34.44-34.81) than other sampling stations. The two stations were greatly influenced by the Kuroshio intrusion, as indicated by KI indexes larger than 0.6 ( Figure 1A). Kuroshio water affected stations A3 and A7 in the middle of the transect from 2A1 to SEATS ( Figure 1B). Although station 2A1 is geographically close to the Pearl River Estuary, it was not affected by fresh water during our sampling period. Salinity at surface and DCM layers of 2A1 station were 34.28 and 34.33, respectively, which was close to those of SEATS station. The distinct temperature between surface and DCM (paired Mann-Whitney U-test, p < 0.01) indicated stratification in water columns (Supplementary Table 2). Nutrients at the surface were mostly below the detection limits but were abundant in the DCM. The lower concentrations of Chl a at the surface (average 0.08 mg m −3 ) than at the DCM (average 0.37 mg m −3 ) (paired Mann-Whitney U-test, p < 0.01) was consistent with the different nutrient regimes.

Richness of Nano-and Pico-Eukaryotes in the NSCS
In our study, 852,873 high-quality DNA sequences and 997,272 high-quality RNA sequences were obtained from the nanoeukaryotic samples (Supplementary Tables 5, 6), whereas 798,801 high-quality DNA sequences, and 787,383 high-quality RNA sequences were obtained from the picoeukaryotic samples (Supplementary Tables 7, 8). Note that picoeukaryotic DNA amplicons in A7Sur, SEATS-Sur, and cDNA amplicons in A7Sur and 2A1Sur were too low concentration for library preparation during high-throughput sequencing, and thus these samples were not included. In this study, 3,198 and 3,233 nanoeukaryotic and picoeukaryotic OTUs, respectively, were generated. All the OTUs were affiliated to ten super-groups: Alveolata, Stramenopiles, Hacrobia, Archaeplastida, Rhizaria, Opisthokonta, Picozoa, Amoebozoa, Excavata, and Apusozoa (Figures 2A,B).
Alveolata was the richest super-group (in terms of the number of OTUs) in both nano-and picoeukaryotes, accounting for 61.2 and 61.0% of total richness (observed OTU number), respectively (Figure 2); followed by Stramenopiles (12.8 and 13.6%) and Rhizaria (9.2 and 9.5%). The other seven supergroups contributed only a small proportion of total richness. Dinophyta and Ciliophora were the top two phyla of Alveolata, Dinophyceae was the richest class of nanoeukaryotes while it was MALV-II for picoeukaryotes.
Among the 3,198 nanoeukaryotic OTUs (3,233 picoeukaryotic OTUs), there were 292 (319) unique DNA OTUs in the MKI stations and 215 (215) in the LKI stations  Figures 4A,C); similar numbers were displayed for the RNA samples (Supplementary Figures 4B,D). In addition, there were more unique nanoeukaryotic DNA OTUs at the DCM than at the surface (466 vs. 192), similar to the unique picoeukaryotic DNA OTUs (397 vs. 228).

Relative Read Abundance of Nano-and Pico-Eukaryotes
Dinophyceae was the most sequenced nanoeukaryotic class (Figures 3A,B and Table 1). The most sequenced OTU was Dinophyceae Gymnodinium_03 aureolum, representing an average relative read abundance of 12.0 and 10.0% in DNA and RNA samples, respectively, and there was no significant difference between surface and DCM or between LKI and MKI stations. In contrast, MALV-II and MAST were the most sequenced picoeukaryotic classes in DNA and RNA samples, respectively (Figures 3C,D and Table 1). Radiolaria Astrosphaera hexagonalis (OTU5) which favored the DCM environments (Kruskal-Wallis test, p < 0.05) in MKI stations (Kruskal-Wallis test, p < 0.05) was the most sequenced OTU in picoeukaryotic DNA samples (averagely 8.2%), while Pelagophyceae Pelagomonas calceolata (OTU1) was the most sequenced in RNA samples (averagely 6.7%). The relative read abundance of the picoeukaryotic OTU1 was not significantly different between MKI and LKI stations (Kruskal-Wallis test, p = 0.79) but higher in DCM than in surface for some stations (Kruskal-Wallis test, p = 0.06).
In nanoeukaryotic samples, some of minor classes showed large variations in relative read abundance. For example, The other category refers to groups with relatively lower sequence abundance.
Radiolaria was much more sequenced in DCM than in surface (Kruskal-Wallis test, p < 0.001) in DNA samples, increases of one order of magnitude occurred in DCM in stations such as A7, K4, SEATS, and K8 ( Figure 3A). Increases of one order of magnitude in DCM relative to surface also occurred to RNA samples of Bacillariophyta and MAST in a few stations ( Figure 3B). In picoeukaryotic samples, the distributions of relative read abundance among classes were not as uneven as in nanoeukaryotic samples, therefore some minor classes with large variations could dominate in a few stations. For example, the relative DNA abundance of Radiolaria was up to 67.6 and 45.8% in DCM of MKI stations A3 and A7; and the exceptional case for Dinophyceae was that it contributed 69.5% of RNA samples in DCM of the MKI station D6.

Beta Diversity of Nano-and Pico-Eukaryotic Communities and Controlling Factors
PCoA analysis showed that the beta diversity of both nano-and pico-eukaryotic communities was mainly shaped by sampling depths (Figure 4), and the surface and DCM samples were clearly separated on axis 1 while there was no apparent clustering on axis 2. Moreover, a 3-way PERMANOVA was used to assess the relative importance of environmental and spatial factors in shaping eukaryotic communities ( In contrast, NCM explained large fractions (71-81%) of variabilities in the occurrence of nanoeukaryotes and picoeukaryotes during the study period (Figure 5). The m value (immigration rate) was higher for nanoeukaryotes than for picoeukaryotes. The NST ruzicka index showed that both DNA and RNA samples of nanoeukaryotes were primarily controlled by stochastic processes (average NST ruzicka = 59 and 62%, respectively; there was no significant difference between DNA and RNA samples, Mann-Whitney U-test, p = 0.17), while the value of 38% for DNA samples of picoeukayotes indicated the dominance of deterministic processes (Figure 6; difference between nano-and pico-eukaryotes was significant, Mann-Whitney U-test, p < 0.001). However, for RNA samples of picoeukaryotes in the DCM, deterministic processes were only marginally stronger than stochastic processes (average NST ruzicka = 50.1%; difference between DNA and RNA samples was significant, Mann-Whitney U-test, p < 0.001). The above results indicated that neutral processes (stochastic processes) played a critical role in shaping the eukaryotic communities during the study period.

Variations in Relative Metabolic Activity of Nano-and Pico-Eukaryotes
PCA analysis was conducted on both the nano-and picoeukaryotic groups. For nanoeukaryotes, four groups that showed large contributions to the total variation and RNA: DNA ratios different from 1 were identified: Dinophyceae, Radiolaria, MAST, and Bacillariophyta ( Figure 7A). The average RNA: DNA ratios of Bacillariophyta, MAST, and Radiolaria were significantly higher for the DCM samples than for the surface samples ( Figure 7B, Mann-Whitney U-test, p < 0.05, Supplementary Table 4). Following the same procedure, five groups of picoeukaryotes were identified: Dinophyceae, Radiolaria, MAST, Mamiellophyceae, and Pelagophyceae ( Figure 7C and  Table 3, Mann-Whitney U-test, p < 0.05). The RNA: DNA ratios of four groups showed significant differences between the surface and DCM (Figure 7D). The average RNA: DNA ratio of Radiolaria was slightly higher for the surface than for DCM (Mann-Whitney U-test, p < 0.05). In contrast, Dinophyceae, MAST, and Pelagophyceae showed significantly higher average RNA: DNA ratios at the DCM than at the surface (Mann-Whitney U-test, p < 0.05). The relationships between the RNA: DNA ratios and the measured environmental factors (KI, salinity, water temperature, NO 2 − + NO 3 − , PO 4 3− , SiO 3 2− , and Chl a) were assessed. The RNA: DNA ratios of the nanoeukaryotic groups Dinophyceae and Radiolaria showed positive correlations with nutrients but negative correlations with the Chl a concentration ( Figure 8A and Supplementary Figure 5), although the correlation coefficients were small. The metabolic responses of picoeukaryotes to environmental factors had roughly two patterns ( Figure 8B and Supplementary Figure 6). Dinophyceae, MAST, and Pelagophyceae increased RNA: DNA ratios with salinity, nutrients, and Chl a concentration but decreased with water temperature. Meanwhile, RNA: DNA ratios of Mamiellophyceae and Radiolaria correlated negatively with nutrients and salinity but positively with water temperature.

Nano-and Pico-Eukaryotic Community Structure in the NSCS During the Study Period
In our study, the high rate of assignment for nano-and pico-eukaryotes could be that many assignments have low taxonomic ranks (e.g., only to "Eukaryota"). Dinophyceae was FIGURE 5 | The neutral community model based on the relationship between OTU occurrence and relative read abundance. OTUs with higher frequency than predicted by the model are shown in red, while those with lower frequency are shown in green. OTUs within the predicted range are displayed in black. The dotted line represents the 95% confidence interval (blue dotted line) around the model prediction.
the dominant nanoeukaryotic group and the most diverse picoeukaryotic group at all sampling sites (Figures 2, 3 and Table 1). Previous studies have revealed that many dinoflagellates were mixotrophic (Stoecker, 1999;Jeong et al., 2004), and the mixotrophic species could consume a wide variety of prey items (Seong et al., 2006;Park et al., 2006;Berge et al., 2008;Glibert et al., 2009). The most sequenced nanoeukaryotic OTU was G. aureolum, a common bloom-forming species in temperate waters that was first described by Hulburt (1957). G. aureolum contains a horseshoe-shaped apical groove, a nuclear fibrous connective (NFC), and nuclear chambers (Jeong et al., 2010). Jeong et al. (2010) have shown that G. aureolum is a mixotrophic dinoflagellate that feeds on Synechococcus, heterotrophic bacteria, and some algal species. The Dinophyceae A. sanguinea (OTU2) was the dominant nanoeukaryotic species in the DCM layer of station 2A1, which was located close to the Pearl River Estuary. A. sanguinea is also a common bloom-forming species with a wide distribution, including the Black Sea (Gómez and Boicenco, 2004), the coast of Hong Kong (Lu and Hodgkiss, 2004), the South Sea off the coast of South Korea (Lim et al., 2008), and the Ariake Sea off the coast of Japan (Katano et al., 2011). Large-scale blooms of A. sanguinea have been reported to be a significant threat to seabirds, fish, and shellfish (Jessup et al., 2009;Reis-Filho et al., 2012).
The MALV-I and II groups accounted for a significant proportion of the picoeukaryotic relative read abundance in the DNA survey. However, their relative read abundance  decreased largely in the RNA survey. Most species in the MALV-I and II groups have been characterized as parasites in a range of hosts (Bachvaroff et al., 2012). Similarly, the most sequenced OTU of picoeukaryotes was Radiolaria A. hexagonalis (OTU_5), accounting for 8.5% of total DNA sequences, but decreased sharply to 0.5% in the RNA survey. The different results between DNA and RNA surveys and the underlying reasons have been extensively discussed. Sample collection by the size-fractionated filtration may introduce artifacts, for example, cell breakage and fragments may go through filters and be retained in smaller size fractions (Massana et al., 2015). The Radiolaria A. hexagonalis is supposed to be an organism with the diameter of 150 µm (Haeckel, 1887). Its DNA sequences appearing in a smaller fraction may represent dead individuals or swarmer cells for reproduction (Li and Endo, 2020). In addition, high rDNA gene copy numbers of some species may cause them to be overrepresented in DNA samples. Massana et al. (2015) also reported that the MALV-I and II groups were overrepresented in the DNA survey for all size fractions of eukaryotes in European coastal waters and they attributed this phenomenon to higher rDNA gene copy numbers of these two groups than other eukaryotes. Overall, it is likely that the DNA survey overestimated the relative read abundance of Radiolaria, MALV-I, and MALV-II, while underestimated that of MAST, Pelagophyceae, and Dinophyceae ( Figure 3D). MAST is a remarkably diverse group of uncultured species . MASTs are heterotrophic flagellates preying on bacteria . They are widely distributed in the temperate and tropical oceans and supposed to play an important role in the microbial food web (Seeleuthner et al., 2018). The RNA survey found the average relative read abundance of 21.4% belonged to the MAST group relative to only 10.3% by the DNA survey (Figure 3). This large bias reinforced the argument that the RNA survey is a more reliable method than the DNA survey to investigate the eukaryotic community structure Numbers are the correlation coefficients, background colors (from dark red to white) indicate significance levels associated to the p values (0.001, 0.01, 0.05) and no significance. (Massana et al., 2015). The RNA-based results showing the abundant MAST reflected active functions of the microbial food web in the NSCS.
In contrast, the DNA and RNA surveys reported roughly the same community structure for nanoeukaryotes in this study (Figure 3). However, we should not interpret this consistency as an indication of that the DNA survey did not obscure the results with artifacts. The consistency between the DNA and RNA samples happened at the group level where Dinophyceae dominated supremely. But differences emerged at the OTU level, for example, A. sanguinea (OTU2) at the DCM layer of station 2A1 occupied 68.2% of RNA sequences relative to 32.7% of DNA sequences. The findings emphasized again the necessity to conduct the RNA survey to tackle with the technical issues of the DNA survey (Massana et al., 2015).
Filters may collect material belonging to species from other size fractions. We found many pairs of OTUs with identical sequences between the nano-and pico-eukaryotes. This phenomenon has long been an issue for analyzing plankton samples (de Vargas et al., 2015;Massana et al., 2015). de Vargas et al. (2015) reported that 64% of all OTUs spanned two or more size classes and Massana et al. (2015) even found 32-55% of picoeukaryotic OTUs and 20-41% of nanoeukaryotic OTUs in filtrate passed through a 0.8 µm-pore-size membrane filter. We found a significant proportion of shared OTUs were MALV-II and MALV-I (also known as Syndiniales). Syndiniales are marine parasites that infect hosts from protists to fish, they could be detected in the pico-sized fraction because they produce small, flagellated cells (dinospores) to infect new hosts (Guillou et al., 2008;Siano et al., 2011). Moreover, identical sequences do not essentially indicate identical species because OTUs are clusters. In spite of deficiencies, filtration is the most feasible way to study diversity for size characteristics of plankton, and the method did reveal distinctness of size classes in terms of relative abundance of OTUs (Table 1 and Figure 3). Besides, BLAST identifies sequences with the tophit taxon in the reference database which cannot guarantee the most accurate result because in fact there always have a series of hits with good scores; although BLAST is the common method, more robust methods for taxonomic annotation shall be further considered.

Differences in Mechanisms Shaping the Community Structure of Nano-and Pico-Eukaryotes in the NSCS
Environmental filtering played an important role in shaping the community structure of both nano-and pico-eukaryotes ( Table 2). Phytoplankton growth in the NSCS is limited by nitrogen in the open water and phosphate and silicate in the river water influenced area (Yin et al., 2000(Yin et al., , 2001. During the season when the study was conducted, the monsoon winds were too weak to break down water-column stratification, which prevents nutrients reentering the surface. Then, nutrients at the surface layer were depleted by light-saturated photosynthesis while nutrients at the DCM layer were abundant due to light limitation on nutrient uptake. Therefore, the community structure at the DCM layer was distinct from that at the surface layer because of the combined effects of photoadaptation and elevated nutrients, and decreased water temperature might play a role as well. Thus, it was not surprising that depth was the governing factor of the beta diversity (Countway et al., 2010;Wu et al., 2017a;Giner et al., 2020). However, most of the variation in community structure were explained by neutral processes (Figures 5, 6), which was consistent with previous studies (Sloan et al., 2006). Furthermore, the NST ruzicka was higher for nanoeukaryotes than for picoeukaryotes, and the VPA results also indicated that nanoeukaryotes were subject to dispersal limitation while picoeukaryotes were not ( Table 2). The latter was consistent with a previous study which suggested increasing dispersal limitation with cell size (de Vargas et al., 2015). Neutral community theory explains the species of a local community as a sample of the meta community regulated by immigration and extinction at the local level, therefore dissimilarities between local communities will be a function of the geographical distance (namely the dispersal limitation). The NCM displayed that the immigration rate for nanoeukaryotes was no less than that for picoeukaryotes ( Figure 5). However, the number of individuals for the same amount of carbon biomass decreases with cell size, then the possibility of extinction at the local level will be higher for nanoeukaryotes than for picoeukaryotes, then the dispersal limitation might be explained. Ecological drift might be another mechanism accounting for the diversity of nanoeukaryotic community in this study. Ecological drift is stochastic changes in population resulted from equal chances of individuals to reproduce or die. Ecological drift drifts community composition from expectations of environmental filtering (Vellend, 2010). In this study, Dinophycease dominated the nanoeukaryotic communities by contributing two-thirds of the relative read abundance (Figures 3A,B), which might limit other taxa from establishing their population, leading to high similarity between local nanoeukaryotic communities.
The Kuroshio intrusion is supposed to mainly affect the eastern NSCS. However, Figure 1B shows the geostrophic currents in the NSCS in early May 2016 (a week before the cruise), which flowed into the NSCS through the northern Luzon Strait, and generated a westward branch at approximately 21 • N. The current reached approximately 114 • E and 20 • N. Therefore, substantial fractions of Kuroshio water were found in the middle of the NSCS (Figure 1A). The VPA analysis showed that variations in the nanoeukaryotic community structure could be partially explained by the KI index. The reason could be different activities of archaea and bacteria that was reported by Xu et al. (2018). That study showed distinctly higher ammonia oxidation rates by archaea and bacteria in the MKI stations than in the LKI stations. Ammonia-oxidizing archaea and bacteria were dependent on dissolved organic matter, which was proportional to the faction of Kuroshio water.

Relative Metabolic Activities in Relation to Community Structure
Some argued that the discrepancy between the DNA and RNA surveys was mainly due to high rDNA copy numbers of some species or dead cells and did not reflect the metabolic activities, which, however, could be assessed by the relative changes of the RNA: DNA ratio between sampling sites (Hu et al., 2016). For example, a previous study in the northwestern Pacific Ocean reported different relative metabolic activities of Mamiellophyceae for different nutrient regimes (Wang et al., 2019). In this study, Mamiellophyceae also showed lower RNA: DNA ratios with increased nutrients (Figure 8B). Dinophyceae and MAST showed significantly higher RNA: DNA ratios at the DCM layers than at the surface layers (Figure 7B), which was consistent with a previous study in the South China Sea (Xu et al., 2017). That study reported that Dinophyceae, MALV-I, and MALV-II had the highest relative metabolic activities in middle layers (200, 300, and 500 m), followed by deep layers (1,000, 1,500, 2,000, and 3,900 m) and shallow layers (5, 25, and 75 m). Our analysis indicated that the relative metabolic activities of Dinophyceae and MAST increased with nutrients and the Chl a concentration (Figure 8).
In contrast, the correlations between the metabolic activities of nanoeukaryotes and environmental factors were not apparent (Figure 8). Nanophytoplankton were expected to have lower nutrient affinity than picophytoplankton, therefore nanophytoplankton would show larger changes of the relative metabolic activities than picophytoplankton for the ranges of nutrients in this study. However, our results contradicted this expectation. The metabolic activities of nanoeukaryotes might change parallelly, thus the changes were not reflected by the relative values. The relatively stable RNA: DNA ratios might account partially for the dominance of nanoeukaryotic Dinophyceae at all sampling sites ( Figure 7B), but the underlying mechanism was not clear.

CONCLUSION
In the NSCS during the study period, Dinophyceae was the most diverse and sequenced nanoeukaryotic group revealed by both DNA and RNA surveys. In contrast, the relative read abundance of MAST, Pelagophyceae, and Dinophyceae in the size fraction of picoeukaryotes might be largely underestimated by the DNA survey. Environmental filtering played an important role in shaping the community structure, and the sampling depth became the governing factor of the beta diversity under the environmental setting of stratification during the study period. The spatial variations in the diversity of nanoeukaryotes were subject to the dispersal limitation under the size rule. The effects of the Kuroshio intrusion on the nanoeukaryotic community structure might also be explained by the dispersal limitation. Overall, neutral processes are critical in shaping the community structure of nanoeukaryotes. The relative metabolic activities of nanoeukaryotes were relatively stable in accordance with the high similarity of community structure between sampling sites. The responses of the relative metabolic activities of picoeukaryotes to environmental factors displayed two distinct patterns. Our findings improved the understanding of the nano-and picoeukaryotic communities in the NSCS and the mechanisms of their assembly.

DATA AVAILABILITY STATEMENT
The data presented in the study are deposited in the GenBank Sequence Read Archive repository, accession number PRJNA684778, https://www.ncbi.nlm.nih.gov/.

AUTHOR CONTRIBUTIONS
FW, YX, and BH designed the experiments. FW carried out the experiments and wrote the manuscript. SC helped analyze the data. YX and BH provided critical feedback and revised the manuscript. XW and JM helped analyze the data and revised the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by grants from the National Natural Science Foundation of China (Nos. U1805241 and 41776146), Monitoring and Protection of Ecology and Environment in the East Pacific Ocean (No. DY135-E2-5-5), and the study was also partially supported by NSFC (No. 41890803) and Xiamen University Presidential Fund (20720180102).

ACKNOWLEDGMENTS
We thank JianYu Hu for providing the CTD data and Minhan Dai for proving the nutrient data. The nutrients concentration data could be acquired from the supporting information of Xu et al. (2018). This study also benefited from the discussion with Wenxue Wu at Sun Yat-sen University, and Lizhen Lin, Jixin Chen, and Heng Zhu at Xiamen University for their assistances in phytoplankton sample collections and analyses. We thank NSFC Shiptime Sharing Project for supporting the cruises by R/V Dongfanghong 2 (No. NORC2015-05), and we also acknowledge captains and crew for their cooperation during the cruise. We truly thank two reviewers whose comments helped to improve and clarify this manuscript.