Multi-phase ecological change on Indian subcontinent from the late Miocene to Pleistocene recorded in the Nicobar Fan

Modern grasslands on the Indian subcontinent, North and South America, and East Africa expanded widely during the late Miocene – earliest Pleistocene, likely in response to increasing aridity. Grasses utilizing the C₄ photosynthetic pathway are more tolerant of high temperatures and dry conditions, and because they induce less C isotope fractionation than plants using the C₃ pathway, the expansion of C₄ grasslands can be traced through the δ¹³C of organic matter in soils and terrigenous marine sediments. We present a high-resolution record of the elemental and isotopic composition of bulk organic matter in the Nicobar Fan sediments from IODP Site U1480, off western Sumatra, to elucidate the timing and pace of the C₃–C₄ plant transition within the ∼1.5 × 10⁶ km² catchments of the Ganges/Brahmaputra river system, which continue to supply voluminous Himalaya-derived sediments to the Bay of Bengal. Using a multi-proxy approach to correct for the effects of marine organic matter and account for major sources of uncertainty, we recognize two phases of C₄ expansion starting at ∼7.1 Ma, and at ∼3.5 Ma, with a stepwise transition at ∼2.5 Ma. These intervals appear to coincide with periods of Indian Ocean and East Asian monsoon intensification, as well as the expansion of Northern Hemisphere glaciation starting at ∼2.7 Ma. Our data from the deep sea for a multi-phased C₄ expansion on the Indian subcontinent are in agreement with terrestrial data from the Indian Siwaliks.

The Neogene expansion of C₄ vegetation represented a major reorganization of the Earth’s terrestrial biome. In many parts of the world, after ∼20 million years, the late Miocene (∼11–5 Ma) witnessed the establishment of grasslands dominated by grass species using the C₄ photosynthetic pathway (Edwards et al. 2010; Polissar et al. 2021). Distinct carbon isotopic values of C₃ and C₄ plants are widely used to reconstruct past hydroclimate, where wetter conditions are associated with more abundant C₃ plants, and drier conditions by greater C₄ coverage (e.g. Kirkels et al.2022,a, b). Using the premise that plant-derived material retains the carbon isotopic signature of its photosynthetic pathway during transfer from plant to sediment, we analysed the isotopic composition of organic carbon (OC) preserved in deep-marine sediments from the Nicobar Fan as a proxy for palaeo-environmental change. We use these carbon isotopes to investigate likely long-term changes in the hydroclimate of the sediment source area for the Middle Miocene-Pleistocene Nicobar Fan, Indian Ocean. In particular, we analysed the finer-grained sediments (mudstones) from the siliciclastic Nicobar Fan as fluvial transport of fine particles generally leads to the enrichment of OC and their effective offshore transport as suspended load (e.g. Bianchi et al. 2018).

Shortly after 10 Ma, a major palaeoclimate transition appears to have led to a synchronous expansion of grasslands on the Indian subcontinent as well as in East Africa and North and South America (e.g. Cerling et al. 1997). Evidence for this comes from the C isotopic composition of tooth enamel (Cerling et al. 1997), palaeosol carbonate nodules (Quade et al. 1989; Freeman and Colarusso, 2001; Karp et al. 2018), terrestrial sediments (Vögeli et al. 2017) and terrigenous marine sediments (Freeman & Colarusso, 2001). All of these record a shift towards higher δ¹³C values that suggest an expansion of vascular plants utilizing the more complex C₄ photosynthetic pathway rather than the C₃ pathway of woody plants, marine photosynthesizers and many grasses. Because the C₄ pathway produces a lower degree of C isotopic fractionation, the increase in δ¹³C of both bulk organic matter and terrestrial biomarkers is interpreted to reflect the expansion of C₄ plants, most of which are grasses, and are better suited than C₃ plants to warmer and more arid conditions. The expansion of C₄ plant coverage in the late Miocene is therefore interpreted to be an indication of increasing aridity (e.g. Quade et al. 1989; Cerling et al. 1997; Freeman & Colarusso, 2001; Zhisheng et al. 2001; An et al. 2005; Scheiter et al. 2012; Karp et al. 2018; Khim et al. 2020; Morley, 2018), but it should also be noted that because the C₄ pathway involves a CO₂ preconcentration mechanism, low atmospheric pCO₂/pO₂ will additionally favour C₄ grassland expansion (Cerling et al. 1997). In a study of past vegetation patterns in the Indian Siwaliks (middle and lower Gangetic Plain), Basu et al. (2015) concluded that the early appearance of C₄ plants occurred at ∼11 Ma. Also, from an analysis of the δ¹³C values of bulk soil organic matter, n-alkane and n-alkanoic acid in the palaeosols from Late Mio-Pleistocene Siwaliks, Roy et al. (2020) identified two phases of vegetation change between at ∼7 Ma and at ∼3 Ma – an expansion of C₄ plants during 7–3 Ma and a higher abundance of C₃ plants after 3 Ma. They invoked a morpho-tectonic control for these changes in vegetation.

