Influence of Estuarine Water on the Microbial Community Structure of Patagonian Fjords

Introduction

One of the world’s largest fjord systems is located in the Patagonian region of southern Chile (41–56°S) (Iriarte et al., 2010). This fragmented coastal region, which possesses several islands, rivers and channels (Iriarte et al., 2014; Cuevas et al., 2019), is characterized by high primary production due to the occurrence of phytoplankton blooms throughout the year (Iriarte et al., 2016; Montero et al., 2017). During the austral summer in the Patagonian fjords, high primary production results in carbon sequestration (∼5 mol C m^–2 year^–1) that is considered high relative to Patagonia’s Atlantic coast (∼1.6 mol C m^–2 year^–1) (Torres et al., 2011). Rivers and rainfall runoff, together with glacier melting and vertical mixing, explain this high primary production due to nutrient input (Iriarte et al., 2010; Montero et al., 2017; Cuevas et al., 2019). The Patagonian fjord system is also characterized by sharp vertical and horizontal oxygen and salinity gradients (Iriarte et al., 2014; Cuevas et al., 2019), which are produced by advective transport and mixing of Estuarine Water (EW) with Sub-Antarctic Water (SAAW) (Silva and Vargas, 2014; Vargas et al., 2018). This mixing of EW and SAAW will also likely result in the mixing of entire microbial communities and nutrients in a process termed “community coalescence” (Rillig et al., 2015).

The Patagonian fjord region is considered a vulnerable ecosystem that is highly sensitive to climatic change and anthropogenic activity (Iriarte et al., 2010; Iriarte, 2018; González et al., 2019). Climate forcing by the El Niño-Southern Oscillation and the Southern Annular Mode has been linked to inter-annual variations in environmental conditions, such as higher sea surface temperature, changes in wind direction and intensity, and decreased freshwater discharge to the marine system due to drought (Garreaud, 2018; León-Muñoz et al., 2018). These regional variations can affect the marine microbial community, modifying primary production, and they have even been proposed as possible triggers of harmful algal blooms in the area (Iriarte, 2018; Moreno-Pino et al., 2018). For example, a harmful algal bloom of Pseudochattonella cf. verruculosa observed in the summer of 2016 was the largest reported to date in this area (León-Muñoz et al., 2018). Additionally, activities and processes, such as aquaculture, freshwater runoff and increased glacier melt (Iriarte, 2018; González et al., 2019), have had an impact on the carbon and nitrogen nutrient input to these coastal areas (Silva and Vargas, 2014; Vargas et al., 2018), subsequently affecting microbial activity (Olsen et al., 2017; González et al., 2019). However, the overall effect of these processes on the local microbial ecology is still unknown.

Most efforts to study microbial communities in the Patagonian fjords have focused on the northern region (41–45°S) (Gutiérrez et al., 2015, 2018; Iriarte et al., 2016; Montero et al., 2017), Baker Channel (48°S) (Gutiérrez et al., 2015) and Cape Horn area (52–56°S) (Moreno-Pino et al., 2018), while microbial communities in the central area (47–52°S) remain underexplored. Analyses of the microbial ecology in this area have shown high heterogeneity of the phytoplankton community, dominated during spring blooms by the diatom genera Skeletonema, Chaetoceros, Rhizosolenia, and Thalassiosira (Iriarte et al., 2016; Montero et al., 2017). Additionally, oxygen and salinity were recently suggested as drivers of vertical water stratification that affects the bacterial community distribution of the Puyuhuapi fjord in northern Patagonia (44°S) (Gutiérrez et al., 2018). To date, no broader study has focused on the structure and diversity patterns of bacterioplankton across the Patagonian fjord region.

Here we investigate the influence of EW on the microbial community structure and diversity in the central Patagonian fjord system (47–52°S) over a linear span of 450 km. This survey is necessary to identify the general patterns of microbial communities and specific microorganisms that can be used as possible biological indicators of change in the water conditions of this southern ecosystem, which is susceptible to climatic change and anthropogenic activity. The primary objectives of our analyses were to (i) characterize the environmental conditions, such as temperature, salinity, oxygen and chlorophyll-a, in the Estuarine Waters of the central Patagonian fjords; (ii) assess the effect of the environmental conditions and geographic distance on the community structures of both bacterioplankton (bacteria and archaea) and eukaryotic phytoplankton; and finally, (iii) determine the influence of environmental drivers on relevant taxa in the region.

Materials and Methods

Study Area and Sampling

The central Patagonian fjord region sampled in this study stretched 450 km, from the Penas Gulf (47°S) to Nelson Strait (52°S), as shown in Figure 1. The extent to which areas in this region are influenced by Estuarine and Sub-Antarctic Waters (or modified SAAW) differs. The study transect comprised 23 stations from the Penas Gulf (47.7°S, station 2) to Nelson Strait (51.4°S, station 48) and was divided into ten geographic areas according to the information provided by the Chilean Hydrographic and Oceanographic Service of the Navy. From north to south, the sampling stations were as follows: Penas Gulf (station 2); Baker Channel (stations 10, 8, 7, and 5, oriented in an E–W direction); Bernardo (stations 18C and 18A), Iceberg (station 21B), Falcon (stations 30 and 25) and Penguin (station 33) fjords; Messier Channel (stations 19, 24, 35, 43, and 44); Europa (stations 39 and 36) and Peel (stations 72, 75, and 70) fjords; and Nelson Strait (stations 68 and 48) (Figure 1). Sampling was performed onboard the scientific research vessel “Cabo de Hornos” (AGS61), operated by the Chilean navy’s Hydrographic and Oceanographic Service as part of the program “Cruceros de Investigación Marina en Areas Remotas” (CIMAR). Specifically, this CIMAR23f campaign was conducted from October 17 to November 20, 2017.

