Ecological status assessment and non-indigenous species in industrial and fishing harbours of the Gulf of Gabès (central Mediterranean Sea)

Port Biological Baseline Surveys (PBBS) are standardized surveys of the indigenous and non-indigenous marine biodiversity within harbour activities. They provide a baseline for monitoring changes in the structure and function of harbour communities. This study conducted in 12 fishing and industrial harbours from January to December 2018 was the first initiative of a Port Baseline Survey aimed to assess the impact of biological invasions in harbours of the Gulf of Gabès (GG), Tunisia. A total of 174 macrobenthos species were recorded, belonging to eight phyla, with a dominance of crustaceans (32%), molluscs (31%) and polychaetes (20%). Among these species, 57 were non-indigenous species (NIS) for Tunisian waters, while 27 species were recorded for the first time in GG harbours, and three decapods (Dyspanopeus sayi, Hippolyte prideauxiana and Pilumnus minutus) and one amphipod (Hamimaera hamigera) were newly recorded from Tunisian waters. Two main categories of harbours are distinguished according to their macrobenthic communities and environmental conditions. The industrial harbours yield higher richness and abundance of NIS than the fishing harbours. The ALEX metric is used to evaluate the biological invasion status of the Gulf of Gabès harbours and shows that their status ranges from unaffected in fishing harbours to extremely affected in industrial harbours. Three biotic indices (AMBI, BO2A and BENTIX) are applied to assess the ecological status of harbours, which varies from moderate to good. ALEX and the other biotic indices are significantly correlated with harbour characteristics, maritime traffic and edaphic factors (organic matter and chemical contamination). The present study provides a data baseline for the implementation of environmental policies and management plans in the future.