The late Miocene expansion of C₄ vegetation has been considered as globally synchronous, but site-specific studies reveal differences in the timing of this expansion and suggest that local conditions played an important role (e.g. Dunlea et al. 2020; Tauxe & Feakins, 2020). Data from marine sediments, since ∼6 Ma in the Bay of Bengal at IODP Site U1445, led Dunlea et al. (2020) to propose that C₄ vegetation on the Indian Peninsula existed before the end of the Miocene but expanded to substantially greater abundance during the mid-Pliocene to mid-Pleistocene (∼3.5–1.5 Ma). In this paper, we also identify the ∼3.5–1.5 Ma expansion (see below). Tauxe and Feakins (2020) re-evaluated the constraints for the carbon isotopic shifts recorded from the Indian subcontinent and demonstrated a diachronous transition ranging in age from ∼7.8 Ma in Pakistan to as late as ∼6 Ma in Nepal. The record from IODP Expedition 355 Site U1457, drilled on the Indus Fan, shows that the transition in peninsular India began at ∼7.2 Ma, with similar records from the African margin witnessing an earlier shift to C₄ dominance starting at ∼10 Ma, with those from Australia and South America during the Pliocene. They argued that diachroneity in vegetation changes around the globe does not invalidate pCO₂ as a driver but is consistent with it being one of several drivers of the global C₄ expansion.

The Indian subcontinent has received particular attention regarding the C₃–C₄ transition in part because of the potential relationship between vegetation and monsoon intensity. The strength of monsoon winds is a fundamental control on aridity on the Indian subcontinent where southwesterly winds in summer bring moisture-laden air from the Arabian Sea and northern Indian Ocean inland. Humid air masses deflected upwards by the Himalaya cause extreme precipitation compared to the dry winters in which northeasterly winds provide little moisture. The separate East Asian monsoon system, while not as pronounced, leads to similar intra-annual precipitation patterns in Central and East Asia. Even if intensification of either monsoon system does not lead to an overall increase in annual precipitation, an increase in summer precipitation, particularly with the trade-off of greater winter aridity, will tend to favour C₄ grasses (An et al. 2005), which are better suited to warm and dry conditions due to their ability to limit water loss during stomatal gas exchange. The intensity of monsoonal winds has also been inferred from marine productivity, which strengthens upwelling off the Arabian Peninsula during the summer monsoon (Kroon et al. 1991; Prell & Kutzbach. 1992), and dust flux proxies (Rea et al. 1998; Wan et al. 2007).

Terrestrial records of monsoon intensity tend to reflect local conditions, which means they may not be representative of continental-scale ecological shifts (Freeman & Colarusso, 2001; Vögeli et al. 2017). Much of the evidence for grassland expansion has come from the Siwalik palaeosols in northern Pakistan and India (Fig. 1), and despite efforts to sample throughout the spatial extent of this unit (e.g. Freeman & Colarusso, 2001; Vögeli et al. 2017), extrapolating regional results from a single geologic unit can introduce substantial uncertainty. Other efforts to understand long-term climate evolution on the Indian subcontinent have focused on the elemental and isotopic composition of organic matter in terrigenous sediments of the Indus and Bengal Fans (Freeman & Colarusso, 2001; Krishna et al. 2013; Karp et al. 2018; Khim et al. 2020), which integrate ecological signals over entire river catchments, and can therefore better capture large-scale environmental changes with less interference from spatial heterogeneity. Clift et al. (2008) used weathering records from the South China Sea, Bay of Bengal and Arabian Sea to propose that over the past 23 million years, there has been a dynamic coupling between Neogene climate and both erosion and deformation in the Himalaya, that is, a strong tectonic control.

The Bengal–Nicobar Fan system (Figs. 1, 2), which is uniquely suited to address climate evolution in South Asia because of the largely continuous record of terrigenous sediment flux to the deep sea since 10 Ma and the extremely high sediment accumulation rates (McNeill et al.2017,a, b; Pickering et al.2019, 2020), also observed in the Bengal Fan (Lee et al. 2019). Submarine fan deposits at IODP Site U1480 are largely the result of various types of sediment gravity flows (SGFs) of Himalaya-sourced material transported via the Ganges/Brahmaputra river system (Pickering et al. 2019, 2020; Chen et al. 2020) and therefore incorporate terrigenous organic matter from the combined catchment area. Petrographic data, including geochemical analyses, all support a single and sustained major sediment provenance to the north from the Himalaya and surrounding areas for the entire time interval cored in the Nicobar Fan (McNeill, Dugan, Petronotis, and the Expedition 362 Scientists, 2017b; McNeill et al.2017a; also see section on provenance in Pickering et al.2019). Arc-derived ash content in sediments at IODP Sites U1480-U1481 suggests that the Sunda forearc contributed only a very minor component. Geochemical and Sr-Nd isotope analyses on the bulk silicate fraction of the Nicobar Fan muds/mudstones indicate an eastern Himalayan source dominated by the Greater Himalaya, with a significant Gangdese arc contribution (a major geological structure in the southern Tibetan Lhasa terrane to the north of the Himalayas), all primarily carried by the palaeo-Brahmaputra river. Thus, we are confident that the mudstones sampled for this study are part of this terrigenous sediment flux from the Himalaya, and not from a separate source area to the east, such as the Irrawaddy drainage system.

The Bengal–Nicobar Fan system contains large amounts of woody material. Rapid export and burial of woody material is a highly efficient pathway of atmospheric CO₂ sequestration (Lee et al. 2019). Large-magnitude, low-frequency rapid wood transport avoids degradation that can occur in soils, where extensive OC degradation limits the efficiency of OC burial at very high erosion rates (ibid.).

