Low-frequency orbital variations controlled climatic and environmental cycles, amplitudes, and trends in northeast Africa during the Plio-Pleistocene

The eastern Mediterranean sapropels, paced by insolation, provide a unique archive of African monsoon strength over the Late Neogene. However, the longer-term climate of this region lacks characterization within the context of changes in ice volume, sea surface temperature gradients, and terrestrial ecosystems. Here, we examine C 28 n -alkanoic acid leaf wax hydrogen and carbon isotopes in sapropels, sourced from northeast Africa, along with vegetation-corrected precipitation isotopes, derived from astronomically dated sediment cores from ODP 160 Sites 966 and 967 since 4.5 million years ago. Despite sampling only wet-phase sapropels for African monsoon variability, we ﬁ nd a larger range in hydrogen isotopes than previously published data across wet-dry precession cycles, indicating the importance of long-term modulation of Green Sahara phases throughout the Neogene. An in ﬂ uence of orbital properties on regional monsoonal hydroclimate is observed, controlling up to 50% of total hydrogen isotope variance, but large changes outside of these typical frequencies account for at least 50% of the total variance. This secular trend may track changes in ice volume, tropical sea surface temperature, sea surface temperature gradients, or even lower-frequency orbital cycles. Long-term hydroclimate and environmental shifts provide new contexts for milestone events in northeast African hominin dispersal and evolution.

T he climate and environment of Africa have fluctuated substantially over the Pliocene and Pleistocene, drastically altering the landscape where our early human ancestors lived, moved, and evolved.Northeast Africa has served as the terrestrial pathway for migration from the African continent to the rest of the world, particularly during Green Sahara intervals, which allowed for movement across this otherwise barrier-like desert.Global, regional, and local processes have been shown to contribute to the patterns of climate variability near this region in Africa by way of glacial-interglacial cycling 1 , monsoon precipitation 2 , and land surface feedbacks 3 .Records of terrestrial climate processes in northeast Africa are central for contextualizing the environments of our hominin ancestors.
Quantitative, high-resolution records of late Pleistocene climate and environmental change highlight northeast Africa as a potential pathway for Homo sapiens to disperse out of Africa 4 .Other paleoclimate records from this same time interval in the eastern Africa region suggest that periods of specific climate conditions are particularly suited for hominin stability and dispersal 5 .Comprehensive modeling efforts link periods of highamplitude, orbitally driven temperature, rainfall, and primary productivity changes to hominin dispersal 6 , supporting evolutionary hypotheses that focus on orbitally driven climate drivers, such as the variability selection hypothesis 7,8 .Reconstructions of the climate and ecosystem change in this region are needed for disentangling the various influences and potential processes involved.However, climate and environmental variability and trends have been difficult to quantify over long time intervals due to the lack of continuous, geochemical proxy records measured in long, well-dated archives.
New quantitative reconstructions are required to assess patterns of trend, amplitude, and periodicity in order to disentangle potential climate drivers.Geochemical reconstructions from welldated sediment cores allow for continuous measurement, which is ideal for the application of time series analysis to tease apart these patterns.Leaf wax records from both marine 4,9 and lacustrine 10,11 sediment cores document mixed signals with features of both glacial-interglacial cycles and the regional monsoon system in the late Pleistocene, highlighting the need for wider application of these approaches.Prior to the late Pleistocene, leaf wax hydrogen (δD wax ) and carbon (δ 13 C wax ) isotope records have also shown direct evidence for the eccentricity-modulation of precession cycles 10,12 .While monsoon strength has been shown to dominate the orbital-scale variability of the earlier Plio-Pleistocene [10][11][12][13][14] , these terrestrial records are often too short to be able to reconstruct climate variability over millions of years and thus capture long-term trends and a wide range of variability frequencies.Further, near-tropical, terrestrial (including lacustrine) sediment cores are often very difficult to date, limiting the application of advanced spectral analyses to decipher periodicities and climate drivers.Recent compilations of terrestrial climate records have extended this view and highlighted the importance of the 400-kyr eccentricity cycle, but the standardization of the included records limits the temporal resolution and, thus, the range of studied frequencies 15 .
To characterize the climate and environment during monsoon maxima of this region central for hominin evolution and dispersal, we present records of precipitation and vegetation from leaf wax biomarkers preserved in sapropel layers from two adjacent Mediterranean Sea sediment cores (Fig. 1) that together span the last 4.5 million years.The sapropel layering is astronomically tuned, thus providing a long, continuous, precisely dated archive from which to reconstruct climate over the Plio-Pleistocene, while geochemical tracers from terrestrial northeast Africa preserved in the cores allow for the quantification of trends, amplitudes, and spectral properties of climate and environmental change.These novel records, along with comparisons to higher-resolution X-ray fluorescence (XRF) analyses and indices, such as the Wet-Dry Index 16 , will address questions of climate drivers and hominin migration through the Plio-Pleistocene.

