Pliocene expansion of C₄ vegetation in the Core Monsoon Zone on the Indian Peninsula

The expansion of C₄ vegetation during the Neogene was one of the largest reorganizations of Earth's terrestrial biome. Once thought to be globally synchronous in the late Miocene, site-specific studies have revealed differences in the timing of the expansion and suggest that local conditions play a substantial role. Here, we examine the expansion of C₄ vegetation on the Indian Peninsula since the late Miocene by constructing a ∼6-million-year paleorecord with marine sediment from the Bay of Bengal at Site U1445, drilled during International Ocean Discovery Program Expedition 353. Analyses of element concentrations indicate that the marine sediment originates from the Mahanadi River in the Core Monsoon Zone (CMZ) of the Indian Peninsula. Hydrogen isotopes of the fatty acids of leaf waxes reveal an overall decrease in the CMZ precipitation since the late Miocene. Carbon isotopes of the leaf wax fatty acids suggest C₄ vegetation on the Indian Peninsula existed before the end of the Miocene but expanded to even higher abundances during the mid-Pliocene to mid-Pleistocene (∼3.5 to 1.5 million years ago). Similar to the CMZ on the Indian Peninsula, a Pliocene expansion or re-expansion has previously been observed in northwest Australia and in East Africa, suggesting that these tropical ecosystems surrounding the Indian Ocean remained highly sensitive to changes in hydroclimate after the initial spread of C₄ plants in late Miocene.

1 Introduction

The expansion of plants using the C₄ photosynthetic pathway is one of the most dramatic reorganizations of the global biome during the Neogene. A widespread late-Miocene expansion (8 to 6 Ma) is well documented, and many studies have interpreted the broadly synchronous timing as ecosystems adapting to decreasing pCO₂ (e.g., Ehleringer et al., 1991; Ehleringer and Cerling, 1995; Cerling et al., 1993, 1997; Herbert et al., 2016). However, an increasing number of studies have shown that the timing, regional patterns, rate, and drivers of C₄ grassland expansion were much more diverse and complex (An et al., 2005; Behrensmeyer et al., 2007; Huang et al., 2007; Edwards et al., 2010; Zhou et al., 2014). Along with low pCO₂, a C₄ photosynthetic pathway is better adapted to higher temperature, aridity, seasonality, and during disturbances such as flood, droughts, and fires (e.g., Edwards et al., 2010, and references therein). The interplay of these parameters varies amongst regions. Resolving the precise timing and factors leading to major changes in vegetation demands site-specific studies (Strömberg, 2011; Zhou et al., 2014).

Our study provides a novel piece of the puzzle in unraveling the complexities of C₄ expansion by constructing a 6-million-year (Myr) record of C₄ vegetation and aridity on the Indian Peninsula. The marine sediment record is from International Ocean Discovery Program (IODP) Site U1445 (17^∘44.72^′ N, 84^∘47.25^′ E; 2503 m water depth; Fig. 1a) drilled in the Bay of Bengal (BoB) close to the mouths of the Mahanadi River. Lithologies at Site U1445 include calcareous fossils, biosilica, silt, and clays (including glauconite), and are generally described as hemipelagic sediment (Clemens et al., 2016). The Indian monsoon dictates climate patterns in the Mahanadi River drainage basin: rainy summers, dry winters, and an annual reversal of wind direction (Gadgil, 2003; Sarkar et al., 2015). Highly sensitive to the seasonal changes, more than 80 % of runoff from the Mahanadi River into the BoB occurs during the summer (Chakrapani and Subramanian, 1990).

Figure 1 (a) Location of IODP Site U1445 in the Bay of Bengal (red star). Site 231 in the Gulf of Aden and Site 722 in the Arabian Sea are plotted for reference (red dots). The topography and bathymetry are represented in the background map. The Mahanadi River and its main tributaries are traced in dark blue. The region shown in the box outlines the zoomed-in area in panels (b) and (c), which are modified from Ponton et al. (2012). (b) Average annual rainfall (meters per year) and (c) natural ecosystems in the region including the Mahanadi River drainage basin.

