Research ArticleCLIMATOLOGY

A multimillennial climatic context for the megafaunal extinctions in Madagascar and Mascarene Islands

See allHide authors and affiliations

Science Advances  16 Oct 2020:
Vol. 6, no. 42, eabb2459
DOI: 10.1126/sciadv.abb2459

Abstract

Madagascar and the Mascarene Islands of Mauritius and Rodrigues underwent catastrophic ecological and landscape transformations, which virtually eliminated their entire endemic vertebrate megafauna during the past millennium. These ecosystem changes have been alternately attributed to either human activities, climate change, or both, but parsing their relative importance, particularly in the case of Madagascar, has proven difficult. Here, we present a multimillennial (approximately the past 8000 years) reconstruction of the southwest Indian Ocean hydroclimate variability using speleothems from the island of Rodrigues, located ∼1600 km east of Madagascar. The record shows a recurring pattern of hydroclimate variability characterized by submillennial-scale drying trends, which were punctuated by decadal-to-multidecadal megadroughts, including during the late Holocene. Our data imply that the megafauna of the Mascarenes and Madagascar were resilient, enduring repeated past episodes of severe climate stress, but collapsed when a major increase in human activity occurred in the context of a prominent drying trend.

INTRODUCTION

Madagascar and the Mascarene Islands of Mauritius and Rodrigues (Fig. 1 and fig. S1) are recognized as biodiversity “hot spots”—highly threatened biogeographic regions with exceptional levels of endemic species (1). The large-scale biodiversity loss in Mascarenes, which were among the last places on Earth to be colonized by humans, occurred just within the past three or four centuries (24). Historical written accounts suggest that Mauritius lost most of its endemic species of terrestrial vertebrates within ~2 centuries of its initial colonization (c. 1638 CE), including the iconic species of the flightless bird Raphus cucullatus, popularly known as “Dodo” (4). A similar fate awaited Rodrigues, ~600 km east of Mauritius (fig. S1), where the period between the human arrival (c. 1691 CE) and the permanent colonization (c. 1790s CE) was marked by island-wide deforestation and biodiversity loss, including the extinction of Dodo’s counterpart, the Rodrigues solitaire (Pezophaps solitaria) (5). The scale of biodiversity loss in Madagascar is even more staggering. Over the past millennium, Madagascar lost virtually all its endemic megafauna (>10 kg in body mass) such as giant lemurs, elephant birds, pygmy hippopotami, and giant tortoises (69) in close temporal association with marked landscape transformations and wider evidence of increased human activities on the island (1011).

Fig. 1 Climatology and proxy locations.

Spatial pattern of December, January, and February (DJF) mean precipitation amount (millimeters per day) (shaded) between 1980 and 2015 from the European Centre for Medium-Range Weather Forecasts Reanalysis fifth-generation (ERA5) precipitation data. The climatological mean location of the southern boundary of the tropical rain belt (Sb-TRB) is indicated by the 4 mm/day isohyet (black). The precipitation amount of less than 4 mm/day is masked for clarity. The mean latitudinal positions of the northernmost and the southernmost Sb-TRB (4 mm/day isohyet) between 55°E and 75°E (dashed vertical lines) and 10°S to 25°S are shown by the red and blue isohyets, respectively (see Materials and Methods). Locations of Rodrigues (star) and other proxy records (yellow circles with numbers) (also see fig. S1) are marked as follows: 1, Anjohibe and Anjokipoty caves (10, 2629); 2, Lake Mitsinjo (25); 3, Lake Ihotry (20); 4, Lake Andolonomby (23); and 5, Ste-Luce peat bog cores (12).

The ecological devastation and biodiversity loss in Mascarenes are widely believed to have been caused by humans (25); however, whether climate change also played a role is not known with certainty because of a dearth of high-resolution paleoclimate records spanning the past few centuries. In Madagascar, identifying the triggers of its megafaunal extinction has proven challenging. The Malagasy megafaunal extinction is often attributed to human activities, either directly (via overhunting) or indirectly (via the introduction of invasive species or landscape transformation) [e.g., (610)]. In contrast, a few studies have linked the Malagasy megafauna extinction to a severe centennial-scale episode of an island-wide drought between 1000 and 900 years before the present (yr B.P.; the present is 1950 CE) that represented the culmination of a gradual millennial-scale drying trend in the region (12). However, some proxy records from the region such as oxygen isotope (δ18O) values of speleothems (1011) and nitrogen isotope (δ15N) values of extinct and extant subfossil vertebrates (1314) indicate no directional changes in hydroclimate and habitat aridity, respectively, during the late Holocene. Another set of hypotheses relate megafaunal extinctions to a complex sequence of “synergistic” processes where a period of aridification lowered the ecosystem’s resilience and amplified the human impact on the megafauna through overhunting and habitat alterations [e.g., (8)]. More recently, a “subsistence shift hypothesis” (11) posits that the megafaunal extinction may have occurred in response to a shift from a hunting/foraging to an agropastoral human economy with a concomitant increase in the human population as a consequence of expanding Indian Ocean trade networks during the late Holocene (15). The increased human population may have thus created unsustainable pressure on the megafauna population, ultimately causing its extinction (11).

All current models of the Malagasy megafaunal extinctions must contend with three sources of uncertainty. First, the timing of Madagascar’s initial human colonization is not settled, with estimates spanning nearly 9000 years—ranging from 1150 yr B.P. (16) to more than 10.5 thousand years (ka) B.P. (17). A recent critical evaluation of reported radiocarbon (14C) dates from various archaeological sites and butchered animal bones support human presence on the island by at least 2000 yr B.P. (18). Estimates based on palaeoecological, genomic, linguistic, and archaeological evidence, however, converge at 1350 to 1150 yr B.P. (16). Second, the timing and rate of megafaunal extinctions are imprecisely constrained. The current compilations of 14C dates from subfossil vertebrates suggest that most megafaunal species went extinct after ~1000 yr B.P. (9), and there are no credible 14C dates for any of the extinct species younger than ~500 yr B.P. (19). The period from ~1000 to 500 yr B.P. (~950 to 1450 CE) therefore can be viewed as the upper bound of the megafaunal “extinction window.” Estimates for the rate of extinction vary widely since the reported 14C dates on subfossil samples are unevenly distributed through time. In addition, estimates for the rate of megafaunal decline are sensitive to the choice of the “extinction model” and the timing of the initial human presence on the island. For example, if we assume that the negative relationship between megafauna and human activities commenced soon after the arrival of humans on the island at 1150 (or 2000) yr B.P. and the last dated megafaunal species died 500 yr B.P., then the entire extinction comprised 650 (or 1500) years. A recent study using 14C dates of extinct and extant subfossil vertebrates provides a more refined estimate of the megafaunal crash ranging between ~1250 and 1050 yr B.P. (~700 and 900 CE) (11). Third, disentangling the effects of climate change (e.g., aridification) and human activities (e.g., biomass burning) in proxy environmental archives from Madagascar is nontrivial. For example, the temporal variations in charcoal particulates, pollen, diatoms in lake, and peat bog sediments, as well as lake level, have been interpreted either as evidence of widespread aridification (12, 20), increased human activity (10, 21, 22), or both (2325).