Introduction
Harbours are one of the most disturbed coastal ecosystems due to intensive anthropogenic pressures (shipping activities, pollution and dredging: Chan et al. 2016;Dauvin et al. 2017;Chatzinikolaou et al. 2018;Romanelli et al. 2019;Dimitriou et al. 2020). In addition, harbours are a major pathway to non-indigenous species (NIS) due to the maritime activities and artificial structures, such as docks and floating pontoons found in harbours, provide suitable habitats to host opportunistic fouling species and therefore facilitate and accelerate the introduction of NIS (Ardura et al. 2015;López-Legentil et al. 2015;Tempesti et al. 2020). In recent years, several studies have characterized the biota present in harbours to allow proper monitoring programmes of these marine environments. However, this strategic step needs to be preceded by a systematic and thorough study of the biological components of the ecosystems concerned (Mandal and Harkantra 2013;Lehtiniemi et al. 2015). Without baseline data on indigenous habitats, it is practically difficult to apply appropriate management protocols of marine biological invasions (Lehtiniemi et al. 2015).
In the Mediterranean Sea, harbours are the primary regions where NIS settle, because adequate habitat diversity is found in these marine ecosystems (Awad et al. 2014;Tempesti et al. Responsible editor: Lotfi Aleya * Nawfel Mosbahi nawfelmosbahi@hotmail.fr 2020). The successful hosting of NIS largely depends on the compatibility between the settling of species in the new habitats and favourable conditions for establishment, including the availability of trophic resources and the absence of predators (Hayden et al. 2009). In recent years, marine bioinvasions have become one of the greatest global threats to the diversity and integrity of indigenous communities (Molnar et al. 2008;Galil et al. 2018). These invasions are considered to be drivers of irreversible impacts in host environments while also affecting the diversity and the stability of native habitats (Spagnolo et al. 2019). The Mediterranean Sea is one of the regions of the world ocean showing large invasive effects due to the impact of NIS, with the increasing number of introductions representing one of the greatest threats to Mediterranean biodiversity (i.e. Zenetos et al. 2010Zenetos et al. , 2012Zenetos et al. , 2017Zenetos et al. , 2018Galil et al. 2018;Dragičević et al. 2019;Bailey et al. 2020). In fact, the number of NIS recorded in the Mediterranean Sea has been increasing, which affects habitats and ecosystem functioning (Katsanevakis et al. , 2016. With more than 1400 km of coastline, Tunisia occupies a key position located at the crossroads between the western and central parts of the Mediterranean Sea. The Sicilian strait is the passageway from south to north and from east to west and is crucial in the analysis of the spread of NIS introduced into the Mediterranean Sea. A total of 136 marine NIS have been recorded in Tunisian waters, with 60 records from northern coasts, and 76 from central and southern coasts (Ounifi-Ben Amor et al. 2016). In 2019, Ben Souissi et al. (2019) reported, for the first time in Tunisia, two annelids, ten crustaceans, two bryozoans and one hydrozoan NIS. Moreover, a total of 27 non-indigenous marine macrophytes have been recorded in coastal Tunisian waters by Sghaier et al. (2016). The rapid spread of NIS in Tunisian waters is linked to ballast water discharge and biofouling in relation to international maritime transport. The Gulf of Gabès (southern Tunisia) includes the highest number of harbours in Tunisia; it is considered one of the most vulnerable aquatic ecosystems in the Mediterranean Sea, due to biological invasions (Ounifi-Ben Amor et al. 2016). Moreover, intensive anthropogenic pressures such as trawling practices, diverse pollution (e.g. phosphogypsum inputs, industrial and urban wastes) and shipping activities cause the deterioration of Posidonia oceanica meadows and the continual decline of marine fisheries resources (Hattab et al. 2013;Boudaya et al. 2019;Mosbahi et al. 2019).
Due to their position at the sediment-water interface, and their relatively long life spans and sedentary habit, macrobenthic invertebrates have been largely used as bioindicators for aquatic monitoring because they respond rapidly to anthropogenic and natural stressors (Reiss and Kroncke 2005;Blanchet et al. 2008;Dauvin et al. 2012Dauvin et al. , 2017. In order to survey the biodiversity of harbours areas, many studies have focused on macrobenthic communities worldwide (Ingole et al. 2009;Mandal and Harkantra 2013;Dauvin et al. 2017;Chatzinikolaou et al. 2018).
Despite this research effort, biological and ecological information on harbour areas is generally scarce, which is problematic as Port Biological Baseline Surveys (PBBS) can be used to provide data against which future changes in the structure and function of marine communities can be measured (Awad et al. 2014;Spagnolo et al. 2019). In Mediterranean harbours, the majority of studies have reported the prevalence of NIS on hard bottoms, including artificial structures in docks and marinas (Spagnolo et al. 2014(Spagnolo et al. , 2019. However, there are few quantitative studies designed to examine soft-bottom communities and the impacts of NIS in soft sediment habitats (Travizia et al. 2019). In Tunisia, the majority of benthic macrofauna studies have focused on macrobenthic communities of shallow coastal areas (Khedhri et al. 2015;Fersi et al. 2018;Boudaya et al. 2019;Mosbahi et al. 2017Mosbahi et al. , 2019Mosbahi et al. , 2020. The few studies carried out on macrobenthic fauna in Tunisia harbours (Aloui-Bejaoui and Afli, 2012;Chatzinikolaou et al. 2018;Chebaane et al. 2019;Dimitriou et al. 2020) have focused on the structure of macrobenthic faunal communities or the identification of NIS in harbour ecosystems, and the ecological quality status (EQs) of these ecosystems is not yet well studied. The aims of the present study covering the majority of the Gulf of Gabès harbours are as follows: (1) to describe their soft-bottom macrobenthic communities and (2) to assess their EQs linked to different levels of anthropogenic stress and biological invasions, so as to determine the relationship between these polluted zones and the establishment of NIS. Furthermore, this article presents baseline data for determining ecological status and bio-invasion impacts in Tunisian harbours and could serve as background information for harbour authorities and environmental managers for the design and implementation of a harmonized monitoring plan in the central Mediterranean.

Study area
The Gulf of Gabès (GG) is located in the central part of the Mediterranean Sea, delimited by the south-eastern coast of Tunisia (35-33°N and 10-13.5°E) (Fig. 1). The GG is a shallow basin with waters no deeper than 50 m as far as 110 km away from the coast, increasing to 200 m deep at 400 km from the shore. Its climate is dry (average annual precipitation: 210 mm year −1 ) and sunny with strong easterly winds. The tide is semidiurnal, with a maximum range of about 2.3 m near the coast around the Kneiss Islands (Sammari et al. 2006).
The GG supports one of the most productive ecosystems in the Mediterranean Sea, not only associated with significant economic importance (contributing about 65% of the national fish production in Tunisia) but also comprising 18 harbours with numerous activities (fishing, commercial, marinas), and 4 commercial harbours (Zarzis, Gabès, Skhira and Sfax) open to international shipping (Hattab et al. 2013;Ben Rejeb-Jenhani et al. 2018).
In spite of its importance for fisheries and high natural heritage value, the GG is considered one of the most heavily polluted Mediterranean areas owing to the presence of well-developed coastal cities generating various sources of industrial, agricultural and domestic contaminants (El Zrelli et al. 2018;Mosbahi et al. 2019). Due to its position in the central part of the Mediterranean, the GG is subject to increasing biological invasions through the continued introduction of NIS not only from the Red Sea via the Suez Canal as well as from Atlantic waters, but also from human activities such as marine traffic and aquaculture (Ounifi-Ben Amor et al. 2016).

Sampling design and sample treatments
During 2018, surveys were conducted seasonally (January (winter), April (spring), July (summer) and October (autumn) 2018) in 12 GG harbours (Fig. 1). During each seasonal survey, to avoid temporal effects, the sampling campaign was carried out in each harbour during the same period. Two main categories of harbours can be distinguished: three industrial harbours and nine fishing harbours ( Fig. 1; Table 7 in Appendix). For each season, three stations were sampled in each harbour to cover the entire harbour basin from its entrance to the inner part, making a total of 36 stations visited four times during 2018. The station positions were accurately determined using a GPS (Global Positioning System, WGS84). For each station, five replicate samples were taken using a Van Veen grab sampling an area of about 0.05 m 2 (for each replicate) and penetrating approximately 0.1m into the sediment. Four replicate samples were preserved for biological analysis covering a total surface area of 0.2 m 2 . The last sample was used to determine sediment characteristics. Each biological sample was sieved on a 1 mm mesh, and the retained fractions were conserved in 5% formaldehyde saline solution. Some collected samples contained a hard fraction including pebbles and shells, which were preserved and considered for further studies. In the laboratory, after staining with Rose Bengal, samples were sorted; the macroinvertebrates were identified to species level and then counted. The nomenclature of macrobenthic species was updated according to the World Register of Marine Species (WoRMS, last accessed on 07 May 2021).
In addition, temperature (T°C), salinity (Sal), pH, dissolved oxygen concentration (mg.L −1 ) and transparency (m) were measured using a thermometer (WTWLF 196), a salinometer (WTW LF 196), a pH meter (WTW 3110), an WTW oximeter and a Secchi disc, respectively, at each station and for each seasonal campaign. The chlorophyll a concentration (Chl a) was determined on 1 L of marine water, which was collected (at each harbour/season) and transported in the dark and at low temperature to the laboratory and then filtered on GFC filters and extracted using 100% acetone. The absorbance was measured with a spectrophotometer at 630 nm, 647 nm, 664 nm and 750 nm, and the concentration was estimated according to Rodier et al. (1996).
After eliminating the coarsest fraction, the sediment characteristics of the samples were determined by sieving around 1 kg of dry sediment on 2, 1, 0.5, 0.25, 0.125 and 0.063 mm meshes and then weighing each sediment fraction. The mud content (% mud) is expressed as dry weight of <63 μm sediment in relation to the total weight, and the organic matter content (% OM) as weight loss on ignition at 500°C for 4h, using 100 g samples of dry sediment (24h at 60°C). Contents of heavy metals (Zn, Cd and Pb), phosphorus (P), fluorine (F) and nitrogen (N) were determined after digesting the powder sample in aqua regia at 95°C, and analysis by inductively coupled plasma-atomic emission spectrometry (ICP-AES) and mass spectrometry (ICP-MS) (Yoshida et al. 2002).

Data preparation
For each station, we determined the most common macrofauna biodiversity indices defining the EQ of a given station or benthic community. Data were used to calculate the species abundance (A, number of individuals estimated per m 2 ), the taxonomic richness (S, number of species 0.2 m 2 ), the Shannon-Weaver diversity index (H′) in log 2 (Shannon and Weaver 1963) and Pielou's evenness (J') (Pielou 1966) for each station and sampling period. Data analysis was carried out using version 6 of the PRIMER® (Plymouth Routines in Multivariate Ecological Research) software package (Clarke and Gorley 2006).
For the assessment of environmental quality, we make use of three biotic indices: AMBI (AZTI Marine Biotic Index, Borja et al. 2000), BENTIX (Simboura and Zenetos 2002) and BO2A (Benthic Opportunistic Annelids Amphipods index; Dauvin et al. 2016;Dauvin 2018). AMBI and BENTIX were calculated using the software from AZTI and then applied to estimate the proportions of five ecological groups (using the species list published by the AZTI web site on 30 June 2019 http://ambi.azti.es/). The BO2A (Benthic Opportunistic Annelids Amphipods) index was calculated as log10 of the ratio of frequencies for opportunistic annelids and amphipods: i.e. the total number of opportunistic annelids and total number of amphipods +1 divided by the overall abundance counted in a sample (see Dauvin et al. 2016).
Furthermore, ALEX (ALien biotic indEX, Çinar and Bakir 2014) is used to evaluate the impact of NIS on indigenous assemblages (Piazzi et al. 2020). This index is based on the abundance percentages of different groups defined from the level of species establishment and invasiveness within each sample (GI, indigenous species; GII, casual NIS; GIII, established NIS; GIV, invasive NIS). Species are assigned to ALEX groups according to observations of their abundance in the area and data from the relevant literature (see Çinar and Bakir 2014).
The identified species were classified into six trophic groups according to the categories used by many authors (e.g. Afli et al. 2008;Jumars et al. 2014;Mosbahi et al. 2017Mosbahi et al. , 2019Mosbahi et al. , 2020. Non-selective deposit feeders (NSDF) are burrowers that ingest the sediment from which they take their food; selective deposit feeders (SDF) ingest organic particles from the sediment surface); suspension feeders (SF) extract suspended food in the water column; carnivores (C) feed on preys and herbivores/grazers (HG) feed on macrophytes and microalgae) ( Table 7 in Appendix).

Statistical techniques
The original data consists of a 'stations ×species' matrix, which is obtained after averaging data for the four seasons. Macrobenthic abundances are firstly square-root transformed to minimize the influence of the most dominant species. Cluster analysis and non-metric multidimensional scaling (n-MDS) based on the Bray-Curtis similarity were applied to visually assess differences in benthic assemblages between stations of the studied harbours. SIMilarity PERcentages (SIMPER) tests were performed to determine which species contributed most to within-group similarity. One-way ANOSIM (ANalysis Of SIMilarities) permutation test is used to assess if the assemblages differences between sampling stations are statistically significant. The BIOENV procedure was used to identify which combination of environmental variables best explains the differences in macrofauna distribution. These multivariate analyses were applied by PRIMER® -v6 software (Clarke and Gorley 2006). Spearman's correlation coefficients were used to determine the existing relationships between the biotic indicators and the environmental factors of the GG harbours. Spearman's correlation coefficients were calculated by the SPSS Statistics 20 software. Three-way analyses of variance (ANOVAs) were used to test spatial and the seasonal variability of the environmental factors (e.g. OM, salinity, temperature) and the benthic diversity (e.g. taxonomic richness, abundance, trophic groups) and biotic indices (e.g. AMBI, BENTIX) between station, harbour and season. The values of the three parameters (station + season + harbour) were considered independent variables, and each response variable (e.g. environmental variable, benthic and biotic indices) was considered dependent variable. Tukeyadjusted post hoc comparisons were performed to compare the categories for variables included in the regression model. The three-way ANOVAs were performed using R software, and a significance level of < 0.05 was considered in all tests.

Environmental characteristics
The organic matter (OM) and elemental contents show significant variation between the studied harbours (three-way ANOVA; F= 18.06; p< 0.05) for OM and for trace metals p< 0.01 in all cases. The higher values of OM and metal contamination are recorded in the industrial harbours (GAI, SKI and SFI) ( Table 1).
The chemical and physical parameters show seasonal variations. The temperature is maximal in summer (33 C°± 0.1 in SKI and HSF), while minimum values are recorded in winter (11 C°± 0.02 in HSF and SFF; three-way ANOVA; F= 10.08; p< 0.01). Salinity shows a seasonal variation (F= 12.18; p< 0.01), and the highest values are enregistered in summer and autumn in SFF (41.5± 0.05) and GAF (41.2± 0.01). pH shows a spatial variation during the sampling period, with low values in industrial harbours (F= 4.12; p< 0.01). Dissolved oxygen shows a significant spatial (F= 3.08; p< 0.01) and seasonal variability (F= 4.16; p< 0.01). Transparency only displays spatial changes (F= 18.04; p< 0.05), with low values being recorded in industrial harbours (SFI, SKI and GAI) during the all-sampling campaign ( Table 8 in Appendix).

Macrobenthic community composition
A total of 16,102 individuals were identified, belonging to 174 macrobenthic species, 109 families and eight zoological groups unequally distributed among the sampling sites. Crustaceans are the dominant group in terms of number of species (55 species; 32 % of total number of species), followed by molluscs (53 species; 31%) and annelid polychaetes (35 species; 20%). The other five phyla taken together account for 17 % of the total number of species (Table 9 in Appendix). As regards the abundance, crustaceans (34% of the total abundance) and polychaetes groups (30%) have the highest number of individuals in the all prospected harbours.
The taxonomic richness differs significantly between harbour site (three-way ANOVA; F= 14.02; p< 0.05) (Fig. 2). The mean taxonomic richness varies from 26.5 in the Zabboussa fishing harbour (ZBF) to 55 species in the Skhira industrial harbour (SKI). As well, the abundance of macrobenthic species shows a significant spatial difference between the studied harbours (three-way ANOVA; F= 112.08; p< 0.05). The mean abundance varies from 6157 ind.m −2 in the fishing harbour of Skhira (SKF) to 9968 ind.m −2 in ZBF. The Shannon's index H' values range from 2.02 (at SKI) to 3.88 bits.ind −1 (at HSF), and values of Pielou's J' vary from 0.48 (at SKI) to 0.82 (at HSF and MAF).
The species number and abundance of the NIS shows significant differences between the GG harbours. The higher number of NIS (three-way ANOVA; F= 13.1; p< 0.05 for number of species) and abundance (three-way ANOVA; F=62.06; p< 0.05 for abundance of NIS) are recorded in the industrial harbours (Gabès, Skhira and Sfax) compared to the fishing harbours (Fig. 2). Three-way ANOVA highlights a significant influence of harbour site (F= 10.112; r 2 =0.062; p< 0.05) and harbour type (F= 2.042; r 2 =0.048; p< 0.05) on the structure of benthic communities (number of species and abundance) between harbours. The harbour macrofauna communities are strongly dominated by carnivores (43%), detritus feeders (32%) and selective deposit feeders (22%). The diversity and abundance of the trophic groups varies significantly between harbour sites (p < 0.05). During the whole year, the benthic communities in industrial harbours are dominated by carnivore groups. On contrary, the trophic groups show a seasonal variation in the fishing harbours (three-way ANOVA; F=28.11; p < 0.05).

Spatial patterns of the macrofauna assemblages
Cluster analysis based on four replicates from the three stations at each harbour highlights a clear separation between internal and external stations (not shown here). Moreover, the cluster analysis and n-MDS ordination allows us to separate the 36 stations (three stations in 12 harbours) into two main groups at a similarity level of 30 %: the first group includes the nine stations sampled in industrial harbours (Gabès, Skhira and Sfax), while the second group corresponds to stations sampled in fishing harbours (Fig. 3). SIMPER analysis shows that the second group (61.50% contribution  (Table 2). ANOSIM analysis reveals a significant difference between the two harbours groups (R ANOSIM = 0.38; p < 0.01).
The BIOENV procedure indicates that the variations in macrofaunal distribution in harbours of the Gulf of Gabès can be explained by a combination of several variables (Table 3, correlation= 0.60). These correspond to sediment characteristics (OM, mud content and pollutant contamination) and harbour characteristics (depth and surface-area). Organic matter content individually shows the highest correlation with species distribution (correlation= 0.24).

Ecological quality status
The average biotic indices AMBI, BO2A and BENTIX differ significantly between sampling sites (three-way ANOVA; p<0.05), with AMBI values ranging from 2.24 to 4.83, BO2A from 0.011 (at MAF) to 0.182 (at SKI) and BENTIX from 2.40 to 4.22. In the same way, ALEX values are statistically different between harbours (three-way ANOVA; F= 16.2; p< 0.05), with values varying from 0.86 (MAF) to 4.20 (SKI). The three industrial harbours (Sfax, Skhira and Gabès) appear to be have moderate ES and are heavily to extremely affected by biological invasion, being strongly dominated by tolerant (EGIII) and opportunistic polychaete species. Contrariwise, biotic indices classify the majority of fishing harbours in good ecological status, with low ALEX values indicating that these harbours are unaffected to moderately affected by biological invasion (Table 4).
The Spearman's rank correlation coefficients between environmental variables and ALEX, AMBI, BENTIX and Shannon-Weaver diversity show that ALEX is positively and significantly correlated with harbour surface area (p r = 0.04), temperature (p r = 0.46), salinity (p r = 0.03), organic matter (p r = 0.018), sand (p r = 0.22) and mud percentage (p r = 0.38) as well as other environmental factors such as chemical contamination involving Pb, Cd, Zn, phosphorus, fluorine and nitrogen; on the other hand, the results show that ALEX is negatively and significantly correlated with depth (p r = −0.32), transparency (p r = −0.25) and dissolved oxygen (p r = −0.11).
The ALEX values (p r = 0.42), number of NIS (p r = 0.26), the number of individuals of NIS (p r = 0.36) and the ratio between the number of individuals of NIS and indigenous species (p r = 0.46) significantly increase with harbour surface area. AMBI and BENTIX are negatively correlated with harbour surface (p r = −0.44), organic matter content (p r = −0.33), mud (p r = −0.32), phosphorus (p r = −0.36) and cadmium (p r = −0.12). BO2A is positively correlated with organic matter (p r = 0.42) and mud (p r = 0.58) and negatively correlated with trace metals such as Cd (p r = −0.14) and Pb (p r = −0.05) ( Table 5).

Discussion
Semi-enclosed coastal areas with low-energy hydrodynamics, harbours are considered disturbed marine environments where different sources of pollution, combined with global changes, result in complex ecological relations and different processes which often mask the distinction between environmental and natural impacts (Kapsimalis et al. 2014;Chan et al. 2016;Chatzinikolaou et al. 2018;Tempesti et al. 2020). The present study investigates the benthic macrofauna structure, the distribution of benthic assemblages and EQ assessment by the application of some functional and biotic indices in response to multiple stressors in the GG harbour ecosystems.
In this study, the species inventory comprises 29 macrobenthic species recorded for the first time in the GG, including four species new to Tunisian waters that have already been described in the Mediterranean Sea. Therefore, taking into account the 136 species listed by Ounifi-Ben Amor et al. (2016) and the latest inventory added by Ben Souissi et al. (2017Souissi et al. ( , 2019, the number of NIS recorded in Tunisian waters now reaches a total of 157. The number of NIS recorded in Tunisia has been increasing during the last few decades. The introduction of new NIS could be related to certain economic activities such as marine traffic, fisheries, aquaculture and tourism (Streftaris and Zenetos 2006;Ben Souissi et al. 2014;Ounifi-Ben Amor et al. 2016).
In terms of taxonomic composition, the structure of communities in GG harbours is similar to that observed in other harbours of the Mediterranean ecosystem, being dominated mainly by molluscs, crustaceans and polychaetes Chatzinikolaou et al. 2018;Travizia et al. 2019;Dimitriou et al. 2020). The number of species (174)  However, this diversity in Gabès harbours appears low compared to the biodiversity recorded recently by Chatzinikolaou et al. (2018) in the El Kantaoui harbour (Tunisia) (211 species), as well as by Travizia et al. (2019) in the Bari harbour (Italy)(224 species) and by Dauvin et al. (2017) for 10 Algerian harbours (847 species). This difference can be attributed to the sampling strategy adopted for each zone (sampling effort, time and methods of sampling), the number of sampled harbours and the specificity of the harbour environment (e.g. environmental conditions, anthropogenic stressors) ( Table 6).
The macrobenthic fauna composition and benthic diversity in GG harbours show spatial variations between the inner basins and the entrances. These differences could be explained by changes in the levels of organic and metal contamination, which vary according to the hydrodynamic regime corresponding to the water exchange between the inner basins and the entrances of the harbours (Grimes 2010;Dauvin et al. 2017). Equally, the spatial distribution of macrobenthic communities in GG harbours reflects the existence of two distinct macrofaunal assemblages. The first assemblage corresponds to the stations sampled in industrial harbours, which appear severely and extremely affected by NIS such as Pinctada imbricata radiata, Cerithium scabridium, Bursatella leachii, Paracerceis sculpta, Cymadusa filosa, Portunus segnis and Libinia dubia. The second assemblage regroups stations sampled in fishing harbours dominated by indigenous opportunistic species such as Gammarus insensibilis, Maera hirondellei, Dexamine spiniventris, Glycera tridactyla, Perinereis cultifera, Nephtys hombergii and Heteromastus filiformis; this assemblage is unaffected or only slightly affected by biological invasions. The difference between the two assemblages could be explained by the difference of environmental parameters between the GG harbours and especially their sediment type, organic matter content and contamination by pollutants. In fact, industrial harbours are characterized by muddy sediment richer in organic matter and heavy metals. Several previous studies in harbours have indicated that sediment features, water depth, hydrodynamics and pollution (i.e. hydrocarbons) are significant environmental factors which are associated with the clustering of samples from harbours in terms of abundance (Blanchard et al. 2002;Dauvin et al. 2017;Chatzinikolaou et al. 2018). The effect of the frequency of marine traffic and environmental disturbances is highlighted by the high number of NIS recorded in the three industrial harbours of the GG (Gabès, Skhira and Sfax). In the Mediterranean Sea, maritime traffic in commercial and industrial harbours is the main pathway evoked for the introduction of NIS (Molnar et al. 2008;Nunes et al. 2014). The introduced species are attached to the submerged parts of vessel hulls, while organisms transported by water or sediment are contained in ballast water discharges (Simkanin et al. 2009;David and Gollasch 2015;Travizia et al. 2019). For these reasons, industrial harbours (marinas as well) are considered coastal environments especially susceptible and vulnerable to the establishment of NIS due to their favourable abiotic and biotic conditions (i.e. trophic availability, less competition and predation) (Çinar et al. 2012;Awad et al. 2014). Several authors have shown that harbours do not represent good habitats for indigenous communities and are instead populated by highly tolerant and opportunistic introduced species. Likewise, shipping is considered the major pathway for NIS introduction in Mediterranean harbours (López-Legentil et al. 2015;Ulman et al. 2019;Tempesti et al. 2020). These authors also found a significant relationship between harbour size and species richness, indicating that larger harbours tend to contain higher numbers of species associated with higher abundances. This is in agreement with the results of our study, which show that the major harbours of the GG (industrial harbours) are those which shelter more species, especially the NIS. Recent studies also show that recreational boating and fishing vessels might play an important role, especially in the spreading of NIS at the small to medium spatial scale (Ferrario et al. 2015;López-Legentil et al. 2015;Ulman et al. 2017Ulman et al. , 2019Tempesti et al. 2020).
The GG harbours are strongly dominated by carnivores, detritus-feeders and selective deposit feeders related to the availability of trophic resources. Several authors have mentioned a significant relationship between benthic trophic structure, sediment contaminants and environmental variables (Guerra García and García Gómez 2005; Chatzinikolaou . Therefore, changes in trophic structure could be used as an indicator of disturbance. The reduction in trophic complexity is associated with organic-rich and chemically contaminated sediments in polluted environments (i.e. harbours), where the benthic assemblages are dominated by opportunistic species (Rakocinski et al. 2000;Putro 2009;Dauvin et al. 2017).
In Mediterranean harbours, several multivariate biotic indices such as AMBI (Borja et al. 2000), BENTIX (Simboura and Zenetos 2002), BO2A (Dauvin et al. 2016) and BQI (Leonardsson et al. 2015) have been recommended and extensively used to assess and monitor EQs (Riera et al. 2011;Sany et al. 2015;Tomassetti et al. 2016;Dimitriou et al. 2020). In the present study, three biotic indices (AMBI, BO2A and BENTIX) are used to assess EQs, allowing us to classify the GG harbours as having poor to good ecological status. The first group includes stations sampled in fishing harbours, which appear to be in good ecological status. On the other hand, the second group is made up of industrial harbours showing a moderate ecological status, dominated by pollution-tolerant species and opportunistic polychaetes species indicative of stressed environments, such as C. capitata; S. (Scolelepis) squamata, M. fuliginous and C. cirratus. These latter species form part of normal benthic communities found in the open sea, colonizing fine sediment habitats, which could be associated with pollution indicator species showing a slightly polluted system. These opportunist species are also observed in other polluted harbours such as Marseilles along the French Mediterranean coast (Bellan et al. 1980) and the Bethioua and Djendjen harbours in Algeria ). Both of these south Mediterranean harbours are strongly affected by diverse anthropogenic pressures which mean that these environments are subject to heavy metal and hydrocarbon pollution, as well as accumulation of excessive OM, the release of warm waters from power plants and the input of nutrients (Bellan et al. 1980;Dauvin et al. 2017). The relationships between biotic indices and environmental stressors of harbours ecosystems have been studied by many ecologists, who show a highly significant negative correlation between functional diversity, biotic indices and pollution proxies (i.e. chemical pollutants, sediment contamination and plastic litter) (Cole et al. 2011;D'Alessandro et al. 2018aD'Alessandro et al. , 2018bDimitriou et al. 2020). In fact, the lower index values recorded by D' Alessandro et al. (2020) in Maltese harbours suggest a highly disturbed community, by virtue of a heterogeneously occupied functional space with a high level of pollution.
The ALEX index is used here to describe the bio-invasion status in GG harbours, showing that the environmental status of benthic habitats varies from unaffected to extremely affected. The high ALEX score recorded in the industrial harbours shows that some areas in the GG have been modified by   (Piazzi et al. 2018(Piazzi et al. , 2020Tempesti et al. 2020). The correlations between ALEX values and environmental variables indicate that NIS largely prefer shallow water habitats with a high content of organic matter and mud, along with enrichment in heavy metals (Pb, Cd and Zn) and certain other elements such as phosphorus, fluorine and nitrogen. This pattern of enrichment had been observed by Çinar and Bakir (2014) on the coasts of Turkey, where NIS are found coexisting in disturbed environments. It is generally accepted that the majority of Lessepsian migrants have managed to colonize the shallow water habitats of the eastern Mediterranean, while the  shallowness of the canal has acted as a filter for the deep-water biota of the Red Sea (Por 1978;Çinar et al. 2011;Çinar and Bakir 2014). ALEX is positively correlated with the surfacearea of GG harbours. In fact, the high number of NIS recorded in the industrial GG harbours seems to be related to shipping activities and associated biosecurity risks, which include hull biofouling, standing waters from commercial vessels (e.g. in ballast tanks, bilges, anchor chains and engine cooling systems) (Minchin et al. 2009;Lawrence and Cordell 2010;Ferrario et al. 2017;Ulman et al. 2019). The industrial harbours of the GG represent one of the most important economic activities of Tunisia; in 2017, more than 1400 international vessels from several destinations around the world visit the three industrial harbours in the GG (OMMP 2018). These results indicate that marine traffic is a main vector responsible for the introduction of NIS into central Mediterranean harbours.

Conclusion and future perspectives
The present overview provides a valuable baseline database of benthic diversity in the GG harbours, covering a part of the central Mediterranean Sea that is still insufficiently studied. Different environmental parameters affect  Two groups of stations are identified by cluster analysis at a 30% similarity level to the range of environmental conditions found in each studied harbour (sediment type, harbour activity, pollution level). Finally, the results of this study suggest that indigenous biota require some protection against biological invasions in harbour ecosystems. Since these harbours offer hubs for invaders (Floerl et al. 2009), they should be priority targets for control actions through the setting up of Biological Invasions Risk Management Programmes in the GG harbours (Latombe et al. 2017), considering that these areas are very sensitive to anthropogenic pressures and biological invasions. Such programmes should be operated by local port management authorities, local fishermen and research committees to file and implement a monitoring strategy for future management actions. For future research, it appears very important to supplement this work by studies on the macrobenthic species of harbour hard substrates since these habitats are known to be colonized by numerous NIS (López-Legentil et al. 2015;Piazzi et al. 2020).