Palaeogeographic reconstructions of the Bengal–Nicobar Fan System, using the plate-tectonic reconstructions of Hall (2012), with the Bengal Fan morphology from Curray (2014), combined with core data from DSDP/ODP/IODP sites in the Indian Ocean (Fig. 2 modified from Pickering et al. 2019), suggest that the earliest submarine-fan deposits likely were routed along the eastern side of the Indian Ocean, as axial SGFs along the Sunda subduction zone trench that would have been a linear topographic low until it was overfilled to construct the Bengal–Nicobar Fan System. For a comprehensive explanation of SGFs, see Pickering and Hiscott (2016). The latest Eocene and early Oligocene Andaman Flysch, now as accreted and uplifted sedimentary rocks forming part of the Andaman Islands, represents the oldest trench deposits (Pickering et al. 2019, Figure 14). There was a substantially increased delivery of coarser-grained terrigenous sediment to the Bengal–Nicobar Fan System between 13.5 and 8.5 Ma (but beginning at ∼27 Ma), switching to the Nicobar Fan after ∼9.5–9.0 Ma, and then back to the Bengal Fan after ∼2 Ma (Pickering et al. 2020, Figure 14). High elevation in the evolving Tibetan Plateau (>4 km), starting at ∼55–45 Ma, appears to have been established by ∼20 Ma (Ding et al. 2022) and, therefore, this was likely an important sediment source to the Nicobar Fan before the time interval considered in this paper.

Although IODP Site U1480 is ∼2,000 km from the sediment source in the Bay of Bengal, sediment accumulation rates were commonly higher than at IODP sites on the distal Bengal Fan (Fig. 1). Additionally, the IODP Site U1480 record more fully captures the depositional history since ∼10 Ma with fewer hiatuses, intervals of bit advancement without coring and better age constraints than many of the Bengal and Indus Fan sites (IODP sites U1450, U1451, U1455 and U1457; Expedition 354 Scientists, 2016,a, b, c; Khim et al. 2020). Site U1480, therefore, represents a unique opportunity to capture a largely continuous record of ecologic and palaeoclimate conditions during the past ∼10 Ma on the Indian subcontinent.

Our datasets give a unique view of the timing and pace of the C₃–C₄ transition. While previous studies relied primarily on compound-specific isotope analyses of terrestrial plant biomarkers (e.g. Freeman & Colarusso, 2001), we instead used simple analytical techniques to generate a much larger dataset capable of revealing more subtle temporal trends and reduce the reliance on single measurements. Because the incorporation of organic matter from the marine ecosystem can mimic a shift towards terrestrial C₄ plant cover (e.g. Galy et al. 2008), we employed a multi-proxy approach that used mixing models based on δ¹³C, TOC/TN and Br/TOC ratios to correct for contributions from marine organic matter. Uncertainty in end-member composition can dramatically limit the utility of mixing models, so we used Monte Carlo simulations using randomly perturbed end-member compositions to place reasonable quantitative bounds on the contributions of C₄ plants to the terrigenous organic matter through time.

The sample set consists of the residues from shipboard carbonate/TOC analyses and additional dedicated samples that were vacuum sealed and kept frozen. These samples came from the muds/mudstones in the Nicobar Fan sediments (Units I–IIIA in IODP Site U1480 (0–1,310.02 mbsf) and are representative of the siliciclastic sediments in the Nicobar Fan). An analysis of siliciclastic sediments from the Indus Fan (Arabian Sea) has already established that hemipelagic mudstones (rather than sands deposited from various SGFs, such as turbidity currents) provide the most suitable material for evaluating any changes in C₄ grasses (Feakins et al.2020). The carbonate residues were freeze-dried for 12 hours, while additional dedicated samples were dried at 50^oC until consecutive weighing indicated no further mass loss. All samples were ground in an agate mortar and pestle. Approximately 10 mg of sediment was weighed into Ag capsules and decarbonated for TOC and δ¹³C analyses using aqueous SO₂ (∼6–8% by mass). Following dropwise addition of SO₂ solution, the samples were heated to 60^oC until dry, at which point additional SO₂ solution was added. We continued this procedure until a total of ∼1.5 mL had been added, representing at least a 100-fold stoichiometric excess of SO₂/CaCO₃. Pre-combusted glass pipettes and polypropylene 96-well plates were used to minimize any possibility of contamination with organic matter during decarbonation. The Ag capsules containing the decarbonated sediment samples were sealed in Sn capsules to aid in final conversion of organic C to CO₂ in the elemental analyser combustion column. Separate sample aliquots of ∼40 mg were not subjected to decarbonation and were sealed in Sn capsules for TN analyses.

Prior results (Kennedy et al. 2005) suggested that aqueous SO₂ incompletely decarbonates samples for TOC analysis, but the amount of SO₂ solution added provided a less than stoichiometric ratio of SO₂/CaCO₃, so decarbonation should not proceed to completion. Because our samples were likely to contain trace detrital dolomite ((Ca,Mg)CO₃), which is resistant to decarbonation, and would artificially raise δ¹³C_TOC measurements, we tested the decarbonation procedure on a powdered dolomite sample from the Scripps Institution of Oceanography collections. The C content of the untreated dolomite implies an Mg concentration of ∼33%, and decarbonation of ∼1 mg, far in excess of the amount in our samples, resulted in C concentrations below the analytical detection limit of <10 ug C. We are therefore confident that the decarbonation procedure we used quantitatively removes even unrealistically large amounts of (Ca,Mg)CO₃ and is suitable for analysing TOC and δ¹³C_TOC. Decarbonation with aqueous SO₂ has the added benefit of producing Ca/Mg salts that are not hygroscopic, obviating the need for the rigorous drying procedures required after decarbonation by HCl fumigation.

