Beta Diversity of Antarctic and Sub-Antarctic Benthic Communities Reveals a Major Role of Stochastic Assembly Processes

Community assembly is the result of both, deterministic and stochastic processes. The former encompasses niche-based local-scale mechanisms such as environmental filtering and biotic interactions; the latter includes ecological drift, probabilistic colonisation, and random extinctions. Using standardised sampling protocols, we show that the spatial variation in species composition (beta diversity) of shallow subtidal macrobenthic communities of sub-Antarctic (Strait of Magellan and Yendegaia Fjord [Beagle Channel]) and Antarctic (Fildes Bay [King George Island, West Antarctic Peninsula]) localities reflects a high contribution of stochastic processes to community assembly. Null model analyses indicated that random sampling from species pools of different sizes drove the observed among-locality differences in incidence- and abundance-based beta diversity. We analysed a normalised stochasticity ratio (NST), which delimits between more deterministic (<50%) and more stochastic (>50%) assembly. NST was notably larger than 50%, with mean values of 69.5% (95% CI = 69.2–69.8%), 62.5% (62.1–62.9%), and 72.8% (72.5–73.2%) in Strait of Magellan, Yendegaia Fjord, and Fildes Bay, respectively. Accordingly, environmental factors, such as depth, seawater temperature, salinity, and underwater light penetration, accounted for a small fraction of the spatial variation in community composition across the three localities. In this region, therefore, stochastic processes could have stronger effects on community assembly than deterministic niche-based factors. As anthropogenic biotic homogenisation continues apace, our study can give useful insights into the major ecological processes in Southern Ocean’ coastal marine communities.


