letters to nature
running and when the locking loop is closed. The closed-loop frequency ﬂuctuations exhibit a root-mean-square (r.m.s.) frequency deviation of ,310 kHz over a 1-min integration time. The residual noise exhibits a white noise spectrum with an Allan variance of 3 £ 10211/t 1/2 for an averaging time of 1 s # t # 25.6 s, and is limited by the electronic noise of our detection system (PD3 in Fig. 4a). A straightforward improvement would be to use a saturable absorption line (which can have a linewidth as narrow as 1 MHz; ref. 13) as a reference. The improved signal-to-noise ratio in HC-PCF also makes overtone absorptions in the visible and near-infrared (for example, P11 of 12C2H2 at 790.703 nm) accessible to laser frequency metrology. To test for gas leakage, the pressuredependent linewidth of the P9 acetylene absorption was monitored daily over two months. The results gave an average value of ,468 MHz with a standard deviation of 30 MHz, that is, the linewidth is consistent with the Doppler limited value and was constant within experimental error, indicating no measurable leakage of gas. In conclusion, all-ﬁbre, ultra-compact, high performance, easyto-use and unconditionally stable gas-laser devices have been reported. The commercial availability of a wide range of all-ﬁbre components (for example, lasers, phase modulators, power attenuators, isolators, Bragg gratings and beam-splitters) makes complex systems easy to design and construct. It is now possible to imagine miniature laser–gas devices, occupying a tiny volume and containing minute amounts of gas, with everyday applications in ﬁelds such as the colour conversion of laser light (perhaps using a built-in diode pump laser), and the measurement and stabilization of laser frequency. The unique features of HC-PCFs make these gas–laser devices extremely efﬁcient, and their impact in many laser-related ﬁelds is likely to be deep and lasting. A
15. MacAdam, K. B., Steinbach, A. & Wieman, C. A narrow band tunable diode laser system with grating feedback, and a saturated absorption spectrometer for Cs and Rb. Am. J. Phys. 69, 1098–1111 (1992). 16. Wieman, C. E. & Hollberg, L. Using diode lasers for atomic physics. Rev. Sci. Instrum. 62, 1–20 (1991). 17. Onae, A. et al. Optical frequency link between an acetylene stabilized laser at 1542 nm and an Rb stabilized laser at 778 nm using a two-color mode-locked ﬁber laser. Opt. Commun. 183, 181–187 (2000).

Acknowledgements The support of the UK Engineering and Physical Sciences Research Council (EPSRC) is acknowledged. F.B. thanks A. Luiten and M. Maric for discussions on frequency stabilization. F.B. is an EPSRC Advanced Research Fellow. Competing interests statement The authors declare that they have no competing ﬁnancial interests. Correspondence and requests for materials should be addressed to F.B. (pysab@bath.ac.uk).

..............................................................

Obliquity pacing of the late Pleistocene glacial terminations
Peter Huybers1 & Carl Wunsch2
1 Woods Hole Oceanographic Institution, Woods Hole, Massachusetts 02543, USA 2 Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA

.............................................................................................................................................................................

Methods
Splicing procedure
Fusion-splicing was carried out using a commercial ﬁlament-based splicer (Vytran FFS-2000-PM). In this splicer, the splicing region is continuously purged by argon gas to stop the ﬁlament burning. This prevents contamination of the splice by solid deposits and water condensation, and also prevents combustion of ﬂammable gases such as hydrogen and acetylene. For high ﬁll pressures, the splice could be consolidated using heat-curable glue.