Any attribution of megafaunal extinction, whether human or climate-induced, warrants careful consideration of the region’s natural climate variability at high temporal precision. Speleothems, which are secondary carbonate precipitates in caves, are extensively used for reconstructing absolute-dated and high-resolution histories of past climates and environments, typically via interpreting temporal variations of their stable oxygen (δ18O) and carbon (δ13C) isotope values and other geochemical parameters [e.g., (10)]. Several speleothem records from northwest Madagascar (fig. S1) spanning portions of the past 2000 years (10, 2630) are characterized by very large increases [~10 per mil (‰)] in δ13C but no concomitant shifts in δ18O during the late Holocene (figs. S2 and S3). These speleothem isotope records have been interpreted to reflect human-induced rather than climate-driven transformations of the island’s ecology from a dominantly C3 landscape to the present-day C4 savannah grasslands (10). However, a considerable spread in the timing, magnitude, rate, and even in the sign of δ13C shifts in speleothems from the same and a nearby cave (figs. S2A and S3B) (10, 2630) suggests complex interplays of human activities (e.g., vegetation burning/clearing above the cave), karstic processes (e.g., changes in drip water discharge, feeding routes, and ventilation), and climate (e.g., changes in precipitation and evaporation) in producing the observed changes. Critically, some of these processes such as biomass burning above the cave, either human or climate-induced, can also affect the karst vadose zone and plausibly overprint the δ18O of meteoric water in both directions (and consequently in the speleothem) by inducing changes in soil properties, water infiltration, cave ventilation, and CO2 degassing (31, 32).

In this study, we developed a record of regional hydroclimate variability over the past ~8000 years using speleothems from Rodrigues, a small volcanic island located ∼1600 km east of Madagascar (fig. S1). Although Rodrigues also underwent prominent ecological and landscape transformations coinciding with the arrival of humans and their permanent settlement, these changes occurred a few centuries after the Malagasy megafaunal extinction. Consequently, the speleothem δ13C and δ18O records from Rodrigues, at least before ~1790s CE, are not overprinted by human activities. We also report a new speleothem record from northwest Madagascar (fig. S1) that spans (~550 to 2000 yr B.P.) the critical interval of the Malagasy megafaunal extinction. We have used the instrumental, reanalysis, and simulated data from an isotope-enabled general circulation model to provide an interpretative framework of our speleothem records, which illustrates that Rodrigues and Madagascar share similar climatic variability on interannual to multidecadal time scales and that the large past shifts in regional climatic conditions manifested on both islands. We therefore suggest that our hydroclimate records from Rodrigues provide an opportunity to place the human activities and the observed ecological and biotic changes in the context of long-term natural hydroclimate variability in the region, which can be useful for a comprehensive understanding of the triggers of the late Holocene ecological and biotic changes in Madagascar and the Mascarenes.

Regional setting and modern climate

The islands of Rodrigues (~19°42′S, ~63°24′E) and Madagascar (~12°S to 26°S and 45°E to 50°E) are located in the southwest Indian Ocean (SWIO) (Fig. 1 and fig. S1). Rodrigues, a small volcanic island (~108 km2), has a drier climate with a mean annual temperature of ~24°C. The mean annual evaporation (~2085 mm) is nearly twice as high as the mean annual precipitation (~1010 mm) as the island is persistently under the influence of relatively dry southeasterly trade winds (fig. S4, A and C). The bulk of Rodrigues annual precipitation (~70%) falls during the austral summer season (November to April) when the intertropical convergence zone (ITCZ) moves south, and the island comes under the influence of the southern boundary of the tropical rain belt (Sb-TRB) (Fig. 1). Austral summer precipitation at Rodrigues is positively correlated with precipitation over much of the SWIO (fig. S4B). The island of Madagascar, off the southeast coast of mainland Africa, exhibits a complex precipitation pattern. A broad ~1200-m-high plateau, interrupted by several massifs that reach heights of up to ~2600 m, runs across the island from north to south (fig. S1). This topographic barrier creates a strong rain shadow to the west and nearly year-round heavy precipitation in the east along a narrow coastal belt due to the orographic uplift of moisture brought by persistent southeast trade winds (33). During austral summer, the ITCZ slopes southwestward from the equator reaching down to ~20°S and lies across the confluence of the northwest “Malagasy monsoon” and the easterly trade winds, creating a prominent north-to-south rainfall gradient (Fig. 1 and fig. S4A) (33). The year-to-year changes in SWIO rainfall have been shown to arise chiefly from (i) north-south shifts in the mean position of the ITCZ (33), (ii) east-west variations in the zonal position of the eastern boundary of the Mascarene High (Eb-MH) (see Materials and Methods), and (iii) changes in the sea surface temperature (SST)–induced convection associated with the El Niño–Southern Oscillation (ENSO) and the Indian Ocean Dipole (34).

We used various gridded instrumental and reanalysis climatic data to understand the relative importance of the ITCZ, ENSO, and the MH in driving the observed interannual to multidecadal changes in SWIO precipitation and to ascertain the degree to which the precipitation and drought variability at Rodrigues and Madagascar are coherent. Our analysis with the European Centre for Medium-Range Weather Forecasts Reanalysis fifth-generation (ERA5) and the 20th-century (ERA-20C) (35) datasets suggests that the December, January, and February (DJF) rainfall variability over the SWIO (~55°E to 75°E and 10°S to 25°S) is sensitive to year-to-year variations in the mean latitude of the Sb-TRB (represented by the precipitation isohyet of 4 mm/day) (see Materials and Methods). This is highlighted by Fig. 2A, which shows precipitation anomalies associated with the anomalous northward (north of ~18.7°S) minus anomalous southward (south of ~20.3°S) positions of Sb-TRB from its climatological mean location of ~19.4°S (see Materials and Methods). Time series comparisons of DJF rainfall over northern Madagascar and Rodrigues with the estimated Sb-TRB positions show a broad coherence on both interannual and decadal time scales (Fig. 2, B and D). Of note is a particularly extreme northward movement of Sb-TRB (>2σ) during 1998/1999, which was associated with large negative precipitation anomalies at both Rodrigues and Madagascar. In addition, a large decrease in precipitation over Rodrigues and Madagascar occurred between the 1920s and 1940s, in tandem with a sustained northward shift in the Sb-TRB position.

Fig. 2 Changes in Sb-TRB positions and regional precipitation.

(A) Spatial pattern of December through January (DJF) precipitation anomalies for the anomalous northward minus southward positions of the southern boundary of the tropical rainfall belt (Sb-TRB). Locations of Rodrigues (star) and other proxy records (yellow circles with numbers) are marked as in Fig. 1 (also see fig. S1). (B) Year-to-year comparisons between the mean latitudinal positions of DJF Sb-TRB (black) and ERA5 precipitation (z-scored) averaged over northern Madagascar [12°S to 20°S, 45°E to 50°E (blue)] and Rodrigues [17°S to 22°S, 60°E to 66°E (red)] (see Materials and Methods). (C) Same as (B) but for 1900 and 2009 from ERA-20C data. The correlation coefficients (r) with 95% confidence interval (brackets) among three time series are listed below (B) and (C). The mean northward and southward latitudinal positions of Sb-TRB (~17.6°S and 21.2°S) are defined by ±1σ departures (dashed lines) from the mean latitudinal position of Sb-TRB (~19.4°S) (solid line), respectively, in (B) and (C). (D) Same as (C) but for a 30-year low pass–filtered time series (see Materials and Methods). Shaded gray bar in (B) highlights the anomalous northward excursion of Sb-TRB and negative precipitation anomalies during the La Niña event of 1998 to 1999. Vertical bar in (D) highlights a sustained northward shift in Sb-TRB and changes in precipitation over Rodrigues and northern Madagascar.

