Mid-Late Holocene Sub-Millennial Scale Inverse Trends of South Asian Summer and Winter Monsoons in Sri Lanka

The South Asian Monsoon (SAM) brings precipitation crucial for agriculture across the densely populated region of South Asia. Identifying the key long-term drivers of the SAM is essential to improve the predictability of future monsoonal trends in the context of current global climate scenarios and increasingly frequent drought and flooding events in this part of the world. Here, we reconstruct ∼6000 years of climatic and environmental history of the South Asian summer monsoon-fed Bolgoda South Lake and the Horton Plains, and the winter monsoon-fed Panama lagoon, in Sri Lanka to better understand monsoonal operation over this island and its connection to broader climate systems. Multiple proxies (diagnostic biomarkers, hydrogen and carbon isotopes of individual n-alkane, grain size, and Zr/Rb elemental ratio) indicate a sub-millennial scale decreasing trend of summer monsoon rainfall in the wet zone of Sri Lanka alongside an increasing trend of winter monsoon rainfall in the dry zone during the last ∼6000 years. We also observed multi-centennial scale arid events in the Bolgoda South Lake and Horton Plains records at ∼3,500 and ∼1,000 cal years BP. Inverse monsoonal behavior during the mid- and late Holocene seems to be led by the southward migration of the mean latitudinal position of ITCZ, induced by varying solar energy distribution between the Northern and Southern hemispheres due to Earth’s processional cycle. Our observations are broadly supported by existing paleoclimatic records from the Indian sub-continent, but abrupt arid phases are asynchronous in the regional records. In addition, these short-term arid conditions do not show systematic correlations with the different modes of climate variables known to have teleconnections with the Indian Ocean monsoon.


INTRODUCTION
The South Asian Monsoon (SAM) is a major component of the Asian monsoon system that plays a critical role in global climate and the socio-economic conditions of the nearly one quarter of the world's human population that lives in in South Asia (Gupta et al., 2020). The SAM brings rainfall to the Indian subcontinent and even a slight shift in intensity and duration of the SAM can cause extreme events such as droughts, floods, and mudslides that are occurring increasingly frequently and impacting many regions (Mirza, 2011). The seasonal migration of the equatorial low-pressure trough commonly called Inter Tropical Convergence Zone (ITCZ), determines the timing, intensity, and direction of tropical monsoon winds and the accompanying rainfall Mohtadi et al., 2016). This migration is thought to be mainly controlled by the annual solar cycle causing variations in the solar induced cross-equatorial pressure gradient. However, alterations to monsoonal operation in the past had the potential to impact environments and societies across South Asia, as predicted scenarios for change suggest it will in the future.
The SAM consists of two primary rainfall seasons. In the boreal summer, the ITCZ migrates northward and brings summer monsoon rainfall (Southwest monsoon; SWM) to the Indian sub-continent, including Sri Lanka (Figure 1). Based on the wind and moisture distribution, the SWM can be further divided into the Arabian Sea Branch (ASB) and the Bay of Bengal Branch (BBB), and the latter brings precipitation to Sri Lanka and most of north-eastern India (Dixit and Tandon, 2016). In the boreal winter, the ITCZ retreats southward and reaches its southern-most position around January. This retreat results in a reversal of wind direction and provides the winter monsoon rainfall (Northeast monsoon; NEM) to the southern tip of India and the northern and eastern parts of Sri Lanka (Figure 1; Banerji et al., 2020;Gupta et al., 2020). Modern instrumental records and paleoclimate records suggest teleconnections between the SAM and different modes of climate variables, mainly, El Niño Southern Oscillation (ENSO), North Atlantic Oscillation (NAO) and Indian Ocean Dipole (IOD) (Banerji et al., 2020;Gupta et al., 2020). However, large variations in this general pattern are known, especially for regions more distant to the equator, which makes it difficult to identify key patterns and driving forces of change (Achyuthan et al., 2016;Misra et al., 2019).
Looking at SAM operation in the past, a gradual weakening of the summer monsoon is observed in records outside the Indian sub-continent since the mid Holocene, which is attributed to a southward migration of the mean latitudinal position of ITCZ (Haug, 2001;Fleitmann et al., 2007). This migration is likely caused by the Earth's precession cycle distributing solar insolation differently to the northern and southern hemisphere. However, so far, additional support for this hypothesis over the Indian sub-continent is needed as the ISM records are very heterogeneous and asynchronous (Banerji et al., 2020;Gupta et al., 2020). Furthermore, although some records from the South Indian peninsular are significantly influenced by the Indian Winter Monsoon (IWM), (e.g., Veena et al., 2014;Rajmanickam et al., 2017;Sandeep et al., 2017), data from sites that mainly report IWM variability (e.g., Ranasinghe et al., 2013a;Böll et al., 2014;Misra et al., 2019) are scarce. Due to the large spatial (latitudinal) difference between the available summer and winter monsoon records, it is not easy to compare and understand the relationship between these two major components of the SAM system and their key driving forces.
As an island located close to the equator and in a central position within the seasonal ITCZ migration path, Sri Lanka is influenced by both summer (Southwest monsoon; SWM) and winter (Northeast monsoon; NEM) monsoon rainfall. Therefore, migration of the mean latitudinal position of ITCZ can inversely influence the seasonal lengths of SWM and NEM of Sri Lanka. The central highlands act as an orographic barrier, facilitating the tracking of these seasonal rainfall patterns (SWM and NEM) separately with a very low spatial distribution. These unique geographical and geomorphological settings of the island enable comparative reconstruction of South Asian summer and winter monsoon variability by filling the above knowledge gaps about the operation of Indian Ocean Monsoon. However, there is a relative dearth of long-term climate archives such as natural lakes, peat bogs, and speleothems in Sri Lanka. Hence, researchers have turned to available brackish water coastal lakes and lagoons (e.g., Ranasinghe et al., 2013b;Ratnayake et al., 2017;Gayantha et al., 2020) to reconstruct the past climate. Isolated peat deposits in the Horton Plains in Central Highlands have also been used to infer broad changes in rainfall patterns during the Holocene (Premathilake and Risberg, 2003;Premathilake, 2012). However, the coastal water bodies consist of signals from both terrestrial and marine environments that demand more specific and sensitive environmental proxies for identifying the different organic matter sources and ongoing biogeochemical processes. Moreover, these records have yet to reach their full potential in chronological resolution and use multi-proxy data to infer paleoclimate changes.
To understand the relationship between South Asian summer and winter monsoon variability and orbitally-induced mean position migration of ITCZ, we reconstruct climatic and environmental changes, focusing specifically on rainfall variability at three different sites in Sri Lanka-Bolgoda South Lake, Horton Plains, and Panama Lagoon-that also cover the major climate and altitudinal zones of the island (Figure 1). Our multi-proxy data enable comparison of the operation of the ISM and IWM over the last 6,000 years of the Holocene and, through comparison with available regional records, provide insights into climate forcings that drive the operation of the SAM. Historical records available for the last ∼2500 years in Sri Lanka indicate the ancient kingdoms moved towards the southwest of the island (i.e., dry zone to wet zone) with increased human settlement of the wet zone over time (De Silva, 1981). Here we also try to understand the influence of long-term monsoon variability on past human settlement patterns in Sri Lanka. The resulting relatively well-established chronology for ISM and IWM variability in Sri Lanka not only provides insights into changing environments of relevance to past human societies in this part of the world, but also provides a framework for predicting future changes in the SAM in South Asia under different scenarios of global climate change.

