Zooplankton sampling was carried out from January 2013 to December 2013 on a twice per month schedule, over two stations situated on the continental shelf off the Saudi Arabian coast in the central Red Sea (offshore station: N22°31′21.08″, E38°99′48.54″; inshore station: N22°18′46.08″, E38°55′41.70″; Figure 1). The offshore station was located in an open deep water area with a depth of over 500 m, while the inshore station was located in a coastal reef-habitat with a depth of 60 m. Samples from several depth strata were collected along two transects parallel to the coast by using MANTA, Bongo and WP2 nets. The MANTA trawl with an inlet aperture of 50 × 15 cm, was fitted with a net of mesh size 150 μm. The Bongo and WP2 nets, with a 60 and 50 cm diameter aperture respectively, were fitted with nets of mesh size 150 μm with 50 μm cod end and 200 μm with a 50 μm collector respectively. The MANTA net was trawled for 10 min at 2 knots only at the offshore station to collect surface samples. The Bongo net was towed obliquely from 50 m depth to the surface and horizontally at 20 m depth for 10 min both offshore and inshore. The WP2 was hauled vertically to the surface from 100 and 50 m depth respectively at the offshore and inshore stations. Flow rate was monitored using a digital flow meter (model 23.090, KC Denmark, Silkeborg, Denmark) mounted on the mouth aperture of the samplers, while the hydrological parameters were simultaneously recorded using a CTD profiler (Idronaut Ocean Seven 320Plus) mounted 1 m above the net. The mean volume of water filtered was calculated following manufacturer's instructions. All samples were preserved in 98% ethanol (Sigma-Aldrich, St. Louis, MO) and stored at 4°C.

Samples collected using the MANTA and bongo nets were filtered through Whatman Folded filter paper (grade 113v, diam. 24 cm) to eliminate the ethanol. The weight was recorded using a digital weighing scale and the total biomass, defined as the animal density of total zooplankton, was calculated by subtracting the wet filter weight. Biomass of the samples collected using the vertical WP2 net was not estimated.

Aliquots of 5–10% of the total volume were used for the genomic analysis. Samples were digested in ATL lysis buffer (Qiagen, Valencia, CA) with 0.5 mg/ml proteinase K at 56°C overnight. Genomic DNA was isolated using the standard phenol-chloroform method (Sambrook and Russell, 2001).

Amplification was performed using the general eukaryotic primer set targeting the hypervariable region of the 18S rRNA gene designed by Amaral-Zettler et al. (2009). The forward primer 1389F was used in combination with the reverse primer 1510R. This set of primers was developed as a standard marker of marine microbial eukaryotic diversity and has also shown reliable amplification of marine zooplankton assemblages (Pearman et al., 2014; de Vargas et al., 2015; Pearman and Irigoien, 2015; Albaina et al., 2016). All primers had a Roche 454 tag and a 10 bp barcode. PCR reactions were undertaken in 50 μl reaction volumes containing 2.5 U Taq polymerase (Invitrogen, Carlsbad, CA), 1X Taq reaction Buffer, 200 μM dNTPs (Invitrogen), 1.5 mM MgSO4, 0.05 mg Bovine Serum Albumin (BSA) and 0.2 μM of each primer. PCR conditions were as follows: an initial denaturation step of 3 min at 94°C, followed by 30 cycles of 94°C for 30 s, 57°C for 45 s and 72°C for 45 s and a final extension at 72°C for 10 min. All the samples were run in duplicate alongside a negative no template control. PCR products were combined and checked on a 1% agarose gel; the band was then excised and purified with Qiagen's gel extraction kit (Qiagen, Valencia, CA). Using a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA) samples were equalized to 50 ng and pooled according to barcode. A maximum of 14 barcodes were pooled together and purity of multiplexed samples was analyzed using an Agilent 2100 bioanalyzer machine (Agilent, Santa Clara, CA).