INTRODUCTION
Community assembly has been a major topic in fundamental and applied ecological research. In this vein, the analysis of the spatial variation in community composition-beta diversity (Whittaker, 1972)-provides important insights into assembly mechanisms (Chase and Myers, 2011;Legendre and De Caceres, 2013). As a relevant dimension of biodiversity, beta diversity is associated with ecosystem functioning and several aspects of stability, such as resistance, resilience, recovery, and invariability (Mori et al., 2018;Valdivia et al., 2021). Since anthropogenic biodiversity loss affects ecosystems functioning and stability especially at large spatial scales, understanding beta diversity at regional scales is particularly relevant today (Isbell et al., 2017;Gonzalez et al., 2020).
Community composition, i.e., the combination of species incidences and abundances, is controlled by both deterministic and stochastic forces (HilleRisLambers et al., 2012;Vellend, 2016;Thompson et al., 2020). On the one hand, deterministic nichebased processes include environmental filtering when species are adapted to narrow environmental conditions, densityindependent physiological responses to abiotic conditions (e.g., thermal tolerance), and density-dependent biological interactions like competition, predation, and facilitation (MacArthur and Levins, 1967;Keddy, 1992;Belyea and Lancaster, 1999;Somero, 2010). Strong correlations between beta diversity and environmental variation would hint for a primary role of deterministic drivers in community assembly (e.g., Legendre et al., 2009;Menegotto et al., 2019;López-Delgado et al., 2020). On the other hand, stochastic processes encompass random events of reproduction and mortality along with random chance for colonisation (probabilistic dispersal), priority effects, historical contingence, and ecological drift led by random changes in species abundances (MacArthur and Wilson, 1967;MacArthur, 1972;Drake, 1991;Hubbell, 2001;Fukami, 2004;Reijenga et al., 2021). When environmental factors account for small proportion of biological variation, large correlations between beta diversity and geographic distances can be used as evidence for a stronger role of stochastic processes (Gilbert and Lechowicz, 2004;Chase and Myers, 2011). Importantly, a predominant influence of these mechanisms on community assembly can generate beta diversity patterns that are indistinguishable from random chance alone (Ning et al., 2019;Liang et al., 2020).
The relative contributions of deterministic and stochastic processes to community assembly may vary among biogeographic regions. For example, the decrease of beta diversity from tropical to temperate forests has been explained as a result of variations in the relative importance of assembly processes (e.g., Qian and Ricklefs, 2007;Soininen et al., 2007). However, the use of null models revealed that such biogeographic pattern is actually the result of differences in species pool size (gamma diversity) and thus random sampling effects: after correcting for differences in gamma diversity, the apparent decrease in beta diversity with increasing latitude disappears (Kraft et al., , 2012Myers et al., 2012). Therefore, processes that generate differences in species pools could have stronger effects on beta diversity than variations in community assembly mechanisms (Ricklefs, 1987;Hubbell, 2001;Benício et al., 2021).
Antarctic and sub-Antarctic marine subtidal communities provide an opportunity to assess the importance of deterministic and stochastic ecological mechanisms. In these regions, marine subtidal communities represent a unique spot of biodiversity. The spatiotemporal distribution of subtidal communities in the Southern Ocean has been shaped by large-scale geological, paleoclimatic, and dispersal processes (Arntz et al., 2005;De Broyer and Koubbi, 2014). As a consequence of the formation of the Antarctic Circumpolar Current (ACC), Antarctica's 33-million-year isolation and associated dispersal limitations between Antarctic and sub-Antarctic areas have boosted the effect of stochasticity in shaping some local communities (Griffiths and Waller, 2016;Convey and Peck, 2019). Yet, recent evidence suggests that dispersal limitations may not be as strong as previously thought, and local environmental filter could account for differences between Antarctic and sub-Antarctic marine communities (Fraser et al., 2018;Avila et al., 2020). The Southern Ocean, for example, has generally cooled since the opening of the Drake Passage (Peck, 2018). Moreover, glacier melting is intensifying in these regions and elsewhere due to climate change, which can be seen as a selective force in both Antarctic and sub-Antarctic coastal communities (Barnes and Clarke, 2011;Stammerjohn et al., 2012;Cauvy-Fraunie and Dangles, 2019). Thus, to what degree niche and stochastic processes influence the assembly of local communities across Antarctic and sub-Antarctic regions is still an open question.
Here, we analyse beta diversity of subtidal hard-bottom communities dominated by macroinvertebrates and macroalgae in three localities across Antarctic and sub-Antarctic regions (Figure 1): Strait of Magellan (ca., −53.7 • S) and Yendegaia Fjord (Beagle Channel; ca., −54.9 • S) in Chilean South Patagonia, and Fildes Bay (King George Island, ca., −62 • S) in West Antarctic Peninsula. Using identical field protocols, we estimated species incidences and abundances in 50 × 50-cm plots. This plot size was appropriate to capture responses to fine-scale environmental heterogeneity and the outcome of local species interactions (see also Kraft et al., 2011Kraft et al., , 2012Menegotto et al., 2019). Alpha diversity was defined as the number of species in a single 0.25 m 2 plot. Gamma diversity of each locality was estimated as asymptotic species richness calculated as first-order Jacknife, which allowed us to reduce the bias in observed species richness. In this context, observed number of species is assumed to be a biased underestimation of the complete assemblage richness (Gotelli and Colwell, 2011). Following Anderson et al. (2011), beta diversity was estimated for each locality as the amongplot variation in species incidences and abundances (Jaccard and Bray-Curtis dissimilarities, respectively). This framework makes clear distinction between presence/absence vs. relative abundance data, and the exclusion of "double-zeros" or jointabsence data, which usually characterise species abundance datasets. In addition, we did not use multiplicative or additive beta partitions (e.g., Tuomisto, 2010) to avoid the mathematical interdependence between beta, alpha, and gamma diversities . The plots at each of the three localities covered a range of environmental conditions, including depth  Table 1. Colour scales depict maximal sea surface temperatures (monthly averages) obtained between 2010 and 2020. Data were derived from a multi-sensor ultra-high-resolution analysis (https://coastwatch.pfeg.noaa.gov/erddap) and visualised in Ocean Data View (Schlitzer, 2018;http://odv.awi.de).
(between 5 and 20 m), distance from the nearest glacier, and river inputs (Huovinen et al., 2016Valdivia et al., 2020;Palacios et al., 2021). The sampling design allowed us to test two competing hypotheses: (H1) Long-term spatial isolation and differences in abiotic environmental conditions underpin different assembly rules among Strait of Magellan, Yendegaia Fjord, and Fildes Bay. Therefore, among-region differences in beta diversity should hold after correcting for between-locality variation in gamma diversity. Also, the relative contributions of deterministic and stochastic processes to community composition should vary among regions.
(H2) Alternatively, broad-scale factors that generate differences in species pools might drive the observed between-region differences in beta diversity. Consequently, any observed difference in beta diversity among regions should disappear after correcting for gamma diversity. Moreover, stochastic processes will account for the largest proportion of beta diversity across the three localities.

Study Regions
Three localities were analysed (Figure 1): Strait of Magellan, Yendegaia fjord (both in Chilean South Patagonia), and Fildes bay (South Shetland Islands, West Antarctic Peninsula). For brevity, the tree localities are referred to as STRA, YEND, and FILD, respectively, (Figure 1 and Table 1).
STRA is a narrow channel that extends ca. 560 km from the Pacific (ca. −52.6 • S) to the Atlantic (ca. −52.5 • S). The western side of the channel has an annual rainfall between 2,000 and 5,000 mm, which, along with glacier meltdown waters, provide can 6,500 m 3 s −1 freshwater to the channel (Brun et al., 2020). Tidal amplitude (a proxy for tidal currents) in the studied sector of the Strait of Magellan is 126 cm (J. Garcés, unpublished data; Medeiros and Kjerfve, 1988). Mean seawater temperature, salinity, and penetration depth of visible light is 7.3 • C, 28.5 PSU, and 32 m, respectively, (see details in section "Environmental variables, " below). YEND is located in the Beagle Channel, which delimits the southwest border of the Cordillera Darwin Icefield. The fjord receives freshwater inputs from Stoppani glacier through a 12-km-long river coming out onto fjord's head , which generates a gradient of salinity and turbidity Palacios et al., 2021). Tidal amplitude in YEND is 97.6 cm (J. Garcés, unpublished data). In this locality, mean seawater temperature, salinity, and underwater light penetration is 6.6 • C, 31.6 PSU, and 21 m, respectively. FILD is located in King George Island, Shetlands Islands. The bay is 14 km long and 6 to 14 km wide, with a maximal depth of 500 m. In FILD, glacier meltdown affects seawater turbidity, temperature, and salinity, which influences the structure of macrobenthic communities . Tidal amplitude in FILD is 123.2 cm (J. Garcés, unpublished data). FILD exhibits mean seawater temperature, salinity, and light penetration of −0.34 • C, 30 PSU, and 31 m, respectively. Beta diversity of hard-bottom macrobenthic species was estimated at each locality as the spatial variation in community composition among 50 × 50-cm plots. From STRA, we analysed 33 randomly distributed plots in two sites: Punta Santa Ana and Faro San Isidro (18 and 15 plots, respectively). From YEND, the same number of plots were distributed across three sites: fiord's head, middle, and mouth (13, 11, and 9 plots, respectively, see also Huovinen et al., 2020). Similarly, at FILD we analysed 33 plots distributed in four sites: Collins, Artigas, Suffield, and Albatross (8, 9, 8, and 8 plots, respectively, see also Valdivia et al., 2020). Sites spanned ca. 200 m 2 . Between-site distance through waterways was 19.3 km in STRA and ranged between 3.0 and 6.4 in YEND, and between 2.4 and 8.0 km in FILD.

Estimation of Species Abundances
Identical field sampling protocols were used to estimate species abundances at each locality. The study was conducted during August 2016 in STRA, July 2019 in YEND, and between January and February 2017 in FILD. At each locality, plots were haphazardly located within areas of rocky substratum, of similar slopes, and lacking large crevices. Plot slopes ranged between 10 • and 20 • and depths between 5 and 20 m. For each plot, we used suction dredge sampling to collect sessile and mobile macrobenthic organisms (>5 mm length; see also Wahle and Steneck, 1991). We used a portable underwater venturesuction dredge, constructed with an auxiliary SCUBA cylinder connected-through a BCD hose and an air-feed valve-to the wall of a 2.3-m-long, 10-cm-diameter PVC pipe. One of the extremes of the pipe dredged the surface of each plot. The other extreme of the pipe was equipped with a 5-mm pore mesh that received the extracted material. Before sampling, a 0.5 m 2 metallic frame was placed on each plot to standardise the sampled area and to visually estimate the percentage of the plot covered by bare rock and sand. Additionally, each plot was scraped with a spatula to collect seaweeds and invertebrates that remained attached to the rocky substratum. Albeit we originally aimed to rocky substratum, this method has been proven to efficiently sample organisms living associated to boulders, rock, and also the sediments beneath (Wahle and Steneck, 1991). However, the method may underestimate the occurrence and abundance of larger organisms, such as large macroalgae.
All samples were placed into independent, labelled plastic bags on the boat and transported within few hours to the facilities of Centro IDEAL in Punta Arenas (STRA) or Profesor Julio Escudero Research Station of the Instituto Antártico Chileno, INACh (FILD). In YEND, samples were fixed in 70% ethanol and transported to the facilities of the Universidad Austral de Chile in Valdivia (south-central Chile) where organisms were cleaned, sorted, and identified to the lowest taxonomic level possible, usually to species level. Invertebrate and seaweed abundances were expressed as wet weights (g 0.25 m −2 , 0.01 g accuracy).
Before the analyses, taxon data were standardised to proportions of the maximum observed for each taxon across the plots.
Alpha diversity was calculated as the number of species observed in each plot. Gamma diversity was expressed as asymptotic species richness in order to reduce the bias due to underestimation of "true" species richness. To this aim, we used first-order Jacknife, which consists of a non-parametric statistical technique in which subsets of plots are removed from the dataset and the estimator is recalculated (Gotelli and Colwell, 2011).

Beta Diversity and Null Model
A null-model approach was used to compare the observed beta diversity to that expected from random samplings of the regional species pool Kraft et al., 2011). Observed beta diversity was measured as between-plot Jaccard (i.e., incidencebased) and Bray-Curtis (i.e., abundance-based) dissimilarities. Dissimilarities were calculated as: where C ij is either Jaccard or Bray-Curtis resemblance between the ith and jth communities. We then simulated null-model species assemblages in each plot by randomising the assemblage data matrix with fixed regional species richness and speciesoccurrence probability proportional to observed frequencies. The null model was iterated 1,000 times. From the observed and null expected assemblages, we calculated standardised effect sizes as: where beta deviation is standardised effect size, D ij is observed dissimilarity (either Jaccard or Bray-Curtis) between the ith and jth communities, G ij is the average null expected dissimilarity matrices, and V ij is the standard deviation of the null expected dissimilarity . G ij was calculated as: where C ij is the mean null expected resemblance between the ith and jth communities (Ning et al., 2019). A value of beta deviation equal to zero indicates that observed spatial variation in community structure does not differ from random sampling, a positive beta deviation indicates higher beta diversity than expected by chance, and a negative beta deviation indicates lower beta diversity than expected by chance Kraft et al., 2011).
To assess the relative contributions of deterministic and stochastic processes to community assembly, we used Ning et al.
(2019)'s normalised stochasticity ratio (NST). NST is bounded between zero and one, with values below and above 0.5 (50%) representing more deterministic or more stochastic assembly, respectively, (Ning et al., 2019): where NSS is the normalised selection strength and D SS and T SS are extreme values of selection strength, reflecting fully deterministic and stochastic assembly, respectively. D C ij is the similarity between community i and j under extremely deterministic assembly. E ij (k) is a given null expected similarity between community i and j under stochastic assembly. ξ is a generalised function for SS ij under observed, extremely deterministic, or stochastic assembly (Ning et al., 2019). In this way, NST allows us to simultaneously consider the situations when, compared to null expectations, deterministic factors lead to either more similar or more dissimilar observed communities (Ning et al., 2019). NST were recalculated for each of 1,000 random samples of data (bootstrapping). The patterns of incidence-and abundance-based beta diversity were similar across localities (e.g., Supplementary Tables 1, 2). To improve brevity, we report the results of incidence-based beta diversity in the main manuscript, while those of abundancebased beta diversity are fully reported in the supplementary electronic material. The analyses also included a matrix of between-site geographic distances, based on latitude-longitude data of each site.