We also examined the link between zonal shifts of the Eb-MH over the southern subtropical Indian Ocean region (~40°E to 105°E and ~20°S to 40°S) and the SWIO precipitation variability during austral summer. The mean longitude of the Eb-MH was identified by the 1535–geopotential meter (gpm) isopiestic contour of 850-hPa geopotential height (see Materials and Methods). Over the instrumental period (1979–2016), the SWIO exhibits negative precipitation anomalies associated with the anomalous westernmost (i.e., west of ~94°E) shift of the Eb-MH from its climatological mean location of 96.6°E (fig. S4D), which broadly resemble the spatial pattern produced by the meridional shifts in the Sb-TRB. However, the correlations between year-to-year shifts in the Eb-MH position and precipitation over Rodrigues and Madagascar are statistically insignificant (fig. S4, F and G). We further investigated the role of ENSO in influencing the SWIO rainfall variability during the instrumental period. ENSO has been shown to influence austral spring to summer rainfall over eastern Africa and SWIO via its impact on the seasonal migration of the ITCZ [which typically shifts south (north) during the El Niño (La Niña) events because of warm (cold) SST anomalies in the SWIO] (36). The DJF precipitation composites for La Niña minus El Niño years suggest a rather modest increase in precipitation over northern Madagascar and Rodrigues during the El Niño events and vice versa (fig. S4E) (37).

To assess the climatic significance of oxygen isotopes in precipitation (δ18Op) over the SWIO, and, consequently, the speleothem δ18O records, we used monthly δ18Op data from Global Network of Isotopes in Precipitation (GNIP) stations from Madagascar, Réunion, and Mauritius together with simulated δ18Op data from the Experimental Climate Prediction Center’s Isotope-incorporated Global Spectral Model (IsoGSM) (fig. S5) (38). The IsoGSM has been widely used for both modern and past climate simulations and has been shown to be in good agreement with the δ18Op observations from GNIP stations globally (38). The model-simulated δ18Op data reproduce the annual cycle in δ18Op from all three GNIP stations with reasonable skills (fig. S5A) (37). The model-simulated DJF δ18Op values extracted from the grid points nearest to Rodrigues show prominent negative correlations with both simulated and observational [Global Precipitation Climatology Project (GPCP)] precipitation amount over the SWIO region (fig. S5C). These spatial correlation patterns are similar to the patterns obtained from correlations of the GNIP δ18Op data from Mauritius (~600 km west of Rodrigues; see fig. S1) and Madagascar (fig. S5B) with the precipitation over the SWIO (fig. S5B), thus providing an important validation of the model’s data. In addition, the linear regression–derived slopes of DJF simulated δ18Op (from the grid points nearest to Rodrigues) with both model-simulated and GPCP precipitation (averaged over the SWIO region) relationships are nearly identical [m ~ −1.0‰/65 mm (r = −0.68)] (fig. S5D). The prominent inverse relationship between the δ18Op and precipitation over the convective-rich SWIO, both in the model and in observational data (fig. S5D), is consistent with theoretical considerations (39) and empirical observations from other tropical regions [e.g., (40)].

In summary, our analyses with the instrumental, reanalysis, and simulated data suggest that hydroclimate variability at Rodrigues, despite the island’s small size (~108 km2) (fig. S1), is broadly representative of the SWIO region, including Madagascar. By virtue of its location, speleothem records from Rodrigues are particularly well suited to reconstruct the region’s hydroclimate variability, particularly in response to past large-scale meridional and zonal shifts in TRB and the ENSO–MH dynamics. Furthermore, the isolated setting of Rodrigues, far removed from large-sized landmasses, and its low topographic relief (fig. S1, C and D), minimizes precipitation and speleothem δ18O variability stemming from processes such as continentality and altitude effects as well as the mixing of distant water vapor sources with substantially different isotopic compositions (37). Consequently, we have used these observations to interpret temporal variations in speleothem δ18O from Rodrigues in terms of hydroclimate variability over the SWIO.

RESULTS

Paleoclimate records

We used δ18O, δ13C, and trace element measurements from two calcite speleothems (LAVI-15-7 and LAVI-4) (data file S1) from La Vierge cave [19.75°S, 63.37°E; ~32 m above sea level (masl)] (fig. S1) in southwest Rodrigues to develop a regional record of hydroclimate variability that spans the past ~8000 years (~−65 to 8060 yr B.P.) (Fig. 3 and fig. S6). The new record supplements the stable isotope data that we previously reported from LAVI-4 for the period between ~3000 and 6015 yr B.P. (37). We also developed a new speleothem (ABC-1) record from Anjohibe cave (15.53°S, 46.88°E, ~95 masl) from northwest Madagascar (fig. S1 and data file S1) that spans the period from ~550 to 2000 yr B.P. with an average temporal resolution of ~5 years (Fig. 3). The chronologies of LAVI-15-7, LAVI-4, and ABC-1 are constrained by 10, 40, and 6 230Th dates, respectively. The age models were established by using an objective age model algorithm (see Materials and Methods, figs. S7 and S8, and data file S2) (41). Because of larger uncertainties of the U-series dates in the younger portion of the LAVI-15-7, we additionally used 13 accelerated mass spectrometry (AMS) 14C dates to constrain its chronology (see Materials and Methods, fig. S7, and data file S2). Aragonitic portions of ABC-1 were identified by x-ray diffraction (XRD) analysis and petrographic examination (see Materials and Methods and fig. S9). Coherent variations between the LAVI-4 and LAVI-15-7 δ18O records during the coeval portion (~1550 and 1980 yr B.P.) (Fig. 3 and fig. S6) as well as the coeval portions of LAVI-4 and PATA-1 from the nearby Patate cave (fig. S1D) between ~3000 and 6000 yr B.P. (37) indicate that speleothem records from Rodrigues primarily reflect climate-forced changes (see Materials and Methods). The LAVI-15-7 and LAVI-4 isotopic profiles were averaged for the period of overlap (~1550 and 1980 yr B.P.) to produce a composite record (hereafter, LAVI) (fig. S6) that provides an average temporal resolution of ~3 years with an average (2σ) temporal uncertainty of ~80 years over the past 2 ka and ~115 years over the remaining portions of the record. The LAVI δ18O and δ13C profiles were further interpolated to 5-year time steps and z-scored to facilitate their comparisons with other regional proxy records (fig. S6).

Fig. 3 The new speleothem stable isotope records from Rodrigues and Madagascar (this study).

(A) δ13C profiles of Rodrigues speleothems LAVI-15-7 (light green) and LAVI-4 (dark green) and smoothed Mg/Ca spectrum intensity ratio of LAVI-15-7 (black) together with 230Th dates and 2σ error bars (color-coded by sample names). VPDB, Vienna Pee Dee belemnite. (B) δ18O profiles of LAVI-15-7 (dark blue) and LAVI-4 (light blue). (C and D) δ13C (light brown) and δ18O (dark brown) profiles of a Madagascar speleothem ABC-1with 230Th dates and 2σ error bars. The mean isotopic values of each record are marked by dashed lines.