Figure 1. Geographic location and study area of the CIMAR23f campaign. (A) Chilean Patagonian fjords as part of the Subantarctic region (46–60°S). Red denotes the Chilean political division within South America. The black dashed square delimits the Patagonian fjord study area. (B) Study area (47.5–52°S) and sampling station locations. Colored circles and numbers represent the location and identification number, respectively, of each station. South Pacific Ocean (SPO), South Atlantic Ocean (SAO).

Surface water at the 23 stations was sampled at two depths (5 and 10 m) using a Niskin bottle-rosette system (46 samples in total). Water samples were analyzed for dissolved oxygen (DO), chlorophyll-a (Chl-a) and used for microbial biomass concentration. Concurrent with water sample collection, temperature and conductivity profiles for the whole column water were obtained using a Seabird 19 CTD. For Chl-a, 1 L water samples were pre-filtered through 200-μm nylon mesh, filtered through 0.22-μm cellulose membrane filters and then stored at −20°C until analysis. The Chl-a concentration (mg m^–3) was subsequently determined fluorometrically using the method described by Parsons et al. (1984). DO (O₂) concentrations (μM) were determined by the modified Winkler method proposed by Strickland and Parsons (1972). For microbial community analysis, seawater (10–15 L) was subjected to serial prefiltration through 200-μm nylon mesh followed by 20-μm polycarbonate filters using a Cole Parmer peristaltic pump system (Model 7553–70, 6 – 600 rpm) at 50 – 100 mL min^–1 for an average of 1 h. Microbial biomass in the filtrate was then concentrated onto 0.22-μm pore size Sterivex units (Durapore; Millipore) and stored at −20°C until analysis.

DNA Extraction

DNA was extracted according to a modified protocol described in Tillett and Neilan (2000). Briefly, the Sterivex membrane filters were manually removed from the housing before DNA extraction. Next, the biomass on the Sterivex filters was resuspended in xanthogenate buffer [1% potassium ethyl xanthogenate (Sigma-Aldrich, United States), 100 mM Tris-HCl (pH 7.4), 20 mM EDTA (pH 8), 800 mM ammonium acetate] with 1% SDS. The mixture was incubated at 65°C for 2 h (hand-shaken every 30 min). After incubation, the tubes were placed on ice for 30 min. The DNA was then extracted with phenol-chloroform-isoamyl alcohol (25:24:1), and the residual phenol was eliminated with chloroform-isoamyl alcohol (24:1). The extract was cleaned by overnight precipitation with cold isopropanol (−20°C) and subsequently washed with 70% ethanol. DNA was quantified using a Qubit 2.0 Fluorometer (Thermo Fisher Scientific, United States). The quality was assessed by spectrophotometry (A₂₆₀/A₂₈₀ ratio) and the integrity was verified by standard 1% agarose gel electrophoresis.

Amplicon Sequencing and Processing

The V4–V5 hypervariable region of the 16S rRNA gene was amplified using the primers 515F (5′-GTGYCAGCMG CCGCGGTAA-3′) and 926R (5′-CCGYCAATTYMTTTRAGT TT-3′), according to Parada et al. (2016). Sequencing was conducted on the Illumina MiSeq platform at the Argonne National Laboratory (Lemont, IL, United States). Raw sequences of the 16S rRNA gene were demultiplexed using the q2-demux plugin implemented in the QIIME2 pipeline (Bokulich et al., 2018; Bolyen et al., 2019). The paired-end sequences were trimmed and merged using DADA2 (Callahan et al., 2016) to obtain amplicon sequencing variants (ASVs). The taxonomy of the 16S rRNA ASVs was assigned using the q2-feature-classifier against the SILVA132 database (Quast et al., 2013) via the “classify consensus vsearch” method (Rognes et al., 2016). ASVs that were taxonomically annotated as “Mitochondria” and “Chloroplast” were excluded from the subsequent bacterial and archaeal community analysis. To obtain the taxonomic identity of eukaryotic phytoplankton (phytoplankton community), ASVs that were annotated as “Chloroplast” were separated and re-annotated using the same methodology applied to the aforementioned ASVs but against the PhytoRef database (Decelle et al., 2015), as described in Fuentes et al. (2019). Finally, to remove singletons and rare taxa, ASVs were filtered by total count sum (>4) and prevalence in at least 10% of the samples (Supplementary Table 1). All 16S rRNA sequences used in this study have been deposited in the NCBI Sequence Read Archive (SRA) under BioProject ID PRJNA670217.

Exploratory Data Analysis

Amplicon sequencing variants analyses were conducted using R packages (R Core Team, 2020), including phyloseq (McMurdie and Holmes, 2013) and ampvis2 (Andersen et al., 2018). Graphical support for the microbial community was conducted with ggplot2 (Wickham, 2016), dendextend (Galili, 2015) and ggmap (Kahle and Wickham, 2013). The Ocean Data View (ODV) software was used for oceanographic variables (Schlitzer, 2020). Maps of the oceanographic conditions by depth were made with a spatial resolution of 1 km (GEBCO 2015 grid) and implemented with coastlines and fill coastlines layers. Diversity analyses were performed with the stats (R Core Team, 2020) and vegan (Oksanen et al., 2019) packages in R. Alpha diversity was measured as the number of distinct ASVs or taxonomic ranks in each sample (observed richness). The Shannon and Pielou indexes were applied exclusively to clean ASVs. For clustering and ordinate analysis, the double-zero effect and rare taxa weight were standardized by Hellinger transformation in ASV filter count data. Hellinger-transformed counts were then used to calculate Euclidean distances among samples (Hellinger distance matrix) (Legendre and Gallagher, 2001). Oceanographic variables were standardized using the z-score method (mean 0, variance 1). Standardized variables were then used to calculate Euclidean distances among the samples. The Ward2 algorithm (Murtagh and Legendre, 2014), implemented in the vegan package, was used to construct the final dendrogram. Silhouette validation criteria, combined with calculated p-values for hierarchical clustering (multiscale bootstrap) implemented in the pvclust package (Suzuki et al., 2019), were used to select the number of clusters (k) in each dataset (Rousseeuw, 1987). Exploratory analysis of the main factors affecting the microbial community structure (i.e., prokaryotic plankton and eukaryotic phytoplankton) was conducted via principal component analysis (PCA) using the ‘envfit’ test function in the vegan package.