Analyses of TOC, δ¹³C of TOC (δ¹³C_TOC) and TN were conducted at the University of California Santa Cruz Stable Isotope Laboratory. Carbon and nitrogen isotopic and elemental composition were determined by Dumas combustion using a Carlo Erba 1108 elemental analyser coupled to a ThermoFinnigan Delta Plus XP isotope ratio mass spectrometer. Analytical precision of internationally calibrated in-house standards is better than 0.2 ‰ for δ¹³C, which means the δ¹³C error bars in Figs. 3–5 would be smaller than the data points. Sample isotopic values are corrected for size, drift and source stretching effects. Carbon and nitrogen elemental composition is calculated based on standards of known elemental composition. Precision of these known compounds is determined to be better than 1%, again indicating error bars smaller than data points in Fig. 3. All δ¹³C values are reported relative to the Vienna Pee Dee Belemnite (VPDB) C isotope standard.

Major and trace element analyses were conducted on sample aliquots that were washed three times with DI water and centrifuged between washings to ensure seawater-derived bromide would not interfere with sediment Br measurements. The sample aliquots were analysed by Instrumental Neutron Activation Analysis (INAA) at the Oregon State University Radiation Center. Precision is estimated at ∼0.7 ppm at the 2σ level for Br, and we also used this value as the detection limit.

Data were first filtered to address the tendency for samples with low TOC to give low δ¹³C values (Fig. 3). We established a cut-off of 0.125 wt% C and an overall δ¹³C threshold of −28.5 to discard unreliable data while still retaining as many as possible. For Br analyses, we discarded data points for samples with over 400 ppm Cl, which may indicate incomplete removal of seawater Br based on co-variation between Br and Cl for samples with high Cl content. Furthermore, data for which Br concentrations were less than two standard deviations from the detection limit of 0.7 ppm were also discarded. These procedures led to 302 measurements of TOC, δ¹³C_TOC and TN as well as 28 measurements of Br (and other major/trace elements) that we deemed to be reliable (Tables S1; S2). Only these data were used in further analyses.

We used two separate three end-member mixing models to correct the δ¹³C_TOC values for the influence of marine organic matter and place quantitative bounds on the fraction of organic matter from C₄ plants. Both consisted of a system of linear equations of the form | $Af = b$ | where

with | $c_j^i$ | equal to the composition of end-member j with respect to the elemental/isotopic system i, | $f = \;\left[ {\matrix{ {{f_{C4}}} \cr {{f_{C3}}} \cr {{f_M}} \cr } } \right]$ | with the fractional contribution | ${f_j}$ | from end-member j and | $b = \left[ {\matrix{ {d_1^1\;d_2^1 \;\ldots } \cr {d_1^2\;d_2^2 \;\ldots } \cr {\hskip 2pt1\hskip 8pt1 \hskip 5pt\ldots } \cr } } \right]$ | with data | $d_k^i$ | for the measured value of data point k with respect to the elemental/isotopic system i. A total of 10⁵ Monte Carlo simulations were run in which end-member compositions were allowed to vary randomly within distributions centred on literature values (e.g. Galy et al. 2008) and/or extreme measured values and with standard deviations encompassing a range of reasonable values (Table 1). The distributions of | ${{{f_{C4}}} \over {{f_{C4}} + {f_{C3}}}}$ | (the fraction of terrigenous organic matter from C₄ plants) are represented in Fig. 6, and we generated smoothed curves fitting the data based on the distributions produced by bootstrap resampling of the data and the smoothed values at ages corresponding to those of the data.

Bulk TOC δ¹³C data have been shown to correlate well with C isotope ratios of terrestrial plant biomarkers in Bengal Fan sediments (Freeman & Colarusso, 2001; Galy et al. 2008), suggesting that marine sediment δ¹³C_TOC (Fig. 3) presents a reasonable first-order view of landscape evolution and therefore palaeoclimate conditions. However, because the δ¹³C values of marine organic matter are typically intermediate between those of vascular C₃ and C₄ plant material, a shift towards higher δ¹³C values could reflect greater marine organic matter input rather than a true shift towards greater C₄ plant abundance. Attempting to remove the marine organic matter δ¹³C signal using only δ¹³C_TOC results in a mathematically underdetermined system in which the mixing contributions of three end-members (C₃, C₄ and marine organic matter) cannot be uniquely defined. Because correcting for inputs from marine organic matter offers the potential to place at least quantitative bounds on the C₃/C₄ contributions, we applied a multi-proxy approach to reduce the influence of marine organic matter contributions.

