Bacterioplankton Diversity and Distribution in Relation to Phytoplankton Community Structure in the Ross Sea Surface Waters

Primary productivity in the Ross Sea region is characterized by intense phytoplankton blooms whose temporal and spatial distribution are driven by changes in environmental conditions as well as interactions with the bacterioplankton community. However, the number of studies reporting the simultaneous diversity of the phytoplankton and bacterioplankton in Antarctic waters are limited. Here, we report data on the bacterial diversity in relation to phytoplankton community structure in the surface waters of the Ross Sea during the Austral summer 2017. Our results show partially overlapping bacterioplankton communities between the stations located in the Terra Nova Bay (TNB) coastal waters and the Ross Sea Open Waters (RSOWs), with a dominance of members belonging to the bacterial phyla Bacteroidetes and Proteobacteria. In the TNB coastal area, microbial communities were characterized by a higher abundance of sequences related to heterotrophic bacterial genera such as Polaribacter spp., together with higher phytoplankton biomass and higher relative abundance of diatoms. On the contrary, the phytoplankton biomass in the RSOW were lower, with relatively higher contribution of haptophytes and a higher abundance of sequences related to oligotrophic and mixothrophic bacterial groups like the Oligotrophic Marine Gammaproteobacteria (OMG) group and SAR11. We show that the rate of diversity change between the two locations is influenced by both abiotic (salinity and the nitrogen to phosphorus ratio) and biotic (phytoplankton community structure) factors. Our data provide new insight into the coexistence of the bacterioplankton and phytoplankton in Antarctic waters, suggesting that specific rather than random interaction contribute to the organic matter cycling in the Southern Ocean.

Primary productivity in the Ross Sea region is characterized by intense phytoplankton blooms whose temporal and spatial distribution are driven by changes in environmental conditions as well as interactions with the bacterioplankton community. However, the number of studies reporting the simultaneous diversity of the phytoplankton and bacterioplankton in Antarctic waters are limited. Here, we report data on the bacterial diversity in relation to phytoplankton community structure in the surface waters of the Ross Sea during the Austral summer 2017. Our results show partially overlapping bacterioplankton communities between the stations located in the Terra Nova Bay (TNB) coastal waters and the Ross Sea Open Waters (RSOWs), with a dominance of members belonging to the bacterial phyla Bacteroidetes and Proteobacteria. In the TNB coastal area, microbial communities were characterized by a higher abundance of sequences related to heterotrophic bacterial genera such as Polaribacter spp., together with higher phytoplankton biomass and higher relative abundance of diatoms. On the contrary, the phytoplankton biomass in the RSOW were lower, with relatively higher contribution of haptophytes and a higher abundance of sequences related to oligotrophic and mixothrophic bacterial groups like the Oligotrophic Marine Gammaproteobacteria (OMG) group and SAR11. We show that the rate of diversity change between the two locations is influenced by both abiotic (salinity and the nitrogen to phosphorus ratio) and biotic (phytoplankton community structure) factors. Our data provide new insight into the coexistence of the bacterioplankton and phytoplankton in Antarctic waters, suggesting that specific rather than random interaction contribute to the organic matter cycling in the Southern Ocean.