Background
Seasonal monsoon strength fluctuation is driven by the changing distribution of solar radiation, or insolation, linked to the wobble, or precession, of Earth's orbit.Orbital precession drives temperature and pressure gradient differences between ocean and land, strengthening and weakening the monsoon system over 21 kyr cycles.Eccentricity, or the degree of circularity of Earth's orbit around the sun, varies on 100 kyr and 400 kyr cycles and modulates the amplitude of precession variability at a given latitude, but has little effect on total insolation.Obliquity, the tilt of the axis of the Earth, modulates the seasonal distribution and latitudinal gradient of insolation, and is thought to mainly contribute to climate variation at the polar regions.These orbital parameters contribute to the varying strength of the seasonal monsoon system, which can, in turn, drive changes in vegetation and terrestrial ecosystems 10,12,17,18 .
On the other hand, global climate during the Plio-Pleistocene is characterized by the gradual strengthening of the glacialinterglacial cycles [19][20][21] .The cooling trend is thought to trigger changes in the amplitude of the oscillations, as well as changes in the periodicity of the cycles from precession-to obliquity-to eccentricity-paced during this time 19,22 .The benthic foraminiferal oxygen isotope stack 19 documents these deep-water temperature and ice volume patterns, but the characteristic Plio-Pleistocene cycles have far-reaching effects in the climate system.Ice volume has been shown to correspond with ocean circulation, atmospheric CO 2 , temperature, and sea surface temperature gradients, all of which may influence the climate conditions and ecosystems of continental Africa 23,24 .
Canonical records of terrestrial African paleoclimate, such as the soil carbonate carbon and oxygen isotope compilations of eastern Africa 25 , document gradual trends towards more open, arid environments.However, these outcrop measurements often represent very fine temporal and spatial scales due to the formation processes, so while it is possible to resolve long-term trends by compiling many individual measurements, analyses of temporal variability are limited.Other evidence of high latitude processes driving long-term trends in African climate, such as pollen 26 and leaf waxes 27,28 , have been documented in lowresolution reconstructions from marine sediment cores off the coast of Africa, which integrate larger areas and temporal scales.Higher-resolution reconstructions from these relatively long marine sediment cores, such as dust fluxes 29,30 and other lithological proxies 16,31,32 are able to resolve the temporal scale needed for variability analyses and document the periodicity changes mimicking the Plio-Pleistocene glacial-interglacial cycles.However, recent analyses suggest the monsoon cycles do not follow the Plio-Pleistocene pattern observed in ice volume/deep sea temperature 33,34 , adding to the need for further study.
There is ample evidence for regional connectivity of precipitation variability over the Plio-Pleistocene.The eastern African rift system carries evidence of synchronous, pulsed deep lake intervals across the basins that correspond to the 400 kyr eccentricity cycle 35 , which is prominent in low-latitude insolation 36 .The Mediterranean sapropel layers are also a classic indicator of the wet-dry cycles that North Africa experienced over the Plio-Pleistocene.Sapropel deposition is caused by the enhanced discharge from the Nile River Delta and other North African inputs to the Mediterranean following strong monsoonal rains over eastern 37,38 and western 39,40 Africa.This fresh water provides a low-density cap to the Mediterranean Sea that prevents mixing and oxygenation deep water, leading to deposition of organic-rich muds during these wet intervals.These strata have been linked to orbital precession 38,[41][42][43][44] and indicate that it is highly likely that wet-dry cycles in this region were driven by both the East African Monsoon (EAM) and West African Monsoon (WAM) systems, modulated dominantly by low-latitude insolation.However, it is difficult to discern the long-term trends, as well as the amplitude of precession-band monsoon intensity in these binary (i.e., deep/shallow lake and dark/light sediment) archives.
The eastern Mediterranean Sea is bound by land masses encompassing a large range of different climate and environmental zones.Provenance studies using lignin biomarkers 14 and mineralogy 45 have indicated that a vast majority of the depositional material comes from the Nile River Valley and the northeastern Sahara Desert (Fig. 1).Leaf waxes from Northeast Africa are also likely transported by aeolian and fluvial pathways to the Mediterranean Sea, capturing the climate and environmental histories across biomes 46 covering the Ethiopian Highlands, the Nile River Valley, and the northeastern Saharan Desert when vegetated, which are seasonally affected by the EAM, and at its maximum extent, the WAM.
The eastern Mediterranean is notable for sapropel formation due to its dynamic relationship with the Nile River.These sediment layers are visually distinctive, marked by dark, organic rich muds embedded in light, nannofossil oozes.Outcrops of these strata dating back to the early Pliocene occur in Sicily and elsewhere, and sapropels are ubiquitous in deep Mediterranean Sea cores.Ocean Drilling Program (ODP) Leg 160 Sites 966 and 967 47 , located on the northern margin of the Eratosthenes Seamount in the Levantine Basin (34°N, 34°E, 927 m and 2253 m water depth, respectively; Fig. 1), recovered a long record of sapropel deposition covering the past 4.5 Myr.Multiple holes at each site provide continuous composite depth scales that capture the full sapropel record (Supplementary Table 1).Due to the lack of preservation of sapropel sediment in ODP Site 967 prior to 3.2 Ma, as denoted by the Barium to Aluminum ratio (Fig. 2e), cores from both Sites 966 and 967 were necessary for leaf wax analyses.
Linking sapropel deposition to orbital cycles has been the basis of high-precision age models in the region 16,[43][44][45] .Here, we build upon previous detailed absolute 48 and tuned 16,45 age models by refining the characterization of sapropel layers in the bottom sections of cores from Sites 967 Holes A and B and Site 966.We use the Laskar 36 orbital solution for up-to-date ages for insolation cycles.Thorough tuning via AnalySeries 49 of the density and reflectance parameters between the cores provides composite age and depth models going back 4.5 Myr (Fig. S1).All climate time series and subsequent interpretations are based off this ODP 966/ 967 composite record, hereafter ODP 966/967.The updated age model largely follows previous ones 16,45 , with maximum changes on the order of 0.077 Myr compared with ref. 45 in the oldest interval of ODP Site 966.