The LAVI-15-7 and LAVI-4 isotopic profiles coherently display a large range of δ18O (~−2 to −6‰) and δ13C (~−2 to −11‰) values (Figs. 3 and 4 and fig. S6) indicating substantial changes in the regional hydroclimate over the past eight millennia. The presence of a modest covariance between δ18O and δ13C (r = 0.6) suggests that regional changes in hydroclimate were accompanied by changes in the local hydrological balance in the catchment area of the cave’s drip water. Typically, such δ18O-δ13C coherence of stalagmites can arise from precipitation-induced changes in vegetation type, density, and soil microbial productivity such that drier conditions lead to higher isotopic values (typically δ13C lag δ18O by a few decades due to soil carbon turnover) (42). Alternatively, this covariance can also reflect climate-induced changes in karst hydrology such that drier (wetter) conditions lead to slower (faster) water infiltration, increased (decreased) pCO2 degassing, and enhanced (reduced) prior calcite precipitation (PCP)—all resulting in higher (lower) δ18O and δ13C values (43). The importance of PCP can be gauged by trace element concentrations such that strong (weak) PCP gives rise to a high (low) Mg/Ca ratio in the speleothem calcite (43). The positive correlation of the LAVI speleothem δ18O and δ13C profiles with the Mg/Ca ratio (Fig. 3A and fig. S9, E and F), together with no lag between δ18O and δ13C (fig. S9, C and D) as well as sharp drops in speleothem growth rate across major positive δ18O and δ13C excursions (fig. S9, E and F), suggests that the latter mechanism is likely the primary reason for the δ18O-δ13C coherence in our records.

Fig. 4 Comparison between Rodrigues and Madagascar speleothem stable isotope records.

(A to D) δ13C and δ18O profiles of LAVI-15-7 (Rodrigues, this study) (A), ABC-1 (Madagascar, this study) (B), AB2 [Madagascar, (26)] (C), and AB3 [Madagascar, (10)] (D) (see fig. S1 for proxy locations). All records are 5-year interpolated. Linear trends (least-squares fits of speleothem time series) that are significantly different from zero at 99% confidence level are shown with thick solid lines. Statistically insignificant linear trends are delineated by dotted lines. Plum and purple bars delineate the aragonite portions of stalagmites ABC-1 and AB2, respectively. Shaded bars on the top (with varying saturation intensity reflecting uncertainties) mark the timelines of permanent human settlements on Mauritius [red; (2, 3)], Rodrigues [blue; (5)], and Madagascar [green; (16, 18)]. The vertical color bars with varying saturation intensity show correlative features in all δ18O time series. The correlation coefficients (r) with 95% confidence interval and P values for each pair of δ18O profiles are reported on the figure. Statistical significance for each paired δ18O time series was assessed after accounting for autocorrelation and reduced degrees of freedom due to interpolation (see Materials and Methods).

Comparison with regional proxy records

Our analysis of instrumental, reanalysis, and model data allows us to propose a testable hypothesis that speleothem δ18O records from Rodrigues can be used for inferring sustained (i.e., lasting for decades or more) and large-scale (e.g., stemming from variation in the width or the mean position of tropical belt and/or ENSO dynamics) changes in past climatic conditions in the SWIO. One key test of this hypothesis entails robust correlations between the Rodrigues speleothem δ18O record and proxy records of hydroclimate variability from Madagascar as well as across the SWIO. The coeval portions of LAVI-15-7 (Rodrigues) and ABC-1(Madagascar) δ18O profiles are broadly replicated [yielding a correlation coefficient of r = 0.22 with a 95% confidence interval of (0.05, 0.37) and P < 0.01], suggesting coherent hydroclimatic variations on these islands on multidecadal to multicentennial time scales (Fig. 4). This is consistent with the observed pattern of precipitation variability during the instrumental period in response to changes in the mean position of the Sb-TRB (Figs. 1 and 2). Notably, the LAVI-15-7 and ABC-1 δ18O profiles are characterized by statistically significant increasing trends and synchronous shorter-term periods of anomalously high δ18O values (Fig. 4). These observations suggest that the SWIO has experienced a drying trend that was punctuated by shorter episodes of relatively higher aridity over the past 1.5 or 2.0 millennia. To further assess the climate fidelity of the late Holocene drying trend in the SWIO, we compared our LAVI-15-7 and ABC-1 records with published speleothem records from Anjohibe (6) and Anjokipoty (1) caves from northwest Madagascar (Fig. 4 and figs. S1, S2B, and S3A) (10, 2630). These records span various portions of the past 2 ka with highly variable temporal resolution (~1.5 to 30 years) and age constraints. Although the vast majority of δ18O values of the published records, including ABC-1, vary within a narrow range (~−4.0 to −6.5‰), no consistent pattern or trends emerge (fig. S3A). However, a comparison of δ18O profiles of ABC-1 and AB2 (26), currently the best-dated and highest-resolution Madagascar record, reveals broad replication as evidenced by a statistically significant correlation [r = 0.31, with 95% confidence interval of (0.12, 0.48), and P < 0.01] (Fig. 4). The δ18O profiles of LAVI-15-7, ABC-1, and AB2 show several correlative features, and the latter is characterized by a weak drying trend. In addition, a comparison of AB2 and ABC-1 with all other published δ18O records reveals statistically insignificant correlations (fig. S2B).

As noted in the introduction, a salient feature of nearly all published speleothem records from Madagascar (fig. S1) is a large magnitude (~8 to 10‰) increase in their δ13C values, whose onset in AB2 and AB3 speleothems occurred at ~850 to 900 CE (Fig. 4 and fig. S2) (10, 11, 26). The large positive δ13C shift has been interpreted to reflect human-driven transformations of the island-wide ecology from a dominantly C3 landscape to the present-day C4 savannah grasslands (10) and/or a major reduction in plant and soil biomass above these caves (30). This prominent δ13C signal is, however, conspicuously absent in our ABC-1 record, which shows a long-term trend, rather than a shift, toward increasing δ13C values of comparable magnitude (~8‰). A close examination of all published δ13C records from Anjohibe cave, however, shows that the onset and duration of this shift vary by as much as from 250 to 500 years, which is too large to be solely attributed to the dating uncertainties of the individual records (figs. S2A and S3B). In addition, a speleothem record from nearby Anjokipoty cave (fig. S1) (28) also shows a large δ13C shift with a comparable magnitude (~10‰) and onset (~875 CE as in AB2 and AB3 records) but in an opposite direction (fig. S2).

Divergent behaviors of Madagascar isotopic records may point to the idiosyncratic nature of speleothem formation in Anjohibe and Anjokipoty caves. Anjohibe cave, in particular, is a very large cave system containing more than 5 km of passages with at least 24 known entrances (2730). Detailed long-term cave monitoring data are absent but critically needed to assess how environmental conditions and/or karstic characteristics may vary within such a complex cave system. It is plausible that key variables such as ventilation, drip discharge, pCO2, temperature, and relative humidity (RH) may differ markedly within the cave, which may explain differences among the individual speleothem records. In addition, we cannot exclude the possibility that human-induced changes such as vegetation clearing and biomass burning may have masked the original isotopic signals through various above-cave and in-cave processes. For example, fire intensity, duration, and frequency have been shown to overprint the δ18O signature of meteoric water in both directions (and consequently in the speleothem) by inducing changes in evapotranspiration, soil properties, water infiltration, cave ventilation, and CO2 degassing (31, 32).