Environmental Variables
Seawater temperature, salinity, and light penetration were analysed as site-level environmental variables ( Table 2). Seawater temperature influences metabolic rates, phenology, and local population growth rates of marine organisms (Strathmann et al., 2002;Baldanzi et al., 2018;Suarez et al., 2020). Also, seawater temperature correlates with other environmental variables such as nutrient and Chlorophyll-a concentration, which can also have predictable effects on benthic diversity (Witman et al., 2008). Variations in salinity can be associated with glacier meltdown processes and river inputs (e.g., Turner et al., 2017). In addition, light penetration is a useful proxy for seawater turbidity and sedimentation, which influences the physiology of subtidal marine organisms (Sahade et al., 2015;Vause et al., 2019).
In each site at STRA and FILD, we deployed a self-contained thermistor (Star-Oddi DST CT, Garðabaer, Iceland), recording seawater temperature ( • C) and salinity (PSU) every 30 min. Each sensor was encased in a PVC pipe, housed in a concrete block, and deployed by SCUBA divers at ca. 10 m depth from mean low water. At STRA, thermistors were deployed in August 2016 and retrieved in October 2016. At FILD, thermistors were deployed in February 2017 and the data were retrieved every 12 months until January 2019. Due to overgrowth of the sensors' conductivity plates by encrusting organisms, we used only the first 2 months of the salinity timeseries. The temperature records were not affected by biofouling because CT temperature sensors are not as sensitive as conductivity plates. Moreover, the temperatures obtained during the summer months were similar to previously recorded values in the region (Höfer et al., 2019). These data are readily available in Valdivia et al. (2020). At YEND, seawater temperature and salinity were measured in each site with a CTD Seabird 19plus, which was deployed from the surface down to a few metres above the bottom in October 2016. In our study, we used the mean of the records between 5 and 30 m depth. These data are available in Giesecke et al. (2019).
Light penetration was expressed as the depth (m) at which photosynthetic active radiation (PAR) decays to 1% of subsurface conditions (Z 1% ). To that aim, we used a multichannel radiometer (PUV-2500, Biospherical Instruments Inc., San Diego, United States). Underwater light profiles at PAR waveband (400-700 nm) were measured during sunny or partly cloudy days, with calm or moderate wave conditions, and within few hours between 11:00 and 16:00 h to minimize the impact of solar zenith angle on light measurements. Z 1% was calculated as: where Z is depth (m), E d is irradiance, E d0 is the irradiance right below seawater surface, K d is vertical diffuse attenuation coefficient of irradiance, and log is natural logarithm (Huovinen et al., 2016).  et al. (2016) and Palacios et al. (2021). In addition to sitescale factors, we analysed plot depth (m) as potential drivers of community composition.