Climate, Topography, and Vegetation of Sri Lanka
Sri Lanka has a tropical monsoonal climate with seasonally varying monsoon systems (Malmgren et al., 2003). The two principal monsoon rainfall seasons are separated by two short inter-monsoon periods that are associated with changing wind direction as part of the annual/seasonal ITCZ migration ( Figure 1). The summer monsoon (Southwest Monsoon/ Yala; SWM) is active in Sri Lanka between May and September and brings >2500 mm of annual rainfall to the southwestern portion of the island. The central highlands of Sri Lanka, which reach a maximum height of ∼2,500 m a.s.l., act as an orographic barrier to these moisture-laden winds, creating a rain shadow effect to the leeward east and northeast and the prevalence of dry conditions outside of Sri Lanka's southwest during this season. The winter monsoon (Northeast Monsoon/ Maha; NEM) brings rainfall between December and February to the northern and eastern portions of Sri Lanka from the Bay of Bengal. By the time the winter monsoon reaches Sri Lanka, having travelled over the Indian sub-continent, these winds bring less moisture than the summer monsoon winds. Nevertheless, the winter monsoon contributes to the majority of the annual rainfall of ∼1750 mm for the northern and eastern portions of the island. The inter-monsoon rainfall occurs when the ITCZ moves over Sri Lanka bringing thunderstorm-type convectional rain and tropical depressions/cyclones (Malmgren et al., 2003;Hapuarachchi and Jayawardena, 2015; also see Supplementary Figure S1).
These climate conditions (amount of rainfall) lead to the division of Sri Lanka into three main climate "zones" (sub regions)-the wet, dry, and intermediate zones ( Figure 1). The central part of the island is mountainous and consists of complex topographical features such as ridges, peaks, plateaus, basins, Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 789291 valleys, and escarpments. The rest of the island is generally flat except for some small hills in the lowlands ( Figure 1). As noted above, these topographical features can strongly influence the different climate elements in Sri Lanka. The three climatic "zones" and topographical features shape the vegetation cover in Sri Lanka. Tropical rainforests and wet zone grasslands characterize the lowland wet zone (below ca. 900 m). The highland wet zone (above ca. 900 m) consists of tropical montane rainforest and montane grasslands (Patana). The lowland dry zone is characterized by seasonal rainforests (mixed-evergreen forests/dry monsoon forests) and dry zone upland and lowland savannas (Talawa and Damana; Perera, 1975;Erdelen, 1988). In addition, significant amounts of mangrove vegetation can also be found along the coastal water bodies, i.e., lagoons and estuaries, in Sri Lanka (Ranawana, 2017).

The Study Sites
The study areas in this investigation were selected to cover the different climatic, topographic, and vegetation zones in Sri Lanka and reconstruct past trends in both the southwest and northeast monsoons (see Figure 1).

Bolgoda South Lake
Bolgoda South Lake (BGSL) is a part of the Bolgoda Lake system located in the coastal lowland wet zone of western Sri Lanka (6°40′ 57'' -6°41′ 42″N; 79°56′ 20″-79°57′ 54″E) (see Supplementary Figure S2). The South Lake is located around 2.8 km from the coast, and it has a surface area of ∼2 km 2 with an average depth of ∼2 m. Between the lake and the sea, hills, and ridges (height up to 20 m) spread along N-S direction forming a morphological barrier that protect BGSL from direct marine influence. It naturally discharges its water into the Indian Ocean via Bolgoda North Lake. BGSL is fed by a small stream from the eastern side, and the lake also directly connects to the Indian Ocean via a narrow artificial stream (recently made) from the western side of the lake. The lake is currently identified as a slightly brackish water body. The BGSL and its catchment receive an annual rainfall of about 2,550 mm from the southwest (summer) monsoon between May and September. Floating and submerged aquatic plants are common in the lake, with a few mangrove swamps/ marshes observed along its southern margin (Ranwella, 1995).

Panama Lagoon
Panama Lagoon (PN) is a shallow estuarine lagoon (highly brackish water lagoon) located on the southeastern coast of the dry zone (6°45′52″-6°46′29″N; 81°48′ 20″-81°49′31″E). It has a surface area of ∼0.73 km 2 and an average depth of 1.48 m. A tributary brings freshwater into the lagoon from the upland region to the west. The salinity of the lagoon varies between 4.5 and 26.6 ppt with monthly variations (Ellepola and Ranawana, 2015). (see Supplementary Figure S2). Panama lagoon is mainly influenced by northeastern (winter) monsoon rainfall between November and January which brings 90% of rainfall into this region with an annual average of 800 mm (Chandrajith et al., 2014). The catchment area of PN is covered by tropical dry evergreen forests (dry monsoon forests). Mangrove vegetation occurs as fringes around a significant part of the lagoon (Ellepola and Ranawana, 2015).

Horton Plains
The Horton Plains (HP) is a national park and a UNESCO World Heritage site located at the eastern extremity of the Central Highlands of Sri Lanka (6°47-6°50′N, 80°46′-80°51′E). The national park covers 31.6 km 2 and is located at an altitude between 2100-2300 m a.s.l. The HP is a gently undulating highland plateau with diverse landscapes, including mires, plains, forested hilltops, grassy slopes, precipices, brooks, and waterfalls (Premathilake and Risberg, 2003). Tributaries of three of the major river systems of Sri Lanka (Mahaweli, Kelani, and Walawe) start in the Horton Plains, forming an essential component of the hydrological regime of the island. The HP is one of the few areas of the island strongly influenced by both southwest monsoon (SWM/summer monsoon) and northeast monsoon (NEM/winter monsoon) rains. Nonetheless, SWM rainfall has the dominant impact on climate and vegetation in HP today. Although the mean annual rainfall in the wet zone is about 2540 mm, for the HP it can exceed 5000 mm. Due to its higher elevation and precipitation, the HP has a markedly colder and wetter climate than the rest of the island, with an unpredictable mountainous microclimate (Chandrajith et al., 2014). The mean annual temperature is 15°C, and January and February are, on average, the driest months of the year. The vegetation in HP mainly consists of Upper Montane Rain Forests (also referred to as Cloud Forests) and Wet Patana grasslands. In addition, a narrow ecotonal belt consisting of shrubs and herbs can be identified (DWC, 2007). (see Supplementary Figure S2).

MATERIALS AND METHODS
Sampling/Sediment Core Collection 70 and 100 cm long sediment cores were taken from Bolgoda South Lake and Panama Lagoon, respectively, using a UWITEC free-fall gravity corer. A 150 cm long peat core was taken from the Horton Plains peat bog using a Russian peat corer. Water quality parameters (i.e., temperature, pH, conductivity, and dissolved oxygen) were measured in situ during the sampling campaign in Sri Lanka. δD and δO values of collected water samples from the sites were measured later at the Stable Isotope Laboratory, MPI-BGC, Jena using a High-Temperature Conversion-Isotope Ratio Mass Spectrometer (HTC-IRMS). The sediment cores retrieved from Bolgoda South Lake (BGSL) and Panama lagoon (PN) were split into two identical longitudinal sections. One section was archived at the German Research Centre for Geosciences (GFZ), Potsdam, Germany, after visually recording the lithological characteristics. The remaining split section was used for XRF core scanning in German Research Centre for Geosciences (GFZ), Potsdam, Germany followed by sub-sampling into 1-cm intervals for geochemical analyses (see below). Only half cores were retrieved with the Russian peat corer, so no archived sections are available for the Horton Plains (HP) core. The HP cores were also sliced into 1-cm intervals for geochemical analyses (XRF core scanning was not viable for this core). The sub-samples were freeze-dried and stored in zip-lock bags in the laboratory at room temperature before pre-treatment and analysis.

C Dating and Construction of Age-Depth Models
Mollusk shells and bulk organic carbon were extracted from sediment samples and analyzed for 14 C ages to establish the chronology and produce age-depth models for different sites.
Considering the diverse organic matter/carbon sources of Panama lagoon, we selected only mollusk shells for 14 C dating the PN core. In the BGSL core, shells or plant macrofossils were unavailable for 14 C dating. In BGSL, terrestrial OM sources are dominant (see below), and there is a minimal direct influence from the ocean. Hence, bulk sediment (bulk organic carbon) was selected for 14 C dating. Ten mollusk shells were selected from the PN core, and 12 bulk sediment samples were selected from the BGSL core for 14 C dating. In the HP core, 15 peat/organic-rich bulk sediment samples were used for 14 C dating after sieving the freeze-dried material with a 1-mm sieve to remove rootlets that could be observed primarily in the top ∼50 cm of the core (also see the Supplementary Tables S1, S2, and S3).
14 C dating was carried out using an Accelerated Mass Spectrometer MICADAS (Ionplus, Dietikon, Switzerland) at the 14 C Laboratory, Max Planck Institute for Biogeochemistry, Jena, Germany (see details about sample preparation in Steinhof et al., 2017). Age-depth models were developed using the R-software package BACON (rbacon) (Blaauw and Christen, 2011). BACON uses Bayesian statistics to reconstruct age-depth models by dividing a core into many vertical sections (of default thickness 5 cm), and through millions of Markov Chain Monte Carlo (MCMC) iterations estimates the accumulation rate (in years/cm) for each of these sections. Using the starting dates for the first section, the accumulation rates for each section are compiled together to form the age-depth model (Blaauw and Christen, 2011).
Due to its direct connection to the sea and high marine influence the 14 C dates of shells in the PN core were calibrated using the Marine20 curve (Reimer et al., 2020) and the post-bomb curve NH Zone 3 was applied for negative 14 C ages (Hua et al., 2013). In addition, 14 C ages in shells were corrected for the local marine reservoir age (ΔR 133 ± 65) based on data from the online Marine Reservoir Database at Queen's University, Belfast, United Kingdom. Bulk organic carbon 14 C dates in the BGSL and HP core were calibrated using the IntCal20 curve (Reimer et al., 2020).