Previous reconstructions of Neogene C₄ expansion in regions affected by the Indian monsoon use deposits originating from the Himalayas and their piedmont regions (France-Lanord and Derry, 1994; Quade and Cerling, 1995; Quade et al., 1995; Cerling et al., 1997; Freeman and Colarusso, 2001; Sanyal et al., 2004; Behrensmeyer et al., 2007; Galy et al., 2010; Ghosh et al., 2017). The Mahanadi River drains a relatively low-elevation region of the Indian Peninsula distinct from the nearby mountain ranges (e.g., the Western Ghats, the Himalayas, the Indo-Burman ranges; Fig. 1; Xie et al., 2006). With minimal orographic precipitation in the Mahanadi River basin, rainfall in this “Core Monsoon Zone” (CMZ) represents the mean behavior of the Indian monsoon (Fig. 1; Ponton et al., 2012; Sarkar et al., 2015; Giosan et al., 2017, and references therein).

Although agriculture dominates present-day vegetation, models of seasonal climate predict that the natural flora of the Mahanadi Basin would be closed-canopy, moist deciduous forests and moist-to-dry woodlands with rare open spaces (Fig. 1c; Zorzi et al., 2015, and references therein). Today, the region encompasses a range of C₃ and C₄ vegetation, but proxies and models suggest that the plant communities are highly sensitive to glacial–interglacial changes with nearly all flora utilizing a C₄ pathway during the last glacial maximum (Galy et al., 2008; Phillips et al., 2014; Zorzi et al., 2015, and references therein). The behavior of vegetation in the CMZ over million-year timescales is unknown.

Here, we use inorganic bulk geochemical analyses to fingerprint the origin of sediment at Site U1445 as being from the Mahanadi River. Then, we use bulk organic and compound-specific biomarkers at the same site, including carbon and hydrogen isotope measurements of leaf wax fatty acids, to reconstruct the changes in C₄ vegetation and rainfall in the CMZ of the Indian Peninsula over the last ∼6 Myr (Fig. 2).

2 Methods

Over million-year timescales, Site U1445 had a constant sedimentation rate ($\sim \mathrm{115}\pm \mathrm{15}$ m Myr⁻¹; Clemens et al., 2016). We fit a locally weighted spline to biostratigraphic and magnetostratigraphic age constraints from Hole U1445A (Clemens et al., 2016) using CLAM software in R (Blaauw, 2010) to estimate the age of our samples (Fig. A1). Our samples and the age constraints were from the same hole. Uncertainty on the ages is estimated to be less than ±0.2 Myr. Shipboard scientists observed thin turbiditic sequences (∼2–20 cm thick) throughout Site U1445 and the expansion and dissociation of gas hydrates upon recovery that may muddle a higher-resolution record (Clemens et al., 2016). However, Site U1445 has fewer and smaller turbidite deposits relative to other sites drilled in this region, and the records spanning million-year timescales are likely relatively undisturbed.

To determine sediment provenance, we measured major, trace, and rare earth element concentrations on 30 bulk sediment samples spanning 0 to 6 Myr including light and dark layers of sediment (Sect. 2.1). To reconstruct hydrological and vegetation changes, we analyzed bulk organics (Sect. 2.2) as well as compound-specific biomarkers (Sect. 2.3) from 57 samples. We constructed the sampling plan to characterize the differences in organic content and isotope composition between the lighter and darker layers of sediment. As such, samples for organic and biomarker analysis were collected from Site U1445 in pairs, visually targeting relatively light and dark layers at similar depths to capture the variability range on shorter timescales while characterizing longer trends.

2.1 Inorganic analyses of bulk major, trace, and rare earth element concentrations

