Seasonal development of a tidal mixing front drives shifts in community structure and diversity of bacterioplankton

Bacterioplankton underpin biogeochemical cycles and an improved understanding of the patterns and drivers of variability in their distribution is needed to determine their wider functioning and importance. Sharp environmental gradients and dispersal barriers associated with ocean fronts are emerging as key determinants of bacterioplankton biodiversity patterns. We examined how the development of the Celtic Sea Front (CF), a tidal mixing front on the Northwest European Shelf affects bacterioplankton communities. We performed 16S‐rRNA metabarcoding on 60 seawater samples collected from three depths (surface, 20 m and seafloor), across two research cruises (May and September 2018), encompassing the intra‐annual range of the CF intensity. Communities above the thermocline of stratified frontal waters were clearly differentiated and less diverse than those below the thermocline and communities in the well‐mixed waters of the Irish Sea. This effect was much more pronounced in September, when the CF was at its peak intensity. The stratified zone likely represents a stressful environment for bacterioplankton due to a combination of high temperatures and low nutrients, which fewer taxa can tolerate. Much of the observed variation was driven by Synechococcus spp. (cyanobacteria), which were more abundant within the stratified zone and are known to thrive in warm oligotrophic waters. Synechococcus spp. are key contributors to global primary productivity and carbon cycling and, as such, variability driven by the CF is likely to influence regional biogeochemical processes. However, further studies are required to explicitly link shifts in community structure to function and quantify their wider importance to pelagic ecosystems.


| INTRODUC TI ON
Bacterioplankton are the largest component of biomass in the ocean where they underpin marine food webs and influence biogeochemical cycles (Arrigo, 2005;Falkowski et al., 2008;Meng et al., 2022;Passow & Carlson, 2012;Pomeroy et al., 2007;Teeling et al., 2012;Zehr & Kudela, 2011). Therefore, understanding drivers of spatial and temporal variability in their distribution and abundance is a critical step in understanding wider biogeochemical processes (Bunse & Pinhassi, 2017;Fuhrman et al., 2015). Due to the interconnected, highly dynamic nature of the marine environment, and the prevalence of large-scale ocean currents, shifts in bacterioplankton communities were traditionally assumed to occur over large spatial scales (Hewson et al., 2006).
However, recent studies have demonstrated the importance of smallerscale localized oceanographic features in shaping ecological patterns.
Frontal zones (i.e. physical interfaces between distinct water bodies) represent some of the worlds' most well-studied oceanographic features (Belkin et al., 2009;Sournia, 1994). Their formation is influenced by a range of processes, which in turn mediate their biophysical characteristics and ecological effects (Olson et al., 1994). They are generally characterized by elevated primary and secondary productivity, driven by processes such as nutrient retention (Franks, 1992a;Traganza et al., 1987) and entrainment of zooplankton (Franks, 1992b;Le Fevre, 1987), which permeates through the food web (Polovina et al., 2001;Waggitt et al., 2018). Our understanding of processes occurring at the microbial level is comparatively lacking, but emerging research has shown that permanent shelf-break fronts, caused by differences in pressure between offshore and onshore waters, are clear horizontal "ecotones" for bacterioplankton communities (Allen et al., 2020;Baltar et al., 2010;Baltar et al., 2016;Djurhuus et al., 2017;Hernando-Morales et al., 2017). Current understanding is particularly poor for coastal frontal zones that exhibit greater spatiotemporal variability (but see Lemonnier et al., 2020). This is important as successional patterns of bacterioplankton communities have the ability to regulate the temporal variability in biogeochemical processes (Bunse & Pinhassi, 2017).
Tidal mixing fronts develop at the boundary between mixed and stratified waters and are an interplay between depth, tidal mixing and solar radiation (Pingree & Griffiths, 1978). In deep waters during warmer months, solar radiation overpowers tidal forces resulting in surface stratification, whereas in shallower waters tidal forces are often strong enough to maintain a mixed water column. The Celtic Sea Front (CF) is a tidal mixing front located at the transition between Atlantic waters of the European Shelf and coastal waters of the Irish Sea (Simpson & Hunter, 1974). It is a typical tidal mixing front that develops in spring, intensifies over summer and dissipates in early autumn, with profound effects on coastal marine ecosystems across the region (Scales et al., 2014;Waggitt et al., 2018).
Here, we examined how the development of the CF drives changes in bacterioplankton communities during periods of minimum and maximum front intensity. In doing so, we aimed to comprehensively assess how a tidal mixing front affects bacterioplankton community structure and diversity both spatially (horizontally and vertically) and temporally through a period of front intensification.

