Marine Microbial Community Composition During the Upwelling Season in the Southern Benguela

The microbial communities of the southern Benguela upwelling region were sampled quarterly through 1 year, with sampling for prokaryotes taking place in May, September, November and February, spanning the 2015–2016 upwelling season. Picoeukaryote samples were taken in November and February only. Community dynamics were assessed at stations both inside and outside a typical upwelling site. 16S and 18S rRNA amplicon results, respectively, revealed differences in both bacterioplankton and picoeukaryote communities in both space and time (season). There was a significant difference between sites in picoeukaryote community structure and diversity during the upwelling season in February, but not in November. Prokaryote community structure showed significant changes by water type as well as by sampling time or site. The parasitic dinoflagellate, Syndiniales, dominated February samples, and diatoms (Mediophyceae) mostly occurred in November samples, with nitrate driving community structure. Prokaryote results revealed presence of Nitrosopumillus, an ammonium oxidizer, offshore in February. Nitrospina sp., a nitrite oxidizer, was also present in September in hypoxic and deep water samples. This study reveals significant changes in community variability, leading to shifts within interspecies interactions in this region in response to upwelling events. This has far reaching implications with regard to biogeochemical cycling and ecosystem functioning at the microbial level.


INTRODUCTION
The Benguela upwelling is one of four major global eastern boundary upwelling systems (Hutchings et al., 2009). The Southern Benguela is situated on the west coast of South Africa, separated from the Northern Benguela by the powerful Lüderitz upwelling cell at 26 • S. The Southern Benguela is characterized by pulsed, seasonal, south-easterly wind-driven coastal upwelling at various coastal points resulting in high primary productivity in coastal waters here. As a result of this high primary productivity, the shelf regions are subject to naturally occurring low-oxygen events inshore (Pitcher et al., 2014). These events are well known for their negative effects on resident fish, crustaceans and benthic invertebrate communities, consequently impacting trophic energy flow (Blamey et al., 2013). Bacteria, however, remain active in low-oxygen water, performing essential functions in the remineralization of organic matter and biogeochemical cycling (Letscher et al., 2015).
Existing research on prokaryotes in upwelling regions has revealed significant links between bacterial production and primary production, emphasizing the important role prokaryotes play in the southern Benguela (Painting et al., 1993), Northwest Indian Ocean (Wiebinga et al., 1997), off Concepcíon in central Chile (Cuevas et al., 2004) and in the Northwest Iberian Peninsula (Prieto et al., 2016). Prieto et al. (2016) in particular revealed how, in the winter months, primary production only responded to nutrient amended waters when heterotrophic bacteria were active. This implies that phytoplankton may rely on some sort of secondary metabolite from bacteria in order to become active.
The pico-sized microbial community (defined in this paper as all prokaryotes and eukaryotes between 0.2 and 3 µm in diameter) can contribute significantly to autotrophic carbon biomass and primary production in the ocean, both from a global perspective (Buitenhuis et al., 2012;Worden et al., 2015) and a more localized one such as the Pacific Ocean (Worden et al., 2004) and the suptropical and tropical NE Atlantic Ocean (Jardillier et al., 2010). High nutrients associated with newly upwelled shelf water typically stimulates high primary production, dominated by diatoms, followed by high bacterial abundance (Montecino et al., 2004;Lamont et al., 2014). An overall reduced biomass tends to be present beyond the shelf edge as upwelled waters age and move offshore. Primary productivity and bacterial abundance decrease in aged upwelled water, due to light and nutrient limitations, and grazing pressure.
Few studies have focused on the microbial composition of plankton communities of the world's upwelling regions. These include Painting (1989) who also focused on the Southern Benguela, Kerkhof et al. (1999) investigated the North Atlantic just off the New Jersey coast, Suzuki et al. (2001) investigated an upwelling region in Monterey Bay in the NE Pacific and Teira et al. (2009) studied six key bacterial groups in the Northwest Iberian Peninsula coastal upwelling waters. Only a handful of studies have included next generation sequencing techniques, which are capable of describing community diversity and structure in great detail (e.g., Bergen et al., 2015;Aldunate et al., 2018). These studies revealed a high occurrence of Bacteroidetes, Roseobacter, and SAR86 clade of gammaproteobacteria in upwelling waters. Taxonomic data taken from SSU rRNA (smallsubunit ribosomal RNA) gene sequences have also revealed patterns in microbial community composition within oxygen deficient environments. Common abundant phyla in low oxygen waters consist of Proteobacteria, Bacteroidetes, Marine group A (previously SAR406 clade), Actinobacteria, and Planctomycetes (Wright et al., 2012). Other taxa accompany these groups, depending on the physical and chemical properties of the site in question. Rather than existing on their own, there is increasing evidence to suggest that groups of microbiota work together within specific niches (Koeppel et al., 2008). More specifically, endemic patterns emerge at the OTU (operational taxonomic unit) level, reinforcing the ecotype model, where genetically interrelated populations demonstrate niche-specific ecological or biogeochemical roles (Koeppel et al., 2008).
The microbial food web in the Southern Benguela is exceptionally variable because of spatial and seasonal variability in the physical environment, which controls nutrient and light availability to surface waters (Boyd and Nelson, 1998;Hutchings et al., 2009). Typically, this coastal upwelling region experiences seasonal upwelling as a consequence of south-easterly winds in spring and summer, with resulting Ekman transport. Plumes of upwelled water move offshore and northwards, eventually mixing with offshore, and aged upwelled waters (Andrews and Hutchings, 1980). Inshore in St Helena Bay, the north flowing Benguela jet current, along with wind stress shear, forms a clockwise circulation. This encourages water retention, and therefore an accumulation of biological material in the bay (Lett et al., 2006).
Previous studies in the Southern Benguela region have measured bacterial biomass during upwelling events. Bacteria and picoeukaryotes tend to peak post upwelling as a result of high dissolved nitrogen and carbon concentrations caused by decaying phytoplankton blooms (Painting, 1989;Moloney et al., 1991;Painting et al., 1993). In St Helena Bay, dinoflagellate blooms typically develop during late summer to early autumn. These blooms accumulate inshore because of down-welling, resulting in episodic hypoxic to anoxic conditions (Pitcher et al., 2014). Such events can have a detrimental effect on commercially important rock lobsters Jasus lalandii, which have been forced out the water and stranded inshore due to low oxygen conditions resulting from these blooms (Cockcroft et al., 2008). An in-depth understanding of the spatiotemporal variability of bacterial and picoplankton composition and distribution would help understanding of biogeochemical processes in this region.
In this study, DNA samples were taken in austral autumn (May), early and late spring (September and November 2015) and summer (February 2016) from surface waters and waters at or near the fluorescence maximum Fmax in coastal, shelf and offshore stations from St Helena Bay. The study aimed to describe the microbial community composition through next generation sequencing of the 16S and 18S rRNA genes, and to assess the degree of variability in community composition in relation to changes in hydrography associated with different upwelling stages. We predict that the microbial communities would vary between the various degrees of upwelling and between inshore vs. offshore waters. This study provides insight into the seasonal and spatial (vertical and horizontal) variability in microbial plankton communities in the St Helena Bay region and reveals how community structure responds to its environment.