Statistical Analyses
The differences in beta diversity (observed, expected, and deviation) and NST among localities were tested with permutational linear models. In order to avoid inflating the rate of type-I error, we used a "treatment" contrast, in which the northmost locality, Strait of Magellan (STRA), was compared against the other two localities. The P-value of each parameter was obtained after 5,000 permutations of the data. Beta deviation and NST were compared against zero and 0.5, respectively, with permutational t-tests. The influence of environmental and spatial (latitude and longitude) variables on observed beta diversity and beta deviations was investigated by means of distance-based partial redundance analyses (db-RDA; Borcard et al., 2004;Peres-Neto et al., 2006). These analyses allowed us to partition the variation in beta diversity and deviations into environmental and spatial fractions. db-BDA were computed for both Jaccard and Bray-Curtis dissimilarities. The spatial variables consisted of spatial eigenfunctions computed from Principal Components of Neighbour Matrices (PCNM) of latitude and longitude data of each site (Dray et al., 2012). PCNM from each neighbourhood matrix were selected according to adjusted R 2 from RDA after stepwise model building. Then, the fractions of variation in community composition accounted for by the spatial relationships among sampling sites (selected PCNM) and depth and site-scale environmental factors were estimated as adjusted R 2 from the db-RDA ordinations.
All analyses were conducted in R programming environment (R Core Team, 2021). The R packages tidyverse and cowplot were used to data pre-processing and plotting, vegan to estimate dissimilarity matrices, NST to estimate bootstrapped NSTs, lmPerm for permutational linear models, and MKinfer for permutational t-tests (Wheeler and Torchiano, 2016;Wilke, 2016;Ning et al., 2019;Oksanen et al., 2019;Wickham et al., 2019;Kohl, 2020). (Figure 2). Temperature decreased from Strait of Magellan (STRA) to Yendegaia Fjord (YEND) and to Fildes Bay (FILD). However, temperature exhibited little between-site variation within YEND and FILD (Figure 2A). Salinity showed a similar spatial pattern in YEND, but with larger between-site variation in FILD. Light penetration (Z 1% ) varied at both spatial scales, hinting for spatial gradients in light attenuation and likely turbidity within each locality ( Figure 2C).
Observed beta diversity followed the spatial pattern of alpha and gamma diversity, increasing from STRA to YEND and from YEND to FILD (Figure 4A; predicted means of STRA, YEND, and FILD = 0.76, 0.81, and 0.85, respectively; R 2 = 0.14; Supplementary Table 2). These patterns were expected by chance ( Figure 4B; predicted means of null expected beta diversity at STRA, YEND, and FILD = 0.76, 0.84, and 0.87, respectively; R 2 = 0.35; Supplementary Table 2). Accordingly, beta deviation-the standardised difference between observed beta diversity and that expected by chance-followed a completely different pattern ( Figure 4C). In average, beta deviation in STRA was indistinguishable from zero (mean = 0.03, 95% [Confidence Intervals = −0.12 and 0.16]; Supplementary  Table 3 Normalised stochasticity ratio reached average values above 50% in the three localities (Figure 5 and Supplementary Table 3). Mean NST (95% confidence intervals) were 0.695 (0.692 and 0.698), 0.625 (0.621 and 0.628), and 0.728 (0.725 and 0.731) in STRA, YEND, and FILD, respectively. This indicates that community assembly was more stochastic than deterministic across localities. Incidence-and abundance-based NST showed similar patterns, with larger values for the latter (Supplementary Figure 2 and Supplementary Tables 2, 3).
Across the three localities, the joint influence of depth, site-scale factors, and spatial distances accounted for 7% of observed beta diversity (incidence-based, Jaccard; Figure 6A). The combined effects of site-scale and spatial factors accounted for 43% of observed beta diversity. The joint effects of depth and spatial factors accounted for 2% of observed beta diversity. Each of the "pure" effects of depth, site-scale, and spatial factors explained a maximum of 12% of observed beta diversity ( Figure 6A). Two PCNM were selected as spatial factors in these analyses, and residual variation was 33% ( Figure 6A). After controlling for gamma diversity, the joint effects of depth, sitescale, and spatial factors decayed to 0% (Figure 6B). Similarly, all other combined effects decreased to zero, excepting that of depth and site-scale factors (1%). The pure effects of site-scale and spatial factors accounted for 1 and 7% of beta deviation, respectively. One PCNM was included in the analysis of beta deviation and residual variation reached a 94% (Figure 6B). Abundance-based beta diversity (Bray-Curtis) exhibited a similar partitioning structure to that of incidence-based beta diversity (Supplementary Figure 3), except that no PCNM was selected in the model for beta deviation. In addition, no explanatory matrix fit to the matrix of abundance-based beta deviation (Supplementary Figure 3).

DISCUSSION
Our results suggest that community assembly at both sides of the ACC is driven by differences in pooled species richness rather than local ecological processes. The observed variation in beta diversity simply disappeared after controlling for among-locality differences in species pools (gamma diversity). At each locality, community assembly was more stochastic than deterministic. Indeed, variance-partitioning suggested a weak fit between scaledependent environmental factors and community composition. Compared to environmental filtering and niche differences, therefore, broad-scale and stochastic events of reproduction, mortality, and dispersal likely have a stronger effect on the assembly of these communities.

Processes Generating Differences in Gamma Diversity
The apparent differences in beta diversity across the Strait of Magellan (STRA), Yendegaia Fjord (YEND), and Fildes Bay (FILD) were the result of random sampling from species pools. This finding suggests that processes that influence gamma diversity play a preponderant role in the assembly of these communities, in agreement with previous surveys conducted FIGURE 3 | Species richness observed at the scales of site (A, alpha diversity) and locality (B, gamma diversity estimated as asymptotic Jacknife species richness).
FIGURE 4 | Incidence-based beta diversity (Jaccard's dissimilarity): Observed (A) and null expected (B) beta diversity across the three localities. Beta deviation (C) is beta diversity after controlling for differences in gamma diversity. Similar results were obtained when using abundance-based beta diversity (Bray-Curtis, Supplementary Figure 1). on marine macrobenthic communities at broader spatial scales (Witman et al., 2004). What processes may contribute to differences in gamma diversity? Density-independent responses to abiotic environmental conditions, density-dependent biotic interactions, and dispersal are processes operating at both local and regional scales (Vellend, 2016;Thompson et al., 2020). Among density-independent factors, substrate composition may be a relevant factor influencing diversity in our study: several burrowing species, for instance, were observed almost exclusively in YEND (e.g., Golfingia margaritacea and Acesta patagonica; Supplementary Table 1), which may hint to differences in substratum composition among the localities. Contrarily, the abundance of benthic communities has been shown to be decorrelated from depth and sediment type in the nearby of STRA and fjords and channels off the South Patagonian Icefield (Ríos and Mutschke, 1999). Thus, these results collectively suggest that substrate may have underpinned compositional differences among localities in our study.
Substrate composition and variability can be strongly influenced by glacier meltdown processes in high-latitude ecosystems (Gutt et al., 1999;Cauvy-Fraunie and Dangles, 2019). For example, accelerated glacier retreat can result in more ice-free rocky substrate exposed to colonisation of benthic macroalgae in Antarctica (Jerosch et al., 2019). In addition, glacier meltdown influences morpho-functional attributes of habitat-forming seaweeds in YEND (Palacios et al., 2021), which could influence understorey species and others closely located. Indeed, YEND exhibited the shallowest light penetration (Z 1% ) in our study and a clear spatial gradient from the head to the mouth of the fjord, which can be attributed to glacier meltdownrelated turbidity. Enhanced turbidity and sedimentation in the nearby of glaciers has strong effects on diversity and composition of subtidal macrobenthic communities in FILD . Depending on substrate slopes and depth (Cardenas and Montiel, 2015), increased sedimentation could well cover denuded rocky substrata at sites located in the nearby FIGURE 5 | Incidence-based normalised stochasticity ratio (NST): The index is bounded between zero and one, depicting fully deterministic or stochastic community assembly, respectively. The 50 % boundary between both extreme scenarios is shown. Similar results were obtained when using abundance-based beta diversity (Bray-Curtis, Supplementary Figure 2). of glaciers, reducing the availability of substrates for settlement (e.g., Khim et al., 2007) and also compromising the fitness of suspension feeders due to appendage clogging and excess of mineral suspension (Husmann et al., 2012;Krzeminska and Kuklinski, 2018;Gutt et al., 2019;Topçu et al., 2019).
In addition to environmental conditions at the locality scale, it is necessary to discuss broad-scale processes that may underpin the structure of these communities. The ACC and Polar Front Zone have isolated Antarctica from other large land masses since ca. 30 million years ago (Barnes et al., 2006). This has limited the long-distance exchange of biota between Antarctica and other regions, which-in addition to a long-term invariability in abiotic environmental conditions-underlies part of a high marine species diversity in Antarctica and probably the differences in gamma diversity observed in our study (Barnes et al., 2006;Poulin et al., 2014;Peck, 2018;Gómez and Huovinen, 2020).
Indeed, long-term isolation likely underpins part of the large proportion of invertebrate species with benthic and non-feeding early developmental stages in Antarctica (Poulin et al., 2002). However, recent evidence indicates that long-distance dispersal of drifting seaweeds across the ACC could well have been taking place after the last glacial maximum (Guillemin et al., 2020). Moreover, anthropogenic global change-related factor, such as increased storminess and human-mediated transport, are boosting the movement of marine organisms toward Antarctic coasts (Chown et al., 2015;Duffy et al., 2017;Fraser et al., 2018;McCarthy et al., 2019;Cardenas et al., 2020). Anthropogenic global change has, therefore, the potential to enhance biotic homogenisation across the Southern Ocean owing the establishment of new species in Antarctic (Olden et al., 2004).
Related to long-distance dispersal limitations, differences in speciation rates and early colonisation times can generate large differences in gamma diversity, even between adjacent biogeographic regions (e.g., Benício et al., 2021). In the Southern Ocean, Milankovitch cycles and concomitant climate change over the last 15 million years have strongly influenced speciation and extinction rates (e.g., Mittelbach et al., 2007;Poulin et al., 2014;Crampton et al., 2016). For instance, the drastic impact of the last glacial maximum triggered recent (ca. 18,000 years ago) macroalgal recolonisation of Antarctic coasts, apparently from open-water adjacent areas ("polynyas"; Thatje et al., 2008;Guillemin et al., 2020). In this line, low population numbers and enhanced generic drift might have contributed to rapid speciation in Antarctica, as suggested for low-latitude regions (Fedorov, 1966;Mittelbach et al., 2007). This, in turn, may explain in part the large observed differences in pooled species richness between the Antarctic and sub-Antarctic regions (Arntz, 1999).

Stochastic Community Assembly Across Localities
The analysis of the normalised stochastic ratio suggested that stochastic processes were more important than deterministic mechanisms in the community assembly of the three analysed FIGURE 6 | Incidence-based multivariate variance partitioning analyses after db-RDA. The analyses were computed for observed dissimilarities (A) and after controlling for differences in gamma diversity (B). Numbers in the Venn diagrams are adjusted R 2 that represent the relative contribution of plot-, site-scale, and spatial factors to observed and corrected dissimilarities. See also Supplementary Figure 3. localities. Correspondingly, physical factors accounted for a modest proportion of community structure after accounting for the among-region differences in gamma diversity. How do we explain the relatively high contribution of stochastic processes, such as priority effects, historical contingency, and ecological drift, to community assembly within the analysed localities? Dispersal limitations, for instance, can be associated with stochastic processes (Lowe and McPeek, 2014;Vellend et al., 2014). Indeed, dispersal limitations can reduce population abundances and thus enhance the effects of demographic stochasticity-random events of birth, reproduction, and deathon community-level properties (Leibold et al., 2004).
On top of that, ecological drift (i.e., random fluctuation in species abundances led by demographic stochasticity) can have stronger effects on assembly in communities with few organisms and when local diversity is very reduced relative to the regional species pool (Shurin and Srivastava, 2005;Myers et al., 2015;Vellend, 2016). Further, recruitment of rare species can be limited by demographic stochasticity and Allee effects; i.e., reduction in organismal fitness at low population size or density (Lande, 1993). In our study, asymptotic gamma diversity was about ninefold alpha diversity across the three localities. Moreover, gamma diversity was comparatively larger than observed number of species in each locality. Therefore, these communities are likely unsaturated and composed by many rare species. Certainly, the prevalence of many species that are rare or singleton indicates that local deterministic processes cannot always account for the spatial variation in biodiversity (Chase and Myers, 2011;Partel et al., 2011).

Caveats and Future Directions
The comparatively high contribution of stochasticity to community assembly would have been the result a low environmental variability within each locality. However, the observed variation in salinity and light penetration have been previously associated with significant responses of the physiology, morphology, and diversity of benthic species in YEND and FILD Palacios et al., 2021). Previous research suggests that habitat heterogeneity determines benthic community assembly in polar regions (Wlodarska-Kowalczuk et al., 2007;Somerfield et al., 2009;Pabis et al., 2015). Albeit those contributions studied soft-bottom communities, they suggest that the assembly of polar benthic communities is actually independent of historical context and alpha and gamma diversities (Pabis et al., 2015). This apparent discrepancy between our and previous results probably stems from the fact that our sampling method, vacuum dredge, allowed us to obtain organisms occurring on both, rocks and sand below the rocks. As a consequence, we were able to identify species that use rocky bottoms (e.g., Aulacomya atra, and several seaweeds), spaces between rocks (e.g., Brachidontes sp.), and species that use sandy bottoms (e.g., G. margaritacea).
On the other hand, our study represents a snapshot of community composition across an extensive region. However, the balance between deterministic and stochastic forces can well vary over time. Rapid warming and increased connectivity in the Southern Ocean could well reduce the effects of both abiotic environmental filters and dispersal limitations on community structure in the future (Cardenas et al., 2020). In terrestrial ecosystems, assembly have been shown to tilt over time toward either stochasticity or determinism, depending on succession, environmental stressors, and temporal covariation in species abundances (Zhou et al., 2014;Ning et al., 2019;Liang et al., 2020;. Contrasting with terrestrial communities, the high prevalence of organisms with seawaterborne propagules in marine communities facilitates the dispersal from regional species pools (Roughgarden et al., 1988;Vermeij and Grosberg, 2010)-this leads to stronger regional effects on marine than terrestrial ecosystems (Cornell and Harrison, 2013). Thus, the effects of climate change on dispersal could have more severe consequences for marine than terrestrial ecosystems (see also Kinlan and Gaines, 2003). Albeit our study did not allow us to detect temporal variability in the deterministic-stochastic balance, our results are well in line with previous studies of marine and terrestrial ecosystems and broader spatial scales (Witman et al., 2004;Kraft et al., 2011Kraft et al., , 2012. The influence of ecological drift on community assembly can be exacerbated in communities with small numbers of individuals (i.e., small community size; Hubbell, 2001;Vellend, 2016). This makes the probability of species extinction to increase as community size decreases (Vellend, 2016). Accordingly, anthropogenic impacts that reduce the number of individuals across species can make the local community more prone to extinctions (Myers et al., 2015). In line to our results, future biodiversity research in the Southern Ocean should target at the effects of community size on assembly processes (Vellend et al., 2014), along with long-term monitoring programmes of species abundances and composition.

CONCLUSION
In summary, differences in species pools seemed to be a major determinant of the observed differences in community composition between and within West Antarctic Peninsula and Chilean south Patagonia. Abiotic environmental factors explained a small fraction of the variation in community composition across the localities. Compared with niche processes, stochastic processes may have a stronger effect on the assembly of these communities. Therefore, conservation policies should aim to protect critical community sizes in order to limit stochastic extinction risk in the Southern Ocean.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
NV: conceptualisation, software, validation, formal analysis, writing -original draft, and visualisation. JG-V: software, visualisation, and writing -review and editing. IgG: investigation and writing -review and editing. IvG, PH, NN, and EM: writing -review and editing. LP: investigation, data curation, writing -review and editing, and supervision. All authors contributed to the article and approved the submitted version.

FUNDING
This study was financially supported by FONDAP grant #15150003, Centro de Investigación en Dinámica de Ecosistemas Marinos de Altas Latitudes (IDEAL). IgG was financially supported by doctoral scholarship #72160109 ("Becas Chile") granted by the Agencia Nacional de investigación y Desarrollo (ANID). While writing, NV was supported by FONDECYT grants #1190529 and #1181300 and IvG and PH by FONDECYT grant #1201069. IM and two reviewers provided useful comments that improved a previous version of this manuscript.