Results
Leaf wax distributions and isotope values.The average chain length (ACL) of our samples is 27.2 and the average Carbon Preference Index (CPI), is 4.7 (minimum 1.7), similar to sedimentary samples with significant leaf wax n-acid contributions.Leaf wax δ 13 C values are significantly correlated between compound homologues (r C26-C28 = 0.73, p C26-C28 < 0.0001, n = 119; r C28-C30 = 0.68, p C28-C30 < 0.0001, n = 121).There are no systematic offsets in isotope value between the long chain compounds, although the amount of variance within C 30 and C 26 are higher than the C 28 δ 13 C wax dataset.Following the correction for δ 13 C CO2 , the n-C 28 δ 13 C wax reconstruction varies between −28.5‰ and −21.5‰, corresponding to 54% and 6% C 3 plant contributions, respectively, although with potentially large, propagated error bounds.The δ 13 C wax record does not contain a significant trend (r = 0.07, p > 0.1, n = 120).δ 13 C wax is correlated to CPI (r = −0.52,p < 0.0001, n = 145), with more C 3 -like values correlating to higher CPI values, but the isotope record is also correlated with various chain length indicators, such as ACL (r = 0.47, p < 0.0001, n = 145) and long:short chain lengths (r = 0.41, p < 0.0001, n = 145), indicating that longer compounds signify more grass-rich environments despite potential variations in the source of the waxes.
We observe significant correlations among δD wax values of long-chain n-acids (r C26-C28 = 0.75, p C26-C28 < 0.0001, n = 120; r C28-C30 = 0.63, p C28-C30 < 0.0001, n = 123), suggesting leaf wax homologues in the samples are derived from a common vascular plant source.Similar to δ 13 C wax , the structure of δD wax between C 26 , C 28 , and C 30 are very similar with no systematic offset.Thus, we use n-C 28 to represent the long-chain δD wax , which following Fig. 2 Eastern Mediterranean leaf wax isotope records of hydroclimate and vegetation.Leaf wax isotope data generated for this study from ODP Sites 967 (post 3.2 Ma) and 966 (pre 3.0 Ma), including carbon isotopes with inferred %C 3 and 2 standard error mean (SEM) shaded (a), hydrogen isotopes 2 SEM shaded (b), δD of precipitation derived from the vegetation correction of δD wax (c), the eccentricity of Earth's orbit 36 (d), and the ratio of Barium to Aluminum in the ODP Site 967 core to chemically denote the sapropel deposition 54 (e) .the sea water δD correction has a mean of −150.1‰ and varies between −128.9‰ and −180.9‰(Fig. 2b).The δD wax record contains a weak, but significant, D-enrichment trend (r = 0.22, p < 0.0001, n = 118).δD wax is not correlated with ACL (r = −0.12,p > 0.1, n = 147) nor CPI (r = −0.10,p > 0.1, n = 147), indicating the lack of influence that compound distribution has on the hydrogen isotope values.Using the δ 13 C wax for a vegetation correction of δD wax , the δD precip record has a mean of −41.9‰ and varies between −17.4‰ and −80.3‰ (Fig. 2c).
Spectral analysis of hydrologic timeseries.Both hydrogen and carbon leaf wax isotope datasets show ample variability throughout the last 4.5 Myr.Lomb-Scargle spectral analysis of the δD precip record displays strong spectral power at the 400 kyr and 938 kyr bands (Fig. 3a).Considering the width of the lowfrequency peak at half maximum (~860 kyr to 1.02 Myr; Fig. 3a), hereafter, we refer to this as a 1-Myr cycle.These strong periodicities are fairly consistent throughout the 4.5 Myr (with buffers of ~800 kyr on each end), yet the 1-Myr signal is more prominent before 2.7 Ma, and the 400-kyr signal is more prominent after 2.7 Ma (Fig. 3b).Gaussian 400-kyr-and 1-Myr-band-pass filters of δD precip supports the spectral analysis results and are in phase with similarly filtered records of insolation (Fig. 4).When the dominant 400-kyr and 1-Myr cycles are removed from the δD precip record, the residuals trend towards more D-enriched values with statistical significance (r = −0.34,p < 0.001, n = 114; Figs.5a and 6a).The 400 kyr band comprises 10% of the variance of the δD precip record, and the 1-Myr band comprises of 15%, with the residual representing 75% of the variance derived from non-orbital cycles.
Although we constrain the range of frequencies used as indicators of climate forcing (See Methods), a 2227 kyr cycle (close to the 2.3-Myr eccentricity cycle) appears to contribute to about a quarter of the total variance of the δD precip record, reducing the contribution of non-orbital cycles to 50% (Fig. 6b).When all three of these cycles have been filtered out of the δD precip record, there is no significant trend (Fig. 6b), implying that the low-frequency 2.3-Myr cycle may be the source of the enrichment trend over the last 4.5 Myr (Fig. 6a).