Network Reconstruction

Briefly, microbial network analyses were conducted in R using spieceasi (Kurtz et al., 2015) and plotted with ggnet (Briatte, 2020). The best taxonomic hit grouped the clean prokaryotic and eukaryotic ASVs. For network reconstruction and interaction estimation, the neighborhood selection method (Meinshausen and Bühlmann method) (Meinshausen and Bühlmann, 2006) was applied with an nlambda of 20. The network was processed in igraph (Csardi and Nepusz, 2006) for further attribute analyses (i.e., degree and betweenness). Regression coefficient values calculated by the SPIEC-EASI method were used to estimate the positive or negative correlation of edges.

Statistical Analysis

The relation between environmental variables (i.e., temperature, DO, salinity, and Chl-a) and ecological indexes (i.e., Shannon and Pielou) were studied via permutational analysis of variance (PERMANOVA) (Anderson, 2001) using the adonis2 function in the R vegan package (9,999 permutations). Classification of fjord areas was included as the explanatory factor for environmental (df = 9, n = 46) and ecological (df = 9, n = 45) indexes. Additionally, Spearman’s rank correlation (rho) was performed among environmental variables (p < 0.01) (Kendall, 1948). To study the relationship among the environmental and ecological dissimilarity matrices, a Mantel test with a Spearman’s rank correlation (rho) was performed in R. Mantel tests (9,999 permutations), facilitated by the mantel function in the vegan package, were performed between the Hellinger distance matrix (i.e., prokaryotic plankton and eukaryotic phytoplankton), Euclidean distance matrix (independent and joined variables) or the matrix of the Haversine distance between geographical coordinates (i.e., latitude and longitude) among samples (Mantel, 1967; Legendre and Fortin, 2010). Adjusted p-values were obtained using the stats package and Benjamini–Hochberg (BH) post hoc tests (Benjamini and Hochberg, 1995). Distance decay analysis was performed as a simple linear regression analysis of the Jaccard similarity index log between samples versus Haversine distance between sites. Statistical analyses were performed using R software with the stats, vegan, hmisc (Harrell, 2020) and corrplot (Wei and Simko, 2017) packages.

To explore spatial autocorrelation in the collected data, a Mantel correlogram analysis was performed with the vegan package function “mantel.correlog,” using the stations’ latitude and longitude values (XY option), and 9,999 permutations (Supplementary Figure 1). The spatial autocorrelation observed for the oceanographic conditions between neighboring samples was then considered using two approaches, as suggested by Dale and Fortin (2002), and applied to all statistical analyses in this study. First, the p-values were adjusted to a more restrictive type I error (i.e., p < 0.01) for non-permutational analyses [i.e., Spearman’s rank correlation (rho), linear regression and PCA envfit]. Second, when possible, analyses that allow randomization or permutation methods were applied (i.e., PERMANOVA, Mantel test and bootstrap).

Results

Oceanographic Conditions of the Subantarctic Patagonian Fjords

Oceanographic conditions in the Patagonian fjord region showed spatial autocorrelation (Supplementary Figure 1) and an evident freshwater influence along the study area, as illustrated by salinity, temperature and DO profiles (Supplementary Figure 2). The vertical distribution of salinity evidenced a two-layer structure with a remarkable halocline within the first meters (∼0–30 m), especially at sites directly influenced by freshwater (i.e., stations 10, 18A, 21B, 30, 39, and 72). A similar pattern was observed for the DO and temperature profiles (Supplementary Figure 2), with a clear oxycline and thermocline in the first meters (∼0–30 m) of the water column. In addition, temperature profiles displayed thermal inversions at some stations (Supplementary Figure 2).

At 5 and 10 m depth, where microbial communities were studied, temperature showed an overall similar pattern, with more considerable differences between depths at the Falcon, Iceberg and Penguin fjords (Figure 2A). Along the study area, mean temperature values ranged between 7.05°C (Falcon Fjord) and 10.86°C (Iceberg Fjord) (Table 1), with significant differences among defined areas (PERMANOVA R² = 0.97; p = 0.0398). Salinity values at 5 and 10 m were similar (Figure 2B), except for the Baker Channel area (Figure 2B). The study sites did not show significant differences in salinity (PERMANOVA R² = −10.24; p = 0.9751), with mean values normally below 28 (Table 1). Only stations influenced by oceanic water, specifically those located in the Messier Channel and Penas Gulf areas (Figure 1), showed salinity values above 28 (Figure 2B and Table 1).

Figure 2. Oceanographic conditions of the Chilean Subantarctic Patagonian fjord region. Study area and station locations during the CIMAR23f campaign from October 17 to November 20, 2017 (austral spring). (A) Temperature, (B) Salinity, (C) Dissolved oxygen, and (D) Chlorophyll-a.

Table 1. Variability of oceanographic conditions in the Patagonian fjords at 5 and 10 m depth.