Threshold measurements
The energy threshold measurement for the generation of the ﬁrst Stokes line S00(1) was carried out in the following manner. The input power was varied by manual rotation of a half-wave plate, and the output light was collected and fed to an optical spectrum analyser and power detectors. The detected signal at the Stokes wavelength was monitored as the input power was increased up to a point where the Stokes signal emerged abruptly from the noise on the optical spectrum analyser. With this technique, we found a relative uncertainty varying between 40% and 70%.
Received 3 November 2004; accepted 6 January 2005; doi:10.1038/nature03349.
1. 2. 3. 4. 5. 6. Luiten, A. N. (ed.) Frequency Measurement and Control (Springer, Berlin, 2001). Hall, J. L. Frequency stabilized lasers—a parochial review. Proc. SPIE 1837, 2–15 (1992). Hentschel, M. et al. Attosecond metrology. Nature 414, 509–513 (2001). Steinmeyer, G., Sutter, D. H., Gallmann, L., Matuschek, N. & Keller, U. Frontiers in ultrashort pulse generation: Pushing the limits in linear and nonlinear optics. Science 286, 1507–1512 (1999). Benabid, F., Knight, J. C., Antonopoulos, G. & Russell, P. S. J. Stimulated Raman scattering in hydrogen-ﬁlled hollow-core photonic crystal ﬁber. Science 298, 399–402 (2002). Benabid, F., Bouwmans, G., Knight, J. C., Russell, P. S. J. & Couny, F. Ultra-high efﬁciency laser wavelength conversion in gas-ﬁlled hollow core photonic crystal ﬁber by pure stimulated rotational Raman scattering in molecular hydrogen. Phys. Rev. Lett. 93, 123903 (2004). Russell, P. Photonic crystal ﬁbers. Science 299, 358–362 (2003). Marcuse, D. Loss analysis of single-mode ﬁber splices. Bell Syst. Tech. J. 55, 703–718 (1977). Dushman, S. Scientiﬁc Foundations of Vacuum Technique (John Wiley & Sons, New York, 1962). Hansch, T. W. in The Hydrogen Atom (ed. Hansch, T. W.) 93–102 (Springer, Berlin, 1989). Huber, A. et al. Hydrogen-deuterium 1S–2S isotope shift and the structure of the deuteron. Phys. Rev. Lett. 80, 468–471 (1998). Quinn, T. J. Practical realization of the deﬁnition of the metre, including recommended radiations of other optical frequency standards (2001). Metrologia 40, 103–133 (2003). Nakagawa, K., de Labachelerie, M., Awaji, Y. & Kourogi, M. Accurate optical frequency atlas of the 1.5 mm bands of acetylene. J. Opt. Soc. Am. B 13, 2708–2714 (1996). Swann, W. C. & Gilbert, S. L. Pressure-induced shift and broadening of 1510–1540-nm acetylene wavelength calibration lines. J. Opt. Soc. Am. B 17, 1263–1270 (2000).

7. 8. 9. 10. 11. 12. 13. 14.

The 100,000-year timescale in the glacial/interglacial cycles of the late Pleistocene epoch (the past ,700,000 years) is commonly attributed to control by variations in the Earth’s orbit1. This hypothesis has inspired models that depend on the Earth’s obliquity (,40,000 yr; ,40 kyr), orbital eccentricity (,100 kyr) and precessional (,20 kyr) ﬂuctuations2–5, with the emphasis usually on eccentricity and precessional forcing. According to a contrasting hypothesis, the glacial cycles arise primarily because of random internal climate variability6–8. Taking these two perspectives together, there are currently more than thirty different models of the seven late-Pleistocene glacial cycles9. Here we present a statistical test of the orbital forcing hypothesis, focusing on the rapid deglaciation events known as terminations10,11. According to our analysis, the null hypothesis that glacial terminations are independent of obliquity can be rejected at the 5% signiﬁcance level, whereas the corresponding null hypotheses for eccentricity and precession cannot be rejected. The simplest inference consistent with the test results is that the ice sheets terminated every second or third obliquity cycle at times of high obliquity, similar to the original proposal by Milankovitch12. We also present simple stochastic and deterministic models that describe the timing of the latePleistocene glacial terminations purely in terms of obliquity forcing. To test whether the glacial variability is related to changes in Earth’s astronomical conﬁguration, we adopt a formal nullhypothesis (H0) that glacial terminations are independent of obliquity variations, and the alternative hypothesis (H1) that glacial terminations are paced by it. Our focus on obliquity is motivated by previous indications of nonlinear interactions between obliquity period and quasi-100-kyr glacial variability13, but we also perform identical tests for pacing by precession and eccentricity. The test is focused on glacial terminations because their magnitude and abruptness facilitate accurate identiﬁcation. Several obstacles must be overcome to distinguish between H0 and H1. A major problem is the need to establish time controls on the glacial variability. Many studies estimate age by assuming a relationship between climate proxy variability and orbital
491