Grain Size Analysis
Grain size analysis was performed on 10 sediment samples from BGSL and 11 samples from PN. 4-6 g of freeze-dried sediment from each layer were used, and the analysis was carried out according to Kilmer and Alexander (1949). First, 0.05% sodium hexametaphosphate was added to the sample to separate the grains. Then, sediment grains greater than 63 µm (sand fraction) were separated by wet sieving using a 63 µm mesh. The sieved extracts, consisting of particles <63 µm (i.e., silt and clay), were homogenized using a Sonic Vibra-Cell VC 750 with an ultrasonic stirrer. The samples were analyzed for grain size using a Sedigraph Micromeritics III Particle Size Analyzer at Linköping University, Sweden. The results were expressed as weight percentages of sand, silt, and clay.

X-Ray Fluorescence Core Scanning
XRF core scanning was conducted only for the BGSL and PN sediment cores as high organic matter and water content in peat, such as that from HP, creates high analytical uncertainties. XRF core scanning was conducted with the ITRAX µ-XRF core scanner at GFZ, Potsdam Germany at 5-mm resolution. Elemental intensities (c.p.s.) of Al, Si, S, Cl, K, Ca, Ti, V, Mn, Fe, Ni, Cu, Rb, Sr, Zr, and relative variations in the coherent and incoherent radiation, were obtained non-destructively. We used the µ-XRF core scanner with a Chromium X-ray source operated at 30kV and 55 mA for 10 S to generate energy-dispersive XRF radiation. The elemental intensities were non-linearly correlated to the elemental concentrations because the direct conversion of the intensity of a single element to its concentration is problematic. This uncertainty in measurement is due to many factors such as matrix effects and varying physical properties of sediments in the core (e.g., water content and grain-size distribution, irregularities of the split core surface, and spatial variations in the thickness of an adhesive pore-water film which forms directly below the protective foil covering the core surface) (Weltje and Tjallingii, 2008). Hence, the results are presented as log ratios of two elemental intensities that are linearly correlated with the log ratios of two elemental concentrations to eliminate the analytical problems with the measurements. (Weltje and Tjallingii, 2008;Martin-Puertas et al., 2017).

Total Organic Carbon Analysis
22 samples from BGSL, 19 samples from PN, and 70 samples from HP were selected for total organic carbon analysis. Around 5-10 mg of freeze-dried sediment was weighed into tin capsules and inorganic carbon was removed by the gradual addition of totally 120 µl H 2 SO 3 (Bisutti et al., 2004). All tin capsules were introduced into an Elemental Analyzer (NA 1110, CE Instruments) and combusted with pure O 2 at 1020°C. Each sample was measured in duplicate at stable isotope laboratory, Max Planck Institute for Biogeochemistry, Jena, Germany (see details about analytical methods in Brooks et al., 2003). An inhouse standard (acetanilide) was measured with the samples for quality control. Relative standard deviation (CV) of the acetanilide standard was 2.9%.

Biomarker Extraction and Analysis
The sediment subsamples were freeze-dried and homogenized before geochemical analysis (lipid extractions). In addition, mollusk shells were removed from the sediment by sieving with a 2-mm sieve or through manual extraction using tweezers, where necessary. 5-15 g of freeze-dried sediment were used for extraction, depending on the TOC contents of samples. Higher amounts of samples were used for lipid extraction in the PN core due to relatively low TOC content in the sediments. Nineteen samples from the PN core, 26 samples from the BGSL core, and 21 samples from the HP core were selected for lipid extraction considering the sediment accumulation rates throughout the cores. Lipid extraction was done using dichloromethane (CH 2 Cl 2 ) and methanol (CH 3 OH) (9:1 v/v mixture) on a Büchi Speedextractor (E-916, Büchi Labortechnik AG, Switzerland). The total lipid extracts (TLE) Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 789291 were then concentrated on a Büchi Syncore under reduced pressure.
The TLE was then separated into four fractions: F1 (alkane/ non-polar hydrocarbon); F2 (aromatic hydrocarbon and ketone); F3 (alcohol and sterol); and F4 (fatty acids and other polar lipids) using the Solid Phase Extraction (SPE) technique on silica gel (0.040-0.063 mm mesh; Merck, Darmstadt, Germany) according to the modified method from Rach et al. (2014). The F1 fraction containing alkanes was eluted with 12 ml of hexane and subsequently de-sulphurised by elution through activated copper powder (with 1 M HCl) in a glass pipet column. The F2 fraction containing aromatic hydrocarbons and ketones was eluted on the same column using 12 ml of hexane: dichloromethane (DCM) (1:1, v/v) mixture. Subsequently, the alcohols and sterols (F3 fraction) were eluted with 12 ml of DCM: acetone (9:1, v/v) mixture, and finally, fatty acids and other polar compounds (F4 fraction) were eluted with 12 ml of DCM: methanol (1:1, v/v) mixture. The F3 fraction was derivatized with bis(trimethylsilyl) trifluoroacetamide (BSTFA) and pyridine, and heated at 75°C for 2 h. All extracts were then dried under a gentle stream of nitrogen before being redissolved in an equal volume (1 ml or 1.5 ml) of the solvent/ solvent mixture.
An aliquot of the extracts (180 µl of F1, F2, and F3 fractions) was spiked with androstane (20 µl) as an internal standard for quantification. The samples were analyzed on an Agilent 6890N gas chromatograph (GC) interfaced with a 5973 MSD quadruple mass spectrometer (MS) with a DB-5 (5% phenyl methyl siloxane) fused silica capillary column (30 m length, 0.25 mm inner diameter, 0.25 μm film thickness) at Linköping University, Sweden. The samples were injected in splitless mode (1 µl inlet pressure of 10 psi with a flow rate 54.3 ml/min), and the injector temperature was maintained at 300°C. A constant flow (1.3 ml/ min) of He was used as the carrier gas. The GC oven was initially maintained at 35°C for 1 min. The temperature was increased first to 130°C at 20°C/min and then to 320°C at 6°C/min where it was maintained isothermally for 15 min. The MS was operated at 70 eV under full scan mode (m/z 40-600), and the ion source and MS quadrupole temperatures were maintained at 230°C and 150°C, respectively. Compounds in the samples were identified based on their retention times relative to an external n-alkane standard mixture (n-C 15 to n-C 33 ) and fragmentation patterns in the NIST MS Library (Version 2.0) and Lipid library (2011). Quantification of the biomarkers was done relative to the peak area of androstane (internal standard).