Library construction was performed using the Roche XLR70 kit and each pooled sample was sequenced on a ¼ run either on a Roche GS FLX system at the KAUST core facility or a Roche GS FLX+ system at IMGM Laboratories (Martinsried, Germany). The sff files generated during sequencing were deposited in the National Center for Biotechnological Information (NCBI) Short Read Archive (SRA) under the study accession number SRP076465.

Sequences were firstly demultiplexed based on the barcode and raw reads were filtered based on quality (q = 25 over a window of 25) and length (minlenth = 100, maxlength = 200) within the QIIME (Caporaso et al., 2010) framework. No ambiguities in bases were allowed and only a single mismatch in the forward primer, whilst a homopolymer length up to 6 was permitted. Sequences were subsequently trimmed to the reverse primer. Clustering into operational taxonomic units (OTUs) was undertaken using a two-step clustering process as described in Pearman and Irigoien (2015). Briefly, the trie function within QIIME was used with CD-HIT (Li and Godzik, 2006) as a primarily clustering step before de novo clustering using USEARCH (Edgar, 2010) at a 97% similarity threshold. Representative sequences were obtained and chimera checked against the SILVA 119 database (Pruesse et al., 2007) using UCHIIME (Edgar et al., 2011).

The resulting chimera cleaned reference sequences were taxonomically assigned using the naïve Bayesian classifier rdp (Wang et al., 2007) within QIIME with the SILVA 119 database as reference. To compare the composition between both stations, reads were rarefied multiple times (n = 100) at an even depth (3,000 reads) and the average used (any samples which did not meet this criteria were removed from further analysis).

Taxonomic compositions were constructed using the R package phyloseq (McMurdie and Holmes, 2013) and plotted using ggplot (Wickham, 2009). Non-metric multidimensional scaling (NMDS) was used to graphically summarize and analyze seasonal changes in the structure of zooplankton communities in the Central Red Sea. Graphs were plotted in R using phyloseq and statistical analysis of the distance matrices was undertaken using Analysis of Similarity (ANOSIM) in vegan (Oksanen et al., 2013). Contour and density plots were produced in R using the lattice package.

Canonical Correspondence Analysis (CCA) was used to elucidate the role played by different environmental variables in modulating the structure and seasonal succession of zooplankton assemblages.

The data sets used consisted of bimonthly species abundance represented by number of reads and nine environmental parameters: Sea Surface Temperature and Salinity (SST and SSS), Mean Temperature and Salinity of the water column for the whole CTD cast (T.mean and S.mean), Mean Temperature and Salinity of the water column from 0 to 50 m (T.50 and S.50) (for inshore stations is the same as T.mean and S.mean), Temperature and Salinity at 20 m (T.20 and S20) and UV index (UVI). The significance of these variables was assessed using Monte Carlo permutation tests (with 999 unrestricted permutations). The species data were log (x+1) transformed and rarely occurring taxa (less than 0.5% of the total abundance) were down-weighted in order to prevent them from greatly influencing the analyses (Ter Braak and Šmilauer, 2002).

The mean volume of water filtered by the MANTA, Bongo horizontal tow, Bongo oblique tow and vertical WP2 nets was 70.4 m³, 236/210 m³ (offshore/inshore), 252/200 m³ (offshore/inshore) and 28/15 m³ (offshore/inshore), respectively (Table S1).

The temperatures ranged between 22.83–31.78°C at the offshore station and 25.31–31.36°C inshore. Salinities oscillated between 39.4–40.5 and 39.3–39.8 psu offshore and inshore respectively (Figure 2).

A total of 3,407,473 sequenced reads were obtained from all 154 samples. Of these, 2,033,939 reads passed quality checks, including the removal of chimeras (for further details, see Table S1). The quality-filtered sequenced data was clustered at a 97% similarity level and after rarefaction a total of 630 metazoan OTUs were obtained. Overall the most diverse and dominant phylogenetic group in terms of OTUs was the Arthropoda accounting for 417 OTUs (66%) followed by Cnidaria (67 OTUs), which also contributed more than 10% of the overall number of OTUs. The Arthropoda were comprised of 174 OTUs attributed to Malacostraca and 145 OTUs belonging to the class Maxillopoda, while the majority of the remaining OTUs could not be assigned at class level.