Discussion
Cycles: the importance of low-frequency orbital periodicities.The presence and periodicity of sapropelic mud deposition exposed in Eastern Mediterranean outcrops and sediment cores highlights the importance of orbital cycles in the climate of the region.There is ample evidence for dominant precession control on monsoon intervals and sapropel occurrence, with some records also documenting obliquity and eccentricity signals.Reconstructions of hydroclimate 4,[11][12][13][14]50 and lake level 34,[51][52][53] in eastern Africa also document control by orbital precession and eccentricity. Moeling exercises further highlight the importance of orbital precession to the regional monsoon system in eastern Africa 2 . In theor, low-frequency orbital modulation of precession and obliquity signals should be apparent in long monsoon records.However, it has remained difficult to quantify the effects of low-frequency cycles due to the lack of long, well-dated records.
The δD precip reconstructed here captures low frequency orbital cycles and their effect on monsoon strength during wet intervals.The dominant 400 kyr-band and 1 Myr-band cycles we find in northeast African precipitation (Fig. 3a) demonstrates that these low frequency cycles contribute significantly to modulation of the monsoon system.Spectral analysis demonstrates that these orbital cycles are robust in the δD precip record and impact the magnitude of the wet endmembers found in the sapropel archives.Further, the specific timing and amplitude of these cycles in ODP 966/967 correspond extremely well with the same cycles observed in the low-latitude summer insolation timeseries (Fig. 4).The in-phase relationship of the δD precip and insolation cycles suggests direct insolation forcing of the monsoon system through mechanisms that modify the strength of the largely precession-paced monsoon cycles.X-ray fluorescence (XRF) analyses and the associated ratios, such as Ti/Al, attributed to the relative proportions of aeolian (Ti) to riverine (Al) inputs 54,55 , and indices, like the Wet-Dry Index from principal component analysis (PCA) in ref. 16 , also demonstrate strong, persistent 400-kyr cyclicity (Fig. 3c), supporting the inference of monsoon precipitation variability at this timescale.
The strength and persistence of the 400-kyr eccentricity cycles in the climate record are notable because eccentricity has minimal influence over the absolute amount of insolation received on the Earth's surface, instead acting as an amplifier of the redistribution of energy by orbital precession (Fig. 3e).Through this Fig. 4 Filtered leaf wax dataset compared with insolation and deep lake intervals.Gaussian filtering of δD precip from this study (blue) and 30°N June 21 st insolation (gray) at the 400-kyr (a) and 1-Myr (b) bands, with exact frequencies noted in Methods.Eastern African deep lake intervals from ref. 56 and ref. 35 highlighted in gray.Fig. 5 Potential drivers of non-orbital trend and variability.The residual δD precip variability from ODP 966/967 after removal of 400 kyr and 1 Myr cycles (a).This residual record is compared with potential external drivers including (b) the benthic stack 19 , (c) the tropical sea surface temperature stack 66 , and (d) Pacific SST gradient (west minus east) based off of Mg/Ca (solid lines) and U K' 37 (dashed lines) 67 .1-Myr moving average of records are documented in orange to demonstrate long-term trends and fluctuations.mechanism, 400-kyr modulation of precession maxima leads to the observed negative relationship of δD precip of sapropel wet intervals with low-latitude summer insolation (Fig. 4a).More positive insolation maxima drive stronger monsoon rainfall leading to more negative δD precip values.Supporting this interpretation, similarities between the timing of the 400-kyr wet intervals in the δD precip record generally correspond to intervals of wide-spread deep lakes in eastern Africa 35,56 (Fig. 4).However, the continuous δD precip record brings new detail to these long-term fluctuations and supports the idea that these lake basins likely experienced extreme dry intervals not fully represented in the geologic record.While the precession and 100-kyr eccentricity cycles have been linked with human evolutionary changes 8,12 , there is a lack of statistical correlation between environmental variability in the 400-kyr cycle and speciation or extinction in Africa 15 .
The 1-Myr cycle that also dominates the spectral power within ODP 966/967 δD precip is potentially driven by low-frequency obliquity modulation with a 1.2 Myr periodicity, although this lies outside of the range of frequencies noted by spectral analysis.Obliquity and its modulation have a more direct influence on insolation at polar, rather than tropical latitudes.However, the 1.2 Myr cycle has been recorded in marine paleoclimate records that are of sufficient resolution and length 57 , and other studies find 41-kyr-band obliquity cycle in the Mediterranean sapropel records 42,55 , although it is relatively insignificant within the Wet-Dry Index 16 record (Fig. 3c, d).The presence of obliquity in lowlatitude paleoclimate reconstructions, including δD precip in this study, suggests that other orbital parameters, such as the crossequatorial insolation gradient 58 , could be important modulators of tropical precipitation strength, which may vary on an intracontinental scale.It is possible that the 1-Myr cycle is derived from currently unknown processes.
Through evolutive spectral analysis (Fig. 3b, d), we visualize how and when the strength of orbital cycles within the hydroclimate records increases and decreases through time.Within the δD precip record from ODP 966/967, the 400 kyr and 1 Myr cycles are consistently strong since 4.5 Ma except when there is a shift from a strong 1 Myr cycle to a strong 400 kyr cycle around 2.7 Ma, and when the 400-kyr cycle weakens at 1 Ma (Fig. 3b).These shifts co-occur with major global boundary condition changes, namely the intensification of Northern Hemisphere glaciation in the late Pliocene and the mid-Pleistocene Transition.These changes to the global climate system have been shown to alter the significance and presence of eccentricity, obliquity, and precession in paleoclimate records 19,29,59 .However, the 400-kyr cycle of the Wet-Dry Index departs from this pattern, with a relative decrease in strength during 1.7-1.0Ma (Fig. 3d).
Amplitudes: climate variability across orbital timescales.Previous work utilizing the sapropel deposition in the Mediterranean region has solidified the importance of precession on African climate 38 .Studies using the thickness of the dark sapropel layers have further advanced our understanding of the significant role that obliquity and eccentricity modulation may have on the EAM and WAM systems 60 .Analysis of geochemical proxies from the long, well-dated ODP 966/967 Eastern Mediterranean sediment cores have enabled us to resolve both high-and low-frequency climate fluctuations 55 .
We applied stationary and evolutive spectral analysis to the Wet-Dry Index 16 to derive properties that indicate significant influences of precession, obliquity, and eccentricity modulation (Fig. 3c, d).While precession (19 kyr and 23 kyr), obliquity (41 kyr), and high-frequency eccentricity (95 kyr and 125 kyr) are significant in the Lomb-Scargle (Fig. 3c), both stationary and evolutive spectral analyses show that the eccentricity 400-kyr cycle is dominant over the last 3 Myr (Fig. 3c, d).This indicates an important role for the low frequency orbital cycles, supported by a recent compilation of African paleoclimate records 15 .However, a longer record is necessary to decipher even longer periodicities.Our new 4.5 Myr-long δD precip record does not directly resolve precession or 41-kyr obliquity cycles, but we can examine how longer frequency orbital cycles might modulate monsoon strength.The strong 400-kyr eccentricity power in both the δD precip and Wet-Dry Index records (Fig. 3a, c) highlights their similar forcing mechanism, despite the relatively weak spectral power found at this band in solar radiation (Fig. 3e).In Fig. 6 Potential 2.3-Myr orbital cycle contribution to drying trend.The residual δD precip variability from ODP 966/967 after removal of 400-kyr and 1-Myr periodicities (a) and 400-kyr, 1-Myr, and 2.3-Myr periodicities (b).This residual record exhibits a drying trend (red) when only the relatively higher frequencies are filtered out, but no trend when the 2.3-Myr cycle is also filtered out.The 2.3-Myr cycle contributes 25% of the original δD precip reconstruction, bringing orbital influences up to 50% of the signal.conjunction with the shorter, higher-resolution records, our longer, lower-resolution δD precip record highlights the importance of the low frequency cycles, suggesting that a wide range of frequencies are critical, and potentially related, in tropical precipitation.
Rose et al. 14 presented high-resolution leaf wax isotopes from two 100-kyr-long intervals at 3.05 and 1.75 Ma also from ODP 966/967 (Fig. S2).They found ample precession-scale variability in each interval, with modulation of the variability driven strongly by eccentricity.Due to our sampling of sapropel intervals only, our new δD wax data corresponds with the most depleted values in these high-resolution intervals (Fig. S2), often more negative compared with the precession-resolved data likely due to the targeted sapropel sampling and increased organic matter production during the wettest intervals.This illustrates how our sampling approach captures modulation of the D-depleted wet monsoon intervals, rather than the magnitude of individual precession-band monsoon cycles.The comparison reinforces how the amplitude of precession-scale climate variability can be assessed with our low-resolution sampling of the paleoclimate record.
The direct effects of high-frequency precession and obliquity seasonal insolation forcing on northern African monsoon strength have been well documented.Our record shows that low frequency orbital modulation of these cycles also affects monsoon strength, and that the scaling between low-frequency insolation forcing and monsoon strength has remained largely stationary over the past 4.5 Myr.The 400-kyr cycle likely reflects eccentricity modulation of maxima in precession-band insolation cycles and its effect on monsoon strength.In contrast, the 1-Myr cycles are potentially related to the 1.2-Myr modulation of obliquity and its effect on the cross-equatorial insolation gradient and moisture supply to the northern African monsoon system.This long and astrochronologically dated record enables us to examine specific intervals with high-and low-amplitude climate variability.These high-amplitude precession-band packets of variability were particularly strong at 3.8, 1.8, 1.0, and 0.2 Ma according to the δD precip reconstruction (Fig. 2c), times that have been highlighted in other studies as exhibiting evidence of deep lakes 56 (Fig. 4) and transitions in African climate variability 29 .These packets, which are likely to be common across northern and eastern Africa, have been linked to pulses in mammalian origination and extinction 61,62 , although re-examination show only a single possible origination pulse at ~2.00-1.75Ma 63 .These specific times are when both the 400-kyr and 1-Myr cycles contribute to more D-depleted values (Fig. 4), implying that the superposition of orbital drivers leads to major regional climate changes.
Trends: secular climate change over the Plio-Pleistocene.Sapropel records from the eastern Mediterranean clearly document orbitally driven oscillations, but long-term patterns have been historically difficult to quantify and disentangle from cyclicity with these methods.Trends are potentially more discernable using more quantitative measurements of hydroclimate, like from XRF analyses and associated ratios and indices 16 .However, other contributing factors, such as wind direction and basin integration area, could compound the long-term hydrological information.Eastern African records of Plio-Pleistocene climate using paleosol oxygen isotopes have supported long-term drying trends or shifts 25 , but there is still ample disagreement in terms of the extent, timing, and drivers of aridification 64,65 .Further, low-frequency orbital cycles could potentially obscure these long-term signatures if the reconstruction is not of sufficient length.
In order to examine the non-orbital influences on the climate and environment over the Plio-Pleistocene, we remove the dominant 400-kyr and 1-Myr cycles from the leaf wax isotope record of precipitation to characterize residual fluctuations and identify non-orbital climate drivers over the Plio-Pleistocene, that make up 75% of the variance in the original δD precip record (Fig. 5a).The derived residual δD precip demonstrates clear longterm structure over the last 4.5 Myr with two wet-to-dry transitions occurring from ~4-3 Ma and ~2 Ma-present, with at least 20‰ variation (Fig. 5a).When viewed over the whole 4.5 Myr, precipitation strength appears to trend toward more arid conditions by ~10‰ (Fig. 5a), and there is a statistically significant trend (Fig. 6a).
Disentangling the potential drivers of secular change is furthered by comparison of the residual δD precip record with published records of global ice volume 19 , tropical sea surface temperature 66 (SST), and Pacific SST gradients 67 (Fig. 5).All of these records demonstrate a long-term trend in aridification, increasing ice volume and global cooling, tropical SST cooling, and increasing Pacific SST gradients, respectively.These processes are likely related, as global sea level 68 and SST cooling could lead to aridification in eastern Africa 1,69,70 and the zonal SST gradient in the Pacific has been linked with terrestrial precipitation in Africa 23 .For instance, paleoclimate reconstruction studies ascribe the anomalously low SST gradient in the late Pliocene to global warmth and more El Niño-like conditions 67,71,72 .Characterization of Indian Ocean zonal SST gradients throughout the Plio-Pleistocene will improve the interpretations of the residual precipitation mechanisms, but the Pacific gradient throughout this time shares a similar structure, with a shallow gradient during wet northeast African interval and a gradual steepening trend during the drying starting at the beginning of the Pleistocene 67 (Fig. 5d).Our data suggests that just before 3 Ma may have been anomalously dry in the context of a wetter Pliocene (Fig. 5a), which could potentially be driven by the tropical SST gradient shifts, although they show a maximum a few hundred thousand years earlier (Fig. 5d).Terrestrial hydroclimate in the tropical mid-Pliocene has also been shown to be driven by vegetative greening 69 , yet δ 13 C wax (Fig. 2a) shows no trend throughout the duration of our Plio-Pleistocene study interval, suggesting that far-afield processes, like Indian and Atlantic Ocean SST and their effects on moisture supply and convective processes over eastern Africa, drive the long-term aridification in northeast Africa.However, it is also possible that this structure in the residual record is oscillatory, and thus driven by a longer orbital periodicity, such as the 2.3 Myr eccentricity cycle.We additionally filter out the 2.3 Myr cycle to examine this possibility (Fig. 6b).When the 2.3 Myr power of our δD precip is quantified, it contributes nearly 25% of the entire structure of the original precipitation reconstruction.Thus, if this is a real signature, despite not being a robust measure due to the length of the δD precip record, the non-orbital contributions to δD precip variance are reduced from 75% (Fig. 6a) to 50% (Fig. 6b).Regardless, ample (at least 50%) hydroclimate change is derived from nonorbital sources, which likely drive trends in the mean climate, rather than modulation of a precession signal, as our targeted sampling allows us to characterize.
Hominin evolution and environmental change.The strong, consistent 400-kyr band detected in the new δD precip record and previously published Wet-Dry Index 16 is attributed to eccentricity forcing, which modulates precession-band variability and, therefore, contributes to the amplitude of environmental change.Indeed, in comparison with a higher-resolution record of leaf wax isotopes across multiple organic-rich and organic-poor layers 14 , we do find that our sapropel dataset tracks the wet endmembers of a highly-variable, precession-driven reconstruction (Fig. S2).Further, the sapropel intervals provide additional information on long term variability characteristics when compared with the eastern African soil carbonate records (δ 18 O sc and δ 13 C sc ), which typically form in dry seasons and/or hot conditions.The δ 18 O sc and δ 13 C sc records reveal gradual aridification and C 4 expansion, respectively, throughout the Plio-Pleistocene in eastern Africa, at least since 2 Ma 25 .Other offshore δ 13 C wax records indicate long term fluctuations in northeast African environments in the Pliocene 73 .The stable mean of ODP 966/967 δ 13 C wax along with the enrichment trend in δ 13 C sc may indicate that vegetation variability increased through time, which is consistent with wellresolved records of vegetation from eastern Africa 10,74 .While ODP 966/967 δD precip and δ 18 O sc both trend towards more arid conditions, the increase in carbon isotope-derived environmental fluctuation may have had a significant effect on hominin evolution by selecting for more generalist traits and potentially hosting more dramatic shifts between drastically different ecosystems 7 .
The northeast African region is central to the story of hominin evolution, particularly in terms of dispersal routes of the genus Homo.The northeastern Sahara is a hypothesized dispersal route along fluvial or coastal corridors (e.g., the Nile River Valley and the Red Sea), and climate 6 and agent-based modeling 75 has shown that aridity and open vegetation in the region limit these pathways.It is likely that stronger, expanded, and more northward EAM and eastward WAM precipitation would bring respite and resources for hominins in northeast Africa.The expansion of C 3 plants (i.e., trees and shrubs in this region), as seen in 13 C-depleted intervals (Fig. 2a), likely enhanced the connectivity of potential habitats, which has been suggested to be a main factor in dispersal characteristics 76 .We find that there were many periods of more woody Green Sahara intervals, including a particularly woody vegetation (over 50% C 3 ) interval at ~2.2 Ma in the Nile River catchment (Fig. 2a), generally coinciding with the first dispersal out of Africa 77,78 .A large, vegetated area connecting east and northern Africa may have triggered a pull-type response in hominins that were now able to survive using larger cranial capacity 79 and move in the lush, ecologically connected, region of northeastern Africa.Some of the particularly woody intervals ( 13 C-depleted) coincide with relatively arid times, potentially limiting the likelihood of human migration due to lack of water resources and indicating that each Green Sahara interval may not have had ideal conditions for hominin habitat expansion.Similarly, the coincidence of the grassiest ( 13 C-enriched) intervals with moderate or extreme wetness, as documented by δD precip , could indicate that additional WAM moisture into to the leaf wax catchment area of northeast Africa during precession minima led to a more vegetated northeast Sahara relative to the typical desert environment.Thus, in some of the Green Sahara intervals, C 4 grassland expansion may have created linked corridors through which hominins could have dispersed.