INTRODUCTION
Antarctica and the Southern Ocean are central to Earth's climate and oceanic circulation systems (Convey and Peck, 2019). Primary productivity in this region is characterized by intense phytoplankton blooms whose temporal and spatial distribution are driven by different environmental conditions, although the mechanisms regulating these processes are still poorly known (Moore and Abbott, 2000;Deppeler and Davidson, 2017). In this context, the Ross Sea is one of the most productive sectors in Antarctica (Smith and Gordon, 1997;Arrigo et al., 1999), with an annual productivity averaging ∼180 g C m −2 year −1 (Arrigo et al., 2008). Phytoplankton blooms here have been well-studied, with the dominance of diatoms and haptophytes presenting different temporal and spatial patterns (Smith et al., 2014;Mangoni et al., 2017). In the last years, however, changes in phytoplankton blooms and dynamics that contrast the classical Antarctic paradigm have been observed as for example the dominance of Phaeocystis antarctica in stratified coastal waters during summer, the high levels of biomass in an area usually considered an High Nutrient Low Chlorophyll (HNLC), or the presence of relative higher percentage of minor functional groups contrasting the classical Antarctic paradigm (DiTullio et al., 2000;Phan-Tan et al., 2018;Mangoni et al., 2019;Bolinesi et al., 2020b). All together, these observations require a revaluation of the phytoplankton-bacteria interaction, which could be playing a key role in structuring the trophodynamics in this area (Bertrand et al., 2007(Bertrand et al., , 2015. The assumption that the Antarctic Ocean microbial communities were generally species poor has been re-discussed in recent years, since many results suggested that microbial diversity is significantly higher than previously recognized (Murray and Grzymski, 2007;Wilkins et al., 2013;Silvi et al., 2016). In the last decade, researchers have demonstrated that in Antarctic pelagic food webs micro-eukaryote dynamics significantly contribute to the structuring of the prokaryotic community (Kirchman et al., 2001;Wilkins et al., 2013;Delmont et al., 2014 and references therein) and in turn prokaryotic diversity can influence phytoplankton productivity under certain conditions (Bertrand et al., 2007(Bertrand et al., , 2015. Prokaryotic communities contribute to or dominate several key ecosystem processes, including primary production, the turnover of biogenic elements, the mineralization of the organic matter, and the degradation of xenobiotics and pollutants (Cole, 1982;Azam et al., 1983;Croft et al., 2005;Azam and Malfatti, 2007;Falkowski et al., 2008;Sher et al., 2011).
Early work in Antarctic waters revealed that bacterial biomass might represent up to 30% of total microbial biomass in coastal areas (Fiala and Delille, 1992). In Antarctica, free-living marine microbial community composition can differ significantly between locations at relatively small spatial and temporal scales, responding to environmental variations of temperature, salinity, nutrients, or the presence of oceanic fronts (see for example, the Supplementary Material of Lozupone and Knight, 2007;Raes et al., 2018). For example, in Terra Nova Bay (TNB), substantial differences in terms of bacterial assemblages have been observed between coastal and offshore stations and along the water column (Celussi et al., 2009). Heterotrophic members of the Alphaproteobacteria and Gammaproteobacteria class of the Proteobacteria have been reported as dominant phylotypes in Antarctic waters Murray and Grzymski, 2007;Wilkins et al., 2013), and studies in these and other marine ecosystems indicate that bacterial growth is frequently dependent on phytoplankton-derived DOM (Church et al., 2000;Morán et al., 2001Morán et al., , 2002Piquet et al., 2011;Ducklow et al., 2012;Kim et al., 2014). A broad diversity among the class Flavobacteria has been also reported in different sub-areas of the Southern Ocean (Abell and Bowman, 2005). According to Piquet et al. (2011), melt water stratification and the transition to non-stabilized Antarctic surface waters may have an impact not only on microeukaryotes but also on bacterial community composition, with a shift from an Alpha-and Gammaproteobacteria to a Cytophaga-Flavobacterium-Bacteroides-dominated community under mixed conditions. Some studies found an unusual presence of strictly anaerobic Epsilonproteobacteria (now reclassified as phylum Campylobacterota) in the bottom of sea ice (Gentile et al., 2006), probably as a consequence of the oxygen decay and sulfide accumulation, caused by high degradation rates of sympagic diatoms by aerobic and anaerobic heterotrophs (Brierley and Thomas, 2002).
Shifts in the Eukaryotic plankton communities have been often reported as rapidly followed by a shift in the bacterial community (Billen and Becquevort, 1991;Piquet et al., 2011), with a series of phytoplankton-bacterial interactions resulting in both positive and negative feedback loops (Bertrand et al., 2015). For example, in temperate waters phytoplankton blooms in spring and summer induce changes in bacterioplankton community structure (Fuhrman et al., 2006;Teeling et al., 2016;Chafee et al., 2018). These community dynamics have been shown to be recurrent, indicating a phytodetritus-driven seasonality and suggesting that the phytoplankton-prokaryotic interactions in surface waters are more sophisticated than previously thought (Seymour et al., 2017). Mechanisms addressing the nature of the mutualistic interaction between phytoplankton and bacterioplankton communities have been proposed over time. For example, the exchange of phytoplankton exudates and bacteria-produced cobalt containing vitamin B 12 represents one of the best studied feedback loops (Bertrand et al., 2007(Bertrand et al., , 2011Bolinesi et al., 2020b). The production and release of vitamin B 12 by bacteria, in fact, depends on the degradation of phytoplankton exudates, establishing a complex feedback mechanism between prokaryotic and phytoplanktonic communities by heterotrophic bacteria (Fang et al., 2017). Similar trophic interaction between phytoplankton and bacterioplankton in the marine ecosystem, and especially in polar regions, might be more common than previously known. The bacterial-phytoplankton interaction may thus evolve in a series of complex relationship, affecting directly or indirectly the micronutrient availability and co-limitation (e.g., iron, cobalt, and vitamin B 12 ; Bertrand et al., 2007;Tagliabue et al., 2017;Bolinesi et al., 2020a) as well as nutrient uptake (Amin et al., 2015;Bertrand et al., 2015). Yet, despite this the number of studies reporting the simultaneous diversity of the phytoplankton and bacterioplankton community in the Antarctic waters is Frontiers in Microbiology | www.frontiersin.org 3 January 2022 | Volume 13 | Article 722900 comparatively few (Di Poi et al., 2013;Flaviani, 2017;Richert et al., 2019). In this study, we analyzed the bacterial diversity in relation to phytoplankton in the surface waters of TNB and Ross Sea, providing new insight into the coexistence of the two communities in Antarctic waters.

Sampling Procedure and Study Site
Seawater samples were collected either at the deep chlorophyll maximum or at ~20 m depth in stations where the fluorescence profile showed a rather homogeneous distribution of chlorophyll-a (Chl-a) in the upper layer ( Table 1). Sampling activities were carried out on the R/V Italica during the Austral Summer 2017 (between 13 and 30 Jan 2017), in the framework of Plankton biodiversity and functioning of the Ross Sea ecosystems in a changing Southern Ocean (P-ROSE) and CDW Effects on glacial mElting and on Bulk of Fe in the Western Ross sea (CELEBeR) projects -Italian National Antarctic Programfunded by the Ministry of Education, University and Research (MIUR). The area of investigation falls within two different zones (Figure 1) of the Ross Sea, the coastal area of TNB and the Ross Sea Open Water (RSOW). Water samples were collected using a carousel sampler (Sea-Bird Electronics 32) equipped with 24 12-L Niskin bottles and a conductivitytemperature-depth (CTD) instrument (9/11 Plus; Sea-Bird Electronics), along a transect from the coastal area of TNB to the oper Ross Sea (Figure 1C), crossed by a north-south aligned secondary transect carried out in the RSOW area. For the bacterial diversity analysis, at each station, 500 ml of seawater was collected from the Niskin bottle and filtered shipboard onto a 0.22 μm filter (Whatman, 47 mm diameter) successively stored at −80°C and transported back to the lab. For the analysis of total phytoplankton biomass, 500 ml of seawater was filtered shipboard through 0.45 μm GF/F filter (Whatman, 47 mm diameter). For the determination of size classes of phytoplankton, 500 ml of seawater was prefiltered on board onto a 20 or 2 μm net, and the flow through filtered on a 0.45 μm GF/F (Whatman 47 mm diameter) obtaining two size fractionated subsamples. All GF/F filters were preserved frozen at −80°C and transported back to the lab for analysis (see below). The contribution of the different size classes was calculated as the difference between total phytoplankton biomass the two size fractions to provide micro-(>20 μm), nano-(between 20 and 2 μm), and pico-(<2 μm) phytoplankton biomass. For the determination of phytoplankton functional groups by chemotaxonomic criteria, 2 L of seawater were filtered on board onto 0. ), water samples were taken directly from the Niskin bottles and stored at −20°C in 20 ml low-density polyethylene containers until laboratory analysis. Monthly mean sea surface Chl-a concentrations in log(mg m −3 ) at 4 km resolution derived from the MODIS-AQUA sensor (Satellite remote sensing Ocean color data) were downloaded from the European Data Portal (Melin, 2013).

