Diversity and biogeography of picoplankton communities from the Straits of Malacca to the South China Sea

The study area extends from the Straits of Malacca to the South China Sea (1°N–21°N, 96°E–114°E). Surface seawater samples were collected from 12 locations during a multidisciplinary cruise carried out in the Eastern Indian Ocean by R/V Shiyan 3, the South China Sea Institute of Oceanology, Chinese Academy of Sciences, from 14 to 19 April 2017 (Fig. 1). The surface water for phylogenetic analysis (2 l) was filtered through a polycarbonate filter with a pore size of 0.22 μm (Millipore Isopore, Bedford, USA) after pre-filtration through a polycarbonate filter with a pore size of 3 μm (Millipore Durapore, Bedford, USA). The filters (0.22 μm pore size) were snap frozen in liquid nitrogen and then stored at −80°C in the laboratory until DNA extraction. Seawater for nutrient analysis was filtered through polycarbonate filters with a pore size of 0.45 μm, collected in acid-washed polyethylene containers and immediately frozen (−20°C) until analysis.

Figure 1

Temperature, salinity and dissolved oxygen (DO) were measured at the sampling locations with a Sea-Bird SBE9 Conductivity-Temperature-Depth (CTD) profiler (USA). Chemical characteristics, including concentrations of nitrate (NO₃-N), nitrite (NO₂-N), ammonium (NH₄-N), silicate (SiO₃-Si) and phosphorus (PO₄-P) were analyzed according to Jiang et al. (2014, 2015).

Total genomic DNA was extracted from the filters using the PowerWater DNA Isolation Kit (MoBio, Carlsbad, CA, United States) according to the manufacturer’s protocol. The DNA was stored at −80°C for subsequent analyses. The V4-V5 region of the 16S rRNA gene with the 515f/907r primer set (Xiong et al. 2014) and the V9 region of the 18S rRNA gene with the 1380f/1510r primer set (Amaral-Zettler et al. 2009) were selected for targeting amplicons. Next generation sequencing library preparations and Illumina MiSeq PE250 sequencing were conducted at the Novogene Company (Beijing, China). All sequences obtained in this study were deposited in the NCBI Sequence Read Archive (SRA) with accession numbers PRJNA492365 for the 16S rRNA gene and PRJNA492462 for the 18S rRNA gene.

The software package QIIME 1.9.1 (Caporaso et al. 2010) was used for 16S and 18S rRNA gene data analysis. The forward and reverse reads were derived from the original DNA fragments, which were merged by using FLASH (Magoč & Salzberg 2011) and assigned to a sample based on a barcode and truncated by cutting off the barcode and primer sequence. Sequences were grouped into operational taxonomic units (OTUs) using the clustering program VSEARCH 1.9.6 (Edgar 2010) at 97% sequence identity. Rarefaction analysis was conducted using the original detected OTUs. The taxonomic assignment was performed by the RDP classifier at a confidence threshold of 0.8 (Wang et al. 2007).

The Shannon index, the Chao1 index and Good’s coverage of samples were determined as described by Schloss et al. 2009. Weighted and unweighted UniFrac, Bray–Curtis and principal coordinate analysis (PCoA) were performed by QIIME. The unweighted pair group method with arithmetic mean (UPGMA) clustering was conducted by unweighted and weighted UniFrac based on the protocol published by Caporaso et al. 2010.

Heatmap and redundancy analysis (RDA) was conducted using the R statistical package to determine the correlation between community composition and environmental parameters (Zhang et al. 2014).

Network analysis was performed according to Qin et al. (2010). Forty most abundant OTUs of picoplankton communities were selected and Spearman correlations were calculated. In order to reduce the network complexity, correlations r below −0.60 or above 0.60 (n =162 and n=128 for 16S and 18S rRNA gene data, respectively) were visualized in a graph using Cytoscape 2.8.2 (Shannon et al. 2003).

After quality filtering of the raw reads, a total of 901199 16S rRNA gene sequences (on average 75100 reads per sample) and 1 238 847 18S rRNA gene sequences (on average 103237 reads per sample) from 12 samples were obtained and clustered into 1914 OTUs and 3931 OTUs for 16S and 18S rRNA genes (97% cutoff), respectively. The Shannon index indicated that eukaryotic picoplankton a-diversity was significantly higher (independent t-test and paired t-test; p < 0.05) than that of the prokaryotic picoplankton (Table 1).