Similar numbers of OTUs were observed in the offshore station (558 OTUs) (when the MANTA net samples were removed) compared with the inshore station (572 OTUs). There was a significant level of overlap in OTU composition between both sites with a total of 507 OTUs observed at both stations. Comparison of sampling methods showed that the MANTA net revealed the lowest number of OTUs (288 OTUs) in the offshore station followed by the Bongo horizontal net (371 OTUs) which accounted for the lowest number of OTUs (389 OTUs) in the inshore station (where the MANTA net was not used). Inshore, 41% of the OTUs (239) were shared between all sampling methods whilst only 28% (162) were shared between the 4 sampling methods in the offshore station. Removing the MANTA net samples revealed that 45% of the OTUs (251) were shared between the three sampling methods utilized in the inshore station. In the inshore station the vertical net had the most unique OTUs (77 OTUs) while the oblique net revealed the most unique OTUs in the offshore station (45 OTUs).

In order to assess statistical differences in the zooplankton communities under study, an Adonis analysis was performed to investigate the effects of the net and the location. Significant interactions were observed between the factors (Net and Location), indicating that there were variations in the trends for both the weighted (F = 1.67; p = 0.03) and unweighted (F = 1.41; p = 0.018) UniFrac analyses (Figure 6).

To investigate the effects of environmental variables (temperature, salinity, UV index) on zooplankton seasonal composition in the central Red Sea, we measured changes in community assemblages for each of the sampling methods separately (MANTA, Bongo horizontal, Bongo oblique, and WP2 vertical) using Canonical Correspondence Analysis (CCA) in R. The CCA analysis revealed clear differences between the cold and hot season communities, reinforcing the results described above. Temperature was found to be the main factor shaping the seasonal structure of the zooplankton communities.

The first CCA axis clearly separated the cold and hot season communities for all sampling methods, showing the highest correlation with temperature followed by salinity (Figure 7). Thus it appears that temperature was the most important factor in determining the structure of the zooplankton communities over the seasonal 1 year-period of the present study. Both the cold and hot season communities were rather scattered along the second axis, which only separated inshore from offshore community composition for the horizontal samples.

For the MANTA samples (Figure 7A), the Monte Carlo permutation test showed that only temperature at 20 m (T20) was significant (p < 0.05) in explaining zooplankton seasonal distribution, being correlated with axis 1 (0.73). Only the first canonical axis was significant (p < 0.01), explaining 56.2% of the variance of the species—environment relationship. The species—environment correlation coefficient for the first CCA axis was 0.76.

For the samples collected with the bongos towed horizontally (Figure 7B), the Monte Carlo permutation test showed that temperature at 20 m (T20) and sea surface salinity (SSS) were the only significant (p < 0.05) environmental variables in explaining zooplankton spatio-temporal distribution. Both were positively correlated with axis 1 (0.63), but T20 was negatively correlated with axis 2 (−0.37) and SSS positively (0.37). Only the first 2 canonical axis were significant (p < 0.01), explaining 33.5% of the variance of the species—environment relationship. The species—environment correlation coefficients for the first CCA axes were 0.82 and 0.66 respectively.

The Monte Carlo permutation test of the samples collected with the bongo nets obliquely towed from 50 m depth to the surface (Figure 7C), revealed that only T20 was significant (p < 0.01) in explaining zooplankton seasonal distribution. T20 was positively correlated with axis 1 (0.79). Only the first canonical axis was significant (p < 0.01), explaining 40% of the variance, and the species—environment correlation coefficient for this axis was 0.82.

Finally, for the samples collected using the WP2 vertical net (Figure 7D), the Monte Carlo permutation test showed that only 2 environmental variables (T20 and S50) were significant (p < 0.05) in explaining zooplankton seasonal distribution. T20 was positively correlated with axis 1 (0.65), whereas S50 was negatively correlated with axis 2 (−0.41). Only the first 2 canonical axes (axis 1 and axis 2) were significant (p < 0.01), explaining 33.3% of the variance, and the species—environment correlation coefficients for these axes were 0.80 and 0.56 respectively.