Sample Collection
Seasonal samples were collected during four cruises of the Integrated Ecosystem Program (IEP) of the South African Department of Environmental Affairs (Verheye et al., 2018 (Figure 1). Stations 3, 5, and 13 were located over the continental shelf (<200 m depth) and stations 7 and 8 were over the continental slope (>200 m depth). Station 3 is assumed to represent inshore waters, stations 5 and 13 shelf waters, and stations 7 and 8 offshore waters.
DNA samples for 16S and 18S sequencing analysis were collected through size-fractionated filtration. A volume of 1 L of water was collected from the Niskin bottles, pre-filtered through a 200 µm mesh, then filtered through a 2 µm filter followed by a 0.2 µm polycarbonate filter. Both filters were immediately flash frozen in liquid nitrogen and stored at -80 • C until processing.
Samples for flow cytometry analysis were collected in triplicate from the Niskin bottles in 2 mL cryovials, preserved with 0.25% glutaraldehyde, flash frozen in liquid nitrogen and stored at -80 • C until analyzed.

DNA Extraction and Amplification
For logistical and funding reasons, two different methods were used for 16S DNA extraction and amplification. We collaborated in 2015 with the Liu lab, which sequenced our May-September 16S samples using 454 pyrosequencing. For our November and February samples, funding became available to send our samples for Illumina sequencing (16S and 18S). As a result of these logistical issues, it must be noted that no 18S data were sequenced for the May and September 2015 samples.

Method 1: May and September 2015 Samples
DNA was extracted using the DNA blood and cell culture Micro Kit (QIAGEN R ) according to their protocol. Extracted DNA was amplified using the general prokaryotic SSU rRNA gene V3-V4 hypervariable region primers, U341F (5 -CC TAC GGG RSG CAG CAG -3 ) (Watanabe et al., 2001) and 787R (5 -CTA CNR GGG TAT CTA A-3 ) (Anna et al., 2010). In brief, the polymerase chain reaction (PCR) was carried out with a 25 µL Master Mix containing 2.5 µL of 10x buffer (100 mM Tris-HCl, pH 8.3 at 25 • C; 500 mM KCl; 15 mM MgCl 2 ; and 0.01% gelatin), 1 µL of MgCl 2 (25 mM) and 0.5 dNTPs (10 mM), 1 µL of forward and reverse primers (10 µM), 0.25 µL of Platinum Taq polymerase (Invitrogen), and 1 µL of DNA template. The PCR program consisted of an initial incubation step at 94 • C for 4 min, followed by 29 cycles with a denaturation step at 94 • C for 30 s, an annealing step at 55 • C for 30 s, and an extension step at 72 • C for 30 s. These cycles were followed by a final extension step at 72 • C for 7 min. Amplifications were verified on a 0.8% agar gel stained with ethidium bromide. MID adaptor linked primers were used for sequencing.

Method 2: November 2015 and February 2016 Samples
DNA was extracted using the DNA blood and cell culture Micro Kit (QIAGEN R ) according to their protocol. Gene amplification of 16S rRNA was carried out using the 16S V3-V4 primer set 341F (5 -CCTAYGGGRBGCASCAG-3 ) and 787R (5 -CTA CNR GGG TAT CTA A-3 ) (Anna et al., 2010). The V4 region of the 18S rRNA gene was amplified using TAREUK454FWD1 (5 -CCAGCASCYGCGGTAATTCC-3 ) and TAREUKREV3 (5 -ACTTTCGTTCTTGATYRA-3 ) (Stoeck et al., 2010) with the barcode. All PCRs were carried out with Phusion R High-Fidelity PCR Master Mix (New England Biolabs). Amplifications were verified on a 2% agarose gel stained with ethidium bromide. MID adaptor linked primers were used for sequencing. PCR products were mixed in equidensity ratios. This PCR product mixture was purified with Qiagen Gel Extraction Kit (Qiagen, Germany). The libraries were generated with TruSeq R DNA PCR-Free Sample Preparation Kit and quantified via Qubit and Q-PCR, to be analyzed by the HiSeq Illumina platform.

May and September 2015 Samples: 454 Sequencing and Bioinformatic Analysis
All nine samples were sequenced in replicate on a 454 GS Junior system (Roche) according to Daigle et al. (2011). All the resulting sequences were screened by removing the low-quality reads, such as reads shorter than 300 nucleotides, incomplete or inaccurate primer sequences, chimeras, and reads with unidentified nucleotides present. Shhh.seqs (a MOTHUR implementation of Amplicon noise) (Quince et al., 2011) was performed on all sequences prior to analysis. In order to minimize the computational burden of all downstream analyses, all unique sequences were identified using MOTHUR (Schloss et al., 2009) and a taxon was assigned to each sequence. A 97% cutoff was applied to all 16S sequences for subsequent analysis using MOTHUR, and the number of sequences per raw sequence was indicated in the OTU table (Table 2). To test for significant differences in bacterial composition between sites and water type, we used permutational multivariate analysis of variance (PERMANOVA) with the function "adonis" in the "vegan" R package. The microbial community diversity ( Table 1) was estimated using the summary single command in MOTHUR. Graphs of species richness, inverse Simpson and constrained ordination plots were generated in R using the packages RAM (McMurdie and Holmes, 2013; Wen Chen and Levesque, 2018), vegan (Jari Oksanen et al., 2017), and ggplot2 (Wickham, 2009). The R scripts are available upon request. Differences between platforms was tested in R using the ADONIS function in vegan.

November 2015 and February 2016 Samples: Illumina Sequencing and Bioinformatic Analysis
Paired-end reads were assigned to samples based on their unique barcodes and truncated by cutting off the barcode and primer sequence. Paired-end reads were merged using FLASH (V1.2.7) 1 (Magoč and Salzberg, 2011). Quality filtering on the raw tags was performed under specific filtering conditions to obtain the high-quality clean tags (Nicholas et al., 2012) according to the Qiime (V1.7.0) 2 (Caporaso et al., 2010) quality control process. Subsequent sequence analysis was performed using MOTHUR version 3.9 (Schloss et al., 2009). All the resulting sequences were further screened by removing the low-quality reads, as above. A 97% cutoff was applied to all sequences for subsequent analysis using MOTHUR, and the number of raw sequences as well as OTU number is indicated in the OTU table ( Table 2). Sequences are uploaded on European Nucleotide Archive under accession numbers ERS3233437-ERS3233461.
Once OTU tables were generated, taxonomic affiliations were assigned using SILVA (16S) and PR2 (18S) databases. Heat trees were generated using the R package Metacoder (Foster et al., 2017) and diversity index box plots as well as Figure 4 heatmaps were generated using the Heatplus package (Ploner, 2019). Heatmaps were generated by extracting the top 40 orders from log transformed OTU abundance data. Scripts are available upon request. Bubbleplots (Supplementary Figure S3) were generated using relative abundances.

Network Analysis
Spearman rank order correlations were calculated between the log transformed top 37 bacterioplankton and picoeukaryote families from the November and February cruises, as no 18S data are available for the May and September cruises. Both the correlation matrix (R matrix) and significance matrix (P matrix) were calculated using the "Hmisc" package in R (Harrell and Dupont, 2019) by calculating all possible pairwise Spearmans Rank correlations between all families. Only robust (Spearman's correlation coefficient >0.8) and statistically significant (p < 0.01) correlations were further considered. The resulting values were plotted in Cytoscape v3.7.2. Other data including the relative OTU distribution of each node (Family) were also imported into Cytoscape, and plotted though the enhancedGraphics package (John et al., 2014).

Flow Cytometry
Heterotrophic bacteria were enumerated using a flow cytometer (LSR II, Becton Dickinson, Heidelberg, Germany) using the method of Marie et al. (1997)  Diversity is based on an operational taxonomic unit cutoff of 97% sequence identity. The Shannon index, Simpson's index of evenness and sequencing coverage are shown. abundance was distinguished from picoeukaryotes through their higher PE signal.

RESULTS
Samples collected from all four cruises were typically comprised of surface samples and samples from the fluorescence maximum, which occurred within the top 10-20 m of the water column ( Table 1). Station 13, which is characterized by perennial low oxygen concentrations in the midwater zone, was sampled in September at the surface, in low oxygen waters (2.42 mL O 2 .L −1 ) at 59 m and near the bottom (105 m). For comparison during this cruise, station 5 was also sampled at the surface and near the bottom (153 m) ( Table 1).

Hydrographic Conditions
Of the surface and near-surface samples obtained during the four cruises, all the shelf samples were from mature upwelled (type 2) water except for the summer (February) surface sample from station 5, which had characteristics of aged upwelled (type 3) water ( Table 1). The surface and near-surface offshore samples had characteristics of mature and aged waters, and samples deeper than 20 m had low oxygen concentrations (<4 mg O 2 .L −1 ). The temperature profiles (Figure 2) showed cooling with depth, with most of the water columns being weakly stratified. The surface waters of the stations closest to the coast (stations 3, 5, and 13) generally were cooler than those at offshore stations (7 and 8). Deep samples from station 5 showed evidence of low-oxygen waters ( Table 1). Salinity did not vary, ranging from 34.55-35.25 PSU throughout the dataset.

Nutrients
All surface nutrient concentrations were low to moderate, with water characterized as either mature-upwelled (May, September) or aged upwelled (November, February) ( Table 1), indicating there had not been recent upwelling. The measured nutrient values (Table 1) can be compared with concentrations found in newly upwelled water: 20.8 µmol.L −1 for nitrate, 16.6 µmol.L −1 for phosphate, 1.88 µmol.L −1 for silicate (Brown and Hutchings, 1987). Surface nutrient concentrations in May and November were generally highest near the coast and decreased offshore. In September and February, the inshore stations were nutrientdepleted at the surface compared to the offshore stations for nitrate (Figures 2b,d), phosphate (Figures 2j,l) and silicate (Figures 2n,p), but tended to have elevated surface concentrations of nitrite (Figures 2f,h). Silicate concentrations were generally lower at the surface than at depth for all stations and months (Figures 2m-p), indicating probable diatom presence in surface waters. Nutrient concentrations at the F max (in May, November and February) were generally high. In September the lowest nutrient concentrations were at the surface, with elevated concentrations in the low oxygen water and near the bottom (Figure 2, second panel).
Surface waters were generally well oxygenated (Table 1), but oxygen concentrations decreased to <4 mg.mL −1 in September in the mid-water low oxygen area and near the bottom at stations 5 and 13 in September, which were the deepest stations sampled.

Bacterioplankton Diversity
Bacterioplankton species diversity was greatest in September shelf surface waters ( Table 2). May offshore waters had lower OTU numbers and species diversity. September samples had similar diversities at stations 13 (Shannon index: 6.15-6.45) and 5 (Shannon index: 5.68-6.31). Taxon diversity in hypoxic waters in September (OMZ and Deep, Table 2) did not differ much from that in the shelf surface waters. In November and February community richness was highest at the surface at the offshore station ( Table 2) and lowest at the surface at the shelf station. Diversity did not vary much among seasons, but there were significant differences in taxon diversity between shelf and offshore samples (Figure 3).
The different platforms used must be considered when comparing May and September 16S data with the February and November results (Loman et al., 2012;Smith and Peay, 2014;Sinclair et al., 2015). All three of these studies have revealed that differences between the two platforms are insignificant, given that similar PCR and bioinformatic protocols were followed. In this study, despite the fact that the primary bias lies within the 454 results which yielded a lower throughput compared to the Illumina data, no significant difference was found between the 454 and Illumina datasets (AMOVA, p > 0.5). As a result, these results have been reported in a single heatmap, however, separate ordination plots were shown in order to properly visualize the variability within hypoxic samples vs. non-hypoxic samples from the rest of the dataset.

Picoeukaryote Diversity
The picoeukaryote community richness was highest in February at the surface at station 5 (Shannon index: 4.07) and lowest at the F max at station 8 in February (Shannon index: 2.38) ( Table 3). Diversity in November was highest at station 5 F max (Shannon index: 3.81), and lowest at station 5, surface waters (Shannon index: 2.85).

Bacterioplankton Community Composition
For May and September samples, 454 sequencing yielded a total of 22,956 sequence reads after denoising, removal of reads present only once and chimeric sequences. The sequence dataset was clustered into 11,189 OTUs at a 97% cutoff.
Abundance data of the top 40 orders (Figure 4), represent almost 90% of the dataset. Species community composition in all samples was dominated by Bacteroidetes (Flavobacteriaceae) and Proteobacteria (Rhodobacteraceae, Oceanospirillales, Rhodobacterales) ( Figure 4A). May surface waters (Station 3) contained the most OTUs but the lowest abundance of the SAR clades, Desulfobacterales and Rhizobiales. May and September samples clustered by water type. September hypoxic waters contained a large percentage of ZD0405 clade of SAR11 (8 and 10% of station 13 OMZ and bottom, 21% of Station 5 bottom). Other notable groups included a marine group B clade, Nitrospina and the Deep I SAR11 clade. In September surface samples, the "core" Rhodobacteriaceae and Flavobacteriaceae dominated the samples overall (20-23%). Further inshore in May, station 3 exhibited higher abundances of Verrucomicrobiaceae, Oceanospirillaceae, Methylophilaceae, and Thiotrichales. In September, shelf surface waters were occupied by few OTUs but included higher numbers of Rickettsiales and SAR406 compared to all other samples.
Two percent of the May and November samples contained the nitrite oxidizing bacterium, Nitrospina, and 1.5% of the November-February samples contained the ammonium oxidizing archaea, Nitrosopumillus (Figure 4).

Ordination
The constrained ordination using log-transformed data ( Figure 5) illustrates how spatial, temporal and environmental variables within the bacterioplankton samples influenced the distribution patterns of the OTUs. For the 16S community structure, 58% of the variability in community composition was explained by the CCA plot ( Figure 5A). The samples grouped by month and by distance offshore for February samples. Nitrate, nitrite, silicate. Phosphate and dissolved oxygen concentrations  Diversity is based on an operational taxonomic unit cutoff of 97% sequence identity. The Shannon index, Simpson and Shannon's indices of evenness and sequencing coverage are shown. all significantly affected the bacterioplankton community variability in November 2015 and February 2016 16S samples, and silicate, phosphate, dissolved oxygen and temperature significantly affected the May and September 16S samples ( Table 4, p< 0.005). 49% of the picoeukaryote community structure was explained by the first two axes (Figure 5C). The samples grouped according to month but not by distance from shore. Nitrate, nitrite, silicate, dissolved oxygen, and Synechococcus abundances were found to significantly affect community variability of the picoeukaryotes (p < 0.01).

Network Analysis
A trend is immediately visible between the two main clusters in Figure 6, with the top cluster having large contributions from families from offshore samples, and the bottom cluster including mostly shelf samples. Noteworthy associations in the top (offshore) cluster include Dinoflagellate group II clades 10 and 11, 8 and 6, which were highly connected with each other, as well as resident picoeukaryotes such as Pelagomonadaceae

DISCUSSION
This study is one of the first of its kind to investigate both prokaryotic and picoeukaryotic community composition and structure in the coastal upwelling region of the Southern Benguela using next generation sequencing methods. Stations 3, 5, and 13 lie in a region that typically experiences active upwelling (Shillington et al., 2006;Hutchings et al., 2009) in the spring and summer months, and stations 7 and 8 are located over the continental slope, receiving modified or aged upwelled water that has moved offshore. Upwelling stages usually start in September and last until the equatorward southeasterly winds stop blowing, around March. Results revealed a shift in community compositions in response to changes in hydrography associated with upwelling stages, with offshore and shelf stations presenting significantly different microbial community compositions in the February bacterioplankton data, and with markedly different microbial communities found in deep, hypoxic (dissolved oxygen <2mg L −1 ) water samples. As wind data only indicate a strong south-easterly wind in the days immediately leading up to the September sampling date (Supplementary Table S1), shifts in community distributions are more likely to be responding to less recent upwelling events shaping the water column (i.e., mature and aged upwelled waters).

Physical Environment
Temperature and salinity profiles suggested that May samples were from a stratified environment with evidence of a dying dinoflagellate (Ceratium sp.) bloom, which took place at station 3. Nitrate profiles provided additional evidence of stratification, with low surface values indicating nutrient consumption by a similar bloom at station 5 (Heisler et al., 2008). September nutrient data indicated different stages of upwelling at stations 5 and 13. Station 13 exhibited 6-fold higher surface nitrate concentrations compared to station 5 (Figure 2) suggesting stratified, bloom-like conditions at station 5, and a recent upwelling event at station 13. This is further supported by Cape Colombine wind data (Supplementary Table S1). These two stations are not far from each other (∼1-2 km) within St Helena Bay, where topography and bathymetry influence water retention. This is further encouraged by the equatorward Benguela coastal jet which, along with wind stress shear, promotes a clockwise circulation here (Penven et al., 2000). This phenomenon plays a role in causing physical, chemical and biological variability at the shelf stations (Nicholson, 2010). Shelf surface water temperatures in both November and February were warmer than 11 • , indicating that, for both of these sampling periods, the water column was relaxed, or in between upwelling events (Brown, 1984) which typically take place on a 7-10 days cycle in the austral summer months (Jury et al., 2010). Nutrients such as nitrate and silicate were depleted compared with those at the F max , indicating recent bloom-like activity in a stratified environment.

Bacterioplankton Community Diversity, Composition
Bacterioplankton species diversity was consistently higher offshore during February and November (Figure 3 and Table 2). It is hypothesized that bacterial production inshore would be expected to be high due to increased availability of carbon substrate in the form of dying dinoflagellate blooms and picoeukaryotes (Supplementary Figure S1). However, only a select group of bacteria are specialized copiotrophs (i.e., Rhodobacteriales and Flavobacteria) (Dang et al., 2007). This would result in a low species diversity inshore as these primary colonizers can outcompete their slower-growing counterparts. As these waters mature and move offshore, the dissolved organic carbon (DOC) brought to the surface will undergo photochemical transformation, rendering it more bioavailable (Mopper et al., 1991). This would allow for an increased number of bacterial species to take hold and thrive. These mature upwelled waters therefore represent a transitional state from upwelling to aged upwelled waters, and this transition would introduce increased resource availability for the prokaryotic community. Low biodiversity offshore in May compared to that in stratified shelf waters could be the product of water current variability in this region. Offshore stations that previously received mature upwelled water during upwelling (species succession) would no longer do so, allowing different water masses to move in (Brown and Hutchings, 1987). Deep and hypoxic waters in September did not show a huge difference in biodiversity compared to surface waters. The high nutrients and available dissolved organic carbon (DOC) here favor the development of micro niche environments, allowing growth of different bacterial taxa despite lower available oxygen. Rhodobacteriales and Flavobacteria were dominant everywhere, but especially in shelf waters (Figures 3B, 4A). Rhodobacteriales, in particular, is a known keystone primary colonizer, one of the fastest users of accumulated nutrients and playing important roles in carbon and sulfur cycling and trace metal reduction (Dang et al., 2007). Subsequent colonizers are Oceanospirillales, Bacteroidetes (using primarily detritus as a carbon source; Bauer et al., 2006), Sphingobacteriales and Alteromonadales. Bacteroidetes were not very abundant in this dataset, and Alteromonadales was most abundant in May and September samples, most likely thriving on large amounts of substrate material, especially diatoms, for which they are major colonizers (Bidle and Azam, 2001). SAR92 (Alteromonadales), SAR11 and SAR86 (Oceanospirillales) are ubiquitous, both in this dataset and elsewhere. Their primary mode of metabolism is the heterotrophic assimilation of dissolved organic carbon, of which there is an abundance in this highly productive region (Giovannoni et al., 2005b;Grote et al., 2012;Tripp, 2013;Cameron et al., 2014). Proteorhodopsin, a light dependant proton pump, has been discovered in these same orders, suggesting they are also capable of phototrophy (Sabehi et al., 2004;Giovannoni et al., 2005a;Stingl et al., 2007). This explains their abundance in both nutrient-and carbon-rich surface waters, and in lightdeficient waters at depth. Bacteria containing proteorhodopsin have a fitness advantage compared to their competitors, as they are capable of surviving prolonged starvation in high light conditions (Gomez-Consarnau et al., 2010). SAR86, in particular, which has the potential for proteorhodopsin-based ATP generation (Dupont et al., 2012), was abundant in mature upwelled waters and at depth (Figure 4).
Some of the top 16S OTUs in hypoxic and bottom waters included prokaryotes such as Desulfobacteriales, and Desulfuromonadales ( Figure 4A). All of these bacterial orders are specialist sulfate and iron reducer groups at depth (Miletto et al., 2011), reducing sulfur to sulfide, partly freeing up particulate organic matter at depth. Despite September shelf Only the Spearman rank correlations that were statistically significant (P < 0.01) and strong (r > 0.8) are shown. Edges with dashed lines represent negative associations (r-values), and nodes are sized relative to how many edges are connected to each node.
samples showing drastically different physical and chemical signatures (station 13 upwelling vs. station 5 relaxed water column), their bacterioplankton community structure clustered together, alongside May offshore surface samples. This could be due to similar carbon sources in each of these samples. Further studies are needed to confirm this. May offshore bacterioplankton clustering with these samples is most likely explained by the high surface water current variability in this region (Lutjeharms, 2007). Surface waters are displaced offshore quickly, often circulating back into the bay.
Community composition varied for the less abundant 16S reads, particularly between non-upwelling and upwelling months. Orders such as Chromatiales and the nitrite-oxidizing Nitrospina were present in May and September, but not in November and February. Chromatiales were mostly present in hypoxic or bottom waters and, in May, in surface waters at station 3. The dominant family present in these samples was Granulosicoccacae, a known primary colonizer of kelp, which is also abundant in these coastal waters (Bengtsson et al., 2011(Bengtsson et al., , 2012. A large number of prokaryotic orders also were present in November and February, such as Verrumicrobiales, Bdellovibrionales, KI89A-clade (Gammaproteobacteria) and Clostridiales. Clostridiales in particular are known obligate anaerobes (Zhao et al., 2003), and could have been be kept alive in oxygenated waters through inhabiting upwelled aggregated organic matter (Bristow, 2018).

Nitrite and Ammonium Oxidation
Both nitrite oxidizing bacteria and ammonium oxidizing archaea were present in this dataset ( Figure 4A). These results are most likely skewed due to the fact that the primer set used for the 16S samples did not have sufficient coverage of Archaea, and so caution must be paid to these results. Despite this, Nitrospina have been known to increase activity in late spring at in situ light levels (Levipan et al., 2016). A significant amount of nitrification is known to take place in coastal and upwelling waters (Fernandez and Farías, 2012;Levipan et al., 2014), performed by ammonium oxidizing archaea and bacteria. Nitrosopumillus maritimus, the ammonium oxidizing archaea found in the November and February samples, were previously thought to be inhibited by light. However, new strains (PS0 and HCA) are proving to have a high tolerance to light (Qin et al., 2014). Cheung et al. (2019) recently redefined these previously defined AOA clades, revealing diverse sub lineages containing varying ecophysiological characteristics. Further investigation is needed in this region in order to properly understand how and who is responsible for regenerating nitrogen.
Bacterioplankton clustering was influenced by nitrogen species as well as silicate and dissolved oxygen concentration (MANOVA, p < 0.001) ( Figure 6A). Given the presence of nitrifying bacteria and archaea in September shelf and February offshore samples, respectively, these organisms probably play an important biogeochemical cycling role in this ecosystem. Correlations with silicate suggest that these species are using senescent diatom populations as their carbon source, i.e., Amin et al. (2012). In general, the community structure agrees with previous studies from upwelling regions. Alonso-Gutiérrez et al. (2009) revealed seasonal changes in the bacterioplankton community structure in the coastal upwelling region of "Ría de Vigo, " NW Spain, driven by Roseobacter, Bacteroidetes, and SAR86, as in this study. Bergen et al. (2015) also revealed high abundances of Bacteroidetes and Gammaproteobacteria in productive upwelled waters in the Northern Benguela, in contrast to an abundance of Pelagibacterales in less productive waters, as also found in this study (i. e., Figures 4, 6).

Picoeukaryote Diversity and Composition
Picoeukaryote species diversity was highest in February in surface waters, most likely due to high nutrient and light availability. High biomass in surface waters was associated with a decreased biodiversity at F max at this same station, likely due to light limitation from surface waters. A drop in species diversity at station 5 surface waters in November compared to the F max sample can most likely be attributed to fewer available nutrients and, therefore, lower biomass in surface waters, allowing for higher light availability in F max waters.
The significant split in picoeukaryote community structure from November to February confirms a seasonal shift in the picoeukaryote community as the upwelling season progresses. Likewise, the split in community structure from shelf to offshore in February vs. the lack of one in November (especially for picoeukaryotes) illustrates upwelling-induced spatial variability, as shown in other upwelling systems in the California current and off northwest Africa (Baltar et al., 2007;Alonso-Gutiérrez et al., 2009;Teira et al., 2009;Allen et al., 2012).
In February, the parasitic Syndiniales dominated our dataset (Figures 3A, 4A). With the increasing use of deep sequencing technologies to describe marine microbial communities, these parasites, which are very difficult to identify under the microscope, have been shown to have a large presence in many microbial datasets. Known to be strictly obligate parasites , the presence of Syndiniales in such high biomass sampling sites implies their active involvement in mitigating the spread of algal blooms. Syndiniales split primarily into two groups, Dino gp I and II, each representing 12 and 26% of the whole dataset, respectively. This agrees with the current literature Chambouvet et al., 2011). Further studies are underway investigating their hosts at this place and time. Figure 6 suggests they may be associated with diatoms. No known study to date has documented Syndiniales in such high numbers in this region.
Picoeukaryotes distributed significantly by their trophic associations (Supplementary Figure S3). November waters were occupied primarily by diatom species, mostly Polar centric Mediophyceae, defined as "floating autotrophs" in this figure.
Despite the fact that these samples are coming from aged upwelled waters with reduced nutrient concentrations, it is likely that this species' signature is resulting from a dying bloom that has stripped most nutrients from surface waters. February samples, however, classified as more mature or aged upwelled waters were dominated by parasites as described above. Different species assemblages occupied these waters inshore vs. offshore. Not enough is known about these to further discuss their ecological function here. Interestingly, a higher percentage of swimming autotrophs as well as mixotrophs were found in February waters compared to November waters, suggesting a life strategy that would benefit in bloom like waters where light is quickly a limiting factor.
The additional cercozoa and a number of diatom species found inshore in November are typical indicators of post upwelling conditions (e.g., Pitcher, 1990). Syndiniales again dominated the F max samples, most likely taking advantage of the sinking dinoflagellate bloom Montagnes et al., 2008;Chambouvet et al., 2011).

Network Analysis: Bacterioplankton -Picoeukaryote Interactions
The network analysis revealed some interesting trends with regard to bacterioplankton-picoeukaryote organism distribution within the samples as well as Dino-Gp II co-occurrences. It has already been established that bacterioplankton species diversity is significantly different between shelf and inshore waters, and eukaryotic diversity is significantly different between shelf and inshore in February but not in November (Figure 3). Dino gp II are known to be highly genetically diverse, occupying up to 30% of some picoeukaryote datasets and infecting a diverse range of hosts . They also previously have been reported to have a high number of correlations with other OTUs in Saanich inlet, a seasonally anoxic fjord (Torres-Beltran et al., 2018). Negative associations were less frequent but mostly occurred between diatoms and Dino gp II clades 1, 3, and 6. The bottom cluster in Figure 6 primarily contained families that were mostly found inshore and rely on fresh carbon typically produced from highly productive inshore/shelf waters. These families (Flavobacteriales, Bacteroidetes, Rhodobacteraceae) can be typically found together in productive surface waters (Buchan et al., 2014) showing evidence of a type of mutualistic metabolism, whereby one species facilitates the growth of others resulting in multispecies biochemical dependencies (Prokopenko et al., 2013).

CONCLUSION
This study provides a snapshot comparison of the microbial (bacterioplankton and picoeukaryote) community in a key upwelling region in the productive Southern Benguela upwelling ecosystem. Picoeukaryotes showed clear spatial variability in communities in November and February, indicating a significant adaptation of the picoeukaryote community as the seasons progressed. 16S OTUs clustered by nitrogen concentration, revealing key prokaryotic signatures as upwelling waters aged and moved offshore. Key carbon metabolisers such as Flavobacteriales and Rhodobacteriales give way to slower-growing species such as SAR11 and Clostridiales. Notable results that deserve further investigation include the presence of ammonium oxidizing archaea offshore, which correlate with the nitrite maximum. If proven to be active, this would indicate that regenerated production is playing a significant role in these waters. The overwhelming presence of Syndiniales parasites in offshore and F max samples in February indicate active carbon and nutrient cycling at a different trophic level. Their further significant negative association with certain diatom species seems to be a novel finding not yet investigated in detail. Further investigations are ongoing to confirm these associations. This could have a large impact on carbon export, because the major processes involved in bloom degradation currently are thought to be zooplankton grazing and bacterial degradation.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the ENA accession numbers ERS3223437-ERS3223461.

AUTHOR CONTRIBUTIONS
ER, CM, and HL designed the study. ER and ZG collected the samples. SC and HL sequenced the May and September 16S samples. ER and ND analyzed the data. All authors read and approved the final manuscript.

FUNDING
This work was based upon research supported by the National Research Foundation project # 98967, and the Swiss Polar Institute Antarctic Circumnavigation Expedition Grant, Project 12 (S. Fawcett). The Applied Center for Climate and Earth Systems Science also partly funded this project.