Dynamic prokaryotic communities in the dark western Mediterranean Sea

Dark ocean microbial dynamics are fundamental to understand ecosystem metabolism and ocean biogeochemical processes. Yet, the ecological response of deep ocean communities to environmental perturbations remains largely unknown. Temporal and spatial dynamics of the meso- and bathypelagic prokaryotic communities were assessed throughout a 2-year seasonal sampling across the western Mediterranean Sea. A common pattern of prokaryotic communities’ depth stratification was observed across the different regions and throughout the seasons. However, sporadic and drastic alterations of the community composition and diversity occurred either at specific water masses or throughout the aphotic zone and at a basin scale. Environmental changes resulted in a major increase in the abundance of rare or low abundant phylotypes and a profound change of the community composition. Our study evidences the temporal dynamism of dark ocean prokaryotic communities, exhibiting long periods of stability but also drastic changes, with implications in community metabolism and carbon fluxes. Taken together, the results highlight the importance of monitoring the temporal patterns of dark ocean prokaryotic communities.

A total of 6570 ASVs (amplicon sequence variants) were identified in the whole dataset before rarefaction. 29, 20 and 24% of the ASVs were found only in the LIW, oWMDW and bottom water, respectively, and 9.5% of ASVs were shared between the three water masses (Supplementary Fig. S4). These shared ASVs contributed 73, 82 and 74% in abundance at LIW, oWMDW and bottom water, respectively ( Supplementary Fig. S4). The variation in the contribution of the most abundant prokaryotes at family level between water masses is shown in Fig. 2. Nitrosopumilaceae was the most abundant phylotype at the family level in the meso-and bathypelagic western Mediterranean Sea waters, decreasing its contribution to the community with depth from 40.6% in the LIW to 39.2% in the oWMDW and 35.3% at the bottom water (Fig. 2a). Therefore, Nitrosopumilaceae abundance decreased with depth from 5 ± 0.3 × 10 4 cells mL −1 (mean ± S.E.M) at LIW to 2.5 ± 0.1 and 1.9 ± 0.1 × 10 4 cells mL −1 at the oWMDW and bottom, respectively. SAR11 clade I contributed similarly (5.5-8.8%) than clade II (7-8.7%) to the prokaryotic communities, decreasing their abundance with depth. SAR11 clade I abundance   (Fig. 2b). Families of Subgroup 6 (Acidobacteria), of Rhodospirillales and SAR86 (Gammaproteobacteria) contributed > 1% to the community, however, they did not significantly differ between water masses. Low abundance family phylotypes (i.e., contributing < 1% relative abundance to the community) showed less pronounced changes in relative abundance with depth than the abundant phylotypes (> 1% relative abundance) (Fig. 2). Principal coordinate analysis based on Bray Curtis (considering ASVs abundance) showed clustering of communities according to water masses and cruise (Fig. 3a). These variables explained 36 and 8.3% of community variability (PERMANOVA), respectively. Weighted UniFrac distance metric (considering phylogenetic distance and ASVs abundance) explained a larger fraction of the prokaryotic community variability (58%) than Bray Curtis distance metric (45%). Weighted UniFrac analysis did not show a clear clustering of communities according to water masses, and discrete samples of the three water masses deviated from the general pattern, strongly differing from the rest (Fig. 3b). Cruise, water mass and station explained 14, 8.6 and 7.5% of variance in community composition (PERMANOVA), respectively.
Temporal dynamics of meso-and bathypelagic ocean prokaryotes. The weighted UniFrac distance between the communities present in different dates and the initial community (Feb16) was computed for each station and water mass in order to quantify the community structure variation along the studied period. The results show different variation patterns across water masses and basins (Fig. 4). Communities that exhibited larger distance throughout the study period (i.e., oWMDW in station B at the Mallorca channel, LIW community in station C at the south Balearic sub-basin and the communities from the three water masses in station E at the south Algerian sub-basin, Fig. 4) resulted from strikingly different community composition in Feb16 compared to other months ( Fig. 5 and Supplementary Fig. S5). The community composition showed periods of change, depicted by increasing distance to the initial community, followed by the return to the initial community composition, represented by decreasing distances (Fig. 4). Notably, changes were observed at all water masses, including the bottom waters (Fig. 4).
Community composition from station C oWMDW and bottom water was similar between Feb16 (winter) and Apr16 (spring) and significantly diverged towards Jul16 (summer) ( Fig. 4 and Supplementary Table S2), subsequently returning to a composition similar to Feb16 in Oct16 (autumn) for the oWMDW community and in Feb17 (winter) for the bottom community, and remained similar during Jun17 and Nov17 (Fig. 4c). Similarly, the communities of station D exhibited a pronounced change in Oct16 (autumn) at the three water masses, followed by the return to the initial community composition in Feb17 (winter) for LIW and oWMDW communities and in Jun17 (summer) for the bottom community ( Fig. 4d and Supplementary Table S2). The three water masses in station E experienced a similar temporal patterns, with the more pronounced changes occurring in Apr16 (spring) and Jul16 (summer) as compared to Feb16 (winter), and communities stabilizing afterwards (Fig. 4).
In order to explore whether the observed changes in communities were due to the abundance variation of ASVs always present in a specific community or due to the transport of ASVs from a different community and subsequent replacement in the studied community, the ASVs were categorized as 'core' (i.e., present in ≥ 80% of samples), 'frequent' (i.e., present in ≥ 50% and < 80% of samples) and 'transient' (i.e., present in < 50% of samples). Core ASVs were assigned to few taxa, mainly Nitrosopumilaceae and SAR11 (Supplementary Fig. S6 and Table S3). One ASV was always abundant (contributing ≥ 1% to all samples), putatively assigned to Nitrosopumilaceae. Seven additional core ASVs were classified as frequently abundant, i.e., contributing ≥ 1% in ≥ 50% of samples. From these, five were assigned to Nitrosopumilaceae, one to SAR11 Clade II and one to order SAR324   Fig. S6), e.g., at stations C, D and E in Jul16, Oct16 and Feb16, respectively ( Fig. 5 and Supplementary Figs. S5 and S5). The family Moraxellaceae (Pseudomonadales), specifically genus Psychrobacter, steeply increased throughout the three water masses in Jul16 (summer) at station C (south Balearic sub-basin), contributing up to 60.4%, 30.5% and 15% at the bottom, oWMDW and LIW communities, respectively ( Fig. 5a-c). Concurrently, diversity indexes of bottom water communities decreased pronouncedly in Jul16 ( Supplementary Fig. S7) and the contribution of HNA cells increased throughout the water column, more pronouncedly at bottom waters ( Supplementary Fig. S2). During the other cruises, this family was not detected at the bottom and was considered rare for the overall water column (contributing ≤ 0.1% in abundance), only it was detected contributing more than 0.1% and less than 1% (i.e., categorized as 'low abundant') to the community of LIW in Oct16 and of oWMDW in Apr16. In terms of absolute abundance, this family was present at abundances of 0-4, 0-5 and 0-9 × 10 2 cells mL −   www.nature.com/scientificreports/  Table S2), was principally caused by a marked increase of Flavobacteriaceae (16.6-46.9%) with respect to the other cruises (0-0.9%) ( Fig. 5d-f). This family amounted to 20-85, 1.5-17 and 0-7.6 × 10 1 cells mL −1 in LIW, oWMDW and bottom communities in other cruises, respectively, and reached up to 2.8, 3.2 and 0.9 × 10 4 cells mL −1 in Oct16, respectively. Besides, Sphingomonadaceae and some low abundant Gammaproteobacteria families increased to 11% and 13% in bottom waters, respectively, with respect to 0-10% and 2-3% in other cruises (Fig. 5f). The large increase in contribution of these ASVs resulted in decreased diversity indexes at all depths ( Supplementary Fig. S7).
Communities of station E (south Algerian sub-basin) in Feb16 (winter) were characterized by a substantial contribution of Gammaproteobacteria families to the three water masses communities (33.2-77.7%). In particular, Alteromonadales contributed 28.4% at oWMDW communities, and pronouncedly decreased its contribution in the following cruises to < 5% ( Supplementary Fig. S5). A large Gammaproteobacteria families contribution in station E was still observed in Apr16 (spring) at oWMDW (36.2%), as compared to the following cruises (6.7-15.9%) (Supplementary Fig. S5).
Although the most pronounced community changes observed mainly corresponded to an increase in the contribution of transient taxa, sometimes the increase in contribution of core taxa also caused community shifts. The bottom community changes at station A in Oct16 (autumn) and station D in Feb17 (winter) ( Supplementary  Fig. S6), located at the north of the Balearic and Algerian sub-basins, respectively, were related to changes in core ASVs from the family Nitrosopumilaceae. Members of this family contributed up to 63% at the bottom community of station A in Oct16 as compared to 28-41% in other cruises (Fig. 5). However, no significant change in its absolute abundance was observed in this station. This family also increased its contribution to bottom water communities in station D to 54% in Feb17 as compared to 24-38% in other cruises ( Supplementary Fig. S5). Concurrently, an increase in absolute abundance of Nitrosopumilaceae occurred in this station, reaching 3.6 × 10 4 cells mL −1 in Feb17 as compared to 1-2 × 10 4 cells mL −1 in other cruises.
Prokaryotic community variation following a ventilation event. We aimed at exploring whether the ventilation event in Apr16 inferred from the oxygen profiles in stations C, D and E caused significant variations in particular taxa differential abundance. In order to do that, we used the LEfSe pipeline 43 . This pipeline allows to associate taxonomic groups with particular environmental conditions, such as oxygen concentration levels, at multiple taxonomic levels. The range of oxygen values reported in this study (between 164.2 and 201 μmol O 2 kg −1 ) are characteristic of oxygenated ocean waters. Therefore, we used oxygen concentration as a proxy of water mass variabilities (such as the bottom ventilation event or dissolved oxygen concentration decrease at LIW) and not as a driving parameter of prokaryotic community changes. Consequently, three categories of oxygen concentration were defined to encompass the variabilities in water mass properties: 'low oxygen' , which includes the lowest oxygen concentration values (164-176 μmol O 2 kg −1 ) at LIW to characterize variations in this water mass; 'high oxygen' , which includes the highest oxygen concentration values at the bottom (> 197 μmol O 2 kg −1 ), corresponding to the bottom ventilation observed in three of the stations, and 'mid oxygen' , corresponding to the values between 'low' and 'high' oxygen concentrations.
The differential abundance of taxonomic groups was related to phylotypes from diverse taxonomic levels and phylogenies: 4 phylum, 11 classes, 20 orders and 26 families were significantly related to low or high oxygen categories ( Fig. 6 and Supplementary Fig. S10). Phylotypes associated to high oxygen concentration, i.e., increased their contribution to the community when the bottom ventilation event occurred, at phylum level were Crenarchaeota, Nanoarchaeota, Acidobacteria and Thermus. Classes associated to high oxygen were one unidentified Crenarchaeota class, Thermococci, Subgroups 5, 11 and 21 (Acidobacteria), Anaerolineae, TK10 (Chloroflexi) and Deinococci (Fig. 6). Other related taxa at order and family levels were the orders Aigarchaeales, Methanofastidiosales, Deinococcales and 10 uncultured orders and families of Acidobacteria.
Low oxygen category used in the analysis was defined as the lower oxygen concentration values recorded at LIW throughout the study. Low oxygen associated taxa were the uncultured family FIN625 of Marine Group II, one class, order and family of SAR406 (Marinimicrobia), and the uncultured family Arctic97B of Verrucomicrobia (Fig. 6). No taxon was related to mid oxygen category.

Discussion
Despite the interest in exploring the microbial communities of the deep ocean and their role in global biogeochemical cycles, knowledge on their temporal variability remains scarce. Here, by performing a spatio-temporal survey of meso-and bathypelagic prokaryotic communities, we show that dark ocean prokaryotic communities were temporally dynamic, changing severely their structure. Our results revealed that low abundant (< 1% relative abundance) and rare (< 0.1% relative abundance) phylotypes can increase drastically after a disturbance, suggesting a rapid response of active prokaryotes to changing environmental conditions. Three groups dominate the aphotic water column in the Mediterranean: members of SAR11 clades I and II and the archaeal family Nitrosopumilaceae (Fig. 2), in agreement with their widespread distribution in the dark ocean 44,45 . These three families together comprise ~ 50% of the total dark Mediterranean Sea prokaryotic community, thus, their metabolic activity may profoundly influence the deep ocean biogeochemical processes. The dominant family in the meso-and bathypelagic Mediterranean Sea, the Nitrosopumilaceae (Fig. 2), are ammonia-oxidizing chemoautotrophs 46 . Deep ocean environments are characterized by very low ammonium concentrations, often below detection limit of current techniques 47 . The archaeal ammonia oxidizers that inhabit this environment are suggested to be adapted to different ammonium concentration and supply rates [48][49][50] , and to have the potential to use additional energy sources 49,51 . The significant contribution of SAR11 members in the aphotic water masses (Fig. 2) fits with their heterotrophic oligotroph lifestyle harbouring a plethora of metabolic capabilities 52 , such as preferential use of organic sulphur compounds 53,54 , which requires a lower energy cost than producing reduced sulphur from sulphate 55 . SAR11 is the most abundant bacterial group in the ocean with depth-specialized phylotypes and ecotypes 52,56 . Families from SAR324 and SAR406 (Marinimicrobia) clades Figure 6. Differential prokaryotic taxa associated to changes in water mass properties. Cladogram showing the phylogenetic distribution of prokaryotic taxa with significantly differential abundance at low or high oxygen concentration computed for the different taxonomic levels. Rings indicate taxonomic levels from domain (inner circle) to family (outer circle). Circle and shaded area colours indicate taxa associated to low (green) or high (red) oxygen concentration for a linear discriminant analysis (LDA) score > 2. Taxa with no significant differences between different oxygen categories are represented by yellow circles. Legend shows the taxonomic level of taxa (p: phylum, c: class, o: order, f: family). www.nature.com/scientificreports/ also contribute a considerable fraction to the aphotic zone communities (Fig. 2), and are frequently found in other deep ocean waters 12,57 . SAR324 members have been also identified as chemoautotrophs based on sulphur oxidation metabolism 56 . SAR406 contribution to total prokaryotes increased with depth, in agreement with previous studies in the Atlantic Ocean 58,59 . This clade is associated with a particle attached lifestyle 60 , reported to contribute significantly to dark dissolved inorganic carbon fixation 61 and to participate in sulphur and nitrogen biogeochemical cycles 62 . The LIW exhibits the minimum oxygen concentration of the Mediterranean Sea water column (164-185 μmol kg −1 in this study), and is described as the oldest water mass due to the slow spreading and long renewal time 63 . The reduced oxygen conditions favour some chemolithoautotrophic processes such as ammonia oxidation due to a larger availability of reduced chemical compounds compared to more oxygenated waters 64 , supporting the larger relative abundance of Nitrosopumilaceae observed at the LIW (Fig. 2). The larger contribution of Nitrosopumilaceae at this depth coincides with an increase in the contribution of the nitrite oxidizer family Nitrospinaceae to the prokaryotic community (Fig. 2), supporting the suggested reciprocal feeding between nitrite oxidizers and archaeal ammonia oxidizers 65 and indicating an important role of the prokaryotic community inhabiting this water mass in the nitrogen cycle of the Mediterranean Sea. The oWMDW is the main bathypelagic Mediterranean water, occupying most of the water column. It is a rather hydrographical stable water mass as compared to LIW and exhibits transitional environmental and community characteristics between the mesopelagic and the new bottom water mass. The decrease in the contribution of particular LIW phylotypes to the community (i.e., Nitrosopumilaceae, Nitrospinaceae, Thioglobaceae and the Marine Group II family FIN625, among others) coincides with the relative increase of heterotrophic oligotrophs such as specific uncultured Marine Group II phylotypes and SAR11 clade II (Fig. 2). The higher diversity and variability of nearbottom communities as compared to the other two water masses studied here probably relates to differences in physicochemical properties and to the influence of resuspension of sediments of different characteristics in the different stations 66,67 . Organic and inorganic compounds leakage from the sediments may support a wider variety of metabolic processes in the near-bottom waters compared to the water masses above 68 . However, it should also be noted that the stations depth range is variable and that the stations C, D and E exhibit an additional water mass (Western Mediterranean Transition) at the bottom. SAR202 larger contribution in deeper waters is supported by their role in the final stages of organic matter oxidation 69 and the reported major contribution of refractory organic matter to the dissolved organic carbon pool below 1000 m depth 70 .
Different major community composition changes occurred in all locations spanning the entire aphotic water column (Figs. 4, 5 and Supplementary Fig. S5), with different members of the prokaryotic community leading the changes. Most phylotypes that experienced substantial changes in their abundance and contribution over the studied period were rare or low abundant in the other communities and associated to the transient ASVs community (i.e., present in < 50% of samples). Experiments with natural communities have shown the persistence of rare deep ocean prokaryotes during long periods of starvation and their feast response to a pulse of carbon 41,71 . However, our findings also suggest that mostly a transport and replacement of the local community by a different origin's community rather than a change of the core community occurred, in agreement with previous observations for surface prokaryotic communities 72 . However, sequencing constraints regarding the detection of rare sequences have to be considered, e.g., sequencing depth might hinder the detection of some rare but frequently present ASVs.
Long-term variations and seasonal patterns have been well studied for surface microbial communities, linked to seasonal changes of environmental factors 48,73,74 . However, temporal variations of deep ocean microbial communities are by far less studied. Seasonality and temporal dynamics of meso-and bathypelagic microbial communities have been linked to surface ocean processes 13,14,16 . In this study, we identified a local change in the LIW community potentially linked to particle fluxes from the epipelagic. The change of the LIW community in winter 2016 (Feb16) at station C was associated to some phototrophic phylotypes (Chloroplast, Cyanobiaceae), Flavobacteriaceae and SAR116, reported to contribute more to photic zone communities 75,76 , indicating that epipelagic prokaryotes reached the intermediate waters. A drastic composition shift was detected in station E in winter 2016 (Feb16) as compared to other sampling dates, associated to a high contribution of Gammaproteobacteria phylotypes, mainly assigned to Alteromonadales order ( Supplementary Fig. S5). This station is located at the eastern Alboran Sea, a region characterized by intense hydrological dynamics, including upwelling events and frontal systems, that influence the temporal variability of particle fluxes 77 , likely impacting on prokaryotic communities. The massive contribution of the rare genus Psychrobacter (Moraxellaceae) in summer 2016 (Jul16) at the bottom of station C and its decreasing contribution towards upper layers suggests a community change caused by a local environmental disturbance. The increase of HNA cells at the deep-water masses together with the increase of turbidity and nitrite concentration, particularly at the bottom, supports the environmental and community disturbance at station C in July 2016 (Supplementary Figs. S7, S8 and S9). The most pronounced community composition shift in station D was detected in autumn 2016 (Oct16) with several taxa responding (e.g., Flavobacteriaceae, Sphingomonadaceae or Gammaproteobacteria phylotypes, Fig. 5), suggesting a disturbance of different origin and that local factors might influence the community changes.
The variability of surface biotic processes influencing the formation of organic particles could explain the community changes observed in bathypelagic waters 13 . The major community changes detected, although observed throughout the three water masses, are more pronounced in the oWMDW and bottom communities. This finding would be in agreement with a direct connectivity between the surface and the bathypelagic through fast sinking particles that elude mesopelagic communities, as has been previously suggested 14 . The fast sinking of particles 14 , together with the physical and hydrochemical particularities of the LIW 63 , support the decoupling of the mesopelagic waters with the above and below waters. Other common processes in the area likely influencing water masses connectivity and community changes are mesoscale eddies, formation of intermediate nepheloid layers or downwelling from the shelf-slope 38 www.nature.com/scientificreports/ continental slopes. Bottom trapped waves 38 are occasional phenomena occurring in the area of study that could alter sediment resuspension or enhance material advection from surrounding areas, influencing near-bottom or aphotic zone prokaryotic communities locally. Pasqual et al. 67 reported differences in sedimentation regimes and organic matter origin and composition between two locations in the same study area. Unfortunately, surface biotic processes or specific physical forcing events were not assessed in this study.
In addition, this study shows a prokaryotic community shift related to a bottom ventilation event. To our knowledge, this is the first evidence of near-bottom community changes at basin scale. Open-ocean convective processes and dense shelf water cascading events during winter leading to water ventilation in near to bottom waters are frequently reported in the north-western Mediterranean 36,79 . The newly formed deep water flows from the north-western Mediterranean towards the Atlantic following a cyclonic path. However, Margirier et al. 80 did not report a bottom-reaching water convection in the north-western Mediterranean in 2016 and we lack data to comprehensively relate it with a specific basin event. Further support of the event is provided by the reported association to sediments, benthic communities or particle attachment lifestyle of phylotypes related to the bottom ventilation 57,81 , suggesting that the community change could be linked to resuspension of sediments.
Taken together, our results indicate that the meso-and bathypelagic ocean prokaryotic communities exhibit long periods of stability (e.g., from February to November 2017 in most stations; Figs. 4, 5 and Supplementary  Fig. S5) but also experience drastic and occasional changes. However, it should be mentioned that given our sampling resolution (every 2-3 months), the exact time when the shifts occurred and its extension cannot be precisely described. Regardless of the microbial shifts detected, the community structure returned to its original composition, strengthening the view of an active and dynamic deep ocean community. Although logistically challenging, further studies should focus on the characterization of temporal variations of deep ocean prokaryotic communities at smaller temporal scales (e.g., daily or weekly) to better understand the temporal succession of community changes and to address the ecologically relevant question of resistance, resilience or sensitivity of these communities. Our results also suggest that the blooming of otherwise rare or absent phylotypes, such as Moraxellaceae, is originated by transient phylotypes. However, this interpretation has to be taken with care due to the limitations of sequencing (e.g., sequencing depth), that might prevent the detection of some rare phylotypes from the core community, frequently labelled as microbial seed bank 82,83 . The microbial seed bank provides biodiversity and metabolic flexibility to the community, allowing a relatively rapid response to environmental changes and contributing to the community resistance or resilience 41,84,85 . Further research on the rare and low abundant phylotypes (e.g., by increasing sequencing depth) and on their metabolic potential and activities will help to better understand their role under environmental disturbances and in the biogeochemical cycles.  (Fig. 1a). Temperature and salinity depth profiles were recorded with a SBE911 Conductivity-Temperature-Depth (CTD) sensor equipped with oxygen (SBE43) and turbidity (SeaPoint) sensors. Calibration procedures were applied according to standard protocols 86 .Salinity measurements were calibrated using Portasal Guildline 8410A. Oxygen concentrations obtained with the SBE43 sensor were calibrated for each cruise based on the comparison to Winkler titration measurements 87 . Turbidity values were normalized for each cruise by subtracting the lowest FTU (Formazin Turbidity Unit) value determined during the cruise.

Methods
Prokaryotic community analyses were conducted on water samples collected from the meso-and bathypelagic water masses identified: the Levantine intermediate water (LIW), old western Mediterranean deep water (oWMDW, collected at 1000 m) and bottom (collected at 5-10 m above the seafloor). The core depth of the LIW was identified based on its signature in a T-S diagram (Fig. 1b) during the CTD downcast. Additional samples for inorganic nutrients and prokaryotic abundance were taken at 700 and 1500 m. A total of 93 samples were collected for community composition analysis and 168 samples for inorganic nutrients and prokaryotic abundance analysis (Supplementary Fig. S1).
Water samples for dissolved inorganic nutrients were collected in 12 mL vials and kept frozen at -20ºC until further processing. Nitrate (NO 3 − ), nitrite (NO 2 − ), phosphate (PO 4 3 ) and silicate (SiO 4 2 ) concentrations were determined using a QuAAtro Gas Segmented Continuous Flow Analyser (SEAL Analytical) following colorimetric methods 87-89 . Prokaryotic abundance. Seawater samples (1.5 mL) were fixed 10 min with glutaraldehyde (0.1% final concentration), frozen in liquid nitrogen and stored at -80ºC until analysis at the home laboratory. Prior to analysis, samples were thawed and stained with SYBR Green I (Sigma-Aldrich, 1 × final concentration) for 10 min in the dark. Prokaryotic cells were counted on an ACCURI C6 flow cytometer (BD Biosciences) based on their signature in a cytogram of side-scatter versus green fluorescence. High nucleic acid (HNA) and low nucleic acid (LNA) populations were distinguished by gating based on their relative green fluorescence signal. Fluores- 16S rRNA gene amplification was performed using the primers 515F-Y (5′-GTG YCA GCMGCC GCG GTAA) and 926R (5′-CCG YCA ATTYMTTT RAG TTT) targeting the hypervariable regions V4 and V5 of both Archaea and Bacteria 90 . PCR reaction mixture contained 1 µL of extracted DNA, forward and reverse primers (1 µM final concentration) and 1 × Kapa HiFi HotStart ReadyMix (Kapa Biosystems, Wilmington, MA). Cycling conditions were applied as described by Parada et al. 90 . PCR products were purified using the Qiaquick PCR purification kit (Qiagen, Hilden), following the manufacturer's protocol. PCR-purified amplicons were checked on a 2% agarose gel and sequenced using Illumina MiSeq 2 × 250 bp with V2 chemistry, resulting in 2,860,134 total reads.
Sequence data was analysed using DADA2 pipeline 91 implemented in QIIME2 (http:// qiime2. org) 92 . Pairedend sequence reads were first demultiplexed and quality filtered, applying the following specific parameters: maximum number of expected errors maxEE = 2, removing all sequences with ambiguities (maxN = 0), truncating the length (truncLen) of forward and reverse reads to 245 bp and 230 bp, respectively, and trimming the primers length (trimLeft). Chimeras were removed using the 'denoise-paired' script implemented in QIIME2 with default settings. After quality check, sequences were clustered in amplicon sequence variants (ASVs) 93 . ASVs alignment was performed using MAFFT 94 and used to infer unrooted and rooted phylogenetic trees with Fasttree 95 . ASVs taxonomy was assigned using the SILVA database (release 132). Phylotypes contributing ≤ 0.1% in abundance to the community were considered 'rare' 42 , and phylotypes with relative abundance < 1% were considered 'low abundant' . The absolute abundance of phylogenetic groups was calculated multiplying the total prokaryotic abundance, measured by flow cytometry, by the percentage of contribution of each group to the community.
Diversity metrics and statistics. Shannon diversity, Pielou's evenness and phylogenetic diversity were calculated with QIIME2 using the ASVs table rarefied to 3007 ASVs. ASVs beta diversity was analysed using Bray Curtis and weighted UniFrac distance matrices. Principal Coordinates Analysis (PCoA) was used to represent the variance among communities according to the distance matrices previously mentioned. Statistically significant environmental variables were assessed by analysis of variance (ANOVA) and permutational multivariate analysis of variance (PERMANOVA) using the 'aov' and 'adonis' function of the 'vegan' package of R. A significance level of P < 0.05 based on 999 permutations was used. The 'ampvis2' package of R was used to determine the presence and frequency of ASVs 96 for the whole dataset (93 samples) in the region studied. ASVs consistently reported in ≥ 80% of samples were considered 'core' ASVs, ASVs present in 50-79% of samples were considered 'frequent' and ASVs present in < 50% of samples were considered 'transient' 97 . The weighted UniFrac distance was used to assess the change in community composition throughout the study period and to calculate the distance between community composition of adjacent sampling dates. Then, one-sample t-tests were used to statistically determine the significance of a community change as compared to the mean community value for each water mass and station. Oxygen concentration was used as a proxy of water mass variability to relate changes in water masses properties with prokaryotic community changes. Temporal variations in the local oxygen minimum layer (i.e., at LIW) and bottom ventilation events are indicated by oxygen fluctuations 98 . Linear discriminant analysis (LDA) combined with effect size measurements (LEfSe) was used to associate particular taxonomic groups with fluctuations in these water masses using changes in oxygen concentration at multiple taxonomic levels (http:// hutte nhower. sph. harva rd. edu/ galaxy) 43 . Statistical analysis was computed from phylum to family level to simplify computation and results. A significance level of P < 0.01 was used in Kruskal-Wallis test. Prokaryotic communities were classified in three categories according to oxygen concentration: 'Low_O2' for samples with O 2 < 176 μmol kg −1 , corresponding to the LIW samples with lowest oxygen concentrations; 'Mid_O2' corresponding to samples with oxygen concentration ranging between 176 ≥ O 2 ≤ 197 μmol kg −1 ; and 'High_O2' when O 2 > 197 μmol kg −1 , concentrations determined when a bottom ventilation event was observed.

Data availability
Sequences reported in this paper have been deposited as 16S rRNA raw sequences in NCBI Sequence Read Archive (SRA) under the accession number PRJNA575848.