In summary, although no consistent picture of hydroclimate variability over the past one to two millennia has yet emerged from the available suite of speleothem δ18O records from Madagascar, we tentatively assign higher confidence to our ABC-1 δ18O record because it shows a higher degree of replication not only with the AB2 (26) but also with the LAVI-15-7 δ18O records from Rodrigues, despite the latter being located 1600 km east of Madagascar. We also note that nearly all speleothem records from Madagascar covering this period are marked by the presence of aragonite layers (2629), which typically form under drier conditions (Fig. 5).

Fig. 5 Timelines of climatic/biotic changes and human settlements.

(A) The horizontal bars with varying saturation intensity show timelines of human settlements on the Mascarene Islands (red and blue), Madagascar (green), and the Malagasy megafaunal extinction (orange). Yellow and cyan circles mark the initial dates of human settlement after (16) and (18), respectively. The dotted line shows the temporal range of human settlements on Madagascar based on contentious evidence with the purple square representing the earliest date (17). (B) The z score–transformed (interpolated to 5 years) composite δ13C (black) and δ18O (cyan) speleothem profiles from Rodrigues (this study). Periods of drier conditions are identified by z score values of >1 (vertical brown bars with varying saturation intensity). The megadroughts are identified by the z score values >1 lasting 20 consecutive years or more (asterisks). Drying trends over the past 2 and 1.5 ka are shown by red and yellow arrows. (C) Proxy-inferred drying trends and droughts from southeast (SE) (12), northwest (NW) (25), and southwest (SW) Madagascar (20, 23) are indicated by shaded bars (fig. S1). (D) The z score δ18O profiles of ANJ94-5 (blue) (29) and 60-year smoothed ABC-1 (purple, this study) from Madagascar. Bars indicate the aragonitic portions in AB2 (green) (26), ABC-1 (purple, this study), and ANJ94-5 (blue) (29). (E) Published 14C dates of extinct subfossils from Madagascar for body masses between 10 and 150 kg (orange circles) and >150 kg (green circles) (9, 13). Vertical orange and green lines mark the upper 50% confidence intervals (CIs) for island-wide reported 4C dates for 10- to 150-kg and >150-kg body-mass groups (9). (F) Pace of megafaunal decline represented by maximum likelihood χ2 statistics showing subfossil extinct-to-extant odds ratio, after (11). Estimated changes in the effective population size (purple curve) after (45). Shaded orange bar is the same as in (A).

We further compared our regional reconstruction of hydroclimate variability to two additional sources of paleoclimate information from Madagascar. Nitrogen isotope (δ15N) values in bulk collagen (13) and source amino acids (14) extracted from radiocarbon-dated extinct and extant subfossil bones of lemur species have been used to infer temporal changes in the habitat aridity. In addition to trophic level, environmental factors such as temperature, soil cycling, and land use are known to influence δ15N of soil and, by extension, δ15N of plants and animals. Moisture availability, however, appears to exert the dominant control such that animals living in drier environments typically show higher δ15N values. δ15N values of bulk collagen from several extinct and extant vertebrate genera from various subfossil sites in Madagascar reveal no directional trend in habitat aridity during the Holocene (13). Because the interpretation δ15N of bulk collagen can be confounded by the trophic level and biological processes that discriminate against 14N, a subsequent study examined δ15N of amino acid–specific fractions from one extinct and one extant lemur species and found no directional change in habitat aridity during the late Holocene (14). Trend estimation in the latter study, however, relied on a very small number of radiocarbon dates with highly uneven temporal distribution over the period of interest. In contrast, the late Holocene trend and megadroughts inferred from our LAVI-15-7 and ABC-1 records are broadly consistent with diatom, charcoal, and pollen-based reconstructions of regional climate from lake sediment cores in the southeast (12), southwest (20, 23), and northwest (25) Madagascar (fig. S1). These lake records indicate island-wide desiccation trends during the late Holocene and peak dry conditions at ~950 yr B.P. (Fig. 5). The pollen record from Lake Mitsinjo from northwest Madagascar shows a progressive desiccation trend that culminated in a hiatus, indicating a severe drought at ~1000 yr B.P. followed by evidence of arid savanna conditions between 1000 and 500 yr B.P. (25). A palaeohydrological record from Lake Ihotry suggests a lowering of the water table and salinization of coastal pans and wetlands between roughly 2250 and 700 yr B.P. (20). Palynological changes in the saline coastal basin of Ambolisatra in southwestern Madagascar also suggest aridification during the past several thousand years (23) as do contemporaneous declines in the abundance of aquatic birds at coastal subfossil sites at Belo-sur-Mer (fig. S1) (7).

DISCUSSION

Climatic context of megafaunal extinctions in the Mascarenes and Madagascar

We used the composite z score–transformed Rodrigues speleothem δ18O and δ13C records to establish a multimillennium climatic context to examine the late Holocene megafaunal extinctions and ecological changes in the Mascarenes and Madagascar (Fig. 5). Periods of “drier” and “wetter” conditions are defined by the δ18O z score values higher and lower than 1σ, respectively. The δ18O and δ13C profiles exhibit submillennial-scale hydroclimatic oscillations with a characteristic sawtooth-like pattern marked by prominent multicentury drying trends followed by sharp reversals (typically within ~10 to 100 years) to short-term pluvial conditions (Fig. 5). A salient feature of our record is that nearly all major drying trends either culminated in or were punctuated by several multidecadal to muticentennial episodes of elevated aridity (Fig. 5). We term these periods as “megadroughts” when the δ18O z score values in our record exceeded 1σ for at least 20 consecutive years. Our definition is consistent with the previous usage of this term to describe the tree ring–inferred extended periods (i.e., two decades or longer) of sustained precipitation deficits during the past one to two millennia in North America, Asia, and Europe [e.g., (44)]. Nearly all of the inferred megadroughts—for example, at ~0.2, 0.4, 0.9, 1.6, 2.8, 3.6, 4.6, and 4.9 ka B.P.—are marked by anomalously high δ18O z score values (>2σ), suggesting that they were also unusually severe in the long-term context of our reconstruction. The megadroughts centered at 2.8 and 3.6 ka B.P. were particularly long and lasted for several centuries (37). Between 8.0 and 3.7 ka B.P., nearly all episodes of drier conditions and megadroughts in our record exhibit one-to-one correspondence with the drier intervals inferred from the ANJ94-5 speleothem δ18O record from Madagascar (29) but none after 4.0 ka B.P. when the ANJ94-5 profile exhibits a substantially reduced variance (Fig. 5). In addition, several megadroughts in our reconstruction over the past 5.5 millennia are evident in a pollen-based reconstruction of major droughts from southeast Madagascar (Fig. 5) (12).