Compound Specific Isotope Analysis of n-Alkanes
The volume of the n-alkane fraction was reduced to 50 µl to increase the concentration by drying under a gentle stream of nitrogen before compound specific isotope analysis, due to relatively low concentrations of the compounds in samples. Stable hydrogen and carbon isotope (δD and δ 13 C) analysis of individual n-alkanes (n-C 15 to n-C 33 ) was performed using a coupled gas chromatograph isotope ratio mass spectrometer (GC-IRMS system) equipped with a 7890A gas chromatograph (Agilent Technologies, Palo Alto, United States) linked via GC Isolink and ConFlo IV interface to a Delta V Plus Isotope Ratio Mass Spectrometer (Thermo Fisher Scientific, Bremen Germany) at the Max Planck Institute for Biogeochemistry, Jena, Germany. The GC was equipped with a DB1-MS column (60 m length, 0.25 mm inner diameter, 0.25 mm film thickness) for hydrogen isotope measurements and DB-1 MS column (30 m length, 0.25 mm inner diameter, 0.25 mm film thickness) for carbon isotope measurements. The sample was injected in splitless mode at 280°C and 2 µl of the extract was injected. The He carrier gas flow was maintained at 1.3 ml/min. The GC oven was maintained for 1 min at 110°C, before the temperature was increased to 320°C, at 5°C/min and held isothermally for 9 min.
Each sample was measured in triplicate with a standardized mixture of n-alkanes (C 15 to C 33 n-alkane standard mixture) of known isotopic composition measured after every sample (3 GC injections). Only peaks with an amplitude >150 mV were used for evaluation. The carbon and hydrogen isotopic values were converted to the Vienna Pee Dee Belemnite (V-PDB) and Vienna Standard Mean Ocean Water (V-SMOW) scales, respectively, using the values of the above-mentioned n-alkane standard mixture (offset correction). In addition, drift corrections were applied, determined by the standards run after every sample (Werner and Brand, 2001). In addition, the H 3 + factor was determined daily (5.50 ± 0.5; n 25); it remained constant over the measurement period indicating a stable condition in the ion source.

Statistical Analysis
A hierarchal cluster analysis according to Ward's method was performed for n-C 15 to n-C 33 in each site. A Principal Component Analysis (PCA) was also performed for n-alkane distributions at each location to identify the correlation between different nalkanes within the zones identified by cluster analysis. Only nalkanes with a relative abundance >5% in any sample were shown in the PCA biplots of each site to avoid visual overloading of variables. The statistical analyses were performed using the R statistical software package (R Core Team, 2020) and R package "factoextra" for visualization (Kassambara and Mundt, 2020).

Stratigraphy and Chronology
The Bolgoda South Lake sediments mainly consisted of light and dark greyish color sandy-silty clay. Shells were rare in this core except for a single horn gastropod (∼3 cm long) between 39 and 41 cm depths (identified as a gastropod belonging to the family Potamididae) and some tiny mollusk shell fragments in a brownish sand layer deposited between 5.5 and 2 cm depth. This sand layer deposited between 5.5-2 cm showed characteristics of a high energy depositional event that is most likely linked to the 2004 tsunami, which is known to have significantly impacted the region (Wijetunge, 2006). Shell fragments with modern 14 C ages (see Supplementary Table  S2) and the uneven or inclined lower stratigraphical boundary of this sand layer further support this interpretation (Jackson, Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 789291 6 2008; Goff et al., 2012). The uppermost 2 cm layer was characterized by dark grey clayey sediments without shells (see Supplementary Figure S3). For the BGSL core, ages ranged between 5620 cal yr BP and 1410 cal yr BP (model median ages) according to the BACON age-depth model constructed based on 14 C ages in bulk sediments ( Figure 2A). The age-depth model truncated at ca.1410 cal yr BP just below the sand layer (tsunami deposit) at 5 cm depth. This suggests that there is a depositional hiatus present at 5 cm depth that is possibly due to erosion of the top sediment layers during the tsunami. The 14 C age of the top 1 cm layer (present surface layer) at BGSL indicated that an (old carbon) reservoir age of 599 ± 19 14 C years exists for the lake that was subsequently used for reservoir correction for all 14 C dates, assuming the reservoir age is constant throughout the studied period. The marshlands or mangrove swamps located adjacent to the south of the lake likely cause this old carbon reservoir effect (Last and Smol, 2002) (see Supplementary Figure S2).
The Panama Lagoon sediment core consisted of a greyish sandy-silty clay containing many shells throughout the core. From 90 to 37 cm depth, the core was characterized by many intact and semi-intact mollusk shells, including gastropods and bivalves, with horn gastropods being the most abundant. These gastropods were identified as Cerethedia cingulate, a species found today in brackish water mudflats and mangrove habitats. From ca. 10 to 5 cm depth, a light-colored sand layer was present with some shell fragments and an uneven lower boundary which was similar to that found near the top of the Bolgoda South Lake core. This deposit was likewise linked to the 2004 tsunami in the Indian Ocean with an uneven lower boundary and negative 14 C ages of mollusk shell fragments (modern ages). However, no significant erosion can be identified due to the tsunami waves, unlike at BGSL. The uppermost ca. 5 cm layer was characterized by a dark brown layer rich in partially decomposed organic detritus (see Supplementary Figure S3). The age model of Panama Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 789291 8 were considered as recently deposited sediments. In Panama lagoon, the regional marine reservoir age of the shells (133 ± 65 yr BP) was applied for the reservoir correction of 14 C shell dates.
The Horton Plains peat core consisted of dark brown nonlaminated peat deposits. The coring was stopped at 150 cm depth because it was difficult to penetrate an identified sandy layer. Between ca. 50 cm and the core top, the peat deposit included partially decomposed or undecomposed rootlets (see Supplementary Figure S3). The retrieved HP core revealed a depositional history spanning from 5500 cal yr BP (model median age) to the present based on the bulk 14 C dates in the peat sequence ( Figure 2C). However, rapid accumulation rates can be observed at HP, from ca.75 cm to the top of the core (between ca. 600 cal yr BP-present) with significant age uncertainties during this period.

n-Alkanes
In sediment cores retrieved from the three sites along a transect, we identified and measured the concentrations of n-alkanes between n-C 15 and n-C 35 chain-length (see also Supplementary Figure S4). Each sediment core was divided into two zones (clusters) by hierarchical cluster analyses, based on the relative distribution of n-alkanes that is taken to indicate different OM sources for the cores. Dominant n-alkanes in these clusters were identified using PCA biplots for each core (Figure 3).
In BGSL, PC1 accounted for 62.3% variance and was characterized by high positive loadings of short-chain nalkanes from C 16 to C 22 (predominantly even chain n-alkanes) and high negative loadings of long-chain n-alkanes from C 31 to C 35 (odd chain n-alkanes predominant). Short-chain even carbon number n-alkanes such as C 16 , C 18 , and C 20 commonly occur in inter-tidal hypersaline environments such as lagoons and estuaries. They generally occur in autochthonous bacterial (secondary) lipids as a product of decaying algal mats or marine algae in saline environments (e.g., Nishimura and Baker, 1986;Elias et al., 1997;Ekpo et al., 2005;Aloulou et al., 2010;Aghadadashi et al., 2017). The abundance of these n-alkanes can thus be used as a proxy to trace marine influence (e.g., sea water intrusion or marine transgression) in coastal water bodies. Long-chain odd carbon number predominant n-alkanes such as C 27 C 29 , C 31 , C 33 , and C 35 , are mainly derived from the leaf waxes of terrestrial higher plants (Meyers, 2003;Eglinton and Eglinton, 2008). Cluster 1 in BGSL showed positive scores with PC1, whereas cluster 2 showed negative scores ( Figure 3A). This trend implies that cluster 1 in BGSL is characterized by mainly n-C 16 , C 18 , and C 20 alkanes (bacterial lipids) whereas cluster 2 mainly consisted of n-C 31 , C 33 , and C 35 alkanes (leaf wax).
In PN, an inverse pattern of n-alkane distributions can be identified based on hierarchical cluster analyses and PCA. PC1 in PN core accounted for 60.3% of variance characterized by high negative loadings of short-chain n-alkanes from C 15 to C 20 (predominantly even n-alkanes) and high positive loadings of long-chain n-alkanes from C 27 to C 33 (predominantly odd n-alkanes). Cluster 1 in PN showed positive scores with PC1, whereas cluster 2 indicated negative scores ( Figure 3B). This trend implied that cluster 1 in PN was characterized mainly by n-C 27 , C 29 , C 31 , and C 33 alkanes (leaf wax), whereas cluster 2 consisted of n-C 16 , C 18 , and C 20 alkanes (bacterial lipids). PC1 scores changed rapidly from negative to positive scores in the BGSL core at ∼ 3750 cal yr BP while in the PN core, this occurred relatively smoothly at ∼ 3200 cal yr BP (Figures 4A,B).
A different pattern of n-alkane distributions was identified in the HP peat core based on multivariate statistical analysis. According to the PCA biplot of HP, PC1 and PC2 accounted for 33.6 and 23.7% of the total variance, respectively ( Figure 3C). The loadings of the dominant n-alkanes correlated with the PC2 axis rather than the PC1 axis in the PCA biplot. Longer chain nalkanes (namely C 31 and C 33 ) derived from grass (Eglinton and Eglinton, 2008) correlated positively with PC2, whereas midchain n-alkanes C 23 and C 25 originated from mosses (like Sphagnum) (Naafs et al., 2019) correlated negatively with PC2. The present-day environment and abundance of Sphagnum growing near to the HP coring location further supports this interpretation.
In addition, n-C 27 and C 29 , derived mainly from woody plants (Meyers, 2003), showed a moderate negative correlation with PC2. Therefore, PC2 scores at HP appear to generally represent variation between grass (n-C 31 and C 33 alkane) and moss (n-C 23 , C 25 alkane). The majority of cluster 1 depths showed positive scores with PC2, whereas cluster 2 depths showed negative scores with PC2 ( Figure 3C). At HP, PC2 scores were negative with a general core upward trend towards positive scores between ∼5500 and 1700 cal yr BP. From ∼1700 cal yr BP to present PC2 scores were primarily positive albeit with significant fluctuations.