Phytoplankton Community Structure and Biomass
Frozen samples were processed in Italy for the determination of Chl-a and phaeopigments (Phaeo-a) content (used as proxy for phytoplankton biomass) using a solution of 90% acetone according to Holm-Hansen et al. (1965), with a spectrofluorometer (Shimadzu) checked daily with a Chl-a standard solution (from Sigma-Aldrich). HPLC pigments separation was performed on an Agilent 1100 HPLC according to the method outlined in Vidussi et al. (1996) as modified by Mangoni et al. (2017). The following biomarker pigments were used as chemotaxonomic descriptors: alloxanthin (cryptophytes), chlorophyll b (chlorophytes), prasinoxanthin (prasinophytes), 19′-butanoyloxyfucoxanthin (pelagophytes), fucoxanthin (diatoms), 19′-hexanoyloxyfucoxanthin (haptophytes), peridinin (dinophytes), and zeaxanthin (Cyanobacteria). The contribution of the main phytoplankton groups to the total Chl-a was estimated on the basis of the concentrations of biomarker pigments using the chemical taxonomy software CHEMTAX (Mackey et al., 1996). The phytoplankton species composition and cells abundance were determined following the Uthermöhl method (Utermöhl, 1931), according to which at least 400 cells were counted per sample with an inverted light microscope (LM) Zeiss Axiovert Observerz. One at 400× magnification and used to estimate specific group abundance. In order to better visualize the shifts in the phytoplankton community, we computed a diatoms to haptophytes ratio (D/H ratio), calculated using the relative abundance of diatoms and haptophytes as: (diatoms/haptophytes)/(diatoms + haptophytes). Phytoplankton photosynthetic efficiency was estimated using the maximum photochemical quantum yields of PS II (Fv/Fm), representing the initial maximum efficiency of photons captured by open PSII reaction centers. The maximum photochemical quantum yields of PS II were measured using a Phyto_PAM II compact unit (waltz) as described in Bolinesi et al. (2020b).
All samples were acclimated in the dark before analysis to minimize the non-photochemical dissipation of excitation, and measurements were blank corrected by filtering the sample through a 0.2 μm filters. For determining Fv/Fm, samples were illuminated with a saturating pulse, and the ratio was calculated using the formula Fv/Fm = (Fm − F 0 )/Fm as previously described.

Nutrient Analysis
Inorganic nutrients were analyzed using a five-channel continuous flow autoanalyzer (Technicon Autoanalyser II), according to the method described by Hansen et al. (1983) adapted to the available instrumentation. Briefly, samples flow was controlled by a multi-channel peristaltic pump, which regulates the flow rate of samples and reagents throughout the analytical procedures. The sample flow was segmented with air bubbles to enhance mixing of reagents and sample, and to reduce smearing. The sample-reagent mixture reacts chemically to produce color in proportion to the concentration of the nutrient in the sample, and is analyzed in a flowcell using a phototube as detector.

Community DNA Extraction
Total community DNA was extracted as previously reported (Giovannelli et al., 2013) with slight modifications. Briefly, each filter was washed with 1.7 ml of extraction buffer solution [100 mM NaCl (pH 8.0), 20 mM EDTA, and 50 mM Tris-HCl (pH 8.0)] added with 10 μl Proteinase K (1 mg/ml) and incubated at 37°C for 1 h with occasionally mixing. A 100 μl volume of 10% sodium dodecyl sulfate (SDS) was added to each sample, followed by incubation at 55°C for 1 h. The liquid phase was collected and centrifuged for 15 min at 10,000×g to separate cellular debris from nucleic acids. The supernatant was extracted twice with an equal volume of phenol, followed by a precipitation with 2.5 volume of ethanol 100% and 0.1 volumes of sodium acetate (3 M) for 12 h at −20°C. The DNA was collected by centrifugation, washed in 70% cold ethanol, air dried, and resuspended in 50 μl of sterile distilled water. The integrity of DNA was assessed spectrophotometrically and by PCR amplification of the 16S rRNA gene using the primes Ribo-For (5′-AGTTTGATCCTGGCTCAG-3′) and Ribo-Rev (5′-CCTACGTATTACCGCGGC-3′) (Fakhry et al., 2008). The PCR products were visualized on 1% agarose gel stained with ethidium bromide.