In total, 34 prokaryotic phyla and 55 eukaryotic phyla were detected. The top five bacterial phyla accounted for more than 97.7% of the total abundance. Proteobacteria (40.4%; a – 31.8%, δ – 4.8%) and Cyanobacteria (39.2%; Prochlorococcus – 23.1%, Synechococcus – 15.5%) were the most abundant prokaryotic group, followed by Actinobacteria (8.1%), Bacteroidetes (6.8%), Firmicutes (1.9%), SAR406 (1.4%; Fig. 2a). At sampling sites S3 and S5, Cyanobacteria reached the highest relative abundance, while the relative abundance of Proteobacteria was the lowest among the twelve sampling sites. At sampling sites S3 and S5, Cyanobacteria reached a higher relative abundance, while Proteobacteria were relatively less abundant among the twelve sampling sites. The top 10 eukaryotic phyla were: Dinoflagellata (8.9%), Cnidaria (8.2%), Retaria (8.1%), Tunicata (8.1%), Arthropoda (7.8%), Protalveolata (4.1%), Ciliophora (3.9%), Ochrophyta (2.4%), Prymnesiophyceae (1.5%) and MAST-3 (1.3%). Cnidaria and Retaria were the least abundant at site S3, whereas Tunicata reached the highest abundance at site S9.

Figure 2

The results of the principal coordinate analysis (PCoA) based on weighted UniFrac distances revealed that most differences could be attributed to geographical location (Fig. 3). The first two PCoA axes explained more variance in the abundant prokaryotic communities (84.43%; Fig. 3a) than in the eukaryotic communities (42.36%; Fig. 3b). Samples from the SCS and site S1 were grouped together (upper side and right side of the PCoA plots, respectively), while samples from four other sampling sites of the SM were separated.

Figure 3

Based on the composition of the top 40 OTUs, samples were clustered into two groups (Fig. 4). The top 40 OTUs of the 16S RNA gene at sampling sites S3, S4 and S5 were clustered together into a subgroup and the remaining OTUs were included in another group (Fig. 4a). Cyanobacteria, especially Synechococcus (OTU4, OTU25, OTU59, OTU899 and OTU1805) were abundant in the SM. Different OTUs, even those associated with the same classification, showed different distribution patterns. Within Proteobacteria, the relative abundance of OTU2 and OTU1494 was higher in the SCS, but the relative abundance of OTU19, OTU28 and OTU29 was high at site S2 located in the SM. Sampling sites S4 and S5 were in the same cluster, other sampling sites were in the same group based on the top 40 OTUs of the 18S RNA gene (Fig. 4b). Each site featured several highly abundant OTUs. The relative abundance of OTU7, OTU11, OTU12, OTU13, OTU15 and OTU19 was higher in the SM than in the SCS. The relative abundance of OTU22, OTU25 and OTU35 was higher at site S7.

Figure 4

Redundancy analysis (RDA) was performed to identify possible correlations between environmental factors and picoplankton distribution (Fig. 5). Each environmental variable in the RDA biplot was represented by an arrow and the length of an individual arrow indicated how much variance was explained by that variable. The RDA analysis based on the prokaryotic picoplankton community composition was consistent with PCoA (Fig. 3a). Sampling sites S2, S3, S4 and S5 are located on the right-hand side of the RDA plot. However, the RDA models for eukaryotic picoplankton assemblages indicated that temperature and DO (statistically significant environmental factors; p < 0.05) only partly accounted for the diversification of communities from the SM to the SCS (Fig. 3b). In the two RDA models of prokaryotic picoplankton and eukaryotic picoplankton, the environmental factors explained about 96.6 and 61.6% of the total variance in the community composition, respectively.

Figure 5

The RDA analysis revealed that salinity (p < 0.01), temperature (p < 0.1) and NO₂-N (p < 0.1) were statistically the most significant variables, explaining the pattern of the prokaryotic picoplankton community composition (Fig. 5a). However, temperature and DO concentrations were significant factors determining the composition of the eukaryotic picoplankton community (p<0.05) (Fig. 5b).

As the whole networks of the total OTUs are too complicated to display, only a very limited number of key OTUs with higher abundance were considered. The top 40 OTUs were examined (Fig. 6). All curves of the network connectivity fitted to power-law model (r < −0.60 or r > 0.60). There are similarities in terms of network size and structure between 16S and 18S rRNA-based sequences. Prokaryotic taxa displayed closer interactions than the eukaryotic community (162 and 128 connections for 16S and 18S rRNA-based sequences, respectively). Positive connections (97/162 and 72/128 for 16S and 18S rRNA-based sequences, respectively) prevail in both prokaryotic and eukaryotic communities.

