Ecosystem Variability along the Estuarine Salinity Gradient: A Case Study of Hooghly River Estuary, West Bengal, India

Abstract

Hooghly River, a ~460 km long distributary of the Ganga River, passes through a highly industrialized Metropolis-Kolkata in West Bengal, India, and eventually empties into the Bay of Bengal at Gangasagar. To determine the patterns and drivers of planktonic community, spatiotemporal variations in water quality and micronutrient content and planktic prokaryotic and microeukaryotic abundance and diversity across the salinity gradient (0.1 to 24.6 PSU) in the Hooghly River estuary (HRE) were studied. Plankton and water samples were collected at six sites during October 2017, February 2018, and June 2018. The biotic parameters—phytoplankton (Chlorophyll a), total bacterial abundance (cfu), and copepods—were significantly higher in the downstream estuarine sites than in the upstream riparian sites; conversely, rotifer and cladoceran abundances were significantly higher at upstream stations. The most culturable bacterial strains were isolated from the two freshwater sites and one at the confluence (estuarine) and are characterized as Bacillus subtilis, Pseudomonas songnenesis, and Exiguobacterium aurantiacum. Among zooplankton, rotifers (0.09 ± 0.14 ind L⁻¹) and cladocerans (5.4 ± 8.87 ind L⁻¹) were recorded in higher abundance and negatively correlated with bacterial concentrations at upstream stations. On the temporal scale, February samples recorded lower proportions of bacterivorous zooplankton at the three upstream stations. Cluster analysis separated samples on the basis of seasons and water mass movement. The February samples showed distinct spatial characteristics, as three freshwater (FW) stations grouped together and segregated at second 2nd hierarchical level, whereas the three estuarine stations formed a separate cluster at the 50% similarity level. Samples collected in October 2017 and June 2018 exhibited mixed attributes. June samples recorded higher influence of freshwater discharge. The zooplankton abundance showed significant negative correlation with Chl a. Our results demonstrate the relative role of river continuum, land-driven lateral discharge, and seawater intrusion in shaping community structure, which needs to be considered in management and conservation planning of aquatic ecosystems, especially in highly productive and overexploited HRE.

1. Introduction

Estuaries are highly productive and dynamic semi-enclosed waterbodies linked to the sea either permanently or periodically and fed by freshwater from river inputs, resulting in a distinct salinity gradient and characteristics biota [1,2,3,4,5]. The complexity in a river estuary is determined by the variability in river water mixing with sea water, resulting in salinity [6], turbidity, and nutrient gradients [7,8,9,10]. There has been a long debate about the functioning of estuaries [11,12,13]; the community structure of phytoplankton, bacterioplankton, and zooplankton; and their relationship and their co-occurrence pattern [14,15]. The common consensus is that planktic communities play a key role in maintaining the ecological functioning of an estuary [16].

The Hooghly River is one of the most important estuarine systems in India because of the discharge from a vast river basin with substantial monsoonal precipitation (70,500 m³·s⁻¹ peak flow at Farakka), its origin from the largest montane river (~2600 km), and its long tidal zone (~280 km). Being an active tidal estuary, it has distinct biological and physico-chemical characteristics [17]. The commissioning of Farakka Barrage in 1975 facilitated the adequate quantity of Ganga water in the Bhagirathi-Hooghly River system, improving the ecosystem health and riverine-estuarine biodiversity [18,19,20,21], finally manifesting as seaward pushing the salinity zone in the estuary [22]. Fish species such as Rita rita, Sperata seenghala, Eutropiicthys vacha, Wallago attu, Clupisoma garua, Labeo calbasu, and Catla catla have made their emergence in the upper zone of the Hooghly River estuary (HRE), namely Tribeni and Banlagarh [23], which were reported from this zone before 1975, i.e., prior to the commissioning of the barrage. HRE provides valuable a nursery and recruitment habitat for commercially important species, such as Hilsa, finfish, and shrimp [24,25].

Any short-term or long-term changes are immediately reflected by the change in planktonic community [26,27], as they are self-sustaining, constituting the important components of the microbial loop while channeling carbon and energy from microbes to higher trophic levels by joining the classical food web [14,28,29,30] in aquatic ecosystems. The microbial loop explains pathway of carbon flow through nutritional food web that begins with dissolved organic matter (DOM) and reaches to the highest trophic levels bypassing some and passing through various trophic levels. The main stakeholders of the microbial loop include bacteria, zooplankton, phytoplankton, and other nutrient-cycling organisms [14,30,31]. The relative densities of bacterivorous, herbivorous, carnivorous, and omnivorous zooplankton are a reliable indicator of the functioning of the microbial loop and of ecosystem health on the spatial scale [32,33]. The zooplankton community comprises diverse feeding groups, such as bacterivores, detritivores, herbivores, and carnivores [34,35,36], forming a bridge between the microbial loop and classical food web. Information concerning co-occurrence, distribution, and community composition of the prokaryotic and eukaryotic plankton in the HRE is lacking [37,38]. The spatiotemporal variations of planktonic communities are highly affected by the hydrochemical parameters and physical forces [26,27,39,40,41]. Therefore, major components of microbial loop, i.e., bacterioplankton, phytoplankton and zooplankton, are likely to be affected by these activities. However, their co-occurrence and distributional patterns have not been studied in the HRE at a spatial scale ranging from fresh water to the estuary mouth. The knowledge of bacterioplankton–zooplankton co-occurrence is very essential, as zooplankton might act as a biotic selector for a specific microbial loop [42]. Bacterioplankton in the present study include culturable strains only because isolation of microbes is still necessary for the extraction of bioactive compounds [43], and this is accomplished by culture-based technique. Descriptions of new taxa of prokaryotes and experimental validation of microbial, ecological, and evolutionary processes are reliably based on culture based techniques. Therefore, this study isolated culturable bacterial strains and concentrated on culture-based methods.

