The spatial distribution of the subtidal benthic macrofauna and its relationship with environmental factors using geostatistical tools: a case study in Trapandé Bay, southern Brazil

Modeling the distribution patterns of the estuarine macrobenthic community has revealed itself as a difficult task due to spatio-temporal heterogeneity. This study uses ordinary kriging and Poisson modeling to generate distribution maps of the subtidal benthic macrofauna in the Trapandé Bay (southeastern Brazil). Samples were taken in duplicate from 36 locations distributed along nine transects perpendicular to the main estuarine axis in October 2006 and March 2007. One-hundred and seventy taxa belonging to 12 phyla, were identified, with dominance of Annelida Polychaeta. Distribution maps were prepared to illustrate the total density, the number of species and the six most numerous taxa, as well as abiotic parameters. The general distribution pattern has revealed that the greatest number of species and the highest densities are at the estuary mouth, decreasing towards its inner areas. However the temporal and spatial changes observed at the estuary mouth have clearly shown the impact of environmental variations such as nutrients and freshwater input, attributed to increased rainfall in March. The increased flow in the Cananeia Sea, coming from the drainage basin, produces major changes in sediment and faunal composition. Ordinary kriging associated with Poisson modeling has proved to be a powerful and promising tool for modeling the macrofauna, despite the fact that it is not frequently used due to the scarcity of appropriate software.