NATURE | VOL 434 | 24 MARCH 2005 | www.nature.com/nature

© 2005 Nature Publishing Group

letters to nature
forcing14,15, but this approach assumes the validity of the hypothesis being tested. Instead, we use an age-model devoid of orbital assumptions and apply it to the leading empirical orthogonal function (EOF1) of ten well-resolved marine d18O records13, a proxy for ice volume and ocean temperature (Fig. 1a). Most simple models of the late Pleistocene glacial cycles have at least four degrees of freedom2, and some have as many as twelve3. Unsurprisingly then, the seven observed quasi-100-kyr glacial cycles are insufﬁcient to distinguish between the skill of the various models16. Models with minimal degrees of freedom are necessary. Other requirements for a useful test include the ability to cope with noisy records, age-model uncertainty and (possibly) nonlinear interactions. Here we test for stability in the phase of the orbital parameters during glacial terminations using Rayleigh’s R (see Methods). To proceed, we must estimate the probability distribution function associated with H0. Of the ﬁve estimation methods explored, the one adopted gives the highest critical value and thus makes H0 the most difﬁcult to reject—a modiﬁed random walk8 representing ice-volume variability: V tþ1 ¼ V t þ ht and if V t $ T o ; terminate: ð1Þ This highly simpliﬁed model posits 1-kyr steps in ice volume, V t , of random length, h t , independently drawn from a normal distribution with standard deviation j ¼ 2 and mean m ¼ 1. The nonzero mean biases the Earth toward glaciation. Once V t reaches a threshold, To ¼ 90, a termination is triggered, and ice-volume is linearly reset to zero over 10 kyr. If the model were deterministic with j ¼ 0 and m ¼ 1, glacial cycles would last exactly 100 kyr, but with j ¼ 2, glacial cycle duration is approximately normally distributed at (100 ^ 20) kyr, a spread consistent with observations11. Initial ice volume is randomly set between 0 and To with equal probability. From a Monte Carlo technique (see Methods), we ﬁnd a critical value of R ¼ 0.60. The observed obliquity phases produce R ¼ 0.70, and H0 is rejected (Fig. 1b). This rejection of H0 is robust to all plausible reformulations of the test. Thus, the phase of obliquity has a statistically signiﬁcant relationship with the timing of deglaciation. The mean phase at deglaciation is indistinguishable from zero and is associated with maxima in obliquity. We estimate the H 1 probability density function by assuming that terminations always initiate at the same phase of obliquity, but that termination timing is subject to identiﬁcation and age-model uncertainties (see Methods). The maximum-likelihood value of H1 is R ¼ 0.69, very near the observed value, further supporting the conclusion that glacial terminations are paced by variations in obliquity. Hypotheses not accounting for the obliquity pacing are unlikely to be correct. Analogous hypothesis tests for precession and eccentricity produce different results (Fig. 1c, d). Age-model uncertainty approaches half a precession cycle, so that the power of the precession test is negligible—even if present, precession pacing of the glacial cycles cannot be discerned. In the case of eccentricity, H0 is not rejected using the random-walk probability estimate (equation (1)), but is rejected using weaker formulations of the eccentricity null hypothesis. The discrepancy arises because null