Figure 6

In the present study, the Illumina MiSeq sequencing technology was employed to describe the diversity and biogeography of picoplankton (including prokaryote and eukaryote) communities from the SM to the SCS. The comparison of the community structure based on OTU percentages indicated a divergent distribution pattern of individual picoplanktic taxa.

The results showed that Proteobacteria and Cyanobacteria were the dominant phyla of prokaryotic picoplankton and Alphaproteobacteria accounted for more than 75% of Proteobacteria in the study area. According to previous studies and based on phylogenetic and phenotypic properties as well as strong adaptability to various marine environments (Sunagawa et al. 2015; Shi et al. 2018), Proteobacteria were identified as the most abundant bacterial phylum in all marine water masses (Suh et al. 2014; Sunagawa et al. 2015; Sun et al. 2015; Li et al. 2017; Shi et al. 2018).

Being important autotrophic picoprokaryotes, Cyanobacteria not only regulate the global primary productivity, but also drive carbon and nitrogen cycles in the marine ecosystem ( Xia et al. 2015). Synechococcus and Prochlorococcus are the most abundant groups of picophytoplankton in the world’s oceans and they have different ecological niches (Zwirglmaier et al. 2008; Xia et al. 2015). Prochlorococcus is the dominant genus in oligotrophic areas, while Synechococcus inhabits mainly coastal mesotrophic and eutrophic environments, such as coastal waters and upwelling regions (Zwirglmaier et al. 2008; Xia et al. 2015). This study revealed that within Cyanobacteria, the abundance of Prochlorococcus in the open sea was significantly higher than that of Synechococcus (p < 0.05), while in the strait area the abundance of Prochlorococcus was significantly lower than that of Synechococcus (p < 0.05). Chen et al. (2011) mentioned that picophytoplankton might be restricted by high concentrations of nutrients and (or) heavy metals in the estuary and coastal areas of the northern SCS. In general, the open sea is usually characterized by oligotrophic conditions, while the strait area is highly eutrophicated by nutrients from terrestrial input and transportation. Therefore, the distribution of Synechococcus and Prochlorococcus showed varying characteristics from the SM to the SCS.

Five phyla – Dinoflagellata, Cnidaria, Retaria, Tunicata and Arthropoda – dominated in eukaryotic picoplankton from the SM to the SCS. We distinguished ten most abundant eukaryotic phyla and the results showed that only three of them (Dinoflagellata, Ochrophyta and Prymnesiophyceae) belonged to phytoplankton, while the remaining seven phyla were affiliated with zooplankton. Dinoflagellates are the major group of unicellular protists inhabiting both marine and freshwater environments (Price & Bhattacharya 2017). These aquatic protists account for the main part of the marine plankton and play an important role in the marine primary productivity, while being represented by some toxic species (Price & Bhattacharya 2017). They form the basis of food chains and are a large component of the food web in the marine ecosystem (Taylor et al. 2008; Sze et al. 2018). High diversity of dinoflagellates was detected by DNA metabarcoding using high-throughput sequencing in Singapore’s waters (Sze et al. 2018). Annett et al. (2010) reported that both heterotrophic and autotrophic dinoflagellates contributed as much as 27% to the total phytoplankton biomass in Ryder Bay. The phylum Cnidaria, containing over 10000 species of animals, is usually classified into five classes: Anthozoa, Cubozoa, Hydrozoa, Scyphozoa and Staurozoa (Kayal et al. 2013). They occur only in aquatic environments, mostly in marine ecosystems, and are ecologically important marine invertebrates (Kayal et al. 2013). Cnidarians are found from deep waters near hydrothermal vents to the polar seabed and tropical reefs (Rocha et al. 2015). The SCS features an extraordinary richness of reef corals (571 known species; Huang et al. 2015), due to the large number of Anthozoa (Cnidaria).

In the present study, the weighted UniFrac PCoA indicated that picoplankton communities of prokaryotes and eukaryotes from the SM and the SCS belong to different groups. Samples from the SCS and site S1 were grouped together and placed on the opposite side of the four other sampling sites in the PCoA plots. Jiang et al. (2015) mentioned that geographical proximity is an important factor affecting the structure of phytoplankton communities from the Pearl River estuary to the SCS. Although site S1 was located in the SM mouth facing the Indian Ocean, the picoplankton communities exhibited characteristics of the open sea, such as the SCS in this study. It has been reported that sampling sites located near islands had a similar bacterial community (Ling et al. 2012). The distance between land and a sampling location may play an important role in the distribution of eukaryotic ultraplankton (Jiang et al. 2014). Therefore, the geographical location could influence the structure of picoplankton communities from the SM to the SCS.

