The study area covers the Saudi territorial waters of the Gulf from Khafji to Ras Tanura (latitude: 26°40′–28°30′ N and longitude: 48°10′–51° 00′ E) (Figure 1). The coastline of this region extends in a more or less straight line in the south-easterly direction, and it is highly convoluted with extensive and complex sheltered bay systems (Basson et al., 1977). The most remarkable feature of the western Arabian Gulf is its extremely low relief. The average slope of the seabed is only 35 cm per kilometer. This results in an unusually extensive intertidal zone, width of which in many places can be measured in kilometers. The depth of the Saudi waters of the Gulf is generally below 60 m.

Sediment samples from the area of the open waters were collected onboard the Cruise 1 of R/V Al-Gosaibi, Saudi Arabia during November–December, 2002. A total of 174 samples from 58 stations (three replicates from each station) along five transects in the Saudi waters of the Gulf (Figure 1) were collected using a Van Veen grab sampler (0.1 m²). Sampling depth ranged from 2.75 to 58 m. Two transects were established along the north-south direction in the coastal (Transects A, 21 stations) and offshore (Transects B, 13 stations) zones, and three transects were established perpendicular to the coast corresponding to the northern (Transects C, 8 stations), central (Transects D, 8 stations), and southern (Transects E, 8 stations) sections of the study area.

Macrobenthic samples were also collected from five locations in two bay systems, viz. Manifa-Tanajib System (MTBS) and Abu Ali Bay System (AABS) during January 2003 (Figure 1). Three stations were located in the MTBS (Tanajib—TN, Balbol—BB, and Manifa—MN) and two in the AABS (Musallamiyah—MS and Mardomah—MD). As oil and gas exploitation in the Gulf is intensive, it is difficult to find control sites (i.e., pristine or undisturbed locations). Hence, the sites were selected based on an initial intensive survey to determine the TPH concentration. Hundreds of core samples were analyzed for the determination of the concentration of TPH in sediments from these locations. Macrobenthic samples and sediment samples were collected from the stations with the lowest (stations including C) and the highest (stations including A) TPH concentrations in these five locations. At each location, three replicate macrobenthic samples were collected from each of two sites.

Altogether 204 macrobenthic samples were used for this study: 174 samples from the open waters and 30 from the bay systems. All the sediment samples collected for the macrobenthic study were sieved through a 0.5 mm mesh sieve and fixed in 5% formalin. Separate grab samples were collected for the sediment characteristics study. As the studies undertaken in the open waters and inner bays were conducted as independent studies, there are differences in the approach of selecting of environmental variables. In the open waters, hydrographical parameters (temperature, salinity, and dissolved oxygen; DO) of bottom waters were measured using a Conductivity Temperature Depth (CTD) profiler and samples were collected for the determination of the grain size and total organic carbon (TOC) in sediments. In the inner bays, the main focus was to test the health of the benthic community, vis-à-vis the concentration of TPH, and hence hydrographical parameters and TOC in sediments were not determined. Sediment samples were collected for the determination of TPH and grain size.

The grain size of the sediments was determined by a combined sieving and pipette method (Carvar, 1971). Sediment textural types were identified using Sheppard's classification (Sheppard, 1954). The concentration of TOC in the sediments was determined by combustion at 900°C and detection with a non-dispersive infrared sensor (NDIR), using a Shimadzu TOC analyzer with an ASI-V autosampler. Standard carbon solutions of 0.5–20 ppm (mg L⁻¹) prepared with potassium hydrogen phthalate C₆H₄ (COOK) (COOH) in C-free water were used for the calibration of the analyzer.

TPH in sediment samples was determined according to the Battelle standard operating procedure, SOP 5–203 (NOAA Status and Trends). Briefly, ~2–30 g of the sediments was serially extracted with dichloromethane (DCM) and sodium sulfate (Na₂SO₄) using an Accelerated Solvent Extraction (ASE) procedure. The sample was weighed into an extraction vessel and spiked with a surrogate internal standard (SIS), solvent was added, and the ASE unit was operated according to the Standard Operating Procedure (SOP). The filtered extract was decanted into an Erlenmeyer flask. The extract was then passed through a 2% deactivated alumina column, and concentrated to 2,000 μl using Kuderna-Danish (KD) and nitrogen evaporation techniques. The TPH concentration in this extract was determined using a gas chromatograph equipped with a flame ionization detector (GC/FID).