Figure 1 The Rayleigh test for phase directionality. a, d18O EOF1 normalized so that negative values indicate more ice. Dots indicate onset of a termination and horizontal bars indicate one-standard-deviation age-model uncertainties13. Termination 3 is split between events 3a and 3b. Vertical lines indicate the time of maxima in obliquity. Note that time runs from right to left in the palaeoclimate convention. b, The obliquity phase (dots) sampled at each termination and plotted on a unit circle. The vector average has a magnitude R ¼ 0.70 (cross mark) exceeding the critical value c ¼ 0.60 (ﬁlled circle), so that H0 is rejected. Furthermore, R is near H1’s maximum-likelihood value (dashed circle). The direction is indistinguishable from maximum obliquity (top of the circle). c, d, Analogous tests are made for precession (R ¼ 0.43, c ¼ 0.60) (c) and eccentricity (R ¼ 0.66, c ¼ 0.84) (d), but in neither case can the corresponding H0 be rejected. See the Supplementary Information for more details.
492

Figure 2 Deterministic and stochastic descriptions of the late-Pleistocene glacial variability. a, Deterministic model results (red) with an obliquity-dependent threshold (black) plotted over EOF1 (brown). b, Periodograms of the deterministic model results (red) and EOF1 (brown). Concentrations of energy are centred on the 1/41-kyr obliquity frequency and the 1/100-kyr glacial band; as well as combination tones at 1/70, 1/29 and 1/23 kyr. The approximate 95% conﬁdence interval is indicated by the vertical bar on the right. c, A realization of the stochastic model. d, Histogram of the time between terminations, derived from many runs of the stochastic model. The observed duration between terminations (triangles, using termination 3a not 3b) coincide with the dominant 80- and 120-kyr modes. e, Histogram of Rayleigh’s R from the stochastic model with the observed obliquity value, R ¼ 0.70, indicated by the triangle.
NATURE | VOL 434 | 24 MARCH 2005 | www.nature.com/nature

© 2005 Nature Publishing Group

letters to nature
hypotheses that assume a glacial timescale of roughly 100 kyr (which seem more physical) tend to have higher R values and are more difﬁcult to reject. Because the hypotheses of negligible inﬂuence of precession and eccentricity on the glacial terminations cannot be rejected, we adopt a minimalist strategy, retaining only obliquity to describe the glacial terminations. From a physical standpoint, support for the obliquity control hypothesis also comes from the fact that maxima in obliquity cause annual average insolation anomalies of up to 10 Wm22 at high latitudes. Furthermore, the annual average and seasonal insolation redistributions associated with obliquity are hemispherically symmetric—as is the glacial variability to within a few thousand years17,18. Finally, although obliquity control does not by itself address the question of why termination ﬁve is anomalously large, it does alleviate the marine isotope stage 11 problem2 of explaining the initiation of termination ﬁve when eccentricity and precession amplitude are small. But how does a forcing with a 40-kyr period pace the 100-kyr latePleistocene glacial variability? One suggestion is that the phase and amplitude modulations of obliquity cause the 100-kyr variability4, but it is difﬁcult to see the climatic signiﬁcance of such modulations. Instead, we suggest that the climate state skips one or two obliquity beats before deglaciating, thus giving quantized glacial-cycle durations of either 80 or 120 kyr. A speculative scenario is for increased obliquity to increase high-latitude insolation and cause heating of an ice sheet, eventually warming the ice–bedrock interface. When the ice sheet is thin, basal temperature and pressure are low19, and the obliquity heating has little effect—a skipped beat. But when the ice sheet is thick, basal temperature and pressure are high19, and additional obliquity heating promotes basal melting. Enhanced lubrication of the ice–bedrock interface may trigger deglaciation by increasing ice ﬂux into the ocean or toward lower latitudes, as well as by increasing the thinning rate and causing inward migration of the ablation zone20. Note that ,10 kyr are required for surface heating to penetrate to the base of an ice sheet19. Unlike precession, changes in obliquity are associated with sustained annual average changes in insolation21,22 and are thus more likely to cause basal warming. Furthermore, if climate was warmer during the early Pleistocene, terminations would be triggered more nearly every obliquity cycle, giving the observed smaller and more rapid glacial variability23. A simple system, consistent with the above observations, is obtained by making the threshold in equation (2) dependent on obliquity:  ð2Þ V tþ1 ¼ V t þ ht and if V t $ T o 2 avt ; terminate:  Here vt is obliquity24, normalized to zero mean and unit variance, with amplitude a. The time-variable threshold makes it more likely for a termination to occur when obliquity is large. Both deterministic and stochastic variants of equation (2) are offered, to emphasize that such models are not theories of climate change, but rather attempts at efﬁcient kinematic descriptions of the data, and that different mechanisms can be consistent with the limited observational records. To make the model deterministic, the ice-volume step size in equation (2) is ﬁxed at h ¼ 1. Setting a ¼ 15, To ¼ 105 and (initial ice volume) V (t¼2700) ¼ 30 provides a good description of the late-Pleistocene glacial variability (Fig. 2a). Many other parameterizations yield similar results—we choose this one because it is particularly simple. The model produces the correct timing of the glacial terminations (using termination 3a, not 3b) and has a squared-cross-correlation with d18O EOF1 of 0.65, an excellent ﬁt considering there are only three adjustable parameters. (Parameters other than h could be adjusted instead.) For comparison, tuning other models having four2 or twelve3 adjustable parameters yields a maximum squared-cross-correlation with EOF1 of 0.24 and 0.74
NATURE | VOL 434 | 24 MARCH 2005 | www.nature.com/nature

