Coastal Hypoxia Drives Microbial Diversity: Elucidation through 16S rRNA Amplicon Sequencing

The formation of oxygen-depleted zones in the bottom waters is one of the most widespread phenomena in coastal areas. Upwelling episodes occurring along the west coast of India due to the southwest monsoon lead to an increase in biological productivity which further lowers the dissolved oxygen in the upwelled waters, which intensies annually between June and October. Here, we have determined the changes in the microbial community in response to the varying oxygen levels and other physicochemical parameters at the Candolim Time Series Station using high-throughput sequencing. Amplicon Sequence Variants across all the samples collected in different seasons were mostly aliated to the phyla Proteobacteria, Actinobacteria, Bacteroidetes, Verrucomicrobia, Chloroexi, Firmicutes and Planctomycetes, with the most dominant being Proteobacteria (21-41%). Statistical analysis revealed that microbial diversity differed signicantly with changing DO, ammonia, nitrate and nitrite concentrations during different seasons. The microbial community shift due to seasonal hypoxia results in the differential biogeochemical cycling of essential nutrients with certain years seeing redox conditions up to sulphate reduction, while certain years seeing only nitrogen loss. Future scenario of global warming will serve as a big challenge for understanding the role of microbial diversity and its implications in the cycling of natural elements.


Introduction
Coastal ecosystems play a pivotal role in regulating the global biogeochemical cycles. Some of the major processes in these ecosystems include the exchange of nutrients, carbon in the form of dissolved organic matter, reactive organic and inorganic trace gases, and various other physical and biogeochemical elements. Coastal ecosystems are signi cant as they are the areas with high productivity [Wiggert et al. 2005]. The increased biological productivity due to monsoonal upwelling leads to increased organic matter production, utilised by heterotrophic microorganisms. The intense decomposition of organic matter leads to rapid utilisation of in situ oxygen, leading to its depletion and forming hypoxic conditions in sub-surface waters [Gomes et la. 2019;Naqvi et al. 2000].
The west coast of India experiences coastal hypoxia due to upwelling, which is mainly in uenced by the southward movement of the West Indian coastal current during the South West Monsoon (SWM) [Madhupratap et al. 1996]. The upwelled waters are already hypoxic (<2 mL L −1 ) as these are from the thermocline region off the continental shelf of India [Naqvi et al. 2000]. During the SWM, the bottom waters having low oxygen levels profoundly impact the biogeochemistry and ecological functioning of the biome and the underlying sediments. Hypoxia can also be caused by the anthropogenic inputs of nutrients and organic matter, leading to a drastic decrease in DO due to eutrophication [Rabalias et al. 2002]. Such low oxygen levels threaten marine ecosystems and are also responsible for the loss of benthic microbial habitats. In addition, prolonged hypoxic conditions may change the microbial community [Crump et al. 2007] in response to oxygen depletion. In a case study carried out in the Black Sea, the presence or absence of oxygen in the water column has been suggested as the main reason for variations in microbial community structure and function [Thamdrup et al. 2000]. Beman and Carolan [2013] reported a non-linear relationship between the DO concentration and bacterial richness in the oxygen-de cient waters of the Eastern Tropical North Paci c Ocean.
Hypoxic events are often considered as one of the major global environmental problems [Diaz and Rosenberg 2008;Naqvi et al. 2010;Breitburg et al. 2018]. The seasonal coastal hypoxic zone of the Eastern Arabian Sea is by far the largest of all the known coastal hypoxic systems formed due to natural and anthropogenic activities [Naqvi et al. 2006;Sudheesh et al. 2020]. It is observed that if hypoxic conditions persist for a long time, then the organic matter and nutrients tend to accumulate in the sediments resulting in the expansion of the hypoxic zone. This leads to a further decrease in DO levels, thus establishing anoxic conditions, favouring the release of H 2 S by microorganisms [Diaz and Rosenberg 2008]. The progression of seasonal coastal hypoxia has been shown to follow a predictable pattern that starts with SWM's advent (June to September). Oxygen depletion at times might also begin in April-May and gradually intensify with time [Naqvi et al. 2006]. The CaTS-G5 station located on the West coast of India is being extensively monitored since 1997 to study the changes occurring due to coastal hypoxia. This could be considered a representative station to understand the oxygen de ciency in the coastal regions [Naqvi et al. 2006]. The oxygen-de cient conditions intensify towards the end of the southwest monsoon, whereas the water column remains well-oxygenated during other periods.
Prolonged oxygen-de cient episodes often change the bacterial community composition [Crump et al. 2007], affecting the ecosystem's functioning [Reed and Martiny 2013]. As oxygen is the most favourable electron acceptor in aerobic aquatic systems, its availability could in uence gene expression of the bacterial community and affect respiratory metabolism. A shift in bacterial community composition during the development of hypoxia was reported by Beman and Carolan [2013] in the oxygen minimum zones of the Eastern North Paci c Ocean, where they recorded a non-linear relationship between bacterial richness and DO concentration. On the contrary, a negative correlation between bacterial richness and DO was reported by Spietz et al. [2015] in a hypoxic estuary. Various bacterial lineages belonging to the classes Alphaproteobacteria, Gammaproteobacteria and Deltaproteobacteria have been reported to be abundant in oxygen-de cient marine waters. Microbes reported from oxygen-depleted waters falling into these lineages are considered to be chemolithoautotrophic sulphur oxidisers that play a key role in sulfur and carbon cycling in the ocean ecosystem [Lipsewers et al. 2017]. Previous studies, carried out in the Northern Indian Ocean oxygen minimum zone (OMZ), based on metagenomic analysis, reported taxa belonging to Proteobacteria, Actinobacteria, Nitrospinota, Planctomycetacia and SAR406 clade; however, the distribution varied along The present study was conducted to investigate the effect of hypoxia on the microbial diversity at the CaTS-G5 station in the Arabian Sea. The microbial communities were deduced by 16S rRNA amplicon sequencing from water samples collected during different seasons, i.e. pre-monsoon, monsoon and post-monsoon. The annual variation in the bacterial diversity was also studied, along with the distribution of bacterial communities in response to the physicochemical environmental variables and predictive functional analysis of the microbial communities.

