Spatiotemporal Distribution of Key Pelagic Microbes in a Seasonal Oxygen-Deficient Coastal Upwelling System of the Eastern South Pacific Ocean

The strong seasonal variability in physical-chemical conditions of the Eastern South Pacific Ocean creates an ideal setting to study spatiotemporal distribution of key marine microbial communities. We herein report a nearly 4-year-long time series of the variability in amoA gene counts of ammonia oxidizing archaea (AOA) and bacteria (betaproteobacteria, bAOB) by quantitative PCR, GI.1a Thaumarchaeota and MG-II Euryarchaeota by CARD-FISH, and the picoplanktonic community by flow cytometry for this area. During spring-summer, non-photosynthetic picoplankton such as MG-II Euryarchaeota and GI.1a Thaumarchaeota peaked at the surface and deeper waters, respectively. General AOA and bAOB achieved higher abundances at the oxycline mainly in summer (up to 105–104 amoA copies mL–1). Generalized additive models for location, scale, and shape (GAMLSS) indicated that season and depth account for 19–46% of variations in the abundance of the groups studied, particularly GI.1a Thaumarchaeota and AOA. The oxygen and nitrite concentration were statistically meaningful predictors for the studied groups. GAMLSS models indicate that ammonia oxidizing assemblage’s variability is coupled with ammonia, nitrite, and nitrate variations. Our results indicate that microbial abundances fluctuation is associated with upwelling variability and oxygen-deficient water conditions that shape the substrates availability and metabolic response of marine microbes, including keystone ammonia oxidizing assemblages and their ecological interactions. Overall, our results support planktonic nitrification activity and its contribution to nitrous oxide excess production in the time series off Concepción and the ecological dynamics regarding AOA and bAOB in coastal waters.