respectively25. If skill is measured by squared-cross-correlation divided by the number of adjustable parameters, equation (2) does more than three times as well. A periodogram of the deterministic model results (Fig. 2b) shows narrow-band concentrations of energy at the average 100-kyr period, the 41-kyr obliquity period, and at previously identiﬁed13 combination tones 1/41 2 1/100 ¼ 1/70, 1/41 þ 1/100 ¼ 1/29, and 1/41 þ 2/100 ¼ 1/23 kyr—in good agreement with the d18O EOF1 periodogram. The appearance of 1/23-kyr narrow-band energy, in the absence of precession-band forcing, highlights the ambiguous origins of this energy band22. The model also has energy at 2/100 kyr, not visible in the observations. Note that the model produces an energetic background continuum consistent with the d18O periodogram even though it is deterministic. For the stochastic case, h t is deﬁned to be normally distributed with j ¼ 2 and m ¼ 1. Here h t represents the unpredictable background weather and climate variability spanning all timescales out to the glacial/interglacial. All other parameter settings are kept the same as in the deterministic case. Now, equation (2) resembles an orderone autoregressive process. Such a process is known to describe the glacial variability well, except for runs associated with glacial terminations26, modelled here using the threshold condition in equation (2). The time between terminations in the stochastic model averages 100 kyr but has a tri-modal distribution with maxima at two (80 kyr), three (120 kyr), and four (160 kyr) obliquity cycles (Fig. 2d). The observed durations between terminations are consistent with the dominant 80- and 120-kyr modes, but suggest that one would see 160-kyr glacial cycles, given a larger sample size. The R value obtained by the stochastic model averages 0.85 (Fig. 2e), higher than the observed R ¼ 0.70; this is as expected because of observational age-model error. Deterministic and stochastic variants of equation (2) thus both describe the late-Pleistocene glacial variability. A description of the early Pleistocene variability23 is obtained by lowering the termination threshold to To ¼ 40, giving smaller-amplitude terminations that occur more nearly every obliquity cycle. Alternatively, instead of specifying a parameter change, the mid-Pleistocene transition can be described using a chaotic model25 (not shown) having spontaneous transitions between 40 and 100-kyr modes of glacial variability. With only seven cycles, it is unclear whether adequate data will ever be available to distinguish between stochastic, simple deterministic and chaotic deterministic models of the glacial variability. The simplest interpretation of our results is that during the Pleistocene epoch, the Earth tends to a glacial state (anthropogenic inﬂuences aside) and deglaciates near some, but not all, obliquity maxima. Obliquity control of the timing of deglaciation, probably in concert with a stochastic forcing, has several other consequences. These include inferences drawn from precession-based and eccentricity-based models of glaciation27, the hemispheric symmetry of glacial cycles, and the efﬁcacy of age-model tuning of cores. These issues will be taken up elsewhere. A

