Integrated risk analysis for rare marine species impacted by fishing: sustainability assessment and population trend modelling

Abstract

Zhou, S., Milton, D. A., and Fry, G. C. 2012. Integrated risk analysis for rare marine species impacted by fishing: sustainability assessment and population trend modelling. – ICES Journal of Marine Science, 69: 271–280.

Risk assessment and conservation management of rare species are challenging due to a lack of data. We developed an integrated risk assessment approach to assess human impact on population sustainability of rare species. The approach involved two components: a quantitative sustainability assessment coupled with modelling trends in relative abundance. Both components took nil catches into account through zero-inflated statistical distributions that simultaneously modelled the zero and non-zero catches separately in submodels. The sustainability assessment used detection–non-detection data for population estimation and linked sustainability to easily collected life-history traits. This component provides an assessment of population sustainability at one point in time. The trend modelling applied zero-inflated negative binomial models to temporal trends in density and dispersion of species. It provided a complement to the static sustainability assessment. We applied this integrated approach to assess the risk to 14 species of rare, protected sea snakes incidentally caught in the Australian Northern Prawn Fishery. This approach can be applicable for risk assessment of many species with limited abundance data, a large number of absences and some presence–absence information only.

Introduction

The majority of species in most ecological communities tend to be rare (Cunningham and Lindenmayer, 2005). Risk assessment and conservation management of rare or threatened species are challenging for both terrestrial and aquatic ecosystems. Rare species typically lack historical data on their distribution and abundance. Records of detection may come from surveys designed for other species and little information may exist about their life history. These difficulties constrain detailed risk assessment and so resource managers often choose to introduce measures to conserve rare species in the belief that they are the species most at risk of extinction (Gaston, 1994). However, there is a trade-off between reducing the risk to populations of a rare or threatened species and reducing the social and economic costs (Rice, 2011). For the sustainable use of economically important natural resources, it is important to know whether complete conservation of rare incidentally caught species is needed (Fletcher et al., 2010).

In terrestrial systems, the most well-known process to assess the species risk involves qualitative assessments by experts (IUCN, 2001). Similarly, most risk assessment methods in marine systems assess the relative risk of population collapse (Milton, 2001; Astles et al., 2009; Heino, 2011; Hobday et al., 2011). To improve the assessment of the absolute risk to species, a sustainability assessment for fishing effect (SAFE) method was developed and applied to fisheries with diverse bycatch (Zhou and Griffiths, 2008; Zhou et al., 2009). However, even quantitative approaches such as these provide only a static assessment and do not account for historical trends (Kronen et al., 2010).

Population trend analysis to detect the potential decline in adult abundance over time is the most widely used method of assessing extinction risk. This approach is used by many organizations including CITES and IUCN. However, the population trend alone, even declining, may not indicate a risk of extinction or overexploitation (Dulvy et al., 2004).

There is a need for an integration of trend modelling with sustainability assessment to better quantify ecological risk, especially for rare and/or threatened species. We develop and apply an integrated approach to assess the risk to sea snakes caught by a prawn (shrimp) trawl fishery in Australian Northern Prawn Fishery (NPF). Sea snakes are relatively rare, air-breathing reptiles found in coastal waters in tropical Australia and the Indo-west Pacific. Sea snakes are caught by the inshore prawn trawl fisheries in southern and southeast Asia and northern Australia (Steubing and Voris, 1990; Wassenberg et al., 1994). They may be at greater risk of extinction from trawling than many fish species as they are long-lived and have a low reproductive rate (Heatwole, 1997; Fry et al., 2001). Several species have recently been identified by the IUCN as threatened (C.T. Elfes, pers. comm.). Under Australian legislation, sea snakes are the protected species. All Australian fishing industries are required to demonstrate the sustainability of their activities when they directly or indirectly interact with protected species.

Methods

Data sources

Data were collected from a number of sources:

Scientific surveys: Scientific surveys within the NPF management area were designed to study fish or prawns. All sea snakes incidentally caught during these surveys were identified, measured (total and snout-vent length), and released. Survey and fishery data were collected from historical 6 × 6-nautical mile grids inside and outside the major prawn trawling grounds. We applied pre-existing strata based on five bioregions identified from species assemblages and experts' opinion (Figure 1). About 100 surveys or research studies had been undertaken in 1114 of 6963 grids across the entire NPF area from 1975 to 2007. Some surveys only contained the presence–absence data for sea snakes, whereas the others recorded the number of sea snakes (count data) in the sample.

