Linking bacterial and archaeal community dynamics to related hydrological, geochemical and environmental characteristics between surface water and groundwater in a karstic estuary

1. Introduction

The subterranean estuary (STE) is a reactive mixing zone between fresh groundwater and recirculated seawater (Moore, 1999; Rocha et al., 2021), which is considered to be a significant contributor of bacteria (Boehm et al., 2004; Knee et al., 2008; Ruiz-González et al., 2021) and biogenic elements such as nutrients (Cai et al., 2015; Wang et al., 2018; Chen et al., 2021a; Santos et al., 2021; Zhao et al., 2021), carbon (Liu et al., 2014; Wang et al., 2014; Chen et al., 2018b, 2022; Yau et al., 2022), metals (Wang et al., 2019; Mayfield et al., 2021; Zhong et al., 2022) and greenhouse gases (Chen et al., 2021b, 2023; Reading et al., 2021; Zhu et al., 2022) to coastal waters. Biogeochemical reactions in STEs are often mediated by microorganisms, which can modify the chemical composition of submarine groundwater (Santoro et al., 2008; Adyasari et al., 2019, 2020; Chen et al., 2019, 2020b; Zhang et al., 2021). For example, bacterial processes contribute to the nitrate removal and dissolved inorganic carbon (DIC) addition in a sandy STE (Chen et al., 2020b). However, compared with inland groundwater, the microbial dimension of submarine groundwater discharge (SGD) remains poorly understood (Archana et al., 2021; Ruiz-González et al., 2021).

As the most abundant and diverse group of life on the Earth, microbial communities have an integral function and play a pivotal role in biogeochemical cycling to marine ecosystem, and that has no exception for SGD (DeLong and Karl, 2005; Ye et al., 2016; Adyasari et al., 2020). Several studies of microbial communities in STEs were reported and novel insights involve using 16S rRNA (Ye et al., 2016; Adyasari et al., 2019, 2020; Chen et al., 2019, 2020b; Wu et al., 2021) and functional genes (Santoro et al., 2006, 2008) to reveal microbial communities in STE systems. Large amounts of nutrients, carbon and greenhouse gases entering the coastal waters via SGD may lead to outbreaks of eutrophication (Hwang et al., 2005), algal blooms (Chen et al., 2020a; Zhao et al., 2021), deoxygenation (Guo et al., 2020) and localized coastal acidification (Wang et al., 2014; Cardenas et al., 2020). However, the coupling of microbial communities and environmental problems (e.g., red tides) in SGD is less reported, which limits the understanding on the biogeochemical cycle and environmental management in STE ecosystems.

The Krka River Estuary is a salt-wedge, highly stratified estuary, located in the central part of the eastern Adriatic coast (Croatia) (Liu et al., 2019). As in other Mediterranean karstic shorelines, anchialine environments are common in the Krka River Estuary. These types of environments are characterized by highly stratified water column and connection with the open sea, usually through the karst carbonate rock, showing noticeable marine as well as terrestrial influences (Žic et al., 2008; Kwokal et al., 2014). Previous ²²²Rn mass balance has shown that anchialine caves or springs can deliver point-source, nitrogen enriched groundwater with high N/P molar ratio to the coastal waters in the Krka River Estuary, which trigger and sustain red tide outbreaks (Chen et al., 2020a). Therefore, these anchialine caves or springs are reasonable and representative groundwater end-members in the Krka River Estuary (Liu et al., 2019; Chen et al., 2020a). As one kind of the important STEs, anchialine system has complex and variable hydrological, geochemical and environmental characteristics in the Krka River Estuary and provides direct access to the aquifer source (Bishop et al., 2015). However, the coupling relationship between these complex characteristics and bacterial and archaeal community dynamics is still poorly understood despite their potentially important implications for nutrient and carbon cycling.