Materials And Methods
Sampling locations and sampling methodology Onboard the ship the temperature data of the seawater samples were obtained using a Conductivity-Temperature-Depth (Seabird Electronics SBE911) rosette system tted Niskin bottles. During the eld trips, it was obtained using a portable CTD (SeaBird Scienti c SBE25 plus V2), and water samples were collected using 5 litre Niskin samplers. The DO samples were analysed following the Carpenter modi cation of the Winkler method [Carpenter 1965]. The titration was automated using ODF PC software compiled in LabView, in which the endpoint is UVC detected photometrically at 365 nm. For DNA extraction, ten litres of seawater was ltered through 0.

Dna Extraction And Sequencing
The Sterivex lter cartridge was cracked open aseptically using a clean plyer. The lter paper was cut with a sterile blade and added to a disruptor tube containing SLX Mlus buffer, and DNA extraction was carried out as per the manufacturer's protocol (OMEGA BioTek, USA). The vacuum dried DNA samples were outsourced for sequencing to Euro ns Genomics (India). The QC passed samples were processed for amplicon generation followed by NGS library preparation using Nextera XT Index Kit (Illumina Inc.). The primers used to amplify the bacterial 16S V3-V4 region were 16S rRNA Forw GCCTACGGGNGGCWGCAG and 16S rRNA Rev ACTACHVGGGTATCTAATCC (Euro ns). The raw reads were used for further analysis.

Sequence Analysis
High-quality clean reads obtained from Illumina were selected using Trimmomatic v 0.38 [Bolger et al. 2014], and the paired-end data was stitched into single-end reads using FLASH (v1.2.11) [Magoč and Salzberg 2011]. ASVs were generated using the DADA2 plugin in the QIIME2 pipeline [Caporaso et al. 2012], and the taxonomic identity was assigned using the Silva database (version 13_8) [McDonald et al. 2012]. The resulting biome le was used for further analysis. Abundance estimation along with diversity analysis and rarefaction analysis was carried out using the QIIME program. The sequence data of the 16S rRNA gene was submitted to National Center for Biotechnology Information (PRJNA 706932; SUB9205387).

Statistical Analyses
The variability in the bacterial classes was compared with 5 environmental variables, namely, temperature, DO, nitrate, nitrite and ammonia using Canonical Correspondence Analysis (CCA) using PAST-4 software V4.03 [Hammer et al. 2021]. QIIME2 pipeline was used for diversity analysis. Beta diversity at the phyla level between the depths and stations was calculated using multivariate scaling analysis (MDS) based on Bray-Curtis distance using Primer 6 software [Clarke et al. 2008]. The PCoA analysis was carried out and plots were obtained using Unweighted and weighted distance metrics in QIIME2.

Predictive Functional Analysis
Predictive functional analysis of the microbial communities present in all eight samples was done using Tax4fun2 [Wemheuer et al. 2020]. The 16S rRNA sequences were compared with reference sequences in BLAST using the runRefBlast function. The predictive functions were calculated based on KEGG ortholog reference pro les using the Functional Prediction function. Heatmap and statistical analysis representing the functional predictions within the samples was done using the STAMP V2. Tukey-Kramer test was carried out to check the differences in functional gene composition within the samples.