| Sampling
Sampling of the CF area was conducted aboard the RV Prince Madog on two independent research cruises in 2018. Cruises took place in May when the CF first became established and again in September when it reached its maximum intensity ( Figure 1). The CF is dynamic, moves with the tide and is affected by weather, meaning its position and vertical thermal profile are difficult to determine exactly, a priori. Therefore, we identified its position in situ where a rapid (within 100 m) increase of >1°C Sea Surface Temperature (SST) was observed, using the ships onboard thermal profiler. Five stations were sampled on either side (one hour steam /~ 10 nautical miles) of this transition, and seawater was sampled from the surface (S), 20 metres (M) and seafloor (B). The stratification of each station was then determined using the difference in temperature between surface (S) and bottom (B) seawater. Where this difference was <1°C, stations were categorized as mixed (Celtic Front Mixed -CFM) and >1°C they were categorized as Stratified (Celtic Front Stratified -CFS) ( Figure S1). In May, this resulted in five CFM and five CFS stations and in September there were seven CFS stations and three CFM stations (Figure 1).
At each station, a conductivity-temperature-depth profile (CTD; Seabird, SBE 911) was taken and 10 L of water sampled at S, M and B depths using a rosette with Niskin bottles. A number of oceanographic parameters were measured with each CTD profile. These were oxygen (mg/L), salinity (s/m), transmissometer beam attenuation (proxy for total suspended matter), chlorophyll a (chl-a) concentration (mg/m 3 ) (via chl-a fluorescence), and irradiance (PAR -Photosynthetic Active Radiation). 500 mL of seawater was concentrated by filtering through a 0.22 μm nitrocellulose filter. Filters were then stored at −20°C until DNA extraction.

| DNA extraction and PCR amplification
DNA extraction and PCR amplification were both performed in a dedicated eDNA laboratory at the University of Salford. Genomic DNA was extracted from whole filters using Powersoil kits (Qiagen).
DNA was quantified using a Qubit (Invitrogen) and re-suspended to 2 ng/μL following the manufacturer's instructions. We amplified DNA for the 16S rRNA V4 region using dual indexed forward and reverse primers according to (Kozich et al., 2013) and (Griffiths et al., 2018).
A blank PCR control was included on each PCR plate as a negative control. PCRs were run in 25 μL volumes using Solis BioDyne 5x HOT FIREPol® Blend Master Mix, 2 μM primers and 1 μL of sample DNA.
Cycling parameters were 15 min at 95°C, followed by 28 cycles of 95°C for 20 s, 50°C for 60 s and 72°C for 60 s and a final extension at 72°C for 10 min. All PCRs were run in duplicate. Duplicates were combined into a single PCR and cleaned using HighPrep™ PCR clean up beads (MagBio). Products were quality checked using an Agilent 2200 TapeStation and quantified using a Qubit (Invitrogen). Samples were pooled in equimolar concentrations in order to minimize PCR and sequencing bias. Paired-end (2 × 250 bp) amplicon sequencing was conducted using V2 chemistry on an Illumina MiSeq platform at the University of Salford, UK.