The mean DO concentrations in this study ranged between 281.74 and 402.67 μM (Table 1), with non-significant differences between areas (PERMANOVA R² = 0.80; p = 0.3567). DO concentrations and saturation percentages demonstrated a well-oxygenated surface layer, with higher DO values at 5 m, especially in the Falcon and Penguin fjords (Figure 2C and Table 1), rather than at 10 m depth where several areas (Baker Channel, Bernardo, Iceberg and Falcon fjords) showed undersaturated mean values (oxygen saturation <100%, Supplementary Table 2). Chl-a values over 3 mg m^–3 were detected at both 5 and 10 m in the Messier Channel (stations 35, 43, and 44) and Peel Fjord (stations 70 and 72) areas (Figure 2D). Meanwhile, the lowest Chl-a concentrations (<1 mg m^–3) were observed in the Penas Gulf, Baker Channel and Nelson Strait areas (Figure 2D). Although mean values ranged from 0.12 to 9.60 mg m^–3, non-significant differences between areas were detected for Chl-a (PERMANOVA R² = −1.27; p = 0.8678) (Table 1).

Among the environmental variables (temperature, salinity and DO), non-significant (p > 0.01) correlations were detected. However, Spearman’s correlation analysis showed significant associations between depth and these variables (i.e., temperature, salinity and DO). For instance, depth was negatively correlated with DO (R = −0.463; p < 0.01) and temperature (R = −0.364; p < 0.05), but positively correlated with salinity (R = 0.477; p < 0.001) (Supplementary Figure 3). Temperature, salinity and DO conditions in the region were used to identify the presence of estuarine freshwater (EFW; S < 11), estuarine brackish water (EBW; S = 11–21) and estuarine saltwater (ESW; S = 21–31) in the first ∼24 m, followed by modified Subantarctic water (MSAAW; S = 31–33) between ∼25–69 m and Subantarctic water (SAAW; S > 33) below 70 m (Supplementary Figure 2).

Effect of Oceanographic and Geographic Factors on the Microbial Community Structure

Alpha diversity analyses were performed to determine the microbial community diversity of species within each fjord. The bacterioplankton communities and eukaryotic phytoplanktonic (hereafter referred to as phytoplankton) communities showed similar alpha diversity patterns among the study areas, with the highest observed richness values observed in the Messier Channel, Peel Fjord and Baker Channel, and the lowest observed richness values in the Penas Gulf (Supplementary Table 3). PERMANOVA, which was applied to investigate the difference between fjord areas, showed a significant difference in the diversity (Shannon index) and evenness (Pielou index) of ASVs from the bacterioplankton communities (PERMANOVA: Shannon index R² = 0.56, p = 0.0003; Pielou index R² = 0.40, p = 0.0252). Conversely, a non-significant difference was found for phytoplankton using the same indexes (PERMANOVA: Shannon index R² = 0.36, p = 0.0617; Pielou index R² = 0.29, p = 0.1756).

To determine the variation in species composition (i.e., beta diversity) among samples in the Patagonian fjords and the effect of oceanographic variables, a PCA plus envfit analysis was performed. The PCA results show that the first three axes (i.e., PC1, PC2, and PC3) explained 73.8 and 46.6% of the variance for bacterioplankton and phytoplankton, respectively (Figure 3). DO was identified as the main factor inducing variation in the diversity of both communities in the Patagonian fjords (Figure 3). In the Europa Fjord, bacteria and phytoplankton responded differentially to DO concentrations. Additionally, lower salinity values showed an influence on Baker Channel bacterioplankton (Figure 3A). Accordingly, salinity seemed to induce separation of the bacterioplankton assemblages from the Baker Channel and Penas Gulf (Figure 3A). This effect was more pronounced for the bacterioplankton communities located in the Baker Channel at 5 m depth, where the mean salinity value was 15 (Table 1).

Figure 3. Principal component analysis (PCA) of the 45 samples analyzed for bacterioplankton (A) and phytoplankton (B) communities. The left side of panels (A,B) shows PC1 and PC2, while the right side shows PC1 and PC3. Samples were collected from 5 m (triangles) and 10 m (circles) depth. Only oceanographic variables with values of significance (p < 0.01) via envfit were plotted (black arrow and red labels). The relative contribution (eigenvalue) of each axis to the total inertia in the data is indicated as a percentage on the axis titles.

Beta diversity-based cluster analyses of samples were performed to determine regional diversity (Figure 4 and Supplementary Figure 4). The results reveal two separate clusters of bacterioplankton (B1 and B2) and four clusters of phytoplankton (P1–P4) (Figure 4A). A geographical division of the Patagonian fjords is reflected in the cluster distribution of both communities (Figure 4A). The more northern stations, represented by samples from Penas Gulf (47°S) and Baker Channel (48°S), grouped with samples from south of Messier Channel (51°S) (clusters B2 and P4). The central area from Bernardo Fjord to Falcon Fjord (48.8–49.5°S), instead grouped exclusively as part of clusters B1 and P1. Finally, stations in the southern fjord region (50–52°S; Figure 1) were heterogeneously represented by a high diversity of groups dispersed within the six different clusters.

Figure 4. Cluster analysis of the Patagonian fjord microbial community distribution. (A) Geographical distribution of bacterioplankton clusters (colored areas labeled as B1 and B2) and phytoplankton clusters (colored areas labeled as P1 to P4) throughout the Chilean Patagonian fjords. Baker Channel (BK), Bernardo Fjord (BF), Messier Channel (MC), Iceberg Fjord (IF), Falcon Fjord (FF), Penguin Fjord (PF), Europa Fjord (EF), Peel Fjord (PlF), Nelson Strait (NS). Última Esperanza Sound is also shown due to its substantial fluvial contribution to the region. (B) Beta diversity correlation of communities. Mantel test with Spearman’s correlation (rho) of the Hellinger distance matrix between the microbial communities. (C) Microbial community distance decay was calculated as the linear regression of the log Jaccard similarity index and the linear distance among samples. Bacterioplankton communities (dark blue). Phytoplankton communities (red).