16S rRNA Gene Sequencing
Partial 16S rRNA gene sequences were obtained using primer pair Probio_Uni/Probio_Rev, which target the V3 region of the 16S rRNA gene sequence (Milani et al., 2013). 16S rRNA gene sequencing was performed using a MiSeq platform (Illumina) at the DNA sequencing facility of GenProbio srl 1 according to the protocol previously reported (Milani et al., 2013). The used primers were tested using Silva TestPrime 1.0 against the SSU r132 allowing 0-2 mismatched and showed a good coverage and specificity for the Bacteria domain (96.1% coverage and 73.3% specificity), while having a low specificity for the Archaea (91.3% coverage and 15.2% specificity). Given the low specificity for Archaea of the primer set used in this study, we have obtained results for the bacterial population only.

Bioinformatics and Statistical Analyses
All sequences were imported in R, and analyzed with the DADA2 package (Callahan et al., 2016). Following the package guidelines, quality plots were performed to check the sequences' quality. Post-QC reads were trimmed using the filterAndTrim command [truncLen = c(155,145), maxN = 0, maxEE = c(2,2), truncQ = 2, rm.phix = TRUE, trimLeft = 17, and trimRight = 15]. After this step, a parametric error model, based on the convergence between the estimation of error rate and the inference of the sample composition, was performed. Paired-end reads were merged and exact Amplicon Sequence Variants (ASVs) inferred using the dada algorithm. Chimeric sequences were removed and prokaryotic taxonomy assigned using the naive Bayesian classifier method against the Silva Database (r132; https://www.arb-silva.de/documentation/release-132/ After the Phyloseq object creation, low abundance ASVs (less than three reads across the dataset), Mitochondria, Chloroplast, Eukaryotes sequences, and potential contaminants (Sheik et al., 2018) were removed. The resulting dataset was represented by 703 unique ASVs and 535,009 reads. ASVs counts were normalized to the median library size across the dataset. Diversity analyses were carried out using the Phyloseq package (McMurdie and Holmes, 2013) with relative abundance set to 100% after the removal of sequences described above. Top abundance ASVs and Genera were defined has having a cumulative relative abundance above 0.1% in our dataset. The alpha diversity was investigated using both the Simpson and Shannon diversity index among the two sampled areas. The beta diversity was investigated using the UNIFRAC and Jaccard diversity index as implemented in the vegan package (Oksanen et al., 2012). Both the abundance weighted and unweighted version of the index was used. The resulting similarity matrix was plotted using non-metric multidimensional scaling techniques implemented in the ordination command of Phyloseq. The resulting ordination was used to investigate correlations with environmental and phytoplankton variables using the envfit and ordisurf functions in vegan. Collinearity among the predictors was checked using a Pearson correlation matrix. The linearity of the correlation between the rate of change in the beta diversity and the variables identified as significant by the envfit function was checked by plotting the non-metric multidimensional scaling (nMDS) axis against the variable. Statistically significant differences in the distribution of abundant bacterial genera were tested using the Chi-square test. Co-correlation networks were calculated as a pairwise distribution of each ASV across the entire dataset using Spearman rank correlation and different ρ cutoff selected for network plotting using the igraph package (Csardi and Nepusz, 2006).