The ratio of TOC to total nitrogen (TN) is a widely used metric for the fraction of marine organic matter in marine sediments because while terrestrial organic C preserved in sediment is generally N-poor, largely proteinaceous marine organic matter tends to have much greater TN, leading to high TOC/TN ratios for terrigenous and low TOC/TN ratios for marine organic matter (Fig. 5). However, TOC/TN as an organic matter provenance proxy suffers from several problems, most notably the potential for inorganic N (primarily as NH₄⁺) to adsorb to clay minerals, producing artificially low TOC/TN values (Müller, 1977), and a wide range in particularly the vascular C₃ and C₄ end-member compositions. We relied on three approaches to lessen these effects: (1) a large dataset of δ¹³C and TOC/TN (Fig. 5), which reduces the importance of individual data that might be influenced by NH₄⁺ adsorption, (2) the additional less explored provenance proxy of Br/TOC, and (3) Monte Carlo simulations to reflect the influence of uncertainty in end-member compositions. We chose to use δ¹³C_TOC as the main proxy for C₃–C₄ transition because the ease of the analytical procedure facilitates the gathering of a large dataset, and indeed, the δ¹³C_TOC shows substantial scatter for nearly all time intervals, which suggest that restricted datasets using a more analytically challenging technique such as compound-specific isotope measurements may not fully reflect the range of isotopic variability during a given time interval.

While we used TOC/TN, measured on all samples, as the main proxy for organic matter provenance, we also analysed the trace element composition of a subset of samples (Figs. 4, 6; Table S2) to use Br/TOC as an additional proxy. The Br/TOC ratio in sediments with predominantly marine organic matter has been shown to far exceed that of sediments rich in terrigenous organic matter (Mayer et al. 2007), presumably due to the incorporation of Br-rich secondary metabolites produced by marine organisms (e.g. Gribble, 1998). Therefore, as Mayer et al. (2007) found, we expect a graph of δ¹³C_TOC and Br/TOC to be consistent with mixing between three extreme end-member compositions, and Fig. 5 confirms this, though samples from IODP Site U1480 are not as Br-rich as those Mayer et al. (2007) analysed from continental shelf cores. To evaluate the potential of Br/TOC ratios as an indicator of organic matter provenance in open ocean settings like Site U1480, we estimated the fraction of C₄ plant coverage using δ¹³C and Br/TOC data to compare with results using the δ¹³C and TOC/TN system (Fig. 6). The degree of agreement between these two methods depends heavily on the choice of marine end-member Br/TOC, but using the ∼6–10 range suggested in Mayer et al. (2007) produces results that are remarkably similar to those generated using TOC/TN ratios. While the two methods of correcting for marine organic matter contributions are not entirely independent (both rely on δ¹³C_TOC data), the agreement between the results they produce is encouraging.

We also used Monte Carlo simulations to help account for the uncertainty in C₄ abundance arising from the possible range of end-member composition. Figure 6 shows the results using 10⁵ end-member compositions for δ¹³C_TOC, TOC/TN and Br/TOC taken randomly from distributions centred at our estimated end-member compositions and with generous standard deviations (Table 1). Because of the scatter in the resulting data, we used a bootstrapping technique to produce a range and best estimate (the resulting distribution median) of a smoothed curve of inferred C₄ landcover (Fig. 6). In general, the raw δ¹³C_TOC record and the inferred C₄ coverage after correction for marine organic inputs show similar behaviour, suggesting that the overall influence of marine organic matter on reconstructions of landscape evolution is small, which accords with the correlation between vascular plant biomarker δ¹³C and δ¹³C_TOC that has been established previously (Freeman & Colarusso, 2001; Galy et al. 2008). However, while the δ¹³C_TOC suggests that C₄ expansion (associated with increasing δ¹³C_TOC) may have started ∼9 Ma, inferred C₄ coverage begins increasing just prior to 7 Ma. The fraction of organic matter from marine sources appears to have increased between 9 and 7 Ma, which causes the magnitude of the marine organic matter correction to increase during this interval and the inferred C₄ coverage to remain nearly constant.

Determining a reasonable range for the δ¹³C composition of the marine organic end-member is particularly challenging because the data do not appear to uniformly support a δ¹³C of ∼−18 to −22 ‰ that the sea surface temperatures at this site would imply for phytoplankton (Descolas-Gros & Fontugne, 1990). Samples from IODP Site U1480E cores 1–2 largely post-date Nicobar Fan deposition and would therefore be expected to record the δ¹³C of marine organic matter. Indeed, these samples average −21.8 ‰ (ignoring the −26.2 ‰ value for U1480E-2H-6 139–140 cm as unrealistic for pure marine organic matter). However, sediments that pre-date the Nicobar Fan, below about 1,250 m, and which should also record a relatively pristine marine δ¹³C signal, have an average δ¹³C of −25.3 ‰ with only two samples exceeding −23 ‰. This raises the possibility that the δ¹³C of modern marine organic matter may not be valid as an end-member composition for samples from deeper in the core that underwent early diagenesis. The mechanism by which organic matter may become depleted in ¹³C is not immediately clear, since processes such as methanogenesis would be expected to leave an organic residue highly enriched in ¹³C. However, as Macko et al. (1994) discuss, degradation of organic matter, particularly protein-rich marine material, could lead to ¹³C depletion of the residue due to a diverse set of reactions and associated fractionations available to more functionalized organic compounds. Because of this apparent discrepancy between C isotope compositions of marine organic matter, we tested the effects of using two different δ¹³C ranges on the model outputs (Fig. 6). In one run, there are two models for C₃ and C₄ plant coverage: one using a marine organic matter δ¹³C of −20 to −23 ‰, representing the highest values in sediments that underlie Nicobar Fan material, and a second using a δ¹³C of −18 to −21% to reflect modern value phytoplankton δ¹³C in appropriately warm regions (Descolas-Gros & Fontugne, 1990). The model results using the modern phytoplankton δ¹³C suggest a lower fraction of C₄ coverage at all times during the past and a smaller initial episode of C₄ expansion commencing ∼7 Ma. We are hesitant to ascribe inordinate physical meaning to the apparent decrease in C₄ coverage at ∼5 Ma due to the data sparsity, primarily resulting from low coring recovery, between ∼4 and 5 Ma. Because the main differences in inferred C₄ coverage lie in the earlier part of the record, it is likely that the first model run with δ¹³C of marine organic matter reflecting the δ¹³C of hemipelagic/pelagic sediments underlying the Nicobar Fan is more likely to accurately reflect the true C₄ plant coverage.

