Physical and Biogeochemical Regionalization of the Southern Ocean and the CCAMLR Zone 48.1

The Southern Ocean plays a major role in the Earth’s climate, provides fisheries products and help the maintenance of biodiversity. The degree of correspondence between physical and biogeochemical spatial variability and regionalization were investigated by calculating the main physical factors that statistically explained the biogeochemical variability within the Southern Ocean and the 48.1 zone of the Commission for the Conservation of Antarctic Marine Living Resources (CCAMLR). The mean value of physical and biogeochemical variables was estimated during austral summer within a grid of 1° × 1° south of 50°S. The regionalization was developed using both non-hierarchical and hierarchical clustering method, whereas BIO-ENV package and distance-based redundancy analysis (db-RDA) were applied in order to calculate which physical factors primarily explained the biogeochemical spatial variability. A total of 12 physical and 18 biogeochemical significant clusters were identified for the Southern Ocean (alpha: 0.05). The combination of bathymetry and sea ice coverage majorly explained biogeochemical variability (Spearman rank correlation coefficient: 0.68) and db-RDA indicated that physical variables expressed the 60.1% of biogeochemical variance. On the other hand, 14 physical and 16 biogeochemical significant clusters were identified for 48.1 CCAMLR zone. Bathymetry was the main factor explaining biogeochemical variability (Spearman coefficient: 0.81) and db-RDA analysis resulted in 77.1% of biogeochemical variance. The correspondence between physical and biogeochemical regions was higher for CCAMLR 48.1 zone with respect to the whole Southern Ocean. Our results provide useful information for both Southern Ocean and CCAMLR 48.1 zone ecosystem management and modeling parametrization.


INTRODUCTION
The Southern Ocean plays a major role in the Earth's climate since the opening of the Drake Passage approximately 41 Ma ago (Scher and Martin, 2006). The Antarctic Circumpolar Current is the major current of the planet (with a transport of 137 Sv; Rintoul and da Silva, 2019) and it connects the Pacific, Indian and Atlantic Ocean basins isolating the Antarctic continent from lower latitudes.
The Antarctic Circumpolar Current and the Southern Ocean influence heat and salt balances and biogeochemical cycles (e.g., C, N, P; Sallée, 2018). The combination of the thermohaline gradient, sea-ice cycle and atmospheric patterns creates a meridional overturning circulation in the Southern Ocean (Tomczak and Godfrey, 2003;Morrison et al., 2015) which is a sensitive mechanism able to influence the greenhouse gases budget between the atmosphere and the ocean (Reid et al., 2009;Gruber et al., 2019). The export of Particulate Organic Carbon (POC) south of 50 • S is 1000 Tg C yr −1 corresponding to 10% of the global POC export, with a significant contribution of dissolved organic carbon in the top 500 m of the water column (Schlitzer, 2002). The Southern Ocean has also been identified as a region providing important ecosystem services, particularly fisheries (CCAMLR, 2018;FAO, 2018) and maintenance of biodiversity (Grant et al., 2013;David and Saucède, 2015).
The Southern Ocean is recognized as a High Nutrient Low Chlorophyll region, approximately 48% of the Southern Ocean mean surface Chlorophyll-a (Chl-a) concentration is lower than 0.25 mg m −3 despite high macronutrients concentration (Pollard et al., 2006). Different hypotheses have been proposed, from light limitation due to deep mixing (Venables and Moore, 2010) to zooplankton grazing pressure (Le Quéré et al., 2016) and bioavailable dissolved iron depletion (Martin et al., 1990;Wadley et al., 2014). The iron limitation hypothesis leaded to several expeditions between 1993 and 2011 to artificially fertilize the photic layer of the Southern Ocean, these experiments increased surface Chl-a up to 15-fold (Boyd et al., 2007;Yoon et al., 2018). Pitchford and Brindley (1999) suggested that different explanations should be regarded as complementary, rather than alternative. The oceanic islands of the Southern Ocean (e.g., South Georgia Islands) stand out from High Nutrient Low Chlorophyll productivity patterns, mainly due to local input of dissolved iron and shallower bathymetry (Venables and Moore, 2010).
The 48.1 zone of the Commission for the Conservation of Antarctic Marine Living Resources (CCAMLR) is a key region for monitoring the Antarctic ecosystem. This area is experimenting strong interannual variability and trends in sea ice coverage and water column temperature Jones et al., 2016) that turn this zone in one of the most climate sensitive region of the planet and most variable (Hendry et al., 2018). Recent studies (Murphy et al., 2016;Comiso et al., 2017) suggested that the interannual oceanographic variability is partly linked with atmospheric patterns, such as El Niño-Southern Oscillation and the Southern Annular Mode. This region is subject to strong anthropic influence via krill fisheries (Boyd, 2009;Ainley and Pauly, 2014;McBride et al., 2014) and local contamination by increasing tourist's presence (Bargagli, 2008;Waller et al., 2017). Intense Antarctic krill (Euphausia superba) fisheries within the 48.1 zone (Brooks et al., 2018) might have important repercussion on the trophic structure and the trophic pathways due to the reduced functional redundancy of polar ecosystems (i.e., few species perform most of the ecological functions; Murphy et al., 2016), the key role of krill within the Antarctic food web (Trivelpiece et al., 2011;Murphy et al., 2012;McBride et al., 2014) and the importance of this region for krill spawning and nursery (Piñones and Fedorov, 2016;Henley et al., 2019;Perry et al., 2019). A specific focus was given to CCAMLR 48.1 zone because of the ecological importance of this climate sensitive area and since no other biogeochemical regionalization study has been previously developed for this region according to our knowledge.
Regionalization is an important tool for conservation and management of the marine environment and marine resources, since it can provide useful information for identification of protected areas and fisheries zones (Berline et al., 2014). Marine protected areas are multi-purpose regions that provide protection for the natural resources they contain and can help marine ecosystems be resilient to natural and anthropic-driven climate changes (Roberts et al., 2017). Currently, the Ross Sea and the South Orkney Islands are the only established marine protected areas within the Southern Ocean, while three more proposal have been submitted for East Antarctica, the Antarctic Peninsula and the Weddell Sea (Sylvester and Brooks, 2020). The terrestrial ice-free area of Antarctica was divided in 16 biogeographic regions by Terauds and Lee (2016), whereas various studies of bioregionalization have been previously performed for the Southern Ocean pelagic (i.e., Grant et al., 2006;Spalding et al., 2007Spalding et al., , 2012Cabré et al., 2016) and benthic communities (Griffiths et al., 2009;Pierrat et al., 2013;Douglass et al., 2014;Hogg et al., 2016;Teschke et al., 2016;Fabri-Ruiz et al., 2020). A new physical and biogeochemical regionalization is needed to fill the existing gap in the understanding of the physical forcing on biogeochemical cycles spatial heterogeneity (Hendry et al., 2018) and to generate useful information for marine protected areas proposals and marine ecosystem modeling.
This study will generate new physical and biogeochemical regionalization utilizing the Southern Ocean data available from satellite-derived and model output databases. The degree of similarity between physical and biogeochemical characteristics obtained with this new regionalization will be compared with previously existing regionalization. In addition, the physical variable that primarily control biogeochemical variability in the Southern Ocean and CCAMLR 48.1 zone will be identified. The work will deliver useful information to the scientific community in order to identify highly productive areas and to understand and predict future changes in the context of the current climate crisis.