Mindful of the caveats associated with the choice of start and endpoints in estimating trends in quasi-oscillatory time series, our data suggest that the megafaunal extinctions and ecological changes in Madagascar and Mascarenes occurred in the context of a major drying trend over the past 1.5 to 2.0 millennia that was punctuated by at least four megadroughts centered at 1.6, 0.9, 0.4, and 0.2 ka B.P. (Fig. 5). This late Holocene drying trend fully encompassed the period marked by the presence of humans on Madagascar [see (16) and (18)], the increased pace of megafaunal decline (~1250 to 1050 yr B.P.), a major shift in human subsistence strategy (11), and a concomitant increase in the island population (Fig. 5) (45). These observations are in sharp contrast with the notion that ecological, biotic, and societal changes in Madagascar during the last Holocene occurred under stable climatic conditions (1014). Our data therefore necessitate a view that requires consideration of human activities in the context of a changing rather than a stable climate as the trigger of the Malagasy megafaunal extinction. Returning to the question of relative roles of human versus climate in triggering the megafaunal extinctions, lack of temporally precise records of biotic, environmental, and cultural changes in Madagascar precludes us from assessing causal links between megafaunal extinctions and climate. Our findings, however, suggest that climate change alone was not likely the trigger of Malagasy megafaunal extinctions. When viewed from a multimillennial context of our climate reconstruction, the drying trend over the past one to two millennia and the embedded megadroughts were not a unique feature of the regional climate. Millennial-scale drying trends and megadroughts occurred repeatedly, which were more severe and longer than their late Holocene counterparts (Fig. 5). And yet, the fossil records from Madagascar point to the survival of all megafauna lineages throughout much of the late Pleistocene and into the late Holocene [e.g., (9)].

As in Madagascar, our climate reconstruction suggests that the ecological transformations of the Mascarenes that coincided with the period between the first human arrival and permanent colonization of the islands also occurred during one of the driest spells of the past eight millennia (Figs. 4 and 5). The fossil records (46, 47) and the naturalists’ accounts of the island (24) suggest that all endemic species in Mauritius were extant until the permanent colonization of the island, attesting to their resilience to the past climatic extremes. This is not to say that the climate did not play a contributory role in ecological changes. For example, in the case of Rodrigues, between 1691 (c. first arrival of humans) and 1790 CE (permanent human colonization), many sailing vessels stopped at Rodrigues to collect tortoises as an abundant source of fresh meat. By one estimate, ~280,000 tortoises were taken from the island between 1735 and 1771 CE (4). With the cessation of grazing by tortoises, brush that accumulated during the long arid phase at Rodrigues may have contributed to destructive fires that removed much of the remaining forest (4).

Conclusions

We have developed a high-resolution (~5-year temporal resolution) reconstruction of the SWIO hydroclimate variability over the past ~8000 years using speleothems from the island of Rodrigues, located ∼1600 km east of Madagascar. Our analysis of the modern observational and model-simulated climate data suggests that Rodrigues and Madagascar exhibit similar climatic variability on interannual to multidecadal time scales, particularly in response to the meridional shifts in the Sb-TRB. A close correspondence between our Rodrigues and Madagascar speleothem records during the late Holocene suggests that the past shifts in the SWIO hydroclimate manifested on both islands. Our Rodrigues paleoclimate record reveals a pattern of recurring submillennial-scale aridity trends punctuated by multidecadal megadroughts throughout the past eight millennia. When viewed in a long-term context, our findings imply that the ecosystems and megafauna in Madagascar and the Mascarene Islands were resilient, surviving prior episodes of severe climate stress on the islands but collapsing when a substantial increase in human activities on the islands coincided with the late Holocene drying trend. Our data provide an important multimillennial climatic context to the megafaunal extinctions and a framework against which current and future models of megafaunal extinctions can be assessed.

MATERIALS AND METHODS

Cave locations, sample collection, and preparation

La Vierge cave (19.76°S, 63.37°E; ~32 masl, ~280 m long) is located in the François Leguat Giant Tortoise and Cave Reserve in southwestern Rodrigues (fig. S1). Stalagmite samples LAVI-4 and LAVI-15-7 were collected at ~50 and 140 m from the cave entrance from a poorly ventilated chamber during two field campaigns in August 2013 (LAVI-4) (37) and June 2015 (LAVI-15-7). The heights of LAVI-15-7 and LAVI-4 are 82 and 402 mm (fig. S7), and their mean widths are 75 and 80 mm, respectively. Both stalagmites were cut along their growth axes using a thin diamond blade and polished. The isotopic data for LAVI-4 between 81 and 274 mm (depth is relative to the top of the sample) was previously reported by Li et al. (37). In this work, we expanded the LAVI-4 record to span from 3.5 to 402.4 mm. In addition, we report new stable isotope data of the LAVI-15-7 sample between 0 and 68 mm (depth is relative to the top of the sample) (fig. S7). ABC-1 was collected from Anjohibe cave (15.53°S, 46.88°E; ~95 masl)—a large cave system in northwestern Madagascar with about 5.3 km of passages and multiple openings (2630). The length and width of the ABC-1 stalagmite are ~480 and ~200 mm, respectively. Here, we present stable isotope data corresponding to 0 to 55 mm (depth is relative to the top of the sample).

Cave monitoring

Long-term measurements of La Vierge cave temperature (T) and RH at the site of sample collection inside the cave are available between 2005 and 2015 CE. Continuous T and RH measurements (at 4-hour increments) were made with automated data loggers (HOBO) in 2007 and 2008 and (>30) spot measurements during the fieldwork season in August 2013 and June 2015. Over the period of observation, the mean T and RH values inside the cave were remarkably stable (i.e., T = 25.5 ± 0.55°C and RH > 98%) (fig. S10). Stalagmite ABC-1 was collected from Anjohibe Cave at ~1.3 km from the entrance from a poorly ventilated chamber. Temperatures at the sampling location vary between 24.5° and 26°C, approximating the mean annual temperature in the region. The RH at the time of collection was >90%, and the pCO2 varied between 400 and 500 parts per million.

230Th dating

Dating of the speleothem samples LAVI-15-7, LAVI-4, and ABC-1 was performed at Xi’an Jiaotong University, China by using Thermo-Finnigan NEPTUNE Plus multicollector inductively coupled plasma mass spectrometers. We used standard chemistry procedures (48) to separate uranium and thorium. A triple-spike (229Th–233U–236U) isotope dilution method was used to correct instrumental fractionation and to determine U/Th isotopic ratios and concentrations (49). U and Th isotopes were measured on a MasCom multiplier behind the retarding potential quadrupole in the peak-jumping mode using standard procedures. Uncertainties in U and Th isotopic measurements were calculated offline at the 2σ level, including corrections for blanks, multiplier dark noise, abundance sensitivity, and contents of the same nuclides in the spike solution. 234U and 230Th decay constants of Cheng et al. (49) were used. Corrected 230Th ages assume an initial 230Th/232Th atomic ratio of 4.4 ± 2.2 × 10−6, and the value for material at secular equilibrium with the bulk earth 232Th/238U value is 3.8. The correction for the top sample of LAVI-15-7 is large because U concentration is low [~45 parts per billion (ppb)], and the detrital 232Th concentration is elevated [>1000 parts per thousand (ppt)]. Similarly, low U concentrations (~86 ppb) or too high detrital 232Th concentrations (>47,000 ppt) lead to large errors of bottom ages of ABC-1 (data file S2).

14C dating