The distribution of benthic organisms in unconsolidated sediments is determined both by abiotic factors, such as sedimentary and physicochemical characteristics, and by biological factors, such as predation, competition and bioturbation (DAUER 1993, LEVINTON 1995, SNELGROVE 1998, THRUSH & DAYTON 2002, MAIRE et al. 2010, FANJUL et al. 2011. The active interdependence among multiple potentially important factors causes the distribution of benthic associations to be frequently asymmetrical, non-linear and heterocedastic (ANDERSON 2008). Therefore, the distribution of benthic associations resulting from these complex interactions is very heterogeneous (SCHNEIDER 1994), which renders the recognition of spatial patterns difficult. In addition to the heterogeneity in species distribution, another complicating factor for modeling efforts is that species variability is not constant along environmental gradients, as it happens in estuaries, but is instead dependent upon the average abundances of species (ELLIOT & WHITFIELD 2011).
Several studies have shown that the spatial variation in the structure of benthic communities and the nature of marine sediments occur at various scales (e.g., THRUSH et al. 1997, BELL et al. 1999, ANDERSON 2008, CHAPMAN et al. 2010, RODRIGUES et al. 2011), contributing to a better understanding of the ecological relationships between physical processes and species distribution (ELLIOT & MCLUSKY 2002).Thus, this work describes the spatial distribution of subtidal benthic macrofaunal associations and the environmental factors in Trapandé Bay, a subtropical estuarine system in southeastern Brazil; compares the changes of these variables between two seasons and; evaluates the relationship between environmental variables along an estuarine gradient and the benthic macrofauna. We hypothesize that the response of the benthic macrofauna to environmental variables may be visually detected in spatial distribution maps, using appropriate geostatistical methods.

MATERIAL AND METHODS
The Trapandé Bay is located in the southern sector of the Lagoon-Estuarine Complex of Cananeia-Iguape, southern coast of São Paulo, Brazil. The estuary is geologically constructed by a series of bars, which form a complex system of channels between Comprida, Cananeia and Cardoso islands. The main islands are separated by channels ( Fig. 1) that connect to the ocean through the Cananeia Inlet, in the central sector of the estuary, the Icapara Inlet, in its northern sector, and the Ararapira Inlet, in its southern sector (BERNARDES & MIRANDA 2001 The main channels that connect with the Trapandé Bay are the Cubatão Sea in the Northwest, and the Cananeia Sea in the Northeast, which have average depths of 6 and 10 m, respectively. The opening to the Southwest is approximately 10 m deep. The regional drainage basin has variable seasonal discharges and may carry high concentrations of suspended sediments (MIRANDA et al. 2002). About 60% of the Ribeira River discharge flows by the internal channels of the Cananeia-Iguape estuary, carrying suspended sediments and nutrients to the Trapandé Bay through the Cubatão Sea and the Cananeia Sea (SAITO et al. 2006). The rainy season extends from December to April, and the dry season from May to November (SILVA 1984).
The tidal regimen is semi-diurnal, with inequalities during a diurnal cycle (MIYAO et al. 1986). Currents are more intense on the southern bank of the Trapandé Bay, near Cardoso Island, and less intense on the northern bank, near the Cananeia Island (TESSLER & SOUZA 1998).
Sampling was carried out in October 2006, during the dry period, and in March 2007, during the period of highest average annual rainfall. In each campaign, nine transects were sampled at the main axis, covering the entire extent of the bay, and four sampling points were established on each transect (Fig. 2). Transects were divided into inner (1, 2 and 3), intermediate (4, 5 and 6) and outer (7, 8 and 9). At each point, two Figures 1-2. A map illustrating the study area, Trapandé Bay and the nine sampled transects (labeled 1 to 9). Sampling was performed at four points in each transect. 1 2 samples were taken for benthic macrofauna and one sample for granulometric and chemical analysis using a Van Veen grab with an area of 0.065 m 2 . The faunal samples were fixed with 5% formaldehyde and washed on 0.5 mm mesh sieves. Animals were counted and identified to the lowest possible taxonomic level. The replicates were pooled and the data were converted into individuals per square meter (ind/m²). Granulometric measures were performed through the method described by SUGUIO (1973) and the data was analyzed using the package rysgran (GILBERT et al. 2012) within the R environment (R DEVELOPMENT CORE TEAM 2012). The sedimentary organic matter content was determined after 1 hour at 550ºC. Total nitrogen and phosphorus content were determined following GRASSHOFF et al. (1983).
Geostatistics was originally developed by MATHERON (1971) andJOURNEL &HUIJBREGTS (1978) in mining geology and is currently employed in a variety of scientific fields, such as soil sciences, meteorology, hydrology and ecology (AXIS- ARROYO & MATEU 2004). Although it is most frequently applied to marine fisheries (PETITGAS 2001, MAYNOU et al. 1998, SIMARD et al. 1992, RUFINO et al. 2004, it has also been used to assess the patterns of spatial distribution of the benthic macrofauna (BERGSTROM et al. 2002, COLE et al. 2001. The geostatistical Gaussian models do not properly estimate the response patterns of benthic species counting due to the existence of a correlation between a species' average density at a given point and its variance (ANDERSON 2008). Whenever such a correlation exists, Poisson models are more adequate and generate results that are far more realistic. Here, geostatistical maps were generated by ordinary kriging, and adjusted to the Poisson model for densities of the most representative species (in terms of species abundance and occurrence), the total number of species, and the total species abundance. Maps were generated using the R package geoRglm (CHRISTENSEN & RIBEIRO JR 2002).

RESULTS
In transects 3, 4 and 5, sediments were predominantly silty. Mean grain size was coarser towards the inner bay and the Cananeia Inlet, with higher values around 2 ⌽ (fine sand) (Fig. 3).
Total nitrogen (TN) ranged between 0 and 3.5 mg/g, with higher values in the inner and intermediate regions, decreasing towards the Cananeia Inlet in both campaigns (Figs 4 and 5). Total phosphorus contents (TP) were higher in the October campaign with maximum values of 3.0 mg/g in the outer region and lowest values in transects 4 and 5. In March, TP did not exceed 0.7 mg/g near the Cananeia Inlet (Figs 6 and 7).
Organic matter concentrations (OM) were higher in the northern sector of the bay in October, especially in transect 4, with values up to 20%. The lowest OM values occurred in the southern sector of the outer transects. In March, the values of OM were homogeneous along transects 1 to 8, despite the high value observed in the southern sector of transect 3. The lowest concentrations were observed in transect 9 (Figs 8 and 9).
A total of 5,531 organisms were identified, distributed among 12 phyla and 170 taxa. Annelids were numerically dominant, totaling 87.8% of the overall abundance, together with molluscs and crustaceans. The most representative taxa in terms of abundance and occurrence were Aricidea sp., Capitellidae, Magelonidae, Sipuncula, Scoletoma tetraura (Schmarda, 1861) and Ophiuroidea, which together totaled 2,154 organisms and comprised 39% of the total abundance (Table I).
In October, species numbers gradually increased towards the mouth of the bay and were greatest in transect 9. In March, species numbers were greatest in transects 5, 6 and 7, where approximately 40 taxa per m 2 were observed (Figs 10 and 11).
Littoridina australis (Orbigny, 1835) was the most abundant taxon due to its extremely high density at a single sampled point (691 individuals in one sample alone in March). However, aside from this extreme value, the abundances of L. australis were low in all the other samples in both campaigns. Likewise, Bulla striata (Bruguière, 1792) showed high densities in specific occurrences in transect 2. Both taxa were removed from the analysis.
The total density maps indicated higher densities in the outer transects in October, whereas in March the highest densities were concentrated in the intermediate transects and decreased in transect 9 (Figs 12 and 13). In March, transect 5 had the highest total density on its northern sector, where the to-   (Figs 14 and 15). There was a gradual increase in its density towards the northern sectors of the outer transects, where the highest densities observed ranged between 400 and 600 ind/m 2 in October. However, in March, the average densities in the outer areas were reduced to 70 ind/m².
Members of Capitellidae were more abundant in the outer areas, particularly in the sampling points at the south shore, and displayed variance between the two campaigns only in transect 9, which presented lower densities in March than in October (Figs 16 and 17). Members of the Magelonidae were first found in high densities in transect 4, and tended to increase in density in the outermost transects. In March, the highest densities of Magelonidae were restricted to the intermediate transects (4, 5 and 6). A maximum of 250 ind/m 2 was observed in the northern sectors of transects 5 and 6, and decreasing densities were noted at the mouth of the bay (Figs 18  and 19).
In October, Ophiuroidea increased in density towards the mouth of the bay, with higher density values on the southern banks of transects 8 and 9, where up to 200 ind/m 2 were recorded. In March, the highest densities of Ophiuroidea were recorded in transects 5, 6 and 7, with a maximum of 50 ind/m 2 on the southern bank sampling points (Figs 20 and 21). The highest densities of S. tetraura were recorded in October, mainly concentrated on the southern sector of transect 5 and the northern sector of transect 8, with about 200 ind/m 2 (Fig. 22). In March the highest densities of S. tetraura occurred in transects 5, 6 and 7. In these locations, average densities varied between 80 and 100 ind/m 2 (Fig. 23). Inner transects showed the lowest densities of S. tetraura in both sampling campaigns.
Sipunculans were distributed in a similar pattern in both campaigns, with higher density values consistently observed in the intermediate transects. In October, their highest densities occurred in transects 4 and 5, shifting to transects 5 and 6 in March (Figs 24 and 25).

DISCUSSION
The Trapandé Bay exhibits a well-defined environmental gradient with population densities and number of taxa decreasing towards the inner region. The intermediated sector, which corresponds to transects 3 to 7, contained a higher number and density of species than the inner region, particularly on the north bank of the estuary, which is shallow. The outer region presented the highest faunal differences between campaigns, presumably due to the powerful hydrodynamic effects generated by interactions with the Cananeia Sea and the Icapará Inlet. This region has strong, dynamic currents near its bottom (>1.0 m/s) due to the influence of the mouth of the Cananeia Sea and the Icapará Inlet (MIYAO et al. 1986), as well as the constriction of the lagoon channels (KUTNER 1962, TESSLER & SOUZA 1998. The friction stress caused by the bottom currents tends to resuspend fine sediments, nutrients and organic matter. The Cananeia Sea influences the outer region of the bay due to the freshwater and nutrient input from the Ribeira do Iguape River through the Valo Grande Channel. Although located 40 km from Trapandé Bay, the Valo Grande Channel influences nutrient cycling in the southern sector of the Cananeia-Iguape system. MIYAO et al. (1986) observed higher concentrations of phosphates in the mouth region of the Trapandé Bay that were caused by the resuspension of sediments at the tide entrance. JORCIN (2000) found higher concentrations of total phosphorus and nitrogen in the Cananeia region during the summer, in areas coinciding with the region of transects 8 and 9 that were associated with a seasonal increase in the levels of organic matter.
The circulation channel of Trapandé Bay is located next to the Cardoso Island and there is a transition within the bay from sandy sediments, associated with the channel, to sediments with a high mud content (silt and clay), which are found on the northern sector of the region, near the Cananeia Island (SAITO et al. 2006). The northern margin of this area is less deep and suffers less hydrodynamic influence.
The distribution pattern of benthic organisms observed on estuaries typically fit the ecocline model. The community progressively changes following the gradual variation in at least one major environmental variable. In estuaries, this variable usually is salinity, but hydrodynamics, organic matter and mud content play important roles. Furthermore, estuaries differ from previously defined terrestrial ecoclines because they present two overlapping gradients in the major stressor: from river to mid-estuary for freshwater species and from sea to mid-estuary for marine species (ATTRILL & RUNDLE 2002).
Sediments with higher organic matter contents are likely to present higher abundances of deposit-feeder polychaetes and further opportunists or early colonizing species (PEARSON & ROSENBERG 1978). Capitellid polychaetes are a good example of opportunistic species, characteristically found in areas with heavy input of organic matter (SCOTT et al. 1987, RIZZO & AMARAL 2001. Individuals of Paraonidae family are considered non-selective surface diggers (FAUCHALD & JUMARS 1979). In the October campaign, Aricidea sp. was predominantly found at the shallow waters near the Cananeia Sea, associated with higher values of TP and OM. However, in March, the densities of Aricidea sp. decreased, most likely due to the increased rainfall and hydrodynamic flows at the estuary's bottom, which leaded to a decrease in TP and OM contents.
Magelonid polychaetes, a pooling of Magelona papillicornis (Müller, 1858), Magelona posterelongata (Bolivar & Lana, 1986) and Magelona variolamellata (Bolivar & Lana, 1986), were present in the intermediate and outer zone of the Trapandé Bay, not well correlated with the observed peaks of organic matter content. Higher densities of magelonids were usually found in sediments with about 5-10% of organic matter. Magelonids are mainly non-selective surface deposit feeders, but some species may alternate to suspension feeding (FAUCHALD & JUMARS 1979). This alteration in feeding strategies may make it possible for magelonids to live both in muddy and sandy sediments (ROUSE & PLEIJEL 2001).
The distribution pattern of the omnivorous S. tetraura was intrinsically linked to the total density and number of species. This specie may be more influenced by interspecific relationships than by environmental variables. Ophiuroids are typically marine and are thus restricted to the outer estuarine area in the dry season (October). In March, most likely due to the freshwater inflow that originates from the North through the Cananeia Sea, their occurrence was restricted to the southern bank of the outer area.
The general benthic community pattern observed in the Trapandé Bay, in which highest densities and species numbers were recorded at the estuary mouth and decreased towards the inner areas, is recurrent in other estuaries from southern and southeastern Brazil, as shown in Cananeia (TOMMASI 1970), Paranaguá Bay (LANA 1986) and Patos Lagoon (BEMVENUTI & NETTO 1998). The temporal and spatial variation in faunal parameters clearly demonstrated the impact of environmental variations in either chemical composition or particle size on the distribution of the macrofauna. The marked temporal variations may be mainly attributed to increased rainfall in March, when the increased flow in the Cananeia Sea coming from the Ribeira do Iguape River affects sediment properties.
Interpolation procedures represent a gain in information over unsampled areas, but with the restriction that the results contain some degree of uncertainty (JEROSCH et al. 2007).
Several statistical/mathematical methodologies exist on how to derive spatially continuous distributions from available localized point data (LI et al. 2011 for a review), for instance decision tree models (PESCH et al. 2011), machine learning methods (LI & HEAP 2008) and combined methods that use ordinary kriging with Gaussian linear models (STELZENMÜLLER et al. 2010) or non-Gaussian generalized linear models (GML) (GOTWAY & STROUP 1997, MAXWELL et al. 2009). GLM is more adequate than linear models to analyze spatial structure of the benthic macrofauna because it can fit Poisson data, which is much more suitable for counting data like abundance or richness. Most available software for analyzing spatial data, however, works only with Gaussian methods and that is the reason why the method is rarely used in modeling benthic macrofauna. ANDERSON, M.J. 2008. Animal-sediment relationships re-visited: characterizing species' distributions along an environmental gradient using canonical analysis and quantile regression