Conclusions
Through the generation of long, well-dated, quantitative records of both precipitation and vegetation of northern Africa, we discern long-term trends and a large range of spectral properties in the climate and environment of northeast Africa over the Plio-Pleistocene.We find a mixture of signals from insolation-driven monsoon oscillations and secular variations likely caused by some combination of globally driven ice, temperature, and atmospheric dynamics.Low-frequency orbital cycles, driven by both obliquity and eccentricity, are prominent in northeast African precipitation, despite the minimal direct influence from lowlatitude summer insolation.When these low-frequency cycles are removed from the original precipitation record, we find that at least 50% of variability is caused by non-orbital influences.
Environmental change between C 3 -and C 4 -dominated ecosystems in northeastern Africa would have impacted hominin resource base and dispersal pathways.Connectivity of habitable zones and resource abundance, conditions more likely during eccentricity-modulated humid, C 3 -rich environments, serving as potential migration pathways, could have favored dispersals around and out of Africa at critical time intervals, particularly at 2.2 Ma, when the environment was the most C 3 -plant-dominated in the entire study interval.The northeast African region acted as a connecting pathway out of and into Africa, and ideal climatological and environmental conditions were needed in this area to serve as a reliable habitat for early humans.
Climate and ecosystem fluctuations were large on orbital timescales throughout the Plio-Pleistocene, and likely prior to the period of study, suggesting that there were numerous Green Sahara intervals going back millions of years 80 .These likely occurred nearly every ~21 kyr cycle and would have been particularly strong during high eccentricity intervals on 100 kyr and 400 kyr cycles.Our record also suggests that monsoonal variability is potentially modulated by the 1.2-Myr obliquity cycle and the 2.3-Myr eccentricity cycle.These findings, along with prior documentation of these low-frequency cycles, indicate complex mixtures of orbital forcings across latitudes as well as secular trends and variations.Orbital-scale variability was larger than any long-term trends or shifts, and thus had a larger influence on environmental sensitivity, but the increase in the number of reconstructions that can resolve cycles, amplitudes, and trends of African paleoclimate are central to the robust of hominin evolutionary hypotheses.