Sterols
The BGSL and PN cores included mangrove-derived triterpenols, primarily taraxerol (e.g., Ranjan et al., 2015;Gayantha et al., 2020). Due to the different total organic carbon (TOC) ranges at BGSL and PN, taraxerol contents were expressed as TOC normalized values to compare the two sites. According to the results, the high average value of taraxerol was observed in the PN core relative to BGSL. In BGSL, taraxerol varied from 717-2578 ng/g dw (dry weight) TOC without any clear trends through the core ( Figure 4A). In the PN core, taraxerol varied from 626-4801 ng/g dw TOC with a clear increasing core upward trend until the tsunami layer ( Figure 4B). In the HP core several sterols and stanols were identified, namely cholesterol, cholestanol, campesterol, campestanol, stigmasterol, sitosterol, and sitostanol. The 5α-sterol to Δ 5 stanol (C 28 , C 29 ) ratio in HP, which is suggested as a proxy for peat humification (Routh et al., 2014;Naafs et al., 2019), showed a general core upward decreasing trend ( Figure 4C).

Compound Specific Carbon and Hydrogen Isotope Variation in n-Alkanes
Among the three sites, core average n-alkane δD values of longchain n-alkanes (n-C 33 , C 31 , and C 29 ) in HP showed the lowest values followed by PN. BGSL showed high (most enriched) Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 789291 FIGURE 4 | Trends of Leaf wax (LW) and bacterial lipid (BL) derived n-alkane (PC1), taraxerol, total organic carbon (TOC), Zr/Rb element ratio, grain size (<63 µm) and carbon and hydrogen isotopes of leaf wax derived n-alkane (n-C 33 ) in (A) Bolgoda South Lake and (B) Panama Lagoon; Trends of moss and grass derived n-alkane (PC2), Stanol: Sterol ratio, TOC, carbon and hydrogen isotopes of moss derived (n-C 23 ) and grass derived (n-C 33 ) alkane in (C) Horton Plains.
Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 789291 10 average δD values in long-chain n-alkanes (Table. 1). Generally, winter monsoon rainfall is more negative in D and 18 O in contrast to summer monsoon rainfall in Sri Lanka and India (Srivastava et al., 2014;Edirisinghe et al., 2017). The primary source of moisture for winter monsoon (or NEM) is the Bay of Bengal where several large rivers discharge their water from the Indian sub-continent. Therefore, the isotope values of water (δD and δ 18 O) in Bay of Bengal are relatively negative (Singh et al., 2010;Achyuthan et al., 2013)

in contrast to the western (open) Indian
Ocean, which is the moisture source for summer monsoon (or SWM) rainfall. This seasonal effect (summer and winter monsoon precipitation due to changing the wind direction) is clearly reflected in the sediment leaf wax n-alkane (C 29 and C 33 ) δD average values in the BGSL and PN cores (Table. 1). In addition, the lowest average leaf wax δD values in HP, relative to both BGSL and PN, reflect the altitude effect on precipitation in the central highlands.
Comparison of the average δ 13 C values of long-chain nalkanes at the three sites indicated no significant differences between the core average values. However, at PN, slightly higher average values were found than at BGSL and HP (Table. 1). Short-chain n-alkanes (n-C 16 and C 18 ) were only present in the coastal sites (BGSL and PN). These sites indicated relatively higher average δ 13 C values in BGSL relative to PN, particularly for n-C 18 alkanes. Average δD values were also relatively high in BGSL short-chain n-alkanes, particularly for n-C 16 alkanes ( Table 1). Mid-chain n-C 23 and C 25 alkane δ 13 C values in HP showed slightly lower values than long-chain n-alkanes in the same core (Table. 1).
Long-chain n-alkane carbon and hydrogen isotopic variations showed a generally inverse relationship between the BGSL and PN cores. In the BGSL core, δD and δ 13 C values for n-C 33 alkanes showed a general core upward increasing trend in their values. Particularly, high δD n-C33 values were observed at ca. 3400 cal yr BP (26 cm depth), deviating from the general trend ( Figure 4A). In contrast to the overall core upward increasing trend at BGSL, δD and δ 13 C values of n-C 33 showed a general core upward decreasing trend in PN ( Figure 4B). In contrast to BGSL and PN, trends for δD n-alkane (i.e., n-C 23 and C 33 ) values in the HP core varied somewhat differently from their corresponding δ 13 C n-alkane values ( Figure 4C). Noticeably higher δD n-C33 and δD n-C23 values occurred at ∼1000 cal yr BP and ∼630 cal yr BP, respectively. In addition, peak values of δ 13 C n-C23 and δ 13 C n-C33 can be observed between 350-160 cal yr BP.
Regarding the δD values of the HP n-alkanes, a core upward increasing trend for n-C 23 was observed with maximum values at ∼ 630 cal yr BP followed by a gradual decline. δD n-C33 also showed a core upward increasing trend, reaching its peak value ∼1000 cal yr BP ( Figure 4C). The δD peak thus appeared slightly earlier in long-chain n-alkanes than in midchain n-alkanes. δ 13 C of n-C 23 alkanes showed relatively low (below average) and stable values between ∼5480-440 cal yr BP followed by a sharp peak between 350 and 160 cal yr BP and then a rapidly decreasing trend towards the top of the core. However, δ 13 C n-C33 showed relatively low values (below −30‰) between ∼5330 and 1250 cal yr BP and shifted to relatively higher values (above -30‰) towards the present ( Figure 4C).

Total Organic Carbon
The highest average total organic carbon (TOC wt%) values were observed in the HP peat deposit, followed by the BGSL and PN cores, respectively. In the HP peat deposit, TOC values varied between 2.13-14.20% with an average of 5.49%. The BGSL core had an average TOC value of 3.11%, varying between 0.71 and 5.09%. In both the BGSL and HP sites the TOC showed a core upward increasing trend ( Figures 4A,C). TOC content was relatively low in the PN core with an average of 0.52 wt%. Except for the relatively high value (1.94%) in the organic-rich surface layer at the top of the core, TOC content was <1% and showed a gentle core upward increasing trend ( Figure 4B).

Grain Size Distribution and Zr/Rb Log-Ratio
Grain size and elemental composition data are only available for Bolgoda South Lake and Panama lagoon core. Both the BGSL and PN cores showed very high contributions from sand (>63 µm) that represented 94 and 93% of the mean values, respectively. However, for the finer fraction (<63 μm, i.e., silt and clay), silt dominated at BGSL, whereas clay dominated in the PN core. In the BGSL core, both silt and clay fractions showed core upward increasing trends that are distinct for silt ( Figure 4A). In contrast, at PN, the clay fraction was relatively high and extended from 1 | Average, minimum and maximum δ 13 C and δ D values of dominant n-alkanes in each core.