Three powdered subsamples (~9 mg) and 10 thin wafers (~15 mg) of stalagmite LAVI-15-7 were collected using a 0.5-mm carbide dental drill for the AMS dating. 14C dating was performed at the Beta Analytic Testing Laboratory, USA. The surface carbonate from the wafer samples was first removed by bathing the samples in 0.1 N HCl at room temperature to effectively remove all surface carbonate (potentially containing carbon from secondary sources). Washed carbonate subsamples were dried at 110°C for 14 hours and placed into an acidification vessel in a vacuum and mixed with 85% H3PO4 acid to generate CO2. The CO2 was cryogenically removed and was graphited by cobalt-catalyzed hydrogen reduction following the standard protocols for AMS counting. Oxalic acid II (National Institute of Standards and Technology Standard Reference Materials 4990C) was used as a modern reference standard with ~3‰ precision to identify the age “greater than modern” (i.e., post-1950 CE). If the fraction of modern carbon (fM) was higher than 1 due to the “bomb effect,” the percent modern carbon (pMC) was reported, and the sample age was deemed younger than 1950 CE. We used the “bomb peak” calibration in Southern Hemisphere (SH) zone 1-2 (the area between ~10°S and 90°S over the Indian Ocean) (50) and the Southern Hemisphere calibration curve (SHCal13) (51) in conjunction with Bayesian analysis (52) to estimate the ages of the top three subsamples within a 95.4% probability range (fig. S7A). When fM was lower than 1, the conventional radiocarbon age was calculated by the radiocarbon decay equation. The conventional radiocarbon ages were calendar-calibrated by using the STAlagmite dating by Radiocarbon (STAR) algorithm (53), which considers reservoir effects in karst system (fig. S7B).

Stable isotope analysis

Stable isotope (δ18O and δ13C) records of LAVI-15-7, LAVI-4, and ABC-1 were established by ~662, ~2004, and ~295 measurements, respectively. We used a computer-controlled triaxial micromill (New Wave MicroMill) to obtain subsamples (~80 μg) that were continuously milled along the growth axes of three stalagmites at increments between 100 and 200 μm. The subsamples were measured using a Finnigan MAT 253 mass spectrometer coupled with an on-line carbonate preparation system (Kiel-IV) in the Isotope Laboratory, Xi’an Jiaotong University, and an on-line carbonate preparation system (GasBench II) connected to an isotope ratio mass spectrometer (DeltaplusXL) at Innsbruck University, respectively. All results are reported in per mil (‰) relative to the Vienna Pee Dee Belemnite standard. Duplicate measurements of standards show a long-term reproducibility of 0.1‰ (1σ) or better for C and O isotopes (data file S1).

Age models

The age model of LAVI-15-7 was initially established by three age-depth curves using (i) only 230Th dates, (ii) only 14C (calibrated) dates, and (iii) a combination of 230Th and 14C dates (fig. S7). The first age model was developed with the age modeling algorithm [COnstruction of Proxy Record from Age models (COPRA)] (41) using 10 dates and assuming the top of the sample as modern (i.e., −65 yr B.P. or 2015 AD, the year the stalagmite was collected). The second age-depth curve was established by using only radiocarbon dates using the STAR age model algorithm (fig. S7B) (53). References for bomb peak calibration in SH zone 1-2 (the area between ~10°S and 90°S over the Indian Ocean) (50) and SHCal13 (51), calibrated radiocarbon ages of subsamples from the surface and at 0.4-mm depth, were obtained using BetaCal 3.21 and the high-probability density range method (52). For other 14C dates, an anchor point of known age (720 ± 34 yr B.P.) at 15.5 mm (blue triangle), the LAVI-15-7 δ13C profile, and SHCal13 (51) were used to construct the age model. The third age model used both 230Th dates and calibrated 14C dates that were modeled using COPRA (41). The results of this study are not sensitive to the choice of age model used. Therefore, we used the COPRA-derived age model with only 230Th dates including an estimated age of −65 yr B.P. at 0 mm as the final age model for LAVI-15-7 (fig. S7E). The age model of LAVI-4 was modeled by COPRA. The age model of ABC-1 was initially established by using COPRA and was subsequently refined further by tuning it to Rodrigues composite records within the 2σ age error. Seven peaks in the Rodrigues composite δ18O records were selected as targets, and seven corresponding counterparts were used as anchors for an age model based on linear interpolation (fig. S8).

Composite Rodrigues records

The δ18O and δ13C values of Rodrigues stalagmites LAVI-15-7 and LAVI-4 were used to produce composite long-term records for the interval 8060 to −65 yr B.P. Both stalagmites records were first z scored, then interpolated into 5-year time steps, and averaged over the period of overlap to yield a composite record (Fig. 5 and fig. S6).

Isotopic equilibrium tests

Two lines of evidence suggest that speleothem calcite at La Vierge cave precipitated at or near isotopic equilibrium. First, the coeval portions of LAVI-4 and LAVI-15-7 δ18O records replicate well between 1550 and 1980 yr B.P. (r = 0.48, P < 0.001) (Fig. 3 and fig. S6), suggesting that both stalagmites primarily record climate variability. In addition, the LAVI-4 δ18O record also exhibits a high degree of replication with the PATA-1 δ18O record from the nearby Patate cave between ~3000 and 6000 yr B.P. (37). Second, the theoretical range of calcite δ18O values [e.g., (54)] in La Vierge cave calculated using the observed cave ambient temperature of 25° to 26.5°C and δ18O values of drip waters (−3.6 to −4.1‰, n = 4) are consistent with the δ18O values of modern calcite of LAVI-15-7 (fig. S10). δ18O measurements of drip water were conducted at the Stable Isotope Biogeochemistry Laboratory, Institute of Earth Environment, Chinese Academy of Sciences, China. All water isotope samples were filtered using a 0.22-μm filter and then analyzed on an Isotopic Liquid Water Analyzer Picarro 2130-i equipped with an A0211 vaporizer operated at 110°C. The water isotope values are reported in per mil (‰) relative to the Vienna standard mean ocean water.

XRD and trace element analysis

Five powdered subsamples (~10 mg) from ABC-1 were drilled for XRD analysis using 0.5-mm carbide drill bits. Analyses were performed on an x-ray diffractometer (SmartLab, Rigaku) at the Frontier Institute of Science and Technology, Xi’an Jiaotong University, China, operated at 40 kV and 30 mA with CuKα radiation.

Trace elements

ABC-1 trace element (Mg/Ca and Sr/Ca) records were established by ~202 data. The LAVI-15-7 and LAVI-4 trace element (Mg/Ca) records were established by ~264 and 1030 data, respectively. Trace element measurements were conducted on a laser-induced breakdown spectroscopy system at the Isotope Laboratory, Xi’an Jiaotong University. Briefly, the polished slab of stalagmite was fixed on the x-y-z stage. An yttrium-aluminum-garnet-Nd laser (Beamtech, Dawa-200) was used to generate a 1064-nm laser using a lens focusing the 5-mm beam diameter to 0.1 mm. Emitted plasma from the stalagmite surface was collected by optical fibers and sent to a four-passage spectrometer (Ocean Optics MX500+) to obtain a spectrum within the 200- to 580-nm range. The intensity of emission lines, including Mg (I, 285.2 nm), Sr (II, 407.8 nm), and Ca (II, 373.7 nm), was used to calculate the ratios of elements. Ca was used as an internal standard to correct for system offset. Time series of Mg/Ca and Sr/Ca were obtained by averaging ratios of 20 laser shots at each subsample position after five cleaning shots. The measurements were performed continuously along the stalagmites’ growth axes at 300-μm increments.

ABC-1 aragonite correction

Aragonitic portions of the ABC-1 stalagmite were identified by the antiphase relationship between Mg/Ca and Sr/Ca along the growth axis (fig. S9). Details of this method are described in (26). Identified aragonitic layers were confirmed by petrographic examination and XRD analysis (fig. S9). The calcite-equivalent δ18O values of aragonitic portions were determined by applying corrections of −1.7‰ for δ13C and −0.75‰ for δ18O.