To summarize, our study reveals two marked seasonal phases in terms of species composition of the zooplankton communities living in the upper water column (top 100 m) in the central Rea Sea, driven mainly by temperature. The cold season was characterized by a lower biodiversity being mainly dominated by classes Maxillopoda and Hydrozoa. Biodiversity increased during the hot season and was defined by the appearance of classes Malacostraca, Cephalopoda and Gastropoda, among the most abundant in terms of number of reads. Table 1 displays the taxonomic designations (at genus level where possible) of the most dominant OTUs for the epipelagic zooplankton community inhabiting the top 100/50 m of the water column (vertical hauls offshore and inshore respectively). Copepods were the dominant component of the Arthropoda accounting for on average 99.7% of the maxillopodan reads at all depth strata, with the majority of the reads corresponding to calanoid copepods, while poecilostomatoids and harpacticoids occurred only rarely (Table 2). In total, 15 families and 32 species of copepods were identified, with 11 families and 17 genera within the order Calanoida. Five calanoid genera dominated, contributing 58.8% of the total number of reads for the total copepod community: Acartia, Canthocalanus, Undinula, Candacia, and Paracalanus. Undinula were found to be present in high abundances all year around at all depth strata studied, while the rest occurred regularly throughout the year, however, dominance differed. All three species recorded within the genus Acartia revealed a remarkably higher abundance in terms of number of reads during the hot season, while three and four different species within the genus Candacia and Paracalanus respectively alternate their prevalence during both seasons. Finally Canthocalanus was more prevalent during the hot season.

A significant body of research has been undertaken for reef fishes and other common sessile organisms like corals or sponges in the Red Sea (e.g., Berumen et al., 2013; Loya et al., 2014; Erpenbeck et al., 2016; Roberts et al., 2016), identifying it as a region of high biodiversity and endemism. However, zooplankton biodiversity remains largely uncharacterized.