The effect of oceanographic variables on beta diversity was then assessed using a Mantel test with Spearman’s correlation (rho). The results confirm that the environment modulates microbial community clusters differentially. For instance, the bacterioplankton community structure of cluster B2 correlated much stronger with salinity and DO (R = 0.5875, p = 0.0018) than that of cluster B1 (R = 0.243, p = 0.0283) (Table 2). Meanwhile, phytoplankton of cluster P1 showed a weak but significant correlation with temperature (R = 0.3205, p = 0.0093), and those of cluster P4 showed a moderate correlation with salinity (R = 0.4729, p = 0.0283) (Table 2). Finally, geographical distance (i.e., Haversine distance) between samples, as well as the correlation of diversity between communities, were studied. Geographic distance showed a weak correlation for both communities (bacterioplankton R = 0. 2398, p = 0.0024; phytoplankton R = 0.2265, p = 0.0018), thus their dissimilarity weakly correlated with the distance between samples. However, a high correlation was observed between the bacterioplankton and phytoplankton communities (R = 0.687, p = 0.0018) (Supplementary Table 4 and Figure 4B).

Table 2. Correlation of Patagonian fjord microbial clusters with oceanographic parameters.

Microbial Core Taxa, Indicators and Interactions in the Patagonian Fjords

To identify the microbial taxonomic composition, and the impact of environmental variables on specific taxa, we focused on identifying the microbial core taxa based on the decreasing prevalence of taxa among samples (Supplementary Table 5). The results showed that no bacterioplankton or phytoplankton ASVs were shared by all samples (n = 45), signifying a highly variable microbial community in the study area. However, when analyzing bacterioplankton at the species level, as well as phytoplankton at the phylum level, it was possible to detect core taxa (Supplementary Table 5). Fifteen bacterioplankton families (9 orders), dominated by Flavobacteriaceae (Flavobacteriales), Rhodobacteraceae (Rhodobacterales) and Pelagibacteraceae SAR11 Clade I (SAR11 clades) (Supplementary Table 6), represented the core taxa and the most abundant bacterioplankton taxa in these Patagonia fjords. In contrast, although the phytoplankton composition was dominated by the families Thalassiosiraceae, Skeletonemataceae, and Bathycoccaceae (Supplementary Table 7), none represented the core (Supplementary Table 5). To confirm if the observed community differences respond to geographic barriers, a distance decay analysis was performed. The results showed a non-significant (p > 0.01) association between the shared ASVs (Jaccard similarity index) and geographical distance (km) between sites (ln prokaryotic R² = 0.007, p = 0.03; phytoplankton R² = 0.008, p = 0.04) (Figure 4C). The same dispersion analysis performed at different taxonomic levels shows a higher dispersion of taxa in the Patagonian fjord communities.

Additionally, we analyzed the effect of environmental factors on specific taxa to identify potential microbial indicators of environmental change within the Patagonian fjords. Here, we defined microbial indicators as those families with a moderate to strong correlation (i.e., R ≥ 0.5 or ≤−0.5) with one or more oceanographic factors. The analysis identified 29 families (dominated by bacteria) that were influenced by at least one oceanographic factor, while eight families were affected by more than one factor. As expected, DO was identified as a critical factor, affecting 14 families (bacteria 12; archaea 2) (Figure 5). Additionally, several families including Gimesiaceae, Beijerinckiaceae and the NS9 marine group, showed a conserved pattern of correlations and were significantly associated with low temperature and low DO concentrations. In contrast, the bacterial families SAR11 Clade III, Burkholderiaceae and Microbacteriaceae were negatively correlated with salinity and Chl-a. These bacterial families also showed a positive correlation with the eukaryotic families Coccomyxaceae and Pyrenomonadaceae. Furthermore, Coccomyxaceae and Pyrenomonadaceae, which were negatively correlated with Chl-a, were the only phytoplankton families to show a significant correlation with an abiotic factor (Figure 5).

Figure 5. Indicator families from the Patagonian fjord microbial communities. Correlation analysis between environmental variables and fjord families. Only taxa with a significant correlation (p < 0.01) and a moderate to strong effect (R ≥ 0.5 or ≤–0.5) with at least one oceanographic factor are shown. The color gradient represents the Spearman’s correlation coefficient values. Only values with significance (p < 0.01), after Benjamini–Hochberg (BH) FDR correction, are shown. The labels are colored by domain according to the legend at the bottom of the plot. Clade III refers to Pelagibacteraceae subgroup III.

Finally, the relationship previously observed in the Mantel test (Figure 4C) between bacterioplankton and phytoplankton taxa (612 and 33 ASVs, respectively) was analyzed in more detail via network reconstruction (Figure 6). The analysis yielded high modularity within the community. This network was mainly dominated by generalist bacteria present in more than one cluster (Figure 6A), while most connected nodes belonged to ASVs within the major bacterioplankton groups (clusters B1 and B2). All phytoplankton ASVs were located, with high connectivity, at the network center (Figure 6A). According to the beta diversity analysis, phytoplankton ASVs that were abundant in cluster P1 correlated more with cluster B1 while ASVs from cluster P4 interacted more with cluster B2. Subsequently, to better understand specific inter-domain interactions, we performed a new correlation analysis by selecting all bacterioplankton families that correlated (i.e., R ≥ 0.5 or ≤−0.5) with phytoplanktonic taxa. The results show interactions between 27 bacterioplankton and 10 phytoplankton families (Figure 6B). The Pyrenomonadaceae and Coccomyxaceae families accounted for most of the positive correlations with bacteria, especially Schleiferiaceae, Spirosomaceae, Sporichthyaceae, and TRA3-20 (Figure 6B). At the same time, Noelaerhabdaceae (order Isochrysidales) exhibited stronger negative correlations with bacteria (i.e., Methylophilaceae) (Figure 6B).