TRB and MH analysis

The location of the southern limb of TRB in the SWIO was identified by the mean latitude of the 4 mm/day isohyet of ERA5 (35) within the area of 55°E to 75°E and 10°S to 25°S for DJF between 1980 and 2015. We used linear interpolation to identify the latitude where the mean precipitation was equal to 4 mm/day following the method described in (55). A time series of the mean latitude of Sb-TRB for DJF was derived, and the first and third quartile values (18.7°S and 20.3°S) represent anomalous northward (2010, 1994, 1988, 2001, 2000, 1980, 1987, 1992, and 1998) and southward (2014, 1982, 2006, 1986, 2011, 1989, 1984, 2002, and 2003) positions of the Sb-TRB (Figs. 1 and 2). The location of Eb-MH was identified by the easternmost longitude of the 1535-gpm isopiestic contour at 850-hPa geopotential height for DJF between 1979 and 2016 following the method described in (56). Typically, the Eb-MH is located at ~96.6°E. The first and third quartile values (94.1° and 100.3°E) were used as the markers to find nine anomalous westward (1984, 2011, 1983, 2010, 1982, 1985, 1989, 1991, and 2000) and eastward (1997, 2015, 2009, 1987, 1994, 1979, 2001, 2006, and 1993) locations of the Eb-MH (fig. S4D).

Data used for climate analysis

The annual cycles of precipitation and temperature at Madagascar (16°S to 15°S, 46°E to 47°E), Mauritius (21°S to 20°S, 57°E to 58°E), and Rodrigues (20°S to 19°S, 63°E to 64°E) (fig. S4C) between 1980 and 2015 were derived from the monthly Climatic Research Unit (CRU) TS4.03 dataset (57). This land-only dataset is at 0.5° resolution and covers the period between 1901 and 2018. No CRU data are available for Rodrigues between 1901 and 1950. The mean DJF precipitation time series between 1980 and 2015 over northern Madagascar (12°S to 20°S, 45°E to 50°E) and Rodrigues (17°S to 22°S, 60°E to 66°E) were extracted from the monthly ERA5 reanalysis dataset (Fig. 2B) and ERA-20C (Fig. 2, C and D). The DJF Sb-TRB position is defined by the mean latitude of 4 mm/day isohyet within a region delineated by 55°E to 75°E and 10°S to 25°S. ERA5 and ERA-20C (35) are global reanalysis datasets for atmosphere, land surface, and oceans with a resolution of 0.25° and 1.4°, respectively. The comparisons of time series derived from ERA5 and ERA-20C with the GPCP v2.3 dataset (58) between 1980 and 2015 show statistically significant positive correlations (r = 0.76 and 0.87, respectively). For ERA-20C and GPCP, the correlation coefficients of mean DJF precipitation at northern Madagascar and Rodrigues between 1979 and 2009 are 0.72 and 0.64, respectively. The correlation coefficients (with 95% confidence intervals) between the Sb-TRB positions and rainfall amount over northern Madagascar and Rodrigues for 1980 and 2015 are 0.57 (−0.76, −0.30) and −0.68 (−0.82, −0.45), respectively. The extreme northward shift of the Sb-TRB (11.7°S) in 1998 coincides with a record-low precipitation amount in both northern Madagascar and Rodrigues. The correlation coefficient between the Sb-TRB positions and precipitation amount for the ERA-20C dataset (1900 and 2009) over northern Madagascar and Rodrigues are −0.20 (−0.37, −0.01) and −0.68 (−0.77, −0.56), respectively.

GNIP δ18Op data

We used monthly δ18Op data from three nearby GNIP stations in Madagascar (18.9°S, 47.6°E, 1295 masl), Réunion (20.9°S, 55.5°E, 70 masl), and Mauritius (20.3°S, 57.5°E, 402 masl) (source: www.iaea.org/services) to characterize the regional δ18Op variability. The δ18Op data from Mauritius, Madagascar, and Réunion comprise 7, 8, and 15 years, and each station contains 77, 87, and 180 monthly δ18Op observations with gaps from 2009 to 2015 CE, 2009 to 2016 CE, and 2002 to 2016 CE, respectively. We find strong correlations between the observed and IsoGSM-simulated annual δ18Op cycles at Madagascar, Réunion, and Mauritius with correlation coefficients of 0.89, 0.93, and 0.71, respectively (fig. S5A). Furthermore, the DJF amount-weighted IsoGSM δ18Op over Rodrigues (17°S to 22°S, 60°E to 66°E) and the regional precipitation amount averaged over the area (26°S to 12°S, 43°E to 80°E) are inversely related with both the model-simulated precipitation [r = −0.63 (−0.79, −0.39), n = 37] and GPCP precipitation [r = −0.68 (−0.82, −0.45), n = 37] amount, illustrating the regional amount effect (fig. S5D). The slope of DJF amount-weighted δ18Op-precipitation amount relationship yields a conversion of 1‰ increase in δ18Op to a decrease of ~0.74 mm/day (corresponding to ~66 mm over DJF) in rainfall (fig. S5).

Statistical analysis

We applied a parametric nonlinear regression technique to constrain the timing of initial changes in each δ13C profile of published Madagascar records (fig. S2A). This method uses a continuous function consisting of two linear parts that are joined at the breakpoint. The break model was fitted to data using a weighted least-squares method with a brute-force search for the breakpoint. Statistical uncertainties of the trend, timing, and breakpoint were evaluated using 2000 block bootstrap simulations, which preserved the distribution and serial dependence of the data over the length of a block. Low-pass filtering of the LAVI composite records was conducted by using a sixth-order Butterworth filter (fig. S6). Cross-wavelet analysis was used to find localized intermittent periodicities of each time series and explore their common power and relative phase in time-frequency space (fig. S9, C and D). The statistical significance of linear trends and correlation coefficients (95% confidence intervals and P values) reported in this paper, where appropriate, account for autocorrelation of time series.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/42/eabb2459/DC1

https://creativecommons.org/licenses/by-nc/4.0/

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank the staff at François Leguat Giant Tortoise and Cave Reserve at Rodrigues for help with speleothem collection and logistics. Funding: This work was supported by the National Natural Science Foundation of China grants (41888101, 41731174, and 41561144003 to H.C.; 41472140 to Y.N.; and 41703007 to G.K.), National Science Foundation (grant 1702816 to R.L.E. and H.C.), the State Key Laboratory of Loess and Quaternary Geology, Institute of Earth Environment, CAS (grant SKLLQG1414 to Y.N.), and the Chinese Academy of Sciences “PIFI Program” (grant 2020VCA0019 to A.S.). Author contributions: A.S. designed the project and wrote the manuscript with H.L. C.S., Y.N., H.L., P.D., and J.Z. conducted the stable isotope measurements. H.C., R.L.E., G.K., X.L., L.S., and H.L. conducted the radiometric dating. A.S., A.A.A., H.C., A.M., J.B., Y.N., H.L., and G.K. conducted the field work. All authors discussed the results and provided input on the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper, in the Supplementary Materials, and at the National Climatic Data Center (https://ncdc.noaa.gov/). Request for materials should be addressed to A.S. (asinha{at}csudh.edu) or H.C. (cheng021{at}xjtu.edu.cn).
View Abstract

Stay Connected to Science Advances

Navigate This Article