Figure 1.

Sampling locations (cross) and the commercially fished area (black square) during 2004–2006 in the five bioregions of the NPF area.

Logbook data: Fishery logbooks provided the location and intensity of trawling data. We used the data from 2004 to 2006 to estimate the distribution and intensity of fishing. Grids were classified as fished if they had at least one boat-day of trawling during these 3 years.

Commercial fisheries with scientific observers or crew-member observers: During 2003–2005, scientific observers and crew-member observers collected a large quantity of bycatch data from commercial prawn trawling in the NPF. Images were taken of each sea snake caught and labelled with a scale bar. Species were identified from the images and their size estimated with digitizing software (Image J).

Museum records: The available records of sea snakes were obtained from specimens held in the Queensland, Northern Territory, Western Australia, and Australian Museums.

Life-history and other parameters: Life-history parameters were obtained from previous studies (Fry et al., 2001; Milton, 2001; Wassenberg et al., 2001; Ward, 2001; Milton et al., 2009), including maximum age, age at maturity, and within-trawl and post-trawl survival rates. Additional life-history parameters (e.g. maturation age, maximum age, and length infinity) were also derived from specimens collected in the scientific surveys and observer programmes.

Sustainability assessment for fishing effects

Estimating incidental fishing mortality rate

For each species, annual incidental fishing mortality rate u (the fraction of population killed by fishing) was estimated from the relative abundance of each species within trawled areas relative to the entire NPF management area, the catchability when the trawl passed over seabed where sea snakes live, the mortality rate after capture (including both within-trawl and post-trawl mortality), and fishing effort:

(1)

where N_R,F is the estimated abundance of a species in bioregion R and the commercially fished (F = 1) or unfished (F = 0) area, A_{R,F = 1} the commercially fished area in bioregion R, a the mean swept-area per boat-day, E_R the annual fishing effort (boat-day) in region R, q the catchability defined as the probability of a sea snake on the trawl path being caught, and s the post-capture survival rate. Equation (1) assumes that there would be no local depletion effects from repeat trawls at the same location and that populations rapidly mix between trawled and untrawled areas. The fishing mortality rate is likely to be overestimated under this assumption.

We initially attempted to estimate the sea snake relative abundance of N_R,F with biophysical predictive models that linked density to environment variables (Guisan and Zimmermann, 2000; Pitcher et al., 2007). Preliminary analyses indicated that this approach resulted in unacceptably large estimates of uncertainty for several species. Moreover, the model parameters could not be estimated for other species due to insufficient data. Hence, we opted to use an alternative and more flexible method to estimate N_R,F. This method required less data and used simple detection–non-detection records that were more widely available (Zhou and Griffiths, 2007). The principal concept of this approach is the same as the zero-inflated binomial model (Hall, 2000), which has been widely used for occupancy modelling (MacKenzie et al., 2006). We briefly describe the method as follows.

We assumed that individual sea snakes were randomly distributed within each of the ten strata (five bioregions each divided into commercially fished and unfished areas), which is considered appropriate as the majority of samples only caught one sea snake. Let g be the area of each grid, then 1 − g/A_R,F is the probability that one particular snake is not present in the grid in stratum (R, F). The probability of grid i being occupied (has at least one sea snake of a particular species) is

(2)

Let m_R,F,i denote the number of survey conducted in grid i, n_R,F,i be the number of times a species is detected (m_R,F,i ≥ n_R,F,i), and D_R,F,i be the probability of detecting one or more individuals of a particular species in stratum (R, F) and grid i. Further, we assumed that grids in the same stratum (bioregion and fished or unfished area) had the same probability of occupancy p_R,F (homogeneous distribution within a stratified stratum) and grids in the same bioregion have the same detectability D_R. The assumption of constant detectability may be difficult to hold for each grid cell. However, simulations show that this assumption has little impact on the total abundance estimation (Zhou and Griffiths, 2007). We can then obtain the joint likelihood function for all grids over all strata as:

(3)

where the underline denotes a vector, f_B the binomial density function, and G₁ and G₀ the grids where a species of sea snakes has or has not been detected, respectively. To ensure that D_R is between 0 and 1, we re-parametrized it as and linked θ to covariates of gear type and the area swept.