The present study aims to elucidate the heterogeneity and co-occurrence of planktic community along the salinity gradient ranging from freshwater to the estuary mouth and ecological drivers shaping the planktic community structure in the HRE. To achieve these objectives, the study identified the ubiquitous nature and heterogeneity in distribution patterns of aquatic biological communities, including bacterioplankton, phytoplankton, and zooplankton (Rotifera, Cladocera, and Copepoda), at spatial scales during October 2017, February 2018, and June 2018 in the HRE, India. The samples were collected at six sites along the main salinity axis (0.1 to 25 PSU), from Barrackpore before the metropolis Kolkata to the estuary mouth. To elicit whether the co-occurrence or abiotic parameters are responsible for differential distribution patterns, we estimated water quality and micronutrient concentrations at all the six sites. At stations where higher correlation coefficient values for Bacteria vs. Rotifera abundance ((R = −0.76) were recorded, we isolated the bacterial strain with >50% occurrence for further characterization.

2. Materials and Methods

2.1. Study Site

The HRE is a part of Ganga River system that originates from Bhagirathi (upper stretch), flows southwards through the lower Ganga deltaic plane, and merges with the Bay of Bengal in Sundarbans as the Hooghly River in West Bengal, India (Figure 1). Kolkata city, one of the largest metropolises along the Hooghly River, having a population of about 14.5 million, utilizes the river water for drinking and domestic and industrial purposes and also discharges sewage and sludge into the river [44]. Beginning upstream of the metropolitan city of Kolkata and downstream to the confluence, six sampling sites were chosen comprising the agricultural-industrial-anthropogenic and riverine (salinity: 0.1 to 0.45 PSU)-estuarine (salinity: 4.32 to 24.6 PSU) salinity gradient along the HRE (Figure 1).

Sampling of surface water was performed at three riverine and three estuarine sites during the post-monsoon season, i.e., October 2017 and spring–February 2018, the and pre-monsoon season, i.e., June 2018, in the HRE (Figure 1). Samples were collected in cleaned polypropylene bottles from the surface to 20 cm depth for physio-chemical, bacteriological, phytoplankton, and zooplankton analysis. The details of sampling sites, abbreviations used hereafter, and common stressors at each site are provided in

Table 1. Ecological stressors, abbreviations used, and coordinates of the six stations in the HRE sampled during October 2017, February 2018, and June 2018 for the present study.

Sampling Stations (Code)	Coordinates	Altitude (ft) ASL	Salinity (PSU) Trophic Status N:P Ratio	Ecological Stresses
Barrackpore (BRK)	22.75272° N 88.36212° E (Gandhighat)	13.12	0.1–0.42 13:0.04	Industrial effluents, domestic sewage disposal, boating, bathing, occasional Immersion of idols
Dakshineshwar (DKS)	22.65643° N 88.35682° E (Dakshineshwarghat)	9.8	0.14–0.45 (18:0.092)	This site is 128 km away from the sea mouth of the river and has an estuarine condition due to significant tidal oscillation of ∼3 m. Here, the river flows through the densely populated region in Kolkata city. Mostly untreated sewage disposes into river water near Dakshineshwar ghat. The river water is also accessed for washing, bathing, and for many religious rituals
Kadamtala (KDM)	22.565° N 88.3387° E (Kadamtalaghat)	3.2	0.13–0.45 (14:0.052)	Bathing, Washing clothes, domestic effluents, ferry service, spiritual rituals, immersion of idols, oil leaching, leakage of oil from mechanized boat.
Kakdwip (KDP)	21.87208° N 88.16383° E (Harwood Point Ferry service)	0	4.32–10.64 (43:0.10)	Frequent dredging, boating, fishing, etc.
Kachuberia (KCB)	21.85903° N 88.14433° E (Kachuberiaghat Near govt. Jetty Gangasagar)	0	6.29–17.79 (60:0.18)	Frequent dredging, boating, fishing, etc.
Gangasagar (GS)	21.63307° N 88.07498° E (Gangasagar Mohana sea Beach)	0	13.09–24.6 (74:0.33)	Boating, tourist activities, dredging, fishing

.

2.2. Environmental Parameters

In total, we collected 72 surface water samples in pre-cleaned, acid-washed polythene bottles for water quality analyses during the three sampling cruises. The salinity was determined by Hanna Instruments HI98319 marine salinity tester, as salinity has influence on the demography of zooplankton [45]. For the determination of dissolved oxygen, the Winkler titration method was used. Surface water temperature was measured using a mercury glass thermometer. The flow rate was measured by mechanical flow meter (hydrobios model 438110) The dissolved nutrients (nitrate and phosphate) were estimated by colorimetric methods (following [46]) using a spectrophotometer (PerkinElmer UV/VIS Spectrometer Lambda 25, Waltham, MA, USA) after filtering the water through 0.45 µm filter paper (GF/F-Whatman, Maidstone, UK) within 12 h of sampling, and the filtered waters were stored in 100 mL pre-cleaned, acid-washed polythene bottles at 4 °C. All the parameters were analyzed following the standard procedures for water sampling and examination of water quality [47].

2.3. Bacteriological Analysis

2.3.1. Enumeration, Isolation, and Characterization of Culturable Bacterial Strains

To study the bioactive potential for further prospecting, we estimated culturable bacterial concentrations and isolated the most culturable strains from three different sites. Five replicates of the water samples were collected in sterile polypropylene bottles from the surface to 15–20 cm deep and transported to the laboratory at 4 °C in an icebox (Table 1; Figure 1) during October 2017, February 2018, and June 2018. In the laboratory, water samples were stored at −21 °C until further processing. In the laboratory, bacterial concentrations were estimated by direct plate count method following [48]. The surface water samples, collected from different sites along an anthropic gradient in the HRE, were spread on media plates. Total bacterial density (colony forming unit: CFU mL⁻¹) was enumerated on nutrient agar plate, which was incubated at 37 °C for 24 ± 1.5 h [48]. The bacterial counts obtained were used to estimate the number of bacteria grown on the media plates used for DNA extraction. Colonies with different morphologies were subcultured into pure cultures by inoculating them into freshly prepared agar plates [49]. The most abundant colony at each of the three sites was recorded for subsequent identification and statistics. At two sites, BRK and DKS were recorded for relatively higher abundance of bacterivorous zooplankton and higher strength of association between bacterioplankton and zooplankton. With an aim to identify highly abundant bacterial colony at these three sites (BRK, DKS, and GS), we isolated the bacterial colony with 50–70% occurrence and further cultured for sequencing and phylogenetic analysis [48,50].

2.3.2. Clustering, Alignment, and Phylogenetic Analysis of 16S rRNA Gene Fragments in Most Culturable BRK2, DKS, and GS1 Bacterial Strains

DNA from the bacterial culture—BRK2, DKS, and GS1 strains—was isolated using the bacterial gDNA isolation kit (XcelGen, Gujarat, India). Isolation of DNA was carried out according to manufacturer’s instructions. First, 1.2% agarose gel was used to evaluate the quality of isolated DNA, and a single band of high-molecular band of the PCR amplicon was detected (Figure S1). Amplification of isolated DNA was performed with 16S rRNA-specific primer (8F and 1492R) using Veriti^® 96 well thermal cycler (Model No. 9902, Thermo Fisher Scientific, Waltham, MA, USA). Sanger sequencing was performed using BDT v3.1 Cycle sequencing kit with M13F and M13R primers on ABI 3730xl Genetic Analyzer was performed after the PCR amplicon was enzymatically purified. A consensus sequence of 1284, 1465 bp, and 1487 of 16S rRNA was generated by using aligner software from forward- and reverse-sequence data. The consensus sequence of all the three strains, accession numbers, and origin are shown as Table S1. All nucleotide sequences were deposited at the National Center for Biotechnology Information (NCBI) strain library with accession numbers provided in (Table S1).

The evolutionary history was inferred by using the maximum likelihood method and Tamura–Nei model [51]. The tree with the highest log likelihood (−10,338.41, −10,518.13, −10,395.19) is shown, respectively, for BRK, DKS, and GS strain. Initial tree(s) for the heuristic search were obtained automatically by applying neighbor-joining and BioNJ algorithms to a matrix of pairwise distances estimated using the Tamura–Nei model and then selecting the topology with superior log likelihood value. This analysis involved 61 nucleotide sequences. There were a total of 1284, 1465, and 1484 positions in the final dataset of each strain (BRK, DKS, and GS). Evolutionary analyses were conducted in MEGA11 [52].

2.4. Phytoplankton Analysis

For qualitative analyses, five replicates of 1 L surface water samples were preserved in neutral Lugol’s solution (1% Lugol’s solution and 4% formalin) and brought to the laboratory for species identification. In the laboratory, samples were concentrated 10 times by centrifugation, and algal cells were observed on a Sedgewick rafter cell under a compound microscope (10×–400× magnification), and abundant species were identified using the standard key [53,54,55,56,57] and AlgaeBase (www.algaebase.org (accessed on 12 July 2018)) [56]. The amount of primary productivity was estimated in terms of chlorophyll a (Chl a). Chlorophyll pigment was analyzed through extraction using a mixture of dimethyl sulfoxide and 90% acetone [58] and enumerated by spectrometry using a Turner TD-700 fluorimeter (New York, NY, USA) following the standard method [46,59].

2.5. Zooplankton Analyses

Zooplankton samples were collected by making surface tows (0–20 cm) with a customized plankton net with 53 μm mesh size, 0.25 m mouth diameter, and preserved in 4% (w/v) buffered formalin immediately after collection in a 100 mL transparent bottle. At each site, 100 L water was filtered in five replicates. In the laboratory, the plastic bottles containing preserved zooplankton were thoroughly mixed, and a 1 mL subsample was drawn with a fine pipette to a Sedgwick–Rafter plankton counting cell for enumeration under the compound microscope (model no: Olympus CX21LED, Bartlett, TN, USA). The numbers per liter of each genus was quantified and calculated using the following formula:

$$\mathrm{N}=\frac{\mathrm{A}\times \mathrm{C}}{\mathrm{L}}$$

where N denotes the number of plankton per liter, A is the average number of plankton in all counts, C is the volume of original concentrate in ml, and L is the volume of original water filtered, expressed in liters. Zooplankton species were identified to their lowest possible taxon [60,61,62,63,64,65,66].

Trophic-Based Zooplankton Community Analysis

Based on the published literature (Supplementary Table S5) and our own observations on propensity of feeding [35,67,68,69,70,71], the zooplankton communities identified at each station were characterized on the basis of functional feeding mode. Different fractions of zooplankton representing bacterivorous, herbivorous, carnivorous, and omnivorous types were segregated following standard literature [34,35,36,69,71,72,73,74,75,76,77,78] and analyzed separately.

2.6. Data Analysis

To elicit variations at spatio-temporal scales, the similarities of community composition among the sampling stations for each sampling date and also among sampling dates were compared. We first determined the centroid vector that represents the average composition of the group/species. Spatial heterogeneity was estimated using the mean and standard deviation of the similarities from the estimated similarity vector. We calculated the Bray–Curtis index to characterize the dissimilarities between samples (β-diversity). Square-root-transformed species abundance data were used for constructing the Bray–Curtis matrix of dissimilarity with average linkages group classification [79]. As the Bray–Curtis similarity mixes the differences due to species losses and species turnover, we also partitioned this index to understand both components of dissimilarity. For abiotic parameters, the distance between two samples was measured by Euclidean distance (ED), as ED is more appropriate for a low-dimensional data set [79].

To characterize the zooplankton diversity present in each sample (α-diversity), we calculated the Shannon diversity index (H′) (H′ = −Sum(Pi × Log(Pi)) (Shannon, 1948), evenness index (J′) (J′ = H′/Log(S)) [80], and species richness (d) (d = S − 1/Log(N)) [81]. To determine the variations among samples, non-metric multidimensional scaling (NMDS) ordination was computed based on Euclidean distance [82]. To identify the drivers of species abundance, pairwise correlation of water quality, and biotic parameters, the degree of a linear association between any two of the parameters was measured using Pearson’s correlation coefficient (R). To test the distribution of data, the Shapiro–Wilk test of normality was used, and outliers were detected using scattering plot prior to Pearson’s correlation analysis. Highly correlated parameters that may influence the community structure were identified. Indexes of dissimilarity, Shannon’s index (α-diversity), Pielou’s index (evenness), and Euclidean distance were calculated with PRIMER-version 6.0.

3. Results

3.1. Spatio-Temporal Patterns

All the estimated abiotic (Figure 2) and biotic (Figure 3) parameters except DO level, rotifer (Figure 3B), and cladoceran abundance, showed significant seaward increase; in contrast, the rotifer abundance showed a significantly seaward decreasing trend (R² = 0.6; p < 0.0001).

3.2. Physicochemical Parameters

The surface water temperature ranged from 26–29 °C with an average of 27 °C and recorded a significant (R² = 0.3; p < 0.02) seaward increase (Figure 2A). The highest average concentration of Ca⁺⁺, Na⁺, and K⁺ was found at the mouth (Figure 2A–C; Table S2). The trophic level-related parameters (nitrate, phosphate, bacterial concentration, and Chl a) showed higher values in the estuarine stations (Table S2; Figure 2J,K); however, the mean Chl a level was the highest at KCB station preceding the GS (Figure 3C). The dissolved oxygen concentration ranged from 6 to 8.2 mg·L⁻¹ but had disorderly spatial variation. The flow rate of the Hooghly River was recorded 0.3 m·s⁻¹ in February and 1 m·s⁻¹ in August.

3.3. Phytoplankton

The phytoplankton species recorded in the sampling stretch were Pediastrum, Spirogyra, Coscinodiscus, Cyclotella, Melosira, Ankistrodesmus, Aulacoseira, Coelastrum, Microcystis, Oscillatoria, Anabaena, Aphanocapsa, Coscinodiscus radiatus, Pleurosigma formosum, Coscinodiscus lineatus, Biddulphia sinensis, and Chaetoceros lorenzianus (Table S3). The mean Chl a concentration varied from 29.1 mg L⁻¹ to 219.9 mg L⁻¹, showing significant increase towards the river plume; the mean Chl a concentration was the highest at KCB station preceding the confluence GS (Figure 3C). The Chl a values showed positive correlation with all the abiotic parameters; however, a significant positive correlation was recorded with nitrate (R = 0.79) and Ca⁺⁺ ion concentration (R = 0.76). With biotic components, the significant negative correlation was recorded between Chl a values and zooplankton abundance (Table S6).

3.4. Zooplankton Community Structure

The symmetric map of all the parameters estimated, in rows and columns in principal coordinates, is given in Figure 4, in which the response category points to separate stations on an ordinal scale. Looking at the spatial scale with respect to the horizontal principal axis, all the zooplankton community at all the riverine stations aggregated together on the right side, whereas the last station at the confluence was set aside from other stations and positioned on the left, showing higher variation among sampling seasons (Figure 4).

The right extreme position of all the riverine samples and their unique position in the middle of the second (vertical) axis indicates that the responses are more in the intermediate categories of the scale rather than a mixture of extreme responses at the temporal scale (Figure 4). The scale values of these samples optimally discriminate between the 18 samples (6 stations × 3 seasons), giving maximum between-sample variance. The second dimension then separates out samples on the basis of seasons, and all the estuarine samples collected in February are polarized towards the top. Upstream estuarine stations (KCB and KDP) aggregated inside, where both extremes of the response scale as well as the missing response are located (Figure 4). As a result, samples were arranged in the ordination of downstream confluence to freshwater stations from left to right. All KDP samples and KCB October samples are in a unique position inside, with relatively high polarization of responses and high missing values. Their position in the middle of both axes reveals more responses in the scale’s intermediate categories rather than a mixture of extreme responses at the spatiotemporal dimension. Figure 4 depicts the principal inertias at the positive ends of each axis, which were measured by adding together the percentages of inertia, i.e., 63.1% + 13.3% = 76.4%. This shows a “residual” of 23.6%, which is not shown in the map.

The unique right-side positioning of all the riverine stations may be attributed to the presence of rotifers (0.09 ± 0.14) and cladocerans (5.4 ± 8.87) in dL⁻¹, captured at upstream freshwater stations only, whereas copepods were present at all stations (Figure 5A). Integrating all zooplankton samples (Figure 5A; Table S4) were dominated by the Copepoda (92%) followed by Cladocera (7%) and Rotifera (1%). Total zooplankton abundance was more affected by seasons and showed a disorderly distribution among stations. Zooplankton density showed a peak in February (beginning of spring) at all sampling stations except KDM, where the peak was recorded in June samples (Figure 5A). The peak of the total zooplankton abundance was mainly contributed by the copepods at all stations (Figure 5A). At the three riverine stations, the rotifer densities were significantly lower during the peak of the total zooplankton abundance (Figure 5A). The indices of diversity, richness, and evenness of zooplankton recorded at six selected stations in the HRE are provided in (Figure 5B–D).

The heterogeneity at the spatial scale was estimated by the Bray–Curtis similarity matrix, which was observed to be higher than the heterogeneity at the temporal scale (Figure 6). It shows that the average similarity between all sampling dates was higher than the average similarity between the stations on a single sampling date. The cluster analysis gives further insights into the relative role of seawater and freshwater in shaping the community structure and segregate stations accordingly (Figure 6). The first hierarchical level separates the June samples of DKS and KDM and February samples of the two uppermost stations (BRK and DKS) from the remaining samples at 90% dissimilarity, from which June samples of DKS were separated at the second hierarchical level at 75% dissimilarity. June samples of DKS mainly represented the Cyclopoida adults, copepodites, and the Cladocera Moina macrocopa, whereas February samples of upper two stations and June samples of KDM grouped together and represented copepod nauplii, cyclopoids, and calanoids. The June samples of BRK, KDP, and KCB grouped together and separated at the third hierarchical level, whereas the February samples of all the estuarine and lowermost riverine (KDM) stations grouped together and separated at the fourth hierarchical level. Samples collected in February 2018 at the three estuarine stations and the lower most riverine (KDM) station grouped together and separated from remaining samples at the hierarchical level V (VIIB). The highly indicative zooplankton of cluster VIIIB are copepod nauplii and harpacticoids (Figure 6). This cluster clearly indicates further upward intrusion of the marine community in February. All the October samples aggregated at intermediate position and showed clear separation of riverine and estuarine stations (Figure 6). October samples showed distinct spatial variations, where riverine samples were separated from estuarine samples at the fifth hierarchical level, and the June samples of GS (estuarine mouth) joined this cluster (Figure 6).

The trophic-based structuring of the zooplankton community revealed that the community was dominated by the omnivores followed by the herbivores and the bacterivores, respectively (Figure 7A–C). Variations in trophic-based community structure were more prominent at the spatial scale than the temporal scale. Bacterivorous and detritivorous species were recorded at upper stations and limited to KDP, whereas downstream stations mainly represented the omnivorous copepods (Figure 7A–C). Detritivores were mainly represented by the bdelloid rotifers at BRK and DKS. Differences in bacterivorous fraction of zooplankton community were not significant among the three riverine stations (Figure 7A–C); however, at temporal scale, February samples recorded significantly lower fractions of bacterivorous zooplankton at all the three riverine stations.

3.5. Bacteriological Analyses

Total culturable bacterial density CFU mL⁻¹ varied from 0.06 to 300 × 10⁵ CFU mL⁻¹ at six selected sampling stations during October 2017, February 2018, and June 2018 in HRE. With the lowest bacterial concentration at Barrackpore and the highest at Gangasagar, unique spatial differences and significant seaward (R = 0.66 p < 0.0028) increase in density of culturable bacteria were recorded in the present study (Figure 8).

Overall bacterial densities were higher at downstream stations (KDP, KCB, and GS) than the upper freshwater stations (BRK, DKS, and KDM). Amongst upstream stations, bacterial densities were negatively correlated with total rotifer density and Chl a concentration; however, the correlation was significant for Bacteria vs. Rotifera abundance only (Figure 9; Table S6, p < 0.01). The most abundant culturable strains were Bacillus subtilis, Pseudomonas songnensis, and Exiguobacterium aurantiacum, respectively, at BRK, DKS, and GS (Figure 10).

3.6. Characterization of Highly Abundant Bacterial Taxa at BRK, DKS, and GS

The branch lengths and topology of a phylogenetic tree of the isolated strains were attained by the maximum likelihood technique. The phylogenetic trees were formed based on evolutionary distance data of 16S rRNA gene sequences. The tree shows the integrated relationship of the isolated strains with most similar Gangetic bacterial species genes based on highly similarity 95–100% taken from the 12 different Gangetic samples generated by illumine and Nanopore sequencing platforms Figure 10A–C.

3.7. Interrelationship among Taxa

The pairwise multiple correlations among means are summarized in Table S4. The total zooplankton abundance showed negative correlations with all other parameters, recording maximum strength of association with Chl a, but the values were not statistically significant. Zooplankton taxon-specific correlation with Chl a gives insight, explaining the group-wise differential relationship: Rotifera (R = −0.67; p < 0.01) and Cladocera (R = −0.43; p < 0.05) showed a significant negative correlation with Chl a, whereas Copepoda (R = −0.4; p < 0.1) did not correlate significantly with Chl a. Among freshwater stations, Cladocera and Copepoda exhibited significant positive correlations in all the three sampling cruises. Rotifera concentration correlated negatively with bacterial concentration, pH, EC, salinity, nitrate, phosphate, and Ca⁺⁺ levels (Table S6). Among biotic components, Chl a values were in a core position, showing negative correlation separately with all the three zooplankton groups, whereas the bacterial concentration showed significant negative correlation with rotifer fraction of the zooplankton only (Table S6).

4. Discussion

The observed spatio-temporal variations in zooplankton, phytoplankton, and bacterioplankton concentrations and the interrelationship among taxa and spatial occurrence of the most abundant culturable bacterial strains in the HRE along the salinity level ranging from 0.01 to 25 (PSU) unequivocally confirm the complex nature and dynamicity of the estuary. The local complexity and seaward gradient of nutrient concentrations, salinity, turbidity, and Chl a [7,41,83] have been demonstrated in various estuarine ecosystems globally. The strength of association between zooplankton (separately with Rotifera, Cladocera, and Copepoda) with bacteria and Chl a in the present study provide additional insights into relative role of bacterial and algal carbon in supporting zooplankton community in highly dynamic, tropical estuary such as the HRE.

4.1. Interrelationship among Taxa

At upper stations, the rotifer community, recording lower abundance with higher diversity, exhibited significant negative correlation with bacterial concentrations, albeit bacteria and rotifer showed an overlapping relationship with abiotic parameters (e.g., pH, EC, salinity, TDS, nitrate, phosphate, and Ca⁺⁺). Similar physical requirements though significant negative correlation between Bacteria and Rotifera reflect strong grazing pressure by Rotifera on prokaryotes. On the other hand, the negative correlations of Chl a concentration with total zooplankton abundance in general and separately with all the zooplankton groups in particular suggest that zooplankton placed major grazing pressure on the phytoplankton community in the whole HRE. The prokaryotic community is controlled by the smaller zooplankton, particularly Rotifera, at the upstream stations. The majority of rotifers are filter-feeding and nanophagous which substantially utilize smaller organisms of the microbial web (bacteria and nano-flagellates and -ciliates) and are able to thrive on bacteria in eutrophic waters [32,84]. Earlier studies have reported the key role of autotrophic protists in shaping the zooplankton community structure by forming the base of food chain in the estuarine ecosystem [85,86].

The observed gradient in taxonomic diversity of zooplankton along the river to the sea continuum has also been reported in other estuaries [1,41,87]. The differences between the bacterial community in freshwater and estuarine systems can be explained by the fact that the freshwater microbial community are directly under grazing pressure by the smaller zooplankton such as Rotifera, whereas at and around the confluence, microbial abundance is mainly controlled by the grazing pressure from developing stages of larger zooplankton or trophic cascade effect, where bacterivorous flagellates are grazed by Copepoda [85,88,89]. In line with this, previous studies have related microbial growth rates with rotifer grazing effects [90,91,92,93]. It may be noted that omnivores are a major controlling factor in estuarine stations, and this higher abundance of omnivorous zooplankton in downstream saline stations offsets the trophic cascade effects, resulting in higher bacterial abundance. The weak association of copepods with Chl a indicate bottom-up functioning of trophic cascade in estuaries. Omnivore-driven dampening of trophic cascades has been reported in other estuaries also [89,94].

4.2. Ecological Drivers of Planktonic Community Structure

The dissimilarity between sites (β-diversity) may reflect two different ecological situations: (i) It could reflect the loss of species from one site compared to the other. In this situation, the species present at one site are nested inside a larger set of species present at the other site. Alternatively, the following site records none of the species recorded at preceding site. In the latter case, there is a complete turnover of species between the two sites. In the HRE, both events occur concurrently, and the Bray–Curtis index takes into account both the components of dissimilarity. Further, in a distance-based multivariate model, variability among 72 samples based on riverine and estuarine locations was consistent with the main drivers, with salinity itself explaining 39.6%, 40.9%, and 37.4% variability, respectively, for rotifers, bacteria, and Chl a concentration. The freshwater-seawater salinity range is not continuous but rather subdivided in distinct stages, which is generally manifested in planktonic abundance [95]. The extensively used salinity classification is the Venice system that is 0–0.5, 0.5–5, 5–18, 18–30, and 30–40 (PSU). Further, based on salinity-range data, ref. [96] explained five salinity zones including 0–4, 2–14, 11–18, 16–27, and 24-marine. The present results based on explicit criteria of planktonic abundance along with the salinity range support the salinity classification by [96].

Surface water temperature increased downward from 26 to 29 °C with the average temperature value of 27.3 °C and showed significant positive correlation with bacteria and Chl a concentration. The role of temperature in maintaining planktonic community and invasion of alien species in estuary has been emphasized earlier [95]. Nutrient-related parameters (nitrate and phosphate) and TDS explain the higher abundance of Rotifera in upstream stations, where 56% of the factors determining Chl a concentration is the nitrate level and 52% of the factors contributing higher bacterial concentration is explained by the phosphate levels in the surface water. These microbes can use nitrate reductase, nitrite reductase, glutamine synthetase, and other compounds for nitrogen assimilation [97,98].

Ionic concentrations also show a vital role in unravelling the zooplankton, phytoplankton, and bacterioplankton community structure [99,100,101]. Hence, the discernible spatial distribution of most abundant bacterial strains and higher bacterial concentration at downstream stations can be explained by differential sodium and calcium requirements in bacterioplankton and phytoplankton. Freshwater and halophilic bacteria have overlapping physiological attributes but different sodium requirements. The marine strains require higher sodium and calcium level to grow, while freshwater and terrestrial strains such as E. coli can multiply at a higher rate without sodium [102]. In an estuary, the riverine and marine bacterial species, having ecologically similar physiological abilities but different sodium requirements, favor the halophilic strains (e.g., Exiguobacterium aurantiacum in present study) when freshwater from rivers enter the sea, as the sodium and calcium dependence does not constitute a fundamental ecological difference; rather, it only regulates the locally adapted strains responsible for them [42].

4.3. Spatio-Temporal Pattern

The lowermost riverine station is KDM, which grouped with upper freshwater stations during June 2018 and, in contrast, grouped together with lower estuarine stations during February 2018 (Figure 6). On the other hand, the clustering of June samples of upper two estuarine stations with freshwater stations is suggestive of a monsoon-based regime shift in the Hooghly estuary. A higher volume of freshwater discharge owing to pre-monsoon rainfall results in seaward extension of riverine biota, whereas in the case of lower discharge and reduced river flow during February, it leads to upward extension of marine biota into a lower freshwater station (KDM). The association of October samples of the KDM station was with the estuarine stations and February and June samples with upstream freshwater stations. This gives an insight in the seawater intrusion process in the HRE. This differential grouping of lower freshwater station with saline stations and with freshwater stations depending upon the season highlights the importance of stratification based on salinity in the suppression of turbulent vertical mixing in the estuary. In February and June, weak stratification and strong vertical mixing prevails during this period, and the riverine discharge counters the seawater intrusion. During the monsoon, the strength of the estuary circulation increases as river discharge rises, while the length of seawater intrusion diminishes. However, in October, during post-monsoon season, the suppression of turbulent mixing and the strength of the estuarine circulation is mostly determined by tidal velocity. Seawater incursion associated with estuary circulation reduces as tidal velocity rises but increases when river discharge rises [103]. In fact, Monsoon flows affect all facets of estuarine hydrobiology and community structure. In the light of monsoonal influence, the HRE may be called a “tropical monsoonal estuary” [104]. The monsoonal precipitation-driven rapid decline of salinity and surface runoffs limits the distribution of marine forms, whereas the intermediate condition favors rapid multiplication of the brackish water forms and re-assemblage of the halotolerant groups, thereby resulting in the shifting regimes of transitional stations. This explains the seasonal shifting of the KDM station from freshwater to salt water, as observed in the present study. Other studies have shown the relevance of inshore water zones as zooplankton sources in large rivers such as the Danube [105,106] and St. Lawrence rivers [107]. Estuarine ecosystems in India show two peak periods, and peak time varies from region to region [108]. The present observation corroborates the previous results, in which two peaks were recorded in different months [18,109,110]. The spatial attributes of rivers are not always continuous [111], as proposed by [112]. Local land-driven discharge, changes in drought and flood regimes, and the establishment of diverse hydrological retention zones [113,114] (due to silt deposition) alter flow and river beds differentially during dry and wet seasons [113,115,116]. Therefore, results suggest that the phytoplankton, zooplankton, and bacterioplankton dynamics in the HRE are controlled by the interplay of hydrological regime, nutrient concentrations, and allochthonous inputs.

4.4. Spatial Occurrence of Highly Abundant Bacterial Strains

Bacterial study includes culture-based evaluation because studying morphology, ecology, bioprospecting for human use, and bioremediation for the purpose of culture-based isolation is practiced globally [117]. Consequently, while many environmental studies focus on large descriptions of microbial diversity through whole genome sequencing (WGS) approaches, other environmental studies rely on culturing approaches to estimate the abundance of a given culturable taxa in the environment. Both approaches are usually performed independently, and this does not allow direct comparison of the benefits and constraints of both methods. However, the culture-based approaches select only a subset of culturable bacteria, and it remains unclear to what extent culture-driven enrichments could be compared between different environmental samples. The present study depends on culture-based analysis, in which the upstream stations BRK and DKS recorded Bacillus subtilis and Pseudomonas songnenesis strains, and the lowermost station GS recorded Exiguobacterium aurantiacum as the highly abundant bacterial strain. The presence of soil bacteria in upstream stations indicate the influence of land driven allochthonous discharge as the major contributory factor of bacterial abundance, which enters the aquatic food chain by the microbial degradation and remineralization and its recycling within the pelagic food web. Bacteria are the only organisms capable of recycling DOM, making them an essential component of ecosystem functioning [85]. The differences in bacterial populations could be due in part to differences we observed in DOM and DOC quantity and quality among the three samples, as it has previously been shown that DOM and DOC strongly structure bacterial communities in aquatic environments [118,119,120].

The differences in bacterial abundance can also be attributed to differential abundance and community composition of zooplankton, as they utilize bacteria directly or indirectly through many other bacterivorous organisms and establish the link to the traditional aquatic food web [121]. Consequently, ecological productivity in estuaries is affected by both top-down and bottom-up mechanisms. Top-down controls, such as meso-zooplankton grazing, may decrease micro-zooplankton populations, enabling phytoplankton species to bloom and altering the overall structure of the microbial community [87,121,122,123]. At riverine sites around Kolkata city, strong top-down effects are major regulators of bacterioplankton abundance by bacterivorous organisms [124]. Additionally, the lower rotifer abundance during the peak of copepod abundance at upstream freshwater stations is also suggestive of top-down control of rotifers by copepods [41,67,125,126]. In contrast, the dominance of omnivorous and herbivorous fractions of zooplankton and higher Chl a and bacterial concentrations near the estuary mouth reflect strong bottom-up impacts, where the nutrient-loaded environment favors microbial and phytoplankton growth that supports omnivorous species and all the trophic guilds in the absence of a distinct trophic cascade [89,93,127].

The bacterial abundance and physicochemical parameters, particularly nutrient concentrations and differential abundance of rotifers and cladocerans, are indicative of land-driven allochthonous influence from urban discharge of the metropolis city of Kolkata. The highest Chl a and Na⁺ at KCB preceding the confluence (GS) attest to the established facts that estuaries, particularly the mixing zone, are the most productive ecosystem and constitute an important system that provides valuable nursery and recruitment ground for commercially important species. However, with the development of sequencing technologies, further study is needed to elucidate the potential novel functions and phylogenetic relationships by sequencing the genomes of entire communities to understand relative contribution of bacterial community in HRE.

The instant change in community structure recorded in this study at lower salinity levels suggests the potential for underlying change in the oligohaline or limnetic stretches. In line with successful management of the San Francisco estuary based on isohaline condition [128], the present results indicate management options for the HRE, recommended as limited withdrawals to fixed fraction of total river flow beyond a minimum flow threshold [129]. The concerned administration needs to fix a minimum flow target in accordance with ideal region-specific isohalines in the estuary.

5. Conclusions

The zooplankton, phytoplankton, bacterioplankton, and abiotic parameters reported in this paper elucidate the patterns and drivers of differential community structures across the salinity gradient in the HRE. Among zooplankton, rotifers and cladocerans are numerically dominant and exert strong selection pressure on bacterial community and clear the suspended particles from the water column at upstream stations, whereas copepods play a major role in structuring microbial community at downstream estuarine stations. The negative correlation between Chl a and bacterial abundance, though insignificant, points to the competition for inorganic nutrients between phytoplankton and bacteria. Spatial variations in the trophic-based zooplankton community structure also suggest differential effects of direct bacterivore behavior by rotifers and omnivore-driven, suppressed trophic cascade effects through copepods, both of which concurrently play an important role in shaping the HRE community. The abiotic parameters such as surface water temperature; elemental concentrations of Ca⁺⁺, Na⁺, and K⁺; and the trophic level-related parameters (nitrate and phosphate) record significant seaward increase, which in turn reflects the increased concentrations of bacteria and Chl a (primary production) at downstream estuarine stations.

The three isolated strains of the most culturable bacteria at Barackpore, Dakshineshwar, and Gangasagar, characterized as Bacillus subtilis, Pseudomonas songnenesis, and Exiguobacterium aurantiacum, indicate differential influences of land-driven discharge and spatial heterogeneity in the prokaryotic community structure. The observed alteration in planktonic community structure in the sampled stretch of the HRE points to larger impacts of water extraction and sewage discharge on salinity level, resulting in changes in the riverine community in the sampled limnetic-to-oligohaline stretches.

At the temporal scale, the increased river discharge during pre-monsoon and monsoon season plays an important role in shaping the community structure by upward extension of marine influence during the waning season but downward extension of river influence during the monsoon. Therefore, the complexity of phytoplankton, mesozoopalnkton, and prokaryote communities responding to variable elemental and nutrient concentrations is driven by the differential mixing of freshwater and marine sources. Both bottom-up and top-down effects play a vital role in shaping the community in the HRE.

Increased urbanization with uncontrolled water extraction, discharge of industrial and domestic wastes in coastal waters near the mouth of the Hooghly River, and shoreline development affect the planktic community, consequentially affecting overall ecosystem health.

Furthermore, this study also highlights the role of land discharge at freshwater stations, season-specific seawater intrusion in the river, and abiotic variables including trophic status as drivers of abundance of the prokaryotic and eukaryotic planktonic community. These results suggest limited wastewater discharge and water withdrawals to a fixed fraction of total river flow beyond a minimum flow threshold maintaining isohaline. The concerned administration needs to fix a minimum flow target in accordance with ideal region-specific isohalines in the estuary. Therefore, for any future planning, the volume of water withdrawals and wastewater discharge, monsoon-driven regime shift, and internal trophic-based regulation mechanisms need to be considered.