n-alkane
Panama Lagoon South Bolgoda Lake n-alkane Horton Plains ∼6000 to 4700 cal yr BP. After that, the clay fraction decreased and showed a stable or slightly decreasing trend until the present. It is, by contrast, difficult to see a clear trend in silt ( Figure 4B). The Zr/Rb ratio of sediment supports the low-resolution grain size data in our records (see discussion). The log ratio of Zr/Rb XRF element intensities showed an overall decreasing trend in BGSL with the lowest values around 3800 cal yr BP. By contrast, in PN, the Zr/Rb log-ratio showed an overall increasing trend with the highest values occurring around 3600 cal yr BP.

Behavior of Monsoon Rainfall and Environmental Impacts in Sri Lanka
Along the transect from the west to the east coast, running through the central highlands of Sri Lanka, the three sampled sites indicate clear paleoenvironmental changes related to monsoon variability since the mid-Holocene. Here, we discuss the multi-proxy data in the Bolgoda South Lake, Panama Lagoon, and Horton Plains to understand the long-term variability of the Indian Ocean monsoon rainfall and its related environmental changes in Sri Lanka during the last ∼6000 cal yr BP.

Bolgoda South Lake
This summer monsoon (or SWM) fed coastal brackish water lake demonstrated a depositional history spanning from ∼5600 to 1400 cal yr BP. Terrestrial plants (leaf wax n-alkanes) dominate between ∼5600 and 3850 cal yr BP in BGSL (cluster 2), whereas bacterial lipids are predominant between 3850 and 1400 cal yr BP (cluster 1) ( Figure 3A). As indicated by BGSL PC1 scores, a rapid shift of the dominant OM source, from terrestrial plants to marine bacterial lipids, can be observed ∼3800 cal yr BP ( Figure 4A). Thus, the marine influence (sea water intrusion) appears to have increased in the BGSL catchment from 3850 to 1400 cal yr BP, likely as a product of decreasing freshwater input into the αlake. Taraxerol is a reliable proxy for the presence of mangrove vegetation surrounding coastal water bodies (Ranjan et al., 2015;Gayantha et al., 2020). However, in BGSL, taraxerol does not show a clear trend throughout the core, with minimal variation from its core average value. This trend indicates that the mangrove vegetation in the vicinity of the site may have been stable over time ( Figure 4A). Furthermore, we infer that salinity changes due to varying marine influence were not significant enough to cause notable changes in the mangrove vegetation around this lake.
Grain size variation primarily reflects energy-driven processes in the lake catchment. High rates of precipitation trigger enhanced erosion and stream discharge that brings coarse particles (sand and coarse silt) to the lake bottom. Conversely, weak precipitation or dry climatic conditions lead to the accumulation of more fine clay particles in sedimentary sequences (Peng et al., 2005;Gayantha et al., 2017). In BGSL, increasing trends of clay and silt content (representing a low energy environment) towards the top of the core reflects the weakening of inland precipitation over time ( Figure 4A). This low-resolution grain size distribution trend is further supported by high-resolution Zr/Rb ratios, which are widely utilized as a proxy for grain size variation (Kylander et al., 2011;Davies et al., 2015). Zirconium is typically associated with the coarse fraction in sediments (medium to coarse silt), whereas Rb is enriched in the fine clay fraction. Therefore, higher values of Zr/Rb indicate coarse sediment and lower values demonstrate finer sediment grains (Kylander et al., 2011). The Zr/Rb ratio in the BGSL core shows an overall core upward decreasing trend, indicating that the finer fraction is increasing in BGSL. The lowest values can be observed ∼3750 cal yr BP indicating minimum freshwater discharge and erosion in the catchment. This observation might be connected to a decline in catchment precipitation and/or arid conditions ( Figure 4A). TOC is mainly associated with the clay fraction in sediments. Therefore, the observed core upward increasing trend of TOC can also be explained by the increasing clay fraction in lake sediments ( Figure 4A).
Changes in stable hydrogen and oxygen isotopes (δD and δ 18 O) in meteoric water are recorded in archives such as continental ice cores (Thompson et al., 1985) and speleothems (Liu et al., 2015;Liu et al., 2020). However, these types of archives are very limited to specific geographical regions, and therefore, measuring the hydrogen isotope ratios (δD values) in biomarkers such as sediment leaf wax n-alkane is a promising method that has been developed recently (Sachse et al., 2012). δD in sedimentary OM can provide palaeohydrological information, i.e. amount, trajectory and source of precipitation (Holtvoeth et al., 2019). Hydrogen isotopes of plant water and plant tissues will, however, be significantly modified due to processes impacting transpiration (i.e., aridity) prior to biosynthesis. Therefore, δD in terrestrial plant lipids reflect a combined signal of precipitation and evapotranspiration (Sessions et al., 1999;Sachse et al., 2012). δ 13 C variability of leaf wax n-alkanes can reflect the water stress of higher plants when there is no major shift in vegetation (C 3 vs C 4 ) and atmospheric CO 2 during the period of concern (Diefendorf and Freimuth, 2017). Both δD and δ 13 C values of leaf wax-derived n-C 33 alkanes show an overall core upward increasing trend at BGSL, indicating a weakening of rainfall in the catchment and increasing aridity from ∼5600 to 1450 cal yr BP ( Figure 4A). Relatively higher δD and δ 13 C values deviated from the general trend between 3700 and 3200 cal yr BP, suggesting the prevalence of drought-like conditions due to a weakening of the summer monsoon ( Figure 4A). A shift of δD n-C33 to lower values ∼3800 cal yr BP indicates the rapid decrease of summer monsoon at this time, coinciding with the rapid shift of OM sources in the sediments (i.e., from leaf wax to bacterial lipids) and the lowest Zr/Rb values ( Figure 4A). Overall, the proxy results suggest a decreasing trend of summer (or SWM) monsoon in the area from 5600 to 1400 cal yr BP with apparent failure of rain ∼3700 cal yr BP and droughts between 3700 and 3200 cal yr BP.

Panama Lagoon
Panama Lagoon is predominantly influenced by winter monsoon (or NEM) rainfall, and the obtained core revealed a depositional history from ∼5900 cal yr BP to the present. All the proxy records showed inverse trends/variations in Panama lagoon compared to Bolgoda South Lake, except the taraxerol contents. Bacterial lipids Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 789291 are dominant between 5900 and 3400 cal yr BP (cluster 2), and leaf wax biomarkers are predominant between 3400 and the present (cluster 1) ( Figure 3B). In addition, PC1 scores in PN gradually increase (negative to positive scores) with time, indicating a gradual decrease of marine influence at the site, in stark contrast to BGSL. This is likely due to increasing freshwater input into the lagoon with time. The taraxerol content showed a gradually increasing trend with time, indicating increasing mangrove vegetation at the site. This suggests that, with additional stream water input, the lagoon water reached optimal salinity levels that enabled mangrove vegetation to flourish around the lagoon. This changeover suggests that Panama lagoon turned from a high salinity to a low salinity brackish water body due to enhanced freshwater input. This changing marine influence is also confirmed by the data from microfauna analysis (see Supplementary Figures S5, S6).
Notably, high clay content can be observed between ∼5900 and 4500 cal yr BP, and this is also confirmed by low (negative) Zr/Rb log-ratios during this period. After that, clay content/fine particles showed no significant variation until the present. The higher clay content between ∼5900 and 4500 cal yr BP suggests a low energy environment, probably due to the weakening of stream discharge. This can occur as a result of significant declines in rainfall and/or a sea level high stand that further diminished the input of stream water into the lagoon. TOC values also show a slight decreasing trend until the organic detritus-rich topmost layer, that may relate to the decreasing clay content in the PN core over time.
In contrast to BGSL, δD and δ 13 C of leaf wax-derived n-C 33 alkanes in Panama lagoon show an overall decreasing trend, indicating a gradual increase in the amount of winter monsoon precipitation and decreasing aridity in the area from ∼5900 cal yr BP to present. At ∼ 4600 cal yr BP, relatively low values of δD n-C33 can be observed deviating from the general trend, suggesting high rainfall around this time. The coarse sediment fraction also increased at PN over this period, further supporting this interpretation. However, the rainfall is not notably low between 5900-4600 cal yr BP according to leaf wax hydrogen isotopes values. Therefore, it is probably a period with a high sea level stand that diminished the energy of freshwater discharge into the lagoon, resulting in high clay content in the sediment. Previous studies related to past sea level changes in Sri Lanka have also suggested sea level transgression between 7000 and 4900 cal yr BP, but the sea level generally stabilized after this period (Weerakkody, 1992;Ranasinghe et al., 2013a).
The core average δ 13 C n-C33 value of PN is slightly higher than for BGSL (about 3‰) but varies within a narrow range (Table. 1). This is probably due to the relatively high abundance of C 4 plants in the dry zone in contrast to the wet zone. However, the values still indicate that the area was dominated by C 3 vegetation without a large shift in vegetation types (C 3 -C 4 transition) (Diefendorf and Freimuth, 2017;Holtvoeth et al., 2019).