The samples that we analyzed for major, trace, and rare earth element concentrations were originally collected for moisture and density (MAD) measurements onboard the JOIDES Resolution during IODP Expedition 353. Each sample was collected with a 2 cm diameter plastic syringe that fits into the top of a 10 cm³ glass vial, allowing for the vial to be completely filled with sediment (Clemens et al., 2016). The samples were dried in a convective oven at $\mathrm{105}{}^{\circ }\mathrm{C}\pm \mathrm{5}{}^{\circ }\mathrm{C}$ for 24 h (Clemens et al., 2016). The remaining sample preparation, digestions, and analyses were conducted at Boston University, and a detailed description of the analytical geochemical procedures are presented in Dunlea et al. (2015). In summary, sediment samples were hand-powdered with an agate mortar and pestle. For major elements, sample powders were digested by flux fusion (Murray et al., 2000) and analyzed by inductively coupled plasma-emission spectrometry (ICP-ES). For analysis of trace and rare earth elements, sample powders were dissolved in a heated acid cocktail (HNO₃, HCl, and HF, with later additions of HNO₃ and H₂O₂ after samples were dried down) under clean-lab conditions and analyzed by inductively couple plasma-mass spectrometry (ICP-MS). Three separate digestions of a matrix-matched in-house sediment standard were analyzed with each batch and determined precision [(standard deviation)∕(average)x100] was ∼2 % of the measured value for each element. The international Standard Reference Material BHVO-2 was analyzed as an unknown with each batch, and the results were consistently found to be accurate within precision for each element.

2.2 Analyses of carbon and nitrogen content and isotopes

Analyses of the abundance of total carbon (TC), total inorganic carbon (TIC), total organic carbon (TOC), nitrogen (N), and the δ¹³C of the TOC component were performed at Woods Hole Oceanographic Institution, and the methods are described in Whiteside et al. (2011). In brief, samples for TOC were weighed into tared silver boats and then acidified to remove carbonates in a closed desiccator for 3 d at 60–65 ^∘C over concentrated HCl. All samples were flash-combusted in a Costech 4010 elemental analyzer coupled via a Finnigan-MAT Conflo II interface to a Thermo DELTA V Plus isotope ratio mass spectrometer. Data were recorded and integrated using the Isodat software package. Post-run calculations were performed for blank corrections, quantifications, and final calibrations.

2.3 Isotope analyses of compound-specific biomarkers

The analyses of compound-specific biomarkers were performed at Brown University (e.g., Daniels et al., 2017). Samples were freeze-dried and lipids were extracted from 3.5 to 4.5 g of sediment using a Dionex 350 accelerated solvent extractor (ASE) with dichloromethane : methanol (9:1 v∕v). The fatty acids in the total lipid extract were separated from the neutral lipids using aminopropyl silica gel chromatography, eluting with a dichloromethane : isopropanol solution followed by ether with 5 % acetic acid.

The fatty acids were methylated to form fatty acid methyl ester (FAME) by dissolving the dried down acid fraction in ∼0.3 mL of toluene and ∼1 mL of 5:95 acetyl chloride : methanol with known isotopic composition. Nitrogen replaced the headspace in the vial before they were capped tightly and heated at 60 ^∘C for 12 h. Once the reaction was complete, the FAMEs were separated from the water by-products formed during the methylation reaction. Samples received ∼1 mL of synthetic saline solution (50 g NaCl L⁻¹ of double-distilled water) and ∼1 mL of hexane, were vigorously shaken, and were then allowed to rest until the hexane separated from the water. The hexane fraction was pipetted into a new vial, avoiding the water. Another ∼1 mL of hexane was added to the sample, shaken, and pipetted into the new vial. To clean the solution and isolate the fatty acids, samples were run through a second silica gel column, eluting with hexane to remove unwanted acids and then dichloromethane to acquire the clean FAME fraction.

The FAME fraction was analyzed on an Agilent 6890 gas chromatograph with a flame ionization detector (GC-FID). Sample blanks were analyzed with every batch. The isotope ratios of the FAME fraction (δD_n−acid and δ¹³C_n−acid) were measured on a Thermo Finnigan Delta^plus XL isotope ratio mass spectrometer with a HP 6890 gas chromatograph and a high-temperature pyrolysis reactor for sample introduction. For δD, three injections of each sample were analyzed, and for δ¹³C, two injections of each sample were analyzed. Between every six injections, an in-house lab standard mixture containing known amounts of various n-acids was analyzed to monitor instrument accuracy and precision. The δD and δ¹³C were corrected for the methyl groups added during derivatization using the following equations:

$$\mathit{\delta }{\mathrm{D}}_{\mathrm{corrected}}=\left[(\mathrm{2}n+\mathrm{2})\cdot \mathit{\delta }{\mathrm{D}}_{\mathrm{measured}}+\mathrm{123.7}\cdot \mathrm{3}\right]/(\mathrm{2}n-\mathrm{1})$$