According to some guidelines for pollution levels in the sediments of the Gulf (Massoud et al., 1996), they are classified as TPH <15 mg kg⁻¹ is unpolluted, 15–50 mg kg⁻¹ is slightly polluted, 50–200 mg kg⁻¹ is moderately polluted, and >200 mg kg⁻¹ is heavily polluted.

Initially, the specimens were subjected to group level sorting and then species level identification of the dominant taxa such as polychaetes and major crustaceans such as amphipods, ostracods, mysids, and cumaceans. The taxa were identified down to species level using various identification keys (Day, 1967; Fauchald, 1977; Barnard and Drummond, 1982; Ruffo, 1982, 1989, 1993; Hall, 1985; Jones, 1986; McCain, 1993; LeCroy, 2000). The taxa identified to species level represent about 75% of the total macrobenthos density. The remaining taxa were identified at the genus or family levels.

Macrobenthic species abundance data were analyzed to determine the total density of macrobenthic organisms (N), species richness (S), and Shannon-Wiener (H′) (Log₂) for species diversity using PRIMER 6 (Plymouth Routines in Multivariate Ecological Research) (version 6.1.5) (Clarke and Warwick, 2001).

The following biotic indices, which are extensively used by researchers around the world to assess the ecological status of the macrobenthic communities were calculated: (i) the AZTI's Marine Biotic Index (AMBI; Borja et al., 2000), and (ii) the multivariate-AMBI (M-AMBI; Muxika et al., 2007). M-AMBI is a factorial analysis which includes Shannon's diversity, richness, and AMBI. These two indices have been validated as effective ecological status assessment tools for benthic communities in different geographical areas, ranging from tropics to high latitudes, subjected to a range of human pressures, including hydrocarbon impacts (Borja et al., 2003; Muxika et al., 2005; Sigovini et al., 2013; Spagnolo et al., 2014).