Study Zone
The International Hydrographic Organization (IHO, 2002) defined the northerner limits of the Southern Ocean at 60 • S. For the purpose of this study, the northern limit of the Southern Ocean was set at 50 • S due to the spatial distribution of the Polar and Subantarctic Fronts (Orsi et al., 1995;Sallée et al., 2008;Koshlyakov and Tarakanov, 2011), especially in the Atlantic and western Indian Ocean (Figure 1). The proposed limit allowed the inclusion of the entire Drake Passage and the Southern Patagonia, permitting the analysis of differences and similarities  Orsi et al. (1995), whereas the circulation patterns were reproduced combining the information from Heywood et al. (2004); Thompson et al. (2009), Armitage et al. (2018; Moffat and Meredith (2018), Thompson et al. (2018), and Vernet et al. (2019).
between Antarctic and Subantarctic habitats from a physical and a biogeochemical oceanographic point of view.
The limits of the CCAMLR 48.1 zone (50-70 • W, 60-70 • S with the exclusion of the of the southern part of the Weddell Sea) are in fact very similar to those that define the Western Antarctic Peninsula and Scotia Arc zone identified by the Southern Ocean Observing System 1 . The 48.1 zone include the shelves and pelagic areas off the Western Antarctic Peninsula from Margarite Bay to the tip of the Peninsula (Figure 1), coastal islands such as the South Shetland Islands and the northern part of the Eastern Antarctic Peninsula. This study dedicated a special focus on the CCAMLR 48.1 zone because it is one of the most important spawning regions for Antarctic krill (Perry et al., 2019) that is directly affected by fisheries (Ainley and Pauly, 2014;Brooks et al., 2018) and has experienced one of the fastest temperature increase of the planet (Turner et al., 2005(Turner et al., , 2016. Krill was the main fishery target within the 48.1 zone and represented the 99.9% of the catches in the area since 2000 (CCAMLR, 2018). Krill landings from the 48.1 zone accounted for 45% of total Southern Ocean krill catches and the eastern and western part of the Bransfield Strait were the regions where most of the 48.1 krill catches occurred (62%), followed by the western pelagic area of the Drake Passage (15%). Most of these catches were made with midwater otter trawls (CCAMLR, 2018).

Physical and Biogeochemical Variables
The variables selected for the physical and biogeochemical regionalization, along with their characteristics, are presented in 1 http://www.soos.aq/activities/rwg/wapsa Table 1. The mean value of the physical and biogeochemical parameters during the "productivity season" (December-March; Supplementary Figures 1, 2) was calculated in a 1 • × 1 • latitude/longitude grid south of 50 • S. The study focused on this period based on the phenology of phytoplankton proliferations (Soppa et al., 2016), the higher krill abundances (Atkinson et al., 2017) and the higher proportion of valid satellite data with respect to austral winter. The majority of the analyzed variables covered a time span of 15 years (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018), with the significant exception of iron concentration and mixed layer depth (MLD) estimations that only were available for a period of 5 years ( Table 1). These model-output databases were chosen despite of their limited time frame because they guaranteed a higher spatial coverage with respect to others in situ databases (i.e., Tagliabue et al., 2012;Li et al., 2017;Holte et al., 2017), especially in the coastal areas of the Antarctic continent. Olson et al. (2016) bathymetry data was obtained from the data.gov portal 2 , whereas Titchner and Rayner (2014) sea ice estimations were downloaded from the Met Office Hadley Center page 3 . Monthly composite of satellite-derived variables (sea surface temperature, photosynthetically active radiation and chlorophyll-a; SST, PAR and Chl-a, respectively) were extracted from the CoastWatch project of the National Oceanic and Atmosphere Administration 4 . Model output variables, such as MLD and iron concentration in the surface water, were obtained from the Southern Ocean State Estimation page 5 . Net Primary Production (NPP) estimations calculated with the "standard" Vertically Generalized Production Model (Behrenfeld and Falkowski, 1997) were acquired from the Ocean Productivity page 6 . Finally, Kostadinov et al. (2016) nanophytoplankton concentration was achieved from the PANGEA data center 7 and Particulate Organic Carbon (POC) exports were extracted from the Australian Antarctic Data Center 8 . It is important to consider that satellite-derived Chl-a was directly used as input parameter for NPP, diatoms and POC export estimations, which may cause redundancy in the spatial variability of these variables and the biogeochemical regionalization. In spite of this, the treatment of Chl-a data was different for each product, since only diatoms concentration exclusively relies on Chl-a concentration and NPP and POC estimations used multiple input variables along with Chla data. Soppa et al. (2016) equation was used for diatoms concentration: log 10 (Diatom) = 1.1559·log 10 (Chl-a) − 0.2901. Behrenfeld and Falkowski (1997) equation for NPP estimations might be summarized as:

represents a function of surface
Chlorophyll-a estimation, f (SST) indicates a function of Sea Surface Temperature, D L is the day length and f(PAR) represents a function of Photosynthetically Active Radiation. On the other hand, POC export flux was calculated by Guillaumot et al. (2017) using: POC export = NPP c · f (Z eu ) · f (Z 0 ); where NPP c is the Net Primary Production in the surface water calculated with the Carbon-based Production Model (Westberry et al., 2008), Z eu indicates the depth of the euphotic zone and Z 0 the depth of the sea floor. The simplified equation for NPP c might be summarize as:

Estimates of Regionalization
Regionalization is a useful methodological tool that allows to identify similarities and dissimilarities between adjacent and remote areas according to different criteria and to group regions with analogous properties. Indeed, a regionalization developed using physical variables might provide the degree of similarity (expressed as distance in a dendrogram) between regions from a physical standpoint. On the other hand, the use of biogeochemical variables will result in a regionalization according to a biogeochemical point of view. The number of regions provided by this analysis is an indicator of heterogeneity; whereas the mean value for each group (i.e., cluster) might help the identification of areas with similar temperature and sea ice patterns as well as high and low productivity zones. The comparison between physical and biogeochemical regionalization can provide useful information about physical and biogeochemical patterns within a defined area.
The Southern Ocean regionalization was developed following the methodology of Grant et al. (2006) and Raymond (2014). The main steps for this analysis were: (I) The mean bathymetry for each pixel was transformed in log 10 (bathymetry); (II) the pixels with at least a missing value among all variables were eliminated (i.e., no interpolation was used); (III) all the variables were normalized between 0 and 1; (IV) non-hierarchical Clustering Large Application (CLARA; Kaufman and Rousseeuw, 1990) was applied to group all the pixels within 250 clusters using Manhattan distance for both datasets as indicated by Grant et al. (2006) and Raymond (2014); (V) the mean value for each of the 250 groups was calculated; (VI) an agglomerative hierarchical clustering method (Unweighted Pair Group Method with Arithmetic mean, UPGMA; Sokal and Michener, 1958) was applied to the results of the previous step; (VII) the number of significant clusters within the data was identified using the Clarke et al. (2008) statistical method (1000 permutations and an alpha value of 0.05 with the assumption of no a priori groups). This section slightly modified the original methodology proposed by Grant et al. (2006), which incorporated a cutoff method for dendrograms. Indeed, the cutoff technique strongly depends on authors criteria, whereas other "traditional" methods to determine numbers of clusters within a dendrogram (i.e., elbow, silhouette, gap statistic and Charrad et al. (2014) methods) were unsuccessful when applied to the data (Supplementary Table 1). On the contrary, the regionalization of CCAMLR 48.1 zone was performed exclusively with UPGMA agglomerative cluster due to reduced sample size. The analyses were carried out using R 9 and Matlab 10 .
The regionalization methodology applied in this study is accompanied by different sources of uncertainties in the pixels grouping and boundaries location between areas. Among different factors, uncertainty might be ascribed to imprecision in the data and epistemic uncertainty related to the knowledge of regionalization process (Grant et al., 2006). Data imprecisions comprise satellite or model estimations errors and bias due to incomplete observations, whereas knowledge errors include the data selection and the processing steps, such as the clustering method.

Relationship Between Physical and Biogeochemical Variables and Regions
The Pearson linear correlation coefficient (p-value: 0.05) between physical and biogeochemical variables, along with the geographical coordinates, was estimated. The BIO-ENV package from Clarke and Ainsworth (1993) was used in order to find the best subset of physical variables that have maximum correlation with biogeochemical dissimilarities throughout Spearman rank correlation coefficient. In addition, Distance-Based Redundancy Analysis (db-RDA; Legendre and Andersson, 1999), which is the direct extension of multiple regression to the multivariate scale and it allows non-Euclidean dissimilarity indices (i.e., Manhattan distance), was performed with 1,000 permutations and a p-value of 0.001. Throughout db-RDA analysis the ordination process was constrained by the biogeochemical variables and the ordination sought the axes that were best explained by a linear combination of explanatory physical variables (Kindt and Coe, 2005;Borcard et al., 2018). An analysis of variance (ANOVA) was performed to assess the significance of constraints. Both BIO-ENV and db-RDA analysis were realized using the R package Vegan from Oksanen et al. (2019).
For every biogeochemical region the number of physical clusters in which it was divided (and vice versa) were calculated in order to evaluate the degree of agreement between physical and biogeochemical regionalization. The clusters that represented less than 20% of total pixels of the corresponding region were filtered to identify the zones that majorly contributes to its surface. Finally, a non-parametric Mann-Kendall test for significance (with a p-value < 0.05) was used to assess the statistical significance of sea ice trend.

Physical Regionalization
The Southern Ocean presented very heterogenous physical conditions during the biological productive season (December to March). Sea ice coverage, indeed, ranged from 0 to 99%, whereas MLD varied between 12 and 103 m and SST (with a mean value of 2 • C) oscillated between −1 and 12 • C ( Table 2). Physical variables were grouped in an agglomerative hierarchical cluster tree with a cophenetic correlation coefficient of 0.73 (Figure 2A) and 12 significant clusters were identified within the dendrogram ( Figure 2B). Two oceanic ice-free zones north of 60 • S (clusters n • 11 and 12) characterized by deep bathymetry and MLD (Supplementary Table 2) covered approximately the 53% of the study zone (38 and 15%, respectively). On the other hand, regions south of the Polar Front (clusters from n • 1 to 9) were characterized by low SST (mean values ranged from −0.96 to 0.47 • C), shallow MLD (varying between 23.13 to 34.42 m) and reduced PAR values (from 14.34 to 27.87 E m −2 d −1 ). The areas of Falkland Islands and South America were merged within one cluster (n • 10) that presented few pixels with shallow bathymetry, no ice coverage and highest PAR and SST values. Clear circum-Antarctic patterns were detected in several oceanic regions north (n • 8, 11, and 12) and south (n • 2 and 3) of 60 • S, whereas this circumpolar structure was absent for coastal clusters ( Figure 2B).

Biogeochemical Regionalization
The dendrogram of biogeochemical variables ( Figure 2C; cophenetic correlation coefficient: 0.78) was divided into 18 significant clusters. Four oceanic groups of the Atlantic, Indian and Pacific sectors (n • 5, 6, 7, and 11; Figure 2D) were characterized by low POC export (ranging from 4.15 to 5.07 mg C m −2 d −1 ), diminished diatoms concentrations (varying between 0.26 and 0.44 mg Chl-a m −3 ) and covered approximately 52% of the Southern Ocean surface (Supplementary Table 3). In contrary, coastal Antarctic and Southern Patagonian areas (clusters n • 16, 17, and 18) presented elevated NPP values (ranging from 0.44 to 1.98 g C m −2 d −1 ), diatoms concentrations between 0.55 and 1.06 mg Chl-a m −3 , nanophytoplankton concentrations from 37.85 to 45.45 mg C m −3 and POC export between 39.90 and 53.06 mg C m −2 d −1 . Group n • 1 enclosed remote areas from the Weddell Sea (close to the Ronne ice shelf) and the area around the Balleny Islands (near the Ross Sea; Figures 1, 2D), probably because these regions shared the highest dissolved iron concentration in the surface of the Southern Ocean (0.2 nM) and intermediate productivity.
The pixels belonging to the Weddell Sea (cluster n • 17) also presented elevated iron concentration, but the productivity and POC export were higher compared with Balleny Islands and  Ronne ice shelf (cluster n • 1). The circum-Antarctic pattern identified for oceanic physical regions was less clear during the biogeochemical regionalization. Oceanic island such as South Georgia, Auckland and Kerguelen were grouped in a separate region (n • 2) with much higher POC exports (29.03 mg C m −2 d −1 ) with respect to the oceanic areas of the Indian and Pacific Ocean (n • 6, 7, and 8). The physical properties of the South America coastal region were significantly unique for this zone, whereas the biogeochemical regionalization showed that Southern Patagonia shared biogeochemical properties with coastal Antarctic regions of Ross, Bellingshausen, Amundsen and Cooperation Sea (cluster n • 18).

Physical Variables That Explain Biogeochemical Spatial Variability
The latitude showed the highest Pearson linear correlation coefficient (p-value: 0.05) with physical variables (0.77, −0.77, −0.80, and −0.83 for sea ice concentration, SST, MLD, and PAR, respectively; Table 3). On the other hand, the neritic-pelagic gradient seemed to majorly control NPP, diatoms concentration, nanophytoplankton and POC export (with correlation coefficient of 0.55, 0.42, 0.68, and 0.81, respectively), whereas dissolved iron concentration in the photic layer appeared to be mainly controlled by sea-ice dynamics (Pearson linear correlation coefficient: 0.71).
The BIO-ENV analysis identified the combination of bathymetry and sea-ice coverage as the subset of environmental variables with best correlation to biogeochemical spatial variability, with a Spearman rank correlation coefficient of 0.68. Db-RDA (Supplementary Figure 3A) indicated that physical variables expressed the 60.1% of total biogeochemical variance.

Physical Regionalization
Important heterogeneity was found in physical variables within the 48.1 CCAMLR zone, with sea ice coverage varying between 0 and 82%, and SST ranging from −0.7 and 3.4 • C ( Table 2). The agglomerative hierarchical cluster tree for physical variables showed a cophenetic correlation coefficient of 0.78 (Figure 3A), whereas 14 significant clusters were identified. A region with an average sea ice coverage of 3% (n • 5) occupied approximately 33% of the entire 48.1 zone with a southwest (SW)-northeast (NE) distribution that included Shetland Islands ( Figure 3B). This cluster presented intermediate depth (mean value: 1,168 m), SST (0.55 • C) and MLD (31.8 m) values as compared with other groups (Supplementary Table 4). The SW-NE spatial pattern also appeared for warmer free-ice oceanic clusters located toward the Drake Passage (n • 1, 2, and 4), with the regions n • 2 and 5 that seemed to be separated by the Southern Boundary (Figure 1). The zone of the 48.1 zone closer to the Weddell Sea (clusters n • 9 and 10) exhibited the highest sea ice concentrations (from 64.3 to 69.1%) and coolest waters (between −0.68 and −0.61 • C). Marguerite Bay, Gerlache Strait and Crystal Sound (biological hotspots along the Western Antarctic Peninsula; Costa et al., 2007; Figure 1) were clustered into group n • 14, which was characterized by elevated sea ice coverage (mean value: 55%), the shallowest MLD (21 m) and lowest PAR value (20 E m −2 d −1 ) within CCAMLR 48.1 zone. These zones shared physical conditions with the southern part of the Eastern Antarctic Peninsula shelf. Finally, one singleton cluster (n • 6) was observed in the northern tip of Antarctic Peninsula.

Biogeochemical Regionalization
We identified 16 significant clusters within the dendrogram of biogeochemical variables (Figure 3C; cophenetic correlation coefficient: 0.84). Gerlache Strait, the eastern part of Marguerite Bay and one oceanic pixel were identified as singleton biogeochemical clusters (n • 2, 9, and 13) in the CCAMLR zone 48.1 (Figure 3D). Marguerite Bay, the coastal zone of the Eastern Antarctic Peninsula and the area around Crystal Sound (clusters n • 1, 2, 3, and 16) were the most productive clusters in terms of NPP, diatoms, nanophytoplankton concentrations and POC exports (Supplementary Table 5). Conversely, the oceanic region off the Western Antarctic Peninsula (clusters from n • 4 to 10) were the most oligotrophic and presented the least POC export within the 48.1 zone. Pixels belonging to group n • 12 appeared to accurately follow the continental shelf slope and the southern Antarctic Circumpolar Current boundary along the Western Antarctic Peninsula (Moffat and Meredith, 2018). The division marked by the Southern Boundary in physical regionalization appeared less evident in the biogeochemical, whereas a separation between the eastern (n • 14) and western part (n • 15 and 16) of the Bransfield Strait was found ( Figure 3D). The similarity between Marguerite Bay and the southern shelf of the Eastern Antarctic Peninsula physical conditions was also observed in the biogeochemical properties (n • 1, 2, and 3).

Physical Variables That Explain Biogeochemical Spatial Variability
The longitudinal gradient in CCAMLR zone had a stronger influence in sea ice concentration and SST ( Table 4) as compared with the results of the entire Southern Ocean ( Table 3). On the other hand, the latitudinal gradient strongly influenced both physical (especially MLD and PAR, with −0.81 and −0.87, respectively) and biogeochemical properties in 48.1 zone (NPP and nanophytoplankton concentrations). Pearson linear correlation coefficient suggested that bathymetry was the main physical factor affecting biogeochemical properties as iron (0.65), POC export (0.96) and diatoms concentrations (0.77; Table 4).
Indeed, BIO-ENV analysis indicated bathymetry as the physical variable with higher correlation to biogeochemical spatial variability (Spearman rank correlation coefficient: 0.81). Additionally, db-RDA (Supplementary Figure 3B) showed that physical variables explained the 77.1% of biogeochemical variance.

Physical and Biogeochemical Regionalization Correspondence
The pixels of each Southern Ocean biogeochemical region belonged on average to 6 different physical zones (± a 95% confidence interval of ±1.5), whereas each physical region was divided in 9 ± 1.4 biogeochemical zones (Figures 4A,B). However, by filtering the results it was calculated that every biogeochemical (physical) region of Southern Ocean averagely corresponded to 1.8 ± 0.2 different physical (biogeochemical) zones. On the contrary, the correspondence degree between biogeochemical and physical regions was higher for CCAMLR 48.1 zone compared with the entire Southern Ocean, since the pixels of each biogeochemical (physical) region on average belonged to 2.3 ± 0.7 (2.7 ± 0.9) physical (biogeochemical) regions (Figures 4C,D). The results decreased to 1.5 ± 0.3 and 1.7 ± 0.4, respectively, when only the pixels that represented more than 20% of total pixels of the corresponding region were considered.

Comparison Between Physical and Biogeochemical Regionalization
This study allowed the identification of new physical and biogeochemical regions within the Southern Ocean, enabled Italics numbers indicated the results not statistically significant (p-value: 0.05). Lon, longitude; Lat, latitude; Bat, bathymetry; Ice, ice coverage; SST, sea surface temperature; MLD, mixed layer depth; PAR, photosynthetically active radiation; Iron, iron concentration in the surface layer (0-50 m); Dia, diatoms concentration; NPP, net primary production; Nano, nanophytoplankton biomass; POC, particulate organic carbon export.  Italics numbers shown the results not statistically significant (p-value: 0.05). See Table 3 for acronyms definitions. The white (colored) dots represent the clusters that cover less (more) of 20%. Panels (A,C) show the pixels distribution for every physical cluster within the biogeochemical clusters, whereas panels (B,D) indicated the pixels distribution for every biogeochemical cluster within the physical clusters.
the comparison with previous regionalization and provided more explanatory information with respect to other studies. In addition, setting the northern limits at 50 • S enabled the creation of new clusters and the study of similarities and differences between Antarctic and Southern Patagonian environments. Cluster analysis revealed a higher number of significant biogeochemical regions with respect to physical zones.
The number of Southern Ocean regions identified in the present work (12 and 18 for physical and biogeochemical variables, respectively) was comparable with other studies of bioregionalization. Indeed, the study of Grant et al. (2006) divided Southern Ocean in 14 different regions based on bathymetry, SST and nutrients concentration, whereas 20 regions were established by Raymond (2014) using bathymetry, SST and ice coverage. Instead, the study of regionalization for CCAMLR 48.1 zone is the first of our knowledge and identifies a surprisingly high number of significant zones (14 physical and 16 biogeochemical regions), probably due to absence of a priori non-hierarchical grouping. The observed circumpolar distribution of physical oceanic clusters marked by the oceanic fronts pattern was similar with the findings of previous regionalization studies (Grant et al., 2006;Spalding et al., 2012;Raymond, 2014). Additionally, Falkland Islands and the coastal zone of the Southern Patagonian were grouped into one physical cluster, which was in agreement with previous works. This work divided the Antarctic coastal area in seven different clusters, which coincided with Raymond (2014) and differed with the regionalization of Grant et al. (2006) that grouped the entire Antarctic coastal area into one region. Grant et al. (2006) additionally identified a separated cluster for the oceanic area of the Weddell gyre, in contrast with the division of this area into four different regions shown by this research and supported by the results of Raymond (2014). The inclusion of sea ice spatial variability into the regionalization variables (Raymond, 2014; this study) might originate higher physical heterogeneity with respect to a regionalization (Grant et al., 2006) that incorporated nitrate and silicate concentrations, since elevated and homogenous macronutrient concentrations are generally observed south of the Southern Boundary and along the Antarctic coastal area (Pollard et al., 2006;Sigman and Hain, 2012).
As expected, the Subantarctic and the Polar Fronts spatial variability influenced the Southern Ocean physical regionalization generating circum-Antarctic clustering ( Figure 2B). On the contrary, the spatial variability of oceanic fronts was not reflected in biogeochemical regionalization (Figure 2D), since biogeochemical regions were identified crossing the oceanic fronts in the Indian and Pacific Ocean. In fact, the distinction was not reflected in NPP and phytoplanktonic groups distribution despite the clear difference in SST and MLD at both sides of the fronts (Supplementary Figure 1). Low NPP estimations might be explained by a combination of light limitation (Venables and Moore, 2010), zooplankton grazing (Le Quéré et al., 2016) and dissolved iron depletion (Tagliabue et al., 2014;Wadley et al., 2014).
The coastal zones around South America and Falkland Islands were merged into one physical cluster significantly different from all others in the Southern Ocean ( Figure 2B). Contrarily, the coastal region around the southern part of Patagonia was divided in two regions according to biogeochemical variables (i.e., 16 and 18) and the cluster n • 18 also appeared near the coastal zone of Ross, Bellingshausen, Amundsen and Cooperation Sea (Figures 1, 2D). This suggested that coastal zones of Patagonia and some Antarctic zones shared common biogeochemical features (high NPP estimations and diatoms concentrations; Supplementary Table 3), but not physical dynamics. The high productivity of coastal Antarctic zones might be fueled by marginal sea-ice melting during austral summer, injecting micronutrients (Lannuzel et al., 2016) and increasing the photic layer stability (Taylor et al., 2013), whereas the high NPP around the southern part of the Patagonia might be sustained by high tidal energy dissipation and mixing above the continental shelf (Romero et al., 2006), continental freshwater discharges (Cuevas et al., 2019) and the possible iron supply by shelf transport (Garcia et al., 2008) or atmospheric dust (Jickells et al., 2005). Future modeling studies might test the influence of physical and chemical factors on primary production and carbon export in Antarctic and Southern Patagonian related biogeochemical zones under different climate change scenarios.
Additionally, Southern Ocean waters surrounding oceanic islands as Auckland, Kerguelen and South Georgia (Figure 1) were grouped in separated biogeochemical clusters (n • 2 and 13; Figure 2D), whereas this distinction was not revealed by physical regionalization. The elevated NPP rates (Supplementary Table 3) and the annual bloom around these oceanic islands can be explained by an "island mass effect" primary produced by downstream iron fertilization via advection (Boyd et al., 2004;Robinson et al., 2016). Sediment pore waters of oceanic islands are a key source of iron in the Southern Ocean (Graham et al., 2015), since iron is released into the photic layer via internal waves and mixing (Korb et al., 2008) and can fuel phytoplankton blooms. In spite of this, interannual variability in the magnitude and extension of the island mass effect might also be influenced by macronutrients and light limitation (Robinson et al., 2016).
The spatial distribution of the Southern Boundary and the Polar Front influenced the physical regionalization of the NW zone of the CCAMLR 48.1 region (clusters n • 2 and 5) but this division was not clearly detected in the biogeochemical regionalization (Figures 1, 3). On the other hand, the physical regionalization of the CCAMLR 48.1 zone didn't identified a difference between the eastern and western parts of the Bransfield Strait which, on the contrary, was observed in the biogeochemical regionalization (Figure 3). Intense continental meltwater input and the influence of the Antarctic Slope Current (Figure 1) that brings waters from the Weddell Sea (Sangrà et al., 2011;Moffat and Meredith, 2018) might cause the biogeochemical division between the western and eastern part of the Bransfield Strait. Marguerite Bay and the southern shelf of the Eastern Antarctic Peninsula presented similar physical and biogeochemical conditions, with the shallowest MLD and highest NPP estimations and diatoms concentrations, probably due to the elevated presence of sea ice during the austral summer and its influence on biogeochemical cycles (Hyatt et al., 2011;Siegert et al., 2019).
No strong consistency was observed between the physical and biogeochemical regionalization of the Southern Ocean (Figure 4), whereas this correlation was higher for CCAMLR 48.1 zone. The higher consistency between physical and biogeochemical regionalization for 48.1 zone might be explained by either the reduced sample size and/or the response of both biogeochemical and physical regions to the spatial patterns of the shelf slope (Figure 3), as also pointed out by Heywood et al. (2014). The inconsistency detected for the Southern Ocean might be ascribed to the lower spatial variability of biogeochemical variables explained by physical factors as compared with the CCAMLR 48.1 area (60 and 77%, respectively). Indeed, Southern Ocean regions with different physical conditions might present the same biogeochemical features (and High Nutrient Low Chlorophyll patterns), as resulted from oceanic areas around the Subantarctic and the Polar Fronts (Figures 2B-D). This difference might induce a more comparable number of physical and biogeochemical clusters for CCAMLR 48.1 zone (14 and 16, respectively) with respect to the Southern Ocean (12 and 18).
Neither physical nor biogeochemical Southern Ocean regionalization identified a separated cluster for CCAMLR 48.1 zone. Indeed, the limits of the CCAMLR zone included heterogeneous zones with coastal and open water areas, islands, as well as areas belonging to either the Western or Eastern part of the Antarctic Peninsula. This heterogeneity was reflected in elevated number of biogeochemical and physical clusters within CCAMLR 48.1 zone, suggesting that the 48.1 limits were identified exclusively for statistical purposes 11 and should be revised following biogeochemical criteria.

Main Physical Factors Explaining Biogeochemical Spatial Variability
BIO-ENV analysis revealed that the combination of bathymetry and sea-ice concentration majorly controlled the biogeochemical spatial variability in the Southern Ocean. Both the neriticpelagic gradient and sea ice variability affect nutrients and iron cycles and stoichiometry (Garcia et al., 2008;Torres et al., 2014;Annett et al., 2017;Cuevas et al., 2019;Henley et al., 2020), water column stability (Hellmer, 2004;Romero et al., 2006;Cuevas et al., 2019) and irradiance balance (Smith Jr. and Comiso, 2008). In addition, bathymetry directly affects the POC oxygen exposure time to remineralization (Dunne et al., 2007), whereas high pulses of krill fecal pellets that strongly contribute to POC exports (Henley et al., 2020) have been discovered near the marginal ice zone (Belcher et al., 2019). As a result, sea ice variability is a sensible mechanism able to influence NPP rates, POC exports and phytoplankton species composition (Deppeler and Davidson, 2017). Changes in dissolved iron supply, light availability, stratification and ice-free period are projected to affect the biogeochemical cycles and increase Southern Ocean NPP rates up to 50%, in contrast with global trend of declining NPP (Henley et al., 2020). Sea ice variability also influence the production of dimethyl sulfide, a climate-active gas with high concentration in the marginal ice zone Hendry et al., 2018;Henley et al., 2019), and plays an important role in controlling the amount of organic matter that reaches the benthic community (Hendry et al., 2018;Henley et al., 2019). The strong interdependency between sea ice and phytoplankton abundance and size spectrum implies that its variability is also able to alter higher trophic levels of the food web, changing the structure and functioning of the entire ecosystem (Kerr et al., 2018). Indeed, a reduction in sea ice coverage and a shift from a cold and dry climate toward maritime conditions might disadvantage some species (such as Adélie penguins, Weddell seals and minke whales) and favor others (e.g., gentoo 11 https://www.ccamlr.org/en/organisation/explanation-terms penguins and humpback whales; Henley et al., 2019 and references therein).
Bathymetry was identified as the physical variable that majorly described the biogeochemical spatial variability in the CCAMLR 48.1 zone, contrarily to the Southern Ocean where the combination of bathymetry and sea ice was established as the main feature controlling biogeochemical variability. The exclusion of sea ice variability in the 48.1 area might be ascribed to the low spatial heterogeneity in sea ice concentration during the austral summer, since most of the Antarctic Peninsula exhibited low sea ice concentrations with exception of Marguerite Bay and the southern part of the Eastern Antarctic Peninsula (Supplementary Figure 2).
The number and extension of biogeochemical regions might be affected by future variability since both sea ice and bathymetry will experiment changes (on different time scale). Indeed, changes in continental shelf width and bathymetry might occur on decadal to century time scales because of the balance between glacial isostatic adjustment (+41 mm yr −1 ; Barletta et al., 2018) and sea level rise (+1.2 mm yr −1 ; King et al., 2012). The spatial variability of sea ice trends for the Southern Ocean during 40 years (from 1978 to 2018; Figure 5) showed significant heterogenous trends along the major Antarctic marginal sea ice zones and the CCAMLR 48.1 zone, mainly explained by surface atmospheric circulation patterns (Eayrs et al., 2019) and deep ocean warming . Despite the negative sea ice trend detected in the Western Antarctic Peninsula since 1978, a hiatus in atmospheric temperature and sea ice trend has been recorded in this area during the last two decades as part of the natural interdecadal climate variability Henley et al., 2019). Parkinson (2019) shown that sea ice coverage trends were reversed during the 2014-2018 period with sea ice loss in the Weddell and Ross Sea and in the Indian and Western Pacific Ocean; whereas a sea ice gain was observed for the Bellingshausen and Amundsen Sea. These opposite sea ice trends were mainly explained by the combination of atmospheric events (El Niño-Southern Oscillation, the Southern Annular Mode fluctuations and regional circulation flows) and oceanic variability (subsurface heat anomalies and polynyas creation; Parkinson, 2019;Turner et al., 2020). Oceanic biogeochemical regions of the Southern Ocean might also be affected by climate variability, since a poleward expansion of subantarctic waters and a contraction in biogeochemical regions closer to the Antarctic continent were projected for 2100 (Reygondeau et al., 2020).

Future Challenges
This work provides useful information for future parametrization and calibration of biogeochemical models for both the Southern Ocean and the 48.1 zone. Further subjective clustering methods might allow the identification of macro physical and biogeochemical provinces (Spalding et al., 2012) in the Southern Ocean in order to facilitate modeling purposes. The results of the biogeochemical regionalization suggested the division of the 48.1 zone in two provinces (the region above the continental shelf and the deeper ocean) for future ecosystem parametrizations. The study additionally suggests the exclusion of the northerner part of the Eastern Antarctic Peninsula and Marguerite Bay during the parametrization of the 48.1 zone because the physical and biogeochemical dynamics seemed to be strongly affected by either higher presence of sea ice for Marguerite Bay (Hyatt et al., 2011;Siegert et al., 2019; Figure 1) or inflow from Weddell Sea in the Eastern Antarctic Peninsula (Moffat and Meredith, 2018).
This study identified the most productive biogeochemical regions for both the Southern Ocean and CCMAR 48.1 zone. This valuable information will be useful in the development and support of marine protected area proposals, although additional information, such as the spatial variability of higher trophic levels (i.e., zooplankton, zoobenthos, fish, birds, and marine mammals), is needed to create a more solid proposition. The biogeochemical regionalization was initially designed with the inclusion of intermediate trophic levels (i.e., zooplankton abundances). Unfortunately, even the largest Euphausia superba and Salpa thompsoni abundance database for the Southern Ocean (i.e., KRILLBASE; Atkinson et al., 2017) presented extended regions with missing data (Supplementary Figures 4A,B). A continuous updating of krill and salp abundance databases within the Southern Ocean is a must, in order to develop detailed studies regarding food web spatial variability. Abundance data of copepods such as Calanoides acutus and Rhincalanus gigas were even sparser in the Southern Ocean (Supplementary  Figures 4C,D), despite of merging the Coastal and Oceanic Plankton Ecology, Production, and Observation Database 12 , the Southern Ocean Continuous Plankton Record database 13 12 https://www.st.nmfs.noaa.gov/copepod/ 13 https://www.scar.org/science/cpr/ and the database provided by Cornils et al. (2018). The summary and standardization of published copepods data within a new Southern Ocean database is considered an important challenge ahead for the scientific community. A regionalization for 48.1 zone including zooplankton abundances is also needed in the near future in order to develop and improve ecosystem management, such as regulation of fisheries, tourism, contamination and creation of marine and terrestrial protected areas and sanctuaries.
The dataset containing the gridded mean value for physical and biogeochemical variables from December to March within the entire Southern Ocean will be shared with colleagues, consistent with the principle of free-access data and analysis in order to face the actual social and climate crisis.

CONCLUSION
This regionalization study allowed the identification of new regions within the Southern Ocean and comparisons with previous regionalization developed using physical and biogeochemical variables. The number of the physical and biogeochemical Southern Ocean regions identified in this study (12 and 18, respectively) was comparable with previous studies of bioregionalization. However, southern Patagonian coastal areas were merged into a physical cluster significantly different from all others in the Southern Ocean, whereas biogeochemical patterns were shared with coastal Antarctic areas (i.e., Ross, Bellingshausen, Amundsen, and Cooperation Sea). The combination of bathymetry and sea ice coverage majorly explained biogeochemical variability within the Southern Ocean (Spearman rank correlation coefficient: 0.68). Since both sea ice and bathymetry will experiment changes, although on different time scale, our results suggest the number and extension of biogeochemical regions might be affected in the future.
Fourteen physical and 16 biogeochemical significant clusters were identified for 48.1 CCAMLR zone, where bathymetry was the main factor explaining biogeochemical spatial variability (0.81). The results suggested the division of the 48.1 zone in two main provinces separated by the shelf slope during future ecosystem parametrizations, with the exclusion of the northerner part of the Eastern Antarctic Peninsula and Marguerite Bay, because those regions were affected by either higher presence of sea ice (for Marguerite Bay) or inflow from Weddell Sea (in the Eastern Antarctic Peninsula).
The correspondence between physical and biogeochemical regions was higher for CCAMLR 48.1 zone with respect to the entire Southern Ocean, probably ascribed to the lower spatial variability of biogeochemical variables explained by physical factors in the Southern Ocean as compared with the CCAMLR 48.1 area (60 and 77%, respectively). Both physical and biogeochemical Southern Ocean regionalization failed to identify a separated cluster for CCAMLR 48.1 zone, suggesting that the 48.1 limits should be revised following biogeochemical criteria.
This study provides useful data for both Southern Ocean and CCAMLR 48.1 zone management and ecosystem parametrization purposes.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.