Here, we hypothesize that bacterial and archaeal community dynamics interact with the hydrological, geochemical and environmental characteristics in STEs. We investigated the bacterial and archaeal diversities by utilizing 16S rRNA-based Illumina MiSeq sequencing technology and related salinity, ²²²Rn (a groundwater tracer), nutrients and carbon in submarine groundwater (including anchialine caves and springs) and surface water along the Krka River Estuary. This study aims to (1) explore coupling relationship between bacterial and archaeal community dynamics and related hydrological, geochemical and environmental characteristics (e.g., algal blooms) in karstic STEs, and (2) provide insights on the mechanisms of controlling the input or removal of terrigenous materials affected by bacterial and archaeal communities in STE ecosystems, and (3) find the related key candidates for bioremediation in the polluted SGD.

2. Materials and methods

2.1 Study area

The Krka River Estuary is an oligotrophic stratified karstic estuary, which is sensitive to inputs of external substances such as nutrients and carbon (Liu et al., 2019; Chen et al., 2020a). The mean water flow of the Krka River is between 40 m³/s and 60 m³/s (Bonacci et al., 2006). Algal blooms were observed in the estuary and showed obvious seasonal outbreaks (Chen et al., 2020a). The Krka River Estuary has a micro-tidal nature (tidal range: 0.2−0.5 m) with strong stratification. The area is characterized by mild, wet winters and warm, dry summers, with an annual average precipitation of ~900 mm ( http://klima.hr/klima_arhiva.php).

2.2 Sample collection and measurements of physic-chemical parameters

In this study, seven bacterial and archaeal samples (anchialine cave, spring and surface water) were collected along the Krka River Estuary in April, 2016 (Fig. 1). Each sample (~5 L) was collected on a 0.22 µm pore size polycarbonate filter (Nuclepore Track-Etched Membrane, Whatman, UK), then placed in a sterile 1.5 mL microcentrifuge tube and was immediately placed and kept in the fridge (−20℃). ²²²Rn samples were collected in 250 mL glass bottles and analyzed by RAD7 detector immediately. Nutrient and carbon samples were collected with polyethylene bottles, then stored at −20℃ and kept away from light. The nutrients ( ${{\rm {NO}}_2^-} $ , ${{\rm {NO}}_3^-} $ , ${{\rm {NH}}_4^+} $ and ${{\rm {PO}}_4^{3-}} $ ) were analyzed using the method of Strickland and Parsons (1972). The DIC and dissolved organic carbon (DOC) were determined using TOC-V Analyser (Shimadzu®, Japan). Salinity and temperature were measured directly in the field using a multi parametric probe (Hach Lange HQ40D, USA).

Figure  1.  Study site location and sampling stations (a) in the Krka River Estuary (KRE). Algal bloom occurred at surface water stations SW1 (b) and SW3 (c) during sampling (April 2016). SW: seawater; GW: groundwater.

DownLoad: Full-Size Img PowerPoint

2.3 DNA extraction, polymerase chain reaction (PCR), and sequencing

Total DNA was extracted from each filter using a MoBio PowerWater® DNA Isolation Kit (MOBIO Laboratories, Carlsbad, USA). DNA concentration and purity were quantified spectrophotometrically with NanoDrop ND2000 (USA).

To decrease PCR bias, minimum numbers of PCR cycles were performed and three independent PCR mixtures were pooled for each sample. The bacterial 16S rRNA genes were amplified using the specific barcoded universal primer pairs 515F (5′-GTGCCAGCMGCCGCGG-3′) and 907R (5′-CCGTCAATTCMTTTRAGTTT-3′) spanning the V4-V5 hypervariable regions (Xiong et al., 2012). Cycling conditions were an initial denaturation at 95℃ for 2 min, 25 cycles of 95℃ for 30 s, 55℃ for 30 s, 72℃ for 30 s and a final 5-min extension at 72℃. Meanwhile, the archaeal 16S rRNA genes were amplified using specific barcoded universal primer pairs 524F10extF (5′-TGYCAGCCGCCGCGGTAA-3′) and Arch958RmodR (5′-YCCGGCGTTGAVTCCAATT-3′) (Pires et al., 2012) using the following amplification conditions: 3 min at 95℃, followed by 35 cycles with 30 s at 95℃, 30 s at 55℃, 45 s at 72℃, and a final extension period of 10 min at 72℃. PCR products from each tagged primer were purified using the AxyPreDNA gel extraction kit (Axygen Biosciences, USA) and then quantified by QuantiFluorTM-ST (Promega, USA). Finally, reaction mixtures were pooled in equimolar ratios and paired-end reads were generated on an Ilumina MiSeq PE250 (Majorbio Bio-Pharm Technology Co., Ltd., China, http://www.majorbio.bioon.com.cn/).

Raw Illumina FASTQ files were demultiplexed, quality-filtered, and analyzed using Quantitative Insights into Microbial Ecology (QIIME) (version 1.9.1) (Caporaso et al., 2010) as the criteria described previously (Chen et al., 2020b). Reads that could not be assembled were discarded. The UCHIME algorithm in the USEARCH platform was used to detect and remove chimeric sequences, which performs both de novo chimera and reference-based detection (Edgar et al., 2011). Operational taxonomic units (OTUs) with 97% similarity cutoff were clustered using UPARSE (Edgar, 2013). The OTU Cluster was used to determine OTU abundances. The number of reads from each sample that was assigned to each OTU was generated as an “OTU table” by using the Usearch_global command. The sequence reads of all samples were randomly resampled to the smallest sample size. The taxonomic assignment of representative OTU sequences was analyzed by the Ribosomal Database Project Classifier (version 2.11, http://rdp.cme.msu.edu/).

2.4 Phylogenetic analyses

The sequences of the representative OTUs obtained in this study were compared to those in the National Center for Biotechnology Information (NCBI) nucleotide database by using Basic Local Alignment Search Tool (BLAST, Altschul et al., 1997) searching. The closest sequences and selected reference sequences were downloaded and aligned using Clustal W. Phylogenetic trees were generated in MEGA7 using the neighbor-joining method with a bootstrap test of 1000 replicates and maximum composite likelihood model (Tamura et al., 2013).

2.5 Statistical analyses

Alpha diversity metrics and coverage were measured and calculated using the Mothur Program (Schloss et al., 2009). Hierarchical clustering of the samples was performed with the complete method using the function hclust of stats package in R (R Development Core Team, 2013; http://www.R-project.org/). Principal-coordinate analyses (PCoA) were performed to show if distinct separations in bacterial or archaeal community structures were present between groundwater and surface water. To determine which environmental variables best explained patterns of bacterial and archaeal communities, canonical correspondence analysis (CCA) was applied (ter Braak, 1986).

3. Results and discussion

3.1 Site description and environmental characteristics

Surface water and submarine groundwater (spring and anchialine cave water) samples were collected in this study (Fig. 1). Groundwater (including spring water and anchialine cave water) had lower salinity and higher nutrient concentrations than surface water (Table 1). Based on this, seven samples were separated into two groups: Group G included 2 spring samples (GW1 and GW2) and 2 anchialine cave samples (GW3 and GW4), and the salinity of Group G ranged from 0.7 to 3.4 with a mean of 1.5; Group S included three surface water samples (SW1, SW2 and SW3) with relatively higher salinity, which ranged from 3.1 to 9.2 with a mean of 6.2. Group S had mean values of 396 Bq/m³, 79.67 μmol/L and 0.42 μmol/L for²²²Rn, ${{\rm {NO}}_3^-} $ and ${{\rm {PO}}_4^{3-}} $ . However, Group G had higher mean values of 1103 Bq/m³, 90.82 μmol/L and 1.28 μmol/L for²²²Rn, ${{\rm {NO}}_3^-} $ and ${{\rm {PO}}_4^{3-}} $ , respectively. ${{\rm {NO}}_2^-} $ concentrations ranged from 0.03 μmol/L to 0.12 μmol/L in Group G, lower than those in Group S (0.26−0.67 μmol/L). Mean ${{\rm {NH}}_4^+} $ concentration in Group S (1.89 μmol/L) was slightly higher than that in Group G (1.32 μmol/L). There was no distinct pattern of DIC concentration between Group G and Group S. DOC concentration in Group S (1.13 mg/L) was higher than that in Group G (0.68 mg/L). The temperature values (13.5−15.3℃) were relatively constant at all sites.

Table  1.  The site descriptions and physic-chemical parameters

Station	SW1	SW2	SW3	GW1	GW2	GW3	GW4
Longitude	15.8344°E	15.8644°E	15.8276°E	15.8322°E	15.9492°E	15.8601°E	15.8106°E
Latitude	43.8370°N	43.8166°N	43.7820°N	43.7807°N	43.6781°N	43.7373°N	43.7775°N
Water type	surface water	surface water	surface water	spring	spring	anchialine cave	anchialine cave
Characteristics of sampling site	algal bloom happening	located in the center of the lake	near the Station GW1; algal bloom happening	near the village; in the open air	near the lagoon	natural cave; smelly water; water hypoxia (DO: 0.62 mg/L)	natural cave
Temperature/℃	13.5	15.0	ND	14.7	15.3	ND	ND
Salinity	3.1	6.4	9.2	1.1	3.4	0.7	0.8
222Rn/(Bq·m−3)	918	54	215	1748	ND	881	681
${{\rm {NO}}_3^-} $/(μmol·L−1)	178.01	18.50	42.51	178.57	79.66	14.22	ND
${{\rm {NO}}_2^-} $/(μmol·L−1)	0.67	0.32	0.26	0.03	0.03	0.12	ND
${{\rm {NH}}_4^+} $/(μmol·L−1)	2.79	1.87	1.00	1.31	0.50	2.14	ND
${{\rm {PO}}_4^{3-}} $/(μmol·L−1)	0.91	0.20	0.16	0.95	0.23	2.67	ND
DIC/(mg·L−1)	65.62	50.77	ND	61.10	ND	ND	ND
DOC/(mg·L−1)	1.32	0.93	ND	0.68	ND	ND	ND
Note: ND represents not determined.

 | Show Table

DownLoad: CSV

For surface water, the ²²²Rn and ${{\rm {NO}}_3^-} $ concentrations of Stations SW1 and SW3 were significantly higher than those of Station SW2, indicating that SW1 and SW3 were significantly affected by high concentrations of groundwater ²²²Rn and ${{\rm {NO}}_3^-} $ . Algal blooms occurred at surface water stations SW1 and SW3 during the sampling period, which further indicates that groundwater provides the major nutrients such as ${{\rm {NO}}_3^-} $ for algae blooms due to the lack of other significant nutrient sources in the Krka River Estuary (Chen et al., 2020a).

3.2 Microbial diversity and distribution in groundwater-surface water interaction

Totally seven bacterial and archaeal samples were analyzed by using high throughput sequencing, and a total of 195 417 high-quality bacterial V4-V5 Illumina sequences and 139 040 high-quality archaeal V4-V5 Illumina sequences were retrieved. There were 2830 bacterial OTUs and 680 archaeal OTUs in the complete OTU data set at the 97% similarity cutoff. Good coverage was 99.5%−99.9% for all samples after subsampling (Table 2). Richness index (including ACE and Chao) of bacterial community ranged from 332 to 398 in surface water, significantly lower than those in groundwater (434−797). There was no significant difference in archaeal richness between surface water (59−138) and groundwater (81−136). Bacterial and archaeal diversities (e.g., Shannon index) showed variability in both surface water and submarine groundwater (Table 2).

Table  2.  High throughput sequence information and diversity and richness estimators of bacterial and archaeal communities in groundwater (spring and anchialine cave water) and surface water

Station	Optimized sequence	Observed OTU	ACE	Chao	Shannon index	Coverage
Bacteria
SW1	31351	317	384	398	3.03	0.997
SW2	21363	278	342	332	3.50	0.997
SW3	32245	299	381	377	2.11	0.997
GW1	23958	329	455	434	2.78	0.995
GW2	28391	715	780	797	4.72	0.996
GW3	28485	401	473	482	3.71	0.997
GW4	29624	491	579	592	4.13	0.996
Archaea
SW1	33311	138	138	138	2.66	0.999
SW2	33556	53	58	59	0.55	0.999
SW3	31043	105	106	105	2.17	0.999
GW1	31363	115	125	132	1.85	0.999
GW2	35450	134	138	136	1.70	0.999
GW3	41968	55	68	81	0.17	0.999
GW4	32349	80	85	86	1.87	0.999
Note: OTU: operational taxonomic units; ACE: abundance-based coverage estimator; Chao: Chao1 species richness.

 | Show Table

DownLoad: CSV

PCoA analyses showed that both the bacterial and archaeal communities in the areas where algal blooms occurred (i.e., Stations SW1 and SW3) were significantly separated from those in other stations along the first axis explaining 21.2% and 50.3% of the variations (Fig. 2), respectively, which may indicate that the spatial niche partitioning of the bacterial and archaeal communities at Stations SW1 and SW3 and related ecological function were different from other stations (Chen et al., 2019).

Figure  2.  The principal-coordinate analyses (PCoA) analyses on operational taxonomic units (OTU) levels of bacterial (a) and archaeal (b) communities in surface water (green triangles) and submarine groundwater (cyan dots). Notice the algal bloom happening at Stations SW1 and SW3.

DownLoad: Full-Size Img PowerPoint

Taxonomic distributions showed that there were no significant variations in the proportion of Illumina sequences among the bacterial samples, but there were some differences among archaeal samples (Fig. 3). For all bacterial samples, Gammaproteobacteria (29.0%−78.8%), Alphaproteobacteria (3.6%−33.5%) and Bacteroidia (0.3%−27.8%) were the three main bacterial taxa. Cyanobacteriia and Actinobacteria also accounted for a larger proportion at surface water stations SW1 (13.7%) and SW3 (17.2%), respectively. For archaeal samples, the most dominant genus at spring water stations GW1 (51.7%) and GW2 (59.7%) was Candidatus Nitrosopumilus. The genus of Candidatus Nitrosoarchaeum was found abundantly at anchialine cave stations GW3 (97.7%) and GW4 (54.1%) and surface water station SW2 (91.3%), and Candidatus Nitrosoarchaeum was also the second dominant genus at spring water stations GW1 (21.1%) and GW2 (12.6%). However, the genera of Candidatus Nitrosopumilus and Candidatus Nitrosoarchaeum were seldom detected at surface water stations SW1 (<1%) and SW3 (0.5%−2.2%) where algal blooms occurred. Crenarchaeotic Group archaea was the most abundant group at surface water stations SW1 and SW3, which further indicates that Crenarchaeotic Group archaea may be related to algal bloom outbreaks at these stations.

Figure  3.  The relative composition of bacterial (a) and archaeal (b) taxa across all samples. Bacterial and archaeal taxa represented by less than 1% reads are pooled as “others”.

DownLoad: Full-Size Img PowerPoint

In previous studies, Proteobacteria dominated all bacterial groundwater samples, which belongs to Alphaproteobacteria, Betaproteobacteria and Gammaproteobacteria at the genus level (Archana et al., 2021). In our study, the microbial diversity and distribution in groundwater samples were consistent with these previous results (Fig. 3). Although previous data on groundwater archaea are limited, these archaea results indicate that Euryarchaeota, Crenarchaeota, Bathyarchaeota and Thaumarchaetoa were widely distributed in groundwater samples (Archana et al., 2021), but the dominant groundwater archaea in this study was Candidatus Nitrosopumilus (Fig. 3).

3.3 Ecological niches of microbial community in groundwater-surface water interaction

Surface water and groundwater interactions in coastal areas are active mixing zones with characteristic ecological functions (Lee et al., 2017). The dominant bacterial and archaeal communities were significantly different in submarine groundwater and surface water (Fig. 4). For example, the dominant bacterial communities in groundwater samples were Limnohabitans and Psychrobacter while that in the surface water were Flavobacterium and Rickettsiella. The results suggest that the environmental differences between the subterranean and surface estuaries may affect the microbial diversity and distribution, which in turn affect their biogeochemical and ecological functions.

Figure  4.  Ternary diagram illustrating the dominant bacterial (a) and archaeal (b) communities in surface water, spring water and anchialine cave water.

DownLoad: Full-Size Img PowerPoint

Previous studies emphasized that microbes participate in the carbon and nutrient cycling in STEs and determined the chemistry of the groundwater reaching the ocean (Chen et al., 2019, 2020b; Ruiz-González et al., 2021). In this study, the neighbor-joining method was adopted to construct the bacterial phylogenetic tree (Fig. 5). The sequences of OTU32B were frequently obtained from groundwater (proportion of OTUs: 1.3%−45.5%) and exhibited 99.0% similarity with Limnohabitans parvus strain II-B4 (NR_125542), which may serve as ammonium oxidizer (Kasalický et al., 2010). Some OTU sequences such as OTU515B, OTU593B, OTU352B, OTU193B and OTU672B were detected in some groundwater stations (proportion of OTUs: 3.2%−6.5%), which had 99.5%−100% matches to Simplicispira limi strain EMB325 (NR_043773) (Lu et al., 2007), Novosphingobium naphthae strain D39 (NR_151934) (Chaudhary and Kim, 2016), Pseudorhodobacter aquaticus strain DC2N1-10 (NR_151912) (Li et al., 2016), Marivivens donghaensis strain AM-4 (NR_148632) (Park et al., 2016) and Aquirufa antheringensis strain 30S-ANTBAC (NR_165029) (Pitt et al., 2019), respectively. These strains can strongly reduce nitrate through nitrate reduction. These OTU sequences (except OTU193B) were also frequently found in surface water, indicating the connectivity in groundwater-surface water interaction.

Figure  5.  Neighbor-joining tree showing phylogenetic relationships among the major bacterial operational taxonomic units (OTUs) and reference 16S rRNA gene sequences retrieved from the National Center for Biotechnology Information (NCBI) GenBank. These OTUs were those that ranked top five reads from at least a single sample. The scale bar represents the estimated number of nucleotide changes per sequence position. Percentage on nodes refer to the percentage of recovery from 1000 bootstrap resamplings. Only values larger than 50% are shown. The numbers in parentheses indicate the percentage composition of reads in each station in the following order: (SW1, SW2, SW3) (GW1, GW2, GW3, GW4). Methanococcus maripaludis (U38941) was used as the outgroup.

DownLoad: Full-Size Img PowerPoint

In the archaeal analysis, neighbor-joining tree was used to show phylogenetic relationships among the representative OTUs (Fig. 6). The sequences of OTU199A were frequently obtained from groundwater (proportion of OTUs: 1.1%−67.0%) but rarely in surface water (proportion of OTUs: 0−0.5%), and showed >99% similarity with Nitrosopumilus cobalaminigenes strain HCA1 (NR_159206) (Qin et al., 2017), N. maritimus SCM1 (NR_102913) (Walker et al., 2010), N. ureiphilus strain PS0 (NR_15920) (Qin et al., 2017) and N. oxyclinae strain HCE1 (NR_159207) (Walker et al., 2010), which belong to marine ammonia-oxidizing archaea. For example, N. maritimus performs a far unrecognized pathway of ammonia conversion to N₂ when oxygen is depleted (Kraft et al., 2022; Martens-Habbena and Qin, 2022). OTU136A, also known as the ammonia-oxidizing archaea sequence, was found in both groundwater (proportion of OTUs: 14.3%−97.6%) and surface water (proportion of OTUs: 0.3%−91.2%), which had 99.6% match to Nitrosarchaeum sp. AC2 (CP030847). In addition, the sequences of OTU233A, OTU234A and OTU270A were obtained from surface water (proportion of OTUs: 0.3%−45.6%) but relatively few in groundwater (proportion of OTUs: 0−6.0%). These sequences exhibited 95.8%−97.1% similarity with N. viennensis strain EN76 (NR_134097) (Stieglmeier et al., 2014), which is also considered as an aerobic and mesophilic, ammonia-oxidizing archaea. Although archaeal analysis was mainly dominated by ammonia-oxidizing process, different strains showed obvious site selectivity, suggesting that high spatial heterogeneity of dominant archaeal species occurred in groundwater-surface water interaction (Ruiz-González et al., 2022).

Therefore, these main species can utilize excess nutrients (e.g., nitrate and ammonia) at groundwater-surface water interface and reduce the amounts of nutrients entering the ocean, thereby affecting the ecological function of the interface by changing the biomass and community structure of the phytoplankton (Chen et al., 2018a; Torre et al., 2019). The utilization of nutrients by these bacteria and archaea communities seems to match the oligonutrients in the Krka River Estuary (Liu et al., 2019). Note that although nutrients can be utilized by microorganisms in STEs, significant input of point-source nutrients may still lead to water quality degradation. For example, large nutrient discharge with high N/P molar ratios (190−320) from Litno Cave (i.e., Station GW1) likely trigger and sustain red tide outbreaks (i.e., Station SW3) (Chen et al., 2020a).

3.4 Connection between microbial community and hydrochemical environment

Hydrochemical parameters such as salinity, temperature and nutrients have commonly been identified as potential environmental drivers of the observed compositional variations in microbial communities at coastal groundwater-surface water interface (Ye et al., 2016; Chen et al., 2020b; Ruiz-González et al., 2022; Yang et al., 2022). We used CCA analysis (Fig. 7) to determine the interaction between microbial communities and hydrochemical environment, and to clarify which environmental factors govern the bacterial and archaeal communities, respectively.

Figure  6.  Neighbor-joining tree showing phylogenetic relationships among the major archaeal operational taxonomic units (OTUs) and reference 16S rRNA gene sequences retrieved from the National Center for Biotechnology Information (NCBI) GenBank. These OTUs were those that ranked top five reads from at least a single sample. The scale bar represents the estimated number of nucleotide changes per sequence position. Percentage on nodes refer to the percentage of recovery from 1000 bootstrap resamplings. Only values larger than 50% are shown. The numbers in parentheses indicate the percentage composition of reads in each station in the following order: (SW1, SW2, SW3) (GW1, GW2, GW3, GW4).

DownLoad: Full-Size Img PowerPoint

Figure  7.  Canonical correspondence analysis (CCA) analysis showing the bacterial (a) and archaeal (b) community compositions in relation to environmental characteristics. The values of axes 1 and 2 are the percentages explained by the corresponding axis.

DownLoad: Full-Size Img PowerPoint

In the current study, CCA analysis revealed that salinity and nutrients (including nitrate, nitrite, ammonia and phosphate) were correlated with the changes of bacterial and archaeal communities (Fig. 7), and in agreement with previous studies in STEs (Ye et al., 2016; Adyasari et al., 2019; Chen et al., 2020b; Ruiz-González et al., 2022). Salinity and nutrient variations may affect microbial metabolism, thus affecting the life-sustaining activities and abundance of microbial communities (Santoro et al., 2006; Xie et al., 2018; Wang et al., 2020; Ruiz-González et al., 2022). In the bacterial analysis, although there were algal bloom outbreaks at both Stations SW1 and SW3, the bacterial communities at Station SW1 were associated with high nutrient concentrations, while a positive association was found between the bacterial communities and salinity at Station SW3 (Fig. 7a). However, salinity and nitrite (a nitrogen cycle intermediate) had a positive association with the archaeal communities at both Stations SW1 and SW3 (Fig. 7b). The results not only indicate that salinity and nutrients play critical roles in the distributions of bacterial and archaeal communities, but also further indicate that variable environmental factors may drive the heterogeneity of bacterial and archaeal communities.

3.5 Ecological implication of microbial community in STEs

SGD has been recognized as one of the most significant pathways for terrestrial dissolved solutes such as nutrients and carbon entry into coastal waters (Santos et al., 2021), but the nature and magnitude of SGD-derived solute inputs depend largely on microbial-dominated biogeochemical transformations in STEs (Ruiz-González et al., 2022). The distribution patterns of microbial communities were closely related to hydrological, geochemical and environmental characteristics in STEs (Fig. 8). Previous studies (including our study site) have shown that the significant changes of nutrients or salinity associated with SGD could result in variations in the phytoplankton biomass or ecological impairments such as hypoxia, eutrophication or algal bloom (Lecher and Mackey, 2018; Adolf et al., 2019; Chen et al., 2020a; Guo et al., 2020; Zhao et al., 2021), which are adding its share of pressure on the coastal zone and increasing pressure on these coastal resources, especially in mariculture regions (Wang et al., 2021; Yu et al., 2022).

Figure  8.  A conceptual model illustrating interaction between microbial community patterns and hydrological, geochemical and environmental characteristics in karstic subterranean esturies (STEs).

DownLoad: Full-Size Img PowerPoint

In the present work, we found that some of the bacterial sequences were related to cultured isolates (e.g., Simplicispira sp., Novosphingobium sp., Pseudorhodobacter sp., Marivivens sp. and Aquirufa sp.), and archaeal sequences were close to ammonia-oxidizing archaea (e.g., Nitrosopumilus spp. and Nitrosarchaeum sp.). These microbes could significantly and effectively participate in the removal of groundwater pollutants such as nitrate or ammonia (Table 3). Not only in karstic STEs, several previous studies have also found that microorganisms can significantly degrade pollutants in different STEs such as estuary, bay, sandy beach, submarine spring, volcanic island and marginal sea before they are discharged into coastal waters (Ye et al., 2016; Adyasari et al., 2019; Chen et al., 2019, 2020b; Jiang et al., 2020; Huang et al., 2021) (Table 3). For example, Limnohabitans, Novosphingobium, Marivivens and Nitrosopumilus are widely distributed in different STEs, which can significantly remove the nitrate, ammonia and organic pollutants in coastal groundwater. Therefore, these potential key species likely play important ecological roles in the in situ degradation of terrestrial pollutants in STEs (Fig. 8), which can be selected as excellent candidates for bioremediation of polluted groundwater, especially coastal saline groundwater.

Table  3.  A summary of the case study for in situ degradation of coastal groundwater pollutants

Study site	Country	Ecosystem type	Key bacteria	Key archaea	In situ degradation of pollutants	References
Huntington Beach	USA	beach	nitrite reductase-encoding gene fragments (nirK and nirS)	NA	denitrification	Santoro et al. (2006)
Huntington Beach	USA	beach	Betaproteobacterial ammonia-oxidizing bacteria	ammonia-oxidizing archaea	nitrogen removal	Santoro et al. (2008)
Yellow Sea	China	marginal sea	Comamonas spp.	NA	degrading aromatic compounds	Ye et al. (2016)
Northern Java	Indonesia	beach	Burkholderiaceae, Limnohabitans, Flavobacterium, Novosphingobium, Acidovorax, and Sediminibacterium	NA	performing nitrification and denitrification	Adyasari et al. (2019)
Qinzhou Bay	China	bay	Gallionella spp., Limnohabitans spp., and Novosphingobium spp.	Bathyarchaeota	utilization of dissolved organic matter; degradation of aromatic compounds	Chen et al. (2019)
Shengsi Island	China	volcanic island/beach	Erythrobacter citreu, Halomonas zhaodongensi, Marivivens donghaensi, Pelagibacterium nitratireducen, Novosphingobium naphtha, Pseudomonas pseudoalcaligene, and Nevskia ramosa	NA	utilization of nitrate and organic pollutants	Chen et al. (2020b)
Sanggou Bay	China	bay/beach	Pseudoalteromonas	Nitrosopumilusspp.	degradation of refractory sediment organic carbon and ammonia removal	Jiang et al. (2020)
Pargos	Mexico	submarine spring	Campylobacteria	NA	nitrate reduction	Huang et al. (2021)
Krka River Estuary	Croatia	estuary	Simplicispira sp., Novosphingobium sp., Pseudorhodobacter sp., Marivivens sp. Aquirufa sp.	Nitrosopumilusspp. Nitrosarchaeumsp.	nitrate and ammonia removal	this study
Note: NA represents no data.

 | Show Table

DownLoad: CSV

4. Conclusions

In this study, we provided perspective into the correlations and heterogeneity of representative microbial taxonomic groups and related ecological functions in groundwater-surface water interaction. The nutrients and salinity were the primary factors in determining the diversity and distribution of microbial communities. Taxonomic composition suggests that the key microbial communities, such as Simplicispira sp., Novosphingobium sp., Pseudorhodobacter sp., Marivivens sp., Aquirufa sp., Nitrosopumilus spp. and Nitrosarchaeum sp., can effectively remove the nitrate and ammonia at groundwater-surface water interface in the Krka River Estuary, which can provide a new perspective for coastal groundwater management such as in situ degradation/remediation of pollutants in STEs. As STEs are affected by high intensity human activities such as nutrient pollution, this study is essential to understanding the linkages between the distribution patterns of microbial communities and hydrological, geochemical and environmental characteristics in STEs.