Surface seawater samples were collected from four stations: St. MPA-2, located in the MPA in Sulaibikhat Bay; St. K6, located in Kuwait Bay; St. A, located in northern coastal waters to the northeast of Kuwait Bay; and St. 18, located in southern offshore waters to the southeast of Kuwait Bay (Figure 1 and Table 1). The sampling was conducted for four consecutive months (November and December 2019, and January and February 2020). Samples were collected using 5 L Niskin bottles (General Oceanics, Florida, USA) and stored in sterilized 5 L polypropylene carboys (Azlon cat# BNP05B) for nucleic acid extraction and shotgun metagenomic sequencing. Seawater samples for the analysis of dissolved inorganic nutrients such as nitrate (NO₃^–), nitrite (NO₂^–), phosphate (PO₄³), and silicic acid (SiO₄^4–) were transferred into 1 L acid-washed polyethylene bottles for laboratory processing. Seawater samples for ammonium (NH₄) measurements were filtered using syringe filters with 25 mm 0.7 μm glass fiber filters (GF/F) and transferred into clean 125 ml polyethylene bottles. The sample volume was then divided into a triplicate set of 10 ml each in 14 ml test tubes with airtight caps. Phenol and hypochlorite reagents were added in the field. A total of 6 samples were collected from St. MPA-2, and a total of 4 samples were collected from each of St. K6, A, and 18. Complete sampling details including GPS coordinates, time, date, etc. are provided in Supplementary Table 1.

Data describing different physicochemical aspects of the environment at the time of sample acquisition was previously described (Fakhraldeen et al., 2023). The environmental parameters measured were total sampling location depth (m), dissolved O₂ levels (ml/L), salinity, and seawater temperature (°C), all of which were measured using a portable water quality profiler equipped with pre-calibrated sensors (AAQ-RINKO, JFE Advantech, Hyogo, Japan). Wave height (m) was measured using pressure-based wave and tide acoustic sensors and sea color was assessed using the Forel and Ule scales. The measured parameters related to air were air temperature in shade (°C), wind direction, wind speed (kph), relative humidity (RH%), and barometric/atmospheric pressure (hectopascals), all of which were assessed/measured using a portable weather station (Kestrel, Pennsylvania, USA). A visual assessment of cloud cover was also made. All measurements are tabulated in Supplementary Table 2.

Seawater samples for nutrient analysis were filtered through 0.7 μm glass fiber filters (GF/F, Whatman, Maidstone, UK), transferred to 125 ml acid-washed polyethylene bottles, and stored at −20°C until analysis. A colorimetric-based method was used to measure the amount of nutrients (excluding ammonium) in the samples using an automated analyzer (Skalar San ++, Breda, The Netherlands) (ROPME, 2010). A combination of the Sagi and Solorzano methods (Sagi, 1966; Solorzano, 1969) was applied to measure ammonium levels. Optical density measurements of the levels of indophenol blue that is produced when ammonium reacts with the Berthelot reagent (containing phenol and hypochlorite) at a wavelength of 630 nm using a double-beam scanning spectrophotometer were performed (Shimadzu UV-2600, Kyoto, Japan). Triplicate measurements were carried out for each nutrient in each sample.

DNA isolation from seawater samples was performed as previously described (Fakhraldeen et al., 2023). Approximately 1 L of each sample was subjected to vacuum filtration through a 0.22 μm polycarbonate membrane (Isopore/Millipore cat# GTTP040700, Darmstadt, Germany). Membranes were then cut using sterile razor blades, and the cut membrane pieces were used as the starting material for DNA isolation. DNA extraction was carried out using the PowerSoil^® DNA Isolation Kit (Qiagen, Hilden, Germany) as per the manufacturer’s instructions. DNA yield assessment was performed using a NanoDrop (Thermo Fisher Scientific, Massachusetts, USA). The quality of the extracted DNA was also assessed through gel electrophoresis. The eluted DNA was run on a 1% agarose gel containing 10 mg/ml ethidium bromide at 100 V for 60 min. Gels were then imaged using the ChemiDoc MP gel documentation system (Bio-Rad, California, USA) to ensure the absence of sheared DNA. The average 260/280 ratio measured by Nanodrop was 1.965, and all bands detected on agarose gels were sharp and intact with little to no smearing observed.

Library preparation, DNA sequencing, and data processing were performed as described in Fakhraldeen et al. (2023). A total of 1 ng of DNA was processed using the Nextera^§ XT DNA Library Preparation Kit (Illumina, California, USA). Library quantity and quality were assessed using the Qubit 2.0 system (Thermo Fisher Scientific, Massachusetts, USA) and the Tape Station 2200 system coupled with high sensitivity D1000 Screen Tapes (Agilent Technologies, Inc., California, USA), respectively. Libraries were pooled in equimolar amounts based on quality control values and then sequenced using the HiSeq platform (Illumina, California, USA) to generate 2 × 150 bp reads at a depth of 20 million paired-end reads per sample (10 million in each direction). Trimmomatic (version v0.38) (Bolger et al., 2014) was then used to cut adapter sequences and trim low-quality bases with the following settings: phred score 33, SLIDINGWINDOW of 4:20, LEADING and TRAILING threshold of 20, and drop if length <100 (MNLEN 100). FastQC (version v0.11.8)¹ was applied to check the quality of raw reads. MEGAHIT (v1.2.9) NGS assembler (Li et al., 2015; Li et al., 2016) was used to assemble the trimmed fastq data and generate assembled contigs of average length 642 bp. Average contig lengths generated by MEGAHIT per sample are provided in Supplementary Table 2.

MetaPhlAn (v3.0) was used to profile the composition of microbial communities (Segata et al., 2012). MetaPhlAn’s latest database files were downloaded from http://cmprod1.cibio.unitn.it/biobakery3/metaphlan_databases/. The most recent mpa file for MetaPhlAn3 is mpa_v30_CHOCOPhlAn_201901. This file was used with the “– index” parameter. The trimmed reads generated by MEGAHIT were used as inputs for MetaPhlAn3 to increase the sensitivity of the output. Default values were kept for all other parameters. “Very-sensitive” mode was used when aligning reads. The minimal mapping quality was set to 5. The minimum total nucleotide length for the markers in a clade for estimating the abundance without considering sub-clade abundances was set to 2000. Additionally, the threshold for the percentage of markers with a non-zero relative abundance for misidentified species was established at 33%.

The R package metacoder (v0.3.4) was used to generate the heat trees (Foster et al., 2017). The node size and color are mapped relative to the observed read numbers and Operational Taxonomic Unit (OTU) counts, respectively.

The R package “microbiome 1.13.10” (Lahti and Shetty, 2018) was used to analyze the core microbiome at the family, genus, and species levels. The “plot_core” function was used to generate the core microbiome figure. The minimum prevalence was set to 0.1 and applied to the rows (taxa) and the columns (detections).

Alpha diversity, represented by four indices (Chao1, Shannon, Simpson, and InvSimpson), was estimated using the R package phyloseq (v1.34.0) (McMurdie and Holmes, 2013). The alpha diversity was compared between each pair of stations (St. MPA-2 (MPA in Sulaibikhat Bay) versus St. K6 (Kuwait Bay), St. MPA-2 versus St. A, and St. MPA-2 versus St. 18) through One Way Analysis of Variance (ANOVA) (p < 0.05).

Principal component analysis (PCA) was performed using the R function prcomp in R package stats4 (v4.0.5)² using singular value decomposition of the centered and scaled data matrix (taxonomy count table). The principal components (PC1, 2, and 3) from the samples are shown on 3D scatter plots generated using scatterplot3d (v0.3-41). The proportion of variance is also depicted on the plots.

Differential heat trees were generated using the same tool and approach described above, except for the node color being mapped to the log 2 median ratio of abundance values (OTU counts) between various pairs of samples instead of the raw OTU counts.

Analysis of Similarity (ANOSIM) was performed using the R package “anosim” in the package vegan (v.2.5-7) (Oksanen, 2010). Differential abundance analysis was performed using the Metastats tool (White et al., 2009). The tool outputs the mean value, variance, standard deviation, p-value, and q-value. The top differentially abundant taxa are displayed graphically.

The biotic and environmental datasets were initially screened using the D’Agostino-Pearson Omnibus Normality Test. As the test suggested a non-normal distribution of data points, the datasets were processed using non-parametric statistical methods. In the environmental metadata, water temperature, salinity, and dissolved oxygen datasets showed no distribution skewness. Hence, their values were not transformed. Values of all nutrient variables and turbidity levels were individually transformed by taking the square root of each value. The selectively transformed and untransformed variables were standardized before constructing the Euclidean distance-based similarity matrix required for ordination analyses. The bacterial OTU abundance values were transformed by taking the natural log of each value in the biological dataset before constructing the Bray-Curtis dissimilarity matrix. The dominance and diversity of bacterial taxa (at different taxonomic levels) at different sampling areas were investigated through cumulative k-dominance plots. The constraint ordination procedure, Canonical Analysis of Principal Coordinates (CAP), was used to explore the microbial assemblage patterns at the species, genus, and family levels. Before the CAP analysis, the sample clustering patterns were probed through the Similarity Profile (SIMPROF) analysis and its Pi statistics. The Permutational Multivariate Analysis of Variance (PERMANOVA) explored the differences in microbiome structure between the sampling locations and the significance of clustering patterns through its pseudo-t-statistics and associated permutational and Monte Carlo-based significance tests. The linear combinations of environmental variables bearing correlation of >0.2 (Pearson) with the biotic assemblages and CAP axes were overlaid as vectors on the CAP plot. Additionally, 3D bubble plots were constructed to show the distribution of environmental variables across the sample groups. A colored, interactive heatmap with XY clustering for samples and species explored the sample clustering based on bacterial species association patterns. A Venn diagram was used to explore the species turnover (i.e., beta diversity) between the sample groups and to enlist the endemic and common species. The relationship between biotic and environmental rank matrices was explored using RELATE, a Mantel test based on Spearman rank correlation. The multivariate Biotic-Environmental interactions procedure (BIOENV) was used to find the best combination of environmental variables responsible for the bacterial assemblage patterns.

All multivariate analyses were performed in the ecological statistical software PRIMER v7.0.11 (PRIMER E-Ltd., UK) (Clarke and Gorley, 2015) with PERMANOVA + 1 add-on for PRIMER (Anderson, 2008). The Venn diagram was constructed using the freely available Bioinformatics & Evolutionary Genomics software.³

Libraries were used to generate paired-end sequences (150 × 2) resulting in an average of 28.3 million reads total and an average of 14.2 million paired-end reads per sample (Phred > Q20 for 97% of reads). Reads were subjected to quality filtering and trimming resulting in the removal of 1.3% of reads. The retained reads yielded approximately 27.97 million clean, processed reads per sample for downstream analysis. Reads from each sample were assembled de novo into scafftigs with the largest scafftigs ranging in length from 37,426 to 276,400 bp. The N50 values of the assembled genomes ranged from 571 to 1024. Detailed read statistics are available in Supplementary Table 3.

A number of physicochemical parameters were assessed at the time of sample acquisition from the four different sampling locations throughout the study period (Supplementary Table 4). The depth of the stations was logged, and measurements showed that the sampling locations in order of increasing average depth were St. MPA-2 (1.8 m) < St. A (6.5 m) < St. K6 (10.7 m) < St. 18 (20.6 m). Average surface seawater sample acquisition depth was 1.0 m for St. MPA-2 and 1.3 m for the other three stations. Dissolved oxygen (DO₂) levels were also measured and showed no significant difference between the sampling stations throughout the study period. Average DO₂ levels, in increasing order, were recorded at St. 18 (6.8 mg/l) < St. K6 (7.6 mg/l) < St. A and St. MPA-2 (both 7.8 mg/l).

Significant differences were recorded, however, in seawater temperature and salinity levels between the stations. The highest average seawater temperature was recorded at St. 18 (22.1°C), which was significantly greater than the average temperatures recorded at both Stations A (17.8°C; p = 0.012) and K6 (19.2°C; p = 0.013). Average temperatures recorded at St. MPA-2 (19.6°C) was not significantly different from the other stations (p = 0.241, 0.795, and 0.242 for Stations A, K6, and 18, respectively; Student’s t-test). Stations K6, A, and 18 all had similar salinity levels ranging from 39.1 at St. A to 40.8 at St. K6. Salinity levels at St. MPA-2 were significantly higher, though, at 42.4 (p < 0.03).

Overall, St. MPA-2 in Sulaibikhat Bay was the shallowest of the four stations examined and had the highest average salinity levels. The offshore station located to the southeast of Kuwait Bay, St. 18, had the highest recorded average seawater temperature compared to the other three stations. Dissolved oxygen levels were similar between all four stations throughout the study period.

Taxonomic profiling was performed on the six samples collected from St. MPA-2 in Sulaibikhat Bay to generate a profile of the bacterial and archaeal communities within that water body. Depictions of the twelve most abundant taxa at the family, genus, and species levels are provided (Figure 2).

Results documented 18 families, including 17 bacterial families and one archaeal family, with relative abundances (RA%) of 92.78 and 7.22%, respectively. Of the bacterial families, 50% (9/18) belonged to the phylum Proteobacteria, 22% (4/18) belonged to the phylum Bacteroidota, 11% (2/18) belonged to the phylum Cyanobacteria, and then one family each belonged to the phyla Actinomycetota and Campylobacterota. The Proteobacteria phylum was represented by two classes, Gammaproteobacteria (5/9) and Alphaproteobacteria (4/9). The phylum Bacteroidota was represented equally by the two classes, Flavobacteriia and Bacteroidia. Furthermore, the family Rhodobacteraceae (class Alphaproteobacteria) had the highest RA (60.64%), followed by the families Oceanospirillaceae (class Gammaproteobacteria; 11.66%) > Flavobacteriaceae (class Flavobacteriia; 9.67%) > Erythrobacteraceae (class Alphaproteobacteria; 3.24%) > Microbacteriaceae (class Actinobacteria; 2.84%) > Halieaceae (class Gammaproteobacteria; 2.74%). The families Campylobacteraceae, Crocinitomicaceae, Synechococcaceae, and Pelagibacteraceae all had RAs between 0.2 and 1% while the remaining families all had RAs < 0.1% (Figures 2A, D and Supplementary Table 5).

A total of 26 genera were detected. The two families with the most genera were Rhodobacteraceae [five genera in decreasing order of RA Rhodobacteraceae_unclassified (21.09%) > Phaeobacter (15.4%) > Thalassobacter (13.62%) > Lentibacter (9.16%) > Nereida (1.37%)] and Flavobacteriaceae (five genera in decreasing order of RA: Dokdonia (6.21%) > Gilvibacter (1.68%) > Polaribacter (1.19%) > Formosa (0.58%) > Cloacibacterium (0.01%). All other families were represented by just one genus each. The five most abundant genera from this latter category were Bermanella (11.66%), Candidatus_Heimdallarchaeota_unclassified (7.22%), Altererythrobacter (3.24%), Candidatus_Aquiluna (2.84%), and Halieaceae_unclassified (2.74%) from the families Oceanospirillaceae, Candidatus_Heimdallarchaeota_unclassified, Erythrobacteraceae, Microbacteriaceae, and Halieaceae, respectively (Figures 2B, E and Supplementary Table 5).

A total of 29 species were documented. The sole archaeal species was Candidatus_Heimdallarchaeota_archaeon belonging to the phylum Candidatus_Heimdallarchaeota and an unclassified family. The ten most abundant bacterial species in decreasing order of abundance were Rhodobacteraceae_bacterium_HIMB11 (21.09%) > Phaeobacter_italicus (15.4%) > Thalassobacter_ stenotrophicus (13.62%) > Bermanella_sp (11.66%) > Lentibacter_algarum (9.16%) > Dokdonia_donghaensis (6.21%) > Altererythrobacter_sp_XM_24bin4 (3.24%) > Candidatus_ Aquiluna_sp_XM_24bin5 (2.84%) > Gammaproteo bacterium_IMCC3088 (2.74%) > Gilvibacter_sp_SZ_19 (1.68%) (Figures 2C, F and Supplementary Table 5).

The core microbiome within the six samples collected from St. MPA-2 in Sulaibikhat Bay was assessed. Results showed that the 12 dominant families included Rhodobacteraceae (represented by five genera collectively accounting for an overall prevalence of 27.8%) and Flavobacteriaceae (represented by four genera collectively accounting for an overall prevalence of 15.3%). The following families were represented by one genus each unless otherwise stated: Erythrobacteraceae (6.94%), Campylobacteraceae (6.94%), Microbacteriaceae (5.56%), Oceanospirillaceae (5.56%), Halieaceae (4.17%), Pseudoalteromonadaceae (represented by two genera collectively accounting for an overall prevalence of 2.8%), an unclassified archaeal family (2.78%), Synechococcaceae (2.78), Crocinitomicaceae (2.78%), Pelagibacteraceae (1.39%), and Bacteroidaceae (1.39%) (Supplementary Figure 1A and Supplementary Table 6).

The 12 most prevalent genera in descending order of abundance were Phaeobacter (8.33%), Rhodobacteraceae_unclassified (8.33%), Altererythrobacter (6.94%), Arcobacter (6.94%), Bermanella (5.56%), Candidatus_Aquiluna (5.56%), Thalassobacter (4.17%), Lentibacter (4.17%), Dokdonia (4.17%), Halieaceae_unclassified (4.17%), Gilvibacter (4.17%), and an unclassified archaeal genera Heimdallarchaeota_unclassified (2.78%) (Supplementary Figure 1B).

31 species were detected at a prevalence above 0.1. The top 12 species represented in the core microbiome were Phaeobacter_italicus (8.33%), Rhodobacteraceae _bacterium_HIMB11 (8.33%), Altererythrobacter_sp_XM_24bin4 (6.94%), Bermanella_sp (5.56%), Thalassobacter_stenotrophicus (4.17%), Lentibacter_algarum (4.17%), Candidatus_Aquiluna _sp_XM_24bin5 (5.56%), Dokdonia_donghaensis (4.17%), Arco bacter_cryaerophilus (6.94%), Gamma_proteobacterium_IMCC3088 (4.17%), Gilvibacter_sp_SZ_19 (4.17%), and Candidatus_Heimdallarchaeota_archaeon (2.78%) (Supplementary Figure 1C).

Biodiversity monitoring of Sulaibikhat Bay has been attempted previously, with efforts made to investigate the endemic macrofauna (Al-Zaidan et al., 2003), benthic foraminifera (Al-Enezi and Frontalini, 2015), and higher organisms such as mud skippers (Al-Behbehani and Ebrahim, 2010) within this environment. Microbial communities associated with the micro and macrofauna of Sulaibikhat Bay have been limitedly explored through conventional culture-dependent techniques and microscopic examination (Al-Mohanna et al., 2007; Al-Zaidan et al., 2006). Genome-centric approaches are gaining popularity over such traditional methods for assessment of biodiversity (Doytchinov and Dimov, 2022). Presently high throughput shotgun metagenomic sequencing of DNA isolated from the surface waters of the MPA within Sulaibikhat Bay as well as three other locations distributed throughout neighboring waters in and around Kuwait Bay successfully mapped the bacterial and archaeal communities at these locations. Not surprisingly, bacterial communities were more abundant and diverse compared to archaeal communities. Shotgun metagenomic techniques have been used by numerous research groups for taxonomic and functional annotations of microscopic organisms including bacteria, archaea, viruses, and fungi in aquatic environments (Doytchinov and Dimov, 2022).

Relatively high abundances of the bacterial families Rhodobacteraceae and Flavobacteriaceae were recorded in the surface seawater samples collected from St. MPA-2 in Sulaibikhat Bay. Both these bacterial families were also found in surface waters of Southeast Florida with reported functional roles in hydrolysis of complex molecules, photosynthesis, and degradation of algal products (Campbell et al., 2015). Rhodobacteraceae was also the core bacterial family in surface seawater collected from the South China Sea and the coastal waters of Hong Kong and Xiamen (Xia et al., 2021). Flavobacteriaceae were exclusively abundant in the warm oligotrophic waters of the South China Sea (Xia et al., 2021). Both these families were associated with Synechococcus cultures in the study conducted by Xia et al. (2021). In the present study, the dominant genus within the family Rhodobacteraceae that was detected within the samples collected from Kuwaiti waters was Phaeobacter, and the dominant species was P. italicus. Similarly, the bacterial genera Cognatishimia and Phaeobacter from the family Rhodobacteraceae were found in the surface waters of Bay of Bengal (Akter et al., 2023).

The Sulaibikhat embayment is considered a pristine habitat, biologically diverse and linked with commercial fishing activities in Kuwait. In the recent past, some signatures of physicochemical stresses such as changing water and sediment profiles, and biological activities were noticed. For example, effluents from nearby desalination plants and the Ghazali outfall are affecting the water column and increasing organic content within the water body. Pyrite deposition in the vicinity of the Ghazali outfall is indicative of sewage accumulation (Alshemmari et al., 2011). In view of such anthropogenically induced eutrophication and the subsequent putative harmful impacts they impart, it is essential to continually biomonitor affected environments, especially protected areas. The relative abundances of the bacterial and archaeal taxa differed between St. MPA-2 and the other sampling locations distributed in and around Kuwait Bay. The variations in the relative abundance of various microbial taxa are likely due to variable environmental and physiochemical conditions during the sampling period, which spanned a period of 4 months.

The bacterial richness and distribution recorded within surface waters of St. MPA-2 was significantly greater than those recorded in St. K6, St. A, and St. 18. This is likely attributed to variations in nutrient availability between the locations. Higher levels of inorganic phosphorus and nitrogen were recorded in St. MPA-2 compared to the other locations. The source of these abiotic elements is most likely the effluents received through thirteen discharge points around the Bay (ROPME, 2010). The inorganic and organic load introduced through the effluent ports provide a breeding bed for several microbial communities and likely contribute heavily to the observed richness in distribution and diversity of these communities within these waters. Similarly, greater Shannon and ACE diversity indices were reported for the areas along the tropical seawater habitats near outfalls adjacent to Broward County in the state of Florida (Campbell et al., 2015). Higher Shannon and Faith indices of 18S rRNA operational taxonomic units (eukaryotic microalgae) were also observed at similar locations along Kuwait’s coast (Kumar et al., 2021). Hahnke et al. (2013) reported that bacterial communities are generally more diverse in the phototrophic layers. Sulaibikhat Bay is a shallow, tidal mudflat with a water column spanning ten meters (maximum depth at the Bay entrance). The designated MPA within Sulaibikhat Bay, however, has an average depth of only 2–3 m. The cumulative k-dominance analysis was also suggestive of the dominance of bacterial taxa and increased diversity in St. MPA-2 compared to the other locations distributed throughout the surrounding waters (St. K6 in Kuwait Bay, St. A in northern coastal waters, and St. 18 in southern offshore waters).

Additionally, differential abundance analysis revealed the taxa that are overrepresented in St. MPA-2 as compared to the other three stations. The species Phaeobacter_italicus, family Rhodobacteraceae, and order Rhodobacterales were more abundant in St. MPA-2 (p < 0.05) versus in St. K6, St. 18, and St. A, respectively. These microbes play a role in carbon assimilation and hence their enrichment in eutrophic marine environments is not surprising (Campbell et al., 2015). Rhodobacteraceae is a group of bacteria that exist in marine environments and are functionally involved in biogeochemical cycling of elements such as sulfur and are known to establish symbiotic relationships with other aquatic microorganisms (Pujalte et al., 2014). The species assemblage analysis performed in this study corroborates the differential abundance findings. Namely, 24 species were found exclusively in Group I (St. MPA-2) and 11 species were specific to Group II (St. A, St. 18, and St. K6; Figure 9).

The microbial populations detected in St. MPA-2 present as a distinct group (Group I) on the canonical analysis plots. The other three stations in Kuwait Bay and surrounding waters collectively formed a separate cluster (Group II). Similar clustering was observed through PCA. The spatial distance and environmental gradients between the MPA and other station locations seems to be the driving factor for the observed clustering pattern. Station K6, located within Kuwait Bay, is the most proximal, whereas St. A is located northeast of the Bay between the Kuwaiti coastline and the Kuwaiti island of Bubiyan, and Station 18 is located in offshore waters southeast of the Bay and is the most remote and distant (Fakhraldeen et al., 2023). These findings were further substantiated by the interactive heatmap clustering. In addition to the spatial differences, varying environmental conditions in the MPA in Sulaibikhat Bay are also speculated to play a role in shaping the microbiome, as environmental parameters such as temperature, salinity, and nutrient (NOx, phosphate, and silicic acid) levels affect the bacterial community structure, as reported from the southern coastal waters of Korea (Cui et al., 2019). In the present study, a high correlation between the bacterial taxa and phosphate levels (Rho > 0.700) suggests this inorganic nutrient significantly affects the microbiome in the Sulaibikhat Bay MPA. A study by Xia et al. (2021) reported that the members of the families Rhodobacteraceae and Flavobacteriaceae, which are abundant in the MPA, are rich in ABC transporters involved in importing nutrients including phosphate, amino acids, and other carbon-containing compounds (Xia et al., 2021). Additionally, members of the family Microbacteriaceae, which is significantly enriched in the MPA, have been reported to play significant roles in phosphate solubilization (Rivas et al., 2004). Interestingly, a negative correlation between excess phosphate loading and abundance of members of the family Flavobacteriaceae was reported in oligotrophic coastal waters indicating that perhaps the effect of phosphate levels on aquatic microbiomes, at least as they related to the family Flavobacteriaceae, may vary in a context specific manner (Chen et al., 2016). Increased phosphorus levels were also shown to shape the microeukaryotic communities of the Pearl River estuary, a bell-shaped estuary at the intersection of the Pearl River and the South China Sea (Zhang et al., 2023).

In this study, the sub-clusters assembled within the main clusters in the heatmap analysis demonstrate the intricate nature of the microbiome at these locations. While the combined effect of various environmental variables can influence the overall biotic structure, the taxa-specific environmental preferences or niches along with seasonal differences in environmental settings at each sampling location, can lead to finer-scale differences in the microbial assemblages. Besides environmental gradients, biotic factors such as competition for limiting nutrient resources by the resident phytoplankton community, and predation (viral lysis and microzooplankton grazing) can also alter the microbiome in the northern Arabian Gulf, which is known for its environmental heterogeneity. This complexity underscores the need for further investigation and understanding of the dynamics of the microbial communities residing in this water body.