and

$${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{corrected}}=\left[(n+\mathrm{1})\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{measured}}+\mathrm{36.52}\right]/n,$$

where n is the carbon chain length of the compound, and the methanol added has δD $=-\mathrm{123.7}$ ‰ and δ¹³C $=-\mathrm{36.52}$ ‰. Analytical uncertainty was calculated by [standard deviation∕average] of the injections and is typically less than 3 % for δD and less than 1 % for δ¹³C. The standard deviations are reported in Table S3 in the Supplement. For every instrument run, samples were analyzed in random order.

2.4 Correction of plant physiological effects on δD

We corrected the δD for differences in C₃ versus C₄ plant physiology (Fig. E1; Smith and Freeman, 2006; Chikaraishi and Naraoka, 2007). First, we calculated the fraction of C₃ versus C₄ vegetation with δ¹³C_FA for each sample, estimating that C₃ and C₄ vegetation have a δ¹³C_FA of $-\mathrm{37.7}\pm \mathrm{1.8}$ ‰ and $-\mathrm{21.1}\pm \mathrm{1.4}$ ‰, respectively (Chikaraishi et al., 2004; Ponton et al., 2012). Then, we approximate that C₃ plants have a δD that is 30 ‰ lighter than C₄ plants (Smith and Freeman, 2006; Chikaraishi and Naraoka, 2007; Ghosh et al., 2017) and use the fraction of C₃ versus C₄ to correct the δD of each sample for the differences in plant physiology (Fig. E1). The corrected data show the same overall trend as the uncorrected data, except the values are shifted to be more negative and the overall change in δD is steeper (Fig. E1).

3 Results

To determine the provenance of the aluminosilicate fraction, we examined the proportions of Al, Ti, Sc, Nb, La, and Th concentrations, because other elements (e.g., Fe, K, Mg, Si, Zr, Hf) may be affected by continental weathering, sorting during transport, and post-depositional authigenic processes. The results from 30 samples have almost constant element proportions of the selected elements, indicating that the aluminosilicate fraction of sediment did not significantly vary over the past 6 Myr. The composition of the 30 samples, even amongst the light and dark layers, matches the composition of lithologies that comprise the Mahanadi Basin such as Precambrian granite and gneisses of the Indian Craton and associated sedimentary deposits (Sharma, 2009; Fig. B1; Table S1). Marine sediment deposits in other parts of the Bay of Bengal closer to the Krishna and Godavari rivers or the Ganges and Brahmaputra rivers have a more mafic or highly variable composition that is not observed at Site U1445 (e.g., Tripathy et al., 2014; Fig. B1). As such, we interpret our results as recording terrestrial changes in the CMZ, specifically the Mahanadi drainage basin.

The pairs of samples used for organic analyses are spaced ∼28 m apart (∼260 kyr intervals), and the adjacent light and dark layers within each pair were 0.2 to 4.3 m apart in the sediment core (2 to 46 kyr; Fig. A1). The color difference can be related to the total organic carbon content (wt %; TOC) and total nitrogen content with darker layers having 1.0 to 2.8 times more than adjacent lighter layers (Figs. 2a, C1; Table S2).

Long-chain normal fatty acids of leaf waxes are derived from land plants and are well preserved during transport and burial in marine sediment (Eglinton and Eglinton, 2008). We focus on the C₃₀ chain length to avoid possible contamination from non-terrestrial sources that contribute shorter-chain-length fatty acids (Fig. D1; Table S3). The results of the δ¹³C of C₃₀ fatty acid of leaf waxes (δ¹³C_FA) show a 5 ‰ increase from mid-Pliocene to mid-Pleistocene (∼3.5 to ∼1.5 Ma), after which δ¹³C_FA decreases and becomes more variable from ∼1.5 Ma to the present (Fig. 2b). The hydrogen isotope compositions of the leaf wax fatty acids (δD_FA) increase gradually over the past 6 Myr, but they have a wide range amongst light and dark layers and shorter time intervals (Fig. 2c; Table S3). Comparing the earlier and later time intervals when δ¹³C_FA changes, before the mid-Pliocene (3.5 Ma), δD_FA ranges from −179 ‰ to −147 ‰, and after the mid-Pleistocene (1.5 Ma), the δD_FA increases to between −166 ‰ and −126 ‰ (Fig. 2c).