Part of the benefit of using terrigenous seafloor sediments such as those of the Nicobar Fan at IODP Site U1480 is that they spatially integrate heterogeneous signals over a wide enough area that localized changes are less likely to alter continental-scale interpretations. However, sediment records from the Bengal Fan may not fully archive terrestrial conditions during this time interval as the Nicobar Fan appears to have received a substantially greater fraction of Himalayan-sourced sediment (Pickering et al.2019, 2020); therefore, IODP Site U1480 may be better suited for understanding C₄ expansion. Additionally, the analytical facility of conducting bulk TOC δ¹³C analyses compared to compound-specific isotope measurements means that the average temporal resolution of our dataset is over ten times greater than that of biomarker δ¹³C records from Freeman and Colarusso (2001) and Karp et al. (2018). The δ¹³C_TOC record from IODP Site U1480 may therefore reveal trends and timing that are not apparent in lower-resolution biomarker datasets. Finally, analyses of terrestrial samples from specific formations like the well-studied Siwalik palaeosols (Fig. 1) may fail to capture the true onset of a transition that is spatially heterogenous on a continental scale. Indeed, the range in values of both δ¹³C_TOC and inferred C₄ coverage during any given time interval in IODP Site U1480 record supports the hypothesis that the entire C₃–C₄ transition from ∼7.1–2 Ma was patchy, and that local terrestrial samples or sparsely sampled compound-specific isotope records may not convey the full spatial variability in ecosystem composition (Vögeli et al. 2017).

Additionally, when compared with many studies using organic matter from sands/sandstones, the C₃–C₄ data presented in this study (from muds/mudstones) record a higher-resolution signal. Prior data from sandy sediments represent a more time-averaged response as these sediments can be stored in ‘temporary sinks’ in river and coastal systems for ∼100 kyr or longer (Gaudemer & Metivier, 1999; Blöthe & Korup, 2013), complicating the correlation of sedimentary archives to climate forcing. Sandy beds tend to be much thicker than muds/mudstones and represent relatively infrequent events compared with muddy sediment delivery to ocean basins. Sand samples depend on the availability of coarser-grained sediment, especially compared with the abundance of finer-grained sediment. Even though plant material will also be stored with sands, to be incorporated into (relatively infrequent) SGFs, plant material is also supplied continuously to ocean basins via higher-frequency processes such as hyperpycnal flows, nepheloid layers and other plumes of suspended finer-grained sediment from river and coastal environments. Thus, understanding any climate change signals is best achieved by analysing muddy sediments from, for example, the distal parts of the Nicobar Fan. It is, therefore, perhaps unsurprising that many studies of sandy sediments tend to show more gradual (and temporally aliased) signals.

Our data suggest that the expansion of C₄ plant coverage in the South Asia was punctuated by two distinct periods of C₄ expansion rather than a consistent gradual shift or single episode as inferred from many terrestrial and marine records (Quade et al. 1989; Prell & Kutzbach. 1992; Freeman & Colarusso, 2001; Gupta & Thomas, 2003). We infer an initial episode of C₄ expansion starting at ∼7.1 Ma, which broadly agrees with the timing inferred from earlier studies (e.g. Bengal Fan δ¹³C in Quade et al. 1989; δ¹³C in Freeman & Colarusso, 2001) as well as a period of monsoon intensification (Rea et al. 1998; Wan et al. 2007) but slightly post-dates the inferred intensification of summer monsoon wind patterns from the Indian subcontinent (e.g. δ¹⁸O from Quade et al. 1989) and Indian Ocean (Kroon et al. 1991; Prell & Kutzbach. 1992; Gupta & Thomas, 2003), in addition to the shift towards C₄ plant dominance in North America, South America and East Africa (Cerling et al. 1997). The second period of C₄ expansion beginning at ∼3.5 Ma and punctuated by a rapid increase at ∼2.5 Ma is inferred to be of greater magnitude and does not appear in records from either terrestrial palaeosols or marine sediments from the Bengal Fan, though the temporal resolution of these records is typically of insufficient resolution during this time interval. Using geochemical and magnetic data obtained from a ferromanganese crust in the western Pacific, together with climate models, Zhong et al. (2022) interpret their data to suggest an equatorward shift of the westerly jet and humidification of Central Asia during the gradual transition from a warm Pliocene climate to the cool Pleistocene, between ∼2.5 and 2.7 Ma. This time interval coincides with our inferred increase in C₄ coverage (Fig. 7). Lu et al. (2020) used a temporally highly resolved (∼500 years) sedimentological record from the Qaidam Basin, NE Tibetan Plateau, to show pronounced glacial–interglacial climate variability during the interval from 2.7 to 2.1 Ma.

