ABSTRACT Title of Dissertation: BIOLOGY, LIFE HISTORY AND CONSERVATION OF ELASMOBRANCHS WITH AN EMPHASIS ON WESTERN ATLANTIC SKATES Michael Guise Frisk, Doctor of Philosophy, 2004 Dissertation Directed By: Associate Professor Dr. Thomas J. Miller University of Maryland Center for Environmental Science Chesapeake Biological Laboratory In this dissertation two approaches were used to increase the knowledge of elasmobranch population dynamics and life history: (1), the comparative approach and (2), the species-specific approach. In the comparative approach I constructed standardized three-stage matrix models for 55 species of sharks and rays. Using these models I (1) conducted elasticity analyses to determine how the vital rates of mortality (M) and fertility (f) influence elasmobranch population growth rate r, (2) estimated sensitivity of elasticity to perturbation in vital rates, and (3) examined the taxonomic distribution of model inputs and species vital rates, such as size at maturity (L mat ), and total length (L max ). I found positive relationships between the elasticity of ? (population growth rate) to changes in juvenile and adult stages to longevity and age at maturity; however, the age at maturity and the elasticity of ? to changes in the adult stage relationship appeared to be invariant. Combining vital rates and elasticities, I found similar suites of life histories and demographics within taxonomic groups at various levels. Further I examined where (or if) elasmobranchs fall in the evolved triangular ordination of life history strategies proposed by Winemiller and Rose (1992). My results indicate that when plotted using only the teleost ordination, elasmobranchs appear to be periodic strategists, outside the limits of the teleost ordination. However, when elasmobranch data is included in the ordination they form the extreme range of equilibrium strategists and are grouped by order. In the species-specific approach, I found evidence for a strong latitudinal trend in maximum size (l ? ) and size at maturation (l mat ) in little skate with individuals in northern regions reaching a larger size at maturity and maximum length and growing slower than little skate from more southern regions. No similar trend was found in winter skate. Little skate is smaller, reaches maturity at a younger age is faster growing and shorter lived then winter skate (Little skate: l ? = 56.1 cm, k = 0.19/yr, T max = 12.5, T mat = 7; Winter skate: l ? = 122.1 cm, k = 0.07/yr, T max = 20.5, T mat = 12.5). Winter skate has higher annual fecundity then little skate of 26-101 and 21-57 eggs per year respectively. Using estimated vital rates for winter skate and National Marine Fisheries Service?s survey data an age-structured model was constructed for winter skate from 1963-1998. The model indicated that the western Atlantic population of winter skate was rebuilding in the 1980?s following overfishing in the 1960?s and 1970?s. BIOLOGY, LIFE HISTORY AND CONSERVATION OF ELASMOBRANCHS WITH AN EMPHASIS ON WESTERN ATLANTIC SKATES By Michael Guise Frisk Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2004 Advisory Committee: Dr. Thomas J. Miller, Chair Dr. Michael J. Fogarty Dr. Michael H. Prager Dr. Christopher L. Rowe Dr. Robert E. Ulanowicz ii ? Copyright by Michael Guise Frisk 2004 ii FOREWARD According to the rules and regulations of the Graduate School of the University of Maryland, a graduate student may co-author work with faculty and colleagues that is included in a Dissertation. In such an event the University requires that a forward be included in the Dissertation, as approved by the Dissertation Committee, that clearly defines the substantial contributions the candidate made to the relevant aspects of the jointly authored work included in the Dissertation. In accordance with the above, the Dissertation Committee approved the following statement of the candidate?s intellectual contribution to the dissertation. Chapters 1 and 3 are entirely the candidate?s own work. All of the remaining chapters are co-authored. Dr. Thomas Miller is co-author on all of these chapters. Dr. Miller offered constructive criticism on the design, analysis, interpretation and effective written communication of the research in addition to providing logistic support. The contributions of Dr. Miller to the Dissertation fall within the normal bounds of graduate supervision. Chapter 2 is co-authored with Dr. Miller and Dr. N. K. Dulvy, Center for Environment, Fisheries and Aquacultural Science in Lowestoft, U.K. Dr. Dulvy provided advice on multivariate statistical analysis. Chapter 6 is co- authored with Dr. Miller and Dr. S. J. D. Martell, Fisheries Centre, University of British Columbia, Canada. Dr. Martell provided advice on age-structured population models and guidance in programming in AD-Model Builider. Further, many important editorial suggestions were provided by committee members Drs. Michael J. Fogarty, Michael H. Prager, Christopher L. Rowe and Robert E. Ulanowicz. iii DEDICATION This work is dedicated to my family and friends who have provided guidance and encouragement in the completion of my education. I would especially like to thank my mom and dad. The effort and contents of this work are dedicated to the memory of my Granddad Jean Gaitley Guise, JR, (August 26 th 1908 ? August 14 th 2004). He instilled in me, respect for the natural world, science and the pursuit of knowledge. I was especially lucky to enjoy his presence for so many years. iv ACKNOWLEDGEMENTS This work could not have been completed without the excellent guidance provided from Dr. T.J. Miller. I would like to thank my committee for their dedication to my education and this project. There are too many people to thank than space would allow. However, I would like to thank all the people that worked or volunteered aboard the NOAA RV Albatross IV, Woods Hole, Massachusetts, who helped collect my samples and provided many good times. I would also like to thank all the people in the Miller lab group and in the CBL community who helped edit my work and assisted in the lab and field. The project was funded by the National Marine Fisheries Service/Sea Grant fellowship in population dynamics. Thanks for all the help! v TABLE OF CONTENTS Foreward......................................................................................................................ii Dedication ...................................................................................................................iii Acknowledgements ....................................................................................................iv Table of Contents ........................................................................................................ v List of Tables ............................................................................................................viii List of Figures.............................................................................................................. x Chapter 1: Introduction .............................................................................................. 1 Sustaining elasmobranch resources: The comparative approach...................... 3 Sustaining elasmobranch resources: The species-specific approach................ 5 Summary........................................................................................................... 9 Chapter 2: Life histories and vulnerability to exploitation of elasmobranchs: Inferences from elasticity, perturbation and phylogenetic analyses ........................... 10 Introduction..................................................................................................... 10 Methods........................................................................................................... 14 The Data............................................................................................... 14 The Models .......................................................................................... 14 Elasticity analyses................................................................................ 17 Sensitivity of elasticity......................................................................... 17 Response to exploitation...................................................................... 18 Phylogenetic Analysis.......................................................................... 19 Results............................................................................................................. 20 Variation in elasticity with life history traits ....................................... 20 Sensitivity of elasticity......................................................................... 21 The response of elasticities to exploitation.......................................... 22 Variation of life histories and elasticities with phylogeny................... 22 Discussion....................................................................................................... 25 Tables.............................................................................................................. 32 Figures............................................................................................................. 36 Chapter 3: Ordination of evolved life history strategies in elasmobranch and teleost fishes: In search of common ground........................................................................... 49 Introduction..................................................................................................... 49 Methods........................................................................................................... 52 vi Data....................................................................................................... 52 Statistical analysis................................................................................. 54 Results............................................................................................................. 56 Discussion....................................................................................................... 58 Tables.............................................................................................................. 61 Figures............................................................................................................. 65 Chapter 4: Age, growth, and latitudinal patterns of two Rajidae species in the notherwestern Atlantic: Little skate, Leucoraja erinacea and winter skate, Leucoraja ocellata........................................................................................................................ 67 Introduction..................................................................................................... 67 Methods........................................................................................................... 71 Sampling ............................................................................................... 71 Length-weight analyses ........................................................................ 72 Slide preparation and age determination............................................... 73 Ageing little skate and winter skate...................................................... 74 Precision and reader bias ...................................................................... 76 Centrum and growth analysis................................................................ 76 Results............................................................................................................. 78 Species identification............................................................................ 78 Length-weight relationship ................................................................... 78 Image quality and annuli identification ................................................ 79 Ageing little skate and winter skate...................................................... 79 Precision and reader bias ...................................................................... 80 Centrum and growth analyses............................................................... 80 Discussion....................................................................................................... 83 Tables.............................................................................................................. 90 Figures............................................................................................................. 93 Chapter 5: Maturation, fecundity, latitudinal patterns, and reproductive strategies of two Rajidae species in the northwestern Atlantic: Little skate, Leucoraja erinacea and winter skate, Leucoraja ocellata............................................................................... 105 Introduction................................................................................................... 105 Methods......................................................................................................... 111 Sampling ............................................................................................. 111 Laboratory procedures ........................................................................ 112 Maturation analysis............................................................................. 112 Fecundity analyses.............................................................................. 114 Optimal age of maturation .................................................................. 114 Results........................................................................................................... 116 Maturation........................................................................................... 116 Fecundity............................................................................................. 118 vii Optimal age of maturation .................................................................. 119 Discussion..................................................................................................... 120 Tables............................................................................................................ 127 Figures........................................................................................................... 132 Chapter 6: An age-structured model of winter skate, Leucoraja ocellata, abundance in the western Atlantic: sustainability and uncertainty ???????????150 Introduction................................................................................................... 150 Methods......................................................................................................... 153 Fishery-independent surveys ............................................................... 153 Model data inputs................................................................................. 154 Population dynamics............................................................................ 155 Model initialization.............................................................................. 156 Observation models ............................................................................. 159 Length frequency analysis ................................................................... 159 Likelihoods .......................................................................................... 161 Estimated parameters........................................................................... 162 Trial runs and ?best? model selection.................................................. 162 Compensatory recruitment................................................................... 163 Results........................................................................................................... 165 Base model runs and model selection.................................................. 165 Model 1 ................................................................................................ 165 Model 9 ................................................................................................ 168 Recruitment anomalies......................................................................... 169 Model generalization ........................................................................... 170 Discussion..................................................................................................... 171 Tables............................................................................................................ 180 Figures........................................................................................................... 185 Chapter 7: Summary .................................................................................. 195 Appendix A.............................................................................................................. 200 Appendix B .............................................................................................................. 206 Appendix C.............................................................................................................. 209 References ............................................................................................................... 213 viii LIST OF TABLES Table 2.1. Data for matrix analyses ........................................................................... 32 Table 2.2. Relationships of vital rates of elasmobranchs species with elasticity of model parameters........................................................................................ 34 Table 2.3. The results of varying exploitation levels and elasticity patterns for shark and skate species......................................................................................... 35 Table 3.1. The mean, standard deviation, and number of observations are shown for age of maturation (year) (T mat ), length of maturation (cm) (L mat ), longevity (years) (T mat ), mean clutch size (f mean ), egg size (ES), number of spawning bouts in a year (SB) and parental care (PC)........................... 61 Table 3.2. Pearson correlations used in the decision to remove uninformative variables from the ordination....................................................................... 62 Table 3.3. Principal components procedure for five variable ordinations for marine and freshwater teleost species.......................................................... 63 Table 3.4. Principal components procedure for five variable ordinations for marine, freshwater and elasmobranch species .......................................................... 64 Table 4.1. General linear model results for little skate and winter skate regional weight vs. length relationships..................................................................... 90 Table 4.2. Little skate and winter skate growth parameters ...................................... 91 Table 4.3. Correlations between parameters estimates of the von Bertalanffy growth equation where, l ? is the asymptotic length, t 0 is the theoretical size at birth and k is the growth rate................................................................................ 92 ix Table 5.1. A. Regression summary and parameter estimates for ovary weight vs. length relationships in little skate and winter skate .................................. 127 Table 5.2. Summary of logistic regression analyses for regional maturity comparisons of little skate ........................................................................ 128 Table 5.3. Fecundity estimates for varying mature egg size and spawning season for little skate.................................................................................................. 129 Table 5.4. Fecundity estimates for varying mature egg size and spawning season for winter skate ............................................................................................... 130 Table 5.5. Seasonal mature egg production in little skate and winter skate ............. 131 Table 6.1. Input data, definitions of parameters and symbols used in..................... 180 Table 6.2. National Marine Fishery Service?s annual fall and spring bottom trawl surveys, reported commercial landings and groundfish effort on southern New England and Georges Bank ............................................................ 182 Table 6.3. A. Results of model runs including parameters estimates K, R 0 , B 0 , ? 1 and ? 2 ............................................................................................................... 183 Table 6.4. Correlation to R0 for varying lower bounds of....................................... 184 x LIST OF FIGURES Figure 2.1. The relationship between elasticity of ? to changes in the juvenile stage and longevity (T max ) .............................................................................. 36 Figure 2.2. The relationship between elasticity of ? to changes in the adult stage and longevity (T max ) ..................................................................................... 37 Figure 2.3. The relationship between elasticity of ? to changes in the juvenile stage and age of maturity (T mat ) ..................................................................... 38 Figure 2.4. The relationship between elasticity of ? to changes in the adult stage and age of maturity (T mat ) ............................................................................ 39 Figure 2.5. The relationship between elasticity of fertility and the transition between stages and longevity (T max ) ..................................................................... 40 Figure 2.6. The relationship between elasticity of ? to changes in fertility and the transition between stages and age of maturity (T mat ).............................. 41 Figure 2.7. The sensitivities of elasticity of the adult stage are shown for species in the superorder Galea for each stage of the matrix .................................. 42 Figure 2.8. The sensitivities of elasticity of the adult stage are shown for species in the superorder Batoidea for each stage of the matrix ............................. 43 Figure 2.9. The sensitivities of elasticity of the juvenile stage are shown for species of the superorders Batoidea and Galea for each stage of the matrix .... 44 Figure 2.10. MDS ordination of elasmobranch life histories and elasticities............ 45 Figure 2.11. MDS ordination of elasmobranch life histories and elasticities showing the relative magnitude of each trait........................................................ 46 xi Figure 2.12. MDS ordination of elasmobranch life histories and elasticities split by phylogenetic groupings, superorder, order and family .......................... 47 Figure 2.13. MDS ordination of elasmobranch life histories and elasticities overlaid with the reproductive mode of each species .......................................... 48 Figure 3.1. The principal components and correlations of the first two axes for teleost loadings for the five variable ordination................................................. 65 Figure 3.2. The principal components and correlations of the first two axes for teleost and elasmobranch loadings for the five variable ordination................... 66 Figure 4.1. Sampling locations at which little skate vertebra were sampled along the N.E. coast................................................................................................. 93 Figure 4.2. Locations at which winter skate vertebra were sampled along the N.E. coast ......................................................................................................... 94 Figure 4.3. Allometric relationship between length (cm) and weight (kg) for little skate ......................................................................................................... 95 Figure 4.4. Allometric relationship between length (cm) and weight (kg) for winter skate. Data are shown for N.E. coast ...................................................... 96 Figure 4.5. Sample images of centra for A) the birthmark of a winter skate, B) and C) two examples of winter skate centra sections and D) a little skate centra section ..................................................................................................... 97 Figure 4.6. von Bertalanffy growth model fit to length (cm) and age (years) data for little skate................................................................................................ 98 Figure 4.7. Relationship between von Bertalanffy asymptotic size (l ? ) and Brody growth coefficient (k) shown for little skate from three regions ............ 99 xii Figure 4.8. von Bertalanffy growth model fit to length (cm) and age (year) data for male and female little skate from the N.E. coast .................................. 100 Figure 4.9. Comparison of current estimates of age/growth rates and previous studies ............................................................................................................... 101 Figure 4.10. von Bertalanffy growth model fit to length (cm) and age (years) data for winter skate ......................................................................................... 102 Figure 4.11. von Bertalanffy growth model fit to length (cm) and age (year) for sex specific data for winter skate .................................................................................... 103 Figure 4.12. Comparison of growth models fit to liberal, conservative, and high readability data..................................................................................... 104 Figure 5.1. Regional coverage of little skate ovaries collected along the northeastern United States coast from Cape Hatteras to Canadian waters................. 132 Figure 5.2. Regional coverage of winter skate ovaries collected along the northeastern United States coast from Cape Hatteras to Canadian waters ................................................................................................................ 133 Figure 5.3. The reproductive track of a little skate shown with two ovaries and two shell glands............................................................................................ 134 Figure 5.4. Ovary weight (g) regressed over individual female total length for little skate ...................................................................................................... 135 Figure 5.5. Ovary weight (g) regressed over individual female total length for winter skate ....................................................................................................... 136 Figure 5.6. Number of reservoir and mature eggs for little skate (eggs > 5 mm) by size of skate........................................................................................... 137 xiii Figure 5.7. Number of reservoir and mature eggs for winter skate (eggs > 5 mm) by size of skate............................................................................................ 138 Figure 5.8. Physiological maturity regressed over length for little skate is shown for the mid-Atlantic, southern New England-Georges Bank and the Gulf of Maine .................................................................................................... 139 Figure 5.9. Proportion mature and age for given lengths are shown for physiological and functional maturity in winter skate ................................................. 140 Figure 5.10. Mature eggs and proportion of females with mature eggs are shown for female little skate (eggs > 10 mm)....................................................... 141 Figure 5.11. Mature eggs and proportion of females with mature eggs are shown for female winter skate (eggs > 15 mm)................................................... 142 Figure 5.12. Proportion mature and age for given lengths are shown for physiological and functional maturity for the mid-Atlantic, southern New England- Georges Bank and the Gulf of Maine for little skate.......................... 143 Figure 5.13. A 3-D representation of the number of eggs in female little skate plotted against egg size for mature females..................................................... 144 Figure 5.14. A 3-D representation of the number of eggs in female winter skate plotted against egg size with mature females ...................................... 145 Figure 5.15. Net reproductive rate (R0) is shown for ages of maturity ranging from 0- 12 for little skate ................................................................................. 146 Figure 5.16. Net reproductive rate (R0) is shown for ages of maturity ranging from 0- 16 for winter skate ............................................................................... 147 Figure 5.17. Age specific egg production with fishing and no fishing (F= 0.35) for xiv little skate............................................................................................. 148 Figure 5.18. Age specific egg production with fishing and no fishing (F= 0.35) for winter skate ......................................................................................... 149 Figure 6.1. Survey CPUE are shown for the fall is CPUE for inshore and offshore regions for the Gulf of Maine, Georges Bank, and southern New England ............................................................................................................... 185 Figure 6.2. Survey aggregated CPUE is shown for the fall survey ......................... 186 Figure 6.3. The length frequencies in three centimeter intervals of the fall survey for offshore regions from 1963 to the present ............................................ 187 Figure 6.4. Predicted biomass and survey CPUE for model 1 for 1963 to 1998..... 188 Figure 6.5. Predicted biomass and commercial CPUE for model 1 for 1963 to 1998 ............................................................................................................... 189 Figure 6.6. Model 1 results assuming a lower bound for B0 of 0.5 for 1963 to 1998 ............................................................................................................... 190 Figure 6.7. Estimated fishing mortality for 1963 to 1998........................................ 191 Figure 6.8. Predicted biomass and survey CPUE for model 9 for 1963 to 1998..... 192 Figure 6.9. Predicted biomass and commercial CPUE for model 9 for 1963 to 1998 .............................................................................................................. 193 Figure 6.10. Results of model 1 with process error in recruitment for 1963 to 1998 .............................................................................................................. 194 1 Chapter 1: INTRODUCTION Chondrichthyans first appeared in the Silurian period (350 million ybp) and are believed to have evolved in marine environments (Moyle and Cech, 2000). Chondrichthyan fishes comprise two subclasses: the Elasmobranchii and Holocephali (Moyle and Cech, 2000). The focus of this work is on the elasmobranchs (sharks, skates and rays), which are characterized by a cartilaginous skeleton, the lack of ossified otoliths, the presence of claspers for reproduction in males and internal fertilization. The dissertation has two goals: to understand life history patterns within elasmobranchs generally, and to understand the population dynamics of two representative species of elasmobranchs ? winter skate and little skate. Evolutionary evidence suggests that egg laying is the ancestral form in elasmobranchs with most species transitioning to live-bearing. In elasmobranchs, live-bearing is the form of parity in 60% of species, egg laying is found in 25% and the rest have some combination of the two (Dulvy and Reynolds, 1997). In Rajidae (skates and rays) it appears that egg laying is also an ancestral form, however, skates transitioned to live bearing and then back to egg laying (Dulvy and Reynolds, 1997). In general, elasmobranchs have very low fecundity compared to most teleost species, delayed maturity and long life-spans (Frisk et al., 2001). The life history of elasmobranchs, especially for large species, places them at risk of overexploitation (Frisk at al., 2001). Current research suggests that many of the large predatory shark species of the western Atlantic have dramatically declined due to commercial and recreational 2 fishing (Baum et al., 2003; Myers and Worm, 2003). The recent declines in skate and shark populations in the western Atlantic are believed to have resulted from increases in commercial, recreational and by-catch harvests (NEFSC, 1999). Casey and Myers (1998) described the possible extirpation of the barndoor skate, Dipturus laevis, in the western Atlantic. Additionally, the 1999 National Marine Fisheries Service (NEFSC, 1999) assessment of the Northeast region skate complex highlighted the overfishing of winter skate and barndoor skate and the decline of thorny skate, Amblyraja radiata, (NEFSC, 1999). This report also expressed concern about the paucity of information on the vital rates of skates (Rajidae), a problem common in efforts to understand the resilience of elasmobranch populations to harvest (Frisk et al., 2001). While declines in many shark and skate populations have been well documented, much is left to be accomplished in the development of sound methodology for their conservation. In addition, there are many differences between elasmobranchs and bony fish that make it challenging to fit them into present population dynamic models. Fish populations have evolved the capacity to overcompensate reproductive output in the face of variability in annual mortality (Pitcher and Hart, 1982). This surplus production provides for the sustainable harvests (Pitcher and Hart, 1982). Elasmobranch fishes exhibit egg-laying or live-bearing and produce large offspring, presumably with high survival rates, and have delayed maturity. These facts combine to produce low surplus productions in elasmobranch populations. Using Jenning?s et al. (1998) index of potential reproductive output at the onset of maturation, Frisk et al. (2001) showed that elasmobranchs had much lower potentials for stock 3 replenishment at the onset of maturation compared to teleost species. Moreover, high offspring survival and low fecundity likely means stock recruitment relationships for sharks, skates and rays are strong and less variable than for teleost species (Smith et al., 1998). Nothing in the above paragraphs would eliminate the use of common fishery models from being used for elasmobranchs. After all, in order for a species to persist it must have at least a minimal capacity to recover from population declines (i.e. surplus production). However, what it does suggest is that elasmobranchs are a group that are often not adequately managed and have not been fully examined in terms of life history theory and population dynamics. Importantly, the comparative method that has yielded considerable insight into the ecology and management of teleost species (e.g., Beverton and Holt, 1959; Beverton, 1963; Pauly, 1978; Hoenig, 1983; Beverton, 1992) has yet to be fully developed for elasmobranchs. I suggest that combining insights from comparative analyses with analyses of individual species will lead to the development of new conservation measures for elasmobranchs. The remainder of this introduction and the details provided in this dissertation will highlight the approaches I will use to improve our understanding of life history traits, vital rates, and the population dynamics of elasmobranchs and the development of conservation measures for the subclass. Sustaining elasmobranch resources: The comparative approach Comparative analyses have proven to be a powerful approach for assessing the sustainability of species and the study of life history evolution for diverse taxa 4 (Harvey and Pagel, 1991; Stearns, 1992; Charnov, 1993). For example, ecologists have used the behavior of the logistic equation to examine the life histories and evolutionary dynamics of species from a broad range of animal taxa (Pianka, 1970; Mertz, 1975; Caswell, 1982; Hall, 1988). The use of comparative analysis in fishery science dates back to pioneers of the field including Beverton and Holt (1959). Later, Pauly (1978) and Hoenig (1983) used the approach in their popular analyses for estimating mortality rates. More recently, Winemiller and Rose (1992) provided a classification scheme for teleost taxa, without a prior theoretical model, based on a comprehensive analysis of life history traits of fishes. In my dissertation I will conduct two comparative analyses of elasmobranch life histories to provide guidance to managers in their efforts to conserve elasmobranch stocks. In the second chapter, I present results of an analysis of vital rates of elasmobranchs using matrix-based projection models. The analyses are based on previously published estimates of vital rates of 45 elasmobranch species. Stage-based projection models are used to estimate the intrinsic rate of natural increase, r. Subsequently, I use elasticity and sensitivity of elasticity analyses to quantify the relative susceptibility to exploitation of the species and use elasticity to provide information on important aspects of the species life history. This work has been submitted and accepted for publication in the Journal of the Northwest Atlantic Fisheries Organization. The paper is authored by myself, Thomas Miller and Nicholas Dulvy, a colleague at the Center for Environment, Fisheries and Aquacultural Science in Lowestoft, UK. 5 In the third chapter, I present a second example of the application of the comparative approach to understanding elasmobranch life histories. Using multivariate statistical analyses, Winemiller and Rose (1992) expanded the r-K continuum to a triangular ordination of three evolved strategies in teleost fishes: (1) small, rapidly maturing, short-lived species (opportunistic strategists), (2) larger, highly fecund fishes with longer life-spans (periodic strategists), and (3) fishes of intermediate size that often exhibit parental care and produce fewer but large offspring (equilibrium strategists). Few studies on the life history patterns of elasmobranchs exist in the literature (Cortes, 2000; Frisk et al., 2001). To address the need of further comparative life history analysis for elasmobranchs, I expanded the Winemiller and Rose analysis by the addition of data for skates, sharks and rays. I estimated where elasmobranchs fit into the combined teleost and cartilaginous fish ordination. I anticipate that this work will be submitted as a manuscript to either Transactions of the American Fisheries Society or to Copeia as a sole-authored paper. Sustaining elasmobranch resources: The species-specific approach Frisk et al. (2002) combined data on little skate, barndoor skate and winter skate and parameter estimates from an empirical analysis of elasmobranch life histories (Frisk et al., 2001), and used a matrix model to estimate fishing limits for the three western Atlantic skate species. The approach of Frisk et al. (2002) utilized both known and empirically estimated vital rates for size at maturity, age at maturity and natural mortality; parameters that are necessary to properly set conservation measures. Still, without field-derived data for individual species, sound conservation 6 polices are difficult to form. Comparative analyses are a proven resource; however, nothing is more reliable than species-specific, field-derived data. Little skate and winter skate occur in the western Atlantic from the mid- Atlantic region to Canadian waters (McEachran, 2002). Skates are caught both as by- catch and targeted in bottom trawl fisheries in the Northern Atlantic. Skates experienced changes in the harvesting pressure over the last 50 years with harvest peaks in the 1970?s and during the late 1980?s and 1990?s (Fogarty and Murawski, 1998). During the 1960?s and 1970?s, skates were targeted largely by foreign fishing fleets (Overholtz and Tyler, 1986; Fogarty and Murawski, 1998). The United States extended its jurisdiction to the 200-mile limit in 1976, effectively removing all foreign fishing. The result was lower landings of skates during the late 1970?s and early 1980?s until the domestic fleet began keeping skates for the market during the late 1980?s and 1990?s (NEFSC, 1999). During the 1980?s, following the removal of foreign fishing pressure, the NMFS annual survey indicated that winter skate exhibited a ten-fold increase in abundance (NEFSC, 1999). Little is known of actual harvest rates of skates in the western Atlantic. Length-based estimates (Beverton and Holt, 1959; Hoenig, 1983) indicate recent fishing mortality is approximately 0.3-0.4 for little skate, below the F 0.1 threshold of 0.71, and 0.2-0.3 for winter skate, above the F 0.1 threshold of 0.1-0.2 (NEFSC, 1999). Very little is known of the vital rates in little skate and winter skate. Six age and growth studies have been conducted on western Atlantic skate species. Johnson (1979), Richards et al. (1963), Waring (1984) and Natanson (1990) provide data on little skate. Sulikowski et al. (2003) and Simon and Frank (1996) provide 7 information on winter skate. Gedamke, (on going) is currently working on barndoor skate. Knowledge of egg production in skates in the western Atlantic is also not well developed. Studies by Johnson (1979) and Richards et al. (1963) are the only attempts to estimate annual fecundity for western Atlantic skates. The effects of latitude and the resulting environmental conditions on the vital rates in terrestrial animals have been explained by Allen?s and Bergman?s rules. These rules note general trends apparent in terrestrial animals: as latitude increases, appendage size decreases, while body size increases. In many marine fish, similar trends have been observed with species? populations in higher latitudes exhibiting slower growth, later age at maturity, longer life-span and increased longevity (Taylor, 1958; Beverton and Holt, 1959; Beverton, 1992). Latitudinal variation in vital rates and energetics is not always as straightforward as the examples presented above. Conover (1990) suggests that the potential for growth varies inversely with latitude and the length of the growing season. This counter-gradient variation has been hypothesized in American shad, Alosa sapidissima, striped bass, Morone saxatilis, and the mummichog, Fundulus heteroclitus. In these species adult body size is independent of latitude while the growing season is shorter in higher latitudes. This suggests that populations in higher latitudes grow faster during a shortened growing season whereas populations in lower latitudes grow slower for a longer duration, with both populations reaching the same size at the completion of the first year of life. Evidence suggests that in higher latitudes individuals of little and winter skate reach a larger size, indicating that they may have slower growth, delayed maturity, and longer life-spans. Length at maturity in winter skate varies from 70 to 109 cm, with 8 individuals maturing at larger lengths in higher latitudes (McEachran, 2002). Total length also varies in skates with larger species found in more northern populations. Delayed age/size at maturity, low population growth rate and values of the intrinsic rate of population increase all have been associated with a low resilience to exploitation in elasmobranch species (Frisk et al., 2001). A central question of my dissertation is whether skates exhibit latitudinal gradients in their life histories. In Chapters 4 and 5 I present results of the analysis of age at maturation and growth rates (Chapter 4) and fecundity (Chapter 5) of little skate and winter skate collected from Cape Hatteras, N.C. to Canadian waters. It is anticipated that both papers will be submitted for publication to the Canadian Journal of Fisheries and Aquatic Sciences as jointly authored papers with Thomas Miller. In the final chapter of the dissertation I use the data developed in Chapters 4 and 5 and data derived from the fishery directly (catch/effort data) to understand the population dynamics of winter skate. I use an age-structured model based on National Marine Fisheries Service?s annual fall and spring bottom trawl survey data, landings data, effort data for the groundfishery and vital rates to model winter skate?s relative abundance trends from 1963-1998. The model is used to predict biomass from 1963 to 1998 to observe how winter skate responded to varying fishing levels over the time series. I use the model to explain the historic fishing and abundance patterns of winter skate and address possible compensatory mechanisms in the species. It is anticipated that this chapter will be submitted for publication to the Canadian Journal of Fisheries and Aquatic Sciences as a paper coauthored by myself, Drs. Steve Martell (Uuniversity of British Columbia), and Thomas Miller. 9 Summary I have undertaken two comparative analyses to enhance the conservation of elasmobranchs: (1), using a three stage matrix model, analyzed elasticity trends in 45 elasmobranch species (Chapter 2); and (2), compared the vital rates in elasmobranchs in relation to the Winemiller and Rose (1992) ordination of evolved life histories in teleost species (Chapter 3). I have developed a better understanding of the population dynamics of little and winter skates in the western Atlantic to improve conservation measures for elasmobranchs by: (1), estimating life history parameter estimates for age, growth, maturation (Chapter 4) and fecundity (Chapter 5) for little skate and winter skate: (2), using survey data, landings data, and effort data to develop an age-structured population dynamics model for winter skate from 1963-1998 (Chapter 6). 10 Chapter 2: LIFE HISTORIES AND VUNERABILITY TO EXPLOITATION OF ELASMOBRANCHS: INFERENCES FROM ELASTICITY, PERTURBATION AND PHYLOGENETIC ANALYSES Introduction Several species of elasmobranchs have been shown to be vulnerable to population declines and even local extinction (Casey and Myers, 1998; Stevens et al., 2000; Simpfendorfer, 2000; Dulvy et al., 2000; Frisk et al., 2001; 2002). The key parameter for determining the vulnerability of a species to population declines when exploited is the intrinsic rate of population increase, r. Species exhibiting a high r are more resilient to exploitation and likely recover more rapidly once harvesting ceases than species with a low r. The degree of vulnerability of individual species has also been linked to life history traits, such that species exhibiting the combination of large maximum body size, slow growth, late maturation (at a large size), and long lifespan appear to be most vulnerable (Walker and Hislop, 1998; Dulvy et al., 2000; Stevens et al., 2000; Frisk et al., 2001). Life histories are constrained by trade-offs; slow- growing species tend to be large bodied and mature later in life and have lower annual reproductive output (Charnov, 1993; Reynolds et al., 2001; Frisk et al., 2002; Roff, 2002). It should not be surprising then to find that species with 'slow' life histories also have low r values (Fenchel, 1974; Musick, 1999; Denney et al., 2002). However, the link between life history and population dynamics and specifically the link to the population-level response of the additional mortality resulting from exploitation remains unclear. 11 Demography, the schedule of survival and reproduction of each age class or life stage in a population, links life history and population dynamics. Matrix models are used often to understand demography because they provide both a convenient method for integrating vital rates (survival and fertility) and extrinsic anthropogenic factors such as exploitation or pollution across age or stage classes, and a means to calculate parameters useful for understanding population dynamics, e.g. (? = e r ) (Cort?s, 1998; Walker and Hislop, 1998; Heppell et al., 1999; Brewster-Geisz and Miller, 2000; Caswell, 2001; Cort?s, 2002; Frisk et al., 2002; Mollet and Cailliet, 2002). However, estimates of the rate of population growth vary with population density. Consequently estimates of population growth rate that currently characterize exploited populations are unfortunately not necessarily representative of the performance of virgin populations (Jennings et al., 1998; Smith et al., 1998; Reynolds et al., 2001; Cort?s, 2002). Sensitivity analysis of such matrix models can be used to help identify which life history stages contribute most to variation in the population growth rate (r). Two common forms of analyses are typically used: sensitivity and elasticity (Benton & Grant, 1999). Sensitivity measures the effect of an absolute change in a vital rate upon population growth rate, whereas elasticity measures the effect of a proportional change in a vital rate on population growth rate (Benton & Grant, 1999; Caswell, 2001). Both sensitivity and elasticity elucidate critical aspects of a species life history, provide insight for the focus of natural selection and indicate where management actions may be most successful. When using sensitivity or elasticity analyses, exact estimates of population growth rates are not needed to understand 12 how management actions, which act via manipulation of vital rates, will generally influence population growth. However, the sensitivities and elasticities calculated only apply to the initial vital rates used to define the projection matrix. To understand how elasticity changes as vital rates vary, the sensitivity of elasticity must be estimated (Caswell, 1996). Sensitivities of elasticity allow for analysis of the population-level consequences of dynamic changes in vital rates that result from perturbation of the vital rates. The pattern and extent of responses to external factors exhibited by an individual species may be constrained by its phylogeny. Species that share a subset of traits derived from a common ancestor often exhibit similar suites of life histories even though they may live in very different habitats (Pagel and Harvey, 1988; Harvey and Pagel, 1991). Thus, phylogenetic influences on variations in life histories and demography must be considered. Elasmobranchs are divided between two contrasting superorders: Galea and Batoidea. The Galea is represented in my analyses by predatory, shallow-water species of requiem sharks (Carcharhinidae) and hound sharks (Triakidae). The Batoidea includes both skates (Rajidae) and rays (Myliobatidae). Here I examine links between elasmobranch demography and life history, using comparative analysis of the outputs of a standardized stage-based matrix projection model parameterized for 55 elasmobranch species. Specifically I: (1) examine how vital rates (juvenile and adult survival and fertility elasticities) vary with life history traits (longevity, age of maturity and body size), (2) examine how elasticities respond to varying levels of exploitation, using four representative species, 13 and (3) test whether suites of life histories and demographic traits are linked to phylogenetic relationships. I report results at the superorder, order and family levels. My results are not intended to estimate limits to exploitation, but rather to explore how different species and phylogenetic groups of species are potentially influenced by exploitation. I hope to add to the discussion of where in the elasmobranch life cycle potential compensatory responses may occur under exploitative or environmental changes and how these may differ across species groupings. My primary goal is to link elasmobranch phylogeny, life history, and conservation. 14 Methods The Data I extracted estimates of age of maturity (T mat ), longevity (T max ), and fecundity (F) from the literature (Table 1). All published data I obtained were included in my analyses. Adequate data were available for 55 species from 12 families, namely: Carcharhinidae (22 species), Rajidae (9), Triakidae (9), Alopiidae (3), Lamnidae (3), Sphyrnidae (2), Squalidae (2 stocks), Urolophidae (2), Dasyatidae (1), Odontaspididae (1), Scyliorhinidae (1), and Myliobatidae (1). Reported estimates of age of maturity (T mat ), length of maturity (L mat ), maximum length (L max ) and longevity (T max ) were usually point estimates. However, if a range was available, the mid-point was used. Reported estimates of fecundity were either the average egg production per year, as for many skates, or the mean number of neonates born per year based on the size and frequency of litters for live-bearing species. Natural mortality rates were estimated using Hoenig's (1983) method, an empirical approach that determines total mortality (Z) using species? maximum age (T max ) as a predictor. Data used in my analyses and references are available at http://hjort.cbl.umces.edu/elasmo.htm. The Models Several species had sufficient estimates of vital rates to justify using age- based (Leslie) projection matrix models, but this was not true of all species. Accordingly, I chose to use stage-based models for all species considered here to standardize methodologies among species. All models ran on an annual time step and 15 involved three stages: an egg/neonate stage, a juvenile stage, and an adult stage. The egg or neonate stage lasted for a duration of one time step, while the juvenile and adult stages may have lasted several years. All models were programmed in MATHCAD (v11. Mathsoft Corp. Cambridge, MA). In each model, individuals had three possible fates: they could survive and stay in the same stage, they could survive and grow into the next stage or they could die. The projection matrix took the form: (1) ? ? ? ? ? ? ? ? ? ? = 32 21 1 0 0 0 PG PG fP A Where P was the probability of surviving and remaining in the same stage, G was the probability of surviving and growing to the next stage and ? was fertility. I assumed a post-breeding census for all species; thus, fecundity had to be weighted by the probability of adult survival (Caswell, 2001). Fertilities were calculated with the following function: (2) FPf ?= where F was annual fecundity. G and P values were calculated using estimates of the probability that an individual survived (?) and the probability that an individual grew to the next stage (? i ) (Caswell, 2001). To determine the P?s and G?s to be used in the models, I assumed individuals within a stage have the same probability of survival, regardless of age. Following Caswell (2001), I iterated values of ? int in the equation: 16 (3) 1 int 1 intint ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ? i ii T i T i T i i ? ? ? ? ? ? ? until ? int equaled the value for population growth rate (?) in the eigen analysis of the projection matrix. The resulting value of ? i was used to estimate appropriate values of P i and G i . Egg/neonate mortality was calculated by assuming that every female must, on average, have one female offspring survive to ensure population persistence. This condition was empirically estimated by calculating the level of first year mortality (M 1 ) necessary to sustain a positive growth rate given values of lifetime fecundity and mortality rates (Fogarty et al. 1987). This was calculated by satisfying the constraint that: (4) ?? = = ? = = 1 11 exp1 i j M T i j mat f where M was a vector of stage-specific mortalities and all other inputs were as defined above. In several species, the vital rates of fecundity, age of maturity, and longevity did not allow for any mortality in the egg/neonate stage. This may have resulted from low estimates of longevity or fecundity. In these cases I used a value of M 1 = 0.0. My approach to estimating M 1 yielded estimates of population growth rate 17 (?) that were correlated with estimates of neonate/egg stage mortality. However, as I was interested in the relationship between vital rates and their influence on growth rate and not estimates of population growth per se, this covariation is not a critical concern. Population growth rate (?) was determined from the dominant eigen-value of matrix A (eq. 1). The intrinsic rate of population increase was calculated as r = ln(?). Elasticity analyses Elasticity was calculated as: (5) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? = ji ji ji a a e , , , ? ? where a i,j was the element in the i th row and j th column of the projection matrix and ? was the population growth rate. I developed regressions and multivariate analyses with model inputs, (T max , T mat , L max , L mat ), and elasticity values to elucidate underlying associations within the data to identify key aspects of the species life histories and phylogenic associations. Sensitivity of elasticity Sensitivity of elasticity determines the magnitude and direction of the effect of changes in individual transition elements in A on elasticity (e i,j ). Sensitivity of the elasticity for a i,j is a measure of the rate of change of elasticity to changes in 18 underlying matrix transitions. Sensitivities of e i,j provide an understanding of how life histories have shaped elasticity patterns (Caswell, 1996). For example, a positive sensitivity of the elasticity (e 3,3 ) with respect to a 3,3 (adult survival) would indicate that increasing the probability of remaining in the adult stage would increase the elasticity of the adult stage, whereas increasing the probability of death (i.e., decreasing a 3,3 ) would have a negative effect on the elasticity of the adult stage (e 3,3 ). If values of sensitivity of elasticities with respect to a i,j were negative, then reducing a i,j would increase e i,j and vice versa. The sensitivity was calculated by taking the second derivative of e i,j with respect to the element (a i,j ) in the elasticity matrix which had the greatest contribution to growth rate (Caswell, 1996; 2001), so that: (6) ji ljki lkji ji lkji ji lk ji aaa a aa a a e , ,, ,, 2 , ,, 2 , , , ? ? + ? ? ? ? ? ?? ? = ? ? ? ? ?? ?? ? ? ? where ? i,k ? j,l are Kronecker delta functions. Response to Exploitation Exploitation can change average vital rates over the short-term as a result of gear selectivity favoring particular age classes or life stages and possibly over the long-term if fishing is applied at a constant level. To understand the dynamics of elasticity for varying fishing levels, I ran models in which the fishing mortality rate ranged from 0-3.0 for juveniles and adults. Here I selected little skate, Leucoraja erinacea, common skate, Dipturus batis, dusky shark, Carcharhinus obscurus, and 19 Atlantic sharpnose shark, Rhizoprionodon terraenovae, to be representative of the range of elasmobranch life histories. Phylogenetic Analysis I used multidimensional scaling (MDS) and analysis of similarities (ANOSIM) to test for differences in suites of demographic and life history traits among taxonomic groupings. ANOSIM is a non-parametric permutation test, analogous to multivariate analysis of variance, that computes a test statistic (Global R) reflecting difference between factors (superorder, order, family), contrasted among species within each factor (Clarke & Warwick, 1994). The test was implemented using PRIMER v5 (Clarke & Gorely, 2001). Both analyses were based on a matrix of Bray-Curtis similarities of fourth root transformed and standardized data. In order to maximize the taxonomic breadth of species included in this analysis, I used the following traits: annual fecundity, T mat , T max , L max , juvenile elasticity, adult elasticity and the interstage elasticity. 20 Results Variation in elasticity with life history traits Positive relationships were found between longevity (T max ) and the elasticity of ? to changes in the juvenile and adult stages (survivals) for all elasmobranchs combined (Figures 1 & 2). A significant relationship was found between age of maturity (T mat ) and the elasticity of ? to changes in juvenile survival (Figure 3). The relationship between elasticity of ? to changes in adult survival and age of maturity was not significant (Figure 4). Longer-lived species and later-maturing species tend to have higher elasticities of juvenile and adult survival, although the rate of increase in elasticity decreases as longevity increases beyond 25 years (Figures 1-4). These associations suggest that population growth rates (?) of short-lived, early-maturing elasmobranch species are less sensitive to changes in survival during juvenile and adult stages than longer-lived, later-maturing species. Due to the structure of the matrices used, the estimated elasticities of ? to changes in the inter-stage transitions, including fertility, were equal. The elasticity of ? to changes in inter-stage transitions was negatively related to both longevity and age of maturity (Figures 5 & 6). These relationships indicate that longer-lived species have lower elasticity of ? to changes in the inter-stage transitions while short-lived and early maturing species have higher elasticity of ? values for changes in inter- stage transitions. Taken all together, these relationships (Figures 1-6) indicate a trade-off between survival and reproduction. Long-lived species may be investing more energy for survival in the juvenile and adult stages, while in short-lived species 21 there appears to be selection pressure to advance rapidly through the stages and reproduce. Sensitivity of Elasticity In 25 of 34 species of Carcharhinidae, elasticity of ? to changes in the adult stage contributed most to overall elasticity of population growth rates, while 8 showed larger contributions for the juvenile stage, and 1 for the fertility and transition stages. Of the 7 (Triakidae) species, 3 showed the greatest elasticity of ? to changes in the adult stage and 4 for the juvenile stage. Of the 9 Rajidae species, 6 showed the greatest elasticity of ? to changes in the adult stage and 3 for the juvenile stage. In total for species in the superorder Galea, 28 had the greatest elasticity of ? to changes in matrix elements for the adult stage and 12 for the juvenile stage and 1 for inter- stage transitions. For the species in the superorder Batoidea, 12 had the greatest elasticity of ? to changes in matrix elements for the adult stage and 3 for the juvenile stage. Sensitivity of elasticity was calculated for the adult stage (stage with the greatest elasticity) of species from the superorder Galea. Figures 7 and 8 show the sensitivity of elasticity values on a percent scale. Increases in survival of the adult stage would have the greatest positive effect of e 3,3 for Galea species (Figure 7). Increases in the probability of transition from juvenile to the adult stage or juvenile survival would have large negative effects of e 3,3 . In some cases, increases in the probability of transition to the juvenile stage had large negative effects on e 3,3 , but importantly, changes to fertility would have little effect. 22 Patterns in the sensitivities of e 3,3 to changes in transition elements a i,j for species from the superorder Batoidea are similar to those found for Galea, with increases in the adult stage having large positive impacts and changes in the juvenile stage and the transitional stages and fertility element having negative effects (Figure 8). For Batoidea and Galea species, the sensitivity of elasticity of species for the juvenile stage was positive, indicating that increased survival would have a positive effect on e 2,2 (Figure 9). Changes in the transition to adulthood would have large negative effects on e 2,2 . Smaller negative effects would result from perturbations of fertility, transition to juveniles and the adult stage survival. The response of elasticities to exploitation Estimates of elasticity varied little as (fishing) mortality rates were increased from low to moderate levels (0-0.4). I illustrate the general pattern by showing details for little skate, common skate, dusky shark and the Atlantic Sharpnose shark (Table 3). When considering reasonable fishing mortality rates for little skate (0-0.4), elasticity changed by 2 % or less. Similar results can be seen for the short-lived, Atlantic sharpnose shark and the long-lived common skate and dusky shark. In probable management scenarios with fishing mortality in juvenile and adult stages ranging from (0-0.4), elasticity stays relatively constant. Variation of life histories and elasticities with phylogeny A good MDS ordination was achieved with a low stress value (Figure 10). The ordination clearly suggests a continuous gradation across the life histories- 23 demography constraint space. High L max values are associated with the right of the ordination, exemplified by Lamniformes such as the white shark, Carcharodon carcharias, and the thintail thresher shark, Alopias vulpinus (Figure 11). High values of T max are associated with the top of the ordination, exemplified by the spurdog, Squalus acanthias, the common skate, Dipturus batis, and the dusky shark, Carcharhinus obscurus (Figure 11). High annual fecundities are associated with egg- laying species in the lower left of the ordination such as the lesser spotted catshark, Scyliorhinus canicula, the thornback ray, Raja clavata, and the cuckoo ray, Leucoraja naevus (Figure 11). The elasticity of ? to changes in the juvenile stage does not appear to vary systematically over the ordination. Adult elasticity appeared to be highest in the middle, whereas the elasticity of ? to changes in the interstage transitions was greatest in species at the bottom of this ordination (Figure 11). There were significant differences among superorders, orders and families in the ordination of their life histories and demography (Figure 12A; superorder - Global R = 0.45, P < 0.001; order - Global R = 0.5, P <0.001; family - Global R = 0.58, P < 0.001). Both superorders appear distinct, exhibiting only a small degree of overlap. The Lamniformes form a clear group on the right side of the order-level ordination (Figure 12B). Some of the dorso-ventrally flattened skates and rays (Myliobatiformes and Rajiformes) overlap with Carcharhiniformes in this ordination (Figure 12C). In the order-level analysis there were significant (P < 0.05) pairwise differences between Carcharhiniformes and both Lamniformes and Rajiformes. Lamniformes were significantly different from all other families. There were no significant pairwise differences between Carcharhiniformes and Myliobatiformes or between Rajiformes 24 and Myliobatiformes. There was no significant pairwise difference between Carcharhiniformes and Myliobatiformes and also Rajiformes and Myliobatiformes. There was considerable overlap in the life histories and demography of families at the centre of the ordination particularly among the shark families Carcharhinidae, Sphyrnidae, and the live-bearing ray families Myliobatidae, Dasyatidae, and Urolophidae. The families around the centre of the ordination included the egg-laying skates Rajidae and a catshark (Scyliorhinidae), the long-lived spurdog (Squalus acanthias, Squalidae), and the large pelagic predators of the mackerel (Lamnidae) and thresher sharks (Alopiidae). 25 Discussion Insights into how elasmobranch population dynamics are regulated can be gained from the combination of elasticity, perturbation and phylogenetic analyses. My analyses revealed three fundamental features of elasmobranch demography and population dynamics. Most fundamentally, I found evidence for a trade-off between survival and reproductive investment (elasticity analysis). Additionally, I observed that in the majority of species growth (survival) into adult stages appeared most important in regulating a species response to exploitation (sensitivity of elasticity). Finally, I found that life history and demographic patterns are phylogenetically constrained, such that the population dynamics and responses to exploitation of related species will be more similar than those of distantly related species (phylogenetic analysis). However, my results suggest considerable overlap of families and orders across the life history-demography range. In my models, juvenile elasticity increased (or was invariant) with increasing age at maturity and maximum age. This is in contrast to Cort?s (2002), who, using an age-based Leslie type model, found that adult elasticity tended to decline with generation time. However, the model presented here, which uses collapsed age classes as stages, and the age-structured model of Cort?s (2002) do produce similar results when age-based elasticities are added together to form each stage. Both analyses provide evidence that the trends observed in elasticity patterns apply to sharks and to the additional species of skates and rays that were included in the present analysis. 26 Elasticity analysis provided insights into what aspects of a species? life history will play important roles in understanding population level changes in response to both short-term changes in harvest policies and to longer term evolutionary pressures (de Kroon et al., 2000). I showed evidence for a trade-off between survival and reproductive investment. Generally, short-lived species had higher elasticities of ? to changes in inter-stage transitions (selection pressure on age of maturity and fertility), whereas long-lived species tended to have higher elasticities of ? to changes in adult and juvenile survival. These are relative differences, and it should be noted that the elasticity of ? to changes in inter-stage transitions is usually less than that of survival for short- and long-lived species. My findings broaden the support of a continuum of life histories for elasmobranchs (Cort?s, 2002). While this had previously been described in univariate terms (Smith et al., 2000, Cort?s, 2000, 2002), I show that the pattern is clearly multivariate and may be described as a slow?fast life history- demography continuum. In particular, there is a high degree of overlap or convergence in the life histories and demography of morphologically and phylogenetically distinct taxa such as skates and rays. The elasmobranchs studied here generally do not show high levels of variation in demographic vital rates, despite widely varying life history traits. However, reproductive mode does emerge as a life history trait that differs markedly among elasmobranchs. Live-bearing evolved 9-10 times from egg-laying ancestors and is found in 60% of elasmobranch species (Wourms, 1977; Wourms & Lombardi, 1992; Dulvy and Reynolds, 1997). In addition, elasmobranchs exhibit a broad diversity in the extent of maternal provisioning of neonates, ranging from yolk-only to uterine 27 cannibalism to complex placentation (Wourms, 1977; Wourms & Lombardi, 1992; Dulvy and Reynolds, 1997). The transition from egg-laying to live-bearing possibly reflects a trade-off occurring when the benefits of increased offspring survival exceeds the cost of reduced fecundity (Goodwin et al., 2002). Live-bearing elasmobranchs are larger than their egg-laying relatives (Goodwin et al., 2002). There are significant differences between live-bearing and egg-laying species in their life history and demography (Global R = 0.688, P <0.0001; Figure 13). However, the differences in life history and demography between an ancestral egg-layer and the skates which derived egg-laying from a live-bearing ancestor (Dulvy and Reynolds, 1997) appear to be minimal. There is considerable overlap in my ordination results between live-bearers which provide limited maternal investment, by provisioning the embryo only via the yolk of the ovum (leicithotrophy) and those which provide additional maternal input through placentation, uterine milk, or oophagy (matrotrophy). My finding of a survival-fertility trade-off and the relative importance of juvenile elasticities, particularly in larger-bodied, later-maturing species, would be consistent with the hypothesis that larger-bodied species might have evolved live-bearing to maximize offspring survival. The key question of what has driven the evolution toward large body sizes in elasmobranchs remains largely open. However, a plausible hypothesis would be that body size constrains the maximum internal volume that offspring can occupy during gestation (Qualls and Shine, 1998; Goodwin et al., 2002). The evolutionary history of elasmobranchs, particularly their reproductive modes, suggests the possibility that elasticity patterns may reflect phylogenetic 28 constraints. The multivariate analysis suggested that at the superorder and order level, patterns in elasticity and demographics in elasmobranch species broadly conformed to phylogenetic relationships. At the family level, Rajidae species formed a group separate from other elasmobranchs and were clearly distinguishable from Carcharhinidae. These findings suggest that phylogeny, demographics and population dynamics are indeed linked. Thus, similar conservation efforts may be applied successfully to closely related species. A clear separation between Rajidae and the other elasmobranch families indicated a divide between egg-layers and live-bearers, but may also indicate that skates are a particularly unusual group. Fecundity estimates were higher in Rajidae than in other elasmobranchs, while age at maturity and longevity were not significantly different (Average fecundity: Rajidae, n = 9, mean (? CL 95% ) = 27.2 ? 12.4; Other elasmobranchs, n = 47, 5.4 ? 2.36). With further research, it may be possible to determine if energetic expenses in reproductive effort differ between egg- layers and live-bearers, both at the level of individual offspring and with regard to total investment. With this caveat, my results suggest that energy per offspring is smaller in egg-layers provided total annual reproductive investment is the same for both groups (Smith & Fretwell, 1974; Einum & Fleming, 2000). This suggests that egg-layers may invest less energy to ensure early juvenile survival and more in fecundity, i.e. utilize a bet-hedging strategy (Stearns, 1992). Yet, demography of only a few egg-laying species has been studied. Clearly data and models for egg- layers, from taxa other than the Rajidae, are required to answer these questions. 29 Elasticity is robust for the range of exploitation levels likely in elasmobranchs and using a ?snapshot approach? (relatively fixed vital rates) should suffice in most management schemes unless extreme changes in survival are expected to occur in individual populations. While there were differences in the response to exploitation, elasticity of all species showed less than a 5 % change in values with exploitation levels of fishing mortality < 1.0. While elasticity provided a ?snapshot? view, the sensitivity of elasticity allowed for a more flexible view of the impact of variation in vital rates. I calculated the sensitivity of elasticity for the stage that had the greatest elasticity of ? in order to observe how that stage?s elasticity is affected by perturbation of other stages. This goes beyond the simple observation of elasticity and instead views the consequences of the dynamics in vital rates. I showed both negative and positive effects of perturbing vital rates. However, the magnitudes of change and not the direction are of importance. For all species, the transition from the juvenile to the adult stage was high, indicating the importance of attaining maturity in elasmobranchs. Contributions of fertility to the stage with the greatest elasticity of ? were low for all species. Long-lived species, often have low elasticity of ? to changes in fertility, and the perturbation analyses further indicated that the dynamics of elasmobranchs are not strongly influenced by egg/neonate production. These findings are contrary to what intuitively might be expected, and other authors have suggested that matrix models provide unreasonable elasticity of ? values for changes in fertility (Mollet and Cailliet, 2002). However, there is evidence to suggest there may be little scope or potential for varying fecundity. Firstly, demographic modeling 30 suggests that further increases in egg production would have diminishing returns for three western Atlantic skates (Frisk et al., 2002). Secondly, body cavity limitation and energetic constraints may impose phenotypic canalization on fecundity in elasmobranchs (Brander, 1981). Therefore, I suggest that compensatory responses in changes in fecundity are not likely to occur. This is consistent with a recent age- based comparative analysis (Cort?s, 2002). However, compensatory responses decreasing the age of maturity may increase lifetime egg production, while average fecundity remains constant. It appears that elasmobranch population dynamics are strongly influenced by juvenile and adult survival and the age of maturity but not fertility (Walker and Hislop, 1998; Heppell et al., 1999; Musick, 1999; Brewter-Geisz and Miller, 2000; Frisk et al., 2002). Recent research on population regulation has given supporting evidence that compensatory responses in fertility are not likely. Density-dependence in life histories has been suggested for the sharpose shark and spiny dogfish (Carlson and Baremore, 2004; Sosebee, 2004). In both cases, fecundity values were not density- dependent. Rather, the spiny dogfish showed a decrease in age of maturity and the sharpnose shark showed an increase in growth after a period of high exploitation (Carlson and Baremore this issue; Sosebee, 2004). Elasticity provides a convenient measure of life history trade-offs. Thus, trade-offs between vital rates (for example, reproduction vs. survival) may be reflected in the partitioning of elasticity. Elasticity also provides a measure of the intensity of selection pressure in each stage of a model (Caswell, 2001). I assumed that all species in this analysis can be represented with my three-stage model, which 31 will, to some extent, produce similar patterns in elasticity. A similar approach to mine was successfully used for analyzing elasticity trends and trade-offs in plant species (Silvertown et al., 1993; Silvertown et al., 1994; but see Shea et al., 1994). Caswell (2001) pointed out that the assumption of one model structure representing all species introduces bias, and this limitation should be borne in mind. The central concern centers on annualizing trade-offs among related species using elasticity from a single model structure that may or may not capture the diversity of life histories of the species. However, I feel that trade-offs can be expressed in my analyses and that the life histories of elasmobranchs do not deviate sufficiently from the basic three- stage structure to apply similar criticisms of Shea et al. (op cit.) and Caswell (op cit.). My results provide a method of prioritizing stages of a species life history that will effectively respond to management options, particularly efforts to increase juvenile and adult survival that would have the greatest impact on population protection (see Cort?s, 2002 for similar results using uncertainty in demographic models for 41 shark populations). I agree with Cort?s (2002) and Heppell et al. (1999) that caution should be used when setting management policy with elasticity analysis alone, as they do not adequately show the impact on other life stages of targeting a certain stage. The method I used to calculate elasticities here and used by Cort?s (2002) does not reflect density-dependent dynamics. However, my sensitivity of elasticity results does add some insight into potential interactions between life stages and potential stages for compensatory behavior. 32 Table 1. Data for matrix analyses: minimum and maximum fecundity or batch size = F min and F max , respectively; period between egg production in years = F period ; average fecundity = F ave (which was calculated as the mid point between minimum and maximum fecundity adjusted for yearly reproductive cycle assuming a 50:50 sex ration); longevity = T max ; age of maturity = T mat ; natural mortality = M; and egg/neonate survival = M 1 . NEP and NWA indicate Atlantic and Pacific stocks in reference to Squalus acanthias. Family Genera Species F min F max F per F ave T max T mat M M1 References Alopiidae Alopias pelagicus * 2 1 1 16 8.6 0.280.00Cort?s 2002. Alopiidae Alopias superciliosus * 2 1 1 20 12.80.220.00Cort?s 2002. Alopiidae Alopias vulpinus 2 4 1 1.5 15 5 0.300.00Cort?s 2002. Carcharhinidae Carcharhinus brachyurus 8 20 1 7 25 19.50.180.00Walter & Ebert 1991; Kailola et al. 1993; Cort?s 2000. Carcharhinidae Carcharhinus leucas 6 12 2 2.25 27 15 0.170.00 Hoenig 1979; Branstetter & Stiles; 1987; Cliff & Dudley 1991; Smith et al. 1998. Carcharhinidae Carcharhinus obscurus 3 14 2 2.12540 21 0.110.00Hoenig 1979; Compagno 1984; Smith et al. 1998. Carcharhinidae Carcharhinus brevipinna 6 10 2 2 12 7.5 0.370.00Cort?s 2002. Carcharhinidae Carcharhinus amblyrhynchos 1 6 1 1.75 18 7 0.250.00 Compagno 1984; DeCrosta et al. 1984; Smith et al. 1998; Myers 1999. Carcharhinidae Carcharhinus falciformis 10 12 1 5.5 25 9 0.180.45Branstetter 1987; Bonfil et al. 1993; Smith et al. 1998. Carcharhinidae Carcharhinus galapagensis 4 16 1 5 24 8 0.190.49 DeCrosta et al. 1984; Last & Stevens 1994; Wetherbee et al. 1996; Smith et al. 1998. Carcharhinidae Carcharhinus limbatus 1 11 2 1.5 18 7 0.250.00 Lillam and Parsons 1989; Castro 1996; Smith et al. 1998; Meyer 1999; Cort?s 2000. Carcharhinidae Carcharhinus longimanus 1 15 2 2 22 5 0.200.09 Compagno 1984; Seki et al. 1998; Smith et al. 1998; Myers 1999. Carcharhinidae Carcharhinus plumbeus 5 12 2 2.12 30 15 0.150.00 Compagno 1984; Grove & Lavenberg 1997; Smith et al. 1998. Carcharhinidae Carcharhinus porosus 2 7 1 2.25 12 6 0.370.00Lessa & Santana 1998; Cort?s 2000. Carcharhinidae Carcharhinus sorrah 3 8 1 2.75 8 2.5 0.551.01Davenport & Stevens 1988; Kailola et al. 1993. Carcharhinidae Galeocerdo cuvier 10 82 2 11.5 28 9 0.161.32De Crosta et al. 1984; Randall 1992; Smith et al. 1998. Carcharhinidae Isogomphodon oxyrhynchus 4 8 2 1.5 12 6.5 0.370.00Lessa et al. 2000. Carcharhinidae Negaprion brevirostris 4 17 2 2.62 25 12.70.180.00 Hoenig 1979; Compagno 1984; Smith et al. 1998; Feldheim et al. 2000. Carcharhinidae Prionace glauca 5 135 2 17.5 20 6 0.221.97Stevens 1975; Tanaka et al. 1990; Smith et al. 1998. Carcharhinidae Rhizoprionodon taylori 1 10 1 2.75 7 1 0.621.47Simpfendorfer 1993; Cort?s 2000. Carcharhinidae Rhizoprionodon terraenovae 1 12 1 3.25 10 4 0.440.30 Parsons 1983; Compagno 1984; Cortez 1995; Smith et al. 1998. Carcharhinidae Scoliodon laticaudus 1 14 1 3.75 6 2 0.731.32Compagno 1984 Carcharhinidae Triaenodon obesus 1 5 1 1.5 16 8 0.280.00Compagno 1984; Smith et al. 1998. Carcharhinidae Carcharhinus tilstoni 1 6 1 1.75 12 3.5 0.370.19Davenport & Stevens 1988; Cort?s 2002. Carcharhinidae Carcharhinus acronotus 3 6 1 2.25 4.5 3 0.960.00Compagno 1984; Carlson et al. 1999; Cortez 2000. Dasyatidae Dasyatis vioacea * 6 1 3 10 3 0.440.66Mollet & Caillet 2002. Lamnidae Carcharodon carcharias 8 10 2 2.25 15 9.5 0.300.00Cort?s 2002 Lamnidae L amna nasus 2 5 1 1.75 22 14 0.200.00Cort?s 2002 Lamnidae Isurus oxyrinchus 9 18 2 3.37 17 7 0.260.00Cort?s 2002 Myliobatidae Myliobatis californicus 2 12 1 3.5 24 5.5 0.190.70Martin & Cailliet 1988a; Martin & Cailliet 1988b. Odontaspididae Carcharias taurus * 2 2 0.5 17 6 0.260.00Cort?s 2002. Rajidae Raja brachyura 40 90 1 32.5 15 9 0.301.41Holden 1972. Rajidae Raja montagui 24 61 1 21.2516.3 11.40.270.61Holden et al. 1971; Holden 1972; Holden 1974. Rajidae Leucoraja naevus 71 150 1 55.2511.3 7.5 0.392.06Du Buit 1975; 1976; 1977; Walker 1998. Rajidae Dipturus batis * 40 1 20 50 11 0.092.18Du Buit 1972. 33 Rajidae Raja clavata 52 153 1 51.2523 11 0.192.19Holden 1972; Ryland & Ajayi 1984. Rajidae Leucoraja erinacea * 30 1 15 8 4 0.551.61Johnson 1979; Waring 1984. Rajidae Dipturus laevis * 47 1 23.5 50 12 0.092.25Casey & Myers 1998; Frisk et al. 2001; Frisk et al. 2002. Rajidae Amblyraja radiata * 17 1 8.5 20 5 0.221.47 Templeman 1984; Vinther 1989; Walker 1995; McEachran 2002. Rajidae Leucoraja ocellata * 35 1 17.5 20 9 0.221.30McEachran 2002; Frisk et al. 2002. Scyliorhinidae Scyliorhinus canicula 96 115 1 52.759 4 0.492.99Cort?s 2002. Sphyrnidae Sphyrna lewini 30 40 1 17.5 17 15 0.260.00Cort?s 2002. Sphyrnidae Sphyrna tiburo 3 15 1 4.5 7 2.5 0.621.51Cort?s 2002. Squalidae Squalus acanthias (NEP) 2 17 0.5 9.5 81 35 0.060.39Cort?s 2002. Squalidae Squalus acanthias (NWA)1 15 0.5 8 40 12 0.110.95Cort?s 2002. Triakidae Mustelus griseus 5 30 1 8.75 9 5.650.490.71Wang & Chen 1981; Compagno 1984. Triakidae Mustelus mustelus 2 23 2 3.12 24 13.50.190.00 Compagno 1984; Smale & Compagno 1997; Goosen & Smale 1997. Triakidae Triakis semifasciata 4 36 1 10 24 13 0.190.26Kusher et al. 1992. Triakidae Mustelus antarcticus 1 31 1 8 13 6.9 0.340.72Lenanton et al. 1990. Triakidae Mustelus californicus 2 5 1 1.75 9 2 0.490.56 Compagno 1984; Yudin & Cailliet 1990; Smith et al. 1998; Triakidae Mustelus canis 4 20 1 6 7 2.5 0.621.79 Moss 1972; Worms 1977; Francis 1981; Compagno 1984; Vooren 1992. Triakidae Galeorhinus galeus 6 52 2 7.25 40 12 0.110.85Compagno 1984; Vooren 1992; Smith et al. 1998. Triakidae Mustelus henlei 3 5 1 2 15 2.5 0.300.69 Yudin & Cailliet 1990; Francis & Francis 1992; Smith et al. 1998. Triakidae Mustelus manazo 2 19 1 5.25 10 2.250.441.66 Tanaka & Mizue 1979; Francis & Francis 1992; Yamaguchi et al 2000. Urolophidae Urolophus lobatus 2 6 1 2 14 4 0.320.06White et al. 2001. Urolophidae Urolophus paucimaculatus 2 6 0.5 4 10 3 0.440.95Edwards 1980. 34 Table 2. Relationships of vital rates of elasmobranchs species with elasticity of model parameters. Where (juv) = juvenile stage; (adult) = adult stage and (f,tr1,tr2) = fertility and the transition stages. * indicates that the relationship was not significant. Species groups Equation r 2 N F p Elasmobranchs e(juv) = 0.15?Ln(T max ) - 0.16 0.51 56 55.48 0.000 Elasmobranchs e(adult) = 0.16?Ln(T max ) + 0.17 0.38 56 33.49 0.000 Elasmobranchs e(f,tr1,tr2) = - 0.07?Ln(T max ) + 0.34 0.79 56 204.34 0.000 Elasmobranchs e(juv) = 0.17?Ln(T mat ) - 0.05 0.92 56 611.00 0.000 Elasmobranchs e(adult) = 0.02?Ln(T mat ) + 0.29 0.05 56 2.82 0.099* Elasmobranchs e(f,tr1,tr2) = - 0.06?Ln(T mat ) + 0.26 0.86 56 323.56 0.000 Requiem sharks e(juv) = 0.19?Ln(T max ) ? 0.26 0.57 34 43.34 0.000 Requiem sharks e(adult) = 0.07?Ln(T max ) + 0.13 0.30 34 13.71 0.001 Requiem sharks e(f,tr1,tr2) = - 0.09?Ln(T max ) + 0.37 0.86 34 196.89 0.000 Requiem sharks e(juv) = 0.18?Ln(T mat ) ? 0.07 0.94 34 487.07 0.000 Requiem sharks e(adult) = 0.01?Ln(T mat ) + 0.26 0.02 34 00.60 0.44* Requiem sharks e(f,tr1,tr2) = - 0.06?Ln(T mat ) + 0.26 0.87 34 208.13 0.000 Houndsharks e(juv) = 0.39?Ln(T max ) - 0.76 0.45 7 4.03 0.101 Houndsharks e(adult) = 0.00?Ln(T max ) + 0.34 0.00 7 0.00 0.99* Houndsharks e(f,tr1,tr2) = - 0.13?Ln(T max ) + 0.47 0.57 7 6.59 0.050 Houndsharks e(juv) = 0.22Ln?(T mat ) - 0.12 0.96 7 135.84 0.000 Houndsharks e(adult) = -0.06?Ln(T mat ) + 0.47 0.22 7 1.45 0.28* Houndsharks e(f,tr1,tr2) = - 0.05?Ln(T mat ) + 0.22 0.68 7 10.46 0.023 Skates e(juv) = 0.06?Ln(T max ) + 0.14 0.37 9 4.10 0.08* Skates e(adult) = 0.10?Ln(T max ) + 0.02 0.89 9 54.64 0.000 Skates e(f,tr1,tr2) = - 0.05?Ln(T max ) + 0.28 0.89 9 57.18 0.000 Skates e(juv) = 0.14?Ln(T mat ) + 0.00 0.97 9 263.73 0.000 Skates e(adult) = 0.08?Ln(T mat ) + 0.15 0.23 9 2.08 0.19* Skates e(f,tr1,tr2) = - 0.07Ln?(T mat ) + 0.28 0.74 9 20.34 0.003 35 Table 3. The results of varying exploitation levels and elasticity patterns for shark and skate species, where (juv) = juvenile stage; (adult) = adult stage and (f,tr1,tr2) = fertility and the transition stages. Little skate Common skate Dusky shark Atlantic Sharpnose shark F e(juv) E(adult) e(f,tr1,tr2) e(juv) e(adult)e(f,tr1,tr2)e(juv) e(adult) e(f,tr1,tr2)e(juv) e(adult) e(f,tr1,tr2) 0 0.19 0.24 0.19 0.34 0.40 0.09 0.42 0.42 0.05 0.19 0.28 0.18 0.2 0.19 0.23 0.19 0.34 0.39 0.09 0.42 0.42 0.05 0.19 0.27 0.18 0.4 0.19 0.22 0.20 0.34 0.38 0.09 0.42 0.42 0.05 0.19 0.26 0.19 0.6 0.18 0.22 0.20 0.34 0.38 0.09 0.41 0.41 0.06 0.18 0.25 0.19 0.8 0.18 0.21 0.20 0.34 0.37 0.10 0.42 0.42 0.05 0.18 0.24 0.19 1 0.17 0.20 0.21 0.34 0.37 0.10 0.42 0.41 0.06 0.18 0.23 0.20 1.2 0.17 0.20 0.21 0.34 0.36 0.10 0.43 0.42 0.05 0.17 0.23 0.20 1.4 0.17 0.19 0.21 0.34 0.36 0.10 0.43 0.43 0.05 0.17 0.22 0.20 1.6 0.17 0.19 0.22 0.35 0.36 0.10 0.44 0.44 0.04 0.17 0.21 0.21 1.8 0.16 0.18 0.22 0.35 0.36 0.10 0.45 0.45 0.03 0.17 0.20 0.21 2 0.16 0.18 0.22 0.36 0.37 0.09 0.45 0.45 0.03 0.16 0.20 0.21 2.2 0.16 0.17 0.22 0.37 0.37 0.09 0.46 0.46 0.02 0.16 0.19 0.22 2.4 0.15 0.17 0.23 0.36 0.37 0.09 0.48 0.48 0.02 0.16 0.19 0.22 2.6 0.15 0.16 0.23 0.35 0.35 0.10 0.48 0.48 0.01 0.16 0.18 0.22 2.8 0.15 0.16 0.23 0.40 0.40 0.07 0.49 0.49 0.01 0.16 0.18 0.22 3 0.15 0.16 0.23 0.41 0.41 0.06 0.50 0.50 0.00 0.15 0.17 0.23 36 Figure 1. The relationship between elasticity of ? to changes in the juvenile stage and longevity (T max ). The least-squares relationship is given by e(juvenile) = 0.16?Ln(T max ) - 0.17 (n = 56, r 2 = 0.51, p = 0.00). 0 0.1 0.2 0.3 0.4 0.5 0.6 0 102030405060708090 T max (years) Elasticity of the juvenile stage 37 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 102030405060708090 T max (years) Elas ticity of the adult s t age Figure 2. The relationship between elasticity of ? to changes in the adult stage and longevity (T max ). The least-squares relationship is given by e(adult) = 0.07?Ln(T max ) + 0.13 (n = 56, r 2 = 0.38, p = 0.00). 38 0 0.1 0.2 0.3 0.4 0.5 0.6 0 5 10 15 20 25 30 35 40 T mat (years) El ast i ci ty o f t h e j u v e ni le s t ag e Figure 3. The relationship between elasticity of ? to changes in the juvenile stage and age of maturity (T mat ). The least-squares relationship is given by e(juvenile) = 0.17?Ln(T mat ) - 0.05 (n = 56, r 2 = 0.92, p = 0.00). 39 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 5 10 15 20 25 30 35 40 T mat (years) Elas ticity of the adult s t age Figure 4. The relationship between elasticity of ? to changes in the adult stage and age of maturity (T mat ). The least-squares relationship is given by e(adult) = 0.02?Ln(T mat ) + 0.29 (n = 56, r 2 = 0.05, p = 0.10). 40 0 0.05 0.1 0.15 0.2 0.25 0.3 0 102030405060708090 T max (years) Elasticity of the inter-stage tr ansitions Figure 5. The relationship between elasticity of fertility and the transition between stages and longevity (T max ). Note: in a three-stage model the elasticity of ? to changes of both the transitional stages and fertility stage will be the same. The least- squares relationship is given by e(inter-stage transitions) = - 0.07?Ln(T max ) + 0.34 (n = 56, r 2 = 0.79, p = 0.00). 41 0 0.05 0.1 0.15 0.2 0.25 0.3 0 5 10 15 20 25 30 35 40 T mat (years) Ela s ticity of the inter - stage transitions Figure 6. The relationship between elasticity of ? to changes in fertility and the transition between stages and age of maturity (T mat ). The least-squares relationship is given by e(inter-stage transitions) = - 0.06?Ln(T mat ) + 0.26 (n = 56, r 2 = 0.86, p = 0.00). 42 -100% -80% -60% -40% -20% 0% 20% 40% 60% C. amblyrhynchos C. falciformis C. galapagensis C. limbatus C. longimanus C. plumbeus C. porosus C. sorrah G. cuvier I. oxyrhynchus N. brevirostris P. glauca R. taylori R. terraenovae S. laticaudus T. obesus C. tilstoni M. antarcticus M. californicus M. canis G. galeus M. henlei M. manazo S. tiburo S. canicula I. oxyrinchus A. vulpinus C. taurus Effect of perturbation of vital rate on adu l t stag e f G1 P2 G2 P3 Figure 7. The sensitivities of elasticity of the adult stage are shown for species in the superorder Galea for each stage of the matrix. The sensitivities to elasticity are shown as a percentage. Thus, the relative proportion (percentage) a stage has (a i,j ), indicates the magnitude of the effect of changes of elements (a i,j ) in A, will have on the resulting elasticities (e i,j ). Here the relative proportion (percentage) is of importance not the direction (negative or positive sign). 43 -100% -80% -60% -40% -20% 0% 20% 40% 60% D. vioacea M. californicus D. batis R. clavata L. erinacea D. laevis A. radiata L. ocellatta U. lobatus U. paucimaculatus S. acanthias (NEP) S. acanthias (NWA) Eff ect of pertu r bation of v ital rate on a dult stage f G1 P2 G2 P3 Figure 8. The sensitivities of elasticity of the adult stage are shown for species in the superorder Batoidea for each stage of the matrix. The sensitivities to elasticity are shown as a percentage. Thus, the relative proportion (percentage) a stage has (a i,j ), indicates the magnitude of the effect of changes of elements (a i,j ) in A, will have on the resulting elasticities (e i,j ). Here the relative proportion (percentage) is of importance not the direction (negative or positive sign). 44 -100% -80% -60% -40% -20% 0% 20% 40% 60% C. brachyurus C. leucas C. obscurus M. griseus M. mustelus T. semifasciata C. brevipinna S. lewini C. carchaerias L. nasus A. pelagicus A. superciliosus R. brachyura R. montagui L. naevus Effect of p e rturbatio n of vit a l rate on juvenile stage f G1 P2 G2 P3 Figure 9. The sensitivities of elasticity of the juvenile stage are shown for species of the superorders Batoidea and Galea for each stage of the matrix. The sensitivities to elasticity are shown as a percentage. Thus, the relative proportion (percentage) a stage has (a i,j ), indicates the magnitude of the effect of changes of elements (a i,j ) in A, will have on the resulting elasticities (e i,j ). Here the relative proportion (percentage) is of importance not the direction (negative or positive sign). 45 Figure 10. MDS ordination of elasmobranch life histories and elasticities. This analysis was based on annual fertility, longevity, age at maturity, maximum length, juvenile and adult survival elasticities and the interstage transitional stage elasticity. Species in the centre of the ordination have been moved outward for clarity. Car Car ambl Car Car Car Car Car leuc Car limb Car long Car obsc Car plum Car poro Car sorr Car tils Gao cuvi Iso oxyr Neg brev Pri glau Rhi tayl Rhi terr Sco lati Trd obes Scy cani Sph lewi Sph tibu Gah Mus anta Mus cali Mus cani Mus gris Mus henl Mus mana Mus must Trk semi Alo pela Alo super Alo vulp Cah Isu oxyr Lam nasu Cas taur Das viol Myl cali Amb radi Dip bati Dip laev Leu ocel Leu Leu naev Raj Raj clav Raj mont Uro loba Uro Sqa acan Sqa acan Stress 0.1 46 Figure 11. MDS ordination of elasmobranch life histories and elasticities showing the relative magnitude of each trait. L max T max Fecundity Juvenile elasticity Adult elasticity Interstage elasticity 47 Figure 12. MDS ordination of elasmobranch life histories and elasticities split by phylogenetic groupings, superorder (A), order (B) and family (C). Galea Squalea Superorder Family Carcharhinidae Scyliorhinidae Sphyrnidae Triakidae Alopiidae Lamnidae Odontaspidae Dasyatidae Myliobatidae Rajidae Urolophidae Squalidae Order Carcharhiniformes Lamniformes Myliobatiformes Rajiformes Squaliformes C B A 48 Figure 13. MDS ordination of elasmobranch life histories and elasticities overlaid with the reproductive mode of each species. Reproductive mode follows Dulvy & Reynolds (1997). Egg-laying Egg-laying Live-bearing Live-bearing Reproductive 49 Chapter 3 ORDINATION OF EVOLVED LIFE HISTORY STRATEGIES IN ELASMOBRANCH AND TELEOST FISHES: IN SEARCH OF COMMON GROUND Introduction Life history analyses can provide useful information for both managers and the understanding of evolutionary trends among related species (Beverton and Holt, 1959; Pauly, 1981; Branstetter, 1990; Hoenig and Gruber, 1990; Charnov, 1993; Cortez, 2000; Frisk et al., 2001; Dulvy and Reynolds, 2002; Williams and Shertzer, 2003; King and McFarlane, 2003; Frisk et al., 2004 -- Chapter 2). The r-K continuum has been a long-standing approach in ecology that classifies species by contrasting life history traits such as large, slow growing, low reproductive potential species (K- strategists) against small, fast growing, high reproductive output species (r-strategists) (Pianka, 1970). Although prevalent in the ecological literature, little theoretical rationale or experimental evidence supports an inherent trade-off between the intrinsic rate of population of growth, r and population carrying capacity, K (Roff, 2002). Roff (2002) suggests that one should focus on the pattern of selective pressures on life history traits rather than on traits defined from a logistic growth model. Even though the logistic equation may not represent complex life histories well, the efforts to place species in distinct groups based on similar traits highlights the utility of life history analyses, and the insight into the evolutionary dynamics of species that can be gained by such pursuits (Mertz, 1975; Caswell, 1982; Hall, 1988; Roff, 2002). 50 Winnemiller and Rose (1992) provided a classification scheme for fish taxa based on a comprehensive multivariate statistical analysis of life history traits of fishes, that did not rely strongly on a prior theoretical model. They expanded the r-K continuum to a triangular ordination of three evolved strategies in teleost fishes: (1) small, rapidly maturing, short lived species (opportunistic strategists), (2) larger, highly fecund fishes with longer life-spans (periodic strategists), and (3) fishes of intermediate size that often exhibit parental care and produce fewer but large offspring (equilibrium strategists). Elasmobranchs have long been defined as the pinnacle example of K-selected species exhibiting large size, late age of maturation, long life-span, low fecundity and high parental investment (Hoenig and Gruber, 1990; Cortes, 2004). Elasmobranchs exhibit a wide range of reproductive modes including, oviparity, vivoparity and ovophagy all of which are uncommon among teleost species. Recent analyses have placed elasmobranchs within the ?equilibrium strategy? due to their very high parental investment and reproductive strategies (Frisk et al., 2001). Yet, their large size and long life-span place them equally well within the ?periodic strategy?. Cortes (2000) performed principal components analysis on shark life history data and found alternative life history strategies within elasmobranchs. His results identified three life histories groupings in sharks: (1), species with high fecundity where offspring are born at a low percentage of adult size; (2), species that are large, slow growing and display low fecundity; and (3) species where young are born at a high percentage of adult size, have low fecundity and faster growth. Cortes (2000) related these groupings to trade-offs in juvenile mortality. He argued that species in group one 51 have higher predation rates on young compared to species in group two, while species in group three are characterized by faster growth out of vulnerable early life stages. Here, I expand the Winnemiller and Rose (1992) ordination by including elasmobranch data and thereby better representing the diversity of life history strategies among fishes. Specifically, I examine where (or if) elasmobranchs fall in the Winemiller and Rose proposed evolved triangular ordination of life history strategies both by plotting data for elasmobranch fishes directly on the original Winemiller and Rose ordination, and by conducting a new analysis that includes data for both teleost and elasmobranch fishes. 52 Methods Data Data for elasmobranchs came from the literature and from previous sources that provided compilations of life history information (Smith et al., 1998; Cortes, 2000; Frisk et al., 2001; Cortes, 2002; Frisk et al., 2004). Where necessary vital rates of Rajidae species were developed from dissections and personal observations (Frisk et al. in prep), unpublished data (Dr. N.K. Dulvy, CEFAS, Lowestoft, UK) and data from the National Marine Fisheries Service?s trawl surveys (NEFSC, 1999). Elasmobranch data included estimates from 5 orders, 11 families and 52 species. Life history data for teleost species were obtained from Winemiller and Rose (1992) and included 17 orders, 53 families and 221 species. Following Winnemiller and Rose, several life history traits were analyzed including: age of maturation (T mat ), longevity (T max ), length of maturation (L mat ), maximum size (L max ), maximum clutch size (f max ), mean clutch size (f mean ), egg size (ES), number of spawning bouts per year (SB) and parental care (PC, described below). Variables were calculated following Winemiller and Rose except for SB and PC, which are defined below. To be included, an elasmobranch species had to have estimates for L mat , f mean , ES, SB, and PC. Several variables analyzed by Winemiller and Rose either did not apply to elasmobranchs or there were insufficient data for elasmobranchs to allow their inclusion. Accordingly, the following variables, used by Winemiller and Rose, were not included in my analyses: range of egg sizes, duration of spawning season, time to hatch, larval growth rate, young of the year growth rate, adult growth rate and fractional growth rate. 53 The number of spawning bouts per year (SB) for live-bearing sharks and rays was estimated based on the frequency of reproduction. It should be noted that for some sharks this value will be less than one as many species have extended resting periods after birthing pups. The number of spawning bouts for oviparous species was calculated with the following equation: (1) b a f f SB = where f a is annual fecundity, f b is batch fecundity. In Rajidae species, the number of eggs in a batch is not well known and was assumed to be two for oviparous species (Liby, 1959; Richards et al., 1963; Johnson, 1979). Parental care was calculated by Winemiller and Rose (1992) on an ordinal scale from 0-8 with 8 being the highest possible parental investment. The level of parental investment was estimated based on several factors including the degree of placement of eggs or zygotes, protection of the eggs and zygotes, and the level of nutritional contribution to offspring (Winemiller and Rose, 1992). For example, Winemiller and Rose assigned viviparous teleost species a value of PC = 8. To include the diversity of reproductive strategies exhibited by elasmobranchs, I assigned PC = 8 for ovoviviparous (internal-eggs-live-bearing) and viviparous (live- bearing) elasmobranchs. Oviparous (egg-laying) elasmobranch species do not readily fit into Winemiller and Rose?s classification system as these species produce large eggs that undergo long gestation and developmental periods of several months (Frisk et al., 2002). Oviparous elasmobranchs were assigned PC = 6 due to their large 54 energy input per individual offspring and gestation periods. A value of 8 could be given to oviparous elasmobranch species but the lower value recognizes the lack of protection afforded to offspring compared to that of the young of live-bearing species. Egg size data for elasmobranchs was difficult to find in the literature. In order to estimate egg size I used a general empirical model to predict egg size (ES) from hatch length (HS) to fill in missing data. The following ratios, averaged, were used for elasmobranch orders: Carcharhiniformes, ES/HS = 0.06, CL 95 = 0.008, n = 8; Lamniformes, ES/HS = 0.06, CL 95 = 0.08, N = 3; Rajiformes, ES/HS = 0.21, CL 95 = 0.05, n = 5. Statistical analysis Univariate means and standard deviations of L mat , f mean , ES, SB, and PC were calculated for elasmobranch species (Table 1). All other statistics were estimated on log e -transformed data (Winemiller and Rose 1992). Two separate principal components (PCA) analyses were conducted. PCA is a multivariate statistical method that creates a linear combination of the original predictors (axes), which are not correlated. PCA assumes that the data are distributed according to a multivariate normal distribution. The principal component axes contain all the information of the original predictors (Kleinbaum et al., 1998). Associated with each axis are ?loadings? of the original life history traits that indicate the linear contribution of each life history trait to the overall principal component score. Because variables had different units, PCAs were performed using log e -transformed data and a correlation matrix as in Winemiller and Rose (1992). Data for freshwater and marine teleosts 55 were combined for all analyses as PCAs of data on subgroups were broadly equivalent. Five life history traits: parental care, size of maturity, mean fecundity, egg size and the number of spawning bouts per year explained most of the variation in 12 variable PCAs in the Winemiller and Rose analysis. These five variables were used to determine the ordination of life history strategies in teleost species of opportunistic, periodic and equilibrium strategists and for my analysis. In the first analysis, I quantified how elasmobranchs fit into the original teleost ordination. I re-analyzed Winemiller and Rose?s (1992) data for teleost fish using only 5 life history variables that were common to teleost and elasmobranch fish. Subsequently, principal component scores for each elasmobranch species were estimated by multiplying the appropriate loadings of the teleost PC axis for each trait value. Subsequently, PCA scatter plots were produced showing the distribution of elasmobranch species in the teleost-based ordination. Subsequently, I conducted a second principal component analysis involving both the teleost and the new elasmobranch data. 56 Results Descriptive statistics are shown for teleost and elasmobranch data in Table 1. Maximum size (L max ), size of maturity (L mat ), longevity (T max ) and age of maturity (T mat ) were highly correlated: (Table 2). Mean clutch size (f mean ) was also associated with maximum clutch size. A total of 151 teleost species were analyzed with a five variable ordination. The first three principal components accounted for 91% of the variation (Table 3). A review of the first principal axis associated species with large size of maturation, large mean clutch size, fewer spawning bouts, and little parental care. The second axis associated large size of maturation, large eggs, high parental care, few spawning bouts, and small clutches. The third axis grouped species with large egg size, many spawning bouts, and little parental care. Figure 1 shows the first two teleost principal components with elasmobranchs species added. The elasmobranchs all fall outside the scatter of teleost data. Elasmobranch species occur as an extreme version of the periodic strategy, closest to teleost species with a large size of maturity, large eggs, and few spawning bouts (Figure 1). The second 5-variable PCA involved a total of 203 teleost and elasmobranch species. The first three principal components accounted for 94% of the variation (Table 4). Principal component 1 reflected species with small size of maturity, small clutches, large egg size and high parental care; as opposed to those species on the second axis of large size of maturation, larger clutches, few spawning bouts and less parental care. The first two principal component axes of this ordination of the combined data are shown in Figure 2. Elasmobranchs, when included in the 57 estimation, are associated with teleost orders that have high parental care, large egg size and fall in the equilibrium strategy (Figure 2). These teleost species that form the extreme values included species with high parental investment in the order Percopsiformes (trout-perches, pirate perches and cave fish), that incubate eggs in the gill chamber of females, and Siluriformes (catfish) which are mouth breeders. Within the elasmobranchs grouping is apparent that orders seem to form separate groups. In particular, the skates (Rajidae) are distributed between the opportunistic and equilibrium strategies. 58 Discussion Superficially elasmobranchs appear to be very different from teleost species with their long life-spans, large size and reproductive strategies - reminiscent of reptiles and perhaps mammals (Frisk et al., 2001). However, there may be reason to expect a common ground between the two major classes of fishes. Elasmobranchs and teleost fishes have evolved in the same environmental conditions and face the same environmental constraints as all aquatic organisms. Elasmobranchs follow many of the same life history invariant relationships as a diversity of taxa do (including other fishes), reflected in basic trade-offs of vital rates (Charnov, 1993; Frisk et al., 2001). Certainly there is potential for a common spectrum of evolved life history traits. However, when the ordination is plotted using teleost loadings, elasmobranchs appear outside the ordination and are on the extreme of the periodic strategy. Even more surprising, they are distant from teleost species of similar reproductive strategies. The original ordination indicated elasmobranch fish should be considered periodic strategists. However, when elasmobranch data is included in the ordination they form the extreme range of equilibrium strategists and are grouped by order. The analysis of the combined dataset appears to unify life history strategies by including large species as equilibrium strategists. Because the ordination of Winemiller and Rose (1992) may not have represented the full range of reproductive modes in fish, the elasmobranchs were erroneously placed with periodic strategists. This may be a result of the fact that the data Winemiller and Rose (1992) used were biased towards species of commercial value. In general, many species representing extreme life 59 history traits are tropical or niche specialists that have received little research attention. Elasmobranchs, on the other hand, have unique reproductive strategies and have been fairly well studied. Thus, the ordination that incorporates elasmobranchs may better represent the diverse life histories of all fishes. In the ordination estimated with elasmobranchs, the oviparous species, including all Rajiformes and one Carcharhiniforme, fell in the intermediate zone between equilibrium and opportunistic strategists. This may be a result of the method I used to estimate spawning bouts and parental care for oviparous species. However, the grouping may reflect underlying biology. Rajiformes have higher fecundity compared to other elasmobranchs and are generally smaller than their live-bearing relatives (Frisk et al., 2004; Goodwin, 2002). There may be a life history trade-off reflected in the separation of oviparous species with other, larger elasmobranch species that partition more energy to individual offspring, (ES, PC), and have fewer spawning bouts (see Frisk et al., 2004; Goodwin et al., 2002). The addition of elasmobranchs to the ordination clarifies the three endpoints and suggests that Winemiller and Rose?s (1992) classification of three evolved strategists reflects fish taxa well. Furthermore, inclusion of additional data for species that are underrepresented in these data sets may further elucidate underlying life history strategies. Understanding of life histories of elasmobranchs is still a burgeoning field, but work has elucidated links between taxonomic groups and conservation of these important resources. Here I revised the equilibrium strategy and reinforced the notion of K-selection, as classification schemes from the logistic equation would suggest. 60 However, my findings and those of Winemiller and Rose (1992) suggest that fish life histories are complex and should not be defined merely by the behavior of the logistic equation alone. Although elasmobranchs are on the extreme, these findings imply that a common ordination with three endpoints represents fish life history well, including both elasmobranchs and teleost species. 61 Table 1. The mean, standard deviation, and number of observations are shown for age of maturation (year) (T mat ), length of maturation (cm) (L mat ), longevity (years) (T mat ), mean clutch size (f mean ), egg size (ES), number of spawning bouts in a year (SB) and parental care (PC). Only orders that have more than 3 observations for the five variable ordinations were included in the table. Note the oviparous, Scyliorhinus canicula, skewed the Carcharhiniformes with much higher fecundity than species in the order. T mat L mat L max T max f mean ES SB PC Carcharhiniformes mean 7.71 142.38 196.97 17.01 15.54 27.38 2.34 7.94 std 5.49 74.27 97.06 9.31 20.9 12.68 1.52 0.34 n 34.00 34.00 34.00 34.00 34.00 34.0 34.00 34.00 Lamniformes mean 8.98 303.64 446.43 17.43 5.00 12.49 0.79 8.00 std 3.39 79.51 132.35 2.64 4.50 13.54 0.26 0.00 n 7.00 7.00 7.00 7.00 7.00 7.00 7.00 7.00 Rajiformes mean 8.99 70.12 100.32 21.57 64.31 28.87 32.16 6.00 std 2.59 23.37 32.89 12.35 31.75 13.86 5.61 0.00 n 8.00 8.00 8.00 8.00 8.00 8.00 8.00 8.00 62 Table 2. Pearson correlations used in the decision to remove uninformative variables from the ordination. Age of maturity (T mat ), longevity (T max ), size of maturation (L mat ), maximum size (L max ), mean fecundity (f mean ) and maximum fecundity (f mat ) for elasmobranch data. P indicates the significance level. T mat T max L mat L max f mean f max T mat 1.00 0.84 0.31 0.23 0.08 0.06 P <0.01 0.02 0.10 0.57 0.66 T max 1.00 0.16 0.14 0.07 0.07 p 0.25 0.32 0.60 0.62 L mat 1.00 0.95 -0.19 -0.17 p <0.01 0.19 0.22 L max 1.00 -0.19 -0.16 p 0.19 0.27 f mean 1.00 0.96 p <0.01 f max 1.00 p 63 Table 3(A-C). Principal components procedure for five variable ordinations for marine and freshwater teleost species (n = 151). Tables A-C represent the correlation matrix, eigenvalues of the correlation matrix, and eigenvectors respectively. Variables included in analysis were length of maturation (L mat ), mean clutch size (f mean ), egg size (ES), number of spawning bouts in a year (SB) and parental care (PC). A. Correlation matrix L mat f mean ES SB PC L mat 1.00 0.72 0.25 -0.51 -0.37 f mean 0.72 1.00 -0.31 -0.39 -0.55 ES 0.25 -0.31 1.00 -0.18 0.28 SB -0.51 -0.39 -0.18 1.00 -0.02 PC -0.37 -0.55 0.28 -0.02 1.00 B. Eigenvalues of the correlation matrix . Eigenvalue Difference Proportion Cumulative 1 2.34 0.88 0.47 0.47 2 1.46 0.72 0.29 0.76 3 0.74 0.37 0.15 0.91 4 0.37 0.27 0.07 0.98 5 0.10 0.02 1.00 C. Eigenvectors Prin1 Prin2 Prin3 Prin4 Prin5 L mat 0.57 0.27 0.22 0.39 -0.63 f mean 0.60 -0.17 -0.13 0.39 0.66 ES -0.06 0.71 0.57 -0.08 0.40 SB -0.38 -0.45 0.59 0.55 0.01 PC -0.41 0.43 -0.51 0.62 0.02 64 Table 4. Principal components procedure for five variable ordinations for marine, fresh water and elasmobranch species (203). Tables A-C represent the correlation matrix, eigenvalues of the correlation matrix, and eigenvectors respectively. Variables included in analysis were length of maturation (L mat ), mean clutch size (f mean ), egg size (ES), number of spawning bouts in a year (SB), and parental care (PC). A. Correlation matrix L mat f mean ES SB PC L mat 1.00 0.54 -0.01 -0.49 -0.30 f mean 0.54 1.00 -0.70 -0.17 -0.79 ES -0.01 -0.70 1.00 -0.04 0.74 SB year -0.49 -0.17 -0.04 1.00 -0.04 PC -0.30 -0.79 0.74 -0.04 1.00 B. Eigenvalues of the correlation matrix . Eigenvalue Difference Proportion Cumulative 1.00 2.67 1.21 0.53 0.53 2.00 1.46 0.90 0.29 0.83 3.00 0.56 0.37 0.11 0.94 4.00 0.19 0.08 0.04 0.98 5.00 0.11 0.02 1.00 C. Eigenvectors Prin1 Prin2 Prin3 Prin4 Prin5 L mat -0.33 0.59 0.54 0.15 -0.47 f mean -0.58 0.02 0.13 0.39 0.70 ES 0.48 0.37 0.45 -0.39 0.52 SB 0.13 -0.68 0.69 0.19 -0.10 PC 0.55 0.22 -0.10 0.80 0.03 65 Figure 1. The principal components and correlations of the first two axes for teleost loadings for the five variable ordination. Elasmobranch data was added using teleost loadings. Arrows show the direction and magnitude of individual variables on the principal components where PC is parental care, ES is egg size, L mat is length of maturation, f mean is mean fecundity and SB is the number of spawning bouts per year. The three endpoints of the Winnemiller and Rose (1992) evolved life history strategies are shown, including: opportunistic strategists, that are small, rapidly maturing, short-lived species; periodic strategists, that are larger, highly fecund fishes with longer life-spans; and equilibrium strategists, fishes of intermediate size that often exhibit parental care and produce fewer but large offspring. -5 -4 -3 -2 -1 0 1 2 3 4 5 -5 -4 -3 -2 -1 0 1 2 3 4 5 Principle component 1 Principle component 2 Elasmobranchs Acipenseriformes Anguilliformes Atheriniformes Clupeiformes Cypriniformes Cyprinodontiformes Hexanchiformes Gadiformes Gasterosteiformes Perciformes Percopsiformes Pleuronectiformes Salmoniformes Scorpaeniformes Scomberiformes Siluriformes Osteoglossiformes Equilibrium strategy Opportunistic strategy Periodic strategy L mat PC ES SB f mean 66 -5 -4 -3 -2 -1 0 1 2 3 4 5 -5 -4 -3 -2 -1 0 1 2 3 4 5 Principle component 1 P r i n ci p l e com p o n en t 2 Carcharhiniformes Hexanchiformes Lamniformes Rajiformes Squaliformes Acipenseriformes Anguilliformes Atheriniformes Clupeiformes Cypriniformes Cyprinodontiformes Gadiformes Gasterosteiformes Perciformes Percopsiformes Pleuronectiformes Salmoniformes Scorpaeniformes Scomberiformes Siluriformes Osteoglossiformes Equilibrium strategy Opportunistic strategy Periodic strategy L mat PC ES SB f mean Figure 2. The principal components and correlations of the first two axes for teleost and elasmobranch loadings for the five variable ordination. Arrows show the direction and magnitude of individual variables on the principal components where PC is parental care, ES is egg size, L mat is length of maturation, f mean is mean fecundity, and SB is number of spawning bouts. The endpoints to the Winnemiller and Rose (1992) ordination of evolved life strategies are also shown. 67 Chapter 4 AGE, GROWTH, AND LATITUDINAL PATTERNS OF TWO RAJIDAE SPECIES IN THE NORTHWESTERN ATLANTIC: LITTLE SKATE, LEUCORAJA ERINACEA, AND WINTER SKATE, LEUCORAJA OCELLATA Introduction Adaptation and evolution work in a framework of trade-offs where vital rates are selected to maximize individual fitness. In many species there is a latitudinal gradient of life histories where local populations adapt and evolve to maximize fitness for specific environmental conditions. For instance the rattlesnake, Crotalus viridis, exhibits a latitudinal gradient in both age at maturity and mortality rates (Shine and Charnov, 1992). However, life history gradients and the bioenergetics that underlie them are not always straightforeword. Conover has proposed the ?counter-gradient? hypothesis wherein the potential for growth varies inversely with latitude (Conover, 1990). Based on his work on silversides (Menidia menidia), Conover suggested that the shorter growth season and temperature-induced overwinter mortality at high latitudes selects for faster growth rates. However, many marine fish species have a latitudinal gradient where populations in higher latitudes exhibit slower growth, later age at maturity, increased longevity and likely lower population productivity (Taylor, 1958; Beverton and Holt, 1959; Beverton, 1992). Little skate and winter skate are both in the genus Leucoraja; yet, they have diverging life histories and have the potential to exhibit latitudinal gradients. Little skate L. erinacea is a small species reaching a size of 57 cm (total length, TL) with moderate growth rates, while winter skate L. ocellata is larger, reaching 111 cm TL, and is slower growing (Richards, 1963; McEachran and Musick, 1973; Johnson, 68 1979; Waring, 1984; NEFSC, 1999; Sulikowski et al., 2003; Alvarado Bremer et al., 2004). Little skate and winter skate occur from North Carolina, U.S. to Canadian waters with the centers of abundance in the southern New England and Georges Bank regions; additionally little skate is abundant in the mid-Atlantic (NEFSC, 1999; McEachran, 2002). Both species are common in the Gulf of Maine and off the Canadian coast but are less abundant in these northern portions of their range (McEachran, 2002). Both species are ecologically important and constitute a large proportion of the demersal biomass in the western Atlantic (Link et al., 2002). Skates likely play an important role in the trophic dynamics in western Atlantic ecosystems as a predator on benthic fauna including amphipods, crabs, small fish, shrimp, mollusks and worms (Nelson, 1993; Murdy at al., 1996). Although traditionally not targeted by the domestic fleet, skates were targeted by foreign fleets prior to implementation of the 200-mile limit (Fogarty and Murawski, 1998). Recently, skates have become more valuable as a commercial species and are being landed and removed from the demersal community at greater rates (NEFSC, 1999; McEachran, 2002). McEachan (2002) reported that US landings grew in the Gulf of Maine from 297 metric tones (mt) in 1981 to 15, 000 mt in 1996. Similar trends have been observed on Georges Bank and southern New England (NEFSC, 1999). Skates are sought for their mild tasting meat, characterized by low fat and cholesterol content, and a taste similar to shellfish (McEachran 2002). Skate wings are sold to the European market and increasingly the US domestic market (McEachran, 2002). Current estimates of fishing mortality have indicated that winter skate may be overfished (NEFSC, 1999). 69 Knowledge of skate life histories in the western Atlantic is increasing. However, the lack of data for some species has made the formulation of management policies difficult. While there is limited information about the vital rates of western Atlantics skates compared to teteost fishes, there has been progress in aging many elasmobranch species including several skates (Caillet and Goldman, 2004). Seven Rajidae species have had annual band formation verified with several techniques including lab-based tetracycline tagging studies, marginal increment, back calculation and size frequency analyses (Caillet and Goldman, 2004; Natanson, 1993). Previous aging studies have verified annuli band formation in both little skate and winter skate (Natanson 1993; Sulikowski et al., 2003). Natanson (1993) conducted laboratory experiments under varying temperature conditions and verified that annual bands were formed in little skate. Her experimental evidence suggests that temperature did not have an impact on band formation (Natanson, 1993). Natanson?s (1993) results indicated that previous estimates of age (Richards 1963; Johnson 1979; Waring 1984) were conservative in their interpretation of annuli and may have underestimated ages in little skate. Sulikowski et al. (2003) estimated the age of winter skate in the Gulf of Maine and performed a successful marginal increment analysis to verify annual band formation. In the Gulf of Maine winter skate?s maximum marginal increment occurred in May and the minimal increment occurred in July indicating that band formation likely occurs between June and July (Sulikowski et al., 2003). To better understand the biology and ecology of little skate and winter skate, estimates of vital rates are needed over the species ranges. In this project I collected little skate and winter skate samples over the species range to estimate region-specific 70 age and growth estimates. Specifically, the goals of the study are to perform regional analyses of vital rates covering the ranges of little skate and winter skate from Cape Hatteras to Canadian waters by estimating age, growth and longevity of each species. 71 Methods Sampling Samples were collected during the National Marine Fisheries Service?s (NMFS) annual fall (September and October: 2001), winter (February: 2001), spring (February, March and April: 2001) and summer (June, July and August: 2001) surveys from the NOAA RV Albatross IV. Skates were collected along the Atlantic coast from Cape Hatteras to Canadian waters (Figures 1 & 2). The fall, spring and winter surveys used a bottom trawl equipped with 1.27 cm mesh liner that was towed for 30 minutes at three knots. Additionally, the winter survey is equipped with a flat net and sweep chain. The summer survey used a standard 8-foot New Bedford scallop dredge with a 2-inch ring chain bag and 1.5 inch mesh webbing that was towed for 15 minutes at 3.8 knots. Total length, disk width, weight and sex of individual skates were recorded aboard the ship. Latitude and longitude was recorded for each tow and used to provide location information for specimens. Small specimens (< 20 cm, TL) were frozen whole at sea. Vertebral and tissue samples of larger specimens were removed and frozen at sea. All samples were then shipped to the Chesapeake Biological Laboratory for further analyses. Additional little skate samples were collected on 12-12-00, 3-29-01 and 5-25-01 aboard the F/V Tony & Jane, a 57 foot scalloper located in Ocean City, MD, captained by Mr. J. Eustler. Skates > 30 cm in total length were identified morphologically. However, little skate and winter skate are very difficult to differentiate morphometrically below 25-30 cm TL (McEachran and Musick, 1973). To differentiate between these species I used a rapid PCR-RFLP assay based on the restriction endonuclease Sty I that generates 72 fragments that are easily characterized through agarose gel electrophoresis (Alvarado Bremer et al., 2004). All skates below 30 cm TL were genetically identified prior to further analyses. In total, 675 male and 1034 female winter skate and 1222 male and 1118 female little skate specimens were collected for age, growth and reproductive analyses (reproductive results reported in Chapter 5). The higher number of females likely does not reflect a biased population sex ratio, but rather the sampling protocol that, in addition to vertebrae, also collected ovaries. Limited numbers of mature female (n = 168) and male (n = 172) winter skate (TL > 76 cm) were caught. The sampling covered the area from Cape Hatteras to the upper regions of the Gulf of Maine. For region-specific analyses, the areas were divided into three regions based on potential physical features that may separate stocks. The following three regions were defined: (1) mid-Atlantic; representing areas north of Cape Hatteras to the Hudson River canyon; (2) southern New England-Georges Bank; north of the Hudson River canyon to the outer edges of Georges Bank; and (3) the Gulf of Maine; from the Northeast channel southwest to Cape Cod. When performing analyses of all regions combined, I will refer to the entire region as the northeast coast (N.E. coast). Specimens for age analyses were selected to ensure coverage of species, size, the geographic regions, and inshore and offshore areas. Length-weight analyses The relationship between individual length and weight was analyzed for the mid-Atlantic, southern New England-Georges Bank and the Gulf of Maine. As little 73 skate had ample data and comparisons were made for each region based on the following area designations designed to provide greater separation between regions. For little skate the mid-Atlantic was defined as below the latitude of 38.72, southern New England-Georges Bank was defined as west of longitude 72.0 and south of latitude 41.7 and the Gulf of Maine was defined as the area above latitude 42.45. Additionally, separate models were estimated for male and female little skate. Linear regression models were fit using log 10 -transformed data for length and weight (Proc REG, SAS Corp.). Analysis of covariance was performed in Proc GLM (SAS, Corp.) to test for significant regional and sex-specific trends. The three regions were compared for significant differences by employing orthogonal contrasts. Slide preparation and age determination Vertebral samples were cleaned of adhering tissue and the 4 th and 8 th vertebrae were removed. To aid in the cleaning process vertebral sections were soaked in warm water for several minutes until tissue was softened. Warm water treatment did not interfere with band appearance or composition (James Gelsleichter, Mote Marine Laboratory, Fl, personal communication, 2003). If the first complete vertebral section could not be identified, one of the first ten was removed. Vertebral centra were then air dried, mounted and cut sagittally using a low speed Buehler isomet saw into thin (1-2 mm) sections, glued to a glass slide, wet sanded with graded sand paper (250-1000 ?m) and stained with Toluidine Blue. For small specimens (< 20 cm), the centra were imbedded in Struers epoxy prior to being sectioned. 74 Several techniques were tried to enhance the readability of slides including no treatment, immersion of centra in RDO, a rapid decalcifying agent, prior to drying and sectioning (Natanson, 1993). Additionally, several stains were used post- sectioning including Rose Bengal and Toluidine Blue. I had limited success polishing slides with 0.3 Micron Alpha Alumina powder as contrast was often reduced. The best results were obtained with emersion in 95% ethanol for at least 12- 24 hours (Richards et al., 1963; Waring, 1984; Zeiner and Wolf, 1993) and staining with Toluidine Blue. Ageing little skate and winter skate I followed Natanson?s (1993) interpretation of annual growth bands for little skate, and Sulikowski et al.?s (2003) interpretation for winter skate. Annuli are easily interpreted for winter skate up to approximately 10 years depending on the specimen. In some samples false rings or check marks were difficult to distinguish from true annuli, especially for old individuals. For each specimen only one vertebral section was read. To estimate age and growth relationships, each centrum section was read by two readers (R1 & R2) without knowledge of prior age estimates or skate length. The following criteria were used in reading little skate. For most little skate sections, R1 read each centrum twice and R2 once. If one of the two R1 reads were within 2 years of R2?s read, then the sample was kept and the average of the two reads in agreement was taken. In some cases the average of the two reads R1 made equaled R2?s single read; in this circumstance the average of R1 was used. In some cases only R2 read a centrum. 75 When only R2 read a centrum, is was read twice and if the reads were within 2 years the average was taken. Two approaches exist in estimating the age of winter skate: (1) ?liberal read?, count everything that can be interpreted as an annuli band; (2) ?conservative read?, count only the most obvious bands as annuli. The liberal read was adopted as my standard ageing technique. This decision was based on the lesson of Natanson?s (1993) validation study for little skate and a ?precautionary? management perspective. To ensure that the choice of the ?liberal? standard did not unduely bias results, the two standards were compared. Winter skate specimens were read once assuming a ?liberal? criterion by R1 and R2, except for a few that were read by only R2. If R1 agreed within 3 years of R2 or if R2 read twice and the reads agreed within three years, the specimen was kept for further analyses. Winter skate conservative reads were read by R2 twice and if they did not agree within three years, the specimen was not included in growth estimates. Subsequently, the resulting growth curves estimated by ?liberal? and ?conservative? reads were compared. In addition, the readability or quality of all winter skate slides were ranked from 0-6, with six being very high readability and zero indicating very low readability. High readability scores indicate that problems encountered in identification of annuli were minimal. Comparisons were made between growth curves based on different readability rankings. 76 Precision and reader bias Reader precision was estimated on raw data with no culling of specimens. However, a few samples were damaged in the cutting and mounting process and could not be read. Average percent error (APE) was used to estimate precision within and between readers: ? = ? ?= R i j jij j X XX R APE 1 1 %100 where R is the number of reads, X ij is the i th reading of the j th fish and j X is the mean of readings of the j th fish (Campana, 2001). Precision estimates were made between and within readers R1 and R2. In addition, a third reader (R3) read a sub-sample of centrum and the reads were only used for precision estimates. Reader bias was estimated using mixed model analysis of variance in which readers were treated as random variables (Proc Mixed, SAS Corp.). I tested the null hypothesis of no significant reader effect on estimated mean ages. Centrum and growth analysis A von Bertalanffy growth function was fit to estimated age and length data at capture (von Bertalanffy, 1957). The function took the following form: )1( )( 0 ttk t eLL ?? ? ?= where L t is length at age t, L ? is the asymptotic length, t is age, t 0 is the age at zero length and k is the growth coefficient. I used the average of all reads that met the 77 predetermined standards (above) as the basis for age in my growth models. All growth models were fit in Proc Nlin, SAS. Because of a lack of small samples in some model fits for winter skate, I assumed a hatching size of 16 cm and set the theoretical size at zero (t 0 ) accordingly. In order to test if growth curves differ significantly by region, the non-linear comparison method residual sum of squares of Chen (1992) was employed. Separate growth models were fit to data from each region and for a pooled for all regions combined. The test of significance to indicate that the individual models are a better description of the data than the pooled model is given by: i i iP iP DF RSS DFDF RSSRSS F ? ? ?? ?? = where RSS is the sum of squares, DF is degrees of freedom and individual and pooled values are indicated by the i and p respectively. The procedure is performed by estimating an F statistic with degrees of freedom of 3?(C-1) and comparing it to an F with (N-3?C) degrees freedom (Chen, 1992). C is the number of curves being compared and N is the total or pooled sample size. If significant differences were found between curves, Kimura?s (1980) likelihood ratio test was used to test for significant differences in model parameters. 78 Results Species identification PCR-RFLP analyses indicated that little skate dominated catches of skate < 30 cm TL. All skates less than 15 cm (TL) (n = 55) were little skate, 25 out of 29 skates 16 to 20 cm (TL) were little skate and 66 of 71 skates 21-39 cm (TL) were little skate. Winter skate were rare among collection of juvenile-sized skates and absent below 16 cm TL. All juvenile winter skate less than 30 cm TL were caught on Georges Bank except one 29 cm (TL) winter skate caught in the mid-Atlantic region. In contrast, juvenile little skate were caught over the range of the species. Length-weight relationship Derived relationships and figures are based on species identification by morphological (> 30 cm) and genetic criteria (< 30 cm). The maximum sizes of female and male little skate were broadly similar (Figure 3. Table 1). Region-specific length-weight relationships were significant for little skate (Table 1). Finally, when compared between regions, there were significant differences between the Gulf of Maine and Southern New England-Georges Bank, and the mid-Atlantic; however, there was not a significant difference between the mid-Atlantic and southern New England-Georges Bank (Table 1). There was a significant allometric relationship between weight and length for winter skate for the N.E. coast (Figure 4: Table 1). Females were heavier at size than males (Table 1). However, male winter skate grew to substantially larger and heavier sizes (Figure 4). Winter skate also exhibited a significant length-weight relationship 79 for each region (Table 1). However, this relationship did not differ among regions (Table 1). Image quality and annuli identification Sample images of centra of little skate and winter skate for old and young skates are shown in Figure 5. Annual band interpretation in both species was adequate to good over the geographic range of the species. In little skate the percent of specimens readability rankings were (0-6) as follows respectively: 0.87, 18.8, 37.0, 35.0, 8.3, 1.0. For winter skate, percent of specimens readability rankings (0-6) were as follows respectively: 0.7, 8.4, 24.1, 33.6, 17.5, 13.5, and 0. Ageing little skate and winter skate The largest female little skate aged was 56 cm and 12 years old, while the largest male aged was 57 cm and 12 years old. The oldest aged little skate was a female T max = 12.5 years and 46 cm TL. The oldest individual in the mid-Atlantic was estimated to be 11 and all other individuals 11 years old or older came from southern-New England-Georges Bank and the Gulf of Maine regions. The largest female winter skate aged was 90 cm (TL) and largest male aged was 107 cm (TL). The oldest female winter skate was 76 cm TL and was estimated to be 19.5 years, and the oldest male was 74 cm TL skate that was estimated as 20.5 years old. I note in passing that a male 88 cm (TL) was aged as old as 24 but estimates were not in agreement between readers (R1 age = 24, R2 age = 20). 80 Precision and reader bias Ninty-nine percent of little skate age estimates were within two years between readers R1 and R2. The little skate APE precision estimate across readers was 20.1% (n = 137, R1, 2-reads, R2, 1-read) and within reader precision was 11.6% (n = 35, R2, 3-reads). No significant differences were found in reader bias for little skate (reader- fish interaction: F 1,101 = 1.28, p = 0.10). Eighty percent of R1 and R2 reads of winter skate specimens were aged within 3 years. The winter skate APE precision estimate across readers was 18% (n = 40, R1, 1-read; R2, 1-read; R3, 1-read) and 8.0% (n = 46) within reader (R3, 3-reads). There was no significant reader bias between reader R1 and R2 (reader-fish interaction: F 1,106 = 0.86, p = 0.71). Centrum and growth analyses Little skate growth models were estimated for the N.E. coast, mid-Atlantic, southern New England-Georges Bank and the Gulf of Maine (Figure 6, Table 2). Parameters of the von Bertalanffy model were highly correlated (Table 3). Since no age-0 little skate were caught in the Gulf of Maine, growth curves were fit with a t 0 fixed to correspond to the hatch size estimated for the N.E. coast (11.2 cm TL). Growth rates (k) declined and asymptotic length (L ? ) increased with latitude for little skate (Figure 7). Latitudinal patterns are apparent and pairwise comparisons between the Gulf of Maine and the Mid-Atlantic were significant based on the Chen et al. (1992) residual sum of squared test (Gulf of Maine vs. mid-Atlantic: F 3,143 = 2.87, p = 0.04; Gulf of Maine vs. southern New England-Georges Bank, F 3,137 = 2.30, p = 0.08; 81 southern New England vs. mid-Atlantic: F 3,181 = 1.88, p = 0.90). Based on Kimura?s likelihood ratio test, there were no significant differences between parameters of models for the Gulf and Maine and the Mid-Atlantic (k, p = 0.94; l ? , p = 0.28; t 0, p = 0.43). Comparison of growth curves by Chen et al.?s (1992) method may have estimated significant differences because of very large sample size, fewer estimates of small and large individuals and may not reflect biological differences. Male little skate appear to grow slower and reach a larger size than females (Figure 8, Table 1). However, there were no significant differences between sex- specific curves (residual sum of squares: male vs. female; F 3,181 = 1.78, p = 0.15). Comparison of the results presented here and those of previous studies on little skate are shown in Figure 9. Work conducted in Long and Block Island Sound and the southern New England and Georges Bank regions all show higher growth rates and shorter longevity (Richards et al., 1963; Johnson, 1979; Waring, 1984). However, more recent work, which includes lab validation of growth bands, found ages and growth rates similar to the present study (Natanson, 1993). Winter skate growth was estimated for the N.E. coast (Figure 10; Table 2). Parameter estimates of the von Bertalanffy model were highly correlated (Table 3). Sample coverage was not adequate for separate growth models for specific regions. Males were larger at age than females (Figure 11). Significant differences were found between the growth models of male and female winter skate (residual sum of squares: Male vs. females; F 3,198 = 2.78, p = 0.042). However, further analyses indicated that individual parameters were not significantly different (k, p = 0.46; l ? , p = 0.72; t 0, p = 0.43). Comparison of growth curves by Chen et al.?s (1992) method 82 may have estimated significant differences because of very large sample size, fewer estimates of small and large individuals and may not reflect biological differences. von Bertalanffy estimates for winter skate assuming ?conservative? annuli interpretation provided nearly identical results to estimates from ?liberal? reads (Figure 12, Table 2). Maximum observed conservative age (T max = 15) was lower than the liberal read (T max = 20.5) while the sizes at age were similar. The growth curve for the ?conservative? fit does not flatten out as much as the model for the ?liberal? read. Differences between the conservative fit and the liberal read were significant (residual sum of squares: conservative vs. liberal; F 3,462 = 22.4, p = 0.0001). Additional comparisons using Kimura?s likelihood ratio tests found no significant differences in parameter values (k, p = 0.81; l ? , p = 0.99; t 0 , p = 0.93). Growth curves of specimens ranked as ?high readability? (H.R. = 4-6) were estimated for winter skate and there was significant difference between higher readability and all reads (Figures 12, Table 2; residual sum of square results: H.R. vs. all reads: F 3,301 = 3.65, p = 0.013). Additional comparisons using Kimura?s likelihood ratio tests found no significant differences in parameter values (k, p = 0.054; l ? , p = 0.30; t 0 , p = 0.48). von Bertalanffy model parameters estimated were similar between curve fits for the N.E. coast, high readability (4-6) and all data (0-6). All three interpretations of annuli provide very similar growth and asymptotic size estimates. 83 Discussion I found evidence for a strong latitudinal trend in maximum size in little skate. Little skate from northern regions were larger but grew more slowly than little skate from more southern regions. No similar trend was found in winter skate. Richards et al. (1963) reported that total length in little skate increased with increasing latitude. Their results indicated that little skate from the Miramichi estuary, New Brunswick, averaged 50 cm (TL), while little skate in southern New England and Delaware Bay averaged 47 cm (TL) and 43 cm (TL) respectively (Richards et al., 1963). Little skate follows a latitudinal gradient of increased size, longevity and decreased growth with increasing latitude (Taylor, 1958; Beverton and Holt, 1959; Beverton, 1992). Other elasmobranchs have shown similar variation in vital rates with varying localities and latitudes. Parson (1993) compared populations of the bennethead shark, Sphirna tiburo, from temperate and tropical regions and found that individuals in the temperate population grew slower and reached a larger size. Yamaguchi et al. (1998), studying the starspotted dogfish, Mustelus manazo, reported significant differences in vital rates between regions in the Japan Sea and Pacific Ocean. While starspotted dogfish populations with the greatest temperature differences showed the greatest difference in growth rate, the intermediate sites showed evidence against a temperature effect. Large size, longevity, and slow growth are common indicators that a species is vulnerable to overexploitation and of low population productivity (Frisk et al., 2001; 2002). Within elasmobranchs, low productivity has been associated with population declines. For example, the more productive gummy shark appeared resilient while 84 the less productive school shark experienced a population decline in the same fishery (Stevens, 1999). I suggest that as little skate populations likely have slower growth, increased longevity and lower productivity in higher latitudes, these populations are likely more susceptible to declines under exploitation. Winter skate specimens were relatively few in comparison to little skate and I was unable to conduct regional growth analyses. A comparison of the growth curves estimated for winter skate in this study to those of Sulikowski et al. (2003) provide some evidence that winter skate may grow slower and attain a larger size in higher latitudes. Additionally, McEachran and Martin (1977) reported that larger individuals of winter skate are more common in higher latitudes. However, Simon and Frank (1996) provide conflicting evidence. They report a growth rate of k = 0.14 yr -1 and asomototic size l ? = 114 cm of for winter skate off of Nova Scotia (Simon and Frank, 1996; Sulikowski et al., 2003). Simon and Frank?s (1996) growth estimates are considerably higher than the present study and that of Sulikowski?s et al. (2003). However, Sulikowski et al. (2003) cites personnel communication with Simon that the oldest individuals were likely underaged by four or more years. Presently, there is little evidence to support a latitudinal gradient of variation in vital rates in winter skate. Although the trends are not clear, a review of maps of juvenile and adult winter skate regional abundances reveal that adults and juveniles are fairly mobile (NMFS, 2003). Adult and juvenile winter skate occur, or gather, in great numbers on Georges Bank, while juveniles appear to occur in the mid-Atlantic region at disproportionately higher abundances than adults (NMFS, 2003). Further, all genetically-identified, newly hatched winter skate were caught on Georges Bank. 85 In contrast, recently hatched specimens of little skate were caught in all regions in the N.E. coast. A possible explanation of these patterns is that Georges Bank is a spawning ground for winter skate. However, sample size was not large and the abundance maps are not conclusive; thus these results are speculative in nature. von Bertalanffy growth estimates indicate that male little skate attain a larger size than females. This is in contrast to winter skate and and batoid species including Raja binoculata, R. brachyuran, R. clavata, R. miraletus, R. montagui, and R. rhina, where the female attained a greater length (Zeiner and Wolf, 1993; Abdel-Aziz, 1992; Holden, 1972; Brander and Palmer, 1985; Waring, 1984; Ryland and Ajayi, 1982; Walmsley-Hart et al., 1999). Richards et al. (1963) also found that male little skate attain a larger size. In fact Richards et al. (1963) noted that females greater than 50 cm (TL) were rare in Long and Block Island sounds while larger males were common. Length-weight relationships indicate that female little skate and winter skate are heavier at length than males. However, male winter skate reach a larger size (TL) and weight. Sosebee (2004) estimated physiological maturity of several skate species in the western Atlantic and in her data set the largest individual little skate and winter skate were males. Although Sulikowski et al.?s (2003) growth estimates for the Gulf of Maine predicted a larger L ? for females (Female L ? = 137 cm; Male L ? = 122 cm), they were unable to collect a large number of adult winter skate and may not have adequately sampled old/large winter skate. Two hypotheses for why male winter skates are larger can be developed. First, fishing gear may select heavier, at length, females at a higher rate than males. 86 Thus, allowing males to live longer than females. Second, sexual selection may favor larger males. Elasmobranchs have internal fertilization, and larger males may be able to control and impregnate females more successfully than smaller males in winter skate and little skate. My findings, like those of Natanson (1993), indicate that the observed maximum age for little skate is older than previous estimates of T max = 5, (Johnson 1979), and T max = 8 years (Richards et al., 1963; Waring, 1984). The oldest female and male little skate were 12.5 and 12 years respectively. My growth estimate for little skate was slower than previous work but is in closer agreement with Natanson?s ageing work. Additionally, previous growth curves did not exhibit asymptote behaviour suggesting that larger, older individuals were not adequately sampled (Richards et al., 1963; Waring, 1984). Assuming that longevity can be estimated as the fraction (95%) of L ? attained in a life-span (Taylor, 1958), the theoretical life- span (T max ) in little skate would equal: 13 in the mid-Atlantic, 14 in southern New England-Georges Bank and 15 in the Gulf of Maine. Previous work conducted on winter skate in the Gulf of Maine indicates that the species may grow slightly slower and reach a larger size than my data suggests (Sulikowski et al., 2003). The slightly lower growth of Sulikowski?s et al. (2003) work may explain the higher L ? (i.e. k and L ? correlation). Additionally, the largest individual that Sulikowski et al. (2003) aged was TL = 94 cm. The present analysis included 13 winter skate that were 94 cm (TL) or larger and all samples above 90 cm (TL) were male. The largest individual caught on the National Marine Fisheries Service?s annual surveys (1963-present) was TL = 111 cm (TL). This is in agreement 87 with the size of L ? for southern New England-Georges Bank, high readability and the conservative growth models. It is possible that my data set better represented the older/larger age classes of winter skate than previous work. The oldest aged female and male winter skate were 19.5 and 20.5 respectively. Previous work aged winter skate as old of 16 and 19 years (Simon and Frank, 1996; Sulikowski et al., 2003). While the juvenile winter skate (0-10 years) were relatively easy to read and determine annuli, adult winter skate were harder to age. Difficulty ageing older fish may have produced the growth curves that do not reach fully an asymptote. However, it is relatively common in Rajidae species to see growth curves not flatten out (Sulikowski et al., 2003). Until age of larger individuals can be verified, estimates of longevity based on aged individuals should be interpreted with caution and may be underestimates. Assuming that longevity can be estimated as the fraction (95%) of L ? attained in a life-span, and the uncertainty discussed above, theoretical life-span (T max ) in winter skate is 35 years (Taylor, 1958). In some adult winter skate, the first 8-10 annuli were clear while false rings or check marks became more common in subsequent annuli. One potential confounding factor is if check marks are formed as a result of mating, which has been proposed in the viviparous Australian sharpnose shark, Rhizoprionodon taylori (Simfendorfer, 1993). Defining spawning time in a serial spawner would likely be a difficult process. A second ageing problem that occurred in some individuals is that the distance between annuli increased towards the outer potions of the centra; indicating that I may be undercounting annuli. I ultimately used the liberal read, in which I 88 counted every mark that could be annuli, as the standard. In this method, readers still avoided counting obvious check marks that were not annuli. However, in some cases the task of distinqishing between check marks and annuli remained difficult. Ultimately, the conservative read provided similar growth rates, but a lower maximum age (15, T max ), and the high readability fit was very similar to the liberal fit. However, the life history of winter skate suggests that the conservative read is not accurate. With age at maturity estimates ranging 10-14, (Chapter 5), suggests longevity is greater than 15 years. The high readability fit, specimens with confidence in annuli interpretation, provided nearly identical maximum age and growth rates as the liberal read. However, the high readability fit lacked many of the largest specimens and may have excluded many of the oldest specimens. Here again underestimation of adult ages may be occurring. Validation of old specimens for winter skate is needed to better understand maximum age and growth in older fish. Until then, a precautionary management philosophy and the lessons of Natanson (1993) work on little skate, calls for a liberal interpretation of annuli. Precision between readers was relatively high at 20 and 18 percent for little skate and winter skate. Within reader precision was lower at 11.6 and 8.0 percent for little skate and winter skate respectively. Precision estimates in elasmobranchs tend to be lower than teleost fish and are often not reported (Campana, 2001; Cailliet & Goldman, 2004). Some little skate specimens that were 19-20 cm (TL) were estimated to be age 0. It was unclear if this was an ageing error or a result of date of birth. Little skate is a continuous spawner and hatching can occur at any time of the year while annuli formation is a seasonal process (Natanson, 1993; Johnson, 1979). 89 It has long been known that elasmobranchs, especially large, late maturing, long-lived and slow growing species, are very susceptible to overexploitation and even local extinction (Holden, 1973; Brander, 1981; Casey and Myers, 1998; Walker and Hislop, 1998; Dulvy et al., 2000; Stevens et al., 2000). Further, species in the family Rajidae have been ranked as one of the marine fishes most susceptible to overexploitation (Dulvy et al., 2000). Previous analyses have indicated that moderate to low exploitation rates could lead to the decline of little skate, Leucoraja erinacea, and winter skate, Leucoraja ocellata, populations in the western Atlantic (Frisk et al., 2002). While skates are likely candidates for overexploitation, a lack of knowledge of vital rates of little skate and winter skate have hindered effective management policies (NEFSC, 1999). Here critical life history rates of age and growth have been estimated over large proportions of both species ranges that reinforce concerns over the potential for overexploitation in these species. 90 Table 1. General linear model results for A) little skate and B) winter skate regional weight vs. length relationships. W max is the largest individual weight; l max is the longest skate (TL). Also provided in the table are orthogonal contrasts for tests for regional differences in weight length relationships. A. Little skate Regions N w max (g) l max ( cm) intercept slope F p N.E. coast 2342 1.19 57 0.004 3.1 49697.2 0.0001 Mid 1093 0.84 53 0.004 3.12 23547.30 0.0001 SNE-GB 1209 1.19 55 0.004 3.12 31603.30 0.0001 GOM 54 1.14 57 0.003 3.35 1407.39 0.0001 Male 1119 1.14 57 0.005 3.10 27865.5 0.0001 Female 1223 1.19 56 0.004 3.12 24277.1 0.0001 Contrasts F P GOM vs. MID 9.28 0.002 GOM vs. SNE-GB 9.13 0.002 MID VS. SNE-GB 0.01 0.904 B. Winter skate Regions n w max (g) l max (cm) intercept slope F p N.E coast 1711 10.28 111 0.003 3.33 83922.00 0.0001 Mid 121 7.21 94 0.003 3.33 5236.87 0.0001 SNE-GB 1555 10.28 111 0.003 3.33 80215.50 0.0001 GOM 37 8.81 99 0.003 3.35 881.47 0.0001 Male 676 10.28 111 0.003 2.28 43340.40 0.0001 Female 1035 6.53 93 0.003 3.39 4438.11 0.0001 Contrasts F P GOM vs. MID 0.090 0.762 GOM vs. SNE-GB 0.010 0.909 MID VS. SNE-GB 0.110 0.735 91 Table 2. Little skate and winter skate growth parameters. Coast indicates region- wide estimates, L ? is the asymptotic size, k is the growth rate, t 0 is the theoretical size at hatch, F is the estimated F statistic, p is test probability, n is sample size, Mid is mid-Atlantic, SNE-GB is the southern New England-Georges Bank, GOM is the Gulf of Maine, H.R. is higher readability, Con. is conservative read. ?1? indicates model was fit with a set t 0 to equal age at hatch of 11.2 and 16 for little and winter skate respectively, ?2? indicates that values were not provided. Instead an r 2 of .95 and Standard error of 0.001 was provided (data from Sulikowski et at. 2003). L ? (cm) k yr -1 t 0 yr F p n Little skate NE Coast 56.10 0.19 -1.17 694.18 0.0001 236 Male 60.13 0.17 -1.16 271.98 0.0001 94 Female 53.94 0.20 -1.22 326.61 0.0001 125 Mid 53.26 0.22 -1.04 274.12 0.0001 99 SNE-GB 54.34 0.20 -1.22 233.83 0.0004 90 GOM 1 59.30 0.18 -1.15 3805.98 0.0001 47 Winter skate NE Coast 122.10 0.07 -2.06 289.46 0.0001 229 Male 1 115.35 0.08 -1.83 1348.32 0.0001 75 Female 1 114.09 0.07 -2.10 2604.00 0.0001 126 HR 1 111.43 0.09 -1.80 2016.77 0.0001 75 Conservative 108.29 0.11 -1.40 5111.11 0.0001 239 GOM 2 131.4 0.64 -1.53 * * 209 92 Table 3. Correlations between parameter estimates of the von Bertalanffy growth equation where, l ? is the asymptotic length, t 0 is the theoretical size at birth and k is the growth rate. Little skate l ? to k l ? 1.00 -0.67 -0.96 t 0 -0.67 1.00 0.81 k -0.96 0.81 1.00 Winter skate l ? 1.00 -0.83 -0.99 t 0 -0.83 1.00 0.90 k -0.99 0.90 1.00 93 Figure 1. Sampling locations at which little skate vertebrae were sampled along the N.E. coast. The data are colour coded by region with green dots indicating the mid- Atlantic, orange dots indicating southern New England-Georges Bank and white dots indicating the Gulf of Maine. 94 Figure 2. Locations at which winter skate vertebrae were sampled along the N.E. coast. Orange dots indicate sampling locations. 95 0 0.2 0.4 0.6 0.8 1 1.2 1.4 102030405060 Length (cm) Wei g ht (kg) Males Females Females Males Figure 3. Allometric relationship between length (cm) and weight (kg) for little skate. Data are shown for N.E. coast separated by males and females. 96 0 2 4 6 8 10 12 14 10 20 30 40 50 60 70 80 90 100 110 120 Length (cm) We i g h t ( k g ) Males Females Males Females Figure 4. Allometric relationship between length (cm) and weight (kg) for winter skate. Data are shown for N.E. coast. 97 Figure 5. Sample images of centra for A) the birthmark of a winter skate, B) and C) two examples of winter skate centrum sections and D) a little skate centrum section. C. Winter skate, L = 63 cm, age = 9-10. Readability = 5. A. Birthmark D. Little skate, L = 53 cm, age = 11-12. Readability = 4. B. Winter skate, L = 55 cm, age = 4. Readability = 6. 98 Figure 6. von Bertalanffy growth model fit to length (cm) and age (years) data for little skate (see Table 1 for model statistics). The figure shows data for the mid- Atlantic (circle), southern New England-Georges bank (triangle), and Gulf of Maine (diamond) regions. 0 10 20 30 40 50 60 -2 0 2 4 6 8 101214 Age (yr) L e ngth (cm) N.E. coast Mid-Atlantic Southern New England-Georges Bank Gulf of Maine 99 50 51 52 53 54 55 56 57 58 59 60 mid-Atlantic southern New England- Georges Bank Gulf of Maine Asymptotic size (cm) 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 Grow th rate ( k ) Asymptotic size Growth rate (k) Figure 7. Relationship between von Bertalanffy asymptotic size (l ? ) and Brody growth coefficient (k) shown for little skate from three regions. 100 Figure 8. von Bertalanffy growth model fit to length (cm) and age (year) data (see Table 1 for model statistics). The figure shows data for male and female winter skate for the N.E. Coast. 0 10 20 30 40 50 60 -202468101214 Age (yr) L e ngth (cm) Male Male Females Females 101 Figure 9. Comparison of current estimates of age/growth rates and previous studies. Based on my data growth models are shown for N.E. coast and compared to previous work for southern New England-Georges bank (SNE-GB), Block Island Sound (B.I.S.), Long Island Sound (L.I.S.) and Natansons laboratory data. 0 10 20 30 40 50 60 -202468101214 Age (yr) L e ngth (cm) N.E. coast - present study Southern New England-Goerges Bank Long and Block Island sounds Block Island Sound Natanson data 1993 102 Figure 10. von Bertalanffy growth model fit to length (cm) and age (years) data for winter skate (see Table 1 for model statistics). The figure shows data for the mid- Atlantic (circle) southern New England-Georges bank (triangle), and Gulf of Maine (diamond) regions. 0 20 40 60 80 100 120 -2 0 2 4 6 8 10 12 14 16 18 20 22 Age (yr) Lengt h (cm) N.E. coast Mid-Atlantic Southern New England-Georges Bank Gulf of Maine 103 Figure 11. von Bertalanffy growth model fit to length (cm) and age (year) data for winter skate (see Table 1 for model statistics). The figure shows data for male and female winter skate for the N.E. coast. 0 20 40 60 80 100 120 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 Age (yr) Length ( c m) Males Females Females Males 104 Figure 12. Comparison of growth models fit to liberal (solid line), conservative (circles, dashed line), and high readability (diamonds, smaller dashed line) data. See Table 1 for details. 0 20 40 60 80 100 120 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 Age (yr) L e ngth (cm) Liberal read Conservative read High readability Con. data H.R. data 105 Chapter 5 MATURATION, FECUNDITY, LATITUDINAL PATTERNS, AND REPRODUCTIVE STRATEGIES OF TWO RAJIDAE SPECIES IN THE NORTHWESTERN ATLANTIC: LITTLE SKATE, LEUCORAJA ERINACEA, AND WINTER SKATE, LEUCORAJA OCELLATA Introduction Population growth rates are key parameters for both developing management strategies (Myers et al., 1999) and understanding the evolution of life histories (Roff, 2001). Indeed Myers et al. (1999) suggest that knowledge of the slope of the stock recruitment function at the origin, a measure of the potential population growth rate, is the single most important parameter to estimate to develop sustainable management regimes. In a management and evolutionary context population growth rates (e.g., R o , r) are very sensitive to age at maturation (Stearns, 1992). Thus many life history traits are correlated with age at maturity. Generally species with delayed maturation are characterized by long generation times and low productivity and have been termed, ?slow? species, whereas ?fast? species have early maturation and fast growth (Fenchel, 1974; Musick, 1999; Denny et al., 2002). Ages at maturity that optimize evolutionary fitness and those that arise under fishing pressure can be quite different. Generally, species with ?fast? life histories perform well under exploitation and ?slow? species tend to be vulnerable to population declines. Previous researchers (e.g., Winemiller and Rose, 1992) have considered elasmobranch fishes to be the epitome of ?slow? species. As such, elasmobranchs optimize fitness by large maximum body size, slow growth, late maturation, long life-spans and high parental investment in offspring (Walker and Hislop, 1998; Dulvy et al., 2000; Stevens et al., 2000; Frisk et al., 2001). The ?slow? 106 life histories of sharks, skates and rays means that they are likely susceptible to population decline under exploitation and will display extended recovery times if perturbed (Casey and Myers, 1998; Musick, 1999). Additionally, elasmobranchs have an extended juvenile period which increases the chance of being harvested before the onset of maturation making them among the most vulnerable taxa to overexploitation (Holden, 1973; Brander, 1981; Casey and Myers, 1998; Walker and Hislop, 1998; Dulvy et al., 2000; Stevens et al., 2000). However, this broad generalization regarding the characteristic ?slow? life histories of elasmobranchs should not be taken to imply that reproductive strategies within this order are uniform. Within elasmobranchs, faster and slower species can be identified (Chapter 3). Little skate (Leucoraja erinacea) and winter skate (L. ocellata) are sibling species sharing similar phylogeny, body plan, habitat and life history (McEachran and Musick, 1973; Gabriel, 1992; Alvarado Bremer et al., 2004; Frisk et al., 2004). These two closely related skate species exhibit differing positions on the ?fast-slow? continuum of elasmobranch life histories (Frisk et al., 2004). Little skate is a small (TL = 57 cm), moderately fast growing species (von Bertalanffy grow rate of k = 0.20) and moderately lived (12 years), while winter skate is larger (TL = 116 cm), slow growing (k = 0.08) and long lived (20+ years) (Chapter 4). Both species are the most frequently taken skates in commercial fisheries in the western Atlantic (NEFSC, 1999). Winter skate is taken for consumption of its wings (pectoral fins), while little skate is used mostly for fishmeal and bait (Waring, 1984). Yet, scientific surveys 107 have indicated little skate has maintained a relatively stable population whereas winter skate has gone through fairly substantial population fluxes (NEFSC, 1999). Even though these species are commonly caught, economically valuable and represent a large portion of the demersal fish community in the western Atlantic little is known about their biology. Little skate and winter skate occur along the western Atlantic coast from Newfoundland to North Carolina (McEachran, 2002). Little skate has a more southern distribution and is less common than winter skate above La Have Bank (McEachran, 2002). Both species occur on sandy to gravelly bottoms and are common in depths up to 111 meters (McEachran, 2002). Like all Rajidae, both species have internal fertilization, are oviparous, and produce few large eggs which are protected by a rectangular egg case with long tendrils extending from all corners. These extensions allow the protective case to anchor safely to objects and seaweed on the seafloor. The egg develops in the egg case for several months, eventually hatching as a fully-formed juvenile skate. Two definitions of female maturation in Rajidae can be recognized: physiological and functional. Physiological maturity is recognized by the production of eggs of any size and can be measured by gland width, histologically or by the presence of reservoir eggs (Nolan et al., 2002; Sosebee, 2002). Functional maturity is recognized by the presence of ripe eggs that will likely be spawned in the ensuing spawning season. While physiological maturity has been estimated in little skate and winter skate, (Sosebee, 2002), functional maturity has not been estimated. In some species there may be a significant time delay between ages at physiological and 108 functional maturity that may result in substantial differences on subsequent management reference points. Although unable to determine annual fecundity, Fritz and Daiber (1963), using samples from a Delaware estuary, determined that little skate and clearnose skate, R. eglanteria, demonstrated a peak in spawning although some reproduction occurred over the entire year. Johnson (1979), working with little skate collected in Long Island Sound and Block Island Sound, found that maturing eggs ranged in size from 5.0 mm to 26.8 mm, with mature eggs starting at 15 mm. He approximated annual fecundity at 30 eggs and indicated that egg production peaks in the fall and summer (Johnson, 1979). Spawning season for winter skate is not well understood, but appears to have peak egg production during the summer and fall (McEachran, 2002). Size at hatch is approximately 19 and 15 cm for little and winter skate, respectively (Chapter 4). Previous estimates of maturity in winter skate suggest that size at maturation increases with latitude (McEachran, 2002). It is not uncommon for marine species in northern sections of their range to exhibit slower growth, later age at maturation and increased longevity (Taylor, 1958; Beverton and Holt, 1959; Jennings and Beverton, 1991; Beverton, 1992). While spatial patterns can be complex, such as the ?counter- gradient? hypothesis proposed by Conover and Present (1990), there is commonly a continuous latitudinal pattern in vital rates (Beverton and Holt, 1959; Conover and Present, 1990; Jennings and Beverton, 1991; Beverton, 1992). These trends indicate that over a species range there are trade-offs that are expressed in the plasticity of vital rates in mortality, growth and longevity. It is possible that observed variation in 109 maximum body size and size at maturation indicate that a latitudinal trend may exist in winter skate. Here I examine reproductive patterns in life history traits in little skate and winter skate ovaries from the spring and fall of 1999 and summer, fall, winter and spring of 2001-2 in the western Atlantic from Cape Hatteras to Canadian waters. Size of physiological and functional maturity will be estimated for little skate and winter skate based on the presence of physiologically mature and functionally mature eggs. Additionally estimates of von Bertalanffy growth rates estimated previously by Frisk (chapter 4) will be used to estimate age of physiological and functional maturity. Latitudinal effects on sizes and ages of maturation will be analyzed for both species. Analyses of egg production over the year will be assessed to estimate potential annual fecundity for little skate and winter skate. Previous analyses have indicated that fecundity may not vary in elasmobranchs and may be relatively fixed due to phenotypic canalization (Brander, 1981; Frisk et al. 2002; Frisk et al., 2004). If fecundity is not elastic, survival and maturation rates may potentially show compensatory changes in elasmobranchs. In this paper the net reproductive rate (R 0 ) will be estimated for optimal age at maturation for little skate and winter skate. Additionally the age specific egg production rate (l t ?m t ) will be calculated assuming a population with no fishing and under-exploitation. I hope this will provide insight into potential compensatory changes in vital rates in little skate and winter skate. Specifically the goals of this study are to: (1) estimate physiological and functional length and age at maturation along each species range; (2) estimate annual 110 egg production; (3), predict optimal age at maturation based on mortality and fecundity rates; and (4), compare the reproductive strategies of two sibling species, and ( 5) place results in a management context. 111 Methods Sampling: Samples of little skate and winter skate were collected from North Carolina to Canadian waters (Figures 1 & 2). Most sampling was conducted as part of the National Marine Fisheries Service's (NMFS) annual fall (September and October: 1999, 2001), winter (2001), spring (February, March and April: 1999, 2002) and summer (June, July and August: 2001) surveys conducted from the NOAA R/V Albatross IV. During the winter, spring and fall surveys, samples were collected using a bottom trawl with ? in (1.27 cm) mesh liner towed for 30 minutes at three knots. Summer samples were collected using a standard New Bedford scallop dredge with a 2-inch ring chain bag and 1-1/2 inch mesh lining. A total of 103 winter skate were collected in 1999 during the NMFS fall and spring surveys. During the 2001-2 NMFS surveys, 1,844 female little skate and 1,050 winter skate specimens were collected. To ensure broad geographic coverage, an additional sample of 40 little skate was collected from the mid-Atlantic Bight over three dates in the winter (12-12- 00), spring (3-29-01) and summer (5-25-01) of 2000/2001 from the catch of the F/V Tony & Jane; a 57?foot scalloper registered in Ocean City, Maryland and captained by Mr. J. Eustler. Length and weight of all individual skates were recorded aboard the ship. Ovaries were removed, frozen, and shipped to the Chesapeake Biological Lab. Analyses were conducted on the entire dataset and on regional subsets of the data. For non-spatial analyses, I refer to the entire area sampled as the northeast (N.E.) coast. Regional subsets of the data were identified based on likely 112 biogeographic ecotones within which vital rates could be expected to be more similar than between regions. The three locations used in the analyses were: (1) mid- Atlantic, defined here as areas north of Cape Hatteras to the Hudson river canyon, (2) southern New England-Georges Bank, defined here as north of the Hudson River canon to the northern edges of Georges Bank; and (3) the Gulf of Maine, defined here as the area from the Northeast channel southwest to Cape Cod. Laboratory procedures Little and winter skate specimens were processed to ensure adequate coverage of geographical regions and maturing and mature length classes. Ovaries were thawed and fixed in Bouin?s solution. The fixed weight of each ovary was recorded. Eggs less than 5 mm, although numerous, were difficult to measure and fragile. Eggs greater than 5 mm diameter were enumerated, and their individual diameters estimated as the average of three measurements (nearest 0.1 mm) using vernier calipers. In a few cases (< 1%), eggs larger than 5 mm were broken and the diameter had to be estimated. Fixed ovary weight was regressed on female total length (TL) for both species. Analysis of covariance and planned orthogonal contrasts were used to test for regional differences. Maturation analysis Individual skate were scored as immature, physiologically mature or functionally mature based on the characteristics of eggs when found. For my 113 analyses, female little skate and winter skate were scored as physiologically mature if they were carrying eggs 5 mm in diameter or larger. Functional maturity was characterized by the presence of mature eggs. Mature eggs were defined as eggs that were as large or larger than those found in the developing egg case or located in the shell gland or oviducts immediately prior to release (Figure 3). The diameter of mature eggs was estimated from samples of eggs in the shell glands or within egg cases found in the oviducts of the female. When estimating functional maturity in a serial spawning species such as both skate species studied, some mature individuals will be in a resting phase at any point in the year and carry no mature eggs. Such individuals would be incorrectly scored as physiologically mature but not functionally mature. To reduce this potential source of error in determining maturity status, the relationship between ovary weight and total length was used as a secondary indicator of maturity. Based on these relationships, little skate were defined as mature if they had ovaries > 3 g (weight), and winter skate were defined as mature if they were > 5 g (weight). These ovary weight thresholds correspond to the maturity stage where ovaries rapidly increase in size preparing for mature egg production. Size-dependent patterns in physiological and functional maturity in both species were examined using logistic regression. The maturity status of a sufficient number of little skate was determined to permit analyses to be conducted at the regional level. The presence of regional differences in size at maturity for little skate was examined using pairwise, planned orthogonal contrasts following logistic 114 regression. Constraints on the number of female winter skate collected meant that the size-dependence of maturity was determined coastwide and not at the regional level. To estimate age-at-maturity, size-dependent probability of maturity functions predicted from logistic regression and expected age-at-size curves developed from analysis of vertebral centra were overlaid. Fecundity analyses I estimated annual egg production as the product of fecundity and spawning frequency. Laboratory observations of little skate (Richards et al., 1963; Johnson, 1979) indicate that individual skate spawn every 5-6 weeks. Accordingly, annual egg production rates of little and winter skate were calculated as the average fecundity multiplied by the expected number of annual spawning bouts (=52/5 and 52/6). I varied the size of functional maturity to determine how the estimated size of mature eggs affects the resulting annual fecundity estimates. Optimal age at maturation Optimal age at maturity was estimated from age-specific schedules of survival and reproduction. The net reproductive rate, R 0 , was calculated as: (1) ? ? = ???? ?=?= mat mat Tx TM t tM t m M mR exp)( 1 exp)( 0 where T mat is age of maturity, M is natural mortality and m t is age specific fecundity (Roff, 2001). The theoretical age of maturity is estimated as the T mat that satisfies: 115 (2) .0 0 = ? ? mat T R This condition can be estimated by plotting R 0 against increasing age at maturation (T mat ) and finding the point of the function where R 0 is maximized. Net reproductive rate, R 0 , was calculated based on estimates of M and hence l t , using Hoenig?s method. Longevities of 12 and 20 were used for little skate and winter skate, respectively (Chapter 4). M was estimated to be 0.35 for little skate and 0.22 for winter skate. Length at age and fecundity at age relationships were used to estimate m t for 1 cm intervals and the 52+ cm group for little skate. Length at age and fecundity at age relationships were used to estimate m t for 1 cm intervals with a 86+ cm group for winter skate. Egg production at age data for winter skate was noisy and 4-year averages were applied to smooth extreme values. Age-specific egg production (l t ?m t ) was calculated to estimate age of peak egg production in little skate and winter skate which provide a estimate of theoretical age at maturation. 116 Results Maturation: Little skate ovary weight (OW) was significantly related to total length (Figure 4 -Table 1). Regression relationships suggested a latitudinal gradient in the exponent of the relationship, but these differences were not significant (Table 1). Similarly, ovary weight was significantly related to length in winter skate (Figure 5, Table 1). An ovary weight of ~ 3 grams in little skate and ~ 18 grams in winter skate appear to relate to maturity stage when mature egg formation proceeds (Figures 4 and 5). In little skate, average mature egg size was 18.5 mm (SD = 2.6, range = 13-23 mm, n = 14). Winter skate mature eggs were larger than little skate eggs. Mature egg size of winter skate was 23.52 mm based on four specimens with egg diameters of 29.9 mm and 27.1 mm, 24.6 mm and 24.3 mm, and one specimen with a single egg case with a diameter of 21.7 mm. Based on these data, eggs with diameters 10 and 15 mm and greater were assumed to be mature for little skate and winter skate, respectively. These mature eggs can be expected to be released in the ensuing spawning season. The smallest little skate found to contain eggs > 5 mm diameter was 38 cm (TL); however, most females began egg production at a length of 38-43 cm and larger (Figure 6). The sizes of winter skate found with eggs > 5mm were larger than in little skate (Figure 7). The smallest winter skate with eggs > 5 mm was 60 cm, with most beginning production at a size of > 65 cm (Figure 7). Physiological maturity was significantly related to female length in little skate with trends indicating that size of 117 maturity increases with latitude (Figure 8-Table 2). Size of 50% physiological maturity in little skate was TL = 42.5 cm, 43.0 cm and 46.5 cm in the mid-Atlantic, southern New England-Georges Bank and the Gulf of Maine regions respectively (Figure 8). The overall test of regional differences was significant (Chi-square = 16.74, p = 0.0002). All regional pairwise comparisons were significant (Table 2). In winter skate, physiological maturity was significantly related to length with 50% physiological maturity occurring at a TL = 66 cm for the N.E. coast (n = 212, chi- square = 122.72, p = 0.0001) (Figure 9). There were not adequate data to observe regional trends in winter skate. Functional maturity in little skate was estimated to be a female carrying eggs 10 mm or larger or if combined ovary weight was 3 g or larger. Similarly, functional maturity was defined in winter skate as females containing either an egg 15 mm or larger or the combined ovary weight was 18 g or larger. Based on these definitions, the production of mature eggs appears to increase in a knife-edge function of size for little skate and to a lesser extent for winter skate (Figures 10 & 11). The proportion of females with mature eggs provides an estimate of the average proportion of females spawning at a given time. Proportions were calculated for intervals of 2 and 3 cm and a plus group for little skate and winter skate respectively. The sample size at length (cm) for little skate was: 41-2 cm (n = 21), 43-4(34), 45-6(57), 47-8(72), 48- 49(47), 50-1(29) and the plus group of 52+(11); and for winter skate: 69-70(n = 15), 71-2(19), 73-4(20), 75-6(22), 77-8(14), 79-80(17), 81-2(12), 83-4(9) and the plus group of 85+(12). On average the proportion of mature female little skate that are 118 carrying a mature egg is 0.73 (95% C.I. ? 0.11) and the proportion for winter skate is 0.55 (95% C.I. ? 0.043). Logistic regression indicated that sizes at 50% functional maturity in little skate occurred at larger lengths at higher latitudes (Figure 12). Size of 50% functional maturity was TL = 43 cm, 44 cm and 46 cm for the mid-Atlantic, southern New England-Georges Bank and the Gulf of Maine regions respectively (Figure 12; Table 2). However, these differences were not significantly different (Chi-square = 2.93, p = 0.23). Using growth data from the von Bertalanffy equation estimated previously by Frisk (Chapter 4) approximate age of functional maturity was T mat = 6.5-7 for the N.E. coast. Thus, in little skate, there does not appear to be any delay between physiological and functional maturity. Estimated size at 50% functional maturity in winter skate occurred at TL = 76 cm (n = 212, chi-square = 113.53, p = 0.0001) (Figure 9, Table 2). Using previous age estimates physiological maturity occurs at approximately age of T mat = 9.5 and functional maturity at age of T mat = 12.5 (Figure 9). In winter skate, there is a three- year delay between the onset of physiological maturity and functional maturity. Fecundity: Depending on the size chosen to define mature eggs, estimates of annual fecundity ranged from 21 - 57 eggs for little skate (Table 3) and 26 - 100 eggs for winter skate (Table 4). Little skate egg production varied little seasonally, although higher values were observed in the spring and fall (Table 5A). Winter skate egg 119 production was highest in the fall; however there was insufficient seasonal coverage to estimate season egg production (Table 5B). The distribution of eggs 5 mm and larger are shown in Figures 13 and 14 for functionally mature little skate greater than 47 cm in total length and all functionally mature winter skate. Individuals appear to be resting or actively laying eggs. Winter skate appears to produce more reservoir eggs and mature eggs than little skate. Also, indicated in the figures is that the number of reservoir eggs and mature eggs increased very little with length in both species. However, in both species the very large individual females may have a slightly higher number of reservoir eggs. Optimal age at maturation: R 0 was optimal at T mat = 7 for little skate (Figure 15), which is nearly identical to my observed age of maturity of 6.5-7. In contrast, for winter skate R 0 was optimal for T mat = 12.5 (Figure 16), identical to the observed age of maturity of 12.5. Age-specific egg production was estimated for little skate and winter skate for populations experiencing either no fishing or a fishing mortality rate F = 0.35. Little skate peak egg production occurred at 7.50 and 6.75 with no fishing and fishing, respectively (Figure 17). Winter skate peak egg production occurs at 14.5 and 10 for no fishing and fishing, respectively (Figure 18). In both species exploitation decreases the age of peak egg production. 120 Discussion The reproductive ecologies of little and winter skate differ substantially with both species representing different locations on the ?fast-slow? continuum of life histories in elasmobranchs. This difference is reflected in the vital rates such as length and age of maturity, annual fecundity and the net reproductive rate (R 0 ). Not only have I documented differences in vital rates between sibling species but I have also found a latitudinal gradient in size of maturation in little skate. What features of the life history variation in these two species account for the observed differences? Life history theory predicts that delayed maturity is associated with additional growth leading to a larger size, increased fecundity and/or greater quality or size of offspring (Stearns, 1992). Within species, winter skate and little skate do not show large increases in fecundity with size. Further, age at maturity appears relatively invariant in little skate, while size at maturity increased with latitude. This is not compelling evidence for selective pressure towards delayed maturation. However, when viewed between species, the benefits of delayed maturation hold up well. Winter skate delays maturation, having a longer time to grow and reach a larger size with higher fecundity than little skate. However, the cost of delayed maturity is lower growth rates, longer generation times and likely lower compensatory ability; thus less resilience to exploitation (Frisk et al., 2002; Frisk et al., 2004). I estimated differences in sizes at physiological and functional maturity and fecundity in both species. I also found evidence for latitudinal gradients in reproductive ecologies. McEachran (2002) reported that size of maturity ranges from 35-50 cm for little skate and 70-109 cm for winter skate and generally increases with 121 latitude (McEachran and Martin, 1977). My data indicated that sizes at 50% physiological maturity and functional maturity varied latitudinally for little skate. Sizes at 50% physiological maturity and 50% functional maturity in the mid-Atlantic were 42.5 and 43 cm. Equivalent estimates for the southern New England-Georges Bank were 43 and 44 and 46 and 46 cm for the Gulf of Maine. These sizes are equivalent to ages of 5.75 to 8.2 years. Importantly, there was no evidence to suggest a temporal delay between physiological and functional maturity. In winter skate, there was no evidence of a latitudinal gradient in sizes at maturity as a result of invariance of estimates and a paucity of samples in the Gulf of Maine and mid- Atlantic regions. Physiological maturity and functional maturity occurred at 66 and 76 cm throughout N.E. coast, equivalent to age of 9.5 and 12.5 years. Apparently the onset of functional maturity is delayed three years relative to physiological maturity in winter skate. This difference in the two measures of maturity suggests that in some skate species, estimating physiological maturation only may drastically under- estimate the age and size at which individuals actually contribute to the reproductive effort of the population. Ovary weight and size at maturation increased with latitude in little skate from the mid-Atlantic, to southern New England-Georges Bank to the Gulf of Maine, while age at maturation appeared to remain constant. Frisk (Chapter 4) found some evidence that life-span was higher and growth rate slower north of the mid-Atlantic region for little skate. Commonly, a larger size of maturation is associated with later age at maturation (Stearns, 1992). Here, I have only found a difference in size of maturity and not age, perhaps reflecting imprecision in ageing of skates. While 122 McEachran and Musick (1977) previously reported that winter skate showed evidence for a latitudinal gradient in size of maturation, I found no evidence that size of maturation increases with latitude. It is not surprising that, of the two species, little skate showed a latitudinal trend in vital rates and winter skate did not. Genetic analyses have indicated that little skate has much higher nucleotide diversity than winter skate (Alvarado Bremer et al., 2004, Alvarado-Bremer, personal communication). Winter skate appears to have a much less nucleotide diversity and apparently faced a genetic bottleneck at some time in its history, possibly caused by exploitation (Alvarado Bremer et al. 2004, Alvarado-Bremer, personal communication). McEachran and Musick?s (1977) work was conducted in 1968-69 and it is possible that the population diversity of winter skate declined during periods of high exploitation (i.e. 1960?s and 1970?s NEFSC, 1999; Fogarty and Murawski, 1998). Further DNA analyses such as microsatellite loci work is needed to determine whether and when a bottleneck may have occurred and if there are genetically separate populations in little skate or winter skate. Optimal age at maturation in little skate and winter skate is nearly identical to observed age at maturity, indicating that selection has acted in both species to maximize growth rate. Peak egg production in winter skate occurs at an age of 14.5. However, when the effects of exploitation are included, peak egg production occurs at an age of 10.5. A less dramatic decrease in age of peak egg production is also seen in little skate. In little skate and especially winter skate exploitation favors selection for early maturation. 123 Historically, winter skate has experienced fishing mortality as high or higher than F = 0.35 (NEFSC, 1999). Little skate has been harvested in the same multi- species fishery as winter skate; however, little skate is smaller and may not be fully recruited to the fishery until it reaches mature sizes. Thus, it is likely that winter skate has experienced a stronger selective force to mature earlier compared to little skate. Gear selectivity often removes the faster growing individuals leaving slower growing individuals and potentially changing the demographics of the population. When exploitation is considered, optimize fitness would be achieved at a lower age at maturation for winter skate. However, the actual impacts of exploitation may counter the force of selection for younger age at maturation. Gear selectivity may favor slower growing winter skate that potentially have a larger and older age at maturity. Information on egg production either in the laboratory or the wild in winter skate is absent. Thus, out of necessity I used egg production rates estimated for little skate in order to estimate fecundity in winter skate. Similar laboratory egg production rates were observed for clearnose skate, Raja eglanteria, which produced a pair of eggs every four days (Libby, 1959). My data indicate that the spawning season in both skates is prolonged, and spawning likely occurs year round. Furthermore, my data indicates that egg production is high in the fall and in the spring for winter skate and relatively even for little skate with increased production of mature eggs in the summer. Annual fecundity was estimated to range from 21 - 57 eggs for little skate and 26 - 100 eggs for winter skate. These values are in the range 124 found for other skates species, Dipturus batis (40), Amblyraja radiata (17), Raja brachyura (40), Raja neavus (80), and Raja montagui (24), and Raja clavata (52) (Holden et al., 1971, Walker, 1994). My results for fecundity in little skate are in the high range or higher than Johnson?s estimates of 26 - 33 eggs annually for little skate. Little skate is believed to spawn and undergo a resting period of three to four months (Johnson, 1979). During a reproductive period, little skate is believed to produce eggs for 5-6 weeks, (Johnson, 1979), with a pair of eggs produced in a minimum of 5 days and commonly over a week (Richards et al., 1963). Johnson (1979) indicates that female little skate have two reproductive periods a year. My fecundity analysis averages egg production over the year for little skate and winter skate populations. Thus, any resting period(s) is/are indirectly accounted for in the estimation of fecundity, assuming mature eggs are released in the ensuing spawning season. This was necessary as both species in the current study are serial spawners and do not appear to have strong seasonal patterns previously observed in little skate (Richards et al. 1963; Johnson, 1979). A question of elasmobranch life history that remains unanswered is why evolution has favored larger body size (Frisk at al., 2004)? One hypothesis is that larger body size increases the volume for gestation (Qaulls and Shine, 1998; Goodwin et al., 2002). My results indicate that the larger of the sibling species, winter skate, has higher fecundity than little skate. Frisk et al. (2002) plotted the intrinsic rate of increase against fecundity from results of matrix models and found that increases in fecundity beyond approximately 17, 18, and 23 for little skate, winter skate and barndoor skate, respectively, had diminishing returns (assuming a 50/50 sex ratio). 125 These values roughly correspond to annual fecundity estimates indicating that there is likely little benefit in increasing fecundity for these three Rajidae species. Bonfil (1996) provided theoretical evidence, based on simulations of an age- structured model, that changes in fecundity have little potential for compensatory mechanisms and more likely candidates are mortality and age at maturation. If, as the evidence indicates, that fecundity is fixed by body cavity volume (i.e. body size), in Rajidae, and perhaps elasmobranchs, compensatory changes may occur with changes in age at maturation as has been observed in the sharpnose shark, Rhizoprionodon porosus, and spiny dogfish, Squalus acanthias, (Carlson and Baremore, 2002; Sosebee, 2002). My prediction of optimal age at maturation in winter skate provides evidence that exploitation could be compensated for by lowering age at maturation. Although, our knowledge of the population dynamics of skates is increasing, (Walker and Hislop, 1998; Stevens et al., 2000; Frisk et al., 2000; Frisk et al., 2002), much work is needed to understand the varying resiliency to exploitation of Rajidae species (Walker and Hislop, 1998; Frisk et al., 2002). In addition to management concerns, work remains in understanding life history and evolutionary trends in elasmobranchs, of which many species display biological demographics on the extreme of the equilibrium strategy in the ordination of life history strategies in fishes (Frisk Chapter 3). Here I analyzed the contrasting life histories of two closely related species occupying different places on the ?fast-slow? life history continuum. Further, understanding these differences between sibling species will help to elucidate the underlying biological trade-offs that result in species-specific resiliency to 126 exploitation. 127 Table 1. A. Regression summary and parameter estimates for ovary weight vs. length relationships in little skate and winter skate. Data for little skate and winter skate where (n) is sample size, (L max ) is maximum length observed and (W max ) is maximum ovary weight observed. Parameters of the ovary weight vs. length relationship are provided where a is the intercept and b is the shape parameter and were fit to log transformed data. SE is standard error and p is the models significance probability. B. Results from orthogonal contrasts of regional ovary weight vs. lengths models. No significant contrasts were found. A. B. Little skate n L max (cm) W max (g) a b SE (b) p Mid 92 52 26.7 -14.44 9.2 0.95 0.0001 SNE-GB 142 53 29.1 -13.27 8.4 0.65 0.0001 GOM 16 56 21.1 -15.84 9.8 1.60 0.0001 Winter skate NE Coast 207 93 192 -9.68 5.8 0.37 0.0001 Little skate Regions F p GOM VS. MID 0.16 0.69 GOM VS. SNE-GB 0.61 0.43 MID VS. SNE-GB 0.39 0.53 128 Table 2. Summary of logistic regression analyses for regional maturity comparisons of little skate. A. The logistic model statistics for physiological and functional maturity in little skate is listed for the mid-Atlantic (MID), southern New England- Georges Bank (SNE-GB) and the Gulf of Maine (GOM). B. the results of orthogonal planned contrasts for physiological maturity for little skate. All pairwise regional contrasts were significant for physiological maturity. Note that there were no significant pairwise regional differences in functional maturity for little skate. A. Physiological maturity Region n chi-square p MID 120 52.96 0.0001 SNE-GB 144 96.54 0.0001 GOM 16 9.99 0.0017 Functional maturity Region n chi-square p MID 92 43.6 0.0001 SNE-GB 144 71.35 0.0001 GOM 16 13.86 0.0002 B. Contrasts for physiological maturity chi-square p SNE-GB VS GOM 4.83 0.03 GOM vs MID 14.34 0.00 MID VS SNE-GB 11.14 0.00 129 Table 3. Fecundity estimates for varying mature egg size and spawning season for little skate. Length is total length of female little skate, n is sample size, > X mm is the number of egg greater then X. Fecundity estimates are shown assuming release periods of 5 and 6 weeks and varying size of mature eggs. Length n >8 mm >9 mm >10 mm > 11 mm > 12 mm 46 40 4.3 3.3 3.2 2.6 2.0 47 32 4.9 3.7 3.0 2.0 1.6 48 34 4.8 3.7 3.2 2.6 1.9 49 13 5.1 4.2 3.9 2.9 2.3 50 15 5.0 3.9 3.2 2.3 2.0 51 14 3.9 3.1 2.5 2.2 1.7 52 8 7.8 5.9 5.4 4.5 3.9 53 2 3.0 2.5 2.0 1.5 1.0 56 1 11.0 11.0 8.0 8.0 5.0 Average 5.5 4.6 3.8 3.2 2.4 Fec. 5 wks 57.37 47.65 39.76 33.04 24.79 Fec. 6 wks 47.81 39.71 33.13 27.53 20.66 130 Table 4. Fecundity estimates for varying mature egg size and spawning season for winter skate. Length is total length of female winter skate, n is sample size, > X mm is the number of egg greater then size X mm. Fecundity estimates are shown assuming release periods of 5 and 6 weeks and varying size of mature eggs. Length of >80 are samples for which length was not recorded but ovary weight indicates they were larger than 80 cm. Length n >12 mm >13 mm >14 mm >15 mm >16 mm >17 mm >18 mm 79 8 7.13 5.38 4.50 3.63 2.88 2.50 2.00 80 9 4.89 4.00 3.11 2.67 1.89 1.56 1.22 81 8 4.25 3.88 3.25 2.75 2.25 2.25 2.00 82 4 14.25 13.00 10.00 8.75 7.00 5.00 3.50 83 1 17.00 13.00 11.00 9.00 9.00 8.00 8.00 84 8 9.50 8.13 6.88 6.63 5.13 4.75 4.50 85 4 6.50 5.75 5.75 5.00 4.75 4.00 3.75 86 1 5.00 2.00 0.00 0.00 0.00 0.00 0.00 87 3 14.33 12.00 9.67 8.00 6.67 5.67 5.67 88 1 5.00 3.00 1.00 1.00 0.00 0.00 0.00 89 1 11.00 10.00 7.00 5.00 3.00 2.00 0.00 91 1 8.00 7.00 5.00 4.00 3.00 1.00 1.00 93 1 10.00 7.00 7.00 5.00 3.00 3.00 3.00 >80 4 18.75 17.75 15.50 13.00 9.75 9.00 7.25 Average 9.69 7.99 6.40 5.32 4.16 3.48 2.99 Fec. (5 wks) 100.73 83.11 66.60 55.28 43.31 36.19 31.12 Fec. (6 wks) 83.94 69.26 55.50 46.07 36.09 30.16 25.93 131 Table 5. Seasonal mature egg production in little skate (A) and winter skate (B). Sampling date, sample size (n), average number and confidence intervals (95%) are listed. A. B. Little skate Date n average 95% C.I. 12.12 8 3.38 1.67 2.7-3.2 53 2.55 0.72 3.10-3.29 32 4.25 1.11 4.7-4.25 17 2.24 1.49 5.25 6 4.50 2.37 6.28-8.15 49 4.94 0.89 9.6-10.18 39 2.95 0.83 Winter skate Date n Average 95% C.I. 9.6-10.18 33 9.70 2.38 2.7-3.2 12 1.92 2.02 6.28-8.15 4 2.25 6.15 3.10-4.25 24 1.62 0.91 132 Figure 1. Regional coverage of little skate ovaries collected along the northeastern United States coast from Cape Hatteras to Canadian waters. Regional level analyses were performed on the mid-Atlantic, representing areas north of Cape Hatteras to the Hudson River canyon; southern New England-Georges Bank, north of the Hudson River canyon to the outer edges of Georges Bank; and the Gulf of Maine, north from the Northeast channel southwest to Cape Cod. The data are color coded by region with green dots indicating the mid-Atlantic, orange dots indicating southern New England-Georges Bank and white dots indicating the Gulf of Maine. 133 Figure 2. Regional coverage of winter skate ovaries collected along the northeastern United States coast from Cape Hatteras to Canadian waters. Regional level analyses were performed on the mid-Atlantic, representing areas north of Cape Hatteras to the Hudson River canyon; southern New England-Georges Bank, north of the Hudson River canon to the outer edges of Georges Bank; and the Gulf of Maine, north from the Northeast channel southwest to Cape Cod. Orange dots indicate sampling locations. 134 Figure 3. The reproductive track of a little skate shown with two ovaries and two shell glands. Mature eggs move from the ovaries to the shell glands where they are encased and moved down the oviduct and released on the sea floor. Johnson et al. 1963 reports that eggs are normally released in pairs, even if one ovary contains no eggs, in that case cases may be empty. n Shell glands n Ovaries 135 0 5 10 15 20 25 30 30 35 40 45 50 55 60 Length (cm) Ovary weight (g) Gulf of Maine Southern New England-Georges Bank Mid-Atlantic Gulf of Maine Southern New England-Georges Bank Mid-Atlantic Figure 4. Ovary weight (g) regressed over individual female total length for little skate. Power functions are fit to three geographical regions including the mid- Atlantic, southern New England-Georges Bank and the Gulf of Maine. At approximately 3 grams it appears that most ovaries are mature and are increasing in weight rapidly. 136 0 20 40 60 80 100 120 140 160 180 200 40 50 60 70 80 90 100 Length (cm) Ovary weight (g) Figure 5. Ovary weight (g) regressed over individual female total length for winter skate. A power function is fit to N.E. coast for winter skate. At approximately 18 grams it appears that most ovaries are mature and are increasing in weight rapidly. 137 0 5 10 15 20 25 30 30 35 40 45 50 55 60 Length (cm) Number of reservoir and mature eggs Figure 6. Number of reservoir and mature eggs for little skate (eggs > 5 mm) by size of skate. Figure shows a steep increase in reservoir and mature egg production in little skate. 138 0 10 20 30 40 50 60 70 40 50 60 70 80 90 100 Length (cm) Number of reservoir and mature eggs Figure 7. Number of reservoir and mature eggs for winter skate (eggs > 5 mm) by size of skate. The figure indicates a steep increase in reservoir and mature egg production in winter skate. 139 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 30 35 40 45 50 55 60 Length (cm) Physiological maturity Mid Atlantic Southern New England-Georges Bank Gulf of Maine Figure 8. Physiological maturity regressed over length for little skate is shown for the mid-Atlantic, southern New England-Georges Bank and the Gulf of Maine. Maturation curves were estimated with presence absence data and a logistic model. 140 Figure 9. Proportion mature and age for given lengths are shown for physiological and functional maturity in winter skate. Age of 50% physiological maturity and functional maturity are shown at the intersection of the drop down line with the second y-axis. Maturation curves were estimated with presence absence data and a logistic model. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 40 50 60 70 80 90 100 Length (cm) Proportion mature 0 5 10 15 20 25 30 Age Pysiologically Mature Functionally Mature 95 % Confidence Intervals Age-N.E. coast Age at functionally mature = 12.5 Age at physiologically mature = 9.5 141 0 2 4 6 8 10 12 14 30 35 40 45 50 55 60 Length (cm) Number of mature eggs 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Proportion of individuals with mature eggs Mature eggs Proportion with mature eggs Figure 10. Mature eggs and proportion of females with mature eggs are shown for female little skate (eggs > 10 mm). 142 0 5 10 15 20 25 40 50 60 70 80 90 100 Length (cm) Number of mature eggs 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Proportion of individuals with mature eggs Figure 11. Mature eggs and proportion of females with mature eggs are shown for female winter skate (eggs > 15 mm). 143 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 30 35 40 45 50 55 60 Length (cm) Proportion functionally mature 0 2 4 6 8 10 12 14 Age Mid-Atlantic Southern New England-Georges Bank Gulf of Maine Age-Mid-Atlantic Age-Southern New England-Georges Bank Age-Gulf of Maine Age at functionally mature (Mid & GOM) = 7 Age at functionally mature (SNE-GB) = 6.5 Figure Figure 12. Proportion mature and age for given lengths are shown for physiological and functional maturity for the mid-Atlantic, southern New England-Georges Bank and the Gulf of Maine for little skate. Age of 50% physiological maturity and functional maturity are shown at the intersection of the drop down line with the second y-axis. Maturation curves were estimated with presence absence data and a logistic model. 144 40 40 41 42 42 43 43 44 44 45 45 46 46 47 47 48 48 49 49 50 50 51 51 52 52 53 53 56 0 2 4 6 8 10 12 Number of eggs Length (cm) Egg size (5-25mm) Figure 13. A 3-D representation of the number of eggs in female little skate plotted against egg size for mature females. Two females are shown for each 1 cm length interval. No apparent increase in egg production can be seen as length increases; although the largest female corresponds to the most eggs. 145 69 71 73 75 77 79 81 83 85 87 89 93 0 5 10 15 20 Number of eggs Total length (cm) Figure 14. A 3-D representation of the number of eggs in female winter skate plotted against egg size with mature females. One individual is shown for each 1 cm interval. No apparent increase in egg production can be seen as length increases; although large individuals appear to have a greater number of eggs. Egg size (5-30mm) 146 Figure 15. Net reproductive rate (R 0 ) is shown for ages at maturity ranging from 0-12 for little skate. Peak net reproductive rate occurs at approximately 7 years in little skate. 0 1 2 3 4 5 6 7 024681012 Age at maturity N e t reproductive rate (R o ) 147 0 1 2 3 4 5 6 7 8 9 4 6 8 10 12 14 16 18 20 22 Age at maturity Net reproductive rate R o Figure 16. Net reproductive rate (R 0 ) is shown for ages of maturity ranging from 0- 16 for winter skate. Peak net reproductive rate occurs at approximately 13 years in winter skate. 148 0 0.5 1 1.5 2 2.5 02468101214 Age A g e s p ec i f i c egg produc t i on w i t h no f i s h i n g (l t m t ) 0 0.05 0.1 0.15 0.2 0.25 E g g p r od u c t i o n w i th fi sh i n g (l t m t ) Unfished Fished Figure 17. Age specific egg production with fishing and no fishing (F= 0.35) for little skate. Note that peak egg production occurs slightly earlier when fishing mortality is assumed for little skate. 149 0 0.5 1 1.5 2 2.5 3 3.5 4 024681012141618 Age Age s p eci f ic egg producti on wi th no fi shing (l t m t ) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Age specific egg produ c tion with no fishing (l t m t ) Unfished Fished Figure 18. Age specific egg production with fishing and no fishing (F= 0.35) for winter skate. Note that peak egg production occurs several years earlier when fishing mortality is assumed for winter skate. 150 Chapter 6 AN AGE-STRUCTURED MODEL OF WINTER SKATE, LEUCORAJA OCELLATA, ABUNDANCE IN THE WESTERN ATLANTIC: SUSTAINABILITY AND UNCERTAINTY Introduction Over the last 30 years, exploitation in the northwestern Atlantic has caused large changes in biomasses of component species biomass and properties of the ecosystem (Fogarty and Murawski, 1998; Link et al., 2002). During the 1970?s and 1980?s many species were declared overfished including haddock (Melanogrammus aeglefinus,), silver hake (Merluccius bilinearis), herring (Clupea harengus), in addition to the well-documented decline of cod (Gadus morhua) (Hutchings & Myers 1994; Fogarty and Murawski, 1998). During this period, landings in the northwestern Atlantic fishery have changed substantially (Fogarty and Murawski, 1998). In part, this change reflects changes in composition of the fishing fleet. During 1960?s and 1970?s, the New England fishery was characterized by the presence of foreign factory trawlers that landed a larger amount and a greater diversity of species than the domestic fleet (Fogarty and Murawski, 1998). Elasmobranchs (skates and dogfish sharks) are one taxonomic group that experienced considerable changes in abundance. Analyses of survey data indicate an increase in elasmobranch biomass in the 1980?s with skates and dogfish accounting for most of the increase (Link et al., 2002, Garrison and Link, 2000). For example, the National Marine Fisheries Service?s (NMFS) annual spring bottom trawl survey indicates that winter skate?s (Leucoraja ocellata) population increased approximately ten-fold from 1979 to 1988 (NEFSC, 1999). Similarly, the NMFS fall survey also 151 shows an approximately a 6-fold increase for this species over the same time period (NEFSC, 1999). In the 1980?s, as commercial groundfish stocks declined due to overfishing (Fogarty and Murawski, 1998), a domestic market developed for skates and landings of these species increased (NEFSC, 1999). In particular, larger skates such as winter skate, thorny skate (Amblyraja radiata) and barndoor skate (Dipturus laevis) were landed and often sold to European markets (McEacran, 2002). Smaller skates such as little skate (L. erinacea) were landed and sold as bait (NEFSC, 1999). Even though their life histories are similar, not all elasmobranch species increased. A review of the abundance trends in the surveys indicates that thorny and barndoor skate remained below their overfishing thresholds set by the federal regulation, that correspond to F = M, for most and in the latter case all of the 1970?s and 1980?s (NEFSC, 1999). In fact, barndoor skate was very rarely caught by NMFS?s surveys in the 1980?s (NEFSC, 1999). Winter skate also dipped below its threshold in the 1990?s and smooth skate (Malacoraja senta) appears to have declined in abundance from highs in the 1970?s (NEFSC, 1999). It has been suggested that skates may have been in greater abundance prior to the beginning of the NMFS survey, confounding the choice of appropriate biomass thresholds for western Atlantic skates (Casey and Myers, 1998). The rapid increase in biomass observed in some elasmobranch species appears to be at odds with the slow growth, delayed maturity, low fecundity and long life- spans that characterize the Class. For example, based on results of a Leslie matrix model for winter skate that assumed a longevity of 20 years, an age of maturity of 9 years, fecundity of 35 and no density-dependence, Frisk et al. (2002) estimated that 152 the rate of population increase (r) was 0.13. Given this rate of increase, a winter skate population would take 19 years to experience a ten-fold increase, assuming exponential growth and no exploitation. Yet, the NMFS?s annual bottom trawl surveys indicate that winter skate abundance, as indexed by survey catch per unit effort (CPUE), increased ten-fold in the spring and six-fold in the fall in less than a decade (NEFSC, 1999). Accordingly, it would seem timely to conduct a more detailed analysis of the dynamics of winter skate to determine likely responses of its population to exploitation. The goals of this paper are to estimate the population productivity and compensatory ability of winter skate from 1963 to 1998 using a single species age- structured model. I will use the model to examine whether the rate and magnitude of increase in population abundances in the 1980?s is to be expected given our knowledge of vital rates of winter skate? Additionally, I will assess whether the rate of increase in the population was the result of ?expected? recruitment levels. To address these questions, an age structured model was constructed using the species vital rates and fishery-independent and ?dependent data from 1963-1998. The model assumed that landings provide an index of catch and effort was assumed to follow that of the groundfishery. Various model structures were employed to address the questions posed above including aggregating survey data, using finely resolved length frequency data and observer error only and mixed error models. Possible compensatory mechanisms including recruitment anomalies were explored. 153 Methods An age-structured population model tuned to relative abundance indices and size composition data was used to assess trends in winter skate abundance. The model used data from fishery-independent surveys, commercial catches, and estimates of vital rates (mortality, growth and fecundity). The model included both observation and process error. The model was constructed using the software program AD Model Builder ?, Sindey B.C., Canada. The full model code is provided in Appendix A. Fishery-independent surveys The spring bottom trawl survey has been conducted since 1968 in offshore regions (Figure 1A & 1B). Inshore stations were added for all regions in 1973 (Figure 1A). The fall survey has been conducted since 1963 for offshore regions in southern New England, Georges Bank and the Gulf of Maine (Figure 1C & 1D). In 1972 inshore stations were added to the fall survey (Figure 1C). The surveys were conducted primarily from the RV?s Albatross IV and Delaware II. Both surveys employ stratified random designs with strata defined by depth and latitude. Depending on weather and other circumstances, the spring and fall surveys each sample approximately 300-400 stations per year. Tows were conducted with a bottom trawl, using ? in (1.27 cm) mesh liner, towed for 30 minutes at 3.5 knots. Processing of catches, including recording the total weight of each species for each tow and sub sampling the catch to estimate individual weights and lengths for all species, was conducted onboard. 154 Time series of winter skate abundances have been developed from aggregated survey CPUE for the fall survey (Figure 2). However, these data do not capture the annual changes in the size frequencies evident in the survey. Figure 3 integrates both abundance and size frequency information. Other fishery-independent data exist for winter skate from a variety of sources. A discussion of the relevance of additional surveys (the Massachusetts Division of Marine Fisheries survey, the Connecticut?s Department of Environmental Protection?s trawl survey and a Canadian Department of Fisheries and Oceans survey) to understanding historical patterns in winter skate abundance is presented in Appendix B. Model data inputs The vital rates of winter skate used in the model were derived in previous dissertation chapters (Table 1). Relative abundance estimates were developed from the NMFS?s annual spring and fall bottom trawl survey (Table 2). To utilize the longest available time series, I used fall survey data for offshore regions (strata 1-40) from 1963 to 1998. Use of only offshore stations is justified as survey trends were very similar to trends for all regions combined. The areas covered included off-shore regions of southern New England, Georges Bank and the Gulf of Maine. Commercial catch estimates from 1963 to 1998 were based on reported landings data. Estimates of effort for the groundfish fishery developed by Fogarty and Murawski (1998) were used as surrogates for effort in the skate fishery (Table 2). Note that effort reported in Fogarty and Murawski (1998) included the years 1960- 155 1992, thus for years 1993-1998 a four year running average of effort was used to estimate missing data. Population dynamics The model projects abundance (numbers at age) in the population forward in time. The section includes equations defining biomass (1), fishing mortality (2), recruitment (3), and propagation of numbers at age (4). Definitions of all parameters used in the model are provided in Table 1. Biomass (B) was calculated as the sum of the products of numbers at age (N a,t ) and mean weight-at-age (w a ) where a = age: (1) a a tat wNB ?= ? ? =1 , Animals age 20 and older were accumulated in a plus group, and it was assumed that asymptotic weight was achieved by age 20. Spawner abundance (S t ) was estimated as biomass as above for a = T mat to the maximum age. Fishing mortality was estimated as the observed catch biomass (C) divided by the average population biomass ( B ): (2) t t t B C F = where C t is catch in year t and tB is average biomass in year t. A Ricker stock recruitment function was used to estimate age-1 recruitment: (3) tt wzS tt eeSN ???= ?? + )( 1,1 ? ? 156 where N 1,t+1 is the number of recruits in year t+1, S t is spawning biomass in year t, ? is the density-independent parameter (or maximum age-1 recruits per unit of spawning biomass), and ? is the density-dependent parameter. The second exponential term (e wz ) represents process error that allows for recruitment variation or anomalies. When assuming an observation error model only, the process error terms were not estimated. The numbers at age were projected forward as: (4) )( ,1,1 exp ta FvM tata NN ?+? ++ = where N a,t is number of fish of age a in year t, M is natural mortality, v a is vulnerability of age a fish to the fishing gear and F t is fishing mortality in year t. Note fishing and natural mortality are not age-specific. Model initialization The model was initialized using a leading parameter scheme. Given the dynamic equations above (eq 1-4), it is possible to reduce the number of initial parameters in this age structured model to two key population parameters that represent the population?s recruit production at equilibrium (R o ) and maximum juvenile survival rate (K). To do this, a series of age-specific relationships based on leading (or unknown) parameters were derived. These age-specific relationships include survivorship, maturity, and growth. Given these age-specific relationships and initial estimates of natural mortality, R o and K, parameters for a stock recruitment relationship for unfished conditions were derived. The age-specific relationships are 157 described in equations 5-8 and the alpha and beta terms of the Ricker recruitment function are provided in equations. 9-11. Survivorship to age a was given by: (5) ( ) () () Aafor 1 Aafor 1 1 = ? = <= ? ? ? ? ? M a M a a M a e e L eL where M is the instantaneous natural mortality and A represents the age of the plus group. The maturity-at-age function was represented as a simple logistic function: (6) ))(( 1 1 lhag a e m ??? + = where m a is the fraction of age a individuals that are sexually mature, g is the shape variable and lh is the age of 50% maturity. Growth was assumed to follow the von Bertalanffy equation, given by: (7) )1( )( 0 tak a ell ?? ? ?= where l ? is the asymptotic length, k is the annual rate at which l ? is reached (the Brody growth coefficient) and t o is the theoretical age at zero length. An allometric relationship was used to estimate weight at age: (8) w a = ??l ? where ? and are parameters of the weight to length relationship. Given the survivorship and maturity at age, the scaled net reproductive rate was given by: (9) a a ae mL ? =? . 158 The incidence function (Eq 9) that represents the scaled net reproductive rate over the life-span of an individual, also provided the foundation to calculate the total reproductive effort in a population with a stable age distribution simply by multiplying ? e by the initial number of recruits (or in this case the leading parameter R o ) so that S o = R o ? e . The next step in deriving the parameters for the Ricker function was to note that two points on the stock recruitment curve have been defined corresponding to (0,0) and (S o ,R o ). Accordingly, alpha was defined as: (10) o o S R K=? and noting that eoo RS ?= this reduces to e K ?? /= where K is the recruitment compensation ratio, or the relative improvement in juvenile survival rate at low spawner abundance. Given initial estimates of K and R 0 and ignoring the process error term, i the density-dependent parameter in the stock recruitment function (?) was estimated as: (11) () ?? ? ? for ngsubstituti when /ln toreduces which ln o o oo o o S R KSKS S R = ? ? ? ? ? ? ? ? ?= . The leading parameters were used to initiate the population in model runs. Subsequently, the population was projected forward in time (eqs. 3 & 4). Using the projected population matrix a time series of predicted survey CPUE and commercial catch were calculated. Best fitting model parameters were determined by maximizing the likelihood of the observed data given estimated parameters for observation and process error models described below. 159 Observation models Three separate observation sub-models were constructed from predictions based on the population dynamics model. The first observation model predicted the commercial CPUE (conditioned on observed catch) assuming lognormal errors. The second observation model predicted the fishery-independent survey abundance information and also assumed a lognormal error distribution. The third observation model focused on the size composition data from fishery-independent trawl surveys. The approach assumed that sampling a fish of any size and any age was a random draw from a multinomial distribution, where the probability of sampling an age a fish was determined by the proportions-at-age in the population. Length frequency analyses Survey data represented a time series of abundance at length. To provide an equivalent time series from the model, the abundance at length was estimated from model predicted abundances at age using a multinomial distribution assuming normal errors. The distribution was chosen to capture the 20 age classes of winter skate. The likelihood function took the following generalized form: (12) ??? ? x pppp ....)( 21 ??= where p x is the predicted probability of capturing an individual at length x, and ? is the number of individuals observed at length x, taken from survey data. The observation model compared the observed survey length distributions to length distributions predicted by the population dynamics model. Note that this approach used the population dynamics model to predict a length distribution that can 160 be divided into an arbitrary number of length bins. Changing the number of length bins does not alter the underlying age structure in the model. The observation model multiplied the proportion at age in the population dynamics model by the vulnerabilities at age, based on survey selectivities, and by the predicted probability of observing a skate of length x given it?s age a to yield a relative length distribution which was used subsequently in the likelihood calculation and compared to the predicted numbers at age. The probability density function for an individual skate in the length interval x-d to x+d for a given age a is given in equation 13. This calculation is based on growth rates (k) and mean length-at-age calculated from the von Bertalanffy model developed in previous dissertation chapters. Assuming variation in length-at-age was normally distributed, the probability of fish of age a falling in a length interval x-d to x+d was estimated as: (13) ? + ? ? ? = dx dx Lx a t a ta eaxP 2 2 , 2 )( 2 1 ? ?? where x is the mid-point of the length interval (Fournier, 1991; Martell, in preparation). Standard deviation in length-at-age was a function of age and the growth coefficient k: (14) }) 1 1 {21( 1 1max 1 2 ? ? ? ? +? = T a p p a e ? ?? where )( k ep ? = . The parameters ? 1 and ? 2 explain how the standard deviation of length changes with age (Fournier, 1991; Martell, in preparation). The parameter ? 1 161 is the standard deviation at age 1 and positive values of ? 2 represents how fast ? increases with increasing age. The following equation estimated the probability of sampling a fish of length x from the population: (15) ?? ? == = = x x T a ax a T a x x axPs axPs p 11 1 max max ? ? where s x is selectivity to the survey, P(x|a) is the probability of fish of age a in length interval x, and ? a is the proportion of age a individuals in the population in year t. Size selectivity to the sampling gear was predicted using a logistic function: (16) ))2(2( 1 1 lhlg x e s ??? + = where g2 is a shape parameter, lh2 is the length at 50% selectivity. The two parameters (g2 and lh2) were estimated outside the model. Likelihoods The overall objective function that was minimized was composed of up to five different observation or process error terms. The number of components used depended on the model configuration considered. Three of the objective functions were utilized on all model runs. The first fitted observed abundances based on commercial (CPUE) and survey CPUE based on scaled z scores. The second fitted observed and predicted catch. In the third the observed length frequencies were fit to the predicted probability of a fish of age a in length interval x. Two additional 162 objective functions were used in some model runs: one included a prior on K and the second allowed process error (wz) in the recruitment function. Estimated parameters The five main parameters estimated in this model are K the maximum juvenile survival rate as the population approaches zero, R 0 the unfished recruitment level (population carrying capacity), B 0 the proportion of the virgin population size at the start of the time series, and ? 1 and ? 2 the standard deviations of the length frequency data. Also estimated are ? and ? the density-independent and -dependent parameters of the Ricker stock recruitment relationship, respectively. The model also calculates annual fishing mortality rate F t for each year based on the values of parameters estimated in the objective function minimization. When the model was fit with process error in recruitment, one recruitment error parameter was added for every year. Trial runs and ?best? model selection Base model runs had no process error and parameters were unconstrained with the exception of an informative prior on K. In a meta-analysis, Meyers et al. (1999) found that K (his? ~ ) ranged from 1-7 in most fishes. Preliminary model results indicated that values of K less than 5 resulted in lower model performance. Unconstrained models tended to have too large a value for K. In order to achieve robust results, a moderate prior for K = 5 was set. Six trial model runs were performed on both spring and fall data that each represented the length frequency data 163 in different sized bins: 1, 3, 6, 9,12 cm bins and a 2 bin model where juveniles (< 75 cm (TL)) and adults (> 75+ cm (TL)) were aggregated into two groups in a fashion similar to a delay-difference model (D.D). In subsequent model runs, models will be referred to by their bin size, or by D.D. Five criteria were used to judge ?best? model fits. Criterion 1 was the maximum likelihoods estimates. We did not use penalized information criteria, i.e., Aikaike?s Information criteria, to assess model performance as the number of parameters was equal in all models except for the one that included recruitment process error. Criteria 2, 3 and 4 consist of the sum of values for correlations, slopes, and intercepts from the relationships of predicted vs observed biomass (CPUE spring) and commercial CPUE vs survey CPUE. The logic behind criteria 2, 3 and 4 was that the condition when both slopes = 1 (combined sum = 2) and each intercept = 0 (combined sum = 0) would indicate predicted values match perfectly with observed values. Thus values for criteria 2, 3, and 4 closer to a "perfect? fit (i.e., =2) indicate a ?better? fitting model. The fifth criterion used was the summation of ? 1 and ? 2 where lower values would indicate better fitting models. Each criterion was given a ?score? of 1-5 with 1 being the best of the five models and 5 being the worst of the five models. These scores were summed and the lowest total score indicated the ?best? fitting model. Compensatory recruitment A process error model was used to test the hypothesis that the rate of increase in population was the result of ?expected? recruitment levels. The ?best? fitting model was run with process error in recruitment (wz) bounded between -10 and 10. 164 The model was run assuming base run conditions with a prior of K of 5 and other parameters unconstrained. This model structure has an additional 41 parameters compared to base runs. 165 Results Base model runs and model selection Six base model runs were completed for the fall survey for models that represented the population size structure in 1, 3, 6, 9 and 12 cm bin sizes and an additional formulation that mimicked a delay-difference structure (Table 3). Considering all models examined, initial population sizes of winter skate varied substantially. The smallest initial population size was less then 1 % of the estimated virgin population size. The largest initial population size estimated was 26% of the estimated virgin population size. The virgin recruitment (R o ) ranged from 203,260,710 to 221. The ? and ? parameters of the Ricker stock-recruitment relationship ranged from 18.19 ? 22.60 and from 3.02 x 10 -8 ? 2.61 x 10 -1 respectively. Estimated values of K varied from 3.29 ? 5.50. Based on overall rankings and model structure, models 1 and 9 were chosen to provide two contrasting model structures for further exploration. Model 1 Model 1 represented the population as comprised of 1 cm length classes. The model included four observation models in its objective function, but did not include any process error in its stock and recruitment function. Non-linear optimization algorithms produced parameter values that best predicted the fall fishery-independent survey and the commercial CPUE time series. This model was chosen for further exploration as it provided the most highly resolved representation of the population size structure. Also, the overall likelihood of the model and the summed r 2 and slope 166 values of the relationships between observed and predicted values were the 2 nd highest of all the alternative models considered. The optimum solution had a steepness value of K = 4.42 and unfished recruitment R 0 = 203,260,710 (Table 3). This model also showed relatively low standard deviations for length at age (? 1 = 38.64, the standard deviation of the first age interval, and ? 2 = 3.31, the rate at with standard deviation increases with age). As with all of the model runs, model predictions tended to underestimate survey and commercial biomass as indicated by the slopes of the fits between predicted and observed values (Table 3). The time series of predicted biomass of model 1 exhibited similar behavior to that for the observed survey CPUE (Figure 4). However, the predicted survey CPUE does not fully reflect the rate of increase in the early 1980?s observed in the survey data. Specifically, residual values (observed ? predicted) were consistently positive in the 1982 - 1988. Moreover, predicted biomass declines less quickly than the observed survey biomass after reaching its peak in the 1980?s, so that the model overestimates survey biomass for the period 1989 - 1997. The predicted commercial CPUE also did not match the observed commercial CPUE trend well for years after the mid 1980?s (Figure 5). In particular, the large increase in the commercial CPUE that began in the 1980?s and continued until 1998 was not captured by the model. There was a strong (r = -0.99) negative correlation between estimates of B 0 , the proportion of virgin biomass at the beginning of the model and R 0 , the virgin recruitment in this model. Indeed this was a consistent problem with all model runs for which correlations varied between -0.7 to -0.99. The correlation between R 0 and 167 B 0 indicates a trade-off in model fitting between the initial population size and the equilibrium population size. This results in model runs in which the initial population size as either very large or having been substantially depleted. To explore this pattern, I ran additional simulations in which the possible values of B 0 were bounded. In these runs the model consistently converged to the lower bounds of B 0 . As the lower bound for B 0 was increased, the pattern of correlation between R 0 shifted from one with B 0 to one with K (Table 4). When the model is run with a B 0 with a lower bound of 0.5 parameter correlations disappear all together (Table 4). However, under these conditions, the model loses its ability to reproduce the pattern in the observed survey biomass (Figure 6). The biomass estimate is greatly increased as B 0 is increased and the predicted biomass is well below the observed survey CPUE during the 1980?s (Figure 6). While a higher B 0 allows the model to predict larger overall biomass it restricts the ability of the predicted biomass to fit the survey CPUE. For model 1, normal approximations of the likelihood profiles of parameter estimates were estimated for B 0 , R 0 and K. Likelihood profiles indicate that all parameters are estimated relatively well with 95% confidence intervals of likelihoods for K of 2.38 to 6.39, for R 0 0.91 to 36.6. Likelihood profiles estimates for B 0 were difficult to estimate as the model estimated B 0 was near zero. The normalized likelihoods indicated that the 95% confidence interval of B 0 was -0.00044 to 0.00047. The estimated parameters of the Ricker equation for model 1 lead to a stock- recruitment relationship that was effectively linear over the range of adult biomasses observed. The time series of relative Fs predicted by model 1 are shown in Figure 7. The model suggests that fishing mortality increased in the beginning of the time 168 series, then declined in the late 1980?s, only to increase over the last 15 years of the survey. Model 9 Model 9 was similar to Model 1 in that it comprised four observation models in its objective function but did not include any process error in its stock and recruitment function. However, Model 9 differed from the first model in that it represented the population as comprising 9 cm length classes. Non-linear optimization algorithms produced parameter values that best predicted the fall fishery-independent survey and the commercial CPUE time series. This second model chosen ?scored? 20, ranking among the worst overall in terms of fit. The model was chosen for further analysis as it represented a more aggregated analysis of survey length frequencies and has reasonable values for all parameters (Table 3). Leading parameters for Model 9 were K of 5.46, and unfished recruitment of R 0 = 43,914 (Figure 8 & 9). B 0 was estimated as 0.24 indicating a starting population at one quarter of the equilibrium level, a much higher level than Model 1. Parameter correlation between R 0 and B 0 was 0.78 and between K and B 0 was 0.33. The 95 % confidence interval for K was 2.87 to 7.95 and for R 0 was 9.76 to 11.6. While the confidence estimates were similar for K they were smaller for R 0 when compared to model 1. The model failed to estimate a confidence interval (95%) for B 0 . The survey CPUE predicted by the model exhibited a gradual increase from 1970 ? 1990 (Fig 8). However, while this general pattern is characteristic of the observed survey CPUE, the model captured neither the rapid rate of increase nor the 169 absolute magnitude of the increase in survey CPUE observed in the 1980?s. Similarly, model predictions of commercial CPUE did not adequately describe the observed commercial CPUE (Fig. 9). Estimated parameters of the Ricker recruitment function for this model (Table 3) lead to a stock-recruitment relationship that exhibited strong density-dependence over the range of adult biomasses observed. Peak recruitment levels are estimated to occur at relevant intermediate stock biomasses. The time series of Fs predicted by model 9 are shown in Figure 7. The model suggested similar trends in fishing mortality that were reported for model 1 with peaks in fishing mortality during the 1960?s and 1970?s, followed by a decrease in the 1980?s, and a subsequent increase over the last 15 years of the time series to peak levels. Recruitment anomalies To address whether the rate of increase in population was the result of ?expected? recruitment levels a model that included process errors in the stock and recruitment function was fit. The inclusion of this additional term in the objective function allowed the recruitment predicted in the model to deviate from the underlying Ricker stock and recruitment function, thereby estimating a recruitment anomaly for each year in the model. Process error was assessed for recruitment on fall survey CPUE data. Optimizing model 1 with process error indicated that estimated recruitment anomalies could play a potential role in the population dynamics of winter skate (Figure 10). A review of yearly estimates of the recruitment anomaly (wz) estimated by the model indicated that wz were large and positive 170 throughout the 1970?s. Recruitment anomalies were smaller and variable in sign from the 1980?s ? present. Model generalizations A trade-off between K and R 0 was common in model parameter estimates so that when K was low R 0 was often high. This was not always the case, but it suggests that there was a trade-off between maximum juvenile survival and equilibrium recruitment. There was also a trade-off between either K (or R 0 ) and B 0 , indicating there was another trade-off between the starting biomass and estimates of maximum juvenile survival and/or equilibrium recruitment. When models were initialized with a B 0 near zero, the model tracked observed survey CPUE well, whereas if B 0 was 0.2 or more, the model had difficulty following survey CPUE trends. Overall the models indicated that in 1963 the winter skate population was likely at low abundance and the rapid increase in the 1980?s probably was less drastic than survey data indicates. 171 Discussion Without estimates of the species? vital rates and the development of assessment models, understanding the population dynamics of winter skate has remained elusive. Here I utilized recently estimated vital rates for winter skate to successfully fit an age-structured model to survey CPUE data and catch data for winter skate. A suite of models was developed that differed in the level of resolution with which they resolved the size structure in the population and in how recruitment variability was included. However, none of the models were able to replicate both the rate and magnitude of the increase in biomass in the 1980?s that is indicated by survey trends. Additionally, all the model structures failed to fit commercial CPUE trends. An age-structured population model was at the core of the developed model. The population model was the same in all formulations examined. Non-linear optimization techniques adjusted parameter values to minimize differences between predicted survey CPUE and predicted commercial CPUE and their equivalent observed CPUE time series in all formulations. However, model formulations differed in the degree to which the population model resolved the size structure in the survey CPUE time series. In some formulations, the population model predicted the size distribution of the survey CPUE time series in 1 cm size classes, in other simulations the size structure was less resolved. Thus in the comparison the alternative model formulations, it is important to remember that the underlying population model structure remained constant and involved the same number of parameters. 172 All models predicted an increase in winter skate population biomass indicated in the surveys in the 1980?s. However, the models differed in the degree to which they could resolve the pattern of increase observed in the survey. Models with highly resolved length structures (Model 1) fit the increase in biomass in the 1980?s better than models with more aggregated structure (Model 9). Model 1 captured the magnitude but not the rate of increase in the 1980?s whereas Model 9 captured neither the magnitude nor the rate of increase. Three inferences can be drawn from the discrepancy between observed and predicted abundances: (i) the model was mis- specified, (ii) the survey does not reflect the dynamics of a closed population, or (iii) empirical estimates of vital rates are inappropriate. Below, I will deal with each possibility in turn. Survey data indicated an emergence of high biomass levels for both adults and juveniles in the early 1980?s following a period of very low adult biomass. Given the age and growth estimates of winter skate (Chapter 4; Sulikowski et al.; 2003), this scenario created a difficult challenge for fitting the model to survey data. It is also important to note that all models indicated lower population biomasses in the early 1960?s than indicated in the survey time series. The failure to fit to the early 1960?s observed data resulted because the model was initiated with a population at a stable age distribution whereas the survey data indicated a relative lack of adults in the survey and high levels of biomass of young juveniles in the early 1960?s. Thus, the reproductive capacity of the population in the early 1960?s was low while 1-4 year old skate abundance was high. 173 An important assumption when interpreting the survey time series is that it represents a closed population. If the population the survey indexed was not closed, it would be unreasonable to expect a population model, which is explicitly closed, to replicate the observed abundances. Changes in the distribution of the winter skate population such as seasonal or long term migrations and range restriction would influence the portion of the population effectively sampled, thereby altering the effective catchability of the survey. However, here I used data for the entire western Atlantic south of Canadian waters and there are no data available to suggest a large scale shift in winter skate?s population structure. Yet, winter skate is known to make onshore-offshore seasonal migrations and may use Georges Bank as a primary spawning ground (Chapter 4). For example, all genetically-identified winter skate hatchlings were caught on Georges Bank and most were caught during the fall (4 of 5 ? Alvarado-Bremer et al. 2004). A review of migration patterns of the species shows a juvenile onshore-offshore migration (Appendix C). Adults are spread over the species range in the spring while in the fall they are clustered on Georges Bank (Appendix C). If Georges Bank is a focus for mating or spawning it could, in part, explain the high variability in survey estimates if reproductive patterns follow seasonal or environmental cues. Unfortunately, our understanding of the species? reproductive biology does not allow for a conclusive solution to these issues. Overall, the different model formulations all yielded realistic parameter estimates. The model with the highest likelihood was model 1, which provided a maximum population growth rate of K = 4.42 and a population capacity of R 0 = 203 million. Both values are high considering the species? life history and the reported 174 range in teleost species for K of 1 to 7 (Myers et al., 1999). The R 0 value is much larger than the peak in abundance in the survey CPUE series. R 0 is sensitive to B 0 and in this case is likely overestimated with the peak in predicted biomass of 83,739, which is less the 1 % of equilibrium biomass. The more aggregated model 9 resulted in a K = 5.46 and R 0 = 43, 914 and peak abundance higher than equilibrium biomass. Together, both models indicated that the population has been below carrying capacity for the duration of the time series. Model results are dependent on the quality of vital rate estimates used to parameterize the model. Error in parameters for age of maturity (t mat ), life span (t max ) and the Brody growth coefficient (k) have the potential to change model performance substantially. Values used in the model are similar to previous estimates for winter skate reported by Sulikowski et al. (2003). Parameter estimates used in the model (Chapter 4) indicated a slightly faster growth rate (k = 0.64 vs 0.9) and smaller asymptotic size (l ? = 131 cm vs. 122 cm) than estimated in Sulikowski et al.?s study. Additionally, the model structure presented in this study assumes that all vital rates are constant for the time series. However, winter skate has experienced large changes in fishing effort and directed harvest during the time series (Murawski and Fogarty, 1998) which may have lead to variation in vital rates. Clearly, it is possible that vital rates differed in the 1960?s when the population was at a low biomass level compared to the large population size of the 1980?s. However, there are no historical estimates of vital rates for winter skate. Fecundity was not directly used in the model but was instead formulated using a Ricker stock recruitment function. Model 1 resulted in essentially a linear 175 recruitment function while model 9 showed a higher degree of compensatory dynamics of the relevant range of adult stock sizes. However, for both formulations, the rate of population increase at small population sizes estimated by the model were in line with estimates derived for estimates for teleosts generally (Myers et al., 1999). While very little is known of recruitment in elasmobranchs, there is theoretical evidence to suggest that recruitment is a linear function of stock size in winter skate as estimated in model 1 (Cortes, 2004). More generally, elasmobranch life histories are characterized by low fecundity and high offspring survival (Chapter 3; Frisk et al. 2001), suggesting that recruitment relationships should be strong and closely reflect adult biomass (Cortes, 2004; Smith et al., 1998). However, winter skate has withstood relatively high exploitation levels implying that recruitment relationships likely have compensatory behavior. Estimated recruitment patterns in model 9 may suggest more density effects than would be expected for an elasmobranch species. Clearly, this is an area in need of continued research effort. All models assume a constant catchability for the time series. Survey data indicates large shifts in late juvenile and adult biomass during several periods in the time series. If cacthability increased with size in winter skate, the large changes in adult biomass during the time series would be easier to explain. There is indeed evidence within the survey time series for substantial changes in the size structure of the winter skate population over the time series. As already noted above, the population in the early 1960?s was characterized by a proportionately high abundance of juveniles, where more recently a full age (and hence size) structure has developed in the population. 176 All models failed to adequately replicate commercial CPUE time series. Estimates of catch and effort data are areas of high uncertainty. Reported aggregated mixed skate landings from Cape Hatteras north to Canadian waters were used as catch estimates which may not reflect the historical catch in winter skate. Additionally, inconsistencies and under-reporting of landings and discards add to uncertainty of catch estimates. Use of aggregate landings also assumes that the distribution of skate species represented in the catch has not changed over time. Given knowledge of exploitation in the groundfishery of the western Atlantic, the trends observed in the skate landings seem to follow the pattern of high effort from the foreign fleet in the 1960?s to the mid 1970?s. The exclusion of the foreign fleet lead to an initial decline in effort. However, the subsequent expansion of the domestic fleet lead to an increase in effort from the mid 1980?s to the present. Thus, reported landings in the 1960?s and 1970?s, which occurred during peaks in effort of the groundfishery, may be underestimates as a result of foreign fleets keeping a greater proportion of skates than the domestic fleet at the time. During the period of high activity of the foreign fleets, reported effort for the groundfishery was probably a good indicator of skate effort. However, in the 1980?s and 1990?s, when the domestic fleet began keeping a greater number of skates, the effort reported for the groundfishery probably was an underestimate of effort directed at skates. This is evident in the very high commercial CPUE in the 1980?s and 1990?s. During this period the domestic fleet was increasingly landing skate for the market. It is possible that in years prior to the mid 1980?s skates were caught mainly as by-catch, with some discard survival, while in later years they were kept for the market. Thus, the 177 estimates of effort of the multi-species groundfishery may not follow the changes appropriate for winter skate. All models estimated similar values for fishing mortality (F) during the time series and suggest that fishing mortality was highest in the late 1960?s to the mid 1970?s and during the 1990?s. During the 1970?s estimates of fishing mortality were as high as F = 0.2 and more recently have reached an F = 0.3. The period of lowest winter skate biomass corresponded to estimated fishing mortality ranging F = 0.05- 0.2. Currently, the estimated fishing mortality is higher than the 1970?s raising concern that current fishing mortality may not be sustainable. The model that incorporated process errors in recruitment indicated a period of higher than expected recruitment throughout the 1970?s. Could competitive release from the declines in other groundfish species have provided a competitive advantage leading to increased recruitment rates in winter skate that help explain the sudden increase in biomass in the late 1970?s and 1980?s? It is plausible that inter- and intra-species specific density effects have a role in recruitment in winter skate. Inclusion of process error did allow the model to capture the rate but not the magnitude of the increase in biomass in the 1980?s as the positive recruitment anomalies in the 1970?s and allowed recruitment to increase as much as 3.3 times greater than observed adult biomass. If survey data are accurate, estimates of recruitment process error provides indirect evidence that reproductive potential and/or survival of 0-age skates may have been higher than observed in the 1970?s. Three potential hypotheses can be put forward to explain these patterns. First, female winter skate may have benefited from reduced competition with other fishes during 178 this period. Accordingly, female skates may have been in better condition during the late 1970?s than in other times. In this scenario, an equal spawning biomass might have produced more offspring in the 1970?s, driving the observed increase in abundance, than in subsequent years. However, increased egg production as a compensatory response is unlikely as a result of body cavity volume limitation in elasmobranchs (Holden, 1972). The other two hypotheses would invoke increased survival of juvenile skates. One hypothesis would suggest increased survival of juvenile skates as a result of increased female provisioning or fitness. A second hypothesis would ascribe the increased juvenile survival to reduced predation pressure from other fishes. Irrespective of whether these hypotheses are valid, my results do suggest that if process error in recruitment is accounted for winter skate population could increase at the rate suggested by the NMFS fall survey. However, the model still underestimates the magnitude of the increase in the 1970?s to the mid- 1980?s. Until more field-derived data are estimated, the importance and effect of compensatory recruitment remains unanswered, but potentially important. In summary, I have developed a suite of models to predict the dynamics of winter skate in the northwest Atlantic. All model formulations predicted that the winter skate populations in the 1960?s was substantially below current carrying capacity. This suggests that the skate population could have been larger at earlier times if current carrying capacity is a reliable estimate of former carrying capacity. All models also indicated that winter skate abundance increased in the 1980?s to a peak in mid 1980?s, and suggest population declines or leveling off since mid the 1990?s. However, the models capture increase in biomass, but do not fully replicate 179 rate of biomass increase, suggesting either model mis-specification, biases in survey estimates, or inaccuracies in current estimates of winter skate vital rates. The role of additional factors such as environmental or biological interactions in determining the rate and magnitude of increase of winter skate remain poorly described. The model presented here does not address these issues. However, it is clear that winter skate, like other skates, is susceptible to population decline (Holden, 1973; Brander, 1981; Dulvy et al., 2000; Casey and Myers, 1998; Walker and Hislop, 1998; Stevens et al., 2000; Frisk et al., 2001; Frisk et al., 2002). Thus, in depth analyses of survey trends and the population structure of winter skate are needed to fully understand the dynamics of winter skate in the western Atlantic. 180 Table 1. Input data, definitions of parameters and symbols used in model. Definitions that are italicized provide model input values. Symbol Definition Subscripts t mat age at maturity t max maximum age a age of fish t year in survey x length X index for length interval Population dynamics B biomass N number of individuals B 0 the proportion of the virgin population size at the start of the time series (ranges from 0-1) v vulnerability to fishery (0,0,0.5,1,1?T max ) w weight F fishing mortality C catch R recruitment wz variable for recruitment process error ? density-independent parameter of the Ricker model ? density-dependent parameter of the Ricker model S spawning stock biomass M natural mortality (0.22, based on Heonig?s method) Initial calculations L survival m maturity l length lh age at %50 maturity (12.5) g shape parameter of the maturity ogive (logistic)(3) s selectivity to survey gear %50 selectivity to the survey gear (Fall=56;Spring=47) shape parameter of the selectivity to survey gear (Fall=0.078;Spring=0.1) l ? asymptotic length (111.43 cm) k growth coefficient (0.09/year) t 0 theoretical age at hatch (-1.80) ? intercept of weight-length power curve (0.0000019) ? exponent of weight-length power curve (3.33) ? ? ?fecundity incidence? 0 ? ? unfished eggs per recruit 181 K maximum increase in juvenile survival as a population approaches zero R 0 average recruitment for a unfished population Length frequency analysis ? 1 standard deviation at age 1 ? 2 standard deviation as age approaches T max ? proportion at age a in year t 182 Table 2. National Marine Fishery Service?s annual fall bottom trawl survey, reported commercial landings and groundfish effort on southern New England and Georges Bank. Fall survey data is from offshore strata (1-40) for southern New England, Georges Bank and the Gulf of Maine. Year Fall survey (n/tow) Landings (mt) Effort (days fished) 1963 3.01 33 45 1964 2.33 4083 55 1965 2.16 2363 75 1966 1.94 2844 90 1967 1.12 4898 85 1968 0.95 6483 82 1969 0.76 9462 99 1970 1.55 4128 66 1971 0.77 5905 98 1972 2.46 8823 120 1973 3.15 7963 95 1974 1.18 3651 89 1975 0.53 3968 125 1976 0.82 1212 60 1977 1.98 1418 31 1978 1.61 1353 29 1979 2.05 1360 27 1980 2.21 1581 35 1981 1.95 847 37 1982 3.63 878 39 1983 4.23 3603 37 1984 5.22 4157 39 1985 3.62 3984 41 1986 5.24 4253 38 1987 6.00 5078 39 1988 4.46 7264 39 1989 3.17 6717 39 1990 3.94 11403 40 1991 4.27 11332 39 1992 2.87 12525 39 1993 2.10 12904 39.25 1994 2.60 8829 39.25 1995 2.25 7222 39.25 1996 2.09 14226 39.25 1997 2.25 10952 39.25 1998 3.02 16936 39.25 183 Table 3. A. Results of model runs including parameters estimates K, R 0 , B 0 , ? 1 and ? 2 . Also shown are ? and ? of the Ricker stock recruitment function. Also shown is the likelihood for each model run and the results of regression fits to predicted biomass (cpue) vs survey cpue and a second regression between predicted biomass and commercial cpue. B. Provided are the ?scores? for selection of ?best? fitting models for each criterion and the final ?score?. A.1. A.2. B. Model Combined r 2 Combined slope Combined intercept Likelihood Combined (l1 & l2 ) 1 2 3 4 5 Score 1 0.94 3.87 67673 323.1 41.9 6 3 1 1 5 16 3 0.98 3.92 132947 271.5 37.1 2 4 4 4 4 18 6 0.97 3.99 139322 252.13 32.9 3 5 5 3 3 18 9 0.96 4.03 144152 247.9 32.3 4 6 6 2 2 20 12 0.96 0.54 106735 273.2 17.9 5 1 2 5 1 14 D.D. 0.98 3.71 124553 125.5 87352.7 1 2 3 6 6 18 Model Bin size K LN(Ro) B 0 ? 1 ? 2 ? ? 1 1 cm 4.42 19.13 0.00 38.64 3.31 3.02E-08 18.19 3 3 cm 5.50 10.75 0.21 33.83 3.27 1.50E-04 22.60 6 6 cm 5.49 10.71 0.23 28.94 3.99 1.56E-04 22.59 9 9 cm 5.46 10.69 0.24 27.86 4.40 0.000158 22.46 12 12 cm 3.29 5.48 0.26 14.00 3.93 2.61E-01 22.5 D.D. D.D. 5.33 10.94 0.17 87339 13.68 0.000122 21.94 Likelihood r 2 (survey) slope intercept r 2 (comer.) slope intercept 323.1 0.42 0.77 11956 0.52 3.1 55717 271.5 0.28 0.83 18554 0.70 3.1 114393 252.13 0.29 0.87 16871 0.68 3.13 122451 247.9 0.29 0.88 16242 0.67 3.10 127910 273.2 0.30 0.32 51421 0.66 0.21 55314 125.5 0.25 0.75 22978 0.73 2.96 101575 184 Table 4. Correlation to R 0 for varying lower bounds of B 0 . Note that the Hessian failed with B 0 lower bound set at 0.4 so the results are not reported. B 0 lower bound 0.0 0.1 0.2 0.3 0.5 K -0.01 -0.95 -0.88 -0.68 0.00 B 0 -0.99 -0.005 0.00 0.00 0.00 185 B. A. Figure 1. Survey CPUE are shown for the fall for the Gulf of Maine, Georges Bank, and southern New England. 0 2 4 6 8 10 12 14 16 18 1962 1966 1970 1974 1978 1982 1986 1990 1994 1998 2002 G e o r ge s B a nk ( num b e r / t o w ) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 G u l f of M a i n e ( num be r / t o w ) Georges Bank Gulf of Maine 0 2 4 6 8 10 12 14 16 1962 1972 1982 1992 2002 In s h o r e (n u m b e r/t o w ) 0 0.5 1 1.5 2 2.5 3 Of f s hor e ( numbe r / t o w) SNE-Inshore SNE-Offshore 186 0 0.5 1 1.5 2 2.5 3 3.5 4 1963 1968 1973 1978 1983 1988 1993 1998 2003 CPU E (number/tow ) Fall Figure 2. Survey CPUE is shown for the fall survey. The fall CPUE includes offshore stations 1-40. 187 30 36 42 48 54 60 66 72 78 84 90 96 102 108+ 196 3 1966 19 69 1972 19 75 1978 19 81 1984 19 87 1990 1 993 1996 1 999 200 2 Length (cm) Year 0.95-1 0.9-0.95 0.85-0.9 0.8-0.85 0.75-0.8 0.7-0.75 0.65-0.7 0.6-0.65 0.55-0.6 0.5-0.55 0.45-0.5 0.4-0.45 0.35-0.4 0.3-0.35 0.25-0.3 0.2-0.25 0.15-0.2 0.1-0.15 0.05-0.1 0-0.05 Figure 3. The length frequencies in three centimeter intervals of the fall survey for offshore regions from 1963 to the present. Year is on the x-axis, length is on the first y-axis and survey CPUE is on the second y-axis. 188 0.00E+00 2.00E+04 4.00E+04 6.00E+04 8.00E+04 1.00E+05 1.20E+05 1960 1965 1970 1975 1980 1985 1990 1995 2000 Biomass Predicted biomass Survey cpue Figure 4. Predicted biomass and survey CPUE for model 1 for 1963 to 1998. Note that both biomass and survey CPUE are scaled values. 189 0 50000 100000 150000 200000 250000 300000 350000 1960 1965 1970 1975 1980 1985 1990 1995 2000 Biomass Predicted biomass Comercial cpue Figure 5. Predicted biomass and commercial CPUE for model 1 for 1963 to 1998. Note that both biomass and commercial CPUE are scaled values. 190 0.00E+00 1.00E+15 2.00E+15 3.00E+15 4.00E+15 5.00E+15 6.00E+15 7.00E+15 8.00E+15 9.00E+15 1960 1965 1970 1975 1980 1985 1990 1995 2000 Biomass Predicted biomass Survey cpue Figure 6. Model 1 results assuming a lower bound for B 0 of 0.5 for 1963 to 1998. Note that both biomass and survey CPUE are scaled values. 191 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 1960 1965 1970 1975 1980 1985 1990 1995 2000 Fishing mor t al i t y (F) Model 1 Model 9 Figure 7. Estimated fishing mortality for 1963 to 1998. Results for model 1 and model 9 are shown in the figure. 192 0.00E+00 5.00E+04 1.00E+05 1.50E+05 2.00E+05 2.50E+05 1960 1965 1970 1975 1980 1985 1990 1995 2000 Biomass Predicted biomass Survey cpue Figure 8. Predicted biomass and survey CPUE for model 9 for 1963 to 1998. Note that biomass and survey values are scaled. 193 0.00E+00 5.00E+04 1.00E+05 1.50E+05 2.00E+05 2.50E+05 3.00E+05 3.50E+05 4.00E+05 4.50E+05 1960 1965 1970 1975 1980 1985 1990 1995 2000 Biomass Predicted biomass Commercial cpue Figure 9. Predicted biomass and commercial CPUE for model 9 for 1963 to 1998. Note that biomass and commercial CPUE are scaled values. 194 0 50000 100000 150000 200000 250000 1960 1965 1970 1975 1980 1985 1990 1995 2000 Biomass -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Anomalies recruitment (wz) Biomass Survey CPUE Anomalies recruitment wz Figure 10. Results of model 1 with process error in recruitment for 1963 to 1998. Scaled biomass, survey CPUE and values of recruitment error (wz) are shown. 195 Chapter 7 SUMMARY Sharks, skates and rays are an important ecological and economic resource (NEFSC, 1999; Link et al., 2002). Yet for many species there is a paucity of information on vital rates (Frisk et al., 2001). This paucity of knowledge of elasmobranch life history has been highlighted by the decline of several species in the western Atlantic (NEFSC, 1999). Barndoor skate, thorny skate and winter skate all have shown population declines and in some cases drastic declines (Casey and Myers, 1998; NEFSC, 1999; DFO, 2002). However, managers are left with few conservation options with the basic knowledge of skate species biology such as maturation and growth rates largely unknown. The goals of my dissertation aimed at improving the conservation of elasmobranchs through comparative life history analyses and species-specific approaches. In the comparative approach I analyzed the life history patterns of elasmobranchs in order to better understand evolutionary relationships and conservation of sharks, skates and rays. This was achieved in chapter 2 ?Life histories and vulnerability to exploitation of elasmobranchs: Inferences from elasticity, perturbation and phylogenetic analyses? and chapter 3 ?Ordination of evolved life history strategies in elasmobranch and teleost fishes: In search of common ground?. In the second part of my dissertation I apply the species-specific approach by estimating vital rates of age, growth, maturation, and fecundity for little skate and winter skate. Latitudinal patterns in age, growth and maturation were analyzed over the ranges of little skate and winter skate in chapter 4 ?Age, growth, and latitudinal 196 patterns of two Rajidae species in the northwestern Atlantic: little skate Leucoraja erinacea, and winter skate, Leucoraja ocellata? and chapter 5 ?Maturation, fecundity, latitudinal patterns and reproductive strategies of two Rajidae species in the northwestern Atlantic: little skate Leucoraja erinacea, and winter skate, Leucoraja ocellata?. Finally in chapter 6 ?An age-structured model of winter skate, Leucoraja ocellata, abundance in the western Atlantic: sustainability and uncertaint? I utilize estimated parameters for winter skate to build an age-structured model to analyze the species? exploitation and abundance trends. In chapter 2, I built empirical models based on phylogenic relationships, vital rates, and results from matrix analyses, (elasticity and sensitivity of elasticity), to provide insights into the population dynamics of elasmobranchs. My elasticity analyses provided evidence of a trade-off between survival and reproductive investment in elasmobranchs. Additionally, I found that for most elasmobranch species survival into the adult stage was important in regulating a species response to exploitation. Specifically these findings suggest conservation polices aimed at reducing juvenile mortality would have the greatest impact on overall population growth rate. Phylogenic analyses provided evidence that closely related species will respond to exploitation or conservation polices in a similar fashion. In chapter 3 I continued analyzing life history parameters in elasmobranch species. I expanded the Winnemiller and Rose (1992) ordination of evolved life history strategies for teleost species by adding elasmobranchs to the analysis. Interestingly, elasmobranch species, using the ordination based on teleost life histories, were placed in the periodic strategy. This was in contrast to the notion that 197 elasmobranchs are the panicle K-selective species (Cortes, 2004). However, when the ordination is re-estimated using the combined data sets for elasmobranchs and teleost species, they were placed in the equilibrium strategy due to their large body size, high parental investment, and large egg size. My results indicate that even though elasmobranchs have life histories similar to reptiles and perhaps mammals, they follow many of the same life history invariant relationships shown for a wide diversity of animal taxa including teleost species (Charnov, 1993; Frisk et al., 2001). In the second part of the dissertation, I analyzed the vital rates of little skate and winter skate along the eastern seaboard. Samples were collected in conjunction with the National Marine Fisheries Service?s annual bottom trawl and scallop surveys. In chapter 4 I found strong evidence for a latitudinal pattern in maximum size in little skate with individuals growing slower and reaching a larger size in northern portions of the species? range. However, no regional trends were found for winter skate. My results did not support a latitudinal increase in size of maturation in winter skate reported by McEachran et al. (1977). Perhaps winter skate?s population decline in the 1970?s has resulted in a reduction of the species? range and homogeneous expression of vital rates. Additionally in chapter 4 my data indicated that little skate is a moderately fast growing species with a k = 0.19 and maximum observed age of T max = 12.5. Theoretical estimates of maximum age indicate that longevity in little skate increases with latitude. Theoretical maximum age in little skate in the mid-Atlantic T max = 13, southern New England-Georges Bank T max = 14, and the Gulf of Maine T max = 15. 198 Winter skate is slower growing compared to little skate with a k = 0.07, longer lived T max = 20.5. Theoretical longevity in winter skate was estimated at T max = 35. In chapter 4 I estimated size and age of maturation, annual fecundity and net reproductive rate for little skate and winter skate. Maturation occurs in two levels, physiological, the presence of reservoir eggs, and functional, the presence of mature eggs. There was no difference in age of physiological and functional maturity in little skate (T mat = 7). However, in winter skate functional maturity (T mat = 12.5) occurred three years after physiological maturity (T mat = 9), indicating that previous estimates only considering physiological maturity underestimate age at maturity by three years. Latitudinal patterns were observed in little skate in size of maturation but not in age of maturation, while no regional patterns were observed in winter skate. Fecundity in little skate ranged from 21-57 eggs per year, and in winter skate 101-26 egg cases per year. From chapters 4 and 5, it is clear that the biology of little skate and winter skate represent different locations on the ?fast-slow? continuum of life histories in elasmobranchs. Little skate has a ?fast? life history characterized by a shorter life- span, earlier age of maturation and a higher growth rate compared to the ?slow? life history of winter skate, which is characterized by long life-span, slow growth and delayed maturation. Interestingly the ?slow? life history of winter skate coincides with higher fecundity, larger and more numerous offspring. In skates, increased fecundity may be a benefit to delaying maturation and obtaining a large body size. In chapter 6 parameters estimated in chapters 4 & 5 are used to parameterize an age-structured model for winter skate. My results provided reasonable estimates 199 of the maximum increase in juvenile survival as a population approaches zero (K) and average recruitment of the unfished population (R 0 ). Model predictions in biomass were not able to capture the rapid increase in survey cpue in the 1980?s and instead indicated a more gradual increase (NEFSC, 1999). Additionally, survey estimates of cpue were variable between regions and seasons. My results also indicate that effort for the grounfishery in recent years may not represent the realized effort winter skate is experiencing. Additional fishery dependant data is needed to better manage and understand the population dynamics of winter skate. 200 Appendix A. Model code for AD Model Builder for the age-structured model utilized in this manuscript. Note the model can be conditioned on catch or effort. However, I do not report any results of models conditioned on effort. //****************************************************** // Programmers: Mike Frisk & Steve Martell // Project Name: SkateBoard.tpl // Date: May 26th, 2004 // Version: 1.01 //******************************************************/ DATA_SECTION init_int ft_flag; init_int nyrs; init_int nage; init_number linf; init_number k; init_number to; init_number a; init_number b; init_number g; init_number lh; init_number lh2; init_number g2; init_vector effort(1,nyrs); init_vector obs_ct(1,nyrs); init_vector yt(1,nyrs); init_number lmin; init_number lstp; init_number nbins; init_matrix lf_data(1,nyrs,1,nbins); vector x(1,nbins); !!x.fill_seqadd(lmin+lstp,lstp); //!!cout<5.)ft(i)=5.; } if(!ft_flag) { ft(i)=q*effort(i); } dvariable st=sum(elem_prod(n(i),mat)); n(i+1,0)=so*st*exp(-beta*st)*exp(wt(i)); n(i+1)(1,nage)=++elem_prod(n(i)(0,nage-1),exp(-m-ft(i)*va(0,nage- 1))); } //cout< 0"<