INTRODUCTION
Coastal marine areas hold a diverse microbial community but with the predominance of abundant taxa, such as Gammaand Alpha-proteobacteria, Bacteroidetes among other bacterial phyla, and archaea from Thaumarchaea and MGII Euryarchaea which has been reported to have a spatialtemporal dynamics, e.g., surface versus subsurface, and in response to climatic and oceanographic conditions, e.g., San Pedro Ocean Time series station (Cram et al., 2015;Parada and Fuhrman, 2017); Monterey Bay (Reji et al., 2019;Tolar et al., 2020), based on their ecological and metabolic traits (e.g., Northwest Mediterranean coast, Galand et al., 2018;Monterey Bay, Reji et al., 2019). Marine microbes from coastal zones are also subjected to anthropological pressure, such as, nutrient fertilization and other pollutants, influencing biogeochemical conditions and oxygen-deficiency (Cavicchioli et al., 2019).
Besides photosynthetic communities, chemoautotrophic assemblages associated with nitrification have been reported to play a central role in the functioning of coastal microbial assemblages (e.g., network analyses, Parada and Fuhrman, 2017;Reji et al., 2020). Nitrification is a chemoautotrophic process associated with a two-steps aerobic reaction that oxidizes ammonia into nitrite followed by the conversion of nitrite into nitrate. Two functional groups represented by ammonia-and nitrite-oxidizing microorganisms catalyze complete marine nitrification, with the rate-limited first step of this reaction performed by ammonia-oxidizing bacteria (AOB) and archaea (AOA) (Ward, 2008). Nitrospira bacteria able to catalyze one-step ammonia oxidation to nitrate (Comammox) represents the only known exception to this rule (Daims et al., 2015;van Kessel et al., 2015). However, the biogeochemical impact of Nitrospira-like Comammox bacteria on the marine nitrogen cycle is still unclear, since these microorganisms appear to be absent in marine ecosystems (Daims et al., 2016).
Phylogenetically complex natural assemblages of ammoniaoxidizing microorganisms include Thaumarchaea, groups 1.1a, 1.1b, and Hot Water Crenarchaeota Group (HWCGIII) and Nitrosocaldus group (Pester et al., 2011), with the group 1.1a representing one of the most ubiquitous and abundant microbial lineages on earth (Karner et al., 2001). In contrast, AOB are less abundant and often undetectable in marine ecosystems (Mincer et al., 2007;Beman et al., 2010;Tolar et al., 2013). Nonetheless, according to the quantification of the functional marker amoA genes (encoding the alpha subunit of the ammonia monooxygenase enzyme), AOB could be more abundant than AOA. For example, higher ratios of benthic AOB-to-AOA amoA gene copies were detected in freshwater and eutrophic waters (Bouskill et al., 2012) or in association with increased amounts of substrate and salinity in estuarine sediments (Bernhard and Bollmann, 2010;Li et al., 2015). In addition, the abundance and potential activity of AOB (as amoA gene counts and transcript counts, respectively) may be important in other marine environments characterized by strong physical-chemical gradients such as nitrification zones in oxygen-deficient areas (Lam et al., 2007) and near-surface depths (Tolar et al., 2013), including diel variability in sunlit coastal waters (Levipan et al., 2016).
Marine microbial time series (oceanic and coastal) could help to identify long-term trends in the biogeochemical cycles and predict future environmental sceneries, based on their inherent physical-chemical variability. However, there are still relatively few marine time series studies related to the variability of major marine archaeal lineages, besides ammonia oxidizing assemblages, especially in highly productive, oxygen-deficient and dynamic coastal upwelling areas. For example, exhaustive studies have been carried out along the California coast at San Pedro Ocean time series (SPOT) and Monterey Bay Aquarium Research Institute (MBARI) time series (e.g., Murray et al., 1999;Mincer et al., 2007;Beman et al., 2010Beman et al., , 2011Steele et al., 2011;Robidart et al., 2012;Parada and Fuhrman, 2017;Tolar et al., 2020), at the station ALOHA in Hawai (Karner et al., 2001), at Devil's Hole in Bermuda (Parsons et al., 2015), and in marine ecosystems such as Chesapeake Bay (Bouskill et al., 2011), the Mediterranean Sea (Galand et al., 2010(Galand et al., , 2018Lambert et al., 2019), and the North Sea (Wuchter et al., 2006). These studies concluded that AOA are the main responsible for ammonia oxidation in oceanic and coastal areas in general.
The coastal area off central-southern Chile (∼36 • S) in the eastern South Pacific Ocean is an ideal location to determine the effects of oceanographic dynamics on microbial communities and their biogeochemical activities, since this is a highly productive seasonal upwelling system. Upwelling period occurs during the austral spring-summer modulated by southwesterly winds (Sobarzo et al., 2007). The upwelled water is influenced by the Equatorial Subsurface Waters (ESSW), characterized by cold (∼9-12 • C), high-salinity (>34.5 psu), nutrient-rich (reaching > 30 µM of nitrate), and oxygendepleted conditions (< 1 mL O 2 L −1 ) (Ahumada and Chuecas, 1979). In contrast, strong water mixing events caused by intense northerly winds occur during the austral fall-winter, with water column being characterized by warmer temperatures (10-13 • C) and low-nutrient concentrations (e.g., <20 µM of nitrate). Moreover, winter mixing events ventilate the surface and middle depths (>1-6 mL O 2 L −1 ), although low oxygen concentrations (<0.5 mL O 2 L −1 ) persist at bottom depths almost the entire year (Sobarzo et al., 2007). During fallwinter, high precipitation and river-runoff reduce salinity (between ∼32.5 to <34.1 psu) at the surface generating strong stratification in the first 10-20 m depth (Sobarzo et al., 2007;Saldías et al., 2012).
The oceanographic dynamic in central-southern Chile along with the high solar radiation in spring-summer promotes a high primary productivity (Montero et al., 2007;Hernández et al., 2012) and faster remineralization rates (Daneri et al., 2000, that stimulate nitrification activity. Nitrification is a relevant process in the study area, reaching high net rates in the oxycline during the spring-summer (200-316 nmol L −1 d −1 ) compared with winter, c.a. 21 nmol L −1 d −1 (Fernandez and Farías, 2012). Recently, ammonia and nitrite oxidation has been detected in this area at nanomolar concentrations of oxygen (5-33 nM, Bristow et al., 2016). Furthermore, molecular studies indicate the presence of an abundant and rich AOA community (Molina et al., 2010;Bertagnolli and Ulloa, 2017), as well as of Nitrospina spp., as the main drivers of nitrification under oxygen-deficient conditions at bottom waters during the spring-summer time (Levipan et al., 2014). In general, the water column is characterized by higher bacterioplankton abundances in summer (Galán et al., 2012), and Nitrosopumilus maritimus predominance (an important AOA ecotype at the Sta.18) in winter months based on metagenomic data (Murillo et al., 2014). Additionally, annual variability of light intensity can favor transient co-variability between bAOB and certain AOA ecotypes to catalyze the first step of nitrification (Levipan et al., 2016).
Here, we investigate physical-chemical variability contribution in shaping the spatiotemporal distribution of marine microbial assemblages including ammonia oxidizers at a coastal upwelling station (Sta.18) using data from a nearly 4-year-long time series. We hypothesize that vertical and seasonal environmental fluctuations will significantly modulate the abundance variability of the microbial community in the study area, whereas oxygen and ammonium will be explanatory variables shaping the distribution and abundances of keystone ammonia oxidizing assemblages. Temperature, salinity, and dissolved oxygen concentrations were obtained using a CTD-O probe (Seabird 23B Electronics, Bellevue, United States). Water samples for chemical (NH 4 + , NO 3 − , NO 2 − and Chlorophyll-a) and biological analyses (picoplankton, DNA extractions, and CARD-FISH) were collected using a rosette sampler equipped with 10 L Niskin bottles. Triplicated NH 4 + samples (40 mL) were fixed on board and fluorometrically analyzed (Turner Designs 10-AU fluorometer) using the method proposed by Holmes et al. (1999). Seawater samples for NO 2 − and NO 3 − analyses (20 mL) were filtered through glass filter (0.7 µm of pore-size, Millipore) and stored frozen (−20 • C) until further processing (Parsons et al., 1984). Chlorophyll-a was determined in discrete seawater samples according to the method described by Holm-Hansen et al. (1965). In addition, discrete seawater samples were fixed in duplicated with glutaraldehyde (final concentration, 0.1% v/v) and stored at −80 • C for flow cytometry analysis. Non-fluorescent picoplankton was stained with SYBR-Green I (Molecular Probes) and processed following Marie et al. (1997)

DNA Extraction
Seawater samples (3 L) were collected monthly for DNA extraction at different depths (0, 5, 10, 15, 20, 30, 40, 50, and 80 m) from January 2005 to January 2009. The samples were pre-filtered serially through a 20-µm-mesh (Nitex R nylon) and 3-µm-filters (MCE, mixed cellulose ester, Millipore membrane filter), and then filtered onto cellulose ester membrane filters (47 mm diameter, 0.22 µm pore-size GPWP04700) through peristaltic pumping. Each membrane filter was soaked with 1 ml of DNA buffer (50 mM Tris-HCl [pH 9.0], 0.75 M sucrose, 400 mM NaCl), flash-frozen in liquid nitrogen, and stored at −80 • C. The DNA was isolated from thawed filters using the phenol-chloroform extraction method as previously described (Molina et al., 2007). In order to estimate the contribution of different ecotypes, seawater samples (3 L) collected at 0, 5, 10, 20, 30, 50, and 80 m depth from two oceanographically contrasting months (August versus December 2011) were extracted by using the PowerSoil TM DNA Isolation Kit (MoBio Laboratories, Solana Beach, CA, United States) in accordance to the manufacturer's specifications, with minor modifications described by Levipan et al. (2014). All DNA extracts were resuspended in 100 µl of nuclease-free water and their concentration and quality (A 260 /A 280 ratio) were determined by optical absorption spectroscopy on a Synergy Mx Microplate Reader (BioTek Instruments, Winooski, VT, United States). The DNA extracts were stored at −80 • C until further analysis.
Profiling of amoA gene counts for AOA belonging to the shallow and deep water column clades named as WCA and WCB, and N. maritimus (Group 1.1a) were determined in two oceanographically contrasting months (August vs. December 2011). qPCR assays for these ecotypes were performed with the primer Arch-amoAR along with primers Arch-amoAFA and Arch-amoAFB for WCA and WCB amoA, respectively (Beman et al., 2008). The ecotype N. maritimus-like was quantified with the primer Nmar423F (Levipan et al., 2016) used in combination with the primer Arch-amoAR. These primers combination could also amplify other species such as Nitrosopumilus oxyclinae and Ca. Nitrosarchaeum but not the WCA cultivated representative Nitrosopelagicus brevis (Santoro et al., 2015, see alignment in Supplementary Table 2). qPCR standard curves were generated from 10-fold dilutions of amoA gene clones (2 × 10 7 to 20 copies) available in our laboratory for N. maritimus (GenBank accession number KJ555107) and WCA and WCB (Molina et al., 2010). The copy number of these amoA clones was calculated using the mentioned formula. PCR efficiencies and correlation coefficients of standard curves for WCA and WCB (E = 91 ± 6%, r 2 = 0.989 ± 0.012) and N. maritimus (E = 95.492%, r 2 = 0.986) were determined.
All qPCR reactions were performed in triplicate using a StepOne TM Real-Time PCR System (Applied Biosystems, Lincoln Centre Drive, Foster City, CA, United States); data were analyzed with the StepOne TM software package (v.2.2.2; Foster city, CA, United States). Each reaction was conducted on a volume of 20 µL containing 5-10 ng of template DNA, forward and reverse primers (0.4 µM, final concentration), and 2X Fast SYBR R Green Master Mix (1X; Applied Biosystems, 850 Lincoln Centre Drive, Foster City, CA, United States). The qPCR program consisted of an initial denaturation for 20 s at 95 • C followed by 40 amplification cycles consisting of 95 • C for 3 s, 20 s annealing at a temperature that varied depending on the primer pair used (Supplementary Table 1), and an extension of 20 s at 72 • C. A melting curve was performed at the end of each qPCR experiment. The detection limit for standards was always observed at CT values lower than 30. The specificity of qPCR assays was verified by checking dissociation curves and standard electrophoresis of resulting amplicons.

Statistical Analyses
The abundances of total non-photosynthetic picoplankton, Thaumarchaeota (formerly known as marine GI Crenarchaeota; Massana et al., 1997), AOA, and bAOB, were modeled in response to the spatiotemporal variability (season and depth) and chemical variability (ammonium, nitrite, nitrate, and oxygen) using Generalized Additive Models for Location, Scale, and Shape (GAMLSS) from R 'gamlss' package (Rigby and Stasinopoulos, 2005). MG-II Euryarchaea was not included in this analysis because there were temporal gaps in its abundance, but a Spearman correlation analysis indicates that it is significantly correlated with non-photosynthetic picoplankton in surface waters (R = 0.4; P = 0.0013). Predictive variables showed a high collinearity (up to Spearman's R = −0.79, P < 0.001) due to the strong seasonality in wind-driven upwelling conditions. To solve this problem, for each microbial group, we fitted a spatiotemporal model with season and depth as predictive variables, and independent models for every chemical parameter used as a predictive variable (i.e., nitrite, nitrate, ammonium, and oxygen). Before running the analysis, predictive variables were scaled (Z-transformed) to allow a direct comparison between parameters derived from resultant models (Schielzeth, 2010). In order to count both skewness and leptokurtosis of the response variables, microbial groups were modeled by using the follow distributions (Rigby and Stasinopoulos, 2005): (i) negative binomial type I, (ii) negative binomial type II, and (iii) Box-Cox t ( Table 1). Model diagnostic procedures were conducted (to choose the most appropriate distribution) using (i) the plot.gamlss function which provided a series of plots for checking normalized quantile residuals of the fitted models, and (ii) the Generalized Akaike Information Criterion (GAIC). Microbial groups were modeled as a linear function of variables seasonality, ammonium, nitrite, nitrate, and oxygen. The depth was adjusted using a non-parametric cubic smoothing spline function available in R 'gamlss' package.
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Seasonal Regime of Upwelling-Favorable Winds and Physical-Chemical Variables in the Study Area
Upwelling-favorable winds (i.e., southwesterly) usually extended from September (early austral spring) to April (early austral fall) in the study area, however, inter-annual differences were detected during our temporal framework (Supplementary Figure 1) Figure 1).
During the austral spring-summer, upwelled waters were characterized by low temperatures (∼9-12 • C), and high salinities (>34.5) and densities (≥26 kg m −3 ), influencing surface layer; whereas during non-upwelling periods the densest water was just restricted to bottom depths (Supplementary Figure 2). The upwelled water was also characterized by suboxic conditions (≤22 µM), and its influence was perceptible up to almost 15 m depth, mostly from November to April (Figure 1). Oxic conditions (>200 µM) were restricted in the first 5-10 m depth over the entire study period and at intermediate depths during the austral fall-winter (i.e., upwelling-unfavorable period) (Figure 1). In general, during the upwelling season, high concentrations of ammonium (>1 µM) were detected at the base of the mixed layer (∼10 m depth), except during the 2005-2006 upwelling period, when elevated ammonium concentrations (up to ∼3 µM) were found at almost the entire water column (Figure 1). During this period, also high nitrite concentrations (∼7 µM) were normally found in the water layer below 40 m depth (Figure 1), except during the 2005-2006 upwelling period, when higher nitrite concentrations were registered both in bottom waters and oxycline. Surface waters showed lower concentrations of nitrate (≤10 µM) during the entire study period, while the higher concentrations of this nutrient (up to ∼35 µM) were detected at mid-depth waters (between ∼20 and 50 m) but not just associated with the seasonal upwelling signature (Figure 1).

Spatiotemporal Variability of Chlorophyll-a and Microbial Assemblages
During fall-winter months (upwelling-unfavorable seasons), chlorophyll-a concentrations were always ≤ 5 mg m −3 throughout the water column (Figure 2). In contrast, during spring-summer months, the highest chlorophyll-a concentration (i.e., >10 mg m −3 ) was tightly associated with surface waters (between 0 and ∼20 m depth), except by the 2006-2007 spring-summer period, when the lower chlorophylla concentrations were detected and the span of maximal chlorophyll-a concentration finished earlier than the other years (Figure 2). Cell abundances of the non-photosynthetic picoplankton ranged between 0.23-6.50 × 10 6 cells mL −1 in the water column, reaching their maxima at surface and bottom waters during late summer and fall, and minima in winter months (Figure 2).
Cell abundances of GI.1a Thaumarchaea and MG-II Euryarchaea (referred to as Thaumarchaea and Euryarchaea in the following sections) are shown in Supplementary  Figure 3. The abundance of Euryarchaeota range between zero and 78 × 10 4 cells mL −1 , accounting for 9.6 ± 9.0% of total picoplankton community (determined by counting DAPIstained cells) and showing four maxima in the first approximately 30 m depth (Supplementary Figure 3). Moreover, its abundance decreased with depth, with only one Euryarchaeota maxima detected at depth during late summer 2008 (see march-associated peak). In contrast, cell abundances of Thaumarchaea were usually between ∼0.2 and 85 × 10 4 cells mL −1 and more abundant during late summer and winter. This pattern was observed in deeper waters, where this group was accounted for 25.1 ± 13.8% of DAPI counts detected below the oxycline (Supplementary  Figure 3). In addition, the contribution of Thaumarchaea to the whole picoplankton community (as % of the total number of DAPI-stained cells) was as high as 50% below 50 m depth (Supplementary Figure 3). In total, ammonia-oxidizing assemblages amoA gene copies quantification indicated that AOA amoA outnumbered bAOB amoA by up to an order of magnitude (Figure 2). Moreover, if it is considered that one amoA gene copy corresponds to one cell, AOA and bAOB were accounted for 1.2 ± 2.9% and 0.008 ± 0.015% of total abundance of the non-photosynthetic picoplankton (determined from flow cytometry) in the entire water column over time, respectively. In addition, AOA amoA was detected in the entire water column, but especially at the oxycline (up to 10 5 gene copies mL −1 between ∼20 and 40 m), during the summer 2006, 2006, and 2008, high AOA amoA abundance was determined below the oxycline (40 m depth) during fall and winter (between April and August, Figure 2). The AOA amoA -to-Thaumarchaeota cell ratios ranged between 2.75 × 10 −7 and 0.72, showing higher numbers in spring followed by fall, summer and winter (Supplementary Figure 4). The number of bAOB amoA ranged from undetectable (i.e., <2 copies mL −1 ) during the 2005-2006 upwelling-favorable period up to 10 4 copies mL −1 in the oxycline from mid-2006 onward, following the variability of AOA amoA, but usually centered at the oxycline and intermediate depths, when detected (Figure 2).
The comparison of bAOB and AOA amoA gene copies and with the AOA amoA/Thaumarchaeota cell ratios showed that the relative contribution of bAOB amoA was separated in spatial and seasonal scale; bAOB were favored at surface and intermediate waters but mainly during summer (Supplementary Figure 5A). Summertime was characterized by low ratios of archaeal nitrifying assemblages based on AOA amoA/Thaumarchaeota cells ratios (Supplementary Figure 5B). In fact, a plateau of ∼181 bAOB amoA (i.e., ∼2.26 as log-transformed value) was determined when AOA amoA/Thaumarchaeota cells ratios were >0.2, usually associated with surface and intermediate depths during summer (Supplementary Figure 5B).
The spatiotemporal differences in gene number or cell abundance between the studied microbial groups (except for Euryarchaeota) are shown in Figure 3. Data derived from Generalized Additive Models for Location, Scale, and Shape (GAMLSS) statistical analyses are shown in Table 1. In general, a significant decrease in microbial abundance was observed during winter, more evident for ammonia oxidizing assemblages and spring for Thaumarchaeota. GAMLSS analyses indicated that spatiotemporal variations (i.e., season and depth) account for >19.1% of the fluctations, i.e., in non-photosynthetic picoplankton (24.3%), Thaumarchaeota, (45.6%), AOA amoA (34.8%) and bAOB amoA (19.1%) in the water column (see model 1 in Table 1).
The distribution of different archaeal ecotypes (N. maritimus, WCA and WCB) was studied only for selected months during winter and summer of 2011 in order to decipher the relevance of specific groups within the AOA amoA (Figure 4). The ammoniaoxidizing ecotypes presented similar trends with depth than the average values obtain at the time series distribution, showing a weak peak at the subsurface (∼10 m depth) and higher values below 30 m depth for AOA amoA during both summer and winter months. N. maritimus amoA presented higher magnitudes followed by WCA, WCB, and bAOB amoA. The most relevant difference observed was that during winter, N. maritimus, WCA, and bAOB amoA peaked at intermediate depths (50 m depth), whereas during summer, WCB amoA maximum was determined at intermediate to bottom depths (50-90 m) (Figure 4).

Environmental Forcing of Microbial Communities, Including Ammonia-Oxidizers
The distribution of total non-photosynthetic picoplankton (flow cytometry) was significantly and positively associated with oxygen and nitrite, and negatively associated with nitrate and depth (Figure 5). Nitrite account for 15.1% of non-photosynthetic picoplankton spatiotemporal variability ( Table 1). In contrast, Thaumarchaeota cell abundances were negatively correlated with oxygen, ammonium, and nitrite, and positively associated with nitrate and depth (Figure 5). Oxygen (45.9%,ammonium (25.2%) and nitrate (26.8%) were the explanatory variables for the Thaumarchaeota fluctuations ( Table 1). A similar tendency was found for AOA amoA (Figure 5), but the most explicative variable was nitrate (15.6%) ( Table 1). Except for depth, bAOB amoA were significantly correlated with ammonium, nitrite, and nitrate ( Figure 5). Nitrate (11.3%) was the variable with the higher explanatory percentage to the bAOB amoA distribution in the time series (Table 1).

DISCUSSION
The biological activity in the study area is mainly influenced by seasonal upwelling events and its remote-forcing variability; see examples in the special volume for the COPAS time series (Escribano and Morales, 2012;Farías et al., 2015). In the present study, upwelling-favorable conditions were associated with higher phytoplanktonic biomass (chlorophyll) and abundances of total non-photosynthetic picoplanktonic communities, except during the 2005-2006 upwelling-favorable period for the non-photosynthetic picoplanktonic component (Figure 2). In fact, upwelling favorable winds cumulative intensity denoted an evident interannual upwelling variability during our study timeframe, characterized by a delayed and less intense events triggered during 2005-2006 (Supplementary Figure 1). The influence of upwelling variability on nonphotosynthetic picoplankton has been reported in other seasonal upwelling ecosystems, such as, northwestern Indian Ocean (Wiebinga et al., 1997) and Benguela  and agrees with previous results from the study area (Cuevas et al., 2004;Daneri et al., 2012).
In the study area, upwelling not only fertilizes the surface with high-nutrient content waters fueling primary and secondary productivity (e.g., Montero et al., 2007), but also develops a seasonal oxygen deficiency in the subsurface associated with Equatorial Subsurface Waters (ESSW) which was evident during our study (Ahumada and Chuecas, 1979;Sobarzo et al., 2007). Oxygen has been reported as a controlling factor of microbial community abundance (Eissler et al., 2010), structure (Aldunate et al., 2018), and biogeochemical active processes (Galán et al., 2014(Galán et al., , 2017Srain et al., 2020) in the study area. During our study, GAMLSS analyses indicated that the spatiotemporal variability in non-photosynthetic picoplanktonic cell abundance (studied by flow cytometry) was significant in the study area, and among the variables studied, oxygen and nitrite largely account for the changes observed (Table 1). Both, oxygen and nitrite are associated variables in subsurface waters, since nitrite accumulation depends on the redox conditions in the water column triggered by oxygen deficiency (e.g., Farías et al., 2015). Nitrite is a key intermediary of many biogeochemical processes, i.e., as a product of aerobic ammonia oxidation and   nitrate reduction, and as a substrate for aerobic nitrite oxidation, anammox and nitrite reduction related with heterotrophic and chemoautotrophic denitrification (Ward, 2008). In fact, coupling of keystone OMZ microbial communities such as SUP05 clade with anammox and denitrifying sulfur-oxidizers were determined to be associated to nitrite based in a biogeochemical model in the Saanich Inlet seasonally anoxic fjord (Louca et al., 2016) and in bottom waters during upwelling periods in the study area (Galán et al., 2014).

Variability of Euryarchaea Cells Using CARD-FISH
Euryarchaea and Thaumarchaea presented a differential temporal and spatial distributions in the water column. Euryarchaea was characterized by a higher contribution at surface seawater at fall or spring. This result agrees with other findings on a temporal rather than a spatial scale since Euryarchaeota showed a narrower distribution than previously reported (Levipan et al., 2007a;Quiñones et al., 2009). However, sampling methodological constraints could influence these spatial differences, that is, our study was focused on freeliving fraction cell counts, whereas the mentioned reports involved bulk DNA analyses (7 L samples) that favor the presence of different particle sizes. Particle aggregates seem to enhance the contribution of Euryarchaea off northern Chile (Ganesh et al., 2014) and in the study area (Levipan et al., 2007a), where methylotrophic methanogens were abundant and viable at depth in a large-particle fraction (between 0.22 and 25 µm). Moreover, differences in vertical distribution patterns of Euryarchaea caused by the filter size fraction have also been reported in contrasting areas near ALOHA time series station (Lincoln et al., 2014) and Oregon coast . In addition, our results could indicate that using CARD-FISH distinct phylotypes were detected, associated with euphotic MGII during 2006 -2007, and also at subsurface during 2008, supported by previous findings of a wider distribution and diversity within Euryarchaeota as in other coastal areas, e.g., SPOT (Parada and Fuhrman, 2017).

The Spatiotemporal Distribution of Thaumarchaeota Cell Abundances
Unlike Euryarchaeaota, thaumarchaeal cells were found to be favored in subsurface waters, characterized by peaks at the oxycline during fall and spring (2007) and maxima at the bottom during wintertime and late summer 2008. Thaumarchaeal cells made up the bulk fraction of the prokaryote community at deep waters and represented a large fraction of the archaea considering the sum of Euryarchaeota and Thaumarchaeota (Supplementary Figure 3). These results are supported by previous findings in the Eastern South Pacific (ESP) off Chile, based on rRNA dot blot hybridization (Levipan et al., 2007b;Quiñones et al., 2009) and lipid-based approach (Rossi et al., 2012;Srain et al., 2015). GALMSS analysis indicates that Thaumarchaeota was the microbial group with the greater percentage of predictability associated with its spatial and temporal variability (Table 1), and with all the other variables (oxygen, nitrate, ammonium and nitrate), associated with substrates and products of nitrification. In general, these results support a recent report considering the water column position (euphotic, 50 m, and deeper waters 60 -500 m) and other factors such as hydrographic explanatory variables, including oxygen, nitrite and ammonium Thaumarchaeota contribution to total microbial community (iTag-16S rRNA) variability in the Monterey Bay coastal time series . However, in the case of the time series off Concepción the Thaumarchaeota depth distribution is resolved within a shallow water column depth (92 m depth).

Contribution of Ammonia Oxidizing Archaea in the Study Area
The AOA amoA abundance results indicate that AOA in the study area reached comparable abundances with those determined in other marine ecosystems, e.g., the northern Gulf of Mexico (up to 10 5 copies mL −1 at ≤100 m depth; Tolar et al., 2013), Sargasso Sea or ETNA (up to 10 5 copies mL −1 , Löscher et al., 2012;Newell et al., 2013), North Sea in winter months (from 10 4 up to 10 5 copies mL −1 , Wuchter et al., 2006;Herfort et al., 2007;Pitcher et al., 2011b), and within the oxycline and just over the upper limit of the OMZ in the Arabian Sea (up to 10 5 copies mL −1 , Newell et al., 2011;Pitcher et al., 2011a;Bouskill et al., 2012). However, variable ratios between amoA and 16S rRNA genes or cell counts have been observed in marine ecosystems (e.g., Wuchter et al., 2006;Mincer et al., 2007;Agogué et al., 2008). In our study, a variable contribution of AOA amoA versus Thaumarchaeota cell counts determined through CARD-FISH was observed, reaching higher ratios during spring and fall (0.2 -0.7) compared with <0.2 ratios observed during summer and winter (Supplementary Figure 4). Similar values were obtained in surface waters from Monterey Bay (Mincer et al., 2007) and in subsurface (100 -150 m depth) from subtropical and equatorial regions of the North Atlantic (Agogué et al., 2008). However, lower ratios of AOA amoA to Thaumarchaeota 16S rRNA gene ratios (<0.01) were associated with the presence of non-nitrifying Thaumarchaea from bathypelagic waters (Agogué et al., 2008) or with methodological biases (Mincer et al., 2007). Alternatively, different ecotypes could account for a wide range of AOA amoA versus Thaumarchaeota 16S rRNA gene ratios (e.g., Sintes et al., 2013). GAMLSS modeling indicates that similar predictive variables account for both Thaumarchaeota cell and AOA amoA fluctuation ( Table 1).

Contribution of Ammonia Oxidizing Archaea Ecotypes in a Representative
Upwelling Versus Non-upwelling Season Profile AOA in the study area were found to be associated to the candidate order Nitrosopumilales (Stieglmeier et al., 2014), and to other ecotypes, i.e., Water Column Surface clade (WCA) with one culture Candidate Nitrosopelagicus brevis and Water Column Deep clades (WCB), as in other coastal marine time series (e.g., Beman et al., 2008;Santoro et al., 2010Santoro et al., , 2015Reji et al., 2019;Tolar et al., 2020). The ammonia-oxidizing assemblages were predominantly characterized by N. maritimus-like followed by WCA amoA reaching higher abundances during winter, whereas the deep WCB amoA ecotype showing maxima during summer at the bottom. This result was supported by previous reports in the study area based on AOA amoA survey using clone libraries (Molina et al., 2010), and iTag sequencing (Bertagnolli and Ulloa, 2017) using primers designed by Pester et al. (2012). These previous studies revealed specific predominant contribution of the WCA amoA in both upwelling and non-upwelling season and of the deep WCB amoA during the upwelling period. The WCB was strongly coupled to oxygen deficient conditions, and was suggested to be advected by upwelling waters toward the continental shelf (Bertagnolli and Ulloa, 2017). Indeed, WCA and N. maritimus amoA ecotypes are the main AOA in the study area in terms of the abundances of amoA genes and transcripts, and related to ammonium oxidation (Molina et al., 2010;Bristow et al., 2016;Levipan et al., 2016). In addition, both archaeal ecotypes, but mainly WCA, appear to have an important role in the ammonium oxidation in other marine ecosystems; e.g., the Gulf of California (Beman et al., 2008), Monterey Bay (Smith et al., 2014;Tolar et al., 2020), from the equatorial Pacific to the Arctic (Shiozaki et al., 2016), and the suboxic zone of the central Baltic Sea (Labrenz et al., 2010). Nonetheless, N. maritimus-like ecotypes may be less abundant at lower latitudes of the ETSP (∼12 and 18 • S; Peng et al., 2013).

Ammonia Oxidizing Bacteria and Its Interaction With Their Archaea Counterparts in the Study Area
The abundance of bAOB amoA ranged from the limit of detection (2 copies) to 10 4 copies mL −1 (Figures 2, 4), reaching higher numbers mainly during spring-summer at intermediate depths. This range was lower than that previously described in the ETSP (10 3 to 10 5 copies mL −1 , Bouskill et al., 2012;Löscher et al., 2012), similar to the range reported by Wuchter et al. (2006) in the North Sea, and wider than that found in locations such as Monterey Bay (undetectable to 10 3 copies mL −1 , Mincer et al., 2007), shallow waters in the California Current (undetectable to 10 2 copies mL −1 , Santoro et al., 2010Santoro et al., , 2013, Gulf of California or ETNP (undetectable to 10 2 copies mL −1 , Beman et al., 2008Beman et al., , 2012, and within the oxycline in the Arabian Sea (from ∼10 2 to 10 3 copies mL −1 , Newell et al., 2011;Bouskill et al., 2012).
Spatiotemporal followed by nitrate and ammonium were significant explanatory variables for bAOB amoA, distribution. In addition, the variability in bAOB amoA appears to be constrained by AOA distribution. Shifts in bAOB amoA were more obvious when plotted against the AOA amoA versus Thaumarchaeota cell ratios estimated from CARD-FISH counts (Supplementary Figure 5B). These results show that the number of bAOB amoA was maximum in surface waters during summerfall when low AOA amoA versus Thaumarchaeota cell ratios were found, and reached a plateau at spring-time. This finding suggests that competitive interactions between AOA and bAOB amoA are more intense in surface waters during springtime and were potentially associated with lower ammonium availability and/or oxygen deficiency. The numerical dominance of AOA amoA in the study area could be attributed to the wellknown kinetic differences in ammonium uptake between AOA and bAOB (Martens-Habbena et al., 2009). It is known that ammonium concentration promotes niche separation not only between phylogenetically distant species of AOA and AOB (Schleper, 2010), but also between different clusters of AOA (Sintes et al., 2016). In the present study, ammonia-oxidizing assemblages significantly varied with ammonium (model 5 in Table 1), suggesting niche overlapping between AOA and AOB and hence substrate competition for ammonium. Nevertheless, in contrast to Fernandez and Farías (2012), who reported that bacterial Ammonium oxidation might benefit from low oxygen concentrations, we found that the oxygen concentration was a statistically meaningful predictor for AOA amoA but not for bAOB amoA distribution (models 3 in Table 1), suggesting that oxygen concentration can be involved in niche differentiation between AOA and bAOB amoA (models 1 and 3 in Table 1). Moreover, our results support the key role that oxygen concentration plays in AOA-A and B ecotypes partitioning in this and other areas (Beman et al., 2008;Bertagnolli and Ulloa, 2017), as well as in controlling (at low concentrations) the distribution and activity of marine AOA populations (Qin et al., 2017).
Ammonia oxidizing assemblages are significant contributors to biogeochemical cycling in the study area associated with high rates of nitrification that are correlated with the availability of ammonium mainly in the oxyclines in the study area (Fernandez and Farías, 2012). Moreover, in a 10-year (2002-2012) biogeochemical report associated with nitrous oxide dynamics the COPAS time series, nitrous oxide positive excess, mainly in the oxyclines and subsurface waters, was associated with nitrification (Farías et al., 2015;Galán et al., 2017), as in other oxygen deficient areas of the Pacific Ocean, e.g., Löscher et al. (2012) and Frey et al. (2020). In order to compare potential role of microbial communities studied here, the average monthly inventories of nitrous oxide reported by Farías et al. (2015) for the 2005 -2009 upwelling periods were compared with nitrate, ammonium, AOA, bAOB and non-fluorescent picoplanktonic community estimated inventories in our study (Supplementary Figure 6 and Table 2). This comparison helps us identify the connection between ammonia oxidizing assemblages with Nitrous oxide and nitrate mainly during 2007 -2009, but not for the previous upwelling periods. During 2005During -2007, nitrous oxide and nitrate inventories could be associated with high ammonium accumulation and an enrichment of picoplanktonic communities with unknown identity in the bottom (Supplementary Figure 6). Therefore, changes in microbial assemblages are associated with upwelling variability, generating significant impact in the water column biogeochemistry in this natural oxygen -deficient upwelling ecosystem.

CONCLUSION
In total, our results indicate the physical-chemical oceanographic conditions associated with the seasonal upwelling systems modulate the distribution and dynamics of abundant marine microbial groups including key ammonia oxidizing assemblages. The statistical GAMLSS analyses indicate that spatiotemporal changes associated with biogeochemical conditions trigger by oxygen-deficiency, including substrate and products of nitrification, were significant predictors shaping the marine microbial community distribution. Biological interactions including AOA and bAOB, and possibly with other microbial groups in surface waters influences the temporal dynamics of ammonia oxidizing assemblages. Considering the key role that coastal marine microbial community plays in the cycle and the balance of nutrients and carbon at a global scale, an appropriate understanding of their natural variability will allow to better predict the metabolic and biogeochemical consequences associated with the expected expansion and intensification of oxygen deficient zones due to anthropogenic forcing.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article are given in Supplementary Table 2 and complementary data will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
VM, LB, and HL wrote the manuscript with the support from SR-F, CA, AG, IM, and OU. VM and LB conceived the original idea. VM and OU supervised the project. All the authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We thank the time series station program of the COPAS Center of the University of Concepción for providing environmental data. We acknowledge the help provided by Kay Kay crew and technical support of S. Collado, M. Montoya, F. Santibañez, R. De la Iglesia, and G. Alarcon at the field and laboratory.