The first episode of C₄ expansion ∼7.1 Ma appears to coincide with the timing of a period of monsoon intensification recorded in aeolian dust flux in the North Pacific (Rea et al. 1998) and South China Sea (Wan et al. 2007) and also agrees well with the timing of δ¹³C increase in carbonate nodules from the Siwalik palaeosols. The decrease in median inferred C₄ coverage from ∼5.5–5 Ma also mirrors the decrease in dust flux, though the gap in the IODP Site U1480 record from ∼4.1–5.1 Ma makes us cautious in this interpretation. This episode of C₄ expansion in IODP Site U1480 record appears to pre-date the increase in δ¹³C of vascular plant biomarkers that is apparent around 6 Ma (Freeman & Colarusso, 2001; Karp et al. 2018; Fig. 6), although there are difficulties in dating C₄ expansion using the compound-specific datasets that are discussed in Section 3.2. Interestingly, the initial pulse in dust flux to sites 885/886 at ∼7.7 Ma does not coincide with either a change in soil carbonate δ¹⁸O, an indication of local monsoon intensity, or the first episode of C₄ expansion, although Rea et al. (1998) note that attributing a cause to this peak in dust flux is difficult.

Other records of monsoon intensity inferred from the loess plateaus suggest that aridification progressed northward through time as the East Asian monsoon intensified, increasing the difference in seasonal precipitation (Shen et al.2017, 2018), which is consistent with our results. We detect the C₄ shift beginning ∼7.1 Ma (cf. the 7.4–7.2 Ma expansion of C4 grasslands on the Indian subcontinent as documented from the Indus Fan, Arabian Sea, by Feakins et al. 2020), which lags slightly behind the change in Indian Ocean monsoon intensity inferred from carbonate nodule δ¹⁸O (Freeman & Colarusso, 2001) but is contemporaneous with the increase in soil carbonate δ¹³C of Quade et al. (1989). Additionally, northward C₄ grassland expansion is consistent with the shift towards higher δ¹³C of mammal tooth enamel, which indicates a diet including more C₄ grasses, shortly after 10 Ma in northern Pakistan but later in China (Arppe et al. 2015).

However, we note that two samples (U1480G-42R-1 12–13 cm and U1480F-96X-1 18–19 cm) show a combination of high TOC (above 1 wt% C), high TOC/TN (∼35 and 57, respectively) and δ¹³C values of −19.0 and −21.7‰ at 8.8 and 7.1 Ma (Fig. 5), which strongly suggests localized areas of anomalously high C₄ coverage. This supports the ‘vegetation mosaic’ proposed by Freeman and Colarusso (2001) in which some environments like river floodplains are likely to be particularly sensitive to monsoon intensity as episodic inundation and desiccation would become more dramatic, favouring C₄ grasses that could resist drier winters. Monsoon intensification is likely, they argue, to have a lesser effect on more mountainous areas not susceptible to periodic flooding, and these environments are therefore more likely to retain C₃ plants. The modern landscape (Galy et al. 2008) and analyses of Siwalik palaeosol organic matter from the eastern front of the Himalaya (Vögeli et al. 2017) support spatial heterogeneity in aridity and C₃ and C₄ plant coverage, as do palynological studies (Morley, 2018) and model results based on orographic effects (Boos & Kuang, 2010). These two samples indicating significantly higher C₄ coverage than the bulk of the dataset may then represent localized pockets of C₄ abundance prior to the regional expansion. This would be consistent with the globally synchronous expansion of C₄ plant coverage beginning ∼9 Ma (Kroon et al. 1991; Prell & Kutzbach. 1992; Cerling et al. 1997) and may argue for a more widespread mechanism than intensification of the South Asian monsoon. However, without additional data to corroborate C₄ expansion starting ∼9 Ma, we are cautious not to over-interpret the dataset.

The explanation by Cerling et al. (1997) that decreasing atmospheric pCO₂ in the late Miocene was a fundamental driver of global climate change is supported by later studies (e.g. Pagani et al. 2005; Beerling & Royer, 2011; Rae et al. 2021). Pagani et al. (2005) pointed out that ‘the fall in pCO₂ likely allowed for a critical expansion of ice sheets on Antarctica and promoted conditions that forced the onset of terrestrial C₄ photosynthesis’.

Finally, expansive wildfires have been proposed as either a cause or result of C₄ grassland expansion (Scheiter et al. 2012), but while records of fire frequency have limited temporal resolution, they do not appear to show an increase that is consistent with the first episode of C₄ expansion starting at ∼7.1 Ma (Karp et al. 2018). The complexity and multitude of explanations for increasing aridity in South and Central Asia make us hesitant to definitively attribute the first C₄ expansion to a single mechanism, but a consensus appears to be building around the primacy of monsoon intensification (e.g. An et al. 2005), which is consistent with our observations.

The second main feature of both the δ¹³C data and inferred C₄ coverage is a gradual increase in inferred C₄ coverage beginning at ∼3.5 Ma and punctuated by a sharp increase in δ¹³C_TOC and C₄ fraction at ∼2.5 Ma, effectively coeval with the major expansion of Northern Hemisphere glaciation (e.g. Raymo et al. 1989). Similar changes have been recorded in the abundance of marine primary producers thought to be particularly responsive to upwelling intensity (Kroon et al. 1991; Prell & Kutzbach. 1992) as well as magnetic susceptibility and aeolian dust flux from loess and marine sediment deposits (Fig. 6; Rea et al. 1998). However, other records of C₄ expansion on the Indian subcontinent do not record an additional transition at ∼2.5 Ma (Quade et al. 1989; Cerling et al. 1997; Freeman & Colarusso, 2001; Karp et al. 2018), though several of these records lack the data coverage to identify such a transition (Fig. 6). While the stepwise ∼2.5 Ma transition in the IODP Site U1480 record could conceivably arise from a sudden change in sediment provenance rather than true C₄ expansion, zircon age spectra from IODP Site U1480 are inconsistent with major sediment provenance shifts (McNeill et al. 2017a).