Figure 2 Analyses of 57 samples from Site U1445 in the Bay of Bengal plotted with age (Ma) and (a) total organic carbon (wt %), (b) carbon isotope values of C₃₀ fatty acids from leaf waxes (δ¹³C_FA, ‰), and (c) hydrogen isotope values of C₃₀ fatty acids from leaf waxes (δD_FA, ‰). The black and white dots are pairs of samples from relatively dark and light layers, respectively, at a similar depth. Triangles are samples not in pairs. The black curves are a nine-point moving average of all samples.

4 Discussion

4.1 C₄ expansion on the Indian Peninsula

The δ¹³C_FA of terrestrial plants is primarily a function of the photosynthetic pathway and isotopic composition of atmospheric CO₂ (e.g., Farquhar et al., 1989). In this study, the 5 ‰ increase in δ¹³C_FA is greater than the reconstructed δ¹³C of atmospheric CO₂ (≤1 ‰; Tipple et al., 2010), suggesting that a correction for ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{{\mathrm{CO}}_{\mathrm{2}}}$ would only slightly adjust our results. Thus, we interpret δ¹³C_FA as reflecting the amount of C₄ relative to C₃ vegetation produced in the CMZ.

From ∼6 Ma until ∼3.5 Ma, approximately 51 % to 81 % (mean 69 % ±9 % SD) of the vegetation in the CMZ utilized a C₄ photosynthetic pathway (Fig. E1). Thus, the environmental threshold for the C₄ photosynthetic pathway had already been crossed before the end of the late Miocene. Later in the mid-Pliocene, the reconstruction shows another distinct expansion reaching 64 % to 92 % (mean 81 % ±7 % SD) of C₄ vegetation in the early Pleistocene (Fig. 2b). The change in vegetation from ∼3.5 to ∼1.5 Ma suggests multiple steps of C₄ expansion in the CMZ, rather than a singular late-Miocene expansion. From ∼1.5 Ma to the present, the average proportion of C₄ vegetation decreased and became more variable (58 % to 92 %, mean 76 % ±10 % SD; Fig. 2b), which may reflect the sensitivity of the region to glacial–interglacial variations observed in shorter records from this region (e.g., Zorzi et al., 2015, and references therein).

4.2 Aridification of the Indian Peninsula

After correcting for the effects of plant physiology (Sect. 2.4), the amount of precipitation and mixing of different air masses can each vary the hydrogen isotopic composition of leaf wax fatty acids (δD_FA; e.g., Eglinton and Eglinton, 2008). The mixing of two air masses with unique δD values was recently observed to drive δD of rainfall in New Delhi, India, but, similar to the amount of precipitation, the relatively depleted δD corresponded to wetter conditions (Hein et al., 2017). Thus, we interpreted the δD_FA in this study as a qualitative proxy for aridity or the relative amount of precipitation. The δD_FA results suggest an overall drying of the CMZ on the Indian Peninsula over the past 6 Myr. The shorter-term scatter in the δD_FA record may reflect higher frequency variations in aridity or rainfall.

4.3 Patterns of C₄ expansion surrounding the Indian Ocean

In this section, we compare our record of C₄ expansion with other compound-specific biomarker records of C₄ expansion at sites in the Indian Ocean or the adjacent land and seas. Multiple proxy records document a late-Miocene C₄ expansion in the Ganges or Brahmaputra river basins such as the Siwalik Group in northern Pakistan or BoB sites receiving outflow sediment (Fig. 3c, France-Lanord and Derry, 1994; Quade and Cerling, 1995; Cerling et al., 1997; Freeman and Colarusso, 2001; Sanyal et al., 2004; Behrensmeyer et al., 2007; Ghosh et al., 2017). Collectively, the reported timing of C₄ expansions in the Himalaya region ranges from 9 to 5 Myr, most commonly 8 to 6 Myr (Behrensmeyer et al., 2007). Rather than a uniform timing, detailed sampling of various deposits around the Siwalik regions shows that C₄ vegetation expansion was staggered amongst nearby sub-environments with different local conditions (Ghosh et al., 2017, and references therein). Another biomarker record documents a late-Miocene expansion in a wide continental region north and west of the Arabian Sea (Site 722; Figs. 1a, 3d; Huang et al., 2007). Once C₄ vegetation expanded at each of these sites, the records suggest there is generally little or no systematic change in the amount of C₄ vegetation after the late Miocene.

