BULLETIN OF MARINE SCIENCE. 89(4):000–000. 2013 http://dx.doi.org/10.5343/bms.2012.1041 1Bulletin of Marine Science Ξ2013ZosĞŶsƟĞů^ĐhooůoĨDĂriŶĞΘƚmosphĞriĐ^ĐiĞŶĐĞoĨ ƚhĞhŶiǀĞrsiƚLJoĨDiĂmi EVIDENCE OF EL NIÑO/LA NIÑA–SOUTHERN OSCILLATION VARIABILITY IN THE NEOGENE- PLEISTOCENE OF PANAMA REVEALED BY A NEW BRYOZOAN ASSEMBLAGE-BASED PROXY Beth Okamura, Aaron O’Dea, Paul Taylor, and Anna Taylor ABSTRACT Here we explore how fossil cheilostome bryozoans can demonstrate El Niño/La Niña–Southern Oscillation (ENSO) variability in ancient tropical environments of the tropical eastern Pacific and southwestern Caribbean when used collectively to produce frequency distributions of estimates of mean annual ranges in temperature (MARTs) via zooid-size MART analysis (zs-MART). The approach is based on linking variation in the sizes of constituent zooids in bryozoan colonies with the temperature regimes in which the zooids developed. The shapes of frequency distributions of zs-MART estimates from modern environments are consistent with known upwelling and non-upwelling environments and ENSO variability. Data from fossil colonies provide evidence that ENSO variability characterized Caribbean environments in what is now present-day Panama during the Miocene and Pliocene (thus supporting neither a permanent El Niño nor a permanent La Niña in the Pliocene), but not the Pleistocene after the Isthmus of Panama closed. Western Atlantic data provide further evidence for ENSO variability in the Pliocene of Florida but not Virginia. Our study shows the potential of bryozoan assemblages to infer variation in seasonal regimes thus encouraging the further development of this proxy for inferring ENSO variability. The ability to infer the quasi-periodic climate patterns expected from the El Niño/ La Niña–Southern Oscillation (ENSO) is of great importance if we are to better understand the wider impacts of present-day and future global climate regimes. Extreme weather and changes in ocean circulation linked with the El Niño cycle are responsible for collapses in South American fisheries (e.g., Arntz and Tarazona 1990), droughts and fires, monsoons and tropical cyclone activity in Australasia (Evans and Allan 2006, Nicholls and Lucas 2007), cycles of malaria in India and South America (Kovats et al. 2003), and may increase the risk of civil wars in tropical countries (Hsiang et al. 2011). There is some evidence that ENSO variability is influenced by global warming and that more frequent and stronger El Niño events may pertain in the future (Zhang et al. 2008, Yeh et al. 2009) with potentially far-reaching conse- quences. The existence of the El Niño cycle in ancient environments is therefore of considerable applied interest because these environments may serve as analogues for future conditions. In this respect it is widely appreciated that the climate during the Pliocene Warm Period is of particular relevance because it was a time of sustained global warmth when the land-sea configuration, continental positions, patterns of ocean circulation, and biota were similar to the present day (Dowsett and Robinson 2009). Evidence that the Pliocene Warm Period was characterized by a permanent El Niño (Wara et al. 2005, Fedorov et al. 2006, Ravelo et al. 2006) or a permanent La FastTrack➲ publication BULLETIN OF MARINE SCIENCE. VOL 89, NO 4. 20132 Niña (Rickaby and Halloran 2005) has therefore attracted considerable attention, but paleoclimatologists have thus far failed to concur. Proxies that have been used to characterise El Niño events in past environments include measurements of organismal growth, skeletal isotopic composition, and life histories associated with organisms such as trees, diatoms, bivalves, corals, and fora- minifera (e.g., Conroy et al. 2009, Khider et al. 2011, Li et al. 2011, Welsh et al. 2011). The underlying rationale is to obtain evidence for inter-annual variability in tem- peratures consistent with departures from long-term average temperatures. Thus, El Niño and La Niña events are characterized respectively by elevated or reduced temperature regimes relative to long-term averages in some specified region of the ocean. Despite this, remarkably little evidence of ENSO variation in ancient environments has been gathered. Annual growth increments and geochemical signals in corals (Watanabe et al. 2011), growth increments in bivalves and wood (Ivany et al. 2011), and combined geochemical data from foraminifera (e.g., Rickaby and Halloran 2005, Wara et al. 2005) are among the exceptions. Here we examine whether fossil chei- lostome bryozoans can provide a novel proxy for inferring ancient ENSO variability when used collectively to produce frequencies of estimates of mean annual ranges in temperature (MART). The method we develop offers the advantage of incorporat- ing data from numerous specimens—an approach similar to that recently developed for analyzing distributions of δ18O measures from individual foraminiferal skeletons (Khider et al. 2011). The rationale used here is that an absence of ENSO variability will be exhibited by a normal distribution of MART estimates obtained from bryozoan assemblages, reflecting convergence on a single value for MART in a given environment. However, the inter-annual variation in temperature regimes associated with ENSO variability should be characterised by values that include low MART estimates reflecting peri- ods when either warmer or cooler temperatures predominate throughout the year (associated with El Niño and La Niña events, respectively) and higher MART esti- mates reflecting greater intra-annual variation in temperature for transitional (i.e., years bounding ENSO events) and “normal” periods. Given such a regime, ENSO variability should be reflected by a relatively broad range of MART values arising from different years. Distributions of such MART values are likely to show little cen- tral tendency and instead may be described by relatively large negative kurtosis val- ues (resulting from flat distributions), by relatively large negative or positive values of skewness (resulting from left or right skewed distributions, respectively), or bimodal- ity (resulting from the record of two extreme climates reflecting years affected and unaffected by ENSO events). Here we examine whether the distributions of MART estimates obtained from time-averaged bryozoan assemblages from Panama can be used to provide evi- dence for ENSO variability in material collected from late-Miocene to Pleistocene environments of the tropical eastern Pacific (TEP) and southwestern Caribbean (SWC). We also analyze data available for present-day bryozoan assemblages from these Panamanian environments since ENSO variability is now absent in the SWC (Kaufmann and Thompson 2005), but may be expected in the Gulf of Panama (TEP) given the evidence for periodic upwelling in similar sites in the region (Gonzalez- Silvera et al. 2004, Palacios and Bograd 2005, Karnauskas et al. 2008). In addition to determining whether distributions of bryozoan-derived MART estimates may O K A M U R A E T A L . : B R Y O Z O A N P R O X Y R E V E A L S A N C I E N T E N S O VA R I A B I L I T Y 3 inform on ENSO variability, our data allow us to test the controversial hypothe- ses that the Pliocene Warm Period was characterized by (1) a permanent El Niño (Wara et al. 2005, Fedorov et al. 2006, Ravelo et al. 2006) or (2) a permanent La Niña (Rickaby and Halloran 2005). Methods Using Bryozoans to Estimate MART.—Bryozoans are colonial, suspension feeding invertebrates common in benthic assemblages. There are some 6000 described species at the present day (Gordon et al. 2009) and most belong to the order Cheilostomata. Colonies of cheilostomes are comprised of asexually budded zooids that are reinforced by skeletal walls composed of calcitic and/or aragonitic carbonate (McKinney and Jackson 1989). Zooid sizes are fixed at budding (Fig. 1) and their size varies as a function of seawater tem- perature according to the temperature-size rule (see below). This temperature-dependency of zooid size can therefore provide an indirect record of the seasonal temperature regime experienced by individual colonies. The zooid-size MART (zs-MART) approach uses a predictive model to infer MART (O’Dea and Okamura 2000a). This predictive model is based on the empirically derived relationship of the mean coefficient of variation (CV) of zooid size and MART. The CV of zooid size was calculated from 157 colonies of 29 cheilo- stome species collected from a variety of environments from tropical to polar regions, and the MART associated with each colony calculated by averaging the differences in annual variation in benthic temperatures obtained from direct thermal measurements from loca- tions close to where the colonies were collected. Algebraic rearrangement of the regression that describes the relationship provides a model for predicting MART as: MART = −3 + 0.745(b), where b is the mean intracolonial CV of zooid size. The regression that provides the basis for this model has 95% confidence limits within ±1 °C across the entire tempera- ture range (see fig. 2C in O’Dea and Okamura 2000a). Model predictions have been used successfully to infer MARTs by examining material of Recent (Panama: O’Dea and Jackson 2002; Wales: Knowles et al. 2010) and fossil (Pliocene UK and Miocene France: O’Dea and Okamura 2000a; seven localities in the Pliocene North Atlantic margin ranging from Panama to the UK: Knowles et al. 2009; ancient Panamanian environments during and Figure 1. Part of a Recent colony of Stylopoma sp. from Puerto Rico showing new zooids develop- ing at the growing edge (top) and length (zl) and width (zw) measurements made for zooid-size mean annual ranges in temperature analysis (inset). Once budded, zooids remain fixed in size. BULLETIN OF MARINE SCIENCE. VOL 89, NO 4. 20134 after the closure of the Isthmus seaway: O’Dea et al. 2007a, 2012; early Pliocene Weddell Sea in Antarctica: Clark et al. 2010; Denmark prior to the K-T boundary: O’Dea et al. 2011) provenance (see Okamura et al. 2011 for review). The range of environments so charac- terized demonstrates that the approach can be informative even in, for example, some tropical regions experiencing little annual variation in temperature (e.g., the Caribbean coast of Panama (O’Dea and Jackson 2002). However, zs-MART data have not been used collectively to demonstrate interannual variability in seawater temperature regimes by analyzing frequency distributions as done here. Data for MART inferences are based on measurements of zooid size as estimated by the frontal area of autozooids (normal feeding zooids), calculated as the product of maximum length and width (Fig. 1). It is of utmost importance to adhere to strict criteria in data col- lection (O’Dea and Okamura 2000a, Okamura et al. 2011). Thus, we measured the maximum lengths and widths of 20 randomly-selected autozooids per colony. These zooids were re- quired to be part of the basal series of zooids (not frontally budded), to show unimpeded growth and to be outside of the zone of astogenetic change (avoiding early colony growth stages when zooids are small). Colonies chosen for measurement were required to offer at least 30 ontogenetically complete autozooids and to show growth that was not influenced by substratum irregularities or competition with other sessile organisms. They also needed to have clearly delimited zooid margins for measurement and no evidence of distortion in dried material. The inverse relationship between body size and temperature, known as the temperature- size rule (Atkinson 1994), is widely regarded as a universal phenomenon describing ecto- therms. Although the mechanistic basis is subject to debate (see Okamura et al. 2011 for review), the generality of this inverse relationship enables a means of inferring temperatures based on body size. For colonial animals like bryozoans, the relationship is expressed at the zooid level, and there is a reasonable body of evidence that zooid size is indeed inversely related to temperature. This evidence comes from diverse laboratory-reared material (e.g., Menon 1972, Hunter and Hughes 1994, Atkinson et al. 2006, Amui-Vedel et al. 2007, O’Dea et al. 2007b), diverse field-collected material (e.g., Okamura 1987, O’Dea and Okamura 2000b, 2000c, O’Dea and Jackson 2002, O’Dea 2005, Lombardi et al. 2006), and the growth of bryo- zoans in the field (O’Dea and Okamura 1999; see Okamura et al. 2011 for further review) There are other factors that influence the ultimate size of zooids, the most obvious being species-specific colonial growth rules, such as the production of annual growth check lines (O’Dea and Okamura 2000b, O’Dea 2005, Okamura et al. unpubl data), substrate irregulari- ties, predation, and competition from other sessile encrusters, the influence of which can be minimized with the judicious rejection of randomly selected zooids as the approach requires (O’Dea and Okamura 2000a, O’Dea 2005). Salinity and food availability (Hageman et al. 2009) have also been shown to influence zooid size, yet the effects are minimal in comparison with those of temperature (Okamura 1987, Hunter and Hughes 1994, O’Dea and Okamura 1999, 2000b, O’Dea 2005, O’Dea et al. 2007b, Okamura et al. 2011). Calculating the CV in zooid size standardizes this relationship among taxa, enabling the collection of common zs- MART estimates from a range of species and genera. Moreover, the method developed here is based on frequency distributions, which should help to reduce the effect of noise in obscuring the signal of interest provided samples sizes are large. Study Material.—We analyzed zs-MART data obtained from cheilostome bryozoans collected for use in previous studies as well as de novo for the present study (Table 1). Recent bryozoans from the Gulf of Panama (TEP), the Gulf of Chiriquí (TEP), San Blas (SWC), and Bocas del Toro (SWC) were obtained from dredge and trawl collections between the years 1996 and 2000 (O’Dea et al. 2004, O’Dea et al. 2007a; Fig. 2). Note that such dredge col- lections represent time-averaged death assemblages, akin to fossil assemblages from bulk samples (Best and Kidwell 2000). Fossil bryozoans from Caribbean sites were extracted from O K A M U R A E T A L . : B R Y O Z O A N P R O X Y R E V E A L S A N C I E N T E N S O VA R I A B I L I T Y 5 Ta bl e 1. L oc al it ie s, b as in (f or R ec en t s it es )/ ge ol og ic al fo rm at io ns (f or fo ss il lo ca li ti es ), la ti tu de (L at ), lo ng it ud e (L on g) , m ed ia n ag e (a nd ra ng e) , t im e pe ri od (E po ch ), s am pl e si ze s [N = th e to ta l n um be r o f c ol on ie s sa m pl ed fo r z oo id -s iz e m ea n an nu al ra ng es in te m pe ra tu re (z s- M A R T ) a na ly si s] a nd th e na m es o f t he s pe ci es fr om w hi ch z s- M A R T d at a w er e co ll ec te d (n um be r of c ol on ie s sa m pl ed p er s pe ci es ). T E P = tr op ic al e as te rn P ac ifi c, S W C = s ou th w es te rn C ar ib be an . L oc al it y R eg io n B as in /f or m at io n L at (° N ) L on g (° W ) F or m at io n (d ep th ) M ed ia n ag e (r an ge ) E po ch N S pe ci es s am pl ed [ no . c ol on ie s] G ul f of P an am a T E P G ul f of P an am a 8. 27 79 .4 0 R ec en t (1 2– 55 m ) R ec en t 43 C up ul ad ri a ex fr ag m in is H er re ra -C ub il la , D ic k, S an ne r an d Ja ck so n, 2 00 6 [1 8] ; c f. C on op eu m r et ic ul um ( L in na eu s, 17 67 ) [5 ]; H ip po po di na fe eg ee ns is ( B us k, 1 88 4) [ 5] ; P ar as m it ti na c ro ss la nd ii ( H as ti ng s, 1 93 0) [ 5] ; O ny ch oc el la cf . a ng ul os a (R eu ss , 1 84 7) [ 5] ; S in ofl us tr a an na e (O sb ur n, 19 53 ) [5 ] G ul f of C hi ri qu í T E P G ul f of C hi ri qu í 7. 77 82 .0 1 R ec en t (3 2 m ) R ec en t 13 C up ul ad ri a ex fr ag m in is [ 13 ] S an B la s S W C S an B la s 9. 56 78 .8 0 R ec en t (2 –3 9 m ) R ec en t 24 St yl op om a sp . 2 [ 5] , S te gi ni po re ll a m ag ni la br is ( B us k, 1 85 4) [5 ]; S ty lo po m a sp on gi te s (P al la s, 1 76 6) [ 5] ; C up ul ad ri a sp . [8 ] B oc as d el T or o S W C B oc as d el T or o 9. 36 82 .3 3 R ec en t (3 1 m ) R ec en t 10 C up ul ad ri a sp . [ 10 ] B ur ic a P en in su la T E P G ul f of C hi ri qu í 8. 28 82 .8 7 A rm ue ll es 0. 36 ( 0. 10 ) P le is to ce ne 36 T ha la m op or el la s p. [ 4] ; P ue ll in a sp . [ 2] ; O ny ch oc el la s p. [6 ]; F lo ri di na s p. [ 4] ; ? E sc ha ri na [ 1] ; T ry po st eg a sp . [ 4] ; ?C op id oz ou m s p. [ 1] ; ? H ip po pl eu ri fe ra [ 1] ; M et ra ra bd ot os pa ci fic um O sb ur n, 1 95 2 [3 ]; A ca nt ho de si a sp . [ 1] ; M ic ro po re ll a sp . [ 3] ; S te gi no po re ll a co rn ut a O sb ur n, 1 95 0 [1 ]; P ar as m it ti na s p. [ 1] ; A nt ro po ra le uc oc yp ha M ar cu s, 19 37 [ 3] ; ? A im ul os ia s p. [ 1] S w an C ay S W C B oc as d el T or o 9. 45 82 .3 0 S w an C ay 1. 40 ( 0. 20 ) P le is to ce ne 12 C up ul ad ri a sp . [ 12 ] U pp er L om as d el M ar S W C L im ón 10 .0 0 83 .0 4 M oi n 1. 60 ( 0. 10 ) P le is to ce ne 15 C up ul ad ri a ch ee th am i H er re ra -C ub il la , D ic k, S an ne r an d Ja ck so n, 2 00 6 [1 1] ; D is co po re ll a sc ut el la H er re ra -C ub il la , D ic k, S an ne r an d Ja ck so n, 2 00 8 [4 ] L ow er L om as d el M ar S W C L im ón 10 .0 0 83 .0 4 M oi n 1. 70 ( 0. 20 ) P le is to ce ne 10 C up ul ad ri a ch ee th am i [ 10 ] N or th w es t E sc ud o de V er ag ua s S W C B oc as d el T or o 9. 11 81 .5 7 E sc ud o de V er ag ua s 2. 00 ( 0. 10 ) P le is to ce ne 10 D is co po re ll a sc ut el la [ 10 ] W il d C an e C ay S W C B oc as d el T or o 8. 81 81 .1 0 G ro un d C re ek 2. 05 ( 0. 15 ) P le is to ce ne 6 D is co po re ll a sc ut el la [ 1] ; C up ul ad ri a ch ee th am i [ 5] F is h H ol e S W C B oc as d el T or o 10 .0 1 83 .0 6 G ro un d C re ek 2. 60 ( 0. 40 ) P li oc en e 5 C up ul ad ri a bi po ro sa ( C an u an d B as sl er , 1 92 3) [ 5] N or th C en tr al E sc ud o de V er ag ua s S W C B oc as d el T or o 9. 11 81 .5 7 E sc ud o de V er ag ua s 2. 75 ( 0. 85 ) P li oc en e 5 C up ul ad ri a pa na m en si s H er re ra -C ub il la , D ic k, S an ne r an d Ja ck so n, 2 00 6 [5 ] L a B om ba S W C L im ón 9. 91 83 .0 7 R ío B an an o 3. 05 ( 0. 15 ) P li oc en e 11 D is co po re ll a sc ut el la [ 9] ; C up ul ad ri a bi po ro sa [ 2] BULLETIN OF MARINE SCIENCE. VOL 89, NO 4. 20136 Ta bl e 1. C on ti nu ed . L oc al it y R eg io n B as in / fo rm at io n L at (° N ) L on g (° W ) F or m at io n (d ep th ) M ed ia n ag e (r an ge ) E po ch N S pe ci es s am pl ed [ no . c ol on ie s] B ru no B lu ff S W C B oc as d el T or o 9. 04 81 .7 4 S ha rk H ol e 3. 45 ( 0. 15 ) P li oc en e 10 D is co po re ll a sc ut el la [ 10 ] C ay o A gu a – P t T ib ur ón to P t P ie dr a S W C B oc as d el T or o 9. 15 82 .0 2 C ay o A gu a 3. 55 ( 0. 05 ) P li oc en e 6 C up ul ad ri a bi po ro sa [ 1] ; D is co po re ll a sc ut el la [ 5] Is la S ol ar te S W C B oc as d el T or o 10 .0 1 83 .0 6 O ld B an k 3. 55 ( 0. 05 ) P li oc en e 10 D is co po re ll a bo ca sd el to ro en si s H er re ra -C ub il la , D ic k, S an ne r an d Ja ck so n, 2 00 8 [5 ]; D is co po re ll a sc ut el la [ 5] S an ta R it a S W C B oc as d el T or o 9. 97 83 .1 2 R ió B on an o 3. 55 ( 0. 05 ) P li oc en e 10 C up ul ad ri a bi po ro sa [ 9] ; D is co po re ll a sc ut el la [ 1] C ay o A gu a P un ta N ís pe ro S W C B oc as d el T or o 9. 17 82 .0 3 C ay o A gu a 3. 55 ( 0. 05 ) P li oc en e 10 C up ul ad ri a pa na m en si s [3 ]; C up ul ad ri a bi po ro sa [ 1] ; D is co po re ll a sc ut el la [ 1] ; D is co po re ll a bo ca sd el to ro en si s [5 ] C ay o A gu a P un ta P ie dr a R oj a S W C B oc as d el T or o 9. 14 82 .0 2 C ay o A gu a 4. 25 ( 0. 75 ) P li oc en e 10 D is co po re ll a sc ut el la [ 2] ; D is co po re ll a bo ca sd el to ro en si s [8 ] Is la P op a S W C B oc as d el T or o 10 .0 1 83 .0 6 C ay o A gu a 4. 25 ( 0. 75 ) P li oc en e 10 D is co po re ll a sc ut el la [ 10 ] C ay o A gu a – w es t s id e P un ta S W C B oc as d el T or o 9. 18 82 .0 5 C ay o A gu a 4. 25 ( 0. 75 ) P li oc en e 10 C up ul ad ri a bi po ro sa [ 4] ; D is co po re ll a sc ut el la [ 2] ; D is co po re ll a ch ee th am i [ 1] ; D is co po re ll a bo ca sd el to ro en si s [3 ] S ha rk H ol e P oi nt S W C B oc as d el T or o 9. 04 81 .7 4 S ha rk H ol e P oi nt 5. 65 ( 0. 05 ) M io ce ne 11 C up ul ad ri a bi po ro sa [ 10 ]; D is co po re ll a sc ut el la [ 1] S ou th V al ie nt e w es t S W C B oc as d el T or o 11 .6 5 69 .7 8 S ha rk H ol e P oi nt 6. 29 ( 0. 97 ) M io ce ne 9 D is co po re ll a sc ut el la [ 8] ; C up ul ad ri a ch ee th am i [ 1] R ío T up is a S W C D ar ié n 8. 30 77 .6 1 T iu rà 6. 35 ( 0. 75 ) M io ce ne 8 D is co po re ll a sc ut el la [ 8] R ío C hi co S W C D ar ié n 8. 17 77 .7 0 T iu rà 6. 35 ( 0. 75 ) M io ce ne 8 D is co po re ll a sc ut el la [ 6] ; C up ul ad ri a bi po ro sa [ 2] C ay o P at te rs on S W C B oc as d el T or o 9. 10 81 .8 9 S ha rk H ol e P oi nt 6. 90 ( 1. 30 ) M io ce ne 7 C up ul ad ri a bi po ro sa [ 7] R ío C al zo ne s S W C C an al 9. 04 80 .6 5 U nn am ed 8. 25 ( 2. 95 ) M io ce ne 11 C up ul ad ri a bi po ro sa [ 3] ; D is co po re ll a sc ut el la [ 8] M at tr es s F ac to ry S W C C an al 9. 37 79 .8 3 G at un 9. 00 ( 0. 40 ) M io ce ne 11 D is co po re ll a sc ut el la [ 5] ; C up ul ad ri a bi po ro sa [ 6] S M R Q ua rr ie s, S ar as ot a G ul f of M ex ic o Ta m ia m i 27 .3 7 82 .3 9 Ta m ia m i 3. 15 ( 0. 70 ) P li oc en e 28 Tr yp os te ga v en us ta ( N or m an , 1 86 4) [ 9] ; M ic ro po re ll a af f. c il ia ta ( P al la s, 1 76 6) [ 10 ]; C up ul ad ri a bi po ro sa [ 4] ; D is co po re ll a sp . [ 5] C hu ck at uc k W es te rn A tl an ti c Y or kt ow n 36 .8 9 76 .5 9 Y or kt ow n 3. 40 ( 1. 20 ) P li oc en e 70 P ue ll in a ca pr on en si s W in st on , 2 00 5 [1 0] ; c f. A ca nt ho de si a co m m en sa le ( K ir kp at ri ck a nd M et ze la ar , 1 92 2) [ 8] ; F lo ri di na re gu la ri s C an u an d B as sl er , 1 92 3 [1 0] ; H ip pa li os in a ro st ri ge ra ( S m it t, 18 73 ) [1 0] ; R ep ta de on el la c ol li ns ae C he et ha m , S an ne r an d Ja ck so n, 2 00 7 [1 0] ; H ip po po ri na s p. [1 0] ; M ic ro po re ll a sp . [ 6] ; T ry po st eg a ve nu st a [6 ] O K A M U R A E T A L . : B R Y O Z O A N P R O X Y R E V E A L S A N C I E N T E N S O VA R I A B I L I T Y 7 geological bulk samples from outcrops on land to provide a large data set for the Caribbean over time (O’Dea et al. 2007a, O’Dea and Jackson 2009; Fig. 2). Fossil bryozoans from the Burica Peninsula (TEP) were collected individually from the Armuelles Formation in March 2011. We also incorporate fossil material hand-picked from geological outcrops in Florida and Virginia, USA (Knowles et al. 2009; Fig. 2), to present a broader picture of ENSO variability in the Pliocene, to complement the picture gained for the TEP. Some data sets are based on a number of species, while others are based on a single spe- cies (Table 1). Because the zs-MART approach is taxon independent (see above), variation in representation of species and genera should not influence our interpretations. Ages of all fossil samples are median values of minimum and maximum ages based on mi- crofossils from the sample or interpolated from ages of samples stratigraphically above and below (Coates et al. 1992, Aubrey and Berggren 1999, Cotton 1999, Jackson et al. 1999, Coates et al. 2004, 2005, Leon-Rodriguez 2007, O’Dea et al. 2007a, O’Dea et al. 2012). Data were plotted as frequency distributions of inferred MART values based on zs-MART estimates (rounded to the nearest 1 °C). We plotted these frequencies against a common range for the x-axis (based on the smallest and largest range across the data sets) to enable straightforward comparisons amongst localities. Anderson-Darling tests for normality were conducted to examine whether there is evidence for a central tendency in the distributions. Associated descriptive statistics for frequency distributions were also assembled for compar- ative insights on typical zs-MART values (mean, median), measure of dispersion (standard deviation) and shapes of distributions (skewness, kurtosis). Tests of skewness and kurtosis were employed to infer shapes of the distributions. We did not adopt a common range for the y-axis (e.g., by plotting as relative frequencies) in our figures of frequency distributions to il- lustrate directly the data analysed in statistical tests. zs-MART distributions for ancient environments from the large Caribbean dataset were analyzed by pooling data for zs-MART estimates to generate frequency distributions in sev- eral ways. Since one of our primary interests is in using zs-MART estimates to understand climate regimes during the Pliocene Warm Period, we assigned data according to established dates for this period. Thus, we grouped data (as median values) according to a Pliocene Warm Period of approximately 5–3 Ma, consistent with, for example, 4.5–3.0 Ma (Wara et al. 2005) and 5–4 Ma (Watanabe et al. 2011). Data for other time periods are assigned according to pre- and post-Pliocene Warm Periods. We also assigned Pliocene data according to the latest revi- sion of the geological time scale by the International Committee of Stratigraphy (Gibbard et Figure 2. (A) Context map with localities in Virginia (Yorktown) and Florida (Tamiami) indicated, and (B) detailed map of localities in the Isthmus of Panama. Sampling locali- ties of Recent (filled circles) and fossil (open circles) collections. BULLETIN OF MARINE SCIENCE. VOL 89, NO 4. 20138 al. 2010) into the Zanclean (from 5.33 to 3.5 Ma) and Piacenzian (from 3.6 to 2.58 Ma) stages. During the latter, global mean annual surface temperatures are estimated to have been 2–3 °C higher than the present day (Haywood et al. 2000, Lunt et al. 2010). Because the median date of 3.55 Ma assigned to a number of fossil bryozoan colonies was near the Zanclean/ Piacenzian border of 3.5–3.6 Ma we included analyses assigning data for this date to the Piacenzian and the Zanclean for comparative purposes. Data for other time periods were grouped into the Pleistocene and the Messinian and Tortonian stages of the Miocene. Results Data for Recent Bryozoans.—Figure 3 presents frequency distributions of Recent inferred MART estimates in the SWC (San Blas and Bocas del Toro) and the TEP (Gulf of Panama and Gulf of Chiriquí). The distribution for San Blas shows a relatively small range in inferred MART values (range = 4 °C; Table 2), an associated low standard deviation (SD = 0.95; Table 2) and relatively small mean and median values (mean = 3.08 °C; median = 2.91 °C; Table 2) reflective of the seasonally stable water temperatures in the Caribbean (O’Dea and Jackson 2002a; Fig. 3A). It is consistent with a normal distribution (Anderson-Darling test: AD = 0.484, P = 0.207). The Bocas del Toro distribution is based on a relatively small sample size (n = 10; Table 1) and shows a very narrow range, a low degree of varia- tion, and relatively low MART values (range = 1 °C, SD = 0.61 °C, mean = 3.39 °C, median = 3.34 °C; Table 2), and is also consistent with a normal distribution (AD = 0.328, P = 0.449; Fig. 3B). Inferred MART estimates from both the Pacific sites are distributed across a relatively broad range of values that are typically higher than those character- izing the Caribbean. Thus, the Gulf of Panama distribution shows a range of 8 °C (Table 2), an associated high standard deviation (SD = 1.96 °C; Table 2), rela- tively high mean and median values (mean = 8.01 °C, median = 7.59 °C; Table 2) and is significantly different from normal (AD = 0.951, P = 0.015). The skewness Figure 3. Frequency distributions of inferred mean annual ranges in temperature (MART) values (based on zooid-size MART estimates) from bryozoans collected in the present day from the Caribbean coast of Panama (A: San Blas, B: Bocas del Toro) and the Pacific coast of Panama (C: Gulf of Panama; D: Gulf of Chiriquí). Sample sizes (n) and associated Anderson Darling test statistic (AD) and result are provided with each distribution. O K A M U R A E T A L . : B R Y O Z O A N P R O X Y R E V E A L S A N C I E N T E N S O VA R I A B I L I T Y 9 test demonstrated a significant positive skew in the distribution resulting from a number of high MART estimates in the right-hand tail (Fig. 3C; legend for Table 2). Although based on a small sample size (n = 13), the Gulf of Chiriquí distribu- tion also exhibits a relatively large range, high variability and high inferred MART values (range = 7 °C, SD = 2.60 °C, mean = 5.83 °C, median = 5.55 °C; Table 2). The Anderson-Darling test suggests (with borderline significance) that the Gulf of Chiriquí MART estimates are not normally distributed (AD = 0.687, P = 0.055; Fig. 3D). Caribbean Fossil Data.—Figure 4 presents frequency distributions of MART estimates for the Caribbean fossil material collected during time intervals ac- cording to post-Pliocene, Pliocene, and pre-Pliocene Warm Periods. The material representing the post-Pliocene Warm Period is characterised by a narrow inferred MART range (5 °C), an associated low standard deviation (1.14) and low mean and median values (3.3 and 3.2, respectively) relative to the values for the Pliocene Warm Period and pre-Pliocene Warm Period (Table 2). The distribution is consis- tent with normality (AD = 0.657, P = 0.082; Fig. 4A). Material analyzed from both the Pliocene Warm Period and the pre-Pliocene Warm Period provides broader inferred MART ranges (7 °C in both cases), higher associated standard deviations Table 2. Descriptive statistics for distributions [for inferred mean annual range in temperature (MART) as based on zooid-size MART data] in Figures 3–6 (and age of distributions). Skewness (S) and kurtosis (K) values are accompanied by test statistics, Z S and Z K , which indicate when critical values are exceeded. Z S = S/ SES, where S is the sample skewness value and SES is the standard error of skewness: √6n (n−1)/(n−2)(n+1) (n+3). Z K = K/SEK, where K is the sample kurtosis value and SEK is the standard error of kurtosis: 2(SES) √(n2 − 1)/(n − 3)(n + 5). In both cases critical values are when the test statistics are >2 or < −2. Zs values >2 indicate that the population is very likely skewed positively, those < −2 that the population is very likely to be skewed negatively. Z K values >2 indicate a population with a sharper peak and longer tails than a normal distribution, values < −2 describe a population with a lower and broader peak and shorter tails than a normal distribution (http://www.tc3.edu/instruct/sbrown/stat/shape.htm; and adapted from Cramer 1997). Z S and Z K values that exceed these critical values are indicated by asterisk. n = sample size; SD = standard deviation. The statistics for the Pliocene Piacenzian and Zanclean are according to the alternate assignments of material aged near the boundary (see Methods and ages below). Distribution Age (Ma) Mean Median n SD Range Skewness (Z S ) Kurtosis (Z K ) Gulf of Panama Recent 8.0 7.6 43 1.96 8 0.87 (2.41)* 0.27 (0.38) Gulf of Chiriqui Recent 5.8 5.6 13 2.60 7 0.24 (0.39) −1.66 (−1.39) San Blas Recent 3.1 2.9 23 0.94 4 −0.16 (−0.33) −0.22 (−0.24) Bocas del Toro Recent 3.4 3.3 10 0.61 1 0.24 (0.35) −0.84 (−0.63) Post-Pliocene Warm Period 1.40–2.75 3.3 3.2 63 1.14 5 0.25 (0.83) −0.76 (−1.28) Pliocene Warm Period 3.05–4.25 5.6 6.1 87 2.11 7 −0.06 (−0.23) −1.19 (−2.33)* Pre-Pliocene Warm Period 5.65–9.00 6.5 7.0 43 2.27 7 −0.15 (−0.42) −1.33 (−1.90) Pleistocene 1.40–2.05 3.5 3.5 53 1.08 4 0.25 (0.76) −0.67 (−1.04) Pliocene Piacenzian (1) 2.60–3.55 5.2 4.7 67 2.33 8 0.02 (0.07) −1.25 (−2.17)* Pliocene Piacenzian (2) 2.60–3.45 4.2 3.4 31 2.35 8 0.63 (1.50) −0.99 (−1.21) Pliocene Zanclean (1) 4.25 5.5 5.5 30 2.06 7 0.11 (0.26) −1.22 (−1.47) Pliocene Zanclean (2) 3.55–4.25 5.9 6.3 66 1.99 7 −0.04 (−0.14) −1.16 (−1.99) Miocene Messinean 5.65–6.90 6.5 7.0 43 2.27 7 −0.15 (−0.42) −1.33 (−1.90) Miocene Tortonian 8.25–9.00 6.1 5.6 22 1.20 6 1.55 (3.16)* 3.11 (3.26)* Armuelles Formation 0.26–0.46 5.9 5.8 36 1.93 7 0.29 (0.74) −0.81 (−1.05) Tamiami Formation 2.80–4.00 6.4 6.2 28 1.90 6 0.16 (0.36) −1.06 (−1.20) Yorktown Formation 2.80–3.50 6.5 6.6 70 1.24 7 −0.25 (−0.87) 0.62 (−1.09) BULLETIN OF MARINE SCIENCE. VOL 89, NO 4. 201310 (2.11 and 2.27 °C, respectively), higher mean (5.6 and 6.5 °C, respectively), and median (6.1 and 7.0 °C, respectively) inferred MART values (see Table 2) and non- normal distributions (Pliocene Warm Period: AD = 1.325, P < 0.005; pre-Pliocene Warm Period: AD = 0.902, P = 0.020; Fig. 4B,C). The kurtosis test (Table 2) dem- onstrates significant negative kurtosis (a flatter than normal distribution) for the distribution of the Pliocene Warm Period. However, a skewness value close to zero (−0.06; Table 2) provides evidence that values tend to group as mirror images on either side of the centre of the distribution since the data are not normally distrib- uted. This tendency is clear from inspection of the distribution. The pre-Pliocene Warm Period distribution is similar in shape that of the Pliocene Warm Period but the kurtosis test statistic (−1.9) did not quite exceed the critical value (−2.0). Similar patterns were obtained when the data were assigned to a greater num- ber of time bins (Fig. 5). Thus, Pleistocene material was characterised by a relative- ly narrow range of inferred MART values, an associated low standard deviation, relatively low mean and median values (Table 2) and a distribution consistent with normality (AD = 0.480, P = 0.224; Fig. 3A). Material analysed from all other time intervals demonstrated relatively broader inferred MART ranges and vari- ability, higher mean and median values (Table 2) and distributions that gener- ally showed no central tendency (AD test results associated with time intervals, Fig. 5). The exception to non-normality depended on whether data from the bor- derline date of 3.55 Ma was assigned to the Piacenzian or the Zanclean. In the former case, the Anderson-Darling test for normality provided a non-significant result for the Zanclean (Fig. 5D,E). Negative kurtosis characterises the distribu- tion for the Piacenzian stage of the Pliocene (based on 3.55–2.6 Ma; Fig. 5B) while the distribution for the Zanclean stage of the Pliocene (based on 4.25 Ma; Fig. 5D) was nearly significant with a test statistic value of −1.99 (critical value > −2). The Miocene Tortonian distribution demonstrated significant positive skewness Figure 4. Frequency distributions of inferred mean annual ranges in temperature (MART) values (based on zooid-size MART estimates) from fossil bryozoans from the Caribbean coast of Panama assigned to: (A) the post-Pliocene Warm Period, (B) the Pliocene Warm Period, (C) the pre-Pliocene Warm Period. Range in age of material (in millions of years) noted in parentheses. Sample sizes (n) and associated Anderson Darling test statistic (AD) and result are provided with each distribution. O K A M U R A E T A L . : B R Y O Z O A N P R O X Y R E V E A L S A N C I E N T E N S O VA R I A B I L I T Y 11 (Table 2, Fig. 5G). In addition, skewness values close to zero and non-normality for the Piacenzian stage (based on 3.55–2.6 Ma; Fig. 5B) and the Zanclean stage (based on 4.25–3.55 Ma; Fig. 5E) provide evidence that values group as mirror im- ages on either side of the centres of the distributions, a tendency that can be seen by inspection of the distributions. Other Fossil Data.—Figure 6 provides distributions of inferred MART val- ues that characterize oceanographic conditions for the Pleistocene (Armuelles Formation) of the Burica Peninsula (TEP) and Pliocene deposits of Florida (the Tamiami Formation) and Virginia (Yorktown Formation). The three distribu- tions show similar ranges of MART estimates (7, 6, and 7 °C for the Armuelles, Figure 5. Frequency distributions of inferred mean annual ranges in temperature (MART) values (based on zooid-size MART estimates) from fossil bryozoans from the Caribbean coast of Panama assigned to the (A) Pleistocene, the (B–E) Piacenzian and Zanclean stages of the Pliocene (according to the latest revision of the geological time scale; Gibbard et al. 2010), the (F) Messinian and (G) Tortonian stages of the Miocene. Because a number of fossil specimens are assigned a median date of 3.55 Ma near the Zanclean/Piacenzian boundary, we provide two alternative breakdowns for the data for these time intervals. Range in age of material (in millions of years) noted in parentheses. Sample sizes (n) and associated Anderson Darling test statistic (AD) and result are provided with each distribution. BULLETIN OF MARINE SCIENCE. VOL 89, NO 4. 201312 Tamiami, and Yorktown formations, respectively; Table 2). The distributions for the Armuelles and Tamiami formations are characterised by higher variability (SD = 1.93 and 1.90 °C, respectively) than that of the Yorktown Formation dis- tribution (SD = 1.24 °C). All three distributions share similar mean and median values (mean = 5.9, 6.4, and 6.5 °C; median = 5.8, 6.2, and 6.6 °C for the Armuelles, Tamiami, and Yorktown formations, respectively), and all are consistent with nor- mality (Armuelles Formation: AD = 0.415, P = 0.318; Tamiami Formation: AD = 0.398, P = 0.344; Yorktown Formation: AD = 0223, P = 0.820). Discussion Environmental Variability in Present-day Panama.—Our inferred MART estimates for the modern Gulf of Panama depict an environment which experiences the greatest range of annual temperature variation relative to all other times and places characterised in this study. In some years MART estimates indicate an annual temperature variation in the Gulf of Panama of as little as 5 °C, while other years may experience a temperature range of up to 13 °C. The frequency distribution of MART estimates did not conform to normality (Fig. 1A). These patterns and their magnitude are therefore consistent with what would be expected for an environment subjected to both inter-annual (=ENSO variability) and intra-annual variation in exposure to cold, upwelling, and warm surface waters. Exactly how ENSO affects the strength of upwelling in the Gulf of Panama remains to be determined, but it is likely that the response to ENSO variability is similar to those observed in the Gulfs of Tehuantepec and Papagayo which also experience wind-jet driven seasonal upwelling in the TEP. In these sites, wind and temperature patterns during El Niño and La Niña oscillation suppress and enhance upwelling re- spectively (Gonzalez-Silvera et al. 2004, Palacios and Bograd 2005, Karnauskas et al. Figure 6. Frequency distributions of inferred mean annual ranges of temperature (MART) values (based on zooid-size MART estimates) for fossil bryozoans from: (A) the Armuelles Formation of the Pacific coast of Panama, (B) the Tamiami Formation of Florida, and (C) the Yorktown Formation of Virginia. Estimated ages of the formations provided in parentheses. Sample sizes (n) and associated Anderson Darling test statistic (AD) and result are provided with each distribution. O K A M U R A E T A L . : B R Y O Z O A N P R O X Y R E V E A L S A N C I E N T E N S O VA R I A B I L I T Y 13 2008). If we assume that the Gulf of Panama responds similarly to ENSO variability given that the drivers of upwelling along the three gulfs are the same (D’Croz and O’Dea 2007), the positive skewness associated with inferred MART frequencies in the Gulf of Panama (Fig. 3A) can be interpreted to reflect a system that is character- ised by ENSO dynamics in the majority of years and relatively few years in which very deep cold upwelling and warm surface waters are both experienced. Inferred MART values for the Gulf of Chiriquí are lower than those for the Gulf of Panama as expected given the lack of strong seasonal upwelling. Nevertheless, the wide inter-annual variation in temperature regimes (Fig. 3D), as much as in the Gulf of Panama, can likewise be attributed to ENSO variability. This wide range of zs-MARTs in the Gulf of Chiriquí (3–10 °C) reflects the seasonal vertical migration of the thermocline (D’Croz and O’Dea 2007), which likely occurs closer to surface waters in La Niña years and farther away in El Niño years. This inference requires further exploration given the current dearth of information on seasonal thermal pat- terns in the region. Distributions of inferred MARTs from both San Blas and Bocas del Toro (Fig. 3C,D) reflect the overall low seasonal variation in thermal regimes typical of the Caribbean coast of the Isthmus of Panama (D’Croz and Robertson 1997, D’Croz et al. 2005). Additionally, both distributions reveal little inter-annual variation and con- form to normality as would be expected given the absence of significant ENSO vari- ability in the seas of the SWC (Kaufmann and Thompson 2005) In view of these results we conclude that the wide and non-normal distributions of MART estimates characterize Isthmian coastal seas that are affected by ENSO, while normal distributions with relatively low variance apply in areas not influenced by ENSO events. ENSO Variability in Ancient Tropical American Environments.— Inferred MARTs from the Mio-Pliocene Caribbean coast of present-day Panama shifted from high values similar to those observed in the present day, seasonally upwelling Gulf of Panama, to modern Caribbean-like MARTs during the Late Pliocene as the Isthmus of Panama isolated the TEP and SWC (O’Dea et al. 2007a, Jackson and O’Dea 2013). All inferred MART distributions for the Mio-Pliocene failed to conform to normality, apart from the one exception from the Zanclean stage, and thus provide evidence that ENSO variability was a permanent feature of Tropical American environments through the Mio-Pliocene. In contrast, MART estimates for the Caribbean Pleistocene material were rela- tively low, narrowly distributed and collectively conformed to a normal distribution, suggesting that by this time the Caribbean coast of Panama was warm and stable (O’Dea et al. 2007a) and did not experience ENSO variability. This pattern was not observed for Pleistocene material from the Pacific coast of Panama (the 0.46–0.26 Ma Armuelles Formation) whose inferred MARTs are distributed across a relatively broad range of values similar to those observed in the present-day Gulf of Chiriquí. The differences in shapes of the two distributions (normality for the Armuelles Formation and non-normality for the Gulf of Chiriquí) suggest some differences in temperature regimes, but the sample size for the Gulf of Chiriquí is small. Ancient Western Atlantic Environments.—To examine further the util- ity and inferences that may be gained by focusing on frequency distributions of in- ferred MART values we applied the approach to Pliocene material from the Tamiami BULLETIN OF MARINE SCIENCE. VOL 89, NO 4. 201314 Formation of Florida and the Yorktown Formation of Virginia (Fig. 6B,E). The in- ferred MART distribution for the Tamiami Formation shows a similar magnitude and range to those characterising the Pliocene material from the Caribbean of Panama. This pattern might be expected for the Tamiami Formation if there were inter-an- nual variation in the degree of upwelling and hence ENSO variability. Upwelling has been inferred for the Tamiami Formation on the basis of abundant and diverse molluscs (e.g., Allmon 1988, 1993) and the presence of cormorant skeletons (Emslie 1992). However, recent high resolution isotopic and Sr/Ca analyses of gastropods suggests that the nutrient-rich conditions were derived from freshwater inputs (Tao and Grossman 2010). Salinity may have influenced zooid sizes although available evidence suggests the effect of salinity is roughly half that of temperature (O’Dea and Okamura 1999). Salinity fluctuations alone are therefore unlikely to explain the relatively broad range of MART estimates obtained for the Tamiami Formation. The conformation to normality of these MART values suggests that during deposition, the Tamiami Formation was characterised by exposure to more regular thermal con- ditions than those experienced by Panamanian Caribbean localities in the Pliocene. The unimodal distribution of inferred MART values for the Yorktown Formation of Virginia is in keeping with an environment subjected to regular colder winter tem- peratures and higher summer temperatures. Was There a Permanent El Niño or La Niña in the Pliocene?—Our analy- ses reject the hypotheses for either a permanent El Niño or a permanent La Niña in the Pliocene TEP. The pervasive ENSO variability in the Caribbean that persisted from 9.0–3.05 Ma demonstrated here can be best explained by oceanographic con- tinuity with the Pacific combined with an ENSO cycle similar to that of today. Our study thus complements recent evidence that rejects a permanent El Niño in the Pliocene based on analyses of monthly-resolved fossil coral skeletons using δ18O and salinity proxies >35 yrs of growth in two specimens (Watanabe et al. 2011) and cou- pled climate model simulations that demonstrate ENSO variability in the Pliocene (e.g., Haywood et al. 2007, Scroxton et al. 2011, von der Heydt et al. 2011). These studies contradict those inferring a permanent El Niño in the Pliocene using geo- chemical and faunal-based estimates of sea surface temperatures derived from fora- minifera, alkenones (e.g., Wara et al. 2005, Fedorov et al. 2006), and the presence of warm-water molluscs in northern Chile (Ragaini et al. 2008). Our data, in the same way, reject the lesser supported hypothesis for a permanent La Niña in the Pliocene (Rickaby and Halloran 2005). At face value it is curious that our data appear to conflict with inferences of a per- manent El Niño in the Pliocene of northern Chile based on molluscan assemblage composition in the Peruvian Province (Ragaini et al. 2008). These Pliocene mollus- can assemblages are suggested to provide evidence for permanent warm conditions based on their similarity with molluscan assemblages observed today in the Panamic Province. However, since the TEP of Panama is characterised by El Niño variability in the present day, the presence of molluscs typical of this region in northern Chile in the Pliocene in fact is in keeping with our findings. Thus, together our study and that of Ragaini et al. (2008) suggest that the Peruvian Province was warmer in the Pliocene than at present as reflected by the extension of present day Panamanian molluscs into the Peruvian Province and that ENSO variability broadly characterised both Central and South American regions. Unfortunately, at present we lack bryo- zoan material from the Pliocene of South America to examine this further. O K A M U R A E T A L . : B R Y O Z O A N P R O X Y R E V E A L S A N C I E N T E N S O VA R I A B I L I T Y 15 The Caribbean Environment over Time.—Schematic summaries of frequen- cy distribution data reveal how coastal southwestern Caribbean environments have changed over the past approximately 9 million years (Fig. 7). Pleistocene and pres- ent-day environments both have low MARTs with normal distributions in keeping with the seasonally stable warm temperatures of the Caribbean and lack of ENSO variability after the rise of the Isthmus of Panama (O’Dea et al. 2007a). The slightly wider distribution for the Pleistocene Caribbean may be expected given the much greater time that is almost certainly embraced by our heavily time-averaged fossil samples. The Piacenzian (Pliocene) MART distribution is considerably wider, shows negative kurtosis and evidence of bimodality. This distribution may be explained by (1) the incorporation of samples from times when Pacific upwelling entered the SWC when the Isthmus seaway was open and also from times when modern non-upwell- ing Caribbean environments had become established after Isthmus seaway closure (O’Dea et al. 2007a), or (2) a TEP environment experiencing both ENSO events (re- sulting in low MARTs), but also years that experience both cold upwelling and warm surface waters (resulting in higher MARTs). The wide and non-normal distributions of the Zanclean (Pliocene) and the Messinean (Miocene) stages are indicative of TEP environments consistently subject to ENSO variability since the Isthmus was fully open at these times. The positively skewed distribution of the Tortonian (Miocene) stage suggests an environment in which ENSO events are dominant and relatively few years in which both cold upwelling and warm surface waters are experienced. Further Developments and Caveats.—Our results provide evidence that data from bryozoan assemblages can be used to infer variation in seasonal regimes. This new method merits further development as a proxy for inferring ENSO variability. Figure 7. Schematic representations of patterns of variation in mean annual ranges in temperature (MART) for the Caribbean over time inferred on the basis of frequency distributions of zooid- size MART estimates in the present study. BULLETIN OF MARINE SCIENCE. VOL 89, NO 4. 201316 The technique has the benefit of being based on many replicate estimates, thus mini- mizing the influence of error. Nevertheless, the approach would benefit from a more explicit determination of error and resolving power. This could be achieved by ob- taining a larger data set of modern bryozoans coupled with associated long-term, in situ, oceanographic data with analysis using randomization approaches to determine minimum sample sizes required for accurate predictions. The adoption of additional proxies of ENSO variability (e.g., growth increment and geochemical data: Rickaby and Halloran 2005, Wara et al. 2005, Ivany et al. 2011, Watanabe et al. 2011) to complement data based on zs-MART estimates would of course strengthen studies by providing independent evidence. Furthermore, stable isotope data for bryozoan material would allow assignment of absolute temperatures to the temperature ranges provided by the zs-MART values (Knowles et al. 2010). Finally, sampling sequential patterns of intra-annual variation in zooid size in long- lived bryozoans (e.g., O’Dea and Jackson 2002) could provide a means of characteris- ing complete ENSO cycles. In conclusion, our study provides a preliminary demonstration of how morpho- metric analysis of bryozoan assemblages can provide insights on the presence or ab- sence of inter-annual variation in temperature regimes, including those generated by ENSO, and therefore whether ENSO variability pertained during the Pliocene Warm Period. The approach offers a potential alternative means of characterising ENSO variability to analytical techniques, such as stable isotope analysis, which may be prohibitively expensive for characterising large numbers of individuals. When further refined and combined with other proxies, bryozoan-assemblage based zs- MART analysis should produce robust interpretations about past environmental regimes. A more thorough understanding of the dynamics of El Niño during times of climate change in the past will help improve predictions of future environmental change in our warming world. Acknowledgments The data analyzed here resulted from funding by the Natural Environment Research Council (NER/S/A/2005/13387), the University of Reading, and the Department of Zoology, Natural History Museum, London (to BO), the Natural History Museum, London (to BO and PDT), and the National System of Investigators (SNI) of the National Research of the National Secretariat for Science, Technology and Innovation of Panama (SENACYT), the National Science Foundation (EAR03-45471), and the National Geographic Society (to AO). The Panamanian “Recursos Minerales” granted collecting permits. F Rodriguez sup- ported many aspects of this research and his contribution is greatly appreciated. Literature Cited Allmon WD. 1988. Ecology of Recent turritelline gastropods (Prosobranchia, Turritellidae): current knowledge and paleontological implications. Palaios. 3:259–284. http://dx.doi. org/10.2307/3514657 Allmon WD. 1993. Age, environment and mode of deposition of the densely fossiliferous Pinecrest Sand (Pliocene of Florida): implications for the role of biological productivity in shell bed formation. Palaios. 8:183–201. http://dx.doi.org/10.2307/3515171 Amui-Vedel A-M, Hayward PJ, Porter JS. 2007. Zooid size and growth rate of the bryozoan Cryptosula palliasiana Moll in relation to temperature, in culture and in its natural envi- ronment. J Exp Mar Biol Ecol. 353:1–12. http://dx.doi.org/10.1016/j.jembe.2007.02.020 O K A M U R A E T A L . : B R Y O Z O A N P R O X Y R E V E A L S A N C I E N T E N S O VA R I A B I L I T Y 17 Arntz WE, Tarazona J. 1990. The effects of El Niño on benthos, fish and fisheries off the South American Pacific coast. In: Glynn PW, editor. Global ecological consequences of the 1982–83 El Niño-Southern Oscillation. Elsevier Oceangr Ser. 52:323–360. Atkinson D. 1994. Temperature and organism size: a biological law for ectotherms? Am Nat. 162:332–342. Atkinson D, Morley SA, Hughes RN. 2006. From cells to colonies: at what levels of body or- ganization does the ‘temperature-size rule’ apply? Evol Dev. 8:202–214. PMid:16509898. http://dx.doi.org/10.1111/j.1525-142X.2006.00090.x Aubry MP, Berggren WA. 1999. Appendix 1, New biostratigraphy. In: Collins LS, Coates AG, editors. A paleobiotic survey of Caribbean faunas from the Neogene of the Isthmus of Panama. Bull Amer Paleontol. 357:38–40. Best MMR, Kidwell SM. 2000. Bivalve taphonomy in tropical mixed siliciclastic-carbonate settings: I. Environmental variation in shell condition. Paleobiology. 26:80–102. http:// dx.doi.org/10.1666/0094-8373(2000)026<0080:BTITMS>2.0.CO;2 Clark N, Williams M, Okamura B, Smellie J, Nelson A, Knowles T, Taylor P, Leng M, Zalasiewicz J, Hayward A. 2010. Early Pliocene Weddell Sea seasonality determined from bryozoans. Stratigraphy. 7:199–206. Coates AG, Collins LS, Aubry MP, Berggren WA. 2004. The geology of the Darién, Panama, and the late Miocene-Pliocene collision of the Panama arc with northwestern South America. Geol Soc Amer Bull. 116:1327–1344. http://dx.doi.org/10.1130/B25275.1 Coates AG, Jackson JBC, Collins LS, Cronin T, Dowsett HJ, Bybell LM, Jung P, Obando JA. 1992. Closure of the Isthmus of Panama - the near-shore marine record of Costa Rica and western Panama. Geol Soc Amer Bull. 104:814–828. http://dx.doi. org/10.1130/0016-7606(1992)104<0814:COTIOP>2.3.CO;2 Coates AG, McNeill DF, Aubry M-P, Berggren WA, Collins LS. 2005. An introduction to the geology of the Bocas del Toro Archipelago, Panama. Caribbean J Sci. 41:374–391. Conroy JL, Restrepo A, Overpeck JT, Steinitz-Kannan M, Cole JE, Colinvaux PA. 2009. Unprecedented recent warming of surface temperatures in the eastern tropical Pacific Ocean. Nature Geosci. 2:46–50. http://dx.doi.org/10.1038/ngeo390 Cotton MA. 1999. Neogene planktonic foraminiferal biochronology of the southern Central American isthmus. Bull Amer Paleo. 357:1–80. Cramer D. 1997. Basic statistics for social research: step-by-step calculations and computer techniques using Minitab. London: Routledge. D’Croz L, Del Rosario JB, Góndola P. 2005. The effect of fresh water runoff on the distribu- tion of dissolved inorganic nutrients and plankton in the Bocas del Toro Archipelago, Caribbean Panama. Carib J Sci. 41:414–429. D’Croz L, O’Dea A. 2007. Variability in upwelling along the Pacific shelf of Panama and implications for the distribution of nutrients and chlorophyll. Estuar Coast Shelf Sci. 73:325–340. http://dx.doi.org/10.1016/j.ecss.2007.01.013 D’Croz L, Robertson DR. 1997. Coastal oceanographic conditions affecting coral reefs on both sides of the Isthmus of Panama. Proc 8th Int Coral Reef Symp. 2:2053–2058. Dowsett HJ, Robinson MM. 2009. Mid-Pliocene equatorial Pacific sea surface tempera- ture reconstruction: a mutli-proxy perspective. Phil Trans R Soc A. 367:109–125. PMid:18854303. http://dx.doi.org/10.1098/rsta.2008.0206 Emslie SD. 1992. The avifauna from APAC shell pit, Sarasota County, Florida. The Plio- Pleistocene stratigraphy and paleontology of southern Florida. Florida Geological Survey Special Publication. 36:179–180. Evans JL, Allan RJ. 2006. El Niño/southern oscillation modification to the structure of the monsoon and tropical cyclone activity in the Australasian region. Int J Climatol. 6:611–623. Fedorov AV, Dekens PS, McCarthy M, Ravelo AC, deMenocal PB, Barreiro M, Pacanowski RC, Philander SG. 2006. The Pliocene paradox (mechanisms for a permanent El Niño). Science. 312:1485–1489. PMid:16763140. http://dx.doi.org/10.1126/science.1122666 BULLETIN OF MARINE SCIENCE. VOL 89, NO 4. 201318 Gibbard PL, Head MJ, Walker JC. 2010. Subcommission on Quaternary stratigraphy. Formal ratification of the Quaternary System/Period and the Pleistocene Series/Epoch with a base at 2.58 Ma. J Quat Sci. 25:96–102. http://dx.doi.org/10.1002/jqs.1338 Gonzalez-Silva A, Santamaria-del-Angel E, Millan-Nunez R, Manzo-Monroy H. 2004. Satellite observations of mesoscale eddies in the Gulfs of Tehuantepec and Papagayo (eastern tropical Pacific). Deep-Sea Res II. 51:587–600. Gordon DP, Taylor PD, Bigey FP. 2009. Phylum Bryozoa. Moss animals, seamats, lace corals. In: Gordon DP, editor. New Zealand inventory of biodiversity. Volume One. Kingdom Animalia. Radiata, Lophotrochozoa, Deuterostomia. Christchurch, NZ: Canterbury University Press. p. 271–297. Hageman SJ, Needham LL, Todd CD. 2009. Threshold effects of food concentration on the skeletal morphology of the bryozoan Electra pilosa (Linnaeus, 1767). Lethaia. 42:438– 451. http://dx.doi.org/10.1111/j.1502-3931.2009.00164.x Haywood AM, Vales PJ, Sellwood BW. 2000. Global scale palaeoclimate reconstruction of the middle Pliocene climate using the UKMO GCM: initial results. Glob Planet Change. 25:239–256. http://dx.doi.org/10.1016/S0921-8181(00)00028-X Haywood AM, Valdes PJ, Peck VL. 2007. A permanent El Niño-like state during the Pliocene? Paleoceanography. 22:PA1213. http://dx.doi.org/10.1029/2006PA001323 Hsiang SM, Meng KC, Cane MA. 2011. Civil conflicts are associated with the global climate. Nature. 476:438–441. PMid:21866157. http://dx.doi.org/10.1038/nature10311 Hunter E, Hughes RN. 1994. The influence of temperature, food ration and genotype on zo- oid size in Celleporella hyalina (L.). In: Hayward PJ, Ryland JS, Taylor PD, editors. Biology and paleobiology of bryozoans. Fredensborg, Denmark: Olsen & Olsen. p. 83–86. Ivany LC, Brey T, Huber M, Buick DP, Schöne BR. 2011. El Niño in the Eocene greenhouse recorded by fossil bivalves and wood from Antarctica. Geophys Res Lett. 38:L16709. http://dx.doi.org/10.1029/2011GL048635 Jackson JBC, Todd JA, Fortunato HM, Jung P. 1999. Diversity and assemblages of Neogene Caribbean Mollusca of lower Central America. In: Collins LS, Coates AG, editors. A paleobiotic survey of Caribbean faunas from the Neogene of the Isthmus of Panama. Lawrence, Kansas: Allen. p. 193–230. Jackson JBC, O’Dea A. 2013. Timing of the oceanographic and biological isolation of the Caribbean Sea from the Tropical Eastern Pacific. Bull Mar Sci. ##:##–##. Karnauskas KB, Busalacchi AJ, Murtugudde R. 2008. Low-frequency variability and remote forcing of gap winds over the East Pacific Warm Pool. J Climate. 21:4901–4918. http:// dx.doi.org/10.1175/2008JCLI1771.1 Kaufmann KW, Thompson RC. 2005. Water temperature variation and the meteorological and hydrographic environment of Bocas del Toro, Panama. Carib J Sci. 41:392–413. Khider D, Stott LD, Emile-Geay J, Thunell R, Hammond DE. 2011. Assessing El Niño Southern Oscillation variability during the past millennium. Paleoceanography. 26:PA3222. http://dx.doi.org/10.1029/2011PA002139 Knowles TM, Leng J, Willams M, Taylor PD, Sloane HJ, Okamura B. 2010. Inferring sea- water temperature regimes using stable oxygen isotopes and zooid size variation in Pentapora foliacea (Bryozoa). Mar Biol. 157:1171–1180. http://dx.doi.org/10.1007/ s00227-010-1397-5 Knowles T, Taylor PD, Williams M, Haywood AM, Okamura B. 2009. Pliocene seasonality across the North Atlantic inferred from cheilostome bryozoans. Palaeogeog Palaeoclimat Palaeoecol. 277:226–235. http://dx.doi.org/10.1016/j.palaeo.2009.04.006 Kovats RA, Bouma MJ, Hajat S, Worrall E, Haines A. 2003. El Niño and health. Lancet. 12:917–932. Leon-Rodriguez L. 2007. Benthic foraminiferal record of the Pleistocene uplift of the sed- imentary deposits of the Burica Peninsula (Costa Rica-Panama) as a result of Cocos Ridge subduction beneath the Central American Arc. MS thesis, Florida International University, Miami, Florida, 106 p. O K A M U R A E T A L . : B R Y O Z O A N P R O X Y R E V E A L S A N C I E N T E N S O VA R I A B I L I T Y 19 Li J, Xie S-P, Cook ER, Huang G, D’Arrigo R, Liu F, Ma J, Zheng X-T. 2011. Interdecadal modulation of El Niño amplitude during the past millennium. Nat. Clim Change. 1:114– 118. http://dx.doi.org/10.1038/nclimate1086 Lombardi C, Cocito S, Occhipinti-Ambrogi A, Hiscock K. 2006. The influence of sea- water temperature on zooid size and growth rate in Pentapora fascialis (Bryozoa: Cheilostomata). Mar Biol. 149:1103–1109. http://dx.doi.org/10.1007/s00227-006-0295-3 Lunt DJ, Haywood AM, Schmidt GA, Salzmann U, Valdes PJ, Dowsett HJ. 2010. Earth sys- tem sensitivity inferred from Pliocene modelling and data. Nat Geosci. 3:60–64. http:// dx.doi.org/10.1038/ngeo706 McKinney FK, Jackson JBC. 1989. Bryozoan evolution. Boston: Unwin Hyman. Menon NR. 1972. Heat tolerance, growth and regeneration in 3 North Sea bryozoans ex- posed to different constant temperatures. Mar Biol. 15:1–11. http://dx.doi.org/10.1007/ BF00347433 Nicholls N, Lucas C. 2007. Interannual variations in area burnt in Tasmanian bushfires: relationships with climate and predictability. Int J Wildland Fire. 16:540–546. http:// dx.doi.org/10.1071/WF06125 O’Dea A. 2005. Zooid size parallels contemporaneous oxygen isotopes in a large colony of Pentapora foliacea (Bryozoa). Mar Biol. 146:1075–1081. http://dx.doi.org/10.1007/ s00227-004-1523-3 O’Dea A, Håkansson E, Taylor P, Okamura B. 2011. Environmental change prior to the K-T boundary inferred from temporal variation in the morphology of cheilostome bryozo- ans. Palaeogeogr Palaeoclimatol Palaeoecol. 308:502–512. http://dx.doi.org/10.1016/j. palaeo.2011.06.001 O’Dea A, Herrera-Cubilla A, Fortunato H, Jackson JBC. 2004. Life history variation in cupuladriid bryozoans from either side of the Isthmus of Panama. Mar Ecol Prog Ser. 280:145–161. http://dx.doi.org/10.3354/meps280145 O’Dea A, Hoyos N, Rodríguez F, De Gracia B, De Gracia C. 2012. History of upwell- ing in the tropical eastern Pacific and the paleogeography of the Isthmus of Panama. Palaeogeogr Palaeoclimatol Palaeoecol. 348–349:59–66. http://dx.doi.org/10.1016/j. palaeo.2012.06.007 O’Dea A, Jackson JBC. 2002. Bryozoan growth mirrors contrasting seasonal regimes across the Isthmus of Panama. Palaeogeogr Palaeoclimatol Palaeoecol. 185:77–94. http:// dx.doi.org/10.1016/S0031-0182(02)00278-X O’Dea A, Jackson JBC. 2009. Environmental change drove macroevolution in cupuladri- id bryozoans. Proc Roy Soc B. 276:3629–3634. PMid:19640882. PMCid:PMC2817302. http://dx.doi.org/10.1098/rspb.2009.0844 O’Dea A, Jackson JBC, Fortunato H, Smith JT, D’Croz L, Johnson KG, Todd JA. 2007a. Environmental change preceded Caribbean extinction by 2 million years. PNAS. 104:5501–5506. PMid:17369359. PMCid:PMC1838446. http://dx.doi.org/10.1073/ pnas.0610947104 O’Dea A, Okamura B. 1999. The influence of seasonal variation in temperature, salinity, and food availability on module size and colony growth in the estuarine bryozoan, Conopeum seurati. Mar Biol. 135:581–588. http://dx.doi.org/10.1007/s002270050659 O’Dea A, Okamura B. 2000a. Intracolony variation in zooid size in cheilostome bryozoans as a new technique for investigating palaeoseasonality. Palaeogeogr Palaeoclim Palaeoecol. 162:319–332. http://dx.doi.org/10.1016/S0031-0182(00)00136-X O’Dea A, Okamura B. 2000b. Life history and environmental inferences through retro- spective morphometric analysis of bryozoans: a preliminary study. J Mar Biol Ass UK. 80:3596–3599. http://dx.doi.org/10.1017/S0025315400003210 O’Dea A, Okamura B. 2000c. Cheilostome bryozoans as indicators of seasonality in the Neogene epicontinental seas of Western Europe. In: Proceedings of the 11th International Bryozoology Association Conference: 74–86. O’Dea A, Rodríguez F, Romero T. 2007b. Response of zooid size in Cupuladria exfrag- minis (Bryozoa) to simulated upwelling temperature. Mar Ecol. 28:1–9. http://dx.doi. org/10.1111/j.1439-0485.2006.00144.x BULLETIN OF MARINE SCIENCE. VOL 89, NO 4. 201320 Okamura B. 1987. Seasonal changes in zooid size and feeding activity in epifaunal colonies of Electra pilosa. In: Ross JRP, editor. Bryozoa: past and present. Bellingham, Washington: Western Washington University. p. 197–204. Okamura B, O’Dea A, Knowles T. 2011. Bryozoan growth and environmental reconstruc- tion by zooid size variation. Mar Ecol Prog Ser. 430:133–146. http://dx.doi.org/10.3354/ meps08965 Palacios DM, Bograd SJ. 2005. A census of Tehuantepec and Papagayo eddies in the northeastern tropical Pacific. Geophys Res Lett. 32:L23606. http://dx.doi. org/10.1029/2005GL024324 Ragaini L, Di Celma C, Cantalamessa G. 2008. Warm-water mollusc assemblages from northern Chile (Mejillones Peninsula): new evidence for permanent El Niño- like conditions during Pliocene warmth? J Geol Soc. 165:1075–1084. http://dx.doi. org/10.1144/0016-76492007-039 Ravelo AC, Dekens PA, McCarthy MD. 2006. Evidence for El Niño-like con- ditions during the Pliocene. GSA Today. 16:4–11. http://dx.doi. org/10.1130/1052-5173(2006)016<4:EFENLC>2.0.CO;2 Rickaby REM, Halloran P. 2005. Cool La Niña during the warmth of the Pliocene? Science. 307:1948–1952. PMid:15790852. http://dx.doi.org/10.1126/science.1104666 Scroxton N, Bonham SG, Rickaby REM, Lawrence SH, Hermoso M, Haywood AM. 2011. Persistent El Niño-Southern Oscillation variation during the Pliocene Epoch. Paleoceanography. 26:PA2215. http://dx.doi.org/10.1029/2010PA002097 Tao K, Grossman EL. 2010. Origin of high productivity in the Pliocene of the Florida plat- form: evidence from stable isotopes and trace elements. Palaios. 25:796–806. http:// dx.doi.org/10.2110/palo.2010.p10-058r von der Heydt AS, Nnafie A, Dijkstra HA. 2011. Cold tongue/warm pool and ENSO dynamics in the Pliocene. Clim Past Disc. 7:997–1027. http://dx.doi.org/10.5194/cpd-7-997-2011 Wara MW, Ravelo AC, Delaney ML. 2005. Permanent El Niño-like conditions during the Pliocene warm period. Science. 309:758–761. PMid:15976271. http://dx.doi.org/10.1126/ science.1112596 Watanabe T, Suzuki A, Minobe S, Kawashima T, Kameo K, Minoshima K, Aguilar YM, Wani R, Kawahata H, Sowa K, Nagai T, Kase T. 2011. Permanent El Niño during the Pliocene warm period not supported by coral evidence. Nature. 471:209–211. PMid:21390128. http://dx.doi.org/10.1038/nature09777 Welsh K, Elliot M, Tudhope A, Ayling B, Chappell J. 2011. Giant bivalves (Tridacna gi- gas) as recorders of ENSO variability. Earth Plan Sci Lett. 307:266–270. http://dx.doi. org/10.1016/j.epsl.2011.05.032 Yeh S-W, Kug J-S, Dewitte B, Kwon M-H, Kirtman BP, Jin F-F. 2009. El Niño in a changing climate. Nature. 461:511–514. PMid:19779449. http://dx.doi.org/10.1038/nature08316 Zhang Q, Guan Y, Yang H. 2008. ENSO amplitude change in observation and coupled mod- els. Adv Atmosph Sci. 25:331–336. http://dx.doi.org/10.1007/s00376-008-0361-5 Date Submitted: 25 May, 2012. Date Accepted: 30 April, 2013. Available Online: 8 August, 2013. Addresses: (BO, AT) Department of Life Sciences, Natural History Museum, Cromwell Road, London, SW7 5BD, United Kingdom. (AO) Smithsonian Tropical Research Institute, P.O. Box 0843-03092, Republic of Panama. (PT, AT) Department of Earth Sciences, Natural History Museum, Cromwell Road, London, SW7 5BD, United Kingdom. Corresponding Author: (BO) Email: .