RDA was conducted to further explore the relationships between the environmental factors and the community structure of prokaryotic and eukaryotic picoplankton. With regard to prokaryotes, the distribution of the sampling sites in the UniFrac PCoA analysis was consistent with the RDA results. Prokaryotic samples were clustered into two groups: the strait area (sampling sites S2–S5) and open waters (sampling sites S1, S6–S12). Temperature (p < 0.01) as well as NO₂-N (p < 0.1) and NO₃-N concentrations could affect the prokaryotic community structure in the SM, while salinity significantly affected it in the coastal and open waters. The relationship between the prokaryotic plankton community and environmental parameters varies in different waters (Sun et al. 2015; Suh et al. 2014; Shi et al. 2017; Halsey et al. 2017; Cram et al. 2015; Xia et al. 2015). Previous studies showed that the concentration of NO₃-N is not necessarily a limiting factor for the growth of bacteria in estuaries and coastal waters, but in oligotrophic waters it could be a major limiting factor affecting the bacterioplankton community due to its low concentration (Xia et al. 2015). In this study, sampling sites S2, S3, S4 and S5 located in the MS maintained a relatively higher temperature than other sampling sites due to anthropogenic influence and received a large amount of inorganic and organic nutrients originating from anthropogenic activities, whereas salinity was the key factor shaping the prokaryotic picoplankton communities at other sampling sites located in the open sea.

In terms of eukaryotic picoplankton, there was a discrepancy between UniFrac PCoA and RDA. This indicated that other variables, in addition to the eight selected environmental factors, may also contribute to the community cluster patterns, such as grazing. Further research is needed to identify specific factors. Environmental factors (temperature, DO, salinity, grazing and nutrient concentrations) can significantly affect the eukaryotic plankton community in the natural environment (Bernardi Aubry et al. 2013; Jiang et al. 2014; Suikkanen et al. 2007; Suzuki et al. 1997; Wu et al. 2014; De Vargas et al. 2015). Wu et al. (2014) reported that temperature and irradiance affected the picoeukaryotic distribution at the surface and 60 m sampling depth in the SCS (p < 0.001). Temperature, irradiance, nutrient concentrations and salinity were mainly correlated with changes in the main phytoplankton groups and the phytoplankton community composition (Xiao et al. 2018). Jiang et al. (2014) demonstrated that temperature, salinity, phosphorus and silicate had a significant impact on the community structure of eukaryotic ultraplankton of the northern SCS. This study indicated that temperature and DO played an important role in shaping the eukaryotic picoplankton community composition from the SM to the SCS (p < 0.05). Picophytoplankton (including Synechococcus, Prochlorococcus and picoeukaryotes) are usually dominant in waters characterized by high temperature and low content of nutrients (Agawin et al. 2000; Li 2002; Finkel et al. 2010; Moran et al. 2010; Chen et al. 2014). Temperature and DO affect the growth rate of picoplankton (Chen et al. 2014; Wu et al. 2014). Therefore, both prokaryotic and eukaryotic picoplankton communities were significantly influenced by temperature (p<0.01 or p <0.05).

Complicated networks were constructed by a variety of species interacting with each other in the complex marine ecosystem (Montoya et al. 2006). The marine ecosystem performs the system-level functions (such as biogeochemical cycling, ecosystem stability) owing to these complicated network interactions, and the functions could not be fulfilled by individual populations (Zhou et al. 2010). Community network models were built to demonstrate the interactions between species. The determined competition or cooperation relationships between species were based on nutrients, space, material and information (Zhang et al. 2014). Although ecological networks are very important to the ecosystem, there are few researches on networks of the picoplankton community. Based on Fig. 6, we found that the network size and structure of the prokaryotic community was significantly different from that of the eukaryotic community. In general, the interaction of bacterioplankton was relatively close compared to eukaryotic picoplankton. The results also showed that the cooperation dominated the relationships of both prokaryotic and eukaryotic communities. This is in line with the previous study (Zhang et al. 2014), which has reported that positive connections dominated the interactions between taxa in the DNA-based networks of the total bacterioplankton in the SCS.