Figure 3 (a) Proxy estimates of pCO₂ over the past 6 Myr (green, δ¹¹B, Bartoli et al., 2011; red, alkenone, Pagani et al., 2009; blue, B∕Ca Tripati et al., 2009; yellow, alkenone, Seki et al., 2010; purple, δ¹¹B, Seki et al., 2010). (b) LR04 global stack of benthic δ¹⁸O records as a proxy for global ice volume (Lisiecki and Raymo, 2005) with arrows marking early glaciation events (De Schepper et al., 2014). Blue and red highlight values above and below modern δ¹⁸O, respectively. (c) δ¹³C of C₃₃ alkanes from Siwalik paleosols in northern Pakistan (white dots) and from sediment at Site 717 in the Bengal Fan (black dots; Freeman and Colarusso, 2001). (d) δ¹³C of C₃₁ alkanes at Site 722 in the Arabian Sea (Huang et al., 2007), which integrates vegetation variability from north and east of the Arabian Sea. (e) δ¹³C of C₂₈ fatty acids at Site 231 in the Gulf of Aden, which records vegetation in East Africa (Feakins et al., 2013; Liddy et al., 2016). (f) δ¹³C of C₃₀ fatty acid at Site U1445 in the Bay of Bengal, which records vegetation from the Mahanadi Basin on the Indian Peninsula (this study). (g) δ¹³C of C₃₃ alkanes from northwest Australia (Andrae et al., 2018).

In contrast, the CMZ of the Indian Peninsula and a few other records around the Indian Ocean document an expansion of C₄ vegetation during the Pliocene (Fig. 3). Marine deposits in the Gulf of Aden originate from northeast Africa and record a late-Miocene C₄ expansion, followed by a relapse to predominantly C₃ vegetation ∼4.3 Ma and a re-expansion of C₄ plants in the Pliocene (Site 231; Figs. 1a, 3e; Feakins et al., 2005, 2013; Liddy et al., 2016). The Pliocene re-expansion is consistent with other records from tropical East Africa (e.g., Levin et al., 2004; Cerling et al., 2011). A C₄ expansion in the Pliocene is also observed in northwest Australia (Site 763A; Fig. 3g; Andrae et al., 2018). There is little evidence of significant C₄ vegetation prior to the Pliocene, suggesting a relatively late onset of C₄ vegetation expansion in northwest Australia (Fig. 3). Additionally, there is evidence that East, South, and Central Asia experienced multiple steps of C₄ expansion through the Pliocene (e.g., An et al., 2005; Passey et al., 2009; Zhou et al., 2014; Miao et al., 2017; Koutsodendris et al., 2019). Collectively, a significant regional expansion in the Pliocene, distinctly after the first late-Miocene expansion, is common at least amongst tropical East Africa, northwest Australia, Asia, and the Indian Peninsula.

4.4 Triggers of C₄ expansion in the Pliocene

The adaptations of the C₄ photosynthetic pathway provides a competitive advantage over C₃ vegetation in environmental conditions with low pCO₂, high temperature, high aridity, or extreme seasonality that can lead to frequent floods, droughts, and fires (e.g., Edwards et al., 2010, and references therein). Many studies hypothesize that the global expansion of C₄ vegetation in the late Miocene was triggered by pCO₂ decreasing below a temperature-dependent threshold (e.g. Polissar et al., 2019; Herbert et al., 2016, and references therein). During the Pliocene, decreasing pCO₂ may have also contributed to the expansion of C₄ vegetation (Fig. 3a; Pagani et al., 2009; Tripati et al., 2009; Seki et al., 2010; Bartoli et al., 2011). However, the heterogenous regional response of the Pliocene C₄ expansion suggests pCO₂ cannot be the only driver. Given that the Pliocene C₄ expansion is observed in multiple regions surrounding the Indian Ocean that have highly seasonal rainfall, regional hydrodynamics likely played a role in the Pliocene C₄ expansion.