In Equation (1), swept-area a by one vessel was estimated to be 2.40 km²boat-day⁻¹ (±0.28 s.d.), and the overall post-capture mortality rate (1 − s) from trawling was 48.5 ± 1.4% (Wassenberg et al., 2001). The catchability q is difficult to obtain without experimental trials specifically designed to estimate this parameter. We used a conservative estimate of 1 based on the results of a “relative catchability” study (Pitcher et al., 2002). The relative catchability was estimated from individual catch rates of multiple gear types used at the same site. The gear that had the highest catch was assigned q = 1. In general, this treatment is conservative because it tends to underestimate abundance.

Sea snake life-history parameters

In all, 633 sea snakes of 12 species were collected during scientific surveys and the scientific and crew-member observer programmes. Species and sex of these snakes were identified, their snout-vent length measured in mm, their reproductive biology examined, and their age estimated following the methods of Fry et al. (2001) and Ward (2001), respectively. The sea snakes' average age at maturity and the maximum reproductive age were determined from annuli counts in the vertebrae and gonad staging. These growth rings in the vertebrae have previously been validated as annual for the two most abundant species in our study area (Ward, 2001). As reptiles have similar life-history traits to fish (Charnov et al., 1993), the length–age frequency data for each species were fitted by the von Bertalanffy growth model (Quinn and Deriso, 1999):

(4)

where L_t is the length at age t, L_∞ the asymptotic length, k the von Bertalanffy growth coefficient, and t₀ the age when an individual would have been of length zero. Estimates of k and L_∞ were then used by some methods to derive the estimates of natural mortality rates (see below).

Sustainability reference points