Diversity of the Bacterial Community
Bacterial diversity was investigated using the 16S rRNA gene sequence. After quality check and data filtering, a total of 535,009 reads were obtained and used to identify 703 unique ASVs. Simpson and Shannon diversity index showed higher diversity in the stations of RSOW than TNB (Figure 3), albeit the differences were not statistically significant (Kruskal-Wallis test). Despite these differences, the bacterial community structure at the phylum level was similar among the stations of the two areas. This apparent similarity is still visible all the way to the family level (Figure 4). Sequences belonging to the Bacteroidetes and Proteobacteria represented the most abundant phyla in all the samples, with on average 50.1% of the reads assigned to the Bacteroidetes and 48.4% to the Proteobacteria, respectively. In addition, sequences classified as belonging to the phyla Cyanobacteria, Firmicutes, Actinobacteria, and Verrucomicrobia were detected in almost all samples. The phylum Bacteroidetes was dominated by the class of Bacteroidia, with the Flavobacteriales representing the most abundant order. Among Flavobacteriales, several ASVs were assigned to the genera Polaribacter, Aurantivirga, and Brumimicrobium, representing the top genera in the studied area ( Figure 5A). Within the Proteobacteria, most sequences were affiliated to the class Gammaproteobacteria (47.41%), abundant both in the coastal (TNB) and in the offshore (RSOW) stations, followed by a lower percentage of ASVs assigned to the class Alphaproteobacteria (about 1% of total ASVs). Gammaproteobacteria were mainly represented by the orders Oceanospirillales, more abundant in the stations of the TNB area than in RSOW area (23.59% of ASVs found in TNB vs. 16.74% found in RSOW), followed by Cellvibrionales and Alteromonadales, which instead were more abundant in the offshore area ( Figure 5B). Within the Oceanospirillales, several species were affiliated with the genera Alcanivorax, Oleispira, Halomonas, and Profundimonas, all known hydrocarbon degraders, distributed almost all the sampled stations, while among the Cellvibrionales, the main representatives were members of the clades SAR92 and OM60(NOR5). In addition, among Alteromonadales, we found mainly sequences classified as members of the Pseudoalteromonas, Marinobacter, and Colweilla genera. Alphaproteobacteria were mainly represented by the order SAR11_Clade and Rhodobacterales, uniformly distributed in both the studied areas representing on average 0.55 and 0.33% of the total reads, respectively ( Figure 5C).
While the abundance of possible hydrocarbonoclastic bacteria (identified based on the similarity to known isolates) was relatively low, representing at most 3% of the total community (Supplementary Figure 1), they were distributed differentially among the sampled areas. Members of the genera Alcanivorax, Marinobacter, and Oleispira were present in all the sampled stations, with a higher abundance in the offshore stations of RSOW with, on average, 0.55, 0.46 , and 0.29% of the total reads, respectively. Albeit with a lower abundance, we found also reads classified as Methylophaga and Methylobacillus in several stations of the RSOW area, while completely absent in the TNB stations.  To explore the bacterial community composition between the stations and to identify the environmental drivers responsible for its structuring, we investigated the beta diversity distribution using both abundance weighted and unweighted dissimilarity indices. The principal coordinate analysis (PCoA) using either the abundance weighted and unweighted Unifrac distance index shows a clear separation between the two studied areas in both plots (Supplementary Figures 2A,B). Abundance weighted estimates of the beta diversity show a total of 72.5% of the variance visualized on the plot (Supplementary Figure 2A), with the two study areas clearly separating along the PCoA axis 2 (explaining 21% of the variance). The only exception to this separation is represented by station ST43, which represents the first station of the RSOW area but is geographically located near the TNB stations (Figure 1). Unweighted diversity measures show a less pronounced separation among the two sampled areas, with a lower percentage of the variance accounted for by PCoA axis 1 and 2 (cumulatively 29% of the variance; Supplementary Figure 2B). Direct comparison of the two plots revealed that the stations had on average similar species, and that differences were due variations in the abundant species.
To investigate the role of measured environmental parameters and the phytoplankton community composition in driving the bacterial diversity in the two studied areas, we performed nMDS ordination based on the abundance weighted Jaccard dissimilarity measure followed by the environmental and community composition vector fitting (Figure 6). Similarly to the results obtained with the abundance weighted PCoA analysis, the nMDS showed a clear separation among the two areas with the exception of station ST43. Linear vector fitting against the nMDS ordination Frontiers in Microbiology | www.frontiersin.org 9 January 2022 | Volume 13 | Article 722900 revealed that the main factors explaining the bacterial diversity were the geographic location (represented by the longitude of the sampled station) and the haptophytes relative abundance, with correlation coefficients of R 2 = 0.86 and R 2 = 0.72 against nMDS1, respectively (Figure 6; Table 5). The nitrogen to phosphorus ratio (N/P ratio), was instead strongly correlated with nMDS2 (R 2 = 0.73), together with salinity and the maximum photosynthetic quantum yield (Fv/Fm), showing a correlation coefficient of R 2 = 0.69 and R 2 = 0.46, respectively. Several other variables were identified by the linear vector fitting as significatively correlated with nMDS axis 1 and 2. This was likely due to the high degree of collinearity present  among the environmental parameters and among the phytoplankton composition descriptors (Figure 7). Collinearity was investigated using Pearson moment correlation among the predictors used for the linear vector fitting, together with the results of the vector fitting for both nMDS axis (Figure 7; Table 5).

Differential Distribution of Key Bacterial Genera in the Ross Sea Surface Waters
The overall distribution of the bacterial community in the multivariate ordination shows that the two areas are different in term of community assemblages and that the separation between the two sampled areas is driven by abundant species. We determined unique and shared ASVs between TNB and RSOW area (Figure 8A), revealing that the shared community is composed of 365 ASVs, while 132 and 206 unique ASVs are found in TBN and RSOW samples, respectively. We identified the top bacterial genera shared among the two studied areas ( Figure 8B), and identified those showing a differential distribution. The results show several well-known members of the Antarctic bacterioplankton community as the most abundant genera, some of which are differentially distributed among the two areas. Sequences belonging to the Polaribacter genus were highly abundant in the dataset, representing on average 47.1% of the total reads. Among them, 41.1% were related to Polaribacter irgensii (Table 6), a psychrophilic heterotrophic gas vacuolate bacterium of the Bacteroidetes phylum, previously isolated from Arctic and Antarctic sea ice (Gosink et al., 1998). As shown in Figure 8B, reads classified as P. irgensii are differentially   (Sunamura et al., 2004;Sunamura and Yanagawa, 2015) or in low oxygen areas of the water column (Walsh et al., 2009), follow a trend similar to P. irgensii, with the percentage of ASVs assigned to this taxon higher in the coastal area of TNB with respect to the offshore RSOW stations (4.1 vs. 2.1%, respectively). Although with a lower abundance (less than 1% of the total reads), we found that the ASVs associated to the genera Profundimonas and Brumicrobium, cold-adapted and facultatively anaerobic heterotrophic bacteria generally found in marine environments (Bowman et al., 2003;Margesin, 2017), are also more prevalent in the TNB samples.
An opposite trend was followed by the reads classified as members of the genera Aurantivirga, and, with a markedly lower abundance, by the reads affiliated to Clade Ia, OM60 (NOR5) clade, Alcanivorax, Marinobacter, and Halomonas ( Figure 8B). Aurantivirga is a member of Flavobacteriaceae whose type strain, Aurantivirga profunda, has been isolated in the deep-sea waters of the Pacific Ocean (Song et al., 2015). Our data show that the ASVs belonging to Aurantivirga are abundant in the RSOW area (on average 4.7%) and become substantially lower in the coastal area of TNB (0.7%, adj.p < 0.01).
Our results show that the ASVs related to Clade Ia are significantly more abundant in the RSOW area (0.7%) with respect to the TNB area (0.3%). A similar distribution was found for members of the OM60 (NOR5) clade with a higher abundance in the offshore stations with respect to the coastal stations (0.6% in RSOW vs. 0.2% in TNB). The genera containing obligate or facultative hydrocarbon oxidizers Alcanivorax, Oleispira, and Marinobacter were also differentially distributed among the sampled areas, showing statistically significant higher abundances in the stations of the RSOW compared to TNB (adj.p < 0.001 for Alcanivorax and adj.p < 0.01 for Oleispira and Marinobacter).
Patterns of co-occurrence among the more prevalent ASVs (here defined as occurring in at least 20% of the sampled stations) were investigated using different cut-off levels of Spearman correlation coefficient, between 0.5 and 0.85 (Figure 9;  Supplementary Table 1). The results revealed a phylogenetically heterogeneous cluster of highly co-occurring ASVs present at low threshold of correlation (from 0.5 to 0.65 ρ), progressively breaking down into small (20-30 ASVs) phylogenetically heterogeneous highly correlated clusters.