Physico-chemical characteristics of the water samples
The temperature of the water samples during the study period varied between 21.8 and 32°C, whereas salinity varied between 33.83 and 35.78 PSU. The temperature at the bottom depths was lower than the upper depths during the months from July to Oct 2018, whereas no signi cant differences were observed in the salinity ( Table 1). The DO levels showed a drastic change from July to October 2018, ranging from oxic conditions (2.3  (Table 1). In 2018, higher nitrate concentrations (3.1 -6.9 µM) were observed in September compared to other months. High nitrate was also observed in October 2019. In contrast to nitrate, higher ammonium concentrations (6.13 -7.06 µM) were observed during October 2018, accompanied by low nitrate values (0.01 -0.06 µM) and the absence of nitrite. High ammonium concentration (9.6 µM) was also observed in the near-bottom waters during July 2018. H 2 S (1.66 µM) was detected at the lower depth, indicating anoxic conditions in October 2018 (Table 1). Among the archaeal phyla, Thermoplasmatota showed dominance during July, October 2018 at both depths. However, during October 2019, Thermoplasmatota showed high abundance in the mid-depths. The two euryarchaeal groups present in all the samples belonged to marine group II and marine group III, with the former being dominant in all the samples. In addition, archaeal phyla Crenarchaeota and Nanoarcheota were present in most samples, albeit not so dominant.
A comparative analysis of relative abundance at the class level between mid and near-bottom depths across all the seasons highlighted that the class Acidobacteria belonging to phyla Actinobacteria showed high abundance in all samples. Among the Proteobacteria, the major classes were Alphaproteobacteria followed by Gammaproteobacteria (Fig. 2). It was noted that Alphaproteobacteria were abundant in mid-depth compared to the near-bottom depth in all the samples. Thermoplasmata belonging to the Phylum Thermoplasmatota showed higher abundance in almost all samples except GV3, GV4 and GV8.
At the genus level, Candidatus Actinomarina, NS4 Marine group, NS5 Marine group, Marinimicrobia, clade 1a, AEGEAN-169 marine group, Marine group III, Sva0996 were ubiquitously present in all the samples. In addition, the SUP05 clade was also present in almost all the samples except GV2.
A comparison between the samples from mid-depths GV1, GV3, GV5 and GV7 showed a similar trend. Proteobacteria was the dominant phylum in almost all samples except for GV5, irrespective of sampling time.
The Shannon alpha diversity analysis of surface depths, i.e. GV1, GV3, GV5 and GV7, indicated that GV3 exhibits higher alpha diversity followed by GV1, GV7 and GV5 (Table 2). Similarly, in the lower depths, Shannon alpha diversity analysis revealed that GV6 displayed higher diversity, followed by GV2, GV4 and GV8. MDS and PCoA analysis showed distinct grouping at phyla level at the mid-depths and near bottom depths. The effect of environmental variables on microbial community structure at the genus level was studied at the G5 station along both depths. Seasonal and annual variations in community structure with respect to 5 environmental variables were investigated using Canonical Correspondence Analysis (CCA). The triplot obtained reveals a strong positive correlation between dissolved oxygen, ammonia and the genus NS4_marine group, Clade_II and SAR 116 clade belonging to the class Alphaproteobacteria (Fig. 3). Archaeal marine group II also shows a slight positive correlation with these two parameters. This correlation was essentially observed in samples GV1 and GV2. NS5 marine group, Marinobacterium, clade 1a, PAUC34f and bacteroidetes showed a positive correlation with nitrate and negative correlation to temperature. Candidatus nitrosopelagicus, SAR202 clade, Marinimicrobia, SAR324, AEGEAN-169 marine group were all negatively correlated to nitrite concentration. Among all the environmental parameters analysed, DO, and ammonia signi cantly contributed to the variation in the community in the study area.

Predictive Functional Analysis
The functions predicted by Tax4fun2 highlighted various diverse functions performed by the organisms that are essential for numerous ecological processes. The dominant pathways reported were involved in metabolism, such as carbohydrate, lipid, energy, nucleotide, amino acid metabolism, genetic information processing ( Supplementary Fig. 3), cellular processes, and environmental information processing. Few of the abundant genes predicted in all the samples were sul te reductase (K00392), thiosulfate sulfurtransferase (K01011), Nitrogen regulatory protein Pii1 (K02589), thioredoxin 1 (K03671), DNA repair and replication protein RecF (K03629) etc. (Supplementary Fig. 3).

Discussion
The Western Shelf of India is considered one of the most productive areas of the Arabian Sea due to nutrient enrichment caused during the Southwest monsoonal upwelling [Naqvi et al. 2000[Naqvi et al. , 2006Shetye et al. 1990].
This area experiences considerable changes in the DO concentration throughout the year. It remains welloxygenated from November -May, subsequently reaching hypoxic conditions by July and intensifying further during September-October. Hypoxia begins in June, followed by a decrease in the DO levels and an increase in denitri cation during late August-early September, further progressing towards anoxia associated with sulphidic conditions during September-October [Naqvi et al. 2006[Naqvi et al. , 2010.
The reduction of DO levels at the sampling site is a characteristic feature observed during the SWM [Naqvi et al. 2006] and is fairly evident in the study. Even in the present study, the conditions switched from hypoxic to suboxic to anoxic from July to October 2018. The reduction in the DO values may be attributed to the consumption of DO during remineralisation of the organic matter produced as a result of upwelling [Naqvi et al. 2000]. Concomitant with the reduction in the DO levels, nitrogen is also lost during September -October due to denitri cation. This is evident from the loss of nitrate and nitrite during these months. These results correspond to the studies carried out by Naqvi et al. [Naqvi et al. 2006 , and these are known to encompass e cient denitri ers. The predominant class, Acidobacteria, is known to be involved in nitrate reduction [Kielak et al. 2016].
The archaeal community was dominated by Euryarchaeal class Thermoplasmata, marine group II and marine group III in all the samples. A higher abundance of marine group II Euryarchaea observed in these hypoxic and anoxic waters were parallel to that reported in Arabian Sea OMZ, in the black sea, and anoxic Cariaco Basin wherein this phylum was present throughout the water column [Bandekar et al. 2018;Madrid et al. 2001;Vetriani et al. 2003]. Marine group II archaea have been known to play an important role in carbon cycling [Lincoln et al. 2014;Zhang et al. 2015]. The higher relative abundance of Euryarchaea during October 2018 seems to indicate its a nity to thrive in anoxic waters, as reported previously by Belmar et al. [2011]. It has also been suggested that the concurrent aerobic and anaerobic processes in oxygen-de cient zones may lead to higher archaeal diversity in those regions [Fussel 2013]. The presence of Nanoarchaeato was reported from all the water samples, though at very low abundance.
The various genera present in all the samples have been known to play a signi cant role in biogeochemical cycles. The OM60/ nor 5 clade belonging to the class Gammaproteobacteria are aerobic anoxygenic photoheterotrophs that play an important role in the carbon cycle [Yan 2009 One of the largest naturally occurring oxygen-de cient zone develops during the late southwest monsoon, over the west coast of India, resulting from upwelling, which leads to depletion of DO in that area [Naqvi et al. 2000]. These low oxygen conditions are known to have severe effects on the coastal marine environment resulting in a loss of marine biodiversity. The prokaryotic communities are intensely affected by seasonal changes and respond by shifting community composition based on these environmental factors. Using high throughput sequencing methods, the abundance and distribution of various microbial lineages can be studied. These studies will serve as a useful tool to interpret the importance of various bacterial and archaeal groups in regulating biogeochemical cycling in hypoxic environments.

Conclusion
Coastal ecosystems are greatly in uenced by anthropogenic activities, which can be attributed to growing industrial development and pollution. These activities have the potential to in uence the natural ow of events, such as the development of seasonal hypoxia along the west coast of India. This eventually adds to the declining oxygen concentration in the coastal environments and leads to changes in ecological functioning.
High throughput sequencing techniques can be used to study the change in microbial diversity, which provides insights into the microbial diversity as well as distribution of various microbial communities. Microorganisms contribute signi cantly to the biogeochemical cycling of several essential elements such as nitrogen, carbon, phosphate, and sulfur. Therefore, the in uence of oxygen depletion on microbial community structure affects coastal environmental management and the global environment. Studying the implications of low oxygen concentration and other environmental variables on bacterial community structure provides insights into the processes that might be affected due to the development of hypoxia.

Author contributions
The study was conceptualized by the corresponding author (SD). The sampling and lab work was carried out by the rst author (VN), third author (SS) and fourth author (ABM). SD, VN and SS carried out the amplicon sequence data analysis. The manuscript was written by SD, VN, and the last author (DMS). All authors have read and approved the submission of the manuscript.

Ethics approval
Not applicable Consent to participate Not applicable Consent to publish Not applicable Figure 1 Percentages of relative abundance (Y-axis) and sampling cruise/ eld trip and depth-wise representation of phyla (X-axis) at the Candolim Time-series station (CaTS) G5.

Figure 2
Percentages of relative abundance (Y-axis) and sampling cruise/ eld trip and depth-wise representation of class (X-axis) at the Candolim Time-series station (CaTS) G5.

Figure 3
Canonical Correspondence Analysis ordination Triplot of bacterial communities at phyla level at the Candolim time-series station G5 associated with environmental variables. At the end of the spokes are the environmental parameters, black are the depths sampled in each cruise/ eld trip, and blue dots represent the bacterial phyla.