Fisheries management usually uses two types of reference points: one based on the biomass and the other based on the mortality rate (Quinn and Deriso, 1999). We did not have sufficient data to estimate the biomass-based reference points such as B_msy so we focused on the mortality-based reference points. Theoretically, if the fishing mortality rate is constant each year, the abundance will stabilize to the level corresponding to that fishing mortality. We used the following two biological reference points: u_msm is the fishing mortality rates corresponding to the maximum sustainable mortality (MSM) rate caused by fishing, and u_crash is the minimum unsustainable fishing mortality rate that, in theory, may eventually lead to population extinction if it continues. We refer to this section as “sustainability assessment” because we used these sustainability reference points to gauge the impact on species. MSM is similar to maximum sustainable yield (MSY) for target fish species. The two corresponding reference points take recruitment compensatory processes into accounts and is currently the standard for single-species stock assessment. Estimating biological reference points can be difficult. Reliable estimation requires parameters derived from quantitative stock assessments and often requires time-series data and considerable biological information. An alternative approach is to identify a relationship between reference points and life-history traits, because the two are fundamentally related to each other (Rochet, 2000). Many studies have tried to establish such a relationship. Among a range of life-history parameters, the natural mortality rate is the most widely used one for fisheries management, including bony fish, sharks, and invertebrate fisheries (e.g. see “Status of U.S. Fisheries” at http://www.nmfs.noaa.gov/sfa/statusoffisheries/appendix_3_2008_sdc_all.pdf). Sea snakes have similar life-history characteristics to sharks and produce few live young annually (Charnov et al., 1993; Fry et al., 2001). We adopted the common practice in fishery management of using the estimated natural mortality rate as a surrogate for the optimal fishing mortality rate. However, we chose a more conservative relationship for sea snakes by adopting u_msm = 1 − exp(−0.5 M) and u_crash = 1 − exp(−M). This is half the value commonly used for fish species. Recent analysis suggests this relationship is close to that of chondrichthyans (S. Zhou et al., unpublished data). We used four methods to estimate natural mortality rate M: The average value of u_msm and u_crash from these four methods was used as reference points and the minimum and maximum values as estimates of the range of uncertainty. We obtained the variances of N_R,F from the statistical model based on the detection–non-detection data (Zhou and Griffiths, 2007). The variances of a and s were obtained from previous studies (Wassenberg et al., 2001). The variance of fishing mortality rate u was estimated by the delta method. For the sustainability reference points, we lacked information on the variance of the input parameters themselves. So instead of using confidence intervals to describe uncertainty, we used the range to indicate uncertainty associated with the two reference points.

• ⁠, where T is the average annual water temperature (Pauly, 1980);

• ln(M) = 1.44 − 0.982 ln(t_m), where t_m is the maximum reproductive age (Hewitt and Hoenig, 2005);

• M = 1.50 k (Jensen, 1997); and

• M = 1.65/t_mat, where t_mat is age at maturity (Jensen, 1997).

Modelling temporal catch trend

The sustainability assessment provided a snapshot of effects of fishing on species but does not take into account temporal changes in populations. The second component of the integrated risk assessment approach modelled temporal trends in populations. We used the data from scientific surveys that contained actual counts of sea snakes in each sample. These count data were characterized by a large number of zero-valued samples and high interannual variation. We investigated several approaches [generalized linear models, Hurdle models, negative binomial models, zero-inflated Poisson regression, and zero-inflated negative binomial (ZINB) models; Cameron and Trivedi, 1998] for modelling catch as a function of covariates, including region, commercially fished or unfished area, depth, size of sampling area (gear swept-area), year, and season. We found that the ZINB regression model (Hall, 2000; Minami et al., 2007) generally fitted the data best based on Akaike's information criterion and Vuong tests (Vuong, 1989).

ZINB is a two-component mixture model combining a point mass at zero with a negative binomial count distribution. There are two sources of zeros: from the point mass and from the count component. The probability function of ZINB is expressed as:

(5)

where c_i is the number of a particular sea snake species caught in sample i, p the probability of non-zero catch that can be considered as the probability of occupancy by that species of sea snakes, f_NB the negative binomial probability density function, μ the mean number of sea snakes caught per sample, and κ the shape parameter that can be interpreted as the degree of aggregation. One parameterization of this probability density function is

(6)

where Γ(θ) is the gamma function. p in Equation (5) directly links to the probability of occupancy, one of the IUCN red list criteria, and is modelled as 1 − p by a logit link:

(7)

where X are the covariates and β are the parameters. The mean catch μ is modelled with a log-link:

(8)

where X can be the same covariates as in Equation (7) and γ the parameter (Hall, 2000; Minami et al., 2007). The integrated mean catch for all samples is the product of the two components:

(9)

We used a stepwise procedure by dropping the variables that were not significant and investigated the goodness-of-fit of the ZINB models with and without predictors (the null model with intercept-only). This was done with likelihood ratio tests between the ZINB models.

To detect the temporal catch trend, we estimated the effect of year on p in Equation (7) and the mean catch μ of count component in the log-linear regression model [Equation (8)]. We held the covariates of continuous variables (depth and swept-area) to their mean values for all years, held the categorical variables (season, area of commercially fished or not, region) to their averaged relative frequency for each level across all years, and only allowed the year parameter to vary. The combined year effect on the integrated mean catch E[c] is the product of the two parts of the ZINB model as in Equation (9).

Results

Distribution of sea snake populations from detection–non-detection data

There were 14 sea snake species caught in 5481 of 18 020 trawl samples (Table 1). Hydrophis elegans was the most common species and was caught in 1809 samples. The second-most common species was Lapemis hardwickii, which was caught in 1097 samples. One species, Hydrophis caerulescens was only caught in commercially unfished areas.

Table 1.

Number of times (detections) when one or more sea snakes were caught (in surveys or commercial trawls) and the detection rate (detections/samples).

Species	Detections		Detection rate
	Unfished	Fished	Unfished	Fished
Acalyptophis peronii	56	74	0.007	0.007
Aipysurus duboisii	14	17	0.002	0.002
Aipysurus eydouxii	107	188	0.013	0.019
Aipysurus laevis	41	220	0.005	0.022
Astrotia stokesii	109	271	0.014	0.027
Disteira kingie	57	30	0.007	0.003
Disteira major	128	525	0.016	0.053
Enhydrina schistose	105	2	0.013	0
Hydrophis caerulescens	12	0	0.001	0
Hydrophis elegans	599	1 210	0.075	0.121
Hydrophis mcdowelli	29	22	0.004	0.002
Hydrophis ornatus	122	301	0.015	0.030
Hydrophis pacificus	18	127	0.002	0.013
Lapemis hardwickii	764	333	0.095	0.033

In recent years, prawn fishing occurred in a relatively small proportion of the grids in the NPF area (∼6%). However, 10 of the 14 species had higher densities in these commercially fished areas than unfished areas. As a result, the proportion of abundance in the fished area was typically higher than the proportion of the fished area (i.e. 6%). For example, Hydrophis pacificus had 25.7% of its population in fished areas, the highest of all species, followed by Aipysurus laevis with 23.0% in fished areas (Figure 2). Three species appeared to have a smaller proportion of their population in fished areas than unfished areas: H. caerulescens (0%), Enhydrina schistosa (0.8%), and L. hardwickii (2.3%).

Figure 2.

The estimated proportion of sea snake abundances in the commercially fished area (±95% CI). Hydrophis caerulescens was only caught in commercially unfished areas.

Relevant life-history parameters

Life-history data were incomplete for some species due to the less number of animals caught (Table 2). Mean length-at-maturity ranged from 472 to 1481 mm and mean clutch size from 3.0 to 12.3. Aipysurus duboisii had the shortest lifespan (t_m = 8.7 years) while Astrotia stokesii had the longest (t_m = 22 years). The average ages at maturity t_mat differed less among species, varying between 2.2 and 4.6 years. Generally, species with a large body size matured at an older age and had a longer lifespan than those with a smaller body size. However, body size appeared to have little effect on clutch size.

Table 2.

Relevant life-history parameters of sea snakes caught in the NPF area.

				Clutch size		L∞		k		M
Species	Lmat (mm)	tm	tmat	Mean	s.e.	Mean	s.e.	Mean	s.e.	Mean	Min	Max
Acalyptophis peronii	716	9.2	2.2	4.5	1.7	1 211	92	0.421	0.103	0.554	0.356	0.75
Aipysurus duboisii	910	8.7	3.1	4.5	1	1 096	62	0.47	0.091	0.534	0.393	0.705
Aipysurus eydouxii	472	12.6	2.2	3.6	0.3	813	45	0.631	0.18	0.641	0.351	0.946
Aipysurus laevis	940	16.2	4.7	6.5	1.8	1 427	201	0.178	0.06	0.271	0.193	0.351
Astrotia stokesii	817	22	4.3	9.9	1.7	1 228	44	0.294	0.045	0.327	0.203	0.441
Disteira kingii	1 000	–	–	4.9	0.6	1 165	259	0.446	0.401	0.512	0.316	0.688
Disteira major	710	14	2.4	4.9	0.3	1 048	23	0.537	0.081	0.561	0.316	0.806
Hydrophis elegans	1 183	16.2	3.3	12.3	1.3	2 038	91	0.25	0.037	0.342	0.219	0.5
Hydrophis inornatus	920	–	–	3	0	–	–	–	–	–	–	–
Hydrophis mcdowelli	635	–	–	3.7	0.9	–	–	–	–	–	–	–
Hydrophis ornatus	800	14.8	3.1	6	0.5	1 134	34	0.578	0.111	0.536	0.299	0.867
Hydrophis pacificus	1 481	13.8	4.6	–	–	1 801	92	0.383	0.07	0.388	0.299	0.574
Lapemis hardwickii	718	18.7	3	4.3	0.2	1 020	43	0.423	0.074	0.449	0.238	0.634

L_mat, length-at-maturity; t_m, maximum reproductive age; t_mat, age at maturity; L_∞, asymptotic length (mm); k, growth parameter; and M, natural mortality.

L_∞ estimates for 11 species ranged from 813 mm for Aipysurus eydouxii to 2,038 mm for H. elegans (Table 2). Aipysurus eydouxii also had the highest growth coefficient (k = 0.63), whereas H. elegans had the lowest (k = 0.25). The mean values of M for each species ranged between 0.27 and 0.64 year⁻¹. However, the uncertainty associated with M was large as estimates varied widely among different methods. The maximum estimate was typically about twice that of the minimum value.

Fishing mortality rate

We estimated an annual fishing mortality rate (u) between 0.0 and 0.055 for the 14 species (Figure 3). Hydrophis pacificus had the highest fishing mortality rate (u = 0.055), followed by A. laevis (u = 0.040). Although the estimated fishing impacts were low, uncertainties associated with u were relatively large. For example, the coefficient of variation (CV) of u ranged from 24% for H. elegans to 67% for E. schistosa. The mean CV for all the 14 species was 40%.

Figure 3.

Comparison between estimated fishing mortality rate u and mortality rate corresponding to maximum sustainable fishing mortality u_msm and minimum unsustainable fishing mortality u_crash. The error bars are ±95% CIs for the estimated fishing mortality rates and range (min–max) for the reference points.

Sustainability of sea snake populations

As M was only estimated for 11 species, the estimates of M for the remaining three species were based on those of closely related and of similar size species as surrogates: Disteira major for E. schistosa and the average values of H. elegans, H. ornatus, and H. pacificus for H. caerulescens and H. mcdowelli. This approximation tended to be more conservative as H. caerulescens and H. mcdowelli were smaller species and thus the average of the larger three species may underestimate M and the resulting reference points.

The mean fishing mortality rate corresponding to the maximum sustainable fishing mortality, u_msm, varied from 0.13 for A. laevis to 0.27 for A. eydouxii (mean = 0.20, s.d. = 0.04; Figure 3). The minimum unsustainable fishing mortality rate, u_crash, for these two species ranged from 0.24 to 0.46 (mean = 0.36, s.d. = 0.07). From the comparison between the estimated u and u_msm or u_crash, there was no species that appeared to be at risk of being either unsustainably impacted or risked becoming extinct in the long term at the current fishing intensity.

Temporal catch trend

In all, 3915 sea snakes were recorded and counted in scientific surveys and scientific and crew-member observer programmes between 1993 and 2007. Only one specimen was caught for each of two species, Enhydrina schistose and H. caerulescens, and eight specimens of H. mcdowelli. Therefore, the ZINB model could not be applied to these species. Sea snakes were rare in the surveys or commercial trawls (Figure 4). The maximum catch recorded in one trawl was eight sea snakes.

Figure 4.

Number of sea snakes (in the natural log scale for all species combined) captured in the Australian NPF. There is only one record for both seven and eight snakes in a sample.

Among the 11 species with sufficient data to apply the ZINB, six species were not sufficiently widespread to include bioregion as a covariate. Furthermore, season could not be included as a covariate for one species. For these species, data had to be combined across bioregions or seasons. Likelihood ratio tests indicated that the ZINB models fitted the data significantly better than the null model without predictors. Nearly all variables appeared to have a significant effect (p < 0.05) on sea snake catch. These effects were either on the probability of zero catch or the mean catch μ. The exceptions were that season was not significant for A. laevis, A. stokesii, and H. ornatus and the commercially fished area was not important for A. stokesii, Disteira kingii, D. major, and H. ornatus (all p > 0.05).

The effect of covariates on sea snake catches varied among species and did not show a consistent pattern. Generally, catch increased in shallower water and in commercially fished areas. As the main objective of this catch trend modelling was to investigate temporal trends in sea snake catches, we extracted the coefficient of the year effect for the 11 species that had sufficient data (Table 3). Year affected the catch of all 11 species for at least one component of the ZINB model, particularly the probability of 1 − p. A non-significant effect may be simply due to low statistic power to detect the trend as a consequence of small sample sizes. Observed catch per sample of most species varied widely from year to year (Figure 5). This did not necessarily reflect changes in overall catchability or abundance between years but may simply be because samples were taken at different areas or different seasons in different years. However, when these variables were held constant to separate a year effect, there was no consistent pattern among the 11 species (Table 3). For example, for some species (D. kingii, D. major, H. elegans, H. ornatus, H. pacificus, and all species combined), the probability of non-zero catch p increased over time but the mean catch μ declined, indicating greater dispersion over time. Four species, A. duboisii, A. eydouxii, A. peronii, and A. stokesii, showed the opposite trend. This indicated that individuals of these species became more aggregated over time. The joint effect was a curvilinear temporal trend in catch for most species. Only one species, L. hardwickii showed increases in both mean catch μ and probability of non-zero catch p over time. Overall, no obvious declining trend could be seen in the integrated mean catch from the product of the two components for these species.

Table 3.

Estimated year effect parameter in two components of the ZINB model: non-zero counts (which estimates mean catch μ) and the zero point mass (which estimates probability of zero catch).

Species	Component	Estimate	s.e.	Z-value	Pr (>|z|)
Acalyptophis peronii	Count	0.352	0.143	2.469	0.014
	Zero	0.886	0.233	3.812	0.000
Aipysurus duboisii	Count	0.522	0.335	1.559	0.119
	Zero	1.952	0.809	2.414	0.016
Aipysurus eydouxii	Count	0.153	0.064	2.397	0.017
	Zero	0.278	0.107	2.597	0.009
Aipysurus laevis	Count	0.228	0.073	3.114	0.002
	Zero	1.447	0.320	4.523	0.000
Astrotia stokesii	Count	0.177	0.072	2.437	0.015
	Zero	0.801	0.252	3.176	0.001
Disteira kingii	Count	−0.284	0.085	−3.350	0.001
	Zero	−1.075	0.558	−1.927	0.054
Disteira major	Count	−0.319	0.058	−5.462	0.000
	Zero	−0.719	0.095	−7.569	0.000
Hydrophis elegans	Count	−0.158	0.021	−7.502	0.000
	Zero	−0.931	0.255	−3.652	0.000
Hydrophis ornatus	Count	−0.198	0.048	−4.157	0.000
	Zero	−0.653	0.268	−2.436	0.015
Hydrophis pacificus	Count	−0.090	0.135	−0.667	0.505
	Zero	−0.950	0.396	−2.402	0.016
Lapemis hardwickii	Count	0.011	0.032	0.335	0.738
	Zero	−0.217	0.085	−2.548	0.011
All species	Count	−0.130	0.019	−6.834	0.000
	Zero	−0.476	0.044	−10.799	0.000

A positive estimate in the count component indicates an increase in catch in the non-zero samples over the years, whereas a positive estimate in the point mass at zero component indicates a increase in the probability of the zero catch over the year.

Figure 5.

The temporal catch trend of each sea snake species modelled with ZINB. Cross, mean observed catch per sample; thin solid line, model predicted catch; dashed line, predicted year effect on mean catch μ from the count part submodel; doted line, predicted year effect on the probability of non-zero catch p from the zero part submodel; thick solid line, predicted year effect on integrated mean catch E[c]. For A. laevis and D. kingii, the dashed lines overlap with the thick solid line, whereas the dotted lines are off the figure (>0.6). The observed catch does not necessarily reflect the change in catch rates or abundance because samples were taken at different areas or seasons each year.

Aside from the temporal trend, a test of the shape parameter κ in Equations (5 and 6) indicated overdispersion in spatial distribution (a tendency towards an aggregated distribution) for five out of the 11 species: A. laevis, D. major, H. elegans, H. ornatus, and L. hardwickii. Although the negative binomial model failed to detect a tendency to aggregate among the other six species, overdispersion was significant when all species were combined, noting a caveat of using κ to measure aggregation (Taylor et al., 1979).

Discussion

Our integrated risk assessment provides a new approach for examining the sustainability of rare species of conservation concern. It involves components of fishery stock-assessment and population viability analysis similar to those used for conservation (Akcakaya and Sjogren-Gulve, 2000) and catch trend modelling for species with low detection rates. This integrated approach will provide a more robust assessment of population risk for rarely caught marine species such as sea snakes, many chondrichthyans, and sea turtles. For these species, it is unlikely that sufficient catch data will ever be available to undertake more detailed population dynamics modelling.

The static sustainability assessment (SAFE; Zhou and Griffiths, 2008; Zhou et al., 2009) and population trend modelling complement each other and integrating both methods improves the robustness of the results. The outcomes from the SAFE analysis indicate whether a species can endure external pressure at that point in time. The assessment does not give an indication whether the population may have changed over time. The population trend modelling examines population fluctuations over time. However, it does not indicate whether the species is being unsustainably impacted by human activities. Integrating the two approaches allows the assessment of both fishing pressure and population state. The results show that species can change in their patterns of dispersion over time, independent of abundance. Besides fishing impacts, the inconsistency of catch trends among sea snake species may be due to natural population variations.

Our catch trend modelling suggests that a ZINB model is appropriate to the analysis of catchability of rare species such as sea snakes. Similar models have been recently applied to other rare marine species such as sharks and seabirds (Welsh et al., 2000; Minami et al., 2007; Shono, 2008). Comparisons among alternative models with similar data found two-component models such as ZINB to be appropriate (Potts and Elith, 2006). This is especially the case in situations where the presence and the non-detection of species (“false” zeros) were likely to be the most common source of non-detection (Martin et al., 2005).

Many studies have shown that the life-history traits of fish species have a close relationship with the resilience of populations to external pressure (Rochet, 2000; Denney et al., 2002; Dulvey and Reynolds, 2002). Since sea snakes are aquatic animals and have similar biological characteristics to elasmobranchs, we estimated sustainability reference points based on their natural mortality rates. We chose to be cautious and conservative and did not set maximum sustainable fishing mortality at the level of natural mortality, as is commonly applied to many fish species (Quinn and Deriso, 1999). Instead, we set the maximum sustainable fishing mortality at half that of the natural mortality rate and the minimum unsustainable fishing mortality at the level of natural mortality. This is the first application of this approach to sea snakes and requires further research to derive a more accurate relationship between life-history parameters and sustainability reference points.

We estimated natural mortality by applying four methods that were derived from empirical relationships developed for other marine animals, including fish, molluscs, and cetaceans. As the estimates were derived from a wide range of taxa, the relationships between life histories and M should also be applicable to sea snakes (Charnov et al., 1993). They are intermediate in size and longevity of those species used to derive the relationships. The results are comparable with natural mortality estimates of related tropical water and land snakes (Webb et al., 2003; Madsen et al., 2006). As there are no published natural mortality estimates for sea snakes or studies that link natural mortality with other life-history parameters, this is an interesting area for future research.

Our analysis showed that no species of sea snake was at risk of overfishing or being unsustainable under 2004–2006 levels of fishing effort and distribution. Furthermore, we adopt more precautionary approaches to risk. First, assigning sea snake q = 1 is very likely an overestimate of their catchability. Second, we have assumed no escapement of sea snakes from trawls. Although there are no measurements of the escapement rate, underwater video footage during previous studies has confirmed that sea snakes can escape from the path of a trawlnet and through meshes once caught (Brewer et al., 1998; Pitcher et al., 2007). These two assumptions are likely to have caused overestimation of fishing mortality. Third, the level of prawn fishing effort in the NPF has greatly reduced from ∼35 000 boat-days in the early 1980s to <6000 boat-days in 2007. As a result, the estimated annual catch of sea snakes has also dramatically declined (Wassenberg et al., 1994; Milton et al., 2009). Indeed, the catch trend modelling did not reveal a significant decline in the catch rate for any species. This indicates that the abundance of all species were probably stable. However, one component of the ZINB model outputs (either the probability of non-zero catch or the mean catch of the non-zero part) had a declining temporal trend for most species. Although it is the product of the two components (pμ) that reflects the catch trend, this phenomenon implies a change at least in the distribution pattern and should be monitored.

The estimated fishing mortality rates were specific to each species. Previous assessments of fishing impact from this fishery (Milton, 2001; Stobutzki et al., 2001) developed a relative index of impact based on life-history patterns. Milton (2001) suggested that two species might be of concern due to fishing impacts—D. kingii and H. pacificus. Our results confirm that H. pacificus has the highest fishing mortality rate. However, given the conservative sustainability reference points, the level of recent fishing mortality on H. pacificus does not appear to be unsustainable. The catch trend analysis supports our sustainability assessment. The risk to the H. pacificus populations has been further reduced in 2007 with the reduction in fleet size and fishing effort.

In contrast, the other species Milton (2001) identified to be of potential concern, D. Kingii, only had ∼10% of its distribution within areas of recent trawl effort. The estimated annual fishing mortality rate was much smaller than our conservative estimate of the MSM. This apparent decline in risk is partly due to the shift in fishing effort away from the southeastern part of the Gulf of Carpentaria where D. kingii was more abundant (Milton et al., 2009). If the fleet continues to fish in the same regions of recent high effort, the overall impact on D. kingii populations should remain low.

The integrated risk assessment with catch trend modelling approaches developed here are also directly applicable to many other animal groups such as insects, reptiles, amphibians, and fish with non-determinate growth. Previous assessments of ecological risk to large groups of species have mostly undertaken broad habitat or population assessments based on the expert opinion (Lenders et al., 2001; Burgess et al., 2006; Lewis and Senior, 2011). Many groups, such as insects, may have sufficient detection–non-detection data to undertake species assessments with the approach developed here. Another advantage will be that changes in the status of individual species in response to management actions should be more objectively measured than the qualitative methods using arbitrary high–low risk ranking (Griffiths et al., 2006).

Acknowledgements

We thank all the crew-members and scientific observers and CSIRO scientists who have collected sea snake data. The curators of the Herpetology sections of the Australian, Queensland, Northern Territory, and Western Australian Museums are thanked for providing sea snake records in their collections. We further thank Drs Shane Griffiths, Tony Courtney, Alec MacCall, and Will Le Quesne for their constructive comments on earlier drafts of the manuscript. This project was jointly funded by CSIRO and the Fisheries Research and Development Corporation (Project 2005/051).

Author notes

Handling editor: Rochelle Seitz