Figure 6. Patagonian fjord microbial network and inter-domain interactions. (A) Microbial co-occurrence network. The node color (ASV) was assigned according to the group with the highest relative abundance observed for the node. Vertices are colored according to the positive (green) or negative (red) correlation between nodes. (B) Correlation between microbial families influenced by interactions with phytoplankton (R ≥ 0.5 or ≤–0.5). The color gradient presents the degree of pair Spearman’s correlation, with blue representing positive correlations and dark yellow representing negative correlations. Only values of significance (p < 0.01), after Benjamini–Hochberg (BH) FDR correction, are shown. The labels are colored by domain according to the legend at the bottom of the plot. Clades I and III refer to Pelagibacteraceae subgroups I and III, respectively.

Discussion

Oceanographic Heterogeneity in the Patagonia Fjord

Oceanographic conditions along the study area showed a clear stratification in the water column. This stratification is associated with the presence of the Estuarine Water (EW) in the first 24 m, followed by the Modified Sub-Antarctic Water (MSAAW) and deeper down by the Sub-Antarctic Water (SAAW) as has been previously described in the area (Silva and Vargas, 2014). Additionally, 3 types of EW, previously identified in the area (Silva and Guzmán, 2006), could be distinguished in the first meters of the water column. At 5 and 10 m depth, the Estuarine Salty Water (ESW) clearly dominates in the study area, with only Estuarine Brackish Water (EBW) present at 5 m in the Baker Channel and the Iceberg fjord. In agreement with previous results from central Patagonia (∼47–53°S) (Silva and Vargas, 2014; Torres et al., 2014), a low salinity, high oxygen surface layer was detected throughout the study area, even in stations located far from a direct freshwater input (i.e., stations 2, 43, 48, and 68), denoting a strong influence of freshwater along the region. Substantial fluvial input has been described in the region, mainly related to Baker Channel and Última Esperanza Sound (González et al., 2013; Silva and Vargas, 2014). The freshwater contribution from rivers and glaciers is characterized by low nitrate and phosphate, but high silicate and organic matter, playing a key role in local dynamics through mixing/stratification processes and influencing the phytoplankton community (González et al., 2013; Silva and Vargas, 2014; Vargas et al., 2018). The salinity and DO conditions associated with EW have been shown to influence the establishment of estuarine microbial communities in northern Patagonia (41°S), delimiting their ecological niches (Gutiérrez et al., 2018). This agrees with our observation in the central Patagonian area, where salinity and DO are also significant drivers of the bacterioplankton structure. These results show a different response of the fjord microbial community structure relative to the general pattern observed in other oceanographic regions where temperature is the main driver (Sunagawa et al., 2015). Our results also confirm the presence of Estuarine Waters, mainly Estuarine Salty Water (ESW, S = 21–31) along the study area, highlighting a conserved pattern of estuarine influence along 450 km of the Chilean Patagonian fjord region. Freshwater input with low nutrients (nitrate and phosphate) and high sediment material, which limits light penetration in the column water, probably negatively affects phytoplankton biomass. In this regard, the lowest Chl-a concentrations observed in this study, namely in the Baker Channel and Nelson Strait areas, were correlated with a direct freshwater influence. In contrast, the highest Chl-a values were detected along the Messier Channel area (stations 19, 24, 35, 43, and 44), where local dynamics allow for the arrival of SAAW, which is characterized by nutrient-rich waters (Silva and Vargas, 2014; Cuevas et al., 2019). The relationship between low Chl-a and high freshwater input has been previously reported in several areas of Patagonia (Vargas et al., 2018; Cuevas et al., 2019) and other glacier regions, such as Maxwell Bay (King George Island, Antarctic Peninsula) (Meredith et al., 2018).

Patagonia Fjord Microbial Community Structure

Freshwater and intertidal waters are characterized by a higher microbial richness than marine waters (Wang et al., 2012). This is reflected in the present study, with areas such as Peel Fjord and Baker Channel showing higher richness and sharp environmental gradients in the surface of the water column, which can be explained by the mixing of sediment and microorganisms from freshwater input (Wang et al., 2012; Rocca et al., 2020). In contrast, due to its connection with an enclosed area containing many fjords, the high richness observed in the Messier Channel area was probably associated with the arrival of sediments and EW, thus accumulating the microbial richness of these fjords. The mixing of these microbial communities (i.e., community coalescence) triggers a structural change that shapes the communities of the Patagonian fjords. However, as in other aquatic environments, the mixing process in this region is not random but facilitated by estuarine circulation, specifically driven by the constant mixing of EW and SAAW that can generate sharp gradients, as observed in this study. This prompts constant community coalescence and the selection of species better adapted to dispersion and mixing (Castledine et al., 2020; Mestre and Höfer, 2020). However, heterogeneous water mixing (i.e., differential mixing between water masses along the fjord system) was also noticed between stations in the study area. These different gradients can also induce variations in the community coalescence mixing ratio and microbial structure (Rillig et al., 2015; Castledine et al., 2020). This was demonstrated by our cluster analysis, in that geographically heterogeneous community structures respond differentially to the measured oceanographic conditions, even when several taxa are shared between clusters, emphasizing the importance of biotic interactions in these communities. Furthermore, evidence suggests that the microbial coalescence process favors system function over the prevalence of input members, selecting the best-adapted taxa, even from the rare biosphere, to ensure productivity of the system (Rillig and Mansour, 2017; Sierocinski et al., 2017; Rocca et al., 2020). Further analysis is now needed to determine if the primary productivity pattern is conserved within different communities of the Patagonian fjords.