Methods
Organic geochemical and isotopic analyses.The sampling approach exclusively targets sapropel layers in order to develop a long record of orbital and secular trends of wet Mediterranean intervals.This strategy allowed analysis of trends in the full 4.5 Myr record while greatly reducing the number of samples that would be required to resolve precession and longer orbital variations.Each sample comprises the entire thickness of a sapropel layer in the core, which ranged from 2 to 22 cm.The whole sapropel layer was sampled to avoid aliasing, and instead to characterize the amplitude of the wet signal above the threshold needed for sapropel generation, deposition, and preservation.Thus, sediment samples may integrate varying amounts of time, but provide information on the average amplitude of the wet end member.However, due to necessary levels of organic content for leaf wax isotope quantification, ghost sapropels and oxidated sections of sapropels were not sampled, and thus bias may have been introduced to favor wet and dry periods, depending on the extent of oxidation 45 .Further, ref. 14 (Fig. S2) discovered very low organic matter concentrations in the nannofossil ooze layers between sapropels, meaning integrated samples are strongly biased to the high organic content at the peak sapropel level.
Organic geochemical analysis was performed on 156 samples and resulted in isotope data from 123 sapropels that cover most of the humid intervals from 4.5 Ma to present.Seven sediment samples did not contain sufficient organic material for robust isotope measurements; 20 samples were analyzed from the same sapropel layer in different sediment cores and were averaged for the final datasets used for interpretations; and six samples were from parts of a single sapropel that was deposited across core section breaks and were weight averaged.New isotope analyses were combined with 24 previously measured data initially presented in ref. 81 .
Lipids from freeze-dried sediment were extracted with dichloromethane:methanol (DCM:MeOH; 9:1 v/v) using a DIO-NEX Accelerated Solvent Extractor 350.The lipids were split into neutral, acid, and polar fractions over an aminopropyl silica gel column with DCM:iso-propanol (2:1 v/v), diethyl ether:acetic acid (24:1 v/v), and methanol as eluents, respectively.The acid fractions were methylated with acidified methanol at 60 °C overnight to produce fatty acid methyl esters (FAMEs).FAMEs were then separated from hydroxy FAMEs with a silica gel column eluted with hexane, DCM (FAME), and MeOH (OH-FAME).A sodium phthalate standard with a previously determined non-exchangeable hydrogen and carbon isotope composition was methylated along with samples to determine the isotope composition of the C and H added during methylation.
From the post-methylation DCM fraction, FAME concentrations were measured using an Agilent GC-MSD/FID (Agilent 7890a GC and 5975 C MSD) equipped with a DB-5ms column (30 m 0.25 mm 0.25 µm).Distributions were summarized using the average Chain length 82 (ACL) and tested for degradation using the Carbon Preference Index 83 (CPI), which, when used for n-alkanoic acids, quantifies the even-over-odd chain length preference to estimate degradation of the waxes, as has been done in previous studies involving leaf wax FAMEs 12 .Both ACL and CPI are calculated using the peak areas of even and odd chain lengths from C 22 to C 32 .
Hydrogen carbon isotopes of waxes (δD wax and δ 13 C wax ) were analyzed on a Thermo Delta V Plus isotope ratio mass spectrometer (IRMS) coupled to an Isolink with a Thermo Trace GC Ultra and are reported on the VSMOW and VPBD scales, respectively (Fig. 2).Samples were measured in triplicate for δD wax and duplicate for δ 13 C wax , and δD wax measurements were corrected for the isotopic composition of the added methyl group during methylation using the phthalic acid methyl ester standard.The H 3 + factor was quantified daily 84,85 .An 'in-house' standard was used to monitor instrument performance.Standard mixtures of n-alkanes with known isotope values (Mix A7 from the Schimmelmann Laboratory at Indiana University) were injected periodically for isotopic calibration.Standard error means calculated using the protocol from ref. 86 were 0.09‰ for δ 13 C wax and 4.30‰ for δD wax .These uncertainties incorporate both the analytical uncertainty on the FAME and Mix A7 measurements, and the uncertainty in realizing the VPDB and VSMOW scales.
Framework for interpreting leaf wax isotopes.Leaf wax biomarkers are a novel paleoclimate proxy that are preserved in a range of archives and amenable to compound specific isotope analyses in order to reconstruct the climate and environmental conditions of terrestrial plants.Plants produce epicuticular waxes to shield leaf surfaces from evaporation and physical damage 87 .These waxes may be ablated and transported by aeolian and fluvial processes to lakes and oceans, where they accumulate and are preserved in sediment over geological time.The waxes include long-chain n-alkanoic acids, which we use to reconstruct water isotope compositions.
Leaf wax δD is primarily controlled by δD precip 88,89 , which, in tropical Africa, is dominantly driven by regional atmospheric dynamics that govern rainfall amount 90,91 .A variety of observational 91,92 , modeling 93 , and paleoclimate 4,12,14,58,70,94 studies have revealed δD precip to be highly sensitive to changes in eastern and western African rainfall on orbital timescales.δD precip can be influenced by a variety of other processes such as moisture source and transport, and a variety of convective processes including the location of convective cells 95 .Therefore, we interpret δD precip as a qualitative indicator of rainfall amount, consistent with previous studies in the region 4,11,12,14,17,95 .
Plants utilize several distinct photosynthetic pathways to fix carbon dioxide 96 , and these pathways result in unique carbon isotopic fractionation patterns reflected in δ 13 C of plant tissues including leaf waxes.This allows for the identification of C 3 versus C 4 plant types from sedimentary δ 13 C wax , which is reflected in surface ocean sediments along the African margin 97 .In eastern Africa, C 3 plants are dominantly woody dicots (trees and shrubs), whereas C 4 plants are dominantly grasses and some sedges 98,99 .The relative abundance of trees and grasses is strongly influenced by precipitation, due to a variety of physiological differences related to water stress, although vegetation can also be influenced by rainfall seasonality, growing season temperature, pCO 2 , herbivory, and fire 26,100-104 .δ 13 C wax records from paleolake sediment have documented that eastern African vegetation fluctuates throughout the entire C 3 -C 4 spectrum on orbital timescales 12,17,102,105,106 .
Corrections for changes in the isotopic compositions in the source water and atmospheric CO 2 are applied to the δD wax and δ 13 C wax measurements, respectively.We use a benthic foraminifera δ 18 O stack 20 to estimate past ocean water isotopes to correct the δD wax for ice volume effects on δD (Fig. S3).We normalize the δ 18 O benthic at each insolation cycle to modern, extract the ice volume component using the δ 18 O benthic at the Last Glacial Maximum, and convert it to δD based on the slope and intercept relationship between seawater and the global meteoric water line (Equation S1; Fig. S3).We then apply this anomaly to each study interval to obtain an ice volume-corrected signal of δD wax .Using a 3 Myr-moving average of δ 13 C CO2 from ref. 107 , we correct our δ 13 C wax measurements for changes in the isotopic makeup of atmospheric CO 2 by interpolating the ages to those of our samples.
C 3 and C 4 plants fractionate hydrogen to different degrees during leaf wax synthesis due to differing metabolic pathways and plant physiologies.This causes different apparent fractionations between leaf waxes and precipitation (ε wax-P ), which can affect paleoclimate records based on δD wax if vegetation changes 89 .We calculated a 'landscape fractionation' 108 based upon C 3 -C 4 plant proportions to calculate δD precip values from δD wax .To calculate %C 3 vegetation on the landscape, we develop a linear mixing model using δ 13 C CO2 -corrected 109,110 δ 13 C wax endmember values for C 3 and C 4 plant types 111 , calculated as −35.4‰ for the C 3 n-C 28 acid endmember and −20.5‰ for the C 4 endmember, after ref. 14 .After applying this C 3 -C 4 mixing model to our δ 13 C wax data, we then applied ε wax-P values of −94‰ and −122‰ for C 3 and C 4 vegetation, respectively 89 , to correct for 'vegetation effects' on δD wax and estimate δD precip 112 , which will be used for further time series analyses of hydroclimate change.Time series analyses.We employ a suite of time series analyses to better understand the drivers of climate and environmental change over the last 4.5 Myr in the Northeast African region.We calculated the Pearson correlation coefficient between various climate and orbital parameters in the ODP 966/967 record.Lomb-Scargle spectral analysis of nonuniformly sampled data with the plomb function in the MATLAB Signal Processing Toolbox 113 was performed to elucidate cycles of significant spectral density (Fig. 3).We focus on the band frequencies between 50% of the Nyquist frequency (25% of the mean sampling frequency) and 25% of the total length of the record.Evolutive Lomb-Scargle spectral analysis was also performed on the δD precip record, the Wet-Dry Index 16 , and mean June 21 st insolation at 30°N 36 using a moving (1 data point step) window size of 32 data-points (~566 kyr) for δD precip , 856 data-points (~428 kyr) for the Wet-Dry Index, and 714 data-points (714 kyr) for insolation (Fig. 3).
Remaining time series analytical techniques required linearly interpolated data to achieve evenly sampled records.Filtering exercises were performed on the newly generated δD precip record and mean June 21 st insolation at 30°N 36 using QAnalySeries 114 to isolate (band-pass) periodicities associated with eccentricity (400 kyr and 2.3 Myr) and obliquity (1.2 Myr) (Fig. 4).The specific frequencies and widths used for filtering were targeted using the Lomb-Scargle spectral analysis results of δD precip (2.49743 ± 0.2, 1.06632 ± 0.15, 0.448977 ± 0.14) and insolation (2.475 ± 0.2; 1 ± 0.2; 0.45 ± 0.175).These significant spectra within the predetermined frequency bounds in the ODP 966/967 reconstructions were then subtracted from the original δD precip record to yield a δD precip residual dataset, which had one outlier removed, that captures the non-orbital evolution of δD (Fig. 5a).Further, the 2.3 Myr cycle was also filtered out of δD precip to compare the relative influence of secular trends and low frequency cycles (Fig. 6).By dividing the variance of each filtered band by the variance of the original δD precip record, we then determined the percent contributions of the orbitally driven cycles and the secular residual variability to the measured paleoclimate record.