The increase in C₄ coverage beginning around 3.5 Ma coincides with a pronounced increase in North Pacific dust flux (Fig. 7), perhaps suggesting that monsoon intensification may have been a primary factor in the second episode of C₄ expansion as well. Within the overall trend towards greater C₄ coverage is a sharp transition at ∼2.5 Ma that is notable both for its magnitude and abruptness. While we note that quantitative estimates of C₄ coverage are more affected by changes in end-member composition than trends, the highest inferred C₄ fractions just prior to 2.5 Ma are ∼50% compared with samples implying 80% C₄ coverage just after 2.5 Ma. This increase is of similar magnitude to the change throughout the entire rest of the record but occurs over about 0.5 Myr. The sample density is relatively low just prior to this event, but even the smoothed bootstrapping results – which will be much less sensitive to individual data points – imply a sharp transition at ∼2.5 Ma, suggesting it is unlikely to arise as an artefact of sample density. Such a dramatic, essentially stepwise, change demands a similarly abrupt cause and coincides with the major expansion of Northern Hemisphere glaciation (e.g. Zachos et al. 2001), which made South and East Asian climates more arid (Shen et al.2017, 2018). However, it could be argued that the overall increase in inferred C₄ coverage between ∼3.2 and 1.7 Ma (Fig. 7) may have been more related to uplift of the Himalaya and Tibetan Plateau leading to increase aridity in the source area for the sediments. Our estimates of C₄ fraction agree remarkably well with those of An et al. (2005) based on modelling results assuming multiple episodes of Himalayan uplift and records of East Asian summer monsoon intensity, so we cannot discount the possibilities that regional uplift, a global climate shift and/or spatially complex landscape evolution conspired to produce the ecological changes apparent in the IODP Site U1480 record.

Altogether, the elemental and isotopic composition of organic matter at IODP Site U1480 provides an unusually complete insight of major climate shifts on the Indian subcontinent since ∼10 Ma and suggests two distinct episodes of C₄ expansion in the source area for the sediments. The first, beginning ∼7.1 Ma, coincides with a period of inferred monsoon intensification and agrees well with the timing of a shift in soil carbonate δ¹³C but pre-dates the ∼6 Ma shift inferred from compound-specific isotope studies of limited temporal resolution. Isolated samples prior to this hint at the possibility of earlier C₄ expansion in some areas, which may help reconcile the timing with the slightly earlier C₄ expansion, observed on multiple continents (Cerling et al. 1997). The second episode of C₄ expansion began at ∼3.5 Ma and was punctuated by a rapid shift at ∼2.5 Ma, consistent with a climatic alteration contemporaneous with the expansion of Northern Hemisphere glaciation, a major climatic event whose effects have not been identified in records of C₄ plant coverage on the Indian subcontinent. While it seems reasonable to see a significant change in the relative importance of C₃ and C₄ plants (cf. Sarangi et al. 2021) in response to global cooling at ∼2.5 Ma, which coincides with the dramatic increased sediment accumulation rates (SARs) in the Nicobar Fan (Site U1480 in Pickering et al. 2019, Fig. 9), this shift in SARs appears not to observed in the Bengal Fan, probably because at that time sediment routing was predominantly to the east in the Nicobar Fan. The underlying reasons for both phases of C₄ expansion may be complex and different, for example, with the latter involving at least a contributory tectonic component due to the sustained uplift of the Himalaya and Tibetan Plateau. Our data from the deep sea for a multi-phased C₄ expansion on the Indian subcontinent are broadly in agreement with terrestrial data from the Indian Siwaliks.

The wide array and complexity of potential forcing mechanisms make drawing definitive conclusions about climate vs tectonic drivers challenging (cf. Lenard et al. 2020); however, the uniquely detailed record of δ¹³C_TOC from IODP Site U1480 underscores the potential for both local factors (e.g. monsoon intensification) and global (expansion of Northern Hemisphere glaciation) ones in causing spatially heterogenous landscape evolution. Our approach shows the promise of combining voluminous analytically facile data, multiple proxies and uncertainty estimation techniques to more quantitatively describe geochemical processes that can be linked to environmental change.

The supplementary material for this article can be found at https://doi.org/10.1017/S0016756823000481

Participation on IODP Expedition 362 and post-expedition analyses were supported by the United States Science Support Program and the UK Natural Environment Research Council (NERC). Funding support was provided by the United States Science Support Program Post Expedition Award GG009393-01 to BMH. Dyke Andreasen and Colin Carney of the University of Santa Cruz Stable Isotope Laboratory provided insight into sample preparation and conducted analyses of organic matter. Sarah Feakins of the University of Southern California provided helpful discussion regarding sample preparation and data interpretation. We are particularly grateful to Marta Torres for her guidance and advice during IODP Expedition 362 and for enabling the INAA analyses. All members of the IODP Expedition 362 Scientific Party contributed to sample collection and generation of peripheral datasets. We thank Sarah Feakins, Peter Clift and an anonymous reviewer for their insightful comments in the review process that have led to a much improved manuscript.