Horton Plains
Today, the Horton Plains area in the Central Highlands is influenced by both SWM (summer) and NEM (winter) rainfall but with a dominant impact from SWM. Therefore, the HP peat deposit is expected to preserve mixed signals of both types of precipitation. This high-altitude peat deposit reveals a depositional history of ∼5500 cal yr BP with a rapid accumulation rate in the top ∼70 cm of the core. Moss and woody plants dominate between 5400 and 1800 cal yr BP (cluster 2), whereas grasses dominate from ∼1800 cal yr BP to the present (cluster 1) ( Figure 3C). PC2 scores indicate that the change of moss and woody plant-dominant peat to grassdominant peat occurs around 1700 cal yr BP ( Figure 4C).
The presence of C 28 and C 29 sterols, such as campesterol, stigmasterol and β-sitosterol, can provide additional insights into input from higher plants and mosses (Ronkainen et al., 2014). 5αstanols are the saturated counterparts of Δ 5 -sterols produced by anaerobic microbial hydrogenation (Routh et al., 2014). Hence, the Σ 5α-stanol: Σ Δ 5 -sterol ratio can be used to trace the degree of microbial hydrogenation or humification in peat (Routh et al., 2014;Naafs et al., 2019). The Σ 5α-stanol: Σ Δ 5 -sterol ratios in HP are very low (below 0.5) in the acrotelm (top-most ∼20 cm) that yielded modern ages (negative 14 C values). In the anaerobic cattotelm (below ∼20 cm), the Σ 5α-stanol: Σ Δ 5 -sterol ratio showed a downcore increasing trend, indicating that microbial hydrogenation increases down the core ( Figure 4C). TOC content also decreases downcore, further supporting the observation of post-diagenetic changes after deposition ( Figure 4C). Transformation of plant derived Δ 5 -sterols to their correspondent 5α-stanols preferentially occurs under anoxic conditions (Naafs et al., 2019). Therefore, low values of Σ 5α-stanol: Σ Δ 5 -sterol ratio suggest oxic conditions in the peat, probably due to a lowering of the water table during dry periods. From ∼5500 to 1800 cal yr BP, the Σ 5α-stanol: Σ Δ 5 -sterol ratio varies in a narrow range but shows very low values at ∼1200 cal yr BP, indicating more oxic conditions (low water table) related to the onset of a dry period.
In HP, both terrestrial plants and moss provide information on paleo-hydrological changes in the area in slightly different ways during the last ∼5500 cal yr BP. Grass (n-C 33 ) records the H isotopic ratios of acrotelm water originating from precipitation and is little affected by evaporation, whereas moss (n-C 23 ) records water inside its cells and between its leaves which is strongly affected by local evaporation (Nichols et al., 2010). However, the influence of summer and winter monsoon in the area, with distinct isotopic signatures, and their fluctuations in amount, makes the interpretations complicated.
Overall, a core upward increasing trend of moss-derived δD n-C23 and grass-derived δD n-C33 values from ∼5500 cal yr BP to present indicates overall decreasing trends of precipitation in the area. Both grass and moss hydrogen isotopic variations show very similar behavior from ∼5500 to ∼ 1250 cal yr BP indicating low evaporation. However, from ∼ 1250 cal yr BP to present they show somewhat different behavior and there is an asynchrony between their peak values. This is probably due to the relatively increased contribution of winter monsoonal precipitation (with a concomitant failure of summer monsoon) during this period as well as the different ways in which grasses and mosses record local hydrological conditions. Since, the moss provides more of an evaporation signal rather than source water signal (Nichols Frontiers (Berger and Loutre, 1994), (B) Total Solar Irradiance (TSI) (Steinhilber et al., 2009), (C) Ti content in Cariaco Basin sediments (Haug, 2001), (D) Number of ENSO events per 100 yrs (Moy et al., 2002), (E) Percentage of sand in El Junco Lake, Galápagos (Conroy et al., 2008), (F) Reconstructed North Atlantic Oscillation Index (Olsen et al., 2012), (G) Hematite-Stained Grains (HSG) percentage in North Atlantic deepsea sediment and Bond events (Bond et al., 2001).
Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 789291 14 et al., 2010), we can assume that the peak in aridity should be around 600 cal yr BP based on the values obtained.
The n-alkane isotope trends (δD n-C33 , δD n-C23 and δ 13 C n-C33 ), and relatively low Σ 5α-stanol: Σ Δ 5 -sterol ratios, indicate a failure of rainfall starting at ∼1250 and the onset of a relatively dry climate period between 1000 and 600 cal yr BP at HP. This is further supported by shifting dominant OM sources in the HP peat core, from moss and woody plants to grasses, which occurred around the same time according to the PC2 score variations.
However, we cannot identify the arid phase recorded in Bolgoda South Lake in Horton Plains despite both sites are dominantly influenced by SWM. This could be due to the HP is also significantly influenced by NEM in addition to SWM (see Supplementary Figure S1), which brings relatively isotopically depleted moisture contrast to SWM, that can mask the arid phase caused by the failure of SWM in HP.