Fig. 1
Fig. 1 Vegetation map of study location.Map of study cores from ODP Leg 160 Sites 966 and 967 in the Eastern Mediterranean with continental Northeast Africa and the fluvial component of the leaf wax shed, i.e., the Nile River watershed, highlighted in blue, with biomes 46 , lakes, and rivers.Mapping software, biomes, and watershed shapefiles from opensource Quantum Geographic Information Systems (QGIS): http://qgis.osgeo.org.

Fig. 3
Fig.3Spectral analyses from leaf wax data, Wet-Dry Index, and insolation.Stationary and evolutive Lomb-Scargle spectral analyses of the δD precip record from this study (a, b), the ref.16 Wet-Dry Index (c, d), and mean insolation36 on June 21 st from 30°N (e, f).Confidence intervals (CI) of Lomb-Scargle spectral power denoted (left) in dashed gray lines.
There is no systematic trend in ACL or CPI, or other indicators such as largest compound or relative abundance of long:middle (C 26 + C 28 + C 30 + C 32 / C 22 + C 24 ) or long:short (C 26 + C 28 + C 30 + C 32 / C 16 + C 18 ) chain lengths, throughout the core.Thus, the long-chain n-acid distributions indicate a common leaf wax source and show no indication for the evolution of leaf wax molecular distributions since 4.5 Ma.