DISCUSSION
The Ross Sea is a complex mosaic of subsystems with physical, chemical, and biological features that respond to biotic and abiotic forcing at different temporal and spatial scales (Smith Frontiers Mangoni et al., 2017;Bolinesi et al., 2020b). Although the phytoplankton blooms dynamics in the Ross Sea have been well described, the drivers regulating the temporal and spatial distribution of the blooms are still debated (Arrigo et al., 1999;DiTullio et al., 2000;Mangoni et al., 2017Mangoni et al., , 2019Saggiomo et al., 2021). Among the biotic processes structuring the phytoplankton community, the interactions with the bacterioplankton has been proposed to be an important factor (Richert et al., 2019). For example, vitamin B12 availability, which is produced exclusively by bacterioplankton, has been suggested to be one of the drivers involved in the diatoms-Phaeocystis shift in the Ross Sea waters (Bertrand et al., 2007(Bertrand et al., , 2015. Despite the bacterioplankton composition in the Ross Sea has been previously investigated (Lo Giudice et al., 2012;Silvi et al., 2016), the phytoplanktonbacterioplankton co-occurrence in the area has been poorly characterized. Many studies have suggested the existence of sophisticated ecological interactions between phytoplankton and bacteria in driving marine biological processes, which can range from a mutualistic exchange of biomolecules and nutrients to the competition for the same limiting inorganic nutrients (Amin et al., 2012(Amin et al., , 2015Seymour et al., 2017). In this study, we have investigated the bacterial diversity of the Ross Sea in relation to phytoplankton community structure in 20 stations located between the coastal area of TNB and the offshore waters of the central Ross Sea (RSOW). The two areas were characterized by different amounts of phytoplankton biomass and dominant functional groups (Figure 2), with highest levels of biomass near the coastline. Pseudo-nitzschia spp. and Fragilariopsis spp. were the most abundant species in TNB, while P. antarctica dominated the stations in the RSOW. The bacterioplankton community was representative of Antarctic surface waters in accordance with previous studies (Abell and Bowman, 2005;Ghiglione and Murray, 2012;Grzymski et al., 2012;Wilkins et al., 2013), with Bacteroidetes and Proteobacteria representing the most abundant phyla. Among Bacteroidetes, members of the Flavobacteriaceae, commonly found bacteria in polar environments (Abell and Bowman, 2005;Williams et al., 2012), dominated both the sampled areas. These bacteria are generally classified as the key degraders of phytoplankton derived organic matter (Teeling et al., 2012;Buchan et al., 2014). Their presence in all the sampled stations suggests that these bacteria might play an important role in the organic carbon cycle of the Ross Sea, potentially impacting

A B
FIGURE 8 | Shared and unique Amplicon Sequence Variants (ASVs) among the two sampled areas (A) and density function of the abundance of the top genera shared among the sampled areas (B). Density plots showing the differential distribution of the major genera between the two sampled areas. The x axis represents the relative abundance of each ASV in the genera while the y axis shows the relative density functions. The vertical solid line represents the mean relative abundance for those genera in the area. Differentially abundant genera are marked with an * according to the result of the Kruskal-Wallis test (*adj.p < 0.05; **adj.p < 0.01; and ***adj.p < 0.001).
Frontiers in Microbiology | www.frontiersin.org 15 January 2022 | Volume 13 | Article 722900 organic carbon transfer to higher trophic levels (Azam and Malfatti, 2007). Despite the apparent similarity in the community composition between the two areas (Figure 4), several differences were present at the Genus and ASVs level. Within our dataset, members related to the Polaribacter genus were dominant in all the sampled stations, and showed a statistically higher abundance in the TNB area ( Figure 8B). The closest relative to our sequences was P. irgensii, with an average similarity of 98.7% (Table 6), a known psychrophilic heterotrophic marine bacterium. Previous studies have reported that Polaribacter species thrive during diatom-blooms (Teeling et al., 2012;Williams et al., 2012), potentially suggesting that their presence in the TNB area might be linked with the higher presence of diatom species we have identified. Similarly, a recent analysis revealed that species of the Aurantivirga genus are among the first prokaryotic taxa responding to diatom bloom in the Southern Ocean (Liu et al., 2020). Sequences related to the genus Aurantivirga are highly abundant in our TNB stations, suggesting once again a possible relationship with the high abundance of diatoms (Figures 2, 8B). Among the Proteobacteria, sequences related to the Gammaproteobacteria were detected in all the sampled stations, with a differential distribution of members related to the Alteromonadales, Cellvibrionales, and Oceanospirillales orders among the two studied areas. Within the Alteromonadales, sequences related to Pseudoalteromonas, Marinobacter, and Colweilla genera were dominant and in higher abundance in the RSOW stations ( Figure 8B). These bacterial genera comprise cold-adapted marine bacteria generally detected in the Southern Ocean waters able to degrade simple sugar, amino acids, organic compounds, and hydrocarbons (Médigue et al., 2005;Methé et al., 2005;Singer et al., 2011;Rosenberg, 2013). The abundance of ASVs belonging to these psychrophilic bacterial groups in the offshore area could be explained by specific adaptation to the conditions found in offshore stations. Similarly, members of the Cellvibrionales order dominated the offshore area. These bacteria, generally affiliated to the Oligotrophic Marine Gammaproteobacteria (OMG) cluster (Cho and Giovannoni, 2004), are able to adapt to nutrient depletion and carbon limitation with the potential to harvest light for mixotrophic growth (Stingl et al., 2007;Spring and Riedel, 2013;Courties et al., 2014). We also detected sequences related to the clades SAR92 and OM60 (NOR5).
Previous reports indicate that members of SAR92 clade dominate during the phytoplankton bloom in the Southern Ocean (West et al., 2008) and are able to establish a close relationship with the productive P. antarctica during the austral summer in the Amundsen sea polynya, suggesting an important role of this bacterial taxa during bloom formation and bloom longevity (Delmont et al., 2014). Consistent with this data, our results show a dominance of the SAR92 clade in the RSOW area which is dominated by P. antarctica. The ability to exploit different metabolic pathways based on the conditions found in the sea water, suggest that both Alteromonadales and Cellvibrionales can play an important role in the microbial FIGURE 9 | Co-occurrence network analysis drawn with increasing Spearman correlation cut off and colored according to the phyla classification of each ASVs.
Frontiers in Microbiology | www.frontiersin.org 16 January 2022 | Volume 13 | Article 722900 food web, contributing to the functioning of the Antarctic marine ecosystems. Taken together, these observations suggest a tight coupling between the phytoplankton community structure and the dominant bacteria identified in the Ross Sea surface waters. Among the top genera identified in our study, members belonging to the SAR11 order of the Alphaproteobacteria Clade Ia, were comparatively low in relative total abundance. Our results show that ASVs classified in this group were more abundant in the RSOW area. Species belonging to the SAR11 order comprise aerobic and free-living oligotrophic chemoheterotrophic bacteria globally distributed in marine environments (Morris et al., 2002). They are believed to contribute significantly to the carbon, nitrogen, and sulfur cycling in the Ocean (Malmstrom et al., 2004), and have been previously reported worldwide (Field et al., 1997;Morris et al., 2002;Brown et al., 2012). Members of the SAR11 group have been shown to be more abundant in the Subantarctic and polar fronts compared to Antarctic zones (Wilkins et al., 2013 and references therein). The low global abundance found in our dataset might be due to the lack of competitive advantage of SAR11 members in the presence of high molecular weight organic carbon during phytoplankton blooms. Members of Clade Ia identified in our dataset represent a specific subgroup of the order generally reported in cold waters (Brown et al., 2012;Delmont et al., 2019). The closest relative to our sequences was Pelagibacter ubique, the most common heterotrophic bacteria found in the ocean (Giovannoni, 1990;, with an average similarity of 100% (Table 6). Studies based on genome sequences analysis and in situ hybridization, revealed that P. ubique is an oligotrophic bacterium with a small genome size and a high metabolic activity, able to assimilate either dissolved free amino acids and dimethylsulfoniopropionate (DMSP; Malmstrom et al., 2004;Giovannoni and Stingl, 2005). DMSP is an organosulfur compound produced by several phytoplankton cells, which can perform a double function in polar microalgae, as osmolyte or cryoprotective agent (Kirst et al., 1991;Karsten et al., 1996). Several reports indicate that P. antarctica is one of the leading producers of DMSP in the Ross Sea (DiTullio and Smith, 1995;Ditullio et al., 2003).
Interestingly, our dataset reveals the presence of obligate or facultative hydrocarbon oxidizers in the Ross Sea surface waters, albeit at abundances below ~3%. While the presence of facultative hydrocarbon oxidizers per se is not indicative of hydrocarbon contamination in the environment, our data reveal that the obligate hydrocarbonoclastic genera Alcanivorax and Oleispira (Yakimov et al., 1998(Yakimov et al., , 2003 were present in all the sampled stations, and markedly more abundant in the RSOW stations. Alcanivorax and Oleispira are two bacterial genera found globally in extremely low abundances (Cafaro et al., 2013), but can become transiently dominant, with relative abundances up to 70-90% of prokaryotic cells, in the presence of hydrocarbon spills (Kasai et al., 2002a,b). Their presence in our dataset is interesting and can be explained in several different ways. While their abundance is significantly lower than reported after oil spillage events (Kasai et al., 2002b;Harayama et al., 2004), it is possible that hydrocarbon and exhaust oils released by the large number of ships transiting the Ross Sea every summer are responsible for keeping them above the normal background levels. Antarctic tourism has been steadily increasing over the years, with over half a million tourist landings reported for the 2017-2018 season (Palmowski, 2020), and fishing activities in proximity of the productive Antarctic waters have also increased (Brooks and Ainley, 2017). Alternatively, previous studies have proposed that members of the Alcanivorax genera are able to maintain viable status in uncontaminated marine waters degrading natural lipids of bacterial and phytoplankton origin, released in the water column due to exudation, sloppy feeding, or viral lysis (McGenity et al., 2012;Lea-Smith et al., 2015;Zadjelovic et al., 2020). During our sampling, the phytoplankton community was undergoing a summer bloom, and thus all the mentioned processes, e.g., exudates, sloppy feeding by the zooplankton grazers, and viral induced lysis, might have contributed in increasing lipid concentrations in seawater. In addition to Alcanivorax and Oleispira, the facultative hydrocarbon degraders Marinobacter and Halomonas have been also identified (Supplementary Figure 1). Members of the genus Marinobacter are slightly or moderately halophilic, able to degrade both aliphatic and aromatic hydrocarbons. Marinobacter spp. capable of growing on hydrocarbons as the sole carbon source has been previously isolated from sediments (Gauthier et al., 1992). The genus Halomonas comprises marine halophilic and/or halotolerant bacteria, known to produce large quantity of exopolysaccharides (EPS) with rheological and active-surface properties (Calvo et al., 1998Martínez-Checa et al., 2002;Pepi et al., 2005;Gutierrez et al., 2020). It is possible that Halomonas-producing EPS provides a mechanism to increase the bioavailability of hydrophobic compounds (e.g., hydrocarbons and lipid aggregates) that Alcanivorax, Oleispira, and Marinobacter utilize for growth.
Information regarding the structure of bacterial communities can be also investigated through the use of network analysis. Most microbial network analysis uses a single correlation threshold to identify meaningful interactions from a correlation matrix (Barberán et al., 2012;Williams et al., 2014;Giovannelli et al., 2016;Connor et al., 2017). This approach, while might provide useful information regarding the biological interactions in the community, requires a priori justification or a sensitivity analysis to demonstrate the robustness of the conclusion with respect to the selected threshold. In addition, single threshold co-correlation analysis is considered controversial by some authors (Hirano and Takemoto, 2019;Blanchet et al., 2020) and it is believed to increase the number of detected false positives. An alternative approach to identifying biological interactions is using co-occurrence to identify ecological response to common environmental factors, an approach recently applied with success across broad ecological gradients (Delgado-Baquerizo et al., 2018;Fan et al., 2018; Frontiers in Microbiology | www.frontiersin.org 17 January 2022 | Volume 13 | Article 722900 Fullerton et al., 2021). Microbial co-occurrence networks can be also used to investigate community assembly and dynamics. To this end, we used an increasing correlation threshold to identify the structuring of the bacterioplankton community in the surface waters of the Ross Sea, preserving more ecological signal compared to a fix-threshold approach (Figure 9). With this approach, the degree at which the network breaks apart and changes in structure during thresholding can be used to identify processes influencing community assembly (Connor et al., 2017). As the threshold is increased, the largest network component becomes sparser as edges among nodes are removed (Figure 9). If modularity and other network statistics follow a linear trend during thresholding, this suggests the presence of a dominant core microbiome, often proposed to be critical or keystone components of the community (Faust and Raes, 2012;Mandakovic et al., 2018). Conversely, non-linear trends in network modularity suggest the absence of representative core microbiome and possibly the absence of keystone species. Our data show that the bacterioplankton community is composed of phylogenetically heterogeneous clusters of bacteria showing progressively lower co-occurrence as the correlation threshold is raised (Figure 9). Usually, high levels of co-occurrence in microbial networks are interpreted as increasing biotic interactions among species, including competition for resources (Faust and Raes, 2012;Widder et al., 2014). The resulting sparse network we have obtained at higher level of correlation threshold (>0.60 ρ,  Table 1), suggests that stochastic processes (e.g., neutrality, dispersion, and physical stress, see Holyoak et al., 2005), rather than competition, are responsible for community assemblage in the surface waters of the Ross Sea. Additionally, our results suggest the absence of a defined core microbiome in the studied area. This might be connected to different bacterial assemblages being associated to distinct phytoplanktonic communities identified in the TNB and the RSOW areas. Despite this, the presence of small clusters (20-30 ASVs) of phylogenetically diverse but recurrent ASVs suggests a strong functional redundancy among the most representative bacterial species. Ecological network theory predicts that communities of tightly connected species should be more fragile. Our data suggest that the higher fragmentation identified might be the dynamic response of the bacterial community to phytoplankton related dynamics, which continuously rework and redistribute microbial niches as blooms succession progresses (Luria et al., 2017).

CONCLUSION
The role of bacterioplankton in Antarctic waters has been reevaluated in the last decades, highlighting the importance of the bacterioplankton-phytoplankton interaction in the ocean carbon cycling, biogeochemical processes, and trophic food chain. Data from our study indicate that bacterioplankton diversity was clearly linked to the phytoplankton community structure and distribution in the Ross Sea surface waters in the austral summer 2017. Members of the Bacteroidetes and Proteobacteria phyla were abundant in all the sampled stations, with a dominance of heterotrophic bacterial species (i.e., Polaribacter genus) in presence of diatom blooms in the coastal area of TNB, and a dominance of oligothrophic and mixotrophic bacterial species in presence of haptophytes blooms in the offshore area. The overall picture emerging suggests the existence of specific rather than random interactions between dominant phytoplankton-bacterioplankton species that contributes to the organic matter cycling in the Southern Ocean.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
OM, AC, and PR designed the study. MM, AC, PR, MSa, FB, and OM performed laboratory analyses. MM, AC, DG, MSe, MB, RM, and GD'E performed bioinformatics and statistical analyses. All authors contributed equally to the manuscript writing and revision.

ACKNOWLEDGMENTS
We express our gratitude to Italian Antarctic National Program (PNRA) and the scientific personnel and crew of the research vessel Italica for logistical support.