South Asian Monsoon During the Middle and Late Holocene
Organic geochemical proxies, especially leaf wax hydrogen (and carbon) isotopes, lithogenic metal element ratios (Zr/Rb), and grain size distributions in Bolgoda South Lake indicate an overall decreasing trend of SAM (Southwest monsoon) in Sri Lanka between ∼5600 and 1400 cal yr BP. In contrast, the same proxies in the Panama lagoon indicate an overall increasing trend of South Asian Winter monsoon rainfall (Northeast monsoon) during ∼5900 cal yr BP to the present. Proxy records from the Horton Plains show more similar variations to the Bolgoda South Lake records, indicating an overall decreasing trend of summer monsoon rainfall ( Figure 5). Despite the fact that Bolgoda South Lake is missing records from ∼1400 cal yr BP to present, our previous study from North Bolgoda Lake (Gayantha et al., 2017;Gayantha et al., 2020) is able to confirm a weakening trend of summer monsoon rainfall in the area during this time (between 1400 cal yr BP-present). Therefore, we can reasonably assume that the weakening trend of summer monsoon rainfall identified in Bolgoda South Lake continues until the present.
The results from Bolgoda South Lake, Panama lagoon and the Horton Plains indicate a sub-millennial scale inverse relationship between South Asian summer and winter monsoon rainfall in Sri Lanka from the mid-Holocene onwards. This inverse relationship between summer and winter monsoon trends generally correlates with the Earth's precession (orbital forcing)-controlled inverse distribution of solar insolation between the Northern and Southern hemispheres during the middle and late Holocene ( Figure 5; Berger and Loutre, 1994). These changes in solar energy distribution primarily induced the southward migration of the mean latitudinal position of the ITCZ according to established records in different parts of the globe ( Figure 5; Fleitmann et al., 2003;Haug, 2001;Schneider et al., 2014). These show that the northernmost (in boreal summer) and southernmost (in boreal winter) positions of the ITCZ, moved southward from the middle Holocene. Being an island located in the path of the seasonal migration path of ITCZ, this phenomenon would have influenced the long-term (submillennial scale) inverse variation of the seasonal lengths, and thus the amounts of summer and winter monsoon rainfall, in Sri Lanka. However, our interpretations on this inverse trend between summer and winter monsoon precipitations are specifically for the SAM and not address the East Asian Monsoon. In addition to the broad patterns of the different monsoonal rains, we also observed multi-centennial scale failure in the summer monsoon between ∼3700-3200 cal yr BP, leading to aridity in the lowland wet zone of Sri Lanka (Figures 4, 5). Severe and prolonged droughts have been identified across various portions of Asia between approximately 4500 and 3500 cal yr BP leading to supposed social "collapse" at certain sites such as Harappa and Mohenjo-Daro (Staubwasser and Weiss, 2006;Sinha et al., 2011;Kathayat et al., 2017).
In addition to insolation and ITCZ-related processes, other possible mechanisms affecting the observed rainfall variations have also been suggested. An increased ENSO frequency during the Late Holocene, in contrast to the mid-Holocene (Moy et al., 2002;Conroy et al., 2008), may have significantly influenced the weakening of summer monsoon rainfall during the Holocene ( Figure 5). Generally, El Niño conditions in the Pacific Ocean cause the failure of Indian summer monsoon rainfall and results in drought conditions in South Asia including Sri Lanka (Rasmusson and Carpenter, 1983;Kumar et al., 2006;Zubair et al., 2007;Sinha et al., 2011). In contrast, a positive correlation has been identified between El Niño and NEM in Southern India (Yadav, 2012). However, only the second inter monsoon (SIM) rainfall (October-November) in Sri Lanka show a positive correlation with El Niño conditions in Sri Lanka (Kane, 1998;Malmgren et al., 2003;Hapuarachchi and Jayawardena, 2015). Therefore, increasing El Niño frequency during the Late Holocene ( Figure 5) could have diminished the summer monsoon rainfall in Sri Lanka. In addition, increasing the El Niño frequency can further drag the ITCZ southward, as observed during the LIA, that again causes a weakening of the Indian summer monsoon rainfall (Brown and Johnson, 2005;Banerji et al., 2019). However, according to the reconstructed ENSO records (Moy et al., 2002;Conroy et al., 2008), there is a more abrupt increase in the frequency of El Niño from middle to late Holocene rather than gradual variation. In addition, during the last ∼1000 yrs BP, it shows a decreasing trend ( Figure 5D). Therefore, it is obvious that increasing El Niño frequency could only have a partial influence on the observed inverse rainfall trends between the South Asian summer and winter monsoon. However, rapid increase of El Niño frequency between ∼600 and 1200 cal yr BP (Moy et al., 2002) possibly have led the arid phase recorded in Horton plains around ∼1000 cal yr BP ( Figure 4©).
The peak of the arid phase identified in this study, ∼3500 cal yr BP, coincides with a sharp negative excursion of the NAO index ( Figure 5; Olsen et al., 2012) and shows a partial overlap with Bond event number 2 ( Figure 5; Bond et al., 2001). Similarly, the onset of dry conditions at Horton Plains around 1000 cal yr BP also partially correlates with Bond event 1 ( Figure 5; Bond et al., 2001;Steinhilber et al., 2009). A negative NAO results in cold Frontiers in Earth Science | www.frontiersin.org December 2021 | Volume 9 | Article 789291 16 winters in northern Europe, causing expanding ice sheets and increased snow cover over Eurasia. Similarly, Bond events represent North Atlantic ice rafting events (∼1500 years cycle) related to cold periods in the Northern subpolar regions. These events are generally synchronous with periods of low solar activity during the Holocene (Bond et al., 2001;Banerji et al., 2020). These conditions also lead to failures in the development of a low-pressure trough in the northern and central parts of India, reducing South Asian summer monsoon rainfall over the Indian sub-continent (Bamzai and Shukla, 1999;Gupta et al., 2020). However, only Bond events 1 and 2 show at least partial correlation with the arid phases identified in our records. Whereas the rest of the Bond events or negative NAO conditions do not appear to lead to clear arid phases in our records. Therefore, it is unclear how exactly these arid events are induced by the interaction of different global climate variables and more high-resolution records are needed to answer this question.

Mid-Late Holocene Human-Environment Interactions in Sri Lanka
Many paleoclimate records available in Sri Lanka and India imply a general weakening trend of Indian summer monsoon during the middle and late Holocene (Prasad et al., 1997;Premathilake and Risberg, 2003;Ponton et al., 2012;Leipe et al., 2014;Prasad et al., 2014;Basavaiah et al., 2015;Ankit et al., 2017;Gayantha et al., 2017;Ratnayake et al., 2017;Gayantha et al., 2020, see Figure 6). The notable exception to this trend can be observed in the records from the South Indian peninsular region, which imply a progressive or stepwise increase in rainfall and/fluctuation in intensity in some cases (Veena et al., 2014;Rajmanickam et al., 2017;Sandeep et al., 2017, see Figure 6). However, these regions receive a significant amount of rainfall from the NEM (winter monsoon; between 17 and 49%) in addition to SWM (summer monsoon) rainfall (Sreekala et al., 2012). Indeed, the proxies used to reconstruct monsoon variability in these previous studies cannot distinguish between the source of monsoon rainfall (SWM or NEM). Regional records for winter monsoon rainfall variation are extremely limited. Some only cover the Late Holocene, while others show no clear trends. Nevertheless, in general, the available records show an overall increasing winter monsoon rainfall towards the present (Ranasinghe et al., 2013b;Böll et al., 2014;Mishra et al., 2019; see Figure 6). Our evidence from Sri Lanka allows us to contribute to existing work by isolating the impacts of the SWM and NEM in different sedimentary archives, supporting a decline in SWM precipitation and increase in NEM precipitation over the past 6,000 years in South Asia ( Figure 6).
Prehistoric farming populations seem to have emerged in the lowland dry zone of Sri Lanka around 2900 cal yr BP (Deraniyagala, 1992;Myrdal-Runebjer, 1996). Dense rainforest in the wet zone and mountainous regions, by contrast, seem to have only seen limited agricultural activities at this time . Flat terrain, savannah-type forests, and more stable NE monsoon rainfall may have attracted prehistoric and early historical agricultural societies to the lowlands of the dry zones. Notably, the extensive urban kingdoms of Anuradhapura and Polonnaruwa, associated with rice agriculture as well as the tending of orchards, emerged between ∼2,500 and 800 cal yr BP with sophisticated hydraulic technology in the lowland dry zone (Jayasena et al., 2011;Bebermeier et al., 2017;Abeywardana et al., 2019). Around 900-800 cal yr BP, however, these kingdoms appear to have begun to shift towards the southwest and central highland areas (towards the wet zone through the intermediate zones), with invasions and failure in the hydraulic systems both postulated as potential drivers (Deraniyagala, 1992;Gilliland et al., 2013). Unfortunately, our NEM records are not of sufficient resolution to clearly understand the influence of rainfall on the abandonment of ancient kingdoms in the dry zone during this period (see Lucero et al., 2015), necessitating further work. That said, we can see that SWM rainfall in Sri Lanka has been decreasing continuously until the present in the lowland wet zone (similar to the highland wet zone) due to the key influence of orbital forcing factors. This weakening of SWM rainfall might have led to a decline in rainforest density (shrinking of the rainforest areas) in the wet zone of Sri Lanka, turning the southwestern areas of Sri Lanka into landscapes more suitable for rice agriculture at a time when cities were relocating from the northeastern portions of the island (Lucero et al., 2015).

CONCLUSION
The inverse relationship between South Asian summer and winter monsoons in Sri Lanka supports the hypothesis that the Earth's precession controlled the southward migration of the mean latitudinal position of the ITCZ during the last ∼6000 cal yr BP, changing the distribution of interhemispheric solar insolation. The southward migration of the northernmost and southernmost positions of the ITCZ during this period increased the seasonal length of the winter monsoon in Sri Lanka while decreasing the seasonal length of the summer monsoon, resulting in an inverse relationship between summer and winter monsoon. The more abrupt weakening phases of rainfall is more pronounced in the summer monsoonal records, suggesting that the summer monsoon is more sensitive to other factors/teleconnections relative to the winter monsoon. However, systematic correlation between different climate variables (eg; NAO, Bond events and ENSO) and these arid events cannot always be identified in our records, suggesting these teleconnections are more complex and strongly influenced by regional factors. Nevertheless, the identified overall weakening trend of Southwest monsoon variability could have perhaps played an important role in social, economic, and political change witnessed across Sri Lanka between 2,500 and 9.00 years ago. Further well-dated high-resolution records are necessary to confirm and understand human-environment interactions in this part of the tropics during the Late Holocene, as well as the primary drivers of abrupt aridity events seen in the available records.