In the modern era, a complex interplay of multiple modes of variability dictate the unique seasonal and multi-decadal precipitation patterns in the regions surrounding the Indian Ocean. For example, rainfall on the Indian Peninsula follows quintessential monsoon behavior, but it is also tied to the intertropical convergence zone (ITCZ), the Walker circulation, and the Indian Ocean dipole (IOD; Gadgil, 2003; Wang et al., 2017). The biannual rains of the (semi)arid tropical East Africa are related to the ITCZ, monsoon winds, the sea surface temperature (SST), and the Walker circulation (e.g., Williams and Funk, 2011; Tierney et al., 2015; Yang et al., 2015). Monsoon rains annually quench northern Australia, but the El Niño–Southern Oscillation (ENSO), the Walker circulation, and the amount of Indonesian Throughflow (ITF) better explain the precipitation in other parts of Australia (Ummenhofer et al., 2009, 2011a, b). Previous studies reason that changes in monsoon precipitation or the frequency of fires likely triggered local expansion of C₄ vegetation (An et al., 2005; Passey et al., 2009; Zhou et al., 2014; Miao et al., 2017; Andrae et al., 2018). However, the sampling resolution of our million-year study cannot resolve these seasonal, decadal, or centennial rhythms in precipitation, so we examine the underlying processes driving variations in each of these modes of variability.

The physical mechanism of each of these modes of variability can be related back to atmospheric pressure gradients. Monsoons are classically defined as a seasonal wind reversal induced by the pressure gradient caused by differential heating of land and sea (Wang et al., 2017). Broadly, the ITCZ is formed by the convergence and uplift of air masses near the Equator due to differential heating between low and high latitudes and the formation of Hadley cells. The IOD and ENSO are cyclic changes in the zonal gradient across the Indian and Pacific oceans, respectively, and affect the Walker circulation. Variable Walker circulation affects ocean upwelling and thermocline depth, which in turn will influence sea surface temperature. Sea surface temperatures can interact with the atmosphere and reinforce or dampen the atmospheric pressure gradients affiliated with the monsoon, ITCZ, IOD, or ENSO dynamics. Changes in any combination of these atmospheric gradients may have altered the seasonal precipitation in the regions surrounding the Indian Ocean during the Pliocene (Fedorov et al., 2013; Wang et al., 2017).

What could have changed atmospheric pressure gradients in the Pliocene? Indonesian Throughflow was being restricted from 5 to 3 Ma and there is evidence for strengthening SST gradients during the Pliocene (Cane and Molnar, 2001; Wara et al., 2005; Karas et al., 2009; Ford et al., 2012; Zhang et al., 2014; Burls and Fedorov, 2017; Christensen et al., 2017; White and Ravelo, 2020). Changes in the zonal and meridional SST gradients would have affected atmospheric pressure gradients in and surrounding the Indian Ocean and modified seasonal precipitation patterns. At the same time, the onset of Northern Hemisphere glaciation likely also helped reorganize the atmospheric pressure gradients surrounding the Indian Ocean (e.g., Koutsodendris et al., 2019). Before the intensification of glaciation at 2.7 Ma, there is evidence of four early local glaciation events at 4.9–4.8, 4.0, 3.6, and 3.3 Ma (De Schepper et al., 2014). The latter two events are within the age model uncertainty of 3.5 Myr when C₄ vegetation expands in various regions around the Indian Ocean. Thus, it is possible that the early local growth of ice, not the intensification of glaciation, altered the Siberian high or other atmospheric pressure gradients near the Indian Ocean, changing cyclic precipitation patterns that lead to a regional expansion in C₄ vegetation.

