The ecological interactions and mechanisms that mediate coexistence among species are often poorly understood, despite its relevance to conservation management. A case in point is the Indo-Pacific bottlenose dolphin (Tursiops aduncus) and Indian Ocean humpback dolphin (Sousa plumbea) that occur sympatrically along the south-east coast of South Africa. To improve on our understanding of their coexistence, data on spatial distribution and habitat use were collected using boat-based surveys during 2014-15 along 145 km between Goukamma and Tsitsikamma MPAs. During 235 hours of survey effort, 55 and 42 encounters with T. aduncus and S. plumbea were recorded respectively. We investigated differences in their space use and habitat preferences. There was a strong overlap of the core areas used by both species, located mainly along Goukamma MPA and the eastern section of Plettenberg Bay. For both species, foraging was observed in mornings and travelling in afternoons. T. aduncus showed a preference for rocky areas in the afternoon while S. plumbea preferred estuarine and sandy habitats. We demonstrated that coexistence among dolphins with overlapping distributional patterns can be driven by fine-scale spatial segregation in feeding activities. Spatial conservation management measures such as the protection of sandy and estuarine habitats are recommended.

Study area and data collection 

Boat-based surveys for dolphins were conducted along the south coast of South Africa within the Agulhas Bioregion in 2014–15. When a dolphin or groups of dolphins were encountered, they were identified and photographed, the position of the sighting was recorded and the behaviour of the dolphins were categorised. The surveys occurred between the western boundary of the Goukamma MPA and the eastern boundary of the Tsitsikamma MPA (Figure 1). The area is strongly influenced by the warm Agulhas Current [22] and wind-driven upwelling results in high levels of primary productivity and associated high prey biomass for predators [22]. 

The surveys were designed as a transect line of 145 km distance, running parallel to the coast (approximately 100 m from the coastline), corresponding with the known coastal preferences of both dolphin species [16,23]. A full survey was defined as transect completed within a calendar month between April 2014 and October 2015 (n: 13, Table 1). This ensured that the survey effort was similar throughout the study site. Survey effort was calculated in hours and encounter rate was calculated as the number of encounters per hour searched. Survey effort was discontinued when the Beaufort sea-state was above 3, when the boat was in transit and during encounters. Surveys were performed at a speed of eight knots with at least two or three experienced observers.

A group was defined as two or more animals within a 100-m radius of each other, with coordinated activities [24]. In addition to the GPS coordinates of the encounters, sea surface temperature (SST) depth and time of day were recorded. Five behavioural categories for T. aduncus were defined according to [25]: travelling, foraging, socializing, milling or resting. Four behavioural categories were defined for S. plumbea according to [26]: travelling, foraging, socializing or resting. See [27] for further information on study area, survey design and the full data collection procedures.

Data analysis

Spatial and behavioural analyses

The spatial distribution of T. aduncus and S. plumbea was assessed by means of a kernel density estimator (KDE) analysis [28]. This probabilistic technique provides estimates of the utilisation distribution (UD), a probability density function that describes the relative use of space by an animal within a defined area, based on a sample of animal locations [28]. The geographical coordinates of each group of animals at the time when they were first encountered provided the spatial data for the analysis. A bivariate normal kernel UD was used to determine the home ranges (95% UD) and core areas (50% UD) for each species. The probability of occurrence was calculated using smoothing parameters for the kernel using the least squares cross-validation (LSCV) method [29]. The LSCV method examines various smoothing parameters and selects the bandwidth that gives the lowest mean integrated squared error for the density estimate [30]. KDE analysis was performed using program R 1.0.143 [31] with ‘adehabitatHR’ package, v. 0.4.14 [32].

To assess spatial segregation between T. aduncus and S. plumbea, we calculated a UD overlap index (UDOI) to quantify the degree of overlap between the home ranges of the two species [33]. A randomization approach was then applied whereby the null hypothesis that there is no difference in spatial distribution between the two species, was tested [34]. For the null hypothesis to be upheld, there should be no significant differences between the degree of overlap calculated for the original home ranges and the degree of overlap calculated for home ranges for which species were randomly assigned. Home ranges were generated from 1000 randomly generated datasets in which either species was randomly allocated to each location. P-values were calculated from the proportion of random overlaps that were smaller than the observed overlap. Behaviour and habitat preference were assessed further in relation to the core areas (i.e. inside and outside 50% UD). Pearson’s chi-square homogeneity tests [35] were used to investigate if behavioural states varied spatially in relation to the core areas and time of the day [20].

Habitat preference and temporal segregation

To assess habitat preference (when first encountered) and potential temporal segregation of both species, 2 km² grid cells along the coast were created using QGIS 2.18.4 [36], resulting in 73 grid cells in total. Each cell was characterised according to the benthic substrate, presence or absence of reefs, mosaic (a combination of reef and seafloor) or estuaries, and in relation to the boundaries of MPAs (inside or outside). Within the study area, three broad types of benthic substrate exist, namely rocky, sandy and mixed coast (i.e. rock and sand). The benthic substrate types were obtained from the benthic and coastal habitat map of the National Biodiversity Assessment [37].

Grid cells were also characterized for each survey according to broad time of day (AM/PM), season, sea surface temperature (SST) and the occurrence of T. aduncus and S. plumbea (calculated as the presence-absence of dolphins in the 2 km² grid cell during each survey). Seasons were defined as (1) summer = December-February, (2) autumn = March-May, (3) winter = June-August, and (4) spring = September-November.

Occurrence (presence-absence) data of T. aduncus and S. plumbea as a function of the spatial and temporal characteristics of the grid cells and surveys were first visually assessed for normality and then analysed separately, using a generalised linear mixed-effects models (GLMM) with binomial distribution (link function: logit). Analyses were conducted using the ‘lme4’ package v. 1.1-12 [38] in the freeware R 1.0.143 [31]. For each model, all predictor variables were included in the analysis: relation to MPA boundaries, benthic substrate, presence or absence of reefs, mosaic and estuaries, and also an interaction term between benthic habitat and broad time of day. The 2 km² grid cell was included as a random intercept to account for repeated measures and models were fitted by restricted maximum likelihood estimation (REML).

Collinearity between all covariates was tested. The Variance Inflation Factor (VIF) scores were calculated for each predictor variable using the ‘vifcor’ function of the ‘usdm’ R package v. 1.1-15 [39]. Only uncorrelated covariates (VIF < 3) were included to avoid misinterpretation of the model [40]. Model selection was based on the Akaike information criterion scores corrected for small sample sizes (AICc), whereby all realistic permutations of predictor covariates were fitted to separate models under maximum likelihood (ML) estimation using the ‘dredge’ function in the ‘MuMIn’ package v. 1.15.6 [41]. All models with a ΔAIC < 2 were considered as candidate subsets of the most parsimonious models. Candidate subsets were then used to generate model-averaged coefficients [42]. As a measure of the goodness-of-fit, pseudo R2 values were calculated for the model to explain the proportion of variance accounted for by the fixed and random factors [43].

References

R programme and Microsoft Excel