AMBI was calculated using AMBI 4.1 software (available in the public domain at http://ambi.azti.es) and the February 2010 species list. Before calculation, the dataset was prepared according to the guidelines in using AMBI (Borja and Muxika, 2005) and M-AMBI (Borja et al., 2008). As some identified species were not included in the species list, the procedure described by Borja and Tunberg (2011) was followed when assigning new species. Additionally, M-AMBI was calculated also using AMBI software. As M-AMBI requires reference conditions to be calculated, and the Gulf is heavily impacted by human activities (see comments and references in the Introduction section), the current structural parameters are not suitable for establishing M-AMBI reference conditions (Borja et al., 2012). Hence, the approach proposed by Borja and Tunberg (2011) and Forchino et al. (2011) was used to establish the reference conditions in the absence of pristine areas. This consisted of increasing the highest diversity and richness values of all replicates by15%, and decreasing AMBI to half of the lowest value. As for the “bad” status, the references were based on the azoic situation (diversity and richness equal to 0 and AMBI equal to 6).

The threshold values for the M-AMBI classification according to Borja et al. (2007) are as follows: “high” ecological status, >0.77; “good,” 0.53–0.77; “moderate,” 0.38–0.53; “poor,” 0.20–0.38; and “bad,” <0.20.

Preliminary analysis indicated that environmental data (depth, bottom water temperature, bottom water salinity, bottom water DO, sand, mud, TOC%, and TPH in the open waters, and sand, mud%, and TPH in the inner bays) and macrobenthic data (S, H′, density, AMBI, M-AMBI, %sensitive species, and %opportunistic species in open waters and inner bays) do not meet the homogeneity of variance assumption. Hence, non-parametric tests were used to analyze the data.

Data from the open waters and inner bays were separately treated due to the different environmental conditions in these two regions, as mentioned in Section Study area. In the case of open waters only two factors were considered—distance to the shore (with two states: nearshore and offshore) and latitude region (with three states: northern, central, and southern); and in the case of inner bays, only one factor was considered—level of TPH (with two states: low and high). In the inner bays, stations MN(C) and BB(C) are lower TPH stations (LTS with TPH concentration <70 mg kg⁻¹) and stations TN(C), TN(A), BB(A), MN(A), MS(C), MS(A), MD(C), and MD(A) are higher TPH stations (HTS with TPH concentration ≥100 mg kg⁻¹). Multivariate analyses were performed on the square-root of the transformed abundance data to detect different macrobenthic assemblages. Bray-Curtis index and group average linkage were used for non-metric multidimensional scaling (MDS) ordination using PRIMER 6.

Significant spatial differences in the environmental and community parameters in the open waters and in the inner bays were established by Permutation Multivariate Analyses of Variance (PERMANOVA) with 999 permutations (Anderson, 2001) in PRIMER 7 (Trial Version 7.0.13). Prior to analysis, the environmental variables were fourth-root transformed and a resemblance matrix was calculated using the Euclidean similarity. Benthic community assemblages and community parameters were square-root transformed and a resemblance matrix was calculated using the Bray Curtis similarity. Significant differences in the nearshore-offshore and northern-southern groups were estimated using the data from the open waters. As weaker dissimilarities were obtained when central was added, it was removed from the tests. In order to test the north-south effect in the nearshore-offshore category, nearshore-offshore was set as the first factor and north-south was set as the second factor. Similarly, in order to test the nearshore-offshore effect in the northern-southern category, northern-southern was set as the first factor and nearshore-offshore was set as the second factor. The interactive effects of both these factors were also tested in the open water data. In the inner bays, significant dissimilarities were established between the LTS and HTS.

Correlations between environmental variables and structural parameters of macrobenthos were calculated using the non-parametric Kendall's tau-b test in SPSS 22 (Release 22.0.0.0). Since multiple analyses were conducted on same variables, to protect from Type I error, Bonferroni correction was applied to keep overall significant level to 0.05.

A total of 392 taxa were identified in the open waters of the Gulf. After assignment, many new species were incorporated to the new AMBI species list (November 2014), which is available free of charge at the AMBI web site. After this assignment, the percentage of unassigned individuals ranged from 0 to 19.8%, with a mean value of 5.2%. Thus, all samples can be taken into account in the analysis.

In the region of open waters, according to the AMBI results, all stations are slightly disturbed, and according to the M-AMBI results, all stations show a high or good status (Figure 4, Table 2). The sensitive species group (EG I) varies from 7.1% (station E5080, southern region) to 48.6% (station C3640, northern region), and opportunistic species (EGs IV and V) from 0.8% (station D4350, central region) to 47% (station D4890, central region). Despite the AMBI homogeneity distribution, there are statistically significant nearshore-offshore differences (higher AMBI values in offshore, PERMANOVA results in Appendix F). A similar difference is also shown by opportunistic species (PERMANOVA results in Appendix F). However, north-south variations within the nearshore-offshore regions are absent. A statistically significant north-south difference is not shown by any of the parameters. However, AMBI shows nearshore-offshore variations within the northern-southern regions (PERMANOVA results in Appendix G). Interactive effects of nearshore-offshore and north-south factors are not noticed in any of these parameters (PERMANOVA results in Appendices F,G).

In the inner bays, AMBI values indicate slightly disturbed conditions at all stations except one [LTS, BB(A)], which is moderately disturbed (Figure 4, Table 2). M-AMBI values show good status for LTS, while, high, good, moderate, and poor status for HTS (Table 2). The sensitive species group varies from 0% [HTS, MD(A)] to 33.9% [HTS, TN(A)], and opportunistic species from 3.8% [HTS, TN(A)] to 59.7% [HTS, BB(A)]. The average values of AMBI, M-AMBI, sensitive species, and opportunistic species show worst status for the HTS compared to the LTS (Figure 4). However, none shows statistically significant differences between the two sites (PERMANOVA results in Appendix H).

The results of the Mann-Kendall correlation analyses between environmental variables and faunal parameters are given in Appendix I. In the open waters, after Bonferroni correction, AMBI is positively correlated to depth (p < 0.05) and mud% (p < 0.05) and negatively correlated to sand% (p < 0.05). Ampeliscid density is positively influenced by depth (p < 0.05), and negatively influenced by salinity (p < 0.05) and DO (p < 0.05).

In the inner bays, significant correlations were not recorded after Bonferroni correction (Appendix I).

Parameters of the bottom waters exhibit nearshore-offshore and north-south variations. A general increase in the temperature and decrease in the salinity and DO toward the south are attributed to the inflow of Indian Ocean water through the Strait of Hormuz (Reynolds, 1993; Abdelrahman and Ahmad, 1995). The high evaporation and turbulence in the coastal areas cause higher salinity and DO in the nearshore waters compared to those of the offshore waters. During late fall, winter, and spring, when the coastal waters are comparatively cooler, the temperature increases toward the offshore areas.

In the Gulf, sediments and their characteristics generally follow distributions that align with the northwest-southeast axis of the Gulf (Wagner and van der Togt, 1973). The predominance of sand in the coastal regions, which are closer to the shore, and shallow in nature, can be attributed to the water circulation prevalent in this area. The gradual replacement of sand by mud fractions from the nearshore to offshore regions is likely to be due to the increasing protection from wave action in the center of the basin, resulting in transport of only fine materials to the center of the Gulf basin (Wagner and van der Togt, 1973). In an ideal situation, higher TOC values are associated with a lower sediment particle size (Khalaf et al., 1986). With the increase of the mud fraction with the depth, TOC also increases, as found in the present study.

Based on the published reports on the TPH concentrations of sediments, Massoud et al. (1996) suggested a TPH value of 10–15 mg/kg (in dry, silt/clay sediment fractions), as the unpolluted or natural background level in the Gulf. In the present study, 17% of the stations, which are mainly in the central and southern regions, has the natural background level of TPH. About 64% of the stations has TPH concentrations >50 mg kg⁻¹.

In general, the benthic density (3,181 ± 359 ind. m⁻²) and H′ (4.9 ± 0.1) in the region of the open waters are comparable to the values recorded in earlier studies in the Saudi waters of the Gulf. For example, a density of 4,710 ind. m⁻² and 2,022 ind. m⁻² have been reported by Coles and McCain (1990) and McCain (1984b), respectively. Also, values in the range of 0.5–5.5 have been reported for H′ by McCain (1984b). A comparatively higher density was observed in the nearshore areas compared to that in the offshore areas, which can be due to the rich primary production in nearshore waters. The supply of food to the subtidal benthic environment depends on the proximity to the shore and the water depth (Levinton, 1982). In general, diversity indices recorded are high at many stations. Diversity values >3 were recorded at all stations.

Spatially related structures in the macrobenthic assemblages along the nearshore-offshore and north-south axes are evident from the MDS and PERMANOVA analyses. The hydrographical and sediment parameters which show variations as a function of the depth, also influence the fauna. However, prominent spatial variations in the structural parameters of macrobenthos are not observed. Fauna at the two stations where low DO values were recorded do not show any sign of impoverishment, as indicated by the total density (2,433–3,090 ind.m⁻²), S (57–66 species), and relatively high H′ values (4.0–4.3).

Despite such high TPH concentrations at a large number of stations, “slightly disturbed” conditions, with either “high” or “good” ecological status indicate that the benthic communities have experienced very low stress in the open waters, 12 years after the 1991 Oil Spill. During the spill, the oil has moved from the coast of Kuwait down to Abu Ali Island in Saudi Arabia along the nearshore area (Fowler et al., 1993; Michel et al., 1993; Sadiq and McCain, 1993; Tawfiq and Olsen, 1993; Readman et al., 1996; Price, 1998). Hence, the residual impact, if any, can be expected mostly at the stations located close to the shore (transect A) from the north down to Abu Ali and nearshore stations of the northern (transect C) and central (transect D) regions. The oil spill has not affected the nearshore area south of Abu Ali (Tawfiq and Olsen, 1993; Readman et al., 1996). Our results indicate that the health of the benthic communities in the nearshore area north of Abu Ali (shown by the average values of AMBI, M-AMBI, S, and H′ are 1.8, 0.72, 64, and 4.63, respectively) is better than that of the nearshore area south of Abu Ali (shown by the values of the same parameters of 2.09, 0.68, 59, and 4.38, respectively) even though the concentration of TPH north of Abu Ali (average of 104 mg kg⁻¹) is higher than that south of Abu Ali (average of 35 mg kg⁻¹). Coles and McCain (1990) have reported the status of benthic community based on samples collected during the 1985–1986 period from three stations in the nearshore area north of Abu Ali. They have reported an average S of 45 species per station from the sand/silt habitat, which is lower than that recorded from the same area in the present study (64 species per station). As the benthos samples collected by Coles and McCain (1990) are not from polluted stations, the comparison confirms the absence of a long-term impact from the 1991 Oil Spill.

Recovery of a marine ecosystem from an oil spill can be defined as the re-establishment of a biological community in which the plants and animal characteristics of the community are present and functioning normally (Kingston, 2002). However, in the recovery assessment it is impossible to say whether the benthic communities that have reestablished after the oil spill are the same or are different from those which would have existed prior to the oil spill. The western Gulf has a large number of oil-related facilities in the form of production wells, platforms, underwater pipelines, refineries etc. Periodically, new installations are built and the existing facilities are upgraded. These oil-related activities have the potential to cause disturbances to the benthic communities due to the physico-chemical alterations of the benthic habitat. As such activities existed prior to the 1991 oil spill, it can be assumed that an “undisturbed condition” with “high′' ecological status for the benthic communities may not have existed in the impacted open waters prior to the spill. In that case, the recovery of the benthic communities to the level of a pristine condition as per the definition of Kingston (2002) may not be possible in the western Gulf as these oil-related activities still continue.

The spatial trends in the structural parameters of the macrobenthos observed in the present study are assumed to be due to the oil-related activities in this region. A great majority of the offshore oil production facilities are situated in the northern region and they decrease toward the south. Accordingly, there is a north-south decrease in the concentration of TPH with a corresponding increase in the oil sensitive organisms such as ampeliscids. Although, relatively lower than the nearshore region, the offshore region also have high concentrations of TPH (Table 1), which have the potential to pose a threat to the marine environment. Despite the slightly higher H′, in the offshore region compared to that of the nearshore region, more disturbed conditions (according to AMBI) and a slightly inferior ecological status (according to M-AMBI) are observed in the offshore region compared to the nearshore region. As the Oil Spill has moved toward south through the nearshore region, the high concentration of TPH in the offshore region cannot be attributed to the 1991 Oil Spill. The high concentration of TPH and the relatively poorer ecological status of the benthic communities in the offshore region may be due to the effects of the offshore platforms and the increased vessel traffic.

TOC values directly proportional to the concentration of TPH have been reported in some parts of the world after oil spill incidents (for example, Wattayakorn and Rungsupa, 2011). However, in most cases, such relationships have not been observed in the Gulf (Abdali et al., 1993; Al-Ghadban et al., 1994; Al-Lihaibi and Al-Omran, 1996; Massoud et al., 1996; Al-Lihaibi and Ghazi, 1997; de Mora et al., 2010) and the conclusions of the present study is in accord with the conclusions of the previous studies. The concentration of TOC can be considered as an indicator of hydrocarbon pollution only when the concentration of TPH is very high (Al-Ghadban et al., 1994).

This study clearly shows the presence of residual effects of the 1991 Oil Spill in the sheltered inner bays of the Gulf, where the grain size of the sediment is more or less the same as that in other areas. Of the five locations in the inner bays evaluated in this study, the TPH concentration is very high in the MTBS, specifically at the station located in Balbol. In a study conducted in 2003, a TPH concentration of 1,000 mg kg⁻¹ was recorded in the surface sediments of an impacted station, indicating that the situation was at an alarming level even 12 years after the oil spill. Sauer et al. (1993) has reported the occurrence of sunken oil from the 1991 Oil Spill at a depth of 0–5 cm from the sediment surface inside the sheltered inner bays. A core sample from Balbol Bay analyzed along with the present study shows a concentration of 38,000 mg kg⁻¹ at a depth of about 15 cm from the seabed. This indicates that the rate of accumulation of oil in the sediments was rapid in the sheltered inner bays due to the absence of strong current or wave action. The post-1991 Oil Spill studies have reported higher TPH values in the sediments of the MTBS and AABS compared to those in the sediments in the open water areas (Fowler et al., 1993; Michel et al., 1993; Readman et al., 1996). About seven months after the Spill, Fowler et al. (1993) recorded sediment hydrocarbon concentrations up to 1,140 mg kg⁻¹ from the open water regions and up to 1,400 mg kg⁻¹ from the sheltered inner bays. Michel et al. (1993) also reported sediment hydrocarbon concentrations up to 900 mg kg⁻¹ from the inner bays of the AABS 1 year after the spill. Even though a consistent reduction in the contamination of the bays has been reported in the subsequent years, the recovery has been slow. Readman et al. (1996) reported an increase of the contamination at certain sites on the impacted coasts in the sheltered bays when the same sites were surveyed during the months of August in 1991, 1992, and 1993. They attributed this to the patchiness in the distribution and the re-deposition of a remobilized oil pockets from the original contamination event.

According to the results of the present study, the four stations from the inner bays can be categorized as having a “moderate-poor” ecological status. Even though, the highest concentration of TPH was observed at the station in Balbol, the other two locations in the MTBS, viz., Tanajib and Manifa, can be categorized as having a “good-higH′' ecological status. On the other hand, both sites studied in the AABS, Mussallamiya and Mardomma, can be categorize as having a “moderate-poor” status, indicating a higher stress on the benthic community compared to the MTBS. However, an accurate picture will be revealed only if a large number of samples are collected from both MTBS and AABS and evaluated.

In this study samples were not collected from the intertidal region. However, high levels of polyaromatic hydrocarbons (PAHs) along the shoreline of the sheltered bays of Saudi Arabia, originating from the 1991 Oil Spill, have been reported (Bejarano and Michel, 2010). These authors note that the concentration of the sum of PAHs in these samples are similar to the levels found in the sheltered areas 1year after the oil spill (100,000–300,000 ng g⁻¹). Persisting high concentration of PAHs indicates that the physical removal or microbial degradation of the oil in these heavily oiled intertidal sediments is very low. As alarming levels of oil are present in the intertidal (Bejarano and Michel, 2010) and subtidal (present study) zones more than even a decade after the spill, the complete recovery may take decades more. Chances of the re-distribution of the sunken oil due to dredging/trenching activities are also present. This area, being in close proximity to the oil fields, may host future facilities. Any activities that require dredging/trenching can cause the re-introduction of the oil into the water column, and thereby the re-distribution to the surface of the seabed. This can impact the resident benthic community and cause further delays in the recovery of the inner bay habitats from the 1991 Oil Spill.

The indices AMBI and M-AMBI have not been commonly used in studies of the Gulf waters. These two indices were used in this study because of their extensive use in different habitats worldwide for assessing the ecological status of benthic communities (Borja et al., 2012). The worldwide acceptance of these indices is due to the fact that the benthic status of a given area determined by AMBI and M-AMBI is highly correlated with the benthic status assigned by professional judgment (Borja et al., 2014). However, one difficulty of using such indices in the Gulf is the selection of reference conditions. Commonly, reference conditions are derived from: (i) comparison with an existing “pristine”site; (ii) past data and information; (iii) numerical models; and (iv) best professional judgment (Borja et al., 2004, 2012; Muxika et al., 2007; Forchino et al., 2011). Due to the widespread oil-related activities in the Gulf, and the likelihood of the resulting stress on the marine environment (Bishop, 2002; Khan et al., 2002; Munawar et al., 2002; Jones et al., 2007; Khan, 2007), finding pristine sites is a great challenge. The Gulf is a naturally stressed water body with high salinity and high seasonal temperature variations. In such environments, generally, the biota will be living at their upper tolerance levels and the impacts of any human pressures can be increased many fold compared to their effects in a normal environment. Considering all these facts, Sheppard et al. (2010) recently categorized the Gulf as a young sea in decline. In this region, the majority of the historical data have been generated from a multitude of Environmental Impact Assessment studies. One drawback of this type of data in serving as a good baseline is that, in most cases, the baseline chosen is the already severely degraded marine ecosystem (Sheppard et al., 2010). Hence, a comparison with the previous data is generally of no use. In such circumstances and in the absence of numerical models, expert judgment may be the best option to determine the reference conditions (Borja et al., 2014). In the present study, the approach used for establishing the reference condition is by selecting the highest richness and diversity values determined and increasing them by 15%. This has been considered to be a reasonable approach (Borja and Tunberg, 2011; Forchino et al., 2011; Paganelli et al., 2011). Based on this reference condition, about 25% of the stations were categorized as having a “higH′' ecological status. These stations have an average S value of 82 species and H′ of 5.1. These values are expected to be used as reference data in the future impact assessments on benthic communities in the Gulf.

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