On the Indian Peninsula, C₄ vegetation expands gradually from 3.5 Ma until the mid-Pleistocene when the range of variability in C₄ vegetation increases. The timing of the change coincides within uncertainty with the mid-Pleistocene transition (∼1.3 Ma) where the glacial periods transitioned from 41 to 100 kyr cycles. Although the trigger for the switch in cyclicity is unknown, the processes responsible for the change or the ice sheets themselves may have affected atmospheric gradients and the modes of climate variability that altered the dynamics of C₄ expansion on the Indian Peninsula. The sampling resolution from 1.5 to 0 Ma in our study is too low to resolve the change in cyclicity. However, previous studies documenting variations in C₄ vegetation relative to C₃ since the last glacial maximum demonstrate the sensitivity of plant communities to recent glacial–interglacial cycles (Galy et al., 2008; Phillips et al., 2014; Zorzi et al., 2015, and references therein).

5 Conclusion

Our study provides a piece of the puzzle in unraveling the complexities of C₄ vegetation expansion. Although C₄ vegetation was established in the CMZ on the Indian Peninsula before the end of the Miocene, the results of this study show another significant expansion in the Pliocene (∼3.5 to 1.5 Myr). The latter expansion is not observed in many records from the orographically wet Himalayas, emphasizing the spatial heterogeneities in C₄ vegetation response – even within the same monsoon system. However, other regions adjacent to the Indian Ocean, such as tropical East Africa, Asia, and northwest Australia, corroborate the observed expansion in the CMZ of the Indian Peninsula and show C₄ vegetation patterns sensitive to the changes in hydroclimate during the Pliocene. The heterogeneous response suggests that pCO₂ cannot be the exclusive driver of the expansion of the C₄ vegetation in the Pliocene and changes in regional hydrodynamics likely contributed. Restriction of Indonesian Throughflow and the onset of Northern Hemisphere glaciation may have altered the atmospheric pressure gradients and the modes of variability that determine seasonal precipitation patterns in the continents surrounding the Indian Ocean, which may have caused the regional expansion of C₄ vegetation in the Pliocene.

Appendix A

Figure A1 Age–depth model for Site U1445. To determine the ages of our samples, we fit the biostratigraphic and magnetostratigraphic age constraints (Clemens et al., 2016) with an age–depth model using CLAM software in R (Blaauw, 2010). We ran iterations of the model with different types of fit and levels of smoothing, and we identified a locally weighted spline with 0.4 smoothing to best represent the trends observed in the age constraints. The differences between the age model's iterations are not significant and would not change the interpretations of this study.

Appendix B

Figure B1 La–Th–Sc diagram of 30 sediment samples from the Bay of Bengal. Samples from IODP Site U1445 (blue squares) are plotted as well as IODP Site U1444 (green diamonds), National Gas Hydrate Program (NGHP) Site 19 (purple triangles), and NGHP Site 16 (brown circles) in the Bay of Bengal. Average upper continental crust (black square, Rudnick and Gao, 2014), post-Archean average Australian shale (black dot, Taylor and McLennan, 1985), and average mid-ocean ridge basalt (Gale et al., 2013) compositions are plotted for reference.

Appendix C

Figure C1 Carbonate and bulk organic analyses at Site U1445. Analysis of 57 samples at IODP Site U1445 for (a) bulk calcium carbonate content (wt %) calculated as (total inorganic carbon × (8.33313CaCO₃ wt∕C wt)), (b) total organic carbon concentration (wt %), (c) total acidified nitrogen content (wt %), (d) the ratio of total organic carbon to total nitrogen (TOC∕TN, wt % ∕ wt %), and (e) carbon isotopes of the total organic carbon (‰). Black dots represent visually darker layers relative to a lighter layer (white dot) at a similar depth. TOC∕TN shows a distinct increase in the mid-Pliocene, but it remains within the range of TOC∕TN expected for marine organic material.

Appendix D

Figure D1 Long-chain fatty acids from leaf waxes extracted from Site U1445. Plotted from left to right are hydrogen isotopes and then carbon isotopes of leaf wax fatty acids from chain lengths of C₂₆, C₂₈, and C₃₀.

Appendix E

Figure E1 Correcting δD for differences in fractionation due to plant physiology. (a) Raw, uncorrected δ¹³C_FA plotted for comparison. (b) Calculated fraction of C₄ vegetation with gray lines indicating the maximum and minimum boundaries. (c) Uncorrected δD data, plotted here for comparison. (d) The δD data corrected for differences between C₃ and C₄ plant physiologies.