Methods
Rayleigh’s R is deﬁned as28: R¼    N 1 X   cos fn þ i sin fn    N  n¼1 ð3Þ

Here, f n is the phase of obliquity stroboscopically sampled at the nth glacial termination. The line brackets indicate the magnitude, making R real and non-negative, with a maximum value of one when the phases are all equal. Relative to measures of phase coupling used to investigate other nonlinear interactions29, Rayleigh’s R requires many fewer phases (roughly ﬁve) to test for phase stability28. To measure R between obliquity and the glacial cycles, we use d18O EOF1 and an agemodel independent of orbital assumptions13 (Fig. 1a). Such an independent age-model is important because even so-called minimal tuning strategies—using only a narrow band of a climate record—tend to align the abrupt glacial terminations with a particular phase of the assumed forcing (as indicated by Monte Carlo simulations), thus biasing records towards H1.

© 2005 Nature Publishing Group

493

letters to nature
The rate of change of EOF1 has a unimodal distribution with a long tail, indicative of rapid melting. Terminations are deﬁned to initiate when the rate of change in EOF1 ﬁrst exceeds the two-standard-deviation level. This criterion identiﬁes each of the usual terminations10,11, and two events in the termination 3 deglacial sequence (termed 3a and 3b for the younger and older events, respectively). Additional rules could be added to exclude 3a or 3b, but this rejection seems ad hoc, and so we use all eight termination events. Note, the timing of termination 3 predicted by the Paillard model3 also coincides with either event 3a or 3b, depending on slight changes in the parameterizations. Reassuringly, results are not sensitive to details of the test: either termination 3a or 3b can be excluded; termination times can be deﬁned using the midpoint between the local minimum and maximum bracketing each termination; and individual benthic or planktic d18O records15 can be used in place of EOF1. The probability density function (PDF) associated with H0 is estimated using the modiﬁed random-walk model (equation (1)). A realization of R is obtained by sampling the phase of obliquity at eight consecutive termination initiations, generated from equation (1), and the PDF of H0 is estimated by binning 104 such realizations of R. Other methods are to assume a uniform phase distribution, use surrogate data techniques30, or to derive statistics from ensemble runs of other models, but all of these give PDF estimates which make H0 more easily rejected and are therefore not used. To estimate the PDF associated with H1, we assume that glacial terminations always occur at the same phase of obliquity, but that the phase observations are subject to identiﬁcation and age-model error. A Monte Carlo technique is used where the timing of the glacial terminations are perturbed according to the estimated age uncertainties13 (these average ^ 9 kyr) and identiﬁcation error (^1 kyr, the EOF1 sampling resolution). A realization of R is then computed using the phase of obliquity at the perturbed ages, and 104 such realizations are binned to estimate the PDF of H1. We estimate the likelihood of correctly rejecting H0 (that is, the power of the obliquity test) to be 0.57. See the Supplementary Information for a listing of the other pertinent data and statistics.
Received 6 August 2004; accepted 24 January 2005; doi:10.1038/nature03401.
1. Hays, J., Imbrie, J. & Shackleton, N. Variations in the earth’s orbit: Pacemaker of the ice ages. Science 194, 1121–1132 (1976). 2. Imbrie, J. & Imbrie, J. Modeling the climatic response to orbital variations. Science 207, 943–953 (1980). 3. Paillard, D. The timing of Pleistocene glaciations from a simple multiple-state climate model. Nature 391, 378–391 (1998). 4. Liu, H. Insolation changes caused by combination of amplitude and frequency modulation of the obliquity. J. Geophys. Res. 104, 25197–25206 (1999). 5. Gildor, H. & Tziperman, E. Sea ice as the glacial cycles’ climate switch: Role of seasonal and orbital forcing. Paleoceanography 15, 605–615 (2000). 6. Saltzman, B. Stochastically-driven climatic ﬂuctuations in the sea-ice, ocean temperature, CO2, feedback system. Tellus 34, 97–112 (1982). 7. Pelletier, J. Coherence resonance and ice ages. J. Geophys. Res. 108, doi:10.1029/2002JD003120 (2003). 8. Wunsch, C. The spectral description of climate change including the 100ky energy. Clim. Dyn. 20, 353–363 (2003). 9. Saltzman, B. Dynamical Paleoclimatology: Generalised Theory of Global Climate Change (Academic, San Diego, 2002). 10. Broecker, W. Terminations. in Milankovitch and Climate (eds Berger, A. et al.) Part 2, 687–698 (D. Riedel, Hingham, 1984). 11. Raymo, M. E. The timing of major climate terminations. Paleoceanography 12, 577–585 (1997). 12. Milankovitch, M. Kanon der Erdbestrahlung und seine Andwendung auf das Eiszeiten-problem (Royal Serbian Academy, Belgrade, 1941). 13. Huybers, P. & Wunsch, C. A depth-derived Pleistocene age-model: Uncertainty estimates, sedimentation variability, and nonlinear climate change. Paleoceanography 19, doi:10.1029/ 2002PA000857 (2004). 14. Imbrie, J. et al. in Milankovitch and Climate (eds Berger, A. et al.) Part 1, 269–305 (D. Riedel Publishing Company, 1984). 15. Shackleton, N. J., Berger, A. & Peltier, W. R. An alternative astronomical calibration of the lower Pleistocene timescale based on ODP site 677. Trans. R. Soc. Edinb. Earth Sci. 81, 251–261 (1990). 16. Roe, G. & Allen, M. A comparison of competing explanations for the 100,000-yr ice age cycle. Geophys. Res. Lett. 26, 2259–2262 (1999). 17. Blunier, T. & Brook, E. Timing of millennial-scale climate change in Antarctica and Greenland during the last glacial period. Science 291, 109–112 (2001). 18. Wunsch, C. Greenland-Antarctic phase relations and millennial time-scale climate ﬂuctuations in the Greenland cores. Quat. Sci. Rev. 22, 1631–1646 (2003). 19. Marshall, S. & Clark, P. Basal temperature evolution of North American ice sheets and implications for the 100-kyr cycle. Geophys. Res. Lett. 29, doi:10.1029/2002GL015192 (2002). 20. Zwally, H. et al. Surface melt-induced acceleration of greenland ice-sheet ﬂow. Science 297, 218–222 (2002). 21. Rubincam, D. Insolation in terms of earth’s orbital parameters. Theor. Appl. Climatol. 48, 195–202 (1994). 22. Huybers, P. & Wunsch, C. Rectiﬁcation and precession-period signals in the climate system. Geophys. Res. Lett. 30, doi:10.1029/2003GL017875 (2003). 23. Raymo, M. & Nisancioglu, K. The 41 kyr world: Milankovitch’s other unsolved mystery. Paleoceanography 18, doi:10.1029/2002PA000791 (2003). 24. Berger, A. & Loutre, M. F. Astronomical solutions for paleoclimate studies over the last 3 million years. Earth Planet. Sci. Lett. 111, 369–382 (1992). 25. Huybers, P. On the Origins of the Ice Ages: Insolation Forcing, Age Models, and Nonlinear Climate Change. PhD thesis, MIT (2004). 26. Wunsch, C. Quantitative estimate of the Milankovitch-forced contribution to observed quaternary climate change. Quat. Sci. Rev. 23, 1001–1012 (2004). 27. Ruddiman, W. F. Oribital insolation, ice volume, and greenhouse gases. Quat. Sci. Rev. 22, 1597–1622 (2003). 28. Upton, G. & Fingleton, B. Spatial Data Analysis by Example Vol. 2 (John Wiley and Sons, Chichester, 1989). 29. Rosenblum, M. & Pikovsky, A. Synchronization: from pendulum clocks to chaotic lasers and chemical oscillators. Contemp. Phys. 44, 401–416 (2003). 30. Schreiber, T. & Schmitz, A. Surrogate time series. Physica D 142, 346–382 (2000).