The network analysis suggests high modularity in the Patagonian communities (i.e., bacterioplankton and phytoplankton), observing two divergent modules made up of positive biotic connections within modules and negative connections between modules. Modular communities, such as those observed here, have been proposed to emerge after particular coalescence events. These coalescence events produce mixed communities that are adapted to the new resulting environment, while maintaining most of their original species interactions. However, modularity also can be observed when a structured environment limits the mixing between communities, thereby supporting localized interactions (Castledine et al., 2020). The lack of overlap between the two observed modules and their different ecological potentials could be associated here with the gradients caused by heterogeneous mixing of the water masses in the Patagonian fjords (i.e., EW and SAAW); that is, local conditions limit mixing of the communities. Nevertheless, some strong biotic connections were observed between bacterioplankton and phytoplankton at the community and family levels. In the Baker Channel, some localized correlations associated with its differential and specific ecological conditions were also observed. The low dispersion barriers in the area allow dispersion of microorganisms in conjunction with local selection, limiting microbial interactions. In that sense a few ubiquitous and more generalist (i.e., better-adapted) families were distributed throughout the microbial communities of the different fjords. New experiments are now needed to better address the differential mixing ratios affecting the community structures, interactions and productivity in the Patagonian fjords.

Along the 450 km of the Patagonian fjords studied here, the bacterioplankton community was characterized by the dominance of Flavobacteriales, Rhodobacterales and SAR11 Clade I (Pelagibacteraceae). These bacterial groups have been previously described as abundant in this region (Gutiérrez et al., 2018) and in other fjords, such as Arctic Svalbard (Cardman et al., 2014). Flavobacteriales and Rhodobacterales are commonly attached to phytoplankton cells during blooms (Buchan et al., 2014) and are specialized for algal polysaccharide remineralization (Teeling et al., 2016; Kappelmann et al., 2019). In that sense, the high abundance of these bacteria, which were collected during the spring of 2017, could be related to the typical seasonal occurrence of spring phytoplankton blooms (Gutiérrez et al., 2018; Cuevas et al., 2019; Saldías et al., 2019). In agreement with Gutierrez et al., who studied the surface waters of Puyuhuapi fjord in northern Patagonia (Gutiérrez et al., 2018), our results show that the abundant and global marine generalist SAR11 is ubiquitous in the Patagonian fjords. However, the present study also reveals the presence of different sub-clades of SAR11 that are associated with different niches. In particular, SAR11 Clade I (i.e., ecotype SAR11 Clade Ia.1), which has been mainly associated with cold marine waters at high latitudes, was the most abundant sub-clade present in the Patagonian fjords. Meanwhile, Clade III (i.e., ecotype SAR11 Clade III.b), related to freshwater (Brown et al., 2012; Giovannoni, 2017), was associated with the low salinity (i.e., EW) and low Chl-a waters, such as those observed in the Baker Channel. Nevertheless, the conclusions presented in this study should be taken with a degree of caution given the inherent limitations of current molecular techniques applied to measure the full diversity or taxonomy of microbial communities. Limitations in DNA extraction and library construction that may over- or under-represent estimates of microbial communities members have not yet been solved (Oyola et al., 2012; Vincent et al., 2017). In addition, a major technical limitation in this study was the use of a single molecular marker (i.e., the 16S rRNA gene), which is unable to resolve the finest taxonomic discrepancies and relies heavily on database annotation (Rosselló-Móra and Amann, 2015).

The Patagonian fjord phytoplanktonic community was dominated by the orders Thalassiosirales (Thalassiosiraceae and Skeletonemataceae) and Mamiellales (Bathycoccaceae), which have been previously described in the region (Takeuchi and Kohshima, 2004; González et al., 2013; Iriarte et al., 2013; Fernández et al., 2017; Montero et al., 2017; Fuentes et al., 2019). However, the families Pyrenomonadaceae and Coccomyxaceae exhibited a higher correlation with bacterioplankton families. Pyrenomonadaceae and Coccomyxaceae were especially abundant in areas of low salinity and low Chl-a, such as the Baker Channel. These two families are mainly comprised of genera identified in terrestrial and freshwater environments (Guiry and Guiry, 2020), supporting the introduction of microorganisms into these fjords by freshwater inputs. The distinctive oceanographic conditions found in the Baker area explain the segregated microbial communities when compared to the rest of the Patagonian fjord areas and reveal possible specific interactions between the phytoplankton families Coccomyxaceae and Pyrenomonadaceae and some bacterial families (i.e., Schleiferiaceae, Spirosomaceae, Sporichthyaceae, TRA3-20, and SAR11 Clade III). How these interactions impact the system, and if these interactions are direct (i.e., mutualism) or indirect (i.e., trophic cascades), should be addressed in future studies.

Microbial Indicators of Environmental Changes in the Patagonia Fjords

Microorganisms are crucial for the formation of an ecosystem. Thus, the essential role of the microbiome within the environmental structure can be studied and associated in relation to normal or perturbed ecological conditions (Glasl et al., 2017; Astudillo-García et al., 2019). Consequently, microbiomes have been used to assess the impact of climate change in highly sensitives ecosystems, such as agricultural soils (Schloter et al., 2018), coral reefs (Glasl et al., 2017), and fjords (Gutiérrez et al., 2015; Eregno et al., 2018). Microbial communities respond rapidly to environmental changes, making them a potential tool to complement current monitoring systems (Glasl et al., 2017; Rocca et al., 2020). However, both the environment and microbial communities are highly dynamic, so the use of these communities as biological indicators presents different challenges, such as understanding system stability and how temporal variation impacts the community structure (Astudillo-García et al., 2019). The use of microbial communities or specific taxa to analyze the health of a natural system will largely depend on our knowledge of the microbial diversity, ecology and stability of the system (Astudillo-García et al., 2019). Exploring the diversity of the microbial diversity is important for understanding their complexity, variation and conservation within an specific environment, which will ultimately improve our understanding of how microorganisms interact with the system (Colwell, 1997; Astudillo-García et al., 2019). Efforts to obtain information on the microbial communities of a given system also help to establish a reference for the analysis of other similar regions.