This study takes advantage of high throughput sequencing technology, which allows a rapid characterization of whole communities, to assess the seasonal profiles of zooplankton assemblages in two locations of the central Red Sea through an annual cycle. Our results constitute the first baseline data on the richness and seasonality of the central Red Sea zooplankton. However, it is important to note that the results presented here are based on the molecular data available in public reference databases and are therefore limited by them (Carugati et al., 2015). Reference databases are known to be incomplete in terms of the global inventory of species having a representative barcode sequence (Carugati et al., 2015) limiting the ability to classify species at low taxonomic levels and the marker gene may not have the resolution of species level identifications (Dell'Anno et al., 2015). Understudied areas like the central Red Sea, which is known for housing a highly specialized and partly endemic fauna (Sonnewald and El-Sherbiny, 2017) are expected to harbor a number of undescribed species, which are unrepresented in reference databases and thus, will be classified as unknown. This is likely to have an impact on the success of fine level taxonomic assignments and results in higher level classifications as used in the current study. Our study reports ~15% of a total of 630 metazoan OTUs could only be assigned at the phylum level. This highlights the scarcity of genetic data in public databases for Red Sea zooplankton, and even includes highly abundant taxa. Nevertheless, the 18S region used here, due to its relatively conserved nature is considered to allow accurate classification of novel sequences (Bucklin et al., 2016). In addition, the sensitivity of molecular studies is revolutionizing analysis of marine biodiversity since it allows the correct discrimination and classification of cryptic and sibling species, as well as of early individual life stages (eggs and larvae) which constitute important impediments in classical taxonomic studies. Thus, this emerging technique has revealed previously hidden taxonomic richness (Lindeque et al., 2013) as well as the discrimination of rare and cryptic species of zooplankton communities (Zhan et al., 2013; Zhan and MacIsaac, 2015).

Nonetheless it is important to highlight that high throughput amplicon sequencing analysis applied to diversity studies is still in its infancy and presents several limitations in its current stage. So far, it has not shown good agreement with species abundance data from morphological taxonomic analysis (Lindeque et al., 2013) and there is an urgent need for the development of standardized practices that would allow comparisons between different studies (reviewed in Bucklin et al., 2016). However, it has already shown great value in providing baseline data for marine metazoan plankton diversity (e.g., Fonseca et al., 2010; Pearman et al., 2014).

Our study reports changes in total biomass and in species composition during the annual cycle of the zooplankton communities inhabiting the upper 100 m of the water column in the central Rea Sea. The main abiotic factor driving the biodiversity of the zooplankton communities was found to be temperature, based on the CCA analysis of the seasonal structure of the zooplankton communities, which define two distinct seasonal phases (Figures 2, 7). Despite the recurrent paradigm of the aseasonal nature of the tropical environments (Morales-Nin and Panfili, 2005), we show here that the sea surface temperatures in the Red Sea display high seasonal variations and zooplankton assemblages shape their biomass and biodiversity accordingly. Thus, acclimation capacity of the zooplankton assemblages living here might be higher than expected and closer to that of non-tropical species.

The two seasonal phases correspond well with the two seasonal periods of stratification and vertical mixing described in the Red Sea (Calbet et al., 2015). During the hot season, the water column is stratified and the surface layers are depleted of nutrients (Reiss and Hottinger, 1985). In colder months, the thermocline deteriorates, facilitating the enrichment of surface waters with nutrients brought up from the depths due to deep convective mixing (Wolf-Vecht et al., 1992; Lindell and Post, 1995). It is well established that environmental conditions, particularly temperature and food availability, have a crucial effect on zooplankton abundance and biomass (Webber and Roff, 1995; Christou, 1998). Accordingly, the maximum biomass peak was observed in February (Figure 4), with values that tripled (oblique tows) or even quintupled (horizontal tows) those observed when temperatures and stratification reached their maximum (Figure 2). When temperatures started decreasing in September, zooplankton biomass increased again, doubling the values observed during the summer. This annual trend corresponds to the offshore station, an open deep water area with a standard depth over 500 m. Biomass densities at the inshore station exhibited a more heterogeneous pattern throughout the year, displaying a third biomass peak when seawater temperatures reached their maximum. This pattern might be explained by the fact that this coastal reef zone is a shallow area (60 m) where stratification is not observed (Figure 2) or it could be due to advective variability in the sampling area.

Similar annual trends of biomass densities were described for the zooplankton communities living in coastal areas in the Gulf of Aqaba (Khalil and El-Rakman, 1997), although our results contrast sharply with a study carried out in the Northern Rea Sea (El-Sherbiny et al., 2007). These authors reported consistently high densities of the epipelagic zooplankton community during late spring and early summer, attaining maximum values in August. However, another multiyear study in the Gulf of Aqaba revealed that highest and lowest densities of zooplankton were highly dependent on the year (Cornils et al., 2007). Moreover, as shown in this study, even close geographic areas might show differences in zooplankton biomass densities.

Differences in species composition of the zooplankton assemblages between the communities inhabiting the offshore and the inshore station were detected by Adonis analysis being significant for both the weighted and unweighted UniFrac distance matrices. This observation may reflect fine-scale geographic structuring in terms of biodiversity for the samples collected with the bongo nets; however due to the limited geographic scope of the present study, more stations are needed to confirm this pattern.

The neuston community analyzed in this study revealed the lowest number of OTUs of all the sampling methods. The MANTA net samples the plankton inhabiting the upper 0.15 m of the sea, which is a restricted ecological niche occupied by a unique community (reviewed in Marshall and Burchardt, 2005). Accordingly, the lower biodiversity compared to the samples collected using the two other net systems is expected. On the other hand, the samples collected with the vertical hauls showed the highest OTU richness of all samples as expected, since the strata sampled is the widest in this study.

Throughout the year the dominant phyla found in our study were Arthropoda and Cnidaria, at all depth strata studied. The majority of reads consisted of classes Maxillopoda and Hydrozoa, the latter being especially abundant during the cold season, accounting for on average 48.1 and 25.6% of the reads respectively. In the surface (MANTA samples, offshore) these two classes together contributed on average 85.8% to the total number of metazoan reads, while their contribution was on average 56.0, 66.4, 81.6% offshore and 67.2, 73.3, 85.9% inshore for the samples collected using the horizontal (20 m), oblique (50 to surface) and WP2 (100/50 to surface) nets respectively. The rise in the water temperature was accompanied by an increase in biodiversity with the appearance of classes Malacostraca, Cephalopoda, Gastropoda and Saggitoidea, among the most abundant, accounting for on average 17.1, 9.3, 1.6 and 0.7% offshore and 20.1, 2.5, 6.2 and 1.1% inshore during the hot season. Similar groups were found to be the most abundant in previous studies in the Red Sea using both metagenetics (Pearman et al., 2014; Pearman and Irigoien, 2015) and classical morpho-taxonomy approaches (Cornils, 2005; Cornils et al., 2007; El-Sherbiny et al., 2007; Dorgham et al., 2012). Studies in other tropical marine regions have reported similar compositions of zooplankton assemblages (e.g., Spermonde archipelago: Cornils et al., 2010; Pohnpei island: Machida et al., 2009). Metagenetic studies in temperate regions have reported significantly lower numbers of OTUs (the Western Channel: Lindeque et al., 2013; the North Sea: Mohrbeck et al., 2015; the Northwest Pacific: Hirai et al., 2015), however, differences in the sampling methodology (e.g., the plankton net), the barcode, primer and protocol choice and/or the bioinformatics analysis pipeline make direct comparisons extremely difficult (Hirai et al., 2015; Bucklin et al., 2016).

As expected, copepods were the dominant component of the Arthropoda inhabiting the upper 100 m of the epipelagic zone. This is consistent with previous studies across the northern and central Red Sea that have reported similar findings in terms of copepod abundance (Khalil and El-Rakman, 1997; El-Sherbiny et al., 2007; Dorgham et al., 2012; Pearman et al., 2014) as well as with studies in other marine basins (e.g., Caribbean: Moore and Sander, 1976; Mediterranean: Pancucci-Papadopoulou et al., 1992; the Western Channel: Eloire et al., 2010). In our study, the majority of the reads corresponded to calanoid copepods, while poecilostomatoids and harpacticoids occurred only rarely (Table 2). Remarkably, harpacticoids, specifically the genus Oncaea, have been commonly reported in the northern Red Sea as important contributors to the total diversity of adult copepods across the year (Aamer et al., 2006; El-Sherbiny et al., 2007; Dorgham et al., 2012), while in our study their contribution appears to be low. The exception is at deeper depths during the hot season with a higher presence offshore than inshore. These differences might be due to regional zooplankton patterns but could also be due to a high abundance of adult stages relative to the biomass of juveniles. This would result in a higher proportion of morphologically identified samples compared with molecular techniques (Lindeque et al., 2013). Five calanoid genera dominated the copepod community: Acartia, Canthocalanus, Undinula, Candacia, Paracalanus, showing differing abundances throughout the year. Similar records were reported by morphological and molecular studies in the Red Sea (Aamer et al., 2006; Cornils et al., 2007; Pearman et al., 2014).

In summary, this study has shown that the seasonal fluctuation of the zooplankton communities living in the upper water column (top 100 m) in the central Rea Sea is driven mainly by temperature and the annual change between vertical mixing and stratification. During the cold season, communities are characterized by a lower biodiversity in terms of number of species and a higher biomass than communities found during the hot season. All throughout the year the zooplankton communities are dominated by class Maxillopoda, mainly represented by copepods and class Hydrozoa. The rise in the water temperature favors the appearance of classes Malacostraca, Cephalopoda, Gastropoda, and Saggitoidea, among the most abundant.

We found evidence for fine-scale geographic differences in terms of biodiversity for the two close proximity areas under study, providing an essential foundation for future comparative studies.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.