Supplementary Information accompanies the paper on www.nature.com/nature. Acknowledgements Useful comments were provided by E. Boyle, W. Curry, T. Herbert, J. McManus, F. Ng, M. Tingley and G. Yang. P.H. is supported by the NOAA Postdoctoral Program in Climate and Global Change and C.W. is supported in part by the National Ocean Partnership Program (ECCO). Competing interests statement The authors declare that they have no competing ﬁnancial interests. Correspondence and requests for materials should be addressed to P.H. (phuybers@whoi.edu).

..............................................................

Two episodes of microbial change coupled with Permo/Triassic faunal mass extinction
Shucheng Xie1,2, Richard D. Pancost1, Hongfu Yin2, Hongmei Wang3 & Richard P. Evershed1
1 Organic Geochemistry Unit, Bristol Biogeochemistry Center, School of Chemistry, University of Bristol, Cantock’s Close, Bristol BS8 1TS, UK 2 Faculty of Earth Science and State Key Laboratory of Geological Processes and Mineral Resources, China University of Geosciences, Wuhan 430074, People’s Republic of China 3 School of Environmental Studies, China University of Geosciences, Wuhan 430074, People’s Republic of China

.............................................................................................................................................................................

Microbial expansion following faunal mass extinctions in Earth history can be studied by petrographic examination of microbialites (microbial crusts)1,2 or well-preserved organic-walled microbes3. However, where preservation is poor, quantiﬁcation of microbial communities can be problematic. We have circumvented this problem by adopting a lipid biomarker-based approach to evaluate microbial community changes across the Permo/Triassic (P/Tr) boundary at Meishan in South China. We present here a biomarker stratigraphic record showing episodic microbial changes coupled with a high-resolution record of invertebrate mass extinction. Variation in the microbial community structure is characterized by the 2-methylhopane (2-MHP) index (a ratio of the abundance of cyanobacterial biomarkers to more general bacterial biomarkers). Two episodes of faunal mass extinction were each preceded by minima in the 2-MHP index, followed by strong maxima, likely reﬂecting microbial responses to the catastrophic events that caused the extinction and initiated ecosystem changes. Hence, both cyanobacterial biomarker and invertebrate fossil records provide evidence for two episodes of biotic crisis across the P/Tr boundary. In bacteria, bacteriohopanepolyols (BHPs) are amphiphilic biochemicals that serve regulatory and rigidifying functions in cell membranes. It has been demonstrated that the cyanobacterial lineage currently comprises the only known quantitatively signiﬁcant source of 2-methyl-BHPs4–6, such that the presence of abundant 2a-methylhopanes (2-MHPs) derived from 2-methyl-BHP precursors is believed to be characteristic of cyanobacteria4. This is particularly true for 2-methyl-BHPs containing more than 32 carbon atoms, which have not been detected in organisms other
NATURE | VOL 434 | 24 MARCH 2005 | www.nature.com/nature

494

© 2005 Nature Publishing Group