| Sequence processing
Sequence processing and analysis was conducted in R. Paired-end reads were processed according to the BIOCONDUCTOR workflow for microbiome data analysis (Callahan et al., 2016). Sequences were trimmed and truncated using the "filterAndTrim" function in DADA2 (truncLen, f = 240, r = 160; truncQ = 2; andtrimLeft, f = 20, r = 19), to remove primers and low-quality reads. Amplicon Sequence Variants (ASVs) were resolved using DADA2. Chimeras were removed using the "removeBimeraDenovo" function in DADA2. Sequence taxonomy was assigned using the RDPnaïve Bayesian classifier against the SILVA release 132 database using the "assignTaxonomy" function in DADA2. Sequence read counts, taxonomic assignments and metadata were assembled as a phyloseq object using the R package ""PHYLOSEQ" (McMurdie & Holmes, 2013) and was used in downstream analysis. Samples containing <1000 reads and taxa contributing <0.01% of the reads in the dataset were removed.
ASVs identified from PCR blanks along with those identified as mitochondria, chloroplasts or Archaea were also removed. Rarefaction curves of the processed reads were saturated, indicating good coverage of bacterial diversity ( Figure S2). Sequences are accessible through the EMBL database (accession no. PRJEB61283). ASV table and metadata are available at (https://figsh are.com/s/94f8a 74268 4070a 1e7a4).

| Statistical analysis
The whole dataset was split into two, based on sampling month (May vs September) and all analysis conducted on these separate F I G U R E 1 Map of sea surface temperature (a + b) and stratification (surface minus bottom temperature) (c + d) for May (left) and September (right) 2018, in the Irish Sea. sampling periods. Alpha diversity for each sample was estimated through the Chao1 index (Chao, 1984) implemented through the "estimate_richness" function in the R package "PHYLOSEQ" (McMurdie & Holmes, 2013). The Chao1 index estimates ASV richness, and the standard error surrounding this estimate, based on the observed number of ASVs, the observed number of ASVs occurring only once, and the observed number of ASVs occurring only twice (Chao, 1984). To account for differences in sequence depth between samples in alpha diversity estimates, the dataset was rarefied to the minimum sample depth (3500 reads), using the "rarefy_even_depth" function in PHYLOSEQ. Alpha diversity and environmental parameters were each compared using a two-way Analysis of Variance Multivariate community analysis was conducted on a relative abundance dataset. Whole pattern differences in community structure were visualized through non-Metric Multidimensional Scaling (nMDS) using the "metaMDS" function in the R package "VEGAN".
Permutational Multivariate Analysis of Variance (PERMANOVA) was used to test differences in community structure between Depth categories and Station Stratification and followed the same design as alpha diversity. PERMANOVA was implemented by the function "adonis2" in VEGAN. A Similarity of Percentage (SIMPER) procedure was conducted in VEGAN to determine which taxa contributed the most to observed dissimilarities.
To identify the environmental factors that drive bacterioplankton community structure, a distance-based redundancy analysis was performed using the function "capscale" in VEGAN. Here, models were constructed through forward selection using the "ordiR2step" function in VEGAN based on z-score transformed environmental variables. Multicollinearity between environmental variables was examined by the Variation Inflation Factor (VIF) using the function "vif.cca" (O'brien, 2007). The significance of each constraint was confirmed with ANOVA for dbRDA using the function "anova.cca". An ANOVA-like permutational test (function "permutest") for dbRDA was used to assess the significance of the full model. Significant environmental variables were plotted on the dbRDA ordination using the "envfit" function.

| Temperature
The CF was easily identified through its thermal profile in both May and September. In both months, surface waters (S) were ~1°C greater at CFS stations (i.e., in stratified waters) compared to the CFM stations (i.e., in mixed waters) ( Figure S3). In both May and September, we observed a significant interaction between Station Stratification and Depth (Table S1). In May, there was no temperature difference between surface (S) and seafloor (B) waters at CFM stations compared to ~1.5°C difference at CFS stations. In September, this pattern intensified and surface waters were 4.3°C warmer than the seafloor (B) at CFS stations, whereas the CFM stations the mean vertical temperature difference was 0.4°C ( Figures S1 and S3). The depth of the thermocline shifted between sampling months. This meant that our 20 m (M) sampling depth was below the thermocline (<tc) in May but above the thermocline (>tc) in September.

| Oxygen
In May, there was no significant effect of Station Stratification or Depth on oxygen levels (Table S1). At CFS stations, oxygen levels ranged from 2.42 ± 0.13 mg/L at the surface to 2.26 ± 0.025 13 mg/L at M and B depths. At CFM stations, oxygen levels were generally lower and ranged from 2.30 ± 0.1 at the surface to 2.44 ± 0.05 at M and B depths ( Figure S3). In September, we observed a significant interaction between Station Stratification and Depth. Oxygen levels were 2.26 ± 0.02 mg/L above the thermocline (S and M) at CFS stations and dropped to 1.94 ± 0.17 at the seafloor (B). A less dramatic decrease in oxygen levels was observed at CFM stations where S and M samples were 2.17 ± 0.006 compared to 2.11 ± 0.008 at the seafloor.

| Salinity
In May, we observed a significant interaction between Station Stratification and Depth (Table S1). Salinity above the thermocline of CFS stations (S) was 34.54 ± 0.06 and was lower than all depths at CFM stations. At CFS stations salinity increased to 34.85 ± 0.035 at M and seafloor (B) depths whereas no increase with depth was observed at CFM stations, where salinity was 34.7 ± 0.04. In September, we observed significant effects of both Station Stratification and Depth but no interaction between the two (Table S1). Salinity was greater at CFS compared to CFM stations and significantly higher at the seafloor than at S and M depths. At CFS stations, salinity was 35.12 ± 0.006 at the seafloor compared to 35.06 ± 0.001 at S and M depths. At CFM stations, salinity at the seafloor was 35.03 ± 0.007 at the seafloor compared to 35.07 ± 0.001 at S and M depths.
Fluorescence was variable, with no differences between Station Stratification and Depth (Table S1). Irradiance significantly decreased with depth in both May and September (Table S1). In May, irradiance fell to zero at 20 m (M) and seafloor depths (B). In September, this was only seen on the seafloor (B). Beam attenuation was variable but generally consistent between months and transects. However, it was considerably lower at the seafloor (B) at CFS stations in September.

| Overall biodiversity patterns
In total, we analysed 60 seawater samples, which resulted in

| Alpha diversity
In May, bacteria ASV richness (Chao1 index) did not differ between CFM and CFS stations, but there was a significant Depth effect (Table 1). Here, the surface communities (S) were significantly lower in diversity (107.6 ± 10.4) than those at the seafloor (B) (155.4 ± 12.7) ( Figure 3). Similarly, in September, there was no significant difference between CFM and CFS stations, but there was a Depth effect with the surface communities exhibiting significantly lower diversity (143 ± 13.6) compared to those at the seafloor (225.8 ± 24.8) (Figure 3). There was a clear difference in richness between communities that were above the thermocline (>tc) at CFS stations. In May, richness in communities above the thermocline (CFS-S) was 87 ± 11.3 compared to 138 ± 6.7 at all other stations. In September, these differences were starker, with a richness of 140 ± 8 above the thermocline (CFS-S + CFS-B), compared to 231 ± 17.9 at all other stations (CFS-B and CFM) (Figure 3).

| Community structure (Beta diversity)
Initial comparisons between sampling months showed bacterioplankton communities to be clearly differentiated (PERMANOVA,  (Table S2; Figure S8). The  nMDS ordination suggested communities above the thermocline (CFS-S) are distinctive and do not overlap with other communities in multivariate space (Figure 4).
In September, there was a Transect × Depth interaction indicating community structure shifted with depth at either CFM or CFS stations ( Table 1). This pattern was consistent when the data were analysed at higher taxonomic ranks (Table S2; Figure S8). Post hoc analysis showed communities differed by depth at CFS stations. Here, communities above thermocline (S and M) were distinct from those below it (B) and all depths at CFM stations (Table S3) was similar between depths. SIMPER analysis revealed ASV3 from the genus Synechococcus drove most of this variation between groups. For example, 31% of reads were from ASV3 in CFS-S compared to 17% in CFM-S (Table S4).
We used dbRDA to model the relationship between environmental variables and bacterial community structure, following forward selection. In both months, global tests were significant and forward selection was performed. For May, the final model consisted of temperature and salinity and explained 28% of overall variation. For September, the model consisted of temperature and oxygen and explained 62% of total variation in community composition (Table S5; Figure 4). However, VIF of both (21.5 and 18.9 respectively) indicated there was a high degree of multicollinearity between these two variables.

| DISCUSS ION
Our study demonstrates the clear role a tidal mixing front has in driving shifts in the structure and diversity of bacterioplankton Temperature is the key determinant of bacterioplankton community structure at a global scale (Sunagawa et al., 2015). In our dataset a large proportion of community variation could be explained by a combination of temperature and oxygen, which was predominantly driven by greater temperature variability of CFS-B communities.
However, temperature co-varies with a number of other environmental drivers in frontal systems, which makes its specific role here difficult to conclude. Perhaps most profound is the overlying impact of seawater stratification on phytoplankton succession. In stratified waters, phytoplankton are locked above the thermocline and rapidly proliferate due to elevated temperatures and increased light availability. This elevated productivity depletes nutrients resulting in phytoplankton mortality, decay and a high transfer of organic matter below the thermocline (Bunse & Pinhassi, 2017). The subsequent high-variability in organic matter availability below the thermocline combined with a more thermally variable environment, and less wind mixing at depth (Wu et al., 2018) may also be responsible for greater variability between our CFS-B communities.
A notable proportion of community variation remained unexplained and the underlying drivers horizontally between CFS and CFM communities were less clear. Thermal metabolic scaling and easier diffusion of nutrients in smaller cells result in a shift in phytoplankton size toward smaller picoplankton, especially cyanobacteria, in warm oligotrophic stratified waters (Alderkamp et al., 2006;Bunse & Pinhassi, 2017;Cadier et al., 2017;Flombaum et al., 2013;Morán et al., 2010;Partensky & Garczarek, 2010). Indeed, communities above the thermocline, at CFS stations, were consistently dominated by cyanobacteria, from the genera Synechococcus, and drove a large proportion of observed differentiation between depths and stations. This shift in phytoplankton community structure can also shift that of the underlying bacterioplankton (Camarena-Gómez et al., 2018), largely due to differing quality and bioavailability of Dissolved Organic Matter (DOM) produced by different phytoplankton (Bolaños et al., 2021;Mühlenbruch et al., 2018). Whilst we did not directly measure nutrient profiles, organic matter or phytoplankton communities at the CF, such successional dynamics have recently been shown to be a fundamental driver of community structure of bacterioplankton at the Ushant Front, another tidal mixing front in the Iroise Sea (Lemonnier et al., 2020). Given (i) the thermocline of the wider Celtic Sea is known to act as a boundary between nutrient depleted surface waters and nutrient rich bottom waters (Pingree et al., 1977) and (ii) phytoplankton communities at the CF are similarly structured between stratified and mixed areas across the CF in a similar manor to the Ushant Front (Birrien et al., 1991;Grepma, 1988;Lemonnier et al., 2020;Pemberton et al., 2004;Videau, 1987), it is likely such dynamics are also fundamental at the CF. Moreover, the soft dispersal barrier imposed between stratified and mixed waters may result in differences arising from ecological drift (Hanson et al., 2012). Clearly, future studies encompassing a greater range of environmental drivers, overlying phytoplankton dynamics and examining underlying community assembly mechanisms are needed to fully understand the underlying drivers of bacterioplankton succession.
We observed a considerable decrease (~half) in alpha diversity, in the communities above the thermocline of CFS stations, when compared to surface waters of CFM stations. This is in contrast to shelf-break fronts, where alpha diversity has been to be similar, or increase, across frontal waters (Allen et al., 2020;Morales et al., 2019).
Differences in the processes underpinning different frontal development likely account for this, with shelf-break fronts formed from pressure gradients and often associate with upwelling rather than the balance between tidal mixing and atmospheric heating in tidal mixing fronts. This heating combined with a lack of replenishment of nutrients likely create a stressful environment or need for specialism. This in turn may act as a selective filter, limiting the number of species that can tolerate local conditions, or create intense competition for limited resources (Chase, 2007;Menge et al., 2002).
Moreover, in September, conditions above the thermocline were characterized by elevated dissolved oxygen and a number of studies have observed a negative relationship between oxygen concentration and alpha diversity (Spietz et al., 2015;Walsh et al., 2016;Wang et al., 2015) due to a dominance of a few opportunistic taxa in highly productive areas (Salter et al., 2015;Wemheuer et al., 2014). Indeed, communities above the thermocline were dominated by a few (~30% of sample abundance) opportunistic Synechococcus sp.
The structure of marine bacterioplankton communities are closely coupled to their function (Galand et al., 2018). Therefore, the differences between water bodies observed here, and the mechanisms driving these shifts, may have implications for wider ecosystem processes and biogeochemical cycles. In our study, the overrepresentation of the genus Synechococcus within stratified surface waters may be important as this genus is a large contributor to primary production in marine surface water (Li, 1994;Richardson & Jackson, 2007). Whilst our findings are a critical first step toward linking pattern to regional process, future studies need to move beyond description and explicitly utilize a range of techniques to characterize function. For example, the quantification of functional genes in the bacterioplankton of mangroves has identified key roles of these communities in carbon degradation (Meng et al., 2022). On top of this, future research needs to move toward quantification of both bacterial abundance and ecological processes. This is important as bacterioplankton abundances can vary hugely over small spatial scales. Indeed, in the Celtic Sea abundances of Synechococcus can vary by orders of magnitude over 10s of kilometres (Martin et al., 2005(Martin et al., , 2008. Coupled with high-resolution biogeochemical context, this will us to accurately quantify the contribution that bacterioplankton communities associated with oceanographic features, such as the CF, make toward regional processes and biogeochemical cycles.

AUTH O R CO NTR I B UTI O N S
NGK and SKM conceived the designed the study. NGK, FM, and JMT conducted all fieldwork work. SSB, AH, ADMD, NJM, IS and JMT conducted all lab work. NGK, SBW and RR analysed the data.
NGK lead the manuscript preparation, and all authors contributed equally to subsequent edits. All authors read and approved the final manuscript.

ACK N O WLE D G E M ENTS
We would like to thank the crew and research technicians of the RV Prince Madog without whom this research would not be possible. Leaders Fellowship (MR/S032827/1).

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare no conflicts of interest.