In the Patagonia region, inter-annual variations (e.g., El Niño-Southern Oscillation and Southern Annular Mode) (Garreaud, 2018) can impact the planktonic community, causing increased harmful algal blooms (León-Muñoz et al., 2018). In this study, the negative associations of SAR11 Clade III with salinity and Chl-a conditions highlight this taxon as a potential indicator of freshwater input into the fjords, especially in the Baker Channel area, which features high freshwater runoff, sediment load penetration and decreased light penetration (González et al., 2013). SAR11 Clade III has not only been previously associated with freshwater (Giovannoni, 2017), but it has also been suggested as a tool to monitor climate change mediated by ocean warming (i.e., by following the migration of tropical ecotypes to polar regions) (Brown et al., 2012). Considering that global climate change is influencing precipitation in the Patagonia region (Garreaud, 2018), with precipitation trends decreasing over north-central Patagonia and increasing south of 50°S (Garreaud et al., 2013), changes in the abundance of SAR11 Clade III could also be a potential indicator of possible climate change effects in the region.

In contrast to the conditions observed in freshwater-influenced stations, Beijerinckiaceae, Gimesiaceae (PVC group; Planctomycetes) and the NS9 marine group (Flavobacteriia; Flavobacteriales) were identified as potential indicators of the low temperature and low DO waters upwelled from the MSAAW. These bacteria have been associated with different depths in the water column (Muck et al., 2019). The NS9 marine group has been commonly found below the deep chlorophyll maximum (DCM) (Cram et al., 2015). Planctomycetes are highly diverse (i.e., genetically) and widespread in aquatic environments (Dedysh and Ivanova, 2019), in which the Gimesiaceae family, in particular, has been associated with deep waters, mainly in lakes (Rojas-Jimenez et al., 2019). In this region, due to the potential changes in freshwater input under the current climate change scenario, alterations to the typical salinity and oxygen concentration could disrupt the dominance of certain taxa in the surface of the column water. This could modify the interactions between surface and subsurface microorganisms in the Patagonian fjords. Further analyses are needed to fully understand how these changes in community structure and interactions could impact the high productivity of the Patagonian fjord system.

Conclusion

The microbial dispersion observed in the central Patagonian fjord area could be due to both estuarine circulation and the normal coalescence process that homogenize microbial communities over hundreds of kilometers within this ecosystem. However, particular local conditions seem to influence the interactions of microorganisms that are possibly co-evolving in these communities. Thus, our results suggest that Patagonia fjords are an ideal location to study selective ecological pressure on microbial coalescent communities. In addition, the microbial communities in Patagonia showed that while high dispersion homogenizes the community, established biotic interactions allow bacterioplankton and phytoplankton to respond differentially to environmental factors. Thus, the significant association identified in this study between bacterioplankton and phytoplankton (i.e., beta diversity and family-level microbial composition) suggests that although abiotic factors directly or indirectly affect particular taxa within the microbial assemblages (e.g., indicators), biological interactions can permeate the entire community structure of the Patagonian fjords. However, particular taxa, such as the SAR11 Clade III, Burkholderiaceae and Microbacteriaceae families, should be further analyzed and potentially used as indicators of freshwater input into the Patagonian fjord region, an ecosystem highly susceptible to climate change. Meanwhile, Gimesiaceae, Beijerinckiaceae and the NS9 marine group could be used to monitor the impact of marine waters upwelled from the SAAW or MSAAW in these fjords. Due to changes in freshwater input under the current climate change scenario, alterations to the typical salinity and DO concentration may disrupt common interactions between surface and subsurface microorganisms in the region. Altogether, this study of the marine microbial communities in the Patagonia fjords, as indicators of environmental variations in Subantarctic fjords, may help us to understand the potential impact of climate change on Patagonian biodiversity and productivity, and ultimately in other fjords worldwide.

Data Availability Statement

The datasets generated for this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: National Center for Biotechnology Information (NCBI), Sequence Read Archive (SRA), BioProject accession ID PRJNA670217.

Author Contributions

JT-L collected and processed the samples in the field. JT-L, JC-A, PA-R, JA, and IM processed and analyzed sequence data, analyzed oceanographic data, helped with the discussion, and critically reviewed and edited the manuscript for journal submission. JT-L and BD wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This research was funded by “Cruceros de Investigación Marina en Areas Remotas (CIMAR)” – CIMAR Fiordos 23 (CONA C23F 17-02) of the “Servicio Hidrográfico y Oceanográfico de la Armada” (SHOA), the Center for Climate and Resilience Research (CR)2 (ANID/FONDAP/1511000), and the Center for Genome Regulation (CRG) (ANID/FONDAP/15200002). Additionally, JT-L was supported by the National Agency for Research and Development (ANID) National Doctoral Scholarship (no. 21171048) during the analysis and writing of this manuscript. JA thanks ANID National Doctoral Scholarship (no. 21130515) for additional funding support.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to thank the scientific research ship “Cabo de Hornos” (AGS61) of the Chilean navy, their crew and the “Servicio Hidrográfico y Oceanográfico de la Armada” (SHOA) for their help and logistical support during the CIMAR23f campaign. We also thank Christina Ridley for her help in correcting the English version of the manuscript, and the BLOOM laboratory team, especially Antonieta Jojot, Juan Ayala, and Viviana Catalan, for sampling and sample processing. The comments of the editor and reviewers helped us to improve the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.611981/full#supplementary-material

References