Niche partitioning of bacterial communities along the stratified water column in the Black Sea.

The Black Sea is the largest semi-closed permanently anoxic basin on our planet with long-term stratification. The study aimed at describing the Black Sea microbial community taxonomic and functional composition within the range of depths spanning across oxic/anoxic interface, and to uncover the factors behind both their vertical and regional differentiation. 16S rRNA gene MiSeq sequencing was applied to get the data on microbial community taxonomy, and the PICRUSt pipeline was used to infer their functional profile. The normoxic zone was mainly inhabited by primary producers and heterotrophic prokaryotes (e.g., Flavobacteriaceae, Rhodobacteraceae, Synechococcaceae) whereas the euxinic zone-by heterotrophic and chemoautotrophic taxa (e.g., MSBL2, Piscirickettsiaceae, and Desulfarculaceae). Assimilatory sulfate reduction and oxygenic photosynthesis were prevailing within the normoxic zone, while the role of nitrification, dissimilatory sulfate reduction, and anoxygenic photosynthesis increased in the oxygen-depleted water column part. Regional differentiation of microbial communities between the Ukrainian shelf and offshore zone was detected as well, yet it was significantly less pronounced than the vertical one. It is suggested that regional differentiation within a well-oxygenated zone is driven by the difference in phytoplankton communities providing various substrates for the prokaryotes, whereas redox stratification is the main driving force behind microbial community vertical structure.

INTRODUCTION

The Black Sea is the largest semi‐closed permanently anoxic basin on our planet (Vetriani et al., 2003). The massive vertical redox stratification has been maintained for over 7000 years and has developed as a result of a synergetic effect of active photosynthetic primary production in the surface waters, the intense organic carbon flux, and the characteristic bottom configuration (Fuchsman et al., 2012; Jones & Gagnon, 1994). Such long‐term stratification is a peculiarity of the Black Sea and is the main difference from the other anoxic water bodies like the Baltic Sea (Hannig et al., 2007) or Saanich Inlet (Manning et al., 2010), where periodic oxygenation occurs. The Black Sea contains at least 3 distinct zones: normoxic, suboxic, and euxinic. Normoxic zone includes upper well‐oxygenated layer (O2 = 8.89 ± 1.19 mg/l). The suboxic zone, where the concentration of dissolved oxygen is extremely low (less than 2 mg/l), is about 50 m thick and occupies the boundary between normoxic and euxinic zones (Fuchsman, Staley, et al., 2012). The suboxic zone is followed by the euxinic zone that is completely devoid of oxygen and contains hydrogen sulfide. According to the most recent data collected in the course of the EMBLAS‐II project (Slobodnik et al., 2017), the depth of the top of the euxinic zone varies between stations but generally begins between 54 and 130 m. Vertical stratification of the water column makes the Black Sea a perfect model to study prokaryotic niche partitioning along with the oxic/anoxic interfaces.

Prokaryotes from the Black Sea water column were intensively studied with both culture‐based and culture‐independent techniques due to their unique character, diversity, and functional capacity. Karl and Knauer (1991) evaluated the microbially mediated organic carbon cycle in the Western Basin of the Black Sea. Aytan et al. (2018) studied biomass composition and trophic interactions of the plankton communities of the Southeastern Black Sea coastal area. The dynamics of bacterioplankton community composition during phytoplankton blooms in the shelf surface waters was discussed by Stoica and Herndl in 2007. Analysis of the photic zone archaeal community structure within the continental shelf in the Gelendzhik area was performed by Merkel et al. (2015).

Microbial communities inhabiting the redoxcline and euxinic zone of the Black Sea attracted substantial attention, as they can drive several alternative metabolic pathways involved in carbon (Fuchsman et al., 2011; Michaelis et al., 2002), nitrogen (Fuchsman et al., 2011; Fuchsman, Staley, et al., 2012; Kirkpatrick et al., 2006; Kuypers et al., 2003; Lin et al., 2006; Oakley et al., 2007; Suominen et al., 2019), and sulfur (Fuchsman et al., 2011; Grote et al., 2007) biogeochemical cycles.

Nitrogen cycling in oxygen‐deficient Black Sea zones has been a major research focus during the past 10–15 years because of the globally important transformation processes (e.g., nitrification, anaerobic ammonium oxidation) occurring in these systems. One of the first reports of anammox in an oxygen‐deprived water column was from the Black Sea (Kuypers et al., 2003), which demonstrated the abundance and importance of anammox bacteria in the nitrogen cycle. Anammox peaks in the Black Sea suboxic zone, where NO3 − and NH4 + profiles intersect (Fuchsman, Staley, et al., 2012; Lam et al., 2007). The sinking organic matter (eg. from the phytoplankton blooms) and benthic ammonium release fuels the upward flux of ammonium from the euxinic layer of the Black Sea (Fuchsman et al., 2008; Fuchsman, Staley, et al., 2012). The anammox is coupled with ammonium oxidation by Thaumarchaeota ammonia‐oxidizers (within the lower oxic zone) and γAOB (within the suboxic zone), which support this process with NOx Lam et al. (2007).

The role of Epsilonproteobacteria and Gammaproteobacteria in primary production and sulfur cycling in the anoxic zone of the Black Sea was demonstrated by Grote et al. (2008) and Glaubitz et al. (2010). Namely, bacteria related to Sulfurimonas and gammaproteobacteria sulfur oxidizer cluster (GSO) were shown to be responsible for sulfide oxidation and CO2 fixation in the redoxcline of the Black Sea (Grote et al., 2008 and Glaubitz et al., 2010). Notably, chemoautotrophic production has been shown to comprise up 43%–89% of the total water column production, reaching its maximum in the Black Sea upper euxinic zone (Ediger et al., 2019; Yilmaz et al., 2006). High chemosynthesis levels observed at the Black Sea suboxic/euxinic interface result in the second maximum of S‐POC and S‐PON (Çoban‐Yildiz et al., 2006; Fuchsman et al., 2019; Kirkpatrick et al., 2019).

The first attempt to systematically characterize Black Sea microbial communities with various metabolic pathways inhabiting oxic/ anoxic interface was performed by Vetriani et al. (2003). The major groups identified were Gammaproteobacteria (with a dominance of Pseudoalteromonas like clones), and Epsilonproteobacteria. Members of Deltaproteobacteria related to sulfate reducers and the Archaea related to phylotypes from the ANME groups that anaerobically oxidize methane were found as well. Yet, the majority of clones were not restricted to a specific depth in the water column, and many of the major T‐RFLP peaks remained uncharacterized (Vetriani et al., 2003). A more fine‐scale analysis was conducted by Fuchsman et al., 2012, and the role of SUP05, Sulfurimonas, Arcobacter, and BS‐GS02 in autotrophic denitrification was uncovered at the Black Sea Western Gyre and mixing zone.

The Black Sea attracted scientific attention due to its pronounced stratification and the presence of a well‐defined oxic‐anoxic interface, which allowed for the detailed studies of particular microbial processes. Taking into account significant variation in taxonomic composition and abundance of Black Sea prokaryotes across the water column, it is highly important to get a comprehensive picture of microbial community taxonomic distribution in relation to varying physical and chemical parameters. Previous studies have focused on the microbial community at one or two stations. Here, we examine the microbial community across 12 stations including the central Black Sea, the northwestern shelf, and the southeastern shelf. An attempt has been made to explain microbial niche partitioning in highly stratified water bodies and to address the factors shaping microbial communities in such ecosystems, both in well‐oxygenated and in oxygen‐depleted waters.

MATERIALS AND METHODS

Bacterioplankton collection

Sixty eight seawater samples from 12 stations were collected during the Joint Black Sea Survey, which was conducted under the EU/UNDP EMBLAS‐II project (Improving Environmental Monitoring in the Black Sea: http://emblasproject.org/) in May‐June 2016. The cruise route went through the central part of the Black Sea (Stations from #2 to #10), as well as through the shelf of the eastern (Stations #11 and #12) and northwestern Black Sea (Station #1). The location of the sampling points is presented in Figure 1 and Appendix Table A1.

FIGURE 1

Sampling stations located across the Black Sea

TABLE A1

Geographical location and physico‐chemical parameters of the samples

Station	Sample	Long.	Lat.	Depth, m	Sigma‐theta	Dissolved O2, mg/l	Fluorescence Seapoint	H2S, mg/l	Temp [°C]	NO3, µg/l	NO2, µg/l	Norg, µg/l	Porg, µg/l	PO4, µ/l	Si, µg/l	Ntotal, µg/l	Ptotal, µg/l
1	1.1	31.331	44.831	0	6.9	8.9	0.3	—	21.4	1.9	<LoQ	306	8.2	<LoQ	410	308	8.2
1	1.2	31.331	44.831	5	9.6	8.8	0.7	—	21.2	1.5	<LoQ	266	7.5	<LoQ	400	268	7.5
1	1.3	31.331	44.831	17	13.6	9.6	1.9	—	9.7	11.9	2.3	253.8	5.6	<LoQ	520	268	5.6
1	1.2.1	31.331	44.831	25	14.0	9.3	1.1	—	8.7	27.4	0.6	219.9	<LoQ	<LoQ	630	248	<LoQ
1	1.4	31.331	44.831	60	14.4	6.1	0.6	—	8.7	34.9	<LoQ	99.1	4.4	12.6	460	134	17
2	2.1	31.3883	44.6361	0	11.4	8.7	0.5	—	21.2	3.9	<LoQ	254	6	<LoQ	27.9	258	6
2	2.2	31.3883	44.6361	12	12.6	10	0.6	—	16	0.5	<LoQ	282	<LoQ	<LoQ	35.9	283	<LoQ
2	2.3	31.3883	44.6361	40	13.9	9.7	1.1	—	9.9	2.5	0.6	244.9	<LoQ	<LoQ	51.8	248	<LoQ
2	2.4	31.3883	44.6361	123	15.4	0.7	0.1	—	8.7	52.8	<LoQ	453.2	3.1	40.3	1020	506	43.4
3	3.1	31.5675	44.1582	0	11.5	9.2	0.3	—	19.9	1	<LoQ	451	8.8	<LoQ	39.8	452	8.8
3	3.2	31.5675	44.1582	12	12.7	9.5	0.5	—	16.5	1.5	<LoQ	266	<LoQ	<LoQ	54.5	268	<LoQ
3	3.3	31.5675	44.1582	40	14.0	9.1	0.9	—	9.3	49.8	<LoQ	228	8	<LoQ	130	278	8
3	3.4	31.5675	44.1582	121	15.5	0.7	0.1	0.9	8.6	60.3	<LoQ	212.7	8.1	41.8	1070	273	49.9
3	3.5.1	31.5675	44.1582	135	15.8	—	0.1	1.1	8.6	—	—	—	—	—	—	—	—
3	3.5.2	31.5675	44.1582	144	16.0	—	0.1	1.2	8.7	—	–	—	—	—	—	—	—
4	4.1	31.7105	43.7538	0	12.8	8.9	0.3	—	18	2.5	<LoQ	146.5	<LoQ	<LoQ	80.9	149	<LoQ
4	4.2	31.7105	43.7538	8.2	12.8	9.6	0.3	—	17	3	<LoQ	389	<LoQ	<LoQ	82.2	392	<LoQ
4	4.3	31.7105	43.7538	44.5	14.3	9.4	1	—	8.8	7	0.7	592.4	3.3	6.3	160	600	9.6
4	4.4	31.7105	43.7538	59.8	15.2	0.9	0.4	—	8.6	149	0.8	560.2	5.9	40.6	830	710	46.6
4	4.5	31.7105	43.7538	99.3	16.0	0.2	0.1	—	8.7	4	0.6	1251.4	1.2	1.2	1470	1256	109.2
5	5.1	31.8332	43.3675	0	12.6	9.1	0.2	—	18.1	0.5	<LoQ	718.5	7.2	<LoQ	67.6	719	7.2
5	5.2	31.8332	43.3675	15.3	13.2	9.3	0.3	—	16	1	<LoQ	669	12.8	<LoQ	69	670	12.8
5	5.3	31.8332	43.3675	45.8	14.7	6.9	0.9	—	8.5	99.6	0.7	484.6	1.3	15.6	370	585	16.7
5	5.4	31.8332	43.3675	96.3	16.0	1.1	0.1	1.5	8.7	64.7	<LoQ	412.3	1.2	118.6	1430	477	119.8
5	5.5	31.8332	43.3675	110.9	16.1	—	0.1	0.8	8.7	—	—	—	—	—	—	—	—
6	6.1	32.5731	43.3657	0	12.7	9.4	0.4	—	17.6	2	<LoQ	464	9.8	<LoQ	87.5	466	9.5
6	6.2	32.5731	43.3657	12.5	13.3	9.4	0.5	—	14.8	2	<LoQ	454	0	<LoQ	91.5	456	<LoQ
6	6.3	32.5731	43.3657	55	14.5	6.9	1.1	—	8.4	74.7	0.5	609.8	3.5	16.5	370	685	20
6	6.4	32.5731	43.3657	98.8	15.9	0.7	—	—	8.6	1	0.5	627	13.2	148	1490	628	161.2
7	7.1	32.861	43.4174	0	12.1	8.6	0.3	—	20	1	<LoQ	257	<LoQ	<LoQ	79.5	258	<LoQ
7	7.2	32.861	43.4174	8	12.2	8.6	0.3	—	19.8	1	<LoQ	287	<LoQ	<LoQ	75.6	288	<LoQ
7	7.3	32.861	43.4174	50	14.5	5	1.2	—	8.4	37.4	0.5	250	2.1	21.9	490	288	24
7	7.4	32.861	43.4174	58	15.0	2	0.3	1.3	8.6	66.2	0.6	181.1	4.2	31	750	248	35.2
7	7.5.1	32.861	43.4174	86.5	15.8	—	0.1	1.3	8.6	—	—	—	—	—	—	—	—
7	19.5.2	32.861	43.4174	100	16.0	—	0.1	1.2	8.6	—	—	—	—	—	—	—	—
8	8.1	34.4132	43.3446	0	12.1	8.8	0.4	—	19.9	5.5	<LoQ	222.5	<LoQ	<LoQ	82.2	228	<LoQ
8	8.2	34.4132	43.3446	15	12.7	9.9	0.7	—	13.3	4.9	<LoQ	253	6	<LoQ	84.8	258	6
8	8.3	34.4132	43.3446	48	14.4	8.6	0.8	—	8.6	61.8	1.6	209.6	4.6	12	210	273	16.6
8	8.4.1	34.4132	43.3446	62	15.3	0.6	0.3	1.2	8.6	47.3	0.7	154.9	10.5	41.8	820	203	52.3
8	8.4.2	34.4132	43.3446	82	15.8	0.2	0.1	1.4	8.6	30.4	2.6	239.9	7.7	<LoQ	1210	273	7.7
8	8.5	34.4132	43.3446	106	16.2	—	0.1	1.4	8.6	—	—	—	—	—	—	—	—
9	9.1	36.0697	43.526	0	12.3	8.6	0.4	—	19.4	2	<LoQ	211	6.5	<LoQ	71.6	213	6.5
9	9.2	36.0697	43.526	14	13.1	9	0.5	—	16.2	2	<LoQ	266	<LoQ	<LoQ	70.3	268	<LoQ
9	9.3	36.0697	43.526	45	14.7	4.2	0.8	—	8.4	18.9	<LoQ	378.1	8.5	26.2	540	397	34.7
9	9.4.1	36.0697	43.526	54	15.4	0.4	0.2	0.8	8.6	38.8	0.6	183.6	1.2	42.4	890	223	43.6
9	9.4.2	36.0697	43.526	70	15.8	0.3	0.2	0.9	8.7	22.9	1.2	149.9	0.6	15.3	1230	174	15.9
9	9.5	36.0697	43.526	83	16.0	—	0.1	1.2	8.6	—	—	—	—	—	—	—	—
10	10.1	39.886	42.2345	0	12.2	9.2	0.4		18.2	2	<LoQ	325	5.9	<LoQ	100	327	5.9
10	10.2	39.886	42.2345	12	12.3	9.4	0.7	—	17.6	3.9	<LoQ	527	8	<LoQ	91.5	531	8.1
10	10.3.1	39.886	42.2345	30.5	14.0	9.6	1	—	9.4	89.6	0.9	1298.5	6.5	<LoQ	110	1389	6.5
10	10.3.2	39.886	42.2345	75	14.6	5.8	0.1	—	8.6	79.6	<LoQ	228.4	6	18	410	308	24
10	10.4.1	39.886	42.2345	112.4	15.4	1.9	0.1	—	8.7	14.7	<LoQ	446.3	2.7	39.7	860	521	42.5
10	10.4.2	39.886	42.2345	129.7	15.8	0.6	0.1	0.2	8.7	59.8	<LoQ	287.2	3.6	49.6	1220	347	53.2
10	10.5.1	39.886	42.2345	159	16.2	—	0	2.2	8.7	2540	<LoQ	3061	5.6	188	1420	5601	193.6
10	10.5.2	39.886	42.2345	748	17.1	—	0	9.1	9	2490	<LoQ	1628	60.8	20.4	7280	4118	81.2
10	10.5.3	39.886	42.2345	996	17.1	—	0	11.1	9	3	<LoQ	240	38.6	39.7	7500	243	78.3
11	11.1	40.5787	42.0172	0	12.1	8.9	0.6	—	18.7	2	<LoQ	570	13.1	<LoQ	67.6	572	13.1
11	11.2	40.5787	42.0172	8.9	12.2	9.3	0.6	—	18.4	3.9	0.7	565.3	9.7	<LoQ	71.6	570	9.7
11	11.3	40.5787	42.0172	35.4	14.0	9.6	1.2	—	9.3	59.8	1.3	390.9	13.8	<LoQ	70.3	452	13.8
11	11.4.1	40.5787	42.0172	134	15.4	1.3	0.1	—	8.6	110	0.6	467.4	0.5	42.1	810	578	42.6
11	11.4.2	40.5787	42.0172	165	16.0	0.3	0.1	—	8.7	139	0.7	136.3	2.2	90	1140	276	92.2
12	12.1	41.2136	41.7841	0	7.2	8.9	0.2	—	18.6	2	<LoQ	425	<LoQ	<LoQ	130	427	<LoQ
12	12.2	41.2136	41.7841	7	12.2	9.5	0.3	—	18.6	1.5	<LoQ	315.5	<LoQ	<LoQ	94.1	317	<LoQ
12	12.3	41.2136	41.7841	35	14.1	9.6	1.5	—	9	39.8	0.9	128.3	<LoQ	<LoQ	110	169	<LoQ
12	12.4	41.2136	41.7841	127	16.0	0.3	0.1	—	8.7	32.4	0.8	373.8	5.9	59.8	1430	407	65.7

Abbreviation: LoQ, limit of quantitation.

Samples were taken from the following layers: surface, thermocline, deep chlorophyll‐a maximum (DCM), suboxic, and euxinic zones. Sampling depths were selected according to the profile of the Conductivity, Temperature, and Pressure Recorder (CTD‐SBE 25 plus, Sea‐Bird Scientific, US). The downcast data on temperature and fluorescence were binned to 1 m depth intervals with SBE Data Processing Software version 7.18 (Sea‐Bird Electronics, Bellevue, Washington, USA). Dissolved oxygen measurements were done onboard using the Winkler method (Hansen, 1999), whereas hydrogen sulfide concentration was determined with back‐titration of unused iodine (Ciesielski & Zakrzewski, 2006).

Five liters of seawater were taken from each of the sampling layers in a Niskin bottle inserted in CTD rosette. 2 and 5 L samples were passed through the Millipore Sterivex‐GP 0.22 µm filters (Millipore, USA) using 50 ml sterile syringes. The filtration volume depended on the turbidity of water. The filters were immediately frozen in liquid nitrogen and subsequently stored at −80°C until DNA extraction.

Phytoplankton collection

During the cruise, samples of phytoplankton were collected from the same stations by vertical series from layers of surface, thermocline, and deep chlorophyll‐a maximum. 3 L of water was collected from each sample with 5 L Niskin bottles, attached to the CTD rosette system. Samples were concentrated on board using the funnel of inverted filtration to the volume of 50–100 ml and then fixed with 4% buffered formaldehyde up to the final concentration of 2% in a sample. Samples were concentrated one more time down to 10–20 ml by slow decantation before the visual assessment. Species were identified and quantified under a light microscope LOMO at 600× magnification in the volume of 0.05 ml. The wet biomass (mg/m3) was calculated by the method of geometric similarity, which identifies cells` shapes as corresponding geometrical Figs (Utermohl, 1958).

Species identification was mainly using Schiller (1937), Kisselew (1950), Proshkina‐Lavrenko (1955), Carmelo (2007), Steidinger and Tangen (1997), Cronberg and Annadotter (2006), and the taxonomic nomenclature according to the online database of World Register of Marine Species (WoRMS) and the Black Sea check‐list http://phyto.bss.ibss.org.ua.

Dominant species were identified according to the Brotskaya–Zenkevich index (Zenkevich & Brotskaya, 1937).

DNA extraction and sequencing

Prokaryotic DNA was extracted from 64 water samples using MO BIO PowerSoil® DNA Isolation Kit (MO BIO Laboratories, Inc., USA). Each Sterivex filter was removed from the casing by cracking the housing with sterile pliers and placed in 2 ml sterile tubes and processed according to the manufacturer's protocol. The DNA quantity and quality were double‐checked using NanoDrop Spectrophotometer (Thermo Fisher Scientific, USA) and Qubit 2.0 Fluorometer (Thermo Fisher Scientific, USA). All samples had sufficient DNA concentration (ranging from 4.74 to 226 ng/μl with most above 30 ng/μl) and purity ratio ~1.8 at A260/280 nm.

Sequencing was performed at MR DNA (Shallowater, TX, USA) using the Illumina MiSeq sequencing platform following the manufacturer's guidelines. A previously described MiSeq 16S rRNA gene protocol based on the bTEFAP process (Chiodini et al., 2015; Dowd et al., 2008) with 16S rRNA gene V4 region bacterial primer pair S‐D‐Bact‐0341‐b‐S‐17 and S‐D‐Bact‐0785‐a‐A‐21 (Klindworth et al., 2013) with a barcode on the forward primer was used to evaluate bacterial diversity of samples. In brief, a single‐step 28 cycle PCR using the HotStarTaq Plus Master Mix Kit (Qiagen, USA) was used under the following conditions: 94°C for 3 min, followed by 28 cycles of 94°C for 30 s, 53°C for 40 s, and 72°C for 1 minute, after which a final elongation step at 72°C for 5 min was performed. Subsequently, PCR products were checked in 2% agarose gel to determine the success of amplification and the relative intensity of bands. Multiple samples were pooled together in equal proportions based on their molecular weight and DNA concentrations. The purification of pooled samples was carried out using calibrated Ampure XP beads. Then, the pooled and the purified PCR product was used to prepare the Illumina DNA library.

Bioinformatic and statistical analyses

The reads were merged and reoriented in the 5’‐3’ direction and <200 bp sequences removed by MR DNA (www.mrdnalab.com, Shallowater, TX, USA). Sequences were depleted of barcodes and sequences of <150 bp were removed in QIIME1 (Caporaso et al., 2010). Sequences were denoised, Operational taxonomic units (OTUs) generated, and chimeras removed. OTUs were defined by de novo clustering at 97% similarity. Taxonomic annotation of OTUs was performed with Greengenes 18.5 database. Sequences were uploaded in the NCBI database under accession number PRJNA576012.

Microbial communities functional profile predictions were performed with PICRUSt2 pipeline (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) based on 16S rRNA gene sequencing data (Douglas et al., 2019). PICRUSt2 allows for microbial community functional potential prediction, yet it has certain limitations that should be recognized when interpreting the results. Firstly, it is limited by the availability of sequenced reference genomes, and secondly, it cannot distinguish the strain‐specific functionality being an amplicon‐based analysis (Douglas et al., 2019).

Shannon diversity and Bray–Curtis distance matrix were calculated in QIIME2 2019.1 (https://qiime2.org/, Caporaso et al., 2010) using the “diversity core metrics” plugin. ANOSIM analysis with Bray–Curtis and 9999 permutations was performed in the vegan R package to test both for the vertical, and the regional statistical difference between the microbial communities and NMDS plots was created for visualization. The difference in Shannon index and specific OTUs (pathways) abundance both between layers (surface, thermocline, DCM, suboxic and euxinic layers) and between regional groups (Ukrainian shelf and Offshore) was tested with the Kruskal–Wallis test (p‐value < 0.05). Mantel test was performed in R vegan package with the purpose to infer environmental parameters explaining the variance in microbial community taxonomic and functional profile. The following matrices were used in the Mantel test: the abundance matrix (based on Bray–Curtis measure), the environmental parameter distance matrix (using Euclidean distance), and the geographical distance matrix (based on Haversine distance). All regional comparison analyses were performed based on the data from the surface, thermocline, and DCM layers, as only those were present in the Ukrainian shelf within its shallow depth.

RESULTS AND DISCUSSION

Environmental parameters of the water layers of the Black Sea

Microbial communities of the Black Sea water column were analyzed from five layers: (i) surface; (ii) thermocline; (iii) deep chlorophyll‐a maximum; (iv) suboxic zone; (v) euxinic zone. The sampling depth of the surface layer was 0–1 m. The thermocline was located within the depth range 5–15 m immediately before the temperature dropped down from 17 to 18°C to 8°C (Figure 2). The deep chlorophyll‐a maximum layer was characterized by the increase in fluorescence, which indicates the photosynthetic activity of primary producers. The maximum fluorescence varied from 0.9 to 1.9 at depths 17–55 m. The suboxic zone had an oxygen content of less than 2 mg/l. It occupied the depths of 112–134 m at the near‐shore stations (St # 1, 2, 3, 10, 12) and 54–99 m at the offshore stations (St. # 4–9) of the Black Sea, corresponding to σθ = 15.01–16.02. The position of the euxinic layer located below the transition suboxic zone varied between the stations with the minimum depth of 83 m at Station #9 and the maximum depth of 140 m at Station #10, which is in accordance with the previous data (Jørgensen et al., 1991).

FIGURE 2

Vertical distribution of physico‐chemical parameters across the sampling grid: (a) temperature, (b) fluorescence at Seapoint, (c) dissolved O2 concentration, (d) H2S concentration

Sequencing output and diversity of the community

Illumina MiSeq high‐throughput sequencing was used to analyze the diversity and composition of prokaryotic communities inhabiting the Black Sea water column. Sequencing output counted from 44,349 to 165,791 sequences (Table A2) of partial 16S rRNA gene per sample with the average length 431.5 bp [175.0–545.0]. Based on rarefaction curves outlined from samples (Table A1), it was concluded that deeper sequencing would not have resulted in significantly higher estimates. Therefore, the study covered both abundant and rare species of the communities.

The number of OTUs was in the range of 203–368 (Table A2). Shannon index indicates the high richness of the communities inhabiting all zones analyzed within the Black Sea water column and even distribution of OTUs within the communities (Figure 3b). Unlike the thermocline and surface layers, the diversity of the DCM, suboxic, and euxinic layers had a broader range of variation. Similarly, the diversity of microbial communities inhabiting the Ukrainian shelf had a broader range of variation than at the offshore zone, yet there was no significant difference in alpha diversity between the two regions (Figure 3a).

FIGURE 3

Shannon diversity indexes of prokaryotic communities populating water column of the Black Sea (a) at different Black Sea regions, (b) at different depths

ANOSIM analysis based on the Bray–Curtis measure revealed both the regional (R = 0.89, p = 0.0001) and the vertical difference (R = 0.64, p = 0.0001) in Black Sea microbial community taxonomic structure. Yet, in terms of the regional structure only the samples from the Ukrainian shelf clustered separately (Figure 4a); therefore, the subsequent analysis was performed with Georgian and offshore samples pooled together. Similarly, the microbial community taxonomic structure was not significantly different in surface and thermocline layers; thus, they are discussed together hereinafter.

FIGURE 4

Bray–Curtis regional and vertical differentiation of Black sea microbial communities based on (a) taxonomic data, (b) functional data (KEGG pathways)

Notably, the microbial functional differentiation profile was different from the taxonomic one (Figure 4b). Namely, only vertical differentiation was significant (R = 0.2655, p = 0.0001). Similarly, to the taxonomic structure, the surface and thermocline clustered together and were not significantly different. Thus, their functional profile is discussed together.

Taxonomic composition of microbial communities populating the Black Sea water column

The taxonomic and functional composition of microbial communities changed in accordance with environmental parameters. Though the data on DNA sequences do not indicate the activity of taxa, the relative abundance of the taxa and functional pathways in the community suggested a clear pattern in accordance with the distribution of chemical parameters (Figure 2, Table A1) at the different layers. Below the abundance of the microbial taxa on the class (Figure 5) and genus level (Figure 6) is discussed within the normoxic, suboxic, and euxinic conditio

FIGURE 5

The Black Sea prokaryotic diversity at the class leve

FIGURE 6

The Black Sea prokaryotic diversity at the family level

Normoxic water column

Normoxic water column includes well‐oxygenated (8.89 ± 1.19 mg/l O2) surface, thermocline, and deep chlorophyll‐a maximum layers. The majority of prokaryotes inhabiting normoxic waters (0–55 m) are represented by primary producers and heterotrophic prokaryotes that drive a microbial loop.

Surface and thermocline layers

Heterotrophic bacteria comprised the major part of the communities inhabiting surface and thermocline layers of the Black Sea, degrading the organic matter released by phytoplankton. The high diversity of heterotrophs can be the result of bacterial specialization to carbon substrates. Emiliania huxleyi can be especially involved in the niche partitioning of heterotrophic prokaryotes. The bacterioplankton was sampled at the end of May, when coccolithophores begin to dominate in the Black Sea (Cokacar et al., 2001; Mikaelyan et al., 2011; Oguz & Merico, 2006; Stoica & Herndl, 2007). Emiliania huxleyi (Brotskaya‐Zenkevich Index (BZI) = 116.1) (Zenkevich & Brotskaya, 1937) amounted up to 30% of the phytoplankton community in the surface layer and had a basin‐wide character of distribution.

Diatoms were the most mass‐abundant species at the moment of sampling. According to Brotskaya–Zenkevich index, the most abundant species were diatoms Dactyliosolen fragilissimus (BZI = 1111.5), Pseudo‐nitzschia delicatissima (BZI = 521), Pseudosolenia calcar‐avis (BZI = 469.7), Proboscia alata (BZI = 300.2) followed by dinoflagellates Diplopsalis lenticula (BZI = 521.0), and Tripos furca (BZI = 135.3). Labile algal exudates (amino acids, polysaccharides, lipids, etc) and particulate organic matter are consumed by distinct and diverse taxa, so the number of prokaryotic taxa may reflect the niche partitioning and optimization of resource assimilation (Landa et al., 2016; Williams et al., 2013).

Alphaproteobacteria (32%–37%), Gammaproteobacteria (14%), and Flavobacteriia (14%), which usually outnumber the other taxa in the marine environment (Williams et al., 2013), dominated surface and thermocline microbial communities (Figure 5). Synechococcophycideae comprised about 10% and 14% at the surface layer and the thermocline, respectively (Figure 5).

A more detailed analysis was conducted at the family level. Flavobacteriaceae, Rhodobacteraceae, and Synechococcaceae were the most abundant families within the well‐oxygenated surface and thermocline layers and comprised 11%, 11%, and 10%–13% respectively (Figure 6). Synechococcaceae family is ubiquitous and widespread in temperate waters and is considered to be the most important component of photosynthetic picoplankton (Scanlan & West, 2002). The high abundance of Flavobacteriaceae and Rhodobacteraceae coincides with the availability of growth substrates provided by blooming phytoplankton. The members of the marine Flavobacteraceae family are aerobic chemoheterotrophic bacteria known to proliferate during both early and mid‐phytoplankton bloom stages, colonizing algae cell surfaces (Eckert et al., 2012; Krüger et al., 2019; Zhang et al., 2019). Flavobacteraceae opportunistic nature and high diversity of polymer‐degrading enzymes contribute to their success in occupying econiches provided by phytoplankton particles, and they are known to account for up to 70% of bacterioplankton populations associated with algal blooms (Bowman, 2006).

Marine Rhodobacteraceae are major vitamin suppliers for eukaryotic primary producers (e.g., diatoms, dinoflagellates, and brown algae), and they have been previously shown to constitute the major part of bacterioplankton community during phytoplankton bloom (Haggerty & Dinsdale, 2017; Zhang et al., 2019). The members of the Roseobacter clade belonging to the Rhodobacteraceae family are associated with coccolithophore bloom in the Northwestern Black Sea and to be strongly dependent on DMSP‐producing phytoplankton species (Bakenhus et al., 2017; Stoica & Herndl, 2007).

Halomonadaceae that belong to Gammaproteobacteria accounted for 6%–7% within the surface and thermocline. The microbes belonging to this family are known as aerobic or facultatively anaerobic, chemoorganotrophs (de la Haba et al., 2014). They are associated with primary producers bloom and have been shown to play an important role in the nitrogen cycle being capable of nitrification, denitrification, or simultaneous nitrification and denitrification (Wilson et al., 2017).

Deep chlorophyll‐a maximum (DCM)

The microbial communities of the DCM layer were characterized by reduced quantities of Alphaproteobacteria (22%) compared to surface and thermocline. Gammaproteobacteria (18%) were the second most abundant class followed by Flavobacteriia (15%) and Synechococcophycideae (13%) (Figure 5). The abundance of Acidimicrobiia increased significantly and constituted 11%.

Bacteria potentially involved in the cleavage of algal polymers occupied a significant position in the community. Indeed, Flavobacteriacea were present in high abundance (8%), as well as Rhodobacteraceae (8%) and Halomonadaceae (7%). Additionally, the abundance of OCS155 (6%) and ZA3409c (2%) belonging to Actinobacteria, NS9 (3%), Plactomycetaceae (3%), Verrucomicrobiaceae (3%), increased significantly compared to the upper layers.

Representatives of the OCS155 family are photoheterotrophic free‐living bacteria (Angly et al., 2016) that have been previously shown to be abundant in deep chlorophyll‐a maximum (West et al., 2016) and to be related to high phytoplankton concentrations (Mizuno et al., 2015; Morris et al., 2012), in particular to the diatom‐enriched environment (Nelson et al., 2014).

NS 9 family belongs to Flavobacteria, which are known for their association with phytoplankton blooms (Liu et al., 2019) due to their ability to degrade high molecular weight organic matter and to colonize substrates (Buchan et al., 2014). Likewise, heterotrophic carbohydrate‐degrading Verrucomicrobiaceae contribute to phytoplankton‐derived organic matter degradation at DCM (Cardman et al., 2014). Planctomycetaceae might be involved in the initial aerobic breakdown of complex organic matter into simpler compounds at DCM (Glöckner et al., 2003). Both free‐living and particle‐attached Planctomycetaceae have been previously identified in microbial communities inhabiting both oxic and suboxic Black Sea zones (Fuchsman, Staley, et al., 2012; Kirkpatrick et al., 2006).

Suboxic zone

The oxygen‐depleted water layer was populated by heterotrophic and chemoautotrophic taxa that benefit both from the remaining light and reduced sulfur compounds. Alphaproteobacteria (12%) and Flavobacteria (8%) were characterized by much lower abundance than within DCM and were outnumbered by Gammaproteobacteria (20%), AB16 (15%), and Acidimicrobia (13%) (Figure 5). Deltaproteobacteria (8%), Chlorobia (7%), and Phycisphaerae (5%) were also constituted a significant portion of the suboxic microbial community. The number of unassigned bacteria has increased in the non‐oxygenated part of the Black Sea water column compared to upper layers and constituted 5% in the suboxic zone (Figure 5).

Synechococcophycideae are not able to survive within the suboxic zone of the Black Sea, as they require more intense irradiation for growth (Carey et al., 2012). Hence, the niche is occupied by green sulfur Chlorobia (7%) and green non‐sulfur bacteria Anaerolineae (1%), which grow in dim euxinic environments (Olson, 1998). Earlier, Manske et al. (2005) have shown phototrophic green sulfur bacteria to persist in the Black Sea redoxcline at >100‐m depth at light intensities of <0.001% of surface radiation. Marschall et al. (2010) have shown that Chlorobium BS‐1 was metabolically active at the offshore station located at the center of the Black Sea at the maximum depth of 150 m in contrast to the bacteria from the periphery. Green sulfur bacteria were mainly represented by the Chlorobiaceae (6%) family in our dataset (Figure 6).

Gammaproteobacteria were represented by Piscirickettsiaceae (6%) and Thiohalorhabdaceae (1%), which increased in their abundance compared to the well‐oxygenated water column part and Halomonadaceae (2%) that decreased (Figure 6). Piscirickettsiaceae are known to be adapted to the oxic/anoxic interface and capable of sulfide oxidation (Borin et al., 2009; Kappler et al., 2016). Similarly, Thiohalorhabdaceae (1%) use reduced sulfur compounds as an energy source and O2 or NO3 as electron acceptors (Sorokin et al., 2020).

AB16 group belongs to Marinimicrobia, which had been frequently detected in oxygen minimum zones and euxinic zones worldwide (Allers et al., 2013; Cui et al., 2019; Thrash et al., 2017), including the Black Sea (Fuchsman et al., 2011; Suominen et al., 2019). Even though Marinimicrobia has been mainly found in oxygen minimum zones (OMZs), it has been shown to have the genomic potential for aerobic respiration (cytochrome c oxidases). At the same time, the genes regulating anaerobic metabolism, for example, dissimilatory nitrate reduction, nitrous oxide reduction, and sulfur reduction are co‐expressed (Suominen et al., 2019). This suggests the existence of a putative potential for switching between anaerobic and aerobic metabolism, which might be beneficial in the ecotone environment of OMZs. A714017 family, which belongs to Ab16 class, was the most abundant in the suboxic zone accounting for 12% (Figure 6).

Acidimicrobia are known to be a metabolically versatile group that can utilize various sources of energy (Mizuno et al., 2015). They have been frequently found to represent a significant portion of microbial communities inhabiting DCM and deep photic zone utilizing photoheterotrophic strategies to benefit from sinking organic matter and remaining light (Ghai et al., 2014; Mizuno et al., 2015). ZA3409c family of Acidimicrobia class is the second most abundant on our dataset with 9% (Figure 6).

Deltaproteobacteria that increased to 8% in the suboxic zone of the Black Sea water column were represented by SAR 324 family (5%). SAR 324 includes bacteria that play an important role in the sulfur and nitrogen cycle in oxygen minimum and anoxic environments being capable of nitrate reduction and dissimilatory sulfate reduction (Cao et al., 2016; Sheik et al., 2014). This group has been previously detected across the Black Sea suboxic zone (Fuchsman et al., 2011).

Alphaproteobacteria that decreased to 12% compared to the surface, thermocline, and DCM were represented by the ubiquitous Pelagibacteraceae (4%) family, Rhodospirillaceae (3%) family, and some minor taxa. Rhodospirillaceae are purple non‐sulfur bacteria capable of growing both photoheterotrophically under anoxic conditions in the presence of light and chemoheterotrophically in the dark (Baldani et al., 2014) playing an important role in the marine nitrogen cycle due to their universal capacity to fix molecular N2 (Madigan et al., 1984).

Epsilonproteobacteria increased considerably in their abundance compared to the upper well‐oxygenated layers. Epsilonproteobacteria were shown to be the major dark CO2‐fixing organisms in the Black Sea constituting 75%–100% of CO2 assimilation (Grote et al., 2007). This class accounted for 2% of the total microbial community and was represented by one family—Helicobacteraceae. The role of both Epsilon‐ and Gammaproteobacteria in dark CO2 fixation within both the Black Sea and Baltic Sea redoxclines has been demonstrated earlier (Glaubitz et al., 2010; Grote et al., 2007). It was also suggested that Epsilonproteobacteria are more active, even though their diversity is low (Glaubitz et al., 2010; Kirkpatrick et al., 2019). According to our data, Gammaproteobacteria significantly outnumber Epsilonproteobacteria, yet this does not allow us to speculate about the difference in chemoautotrophic activity of these two groups.

Hydrogen sulfide was absent at the majority of the stations within the suboxic layer. The abundance of sulfur‐oxidizing taxa is likely to be connected with the cryptic sulfur cycle, in which the processes of sulfate oxidation and reduction are linked and coupled with denitrification and anammox (Canfield et al., 2010; Hannig et al., 2007; Kirkpatrick et al., 2012). The presence of the cryptic cycle is argued by the presence of sulfate‐reducing Deltaproteobacteria: SAR324 (5%), Desulfarculaceae (0.6%), and Desulfobacteraceae (0.6%) that can occupy anaerobic niches of organic particulate matter (Fuchsman, Staley, et al., 2012). Together with green sulfur and non‐sulfur bacteria, these taxa are supposed to supply sulfur‐oxidizing bacteria with H2S.

Euxinic layer

Gammaproteobacteria (12%) decreased in the Black Sea euxinic zone being outnumbered by the representatives of AB16 class (17) and Deltaproteobacteria (13%) (Figure 5), which is in line with the prevalence of sulfur reduction process over sulfur oxidation in anoxic conditions. Epsilonproteobacteria increased in their abundance reaching 5%.

The taxonomic composition of the microbial communities indicates a smooth shift between the suboxic and euxinic zones: The upper part of the anoxic zone (83–159 m) seems to be an ecotone between the oxygen‐deprived and deep euxinic zone and is populated with communities similar to those from the suboxic layer. Notably, the deep anoxic water was populated with a different microbial community with a pronounced shift toward fermentation potentially carried out by the representatives of MSBL2 class (Spirochetes) (Daffonchio et al., 2006; Dong et al., 2018), which constituted 15% on average. The process of fermentation results in acetate, ethanol, and hydrogen production, which is beneficial for sulfate reducers, like Deltaproteobacteria that use hydrogen as a reductant in a deep anoxic environment (Dong et al., 2018). Hydrogen‐producing bacteria are also likely to support anaerobic respiration of organohalide respiring Dehalococcoides (Röling et al., 2007), which accounted for 5% within the deep euxinic environment. Similarly, Anaerolineae (Chloroflexi), which constituted a significant portion of the deep anoxic microbial community (14%) (Figure 5), benefit from the hydrogen‐producing bacteria. Even though the members of this genus are typically photoautotrophs and photoheterotrophs, they are capable of growing in the dark using hydrogen and H2S as electron donors, and a wide variety of organic substances, such as sugars or amino acids, as carbon source (Vuillemin et al., 2020).

As in the oxygen minimum zone, phototrophic Chlorobiaceae was among the most abundant families reaching 22% in the upper anoxic layer (Figure 6). This genus is likely to benefit from the remaining radiation at the depths of 80–110 m, with oxygen deprivation and H2S as necessary components no matter whether it is the suboxic or euxinic zone. Its maximum number was observed at the depth range of 82–99 m within suboxic and euxinic zones, and it was much less abundant at lower and higher depths across the non‐oxygenated water column. Indeed, its relative abundance dropped to 0.08% in the deep anoxic zone.

Several prokaryotic guilds that were detected within the suboxic zone were still abundant in the upper part of the euxinic layer despite it does not seem to be a fitting environment for them. It concerns the prokaryotes that require high‐potential electron acceptors like NO3 that are absent in this environment. Yet, the families comprising sulfur‐oxidizing chemoautotrophic bacteria—Piscirickettsiaceae and Helicobacteraceae—constituted the average of 6% of the community (Figure 6). As expected, the average abundance of these families decreased to 0.05% and 0.2% respectively at higher depth.

Sulfate‐reducing bacteria that act as important degraders of organic matter in anaerobic conditions (Jørgensen, 1982; Muyzer & Stams, 2008) were diverse and abundant within the whole euxinic zone, albeit they were in higher abundance at deeper layers. The group harbored a variety of sulfate‐reducing bacteria belonging to Desulfarculaceae (6%), SAR324 (3%), and Desulfobacteraceae (2%). Desulfarculaceae (6%) and Desulfobacteraceae (2%) were among the dominant families in the deep anoxic environment constituting 6% and 4%, respectively. SAR 324 abundance decreased with depth and accounted for 0.02%.

Heterotrophic representatives of the microbial community included A714017 (16%) and SHAS460 (2%) families that both belong to the SAR406 clade and have a potential for the reduction of sulfur compounds (Figure 6).

Notably, only 7% of the microbial community inhabiting the upper anoxic zone remained unidentified at the family level, whereas 30% of microbes from the deep euxinic zone could not be taxonomically assigned. This indicates that a large portion of euxinic microbial communities is still “a dark matter” and highlights the importance of future culture‐based studies targeting them (Suominen et al., 2019).

Spearman correlation analysis resulted in hierarchical clustering of Black Sea microbial communities in 6 groups (Figure 7), which illustrates their co‐occurrence and might be extrapolated to their common ecological niche requirements. The taxa that benefit from the sulfidic environment were as follows: Desulfobacteraceae, Helicobacteraceae, Chlorobiaceae, Piscirickettsiaceae, Desulfarculaceae, SHAS460, A714017 clustered together. The second cluster included prokaryotes that are capable of growing under anoxic or oxygen minimum conditions but require light–photoheterotrophic ZA3409c and Rhodospirillaceae, as well as Thiohalorhabdaceae that depend on O2 or NO3 as electron acceptors. The third group comprised Alteromonadaceaea, Microbacteraceae, Rhodobacteracea, and Saprospiraceae that are mostly particle‐attached, frequently found in algal blooms, and are known to be involved in complex polysaccharides degradation. The fourth cluster encompassed the majority of prokaryotic families, both heterotrophic and autotrophic, both free‐living and particle‐attached, mostly aerobic—Synechococcaceae, Flavobactericeae, Pelagibacteriaceae, and Balneolaceae. Copiotrophs adapted to a high concentration of phytoplankton and large‐molecular weight compounds—OSC155, Planctomycetacea Verrucomicrobiaceae clustered in the fifth group. Pseudomonadaceae and Pseudoalteromonadaceae clustered separately.

FIGURE 7

Hierarchical clustering of Black Sea prokaryotic families based on Spearman correlation

The differentiation of prokaryotic groups that are frequently found in surface water, in particular during algal blooms, might reflect their specialization on different growth substrates provided by algae. Indeed, both the third and the fourth group were found to correlate positively with the abundance of diatoms in the surface layer (r = 0.67, p = 0034 and r = 0.75, p = 0.001, respectively).

Functional composition of microbial communities populating the Black Sea water column

The pathways responsible for the general functions connected to nucleotide and amino acid metabolism, and genetic information processing comprised more than 50% of microbial functional repertoire. A set of orthologs specific for the most important functions in marine ecosystems was chosen for subsequent detailed analysis.

The distribution of microbial community functions was in line with the gradient of chemical parameters and with the taxonomic distribution of prokaryotes (Figure 2). According to PICRUSt, oxygenic photosynthesis (19%) and assimilatory sulfate reduction (40%–47%) dominated the microbial community inhabiting surface and thermocline. Assimilatory sulfate reduction is known to be present in most aerobic marine microorganisms, and it is a vital process as it ensures sulfur supply for the biosynthesis of S‐containing amino acids (Tripp et al., 2008).

Naturally, the PICRUSt ratio of oxygenic photosynthesis decreased with depth, whereas the ratio of anoxygenic photosynthesis increased reaching 7% at oxygen minimum and euxinic zones (Figure 8). Anoxygenic phototrophs typically inhabit a specific econiche, where light reaches sulfidic water layers and include facultatively anoxygenic cyanobacteria, Chlorobi, Chloroflexi, and purple bacteria (Overmann & Manske, 2006). This microbial group is known to benefit both from the remaining light and inorganic electron donors (e.g., nitrite, ferrous iron, molecular hydrogen, and reduced sulfur compounds) (Haas et al., 2018). The green sulfur bacteria, which are known for their special adaptations to the low‐light euxinic environment due to the presence of large photosynthetic antennae and high sulfide tolerance, have been previously detected within the Black Sea suboxic zone (Manske et al., 2005; Overmann & Manske, 2006) and high concentrations of bacteriochlorophyll‐e were observed in the suboxic and upper sulfidic zones (Ediger et al., 2019).

FIGURE 8

Functional pathways distribution in microbial communities inhabiting Black Sea water column

The genes responsible for the variety of processes in the nitrogen cycle have been detected in a high ratio as well (Figure 8). Indeed, according to PICRUSt, assimilatory nitrate reduction and dissimilatory nitrate reduction to ammonium (DNRA) accounted for 9% and 8% respectively within the surface layer and 8% and 7% within the thermocline. There was a trend of assimilatory nitrate reduction to decrease with depth, whereas the ratio of DNRA increased in oxygen minimum and anoxic conditions. This could be explained by the fact that assimilatory nitrate reduction involves using the ammonia from nitrate reduction as a nitrogen source (Shapleigh, 2009), which is an assimilation process strongly dependent on energy. This energy can be acquired via photosynthesis, which yields enough ATP for this process. As light and consequently photosynthesis becomes limited with depth, the energy required for nitrogen assimilation becomes scarce as well. Thus, the DNRA pathway is dominating over the assimilatory one at the suboxic and euxinic zone.

DNRA is a chemoheterotrophic anaerobic process, which results in NH4 + production and is therefore connected to nitrification and anammox (Van den Berg et al., 2017; Pajares & Ramos, 2019). DNRA and denitrification can occur simultaneously and rely on NO3− and carbon concentrations. Therefore, these two processes are assumed to compete for nitrate in the microbial community. It was traditionally proposed that the combination of electron donors (e.g., carbon) and electron acceptors (e.g., NO3 −) availability is the driving factor behind DNRA‐denitrification partitioning (Tiedje, 1988). However, contrasting data on the role of limiting factors in NO3 − to N2 versus to NH4 + are coming from culture‐based studies. Kraft et al. suggested that the NO3/NO2 and C/N ratios are influencing DNRA and denitrification partitioning, whereas Van den Berg et al. proposed that these processes are driven by more complex environmental factors, such as a combination of pH, sulfide concentrations, type, and complexity of electron donor (Van den Berg et al., 2017; Kraft et al., 2014). Contrary to the traditional understanding, DNRA can also be an autotrophic process driven by inorganic nutrients, in particular sulfide (S2−), elemental S, or Fe2+ (Slobodkina et al., 2017). According to our PICRUSt data, the DNRA pathway accounted for 7%–8% in the oxygenated normoxic zone and increased to 11%–12% in the oxygen‐depleted part of the water column. The presence of DNRA, as well as denitrification within the surface, thermocline, and DCM, can be explained by the microaerophilic environment, which could result from the particles coming from dead phytoplankton cells. Our PICRUSt data predicting significant DNRA pathway abundance are in contrast to the data from the oxic–anoxic interface of the Baltic Sea, where in situ measured DNRA rates were rather low (Bonaglia et al., 2016). Such a contrasting pattern should be further investigated in Black Sea microbial communities using incubation experiments, as the PICRUSt data are not comparable to the in situ measurements.

According to PICRUSt, denitrification accounted for 5%–7% in the oxygenated part of the water column and its role increased with depth reaching 12% in oxygen minimum and euxinic zone. PICRUSt predicted the presence of the complete denitrification pathway in the present dataset: narG, narZ, nxrA (nitrate reductase), nirK (nitrite reductase), norB (nitric oxide reductase subunit B), nosZ (nitrous oxide reductase) genes. Yet, no nirS gene was predicted in the studied microbial community. Several studies targeting microbial denitrification in marine OMZs have reported a significant prevalence of the nirK gene over nirS (Fuchsman et al., 2017; Pajares et al., 2019). However, both nirK and nirS genes responsible for the reduction of nitrite to nitric oxide have been previously detected in the Black Sea oxygen minimum zone (Kirkpatrick et al., 2012; Oakley et al., 2007). This discrepancy between our data and the previously obtained estimates might result from the fact that PICRUSt relies on the availability of the reference genomes and 16S rRNA gene sequences in the reference databases, whereas denitrifying bacteria have only been identified by their functional genes both in the Black Sea (Oakley et al., 2007 and Kirkpatrick et al., 2012) and in oxygen minimum zones (Fuchsman et al., 2017).

As predicted by PICRUSt, the nitrification pathway reached its maximum of 27% at the DCM layer and was much lower at the oxygen minimum and euxinic zone. This process had mostly been considered aerobic and involves 3 steps: ammonia oxidation to hydroxylamine (NH2OH), NH2OH oxidation to NO2 − (performed by ammonia‐oxidizing bacteria: Betaproteobacteria (Nitrosomonas and Nitrosospira) and Gammaproteobacteria (Nitrosococcus) (Pajares & Ramos, 2019) and ammonia‐oxidizing archaea: Thaumarchaeota phylum, e.g., Nitrosopumilus maritimus and Cenarchaeum symbiosum) (Lam et al., 2007)), and further conversion of NO2 − to NO3 − (performed by nitrite‐oxidizing bacteria (NOB): Chloroflexi, Nitrospirae, Nitrospinae, and several classes of Proteobacteria (Daims et al., 2016)). PICRUSt predicted the genes responsible for all steps of this pathway in our dataset; thus, the studied microbial community can potentially harbor both the prokaryotes performing certain steps and those capable of complete NH4 oxidation to NO3 − (comammox, complete ammonia oxidation) (Casciotti & Buchwald, 2012; Pajares & Ramos, 2019).

The maximum nitrification observed at DCM could result from the availability of phytoplankton cell lysis products within the deep chlorophyll‐a maximum zone. Phytoplankton cell lysis releases reduced organic N, which can be either assimilated or oxidized to NO3 − through nitrification (Casciotti & Buchwald, 2012). Microaerobic nitrification (ammonia oxidation) has been previously observed in the Black Sea suboxic zone providing the local source of nitrite for the anammox (Lam et al., 2007). The PICRUSt predicted nitrification (nitrite to nitrate oxidation) within suboxic and euxinic Black Sea zones is in line with the previous data from OMZs around the world, where this process was suggested to be performed by specialist nitrifiers with exceptionally high affinity for O2 (Bristow et al., 2016, 2017; Sun et al., 2017).

Intriguingly, according to our data, nitrification contributed much less to the nitrogen cycle within Black Sea surface waters. Photoinhibition might act as an indirect control factor for the nitrification process (Horak et al., 2018; Lam et al., 2007; Peng et al., 2016), and it has been suggested that slow‐growing nitrifying prokaryotes might be outcompeted by phytoplankton at the sunlit surface (Smith et al., 2014).

PICRUSt predicted nitrogen fixation pathway abundance increased with depth and reached its maximum of 11% within the euxinic zone (Figure 9). Diazotrophs capable of biological nitrogen fixation include microorganisms with varying metabolism: cyanobacteria, heterotrophic (e.g., Klebsiella, Vibrio), and phototrophic prokaryotes (e.g., Chlorobium, Chromatium, Rhodospirillum), and strict anaerobes (e.g., Clostridium, Desulfovibrio), all of which have nitrogenase complex (Pajares & Ramos, 2019). The prevalence of biological nitrogen fixation in oxygen‐depleted over surface Black Sea waters can putatively be explained by the sensitivity of nitrogenase complex to O2 (Pajares & Ramos, 2019), while a high ratio of denitrification within oxygen minimum and euxinic zone can enhance this process as well (Deutsch et al., 2007). This finding is in line with the earlier observations of active nitrogen fixation in the Black Sea dark suboxic and uppermost sulfidic zones (Kirkpatrick et al., 2019).

FIGURE 9

Spearman correlation of Black Sea microbial community functional pathways and environmental parameters (blue—positive correlation, red—negative correlation, only significant correlations are shown)

PICRUSt failed to detect any genes within the anaerobic ammonium oxidation pathway. Denitrification is prevailing over anammox in several OMZs (Bulow et al., 2010; Dalsgaard et al., 2003). However, the previous studies of anammox in Black Sea waters indicated that this process exceeded denitrification and could thus be considered the main fixed N2 sink in the central Black Sea (Jensen et al., 2008; Kuypers et al., 2003), and different species of anammox bacteria of Candidatus Scalindua genus have been found both at the upper and the lower suboxic zone (Fuchsman, Staley, et al., 2012).

Candidatus Scalindua has been detected in our data from the Black Sea suboxic zone (Sigma‐theta 15.8–16.0), yet its abundance was rather low (up to 0.5%), and it could have arisen from the typically limited capacity of the universal 16S rRNA gene primers to detect this specific genus (Fuchsman, Staley, et al., 2012; Harhangi et al., 2012; Yang et al., 2020). As a result, Candidatus Scalindua might have been underrepresented in our dataset, which prevented PICRUSt from the efficient mapping of its 16S rRNA gene sequences (Lincy & Manohar, 2020).

Genes responsible for transformations within dissimilatory sulfate reduction adenylylsulfate reductase (AprAB: sulfite <‐> adenylyl sulfate (APS) transformation) and dissimilatory sulfite reductase alpha (DsrAB: sulfide <‐> sulfite transformation) were predicted by PICRUSt in our dataset and comprised between 2% (in well‐oxygenated waters) and 7%–8% in oxygen minimum and euxinic zones. The previous research has already identified the AprAB and DsrAB are the key enzymes mediating dissimilatory sulfate reduction in Black Sea microbial mats (Basen et al., 2011). According to our results, dissimilatory sulfate reduction coincides with increased DNRA and denitrification processes within the oxygen‐depleted water column part (Figure 8), which is in line with the observations of coupled nitrogen and sulfur cycles in the Black sea (Fuchsman, Staley, et al., 2012) and the Gotland Deep of the Baltic Sea (Hannig et al., 2007). Even though dissimilatory sulfate reduction has been considered to be obligatory anaerobic (Jørgensen et al., 2019), we have detected the genes responsible for this pathway within the well‐oxygenated water column part. Yet, given their minor ratio of 2%, it can be suggested that dissimilatory sulfate reduction is occurring within anaerobic econiches (sinking particles) at the surface, thermocline, and DCM of the Black Sea (Fuchsman, Murray, et al., 2012; Rabus et al., 2006).

Microbial community differentiation based on taxonomic and functional profiles

Microbial community differentiation pattern was checked on both regional (Ukrainian shelf vs offshore zone comparison) and vertical scale (comparison of layers within normoxic, oxygen minimum, and euxinic zones) with the purpose to uncover the factors shaping the taxonomic and functional structure of the Black Sea prokaryotic communities.

Regional differentiation

According to the ANOSIM analysis, Ukrainian shelf and offshore samples clustered separately; therefore, these two zones were considered for further analysis. As the Ukrainian shelf is limited by the 100 m isobaths, there is no oxygen minimum and euxinic zone there. Thus, the regional comparison was performed within the surface, thermocline, and DCM layers. The taxonomic differentiation was examined on the class and family levels.

The highest number of significantly different classes, 14, was identified within the surface layer (Kruskal–Wallis, cutoff 0.05). The ratio of Saprospirae, Spartobacteria, and Mollicutes was significantly higher at the Ukrainian shelf, whereas Betaproteobacteria were more abundant in the offshore waters. Six significantly different classes were identified within the thermocline. More BME43 representatives were detected within the Ukrainian shelf and more Oscillatoriophycideae—in the offshore zone. DCM layers harbored 12 classes with abundance significantly different between Ukrainian shelf and offshore. More Verrucomicrobiae and Cytophagia were detected within the Ukrainian shelf, whereas Bacilli, SAR202, and Epsilonproteobacteria were more abundant in offshore waters.

Generally, the differentially abundant taxa were detected in low abundance <5%, whereas the dominant taxonomic group ratio was similar.

Naturally, more significantly different taxa were detected at the family level. 43 differentially abundant families were identified at the surface layer. Ukrainian shelf surface waters harbored more Rhodobacteraceae, Saprospiraceae, Microbacteriacea, and Caulobacteraceae representatives, whereas more Rhodospirillaceae, Piscirickettsiaceae, and Cyanobacteriaceae were present in offshore waters. 38 families with a different ratio in the two regions were detected within thermocline and included Saprospiraceae, Caulobacteraceae, Comamonadaceae, and Actinomycetaceae (more within Ukrainian shelf), as well as Marine Group II, Geobacteraceae, and Microbacteriacea (more in the offshore zone). Similarly, 38 significantly different families were identified within the DCM with more Rhodobacteraceae, Verrucomicrobiaceae, Sphingomonadaceae, and Microbacteriacea present in Ukrainian shelf waters. Notably, the ratio of Rhodobacteraceae, which is one of the dominant families in the studied communities, was significantly different in both surface and DCM of Ukrainian shelf and offshore waters.

Even though ANOSIM analysis did not reveal any significant functional differentiation between Ukrainian shelf and offshore waters, the Kruskal–Wallis (cut‐off 0.05) analysis detected some functional variation between the two regions, when surface, thermocline, and DCM layers were compared separately. Most functional differentiation was observed at the surface with 106 significantly different functional pathways including carbon fixation in prokaryotes, porphyrin and chlorophyll‐a metabolism, amino acid biosynthesis, nitrogen metabolism, carbon fixation in photosynthetic organisms, photosynthesis, carotenoid biosynthesis, and sulfur metabolism. 20 and 8 functional pathways with significantly different ratios were detected within thermocline and DCM layers, respectively.

The differences detected between the Ukrainian shelf and offshore might be attributed to the overall higher phytoplankton diversity and biomass detected at the Ukrainian shelf compared to the offshore zone at the moment of sampling, which influenced the type and availability of substrates for bacterioplankton communities (Slobodnik et al., 2017).

Vertical differentiation

Taking into account the fact that our study was limited with only one shelf station sampled, microbial community vertical differentiation analysis was conducted for the offshore zone only. Kruskal–Wallis analysis identified 132 prokaryotic classes and 292 families with significantly different abundances across the Black Sea water column. A detailed description of microbial community taxonomic structure was provided earlier.

In total, 288 significantly different functional pathways were detected across the Black Sea water column by Kruskal–Wallis test based on PICRUSt analysis results (Table A3). As expected from the ANOSIM results, no significantly different pathways were detected between surface and thermocline. The highest number of significantly different pathways was detected when surface and thermocline were compared against oxygen minimum and euxinic zones. The differentially abundant pathways included both those connected to the general metabolism (amino acids and nucleotide metabolism, transporters, etc.) and those responsible for more specific functions (photosynthesis, N and S metabolism, xenobiotic degradation, etc.).

A Mantel test was conducted with the purpose to corroborate the leading role of environmental factors over geographic distance in shaping microbial communities of the Black Sea. The depth, O2, and H2S concentration and nutrient concentration (PO4, NO3, Ptotal, Ntotal, and Si were shown to be significant in explaining microbial communities’ variation when the Mantel test was performed on taxonomic data (Table 1). Less environmental factors explained the data distribution when functional pathways were considered. Yet, when a specific set of genes responsible for nitrogen, sulfur turnover, and photosynthesis, was analyzed, the depth, O2, H2S, and nutrient concentration (PO4, NO2, Ntotal, and Si) were indicated as significant drivers of microbial community structure (Table 1). NH4 was not measured during the cruise; thus, there was no possibility to include this very important factor in the analysis.

TABLE 1

Environment impact on Black Sea microbial community structure and functions.

Factor	Taxonomy, class	Taxonomy, genus	Functions (all)	Functions (specific)
Geographical distance	0.0065	0.19	0.11	0.01
Depth	0.66*	0.66*	0.14*	0.48*
O2	0.5*	0.54*	0.16*	0.41*
H2S	0.23*	0.22*	0.55*	0.54*
PO4	0.29*	0.3*	0.06	0.25*
NO2	0.0019	0.02	−0.07	0.14*
NO3	0.0023*	0.19*	−0.02	−0.09
Ptotal	0.4*	0.41*	0.007	0.083
Ntotal	0.13*	0.11*	0.05	0.38*
Si	0.41*	0.44*	0.11*	0.22*

Significant impact (p < 0.05).

The insignificant correlation of geographical distance matrix with taxonomic and functional data indicates that the Black Sea microbial communities are influenced primarily by environmental factors that are in turn shaped by pronounced water column vertical stratification.

Spearman correlation analysis of the microbial functional structure and environmental parameters revealed that depth, O2, PO4, and NO3 concentration are the most important factors (Figure 9). Anoxygenic photosynthesis, nitrogen fixation, and dissimilatory sulfate reduction clustered together being abundant within the oxygen‐depleted zone. Anoxygenic photosynthesis is known to be fueled by inorganic electron donors, such as hydrogen sulfide, which is the product of dissimilatory sulfate reduction. There is evidence of sulfate reduction and sulfide oxidation occurring as parallel processes in OMZs as a part of a cryptic sulfur cycle (Canfield et al., 2010; Carolan et al., 2015; Lavik et al., 2009). The sulfide oxidation process is coupled with nitrate reduction resulting in N2 production, which in turn might be responsible for the high nitrogen fixation ratio in oxygen‐depleted waters. This is supported by high nitrate concentrations detected within the suboxic zone (Table A1).

Assimilatory nitrate reduction, nitrification, and assimilatory sulfate reduction pathways clustered together as well (Figure 9). Heterotrophic nitrification of organic matter, which could come from phytoplankton cell lysis, results in nitrate release, which is reduced back to ammonia via the process of assimilatory nitrate reduction. This statement is supported by high NO3 concentration within DCM (Table A1).

Denitrification, dissimilatory nitrate reduction, and oxygenic photosynthesis clustered separately according to their correlation patterns.

CONCLUSIONS

A large‐scale study of Black Sea prokaryotic diversity revealed pronounced taxonomic and functional differentiation of microbial communities in accordance with water column stratification. According to PICRUSt, a shift has been observed from primary producers and heterotrophs inhabiting well‐oxygenated and irradiated waters toward specialized taxa capable of anoxygenic photosynthesis, sulfate reduction, and fermentation and benefiting from the suboxic and euxinic environment. Regional differentiation was detected between microbial communities inhabiting oxygenated water layers in the Ukrainian shelf and offshore zone, even though it was less pronounced than the vertical one. It was suggested that the regional microbial community structure can be driven by the presence of various phytoplankton taxa providing different substrates for the prokaryotes in the two regions. At the same time, vertical differentiation is likely to be mostly driven by oxygen concentration. The data represent a comprehensive 16S rRNA gene dataset for prokaryotes inhabiting the specific Black Sea environment to date and thus it can be a useful resource for further investigations targeting specific taxonomic groups, metabolic pathways, or biogeochemical cycles.

CONFLICT OF INTEREST

None declared.

AUTHOR CONTRIBUTIONS

Mariia Pavlovska: Formal analysis (lead); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing‐original draft (lead). Ievgeniia Prekrasna: Investigation (equal); Methodology (equal); Validation (equal); Visualization (supporting); Writing‐original draft (supporting). Evgen Dykyi: Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal). Andrii Zotov: Formal analysis (supporting); Investigation (equal); Validation (equal). Artem Dzhulai: Investigation (supporting); Validation (supporting); Visualization (supporting). Alina Frolova: Investigation (supporting); Software (supporting); Validation (equal). Jaroslav Slobodnik: Formal analysis (supporting); Project administration (equal); Resources (lead); Writing‐original draft (supporting). Elena Stoica: Formal analysis (equal); Methodology (equal); Project administration (supporting); Resources (equal); Supervision (equal); Validation (supporting); Writing‐original draft (supporting).

ETHICS STATEMENT

None required.

TABLE A2

Sequencing depth and OTU number per sample

Sample	Description	# Sequences	# OTUs
4.1	Surface	112818	234
6.1	Surface	122869	247
11.1	Surface	165791	243
12.1	Surface	120003	251
10.1	Surface	101991	249
9.1	Surface	112285	270
8.1	Surface	66287	279
7.1	Surface	135270	240
3.1	Surface	115907	261
2.1	Surface	97151	256
1.1	Surface	61639	203
5.1	Surface	108026	262
5.2	Thermocline	118665	240
6.2	Thermocline	106002	239
11.2	Thermocline	126090	243
4.2	Thermocline	131693	252
12.2	Thermocline	66894	243
10.2	Thermocline	108356	213
9.2	Thermocline	92497	284
8.2	Thermocline	148730	253
7.2	Thermocline	114979	233
3.2	Thermocline	116231	253
2.2	Thermocline	110929	268
1.2.1	Thermocline	94118	241
1.2	Thermocline	84905	299
6.3	Deep chlorophyll‐a max	103662	326
11.3	Deep chlorophyll‐a max	156306	304
12.3	Deep chlorophyll‐a max	96783	319
12.3	Deep chlorophyll‐a max	77813	335
10.3.1	Deep chlorophyll‐a max	66791	290
10.3.2	Deep chlorophyll‐a max	54802	352
4.3	Deep chlorophyll‐a max	121075	289
9.3	Deep chlorophyll‐a max	87731	270
8.3	Deep chlorophyll‐a max	103168	294
7.3	Deep chlorophyll‐a max	74679	315
3.3	Deep chlorophyll‐a max	81552	299
2.3	Deep chlorophyll‐a max	113777	274
1.3	Deep chlorophyll‐a max	79565	306
1.3	Deep chlorophyll‐a max	82511	274
5.3	Deep chlorophyll‐a max	119752	340
6.4	Suboxic	139774	259
11.4.1	Suboxic	117976	351
11.4.2	Suboxic	100963	312
12.4	Suboxic	119660	253
10.4.1	Suboxic	80871	362
10.4.2	Suboxic	66927	324
9.4.1	Suboxic	99754	304
9.4.2	Suboxic	107021	302
4.4	Suboxic	100684	312
8.4.1	Suboxic	84677	299
8.4.2	Suboxic	92986	258
4.5	Suboxic	108267	230
7.4	Suboxic	127979	350
3.4	Suboxic	81275	299
5.4	Suboxic	142752	259
2.4	Suboxic	44524	324
5.4	Suboxic	121457	307
5.5	Euxinic	103587	339
10.5.1	Euxinic	90805	246
10.5.2	Euxinic	79906	339
10.5.3	Euxinic	91763	326
9.5	Euxinic	116374	249
8.5	Euxinic	100811	274
7.5.1	Euxinic	95154	368
7.5.2	Euxinic	108401	310
3.5.1	Euxinic	90621	341
3.5.2	Euxinic	89964	337

TABLE A3

Functional pathways with significantly different abundance based on the results of the Kruskal–Wallis test (FDR p < 0.05)

	DCM	Suboxic	Euxinic
Surface	166: For example, transporters, amino acid metabolism, carbon fixation pathways in prokaryotes, N and S metabolism, xenobiotics degradation	175 For example, transporters, amino acid metabolism, photosynthesis–antenna proteins, in prokaryotes, xenobiotics degradation, beta‐Lactam resistance	211: transporters, amino acid metabolism, photosynthesis–antenna proteins, in prokaryotes, xenobiotics degradation, beta‐Lactam resistance, N and S metabolism
Thermocline	140: For example, amino acid metabolism, N and S metabolism, xenobiotics degradation	154: For example, amino acid metabolism, photosynthesis–antenna proteins, carbon fixation in photosynthetic organisms, xenobiotics degradation	203: For example, transporters, amino acid metabolism, photosynthesis–antenna proteins, in prokaryotes, N and S metabolism, xenobiotics degradation, beta‐Lactam resistance
DCM	—	37: For example, nucleotide metabolism, photosynthesis–antenna proteins, xenobiotics degradation, beta‐Lactam resistance	88: For example, nucleotide metabolism, photosynthesis–antenna proteins, xenobiotics degradation, beta‐Lactam resistance
Suboxic		—
Euxinic			—