Nearly half of the Pearl River water discharges into the SCS through the PRE via three main tributaries in the lower drainage basin: North River, West River, and East River ( Figure 1 ). A 3-meter gravity core was collected at 21-m water depth off Guishan Island in the PRE (22.1315°N, 113.8055°E) in October 2017. The core was sectioned at an interval of 5cm for 0–100 cm and 10cm for 100–300 cm. The samples were immediately stored at -80°C until further analysis.

The total ²¹⁰Pb (²¹⁰Pb_t), ²²⁶Ra, and ¹³⁷Cs activities were analyzed in dry samples (6∼9g) with a low background high-purity germanium (HPGe) γ-ray detector (EG& G Ortec Ltd., USA) at the State Key Laboratory of Marine Geology, Tongji University. The samples were sealed in polyethylene tubes to allow for radioactive equilibration for 30 days before analysis (Ye et al., 2020).

The activities of ¹³⁷Cs were too low to obtain any confident results, therefore only ²¹⁰Pb dating was used in this study. ²²⁶Ra was used as an index of supported ²¹⁰Pb (²¹⁰Pb_su), and excess ²¹⁰Pb (²¹⁰Pb_ex) activities were calculated by subtracting ²¹⁰Pb_su activities from ²¹⁰Pb_t activities. The sediment accumulation rate (SAR) was then calculated with the constant rate of supply (CRS) model which assumes that the flux of ²¹⁰Pb to the accumulating sediment is constant during a timescale of 100–200 years. The chronology was derived by fitting the exponential ²¹⁰Pb decay profiles with the cumulative dry mass (Appleby and Oldfield, 1978; Zhang et al., 2014; Li et al., 2021). Fitting analysis was done using the ‘Exp2PMod1’ function in the Origin 2021 program.

Sediment samples are first floated in deionized water, dispersed through agitation, and then sieved to <180 microns to remove any root materials. Then the <180 microns organic sediment size fraction is treated with a series of hot (~ 70-90°C) acid leaches with HCl at a concentration of 0.1N for a period of 4-12 hours. Additional applications of HCl are provided until any carbonate presence has been completely removed. Samples are then rinsed to neutral with deionized water and dried at 100°C until dry. The sediments are then homogenized to insure an equal dispersion of the available carbon. A small aliquot is then tested with concentrated HCl to check for completion of carbonate removal. The sample was then measured either in whole or where applicable appropriately sub-sampled for combustion and analysis. The TOC (%), TN (%) content was determined on an elemental analyzer (Vaio EL Cube, Germany) after decarbonation. Stable isotopes of TOC (δ¹³C) were measured using an isotope ratio mass spectrometer (Thermo Science Delta Plus, USA) connected on-line to an elemental analyzer (Carlo Erba Instruments Flash 1112, USA). The C/N ratio was calculated as the atomic ratio of TOC and TN. The ¹⁴C of these samples was measured by accelerator mass spectrometer (AMS) interfaced with an elemental analyzer at the Beta Analytic testing laboratory, USA. Radiocarbon data were expressed as Δ¹⁴C values and fraction modern (Fm).

The grain size of the sediments was analyzed by a laser grain-size analyzer (Mastersizer 3000, Malven Instruments Ltd., UK) following the methodology described by Jiang et al. (2016). Briefly, about 0.2 g of the samples were treated with 15ml 10% (v/v) hydrogen peroxide to remove the organic matter. Carbonates were then removed by gradual addition of 15 ml of 10% HCl. The sample residue was dispersed with 10 ml of 10% (NaPO₃)₆ on an ultrasonic vibrator for 10 min before instrumental analysis. The particle sizes less than 4 μm were defined as clay, 4-63 μm as silt, and > 63 μm as sand. To sort the grain-size distribution into valuable information on geological processes and palaeo-environmental changes, end-member analysis (EMA) was applied to estimate end-member variations according to the methods by Prins et al. (2000). In this study, we used a newly developed GUI software of AnalySize for processing and unmixing grain size data (Paterson and Heslop, 2015) to determine the grain-size distributions of the detrital fraction in our sediment core. Porosity of each sample was calculated from the water content (wet-dry weight) prior to and after freeze-drying.

Lignin analyses were performed using CuO oxidation method of Hedges and Ertel (1982), as modified by Bianchi et al. (2002). Homogenized sediments containing ca. 3 to 5 mg of organic carbon were transferred to stainless-steel reaction vials with 330 ± 4 mg CuO and 3 to 5 ml 2 N NaOH in glove box purged by nitrogen and then digested at 150°C for 3 h. Reaction products were neutralized and extracted with diethyl ether (peroxides removed), filtered through combusted glass-fiber filled with glass wool, dried under N₂, and converted to trimethylsilyl derivatives using bis-(trimethylsilyl)-trifluoroacetamide (BSTFA). Lignin-phenol derivatives were analyzed with an Agilent 7890-Gas Chromatograph/5977-Mass Spectrometric Detector (GC–MS).

Quantification of lignin-phenols was based on a mixed standard calibration curve containing known amounts of 12 lignin reaction products of interest as well as the internal standard ethyl vanillin. Eight lignin-phenol oxidation monomer products (LOPs): C (ferulic acid+cinnamic acid), V (vanillin+acetovanillone+vanillic acid) and S (syringealdehyde+ acetosyringone+syringic acid) were quantified and used as molecular indicators for source and diagenetic state of terrestrial vascular plant tissue (Hedges and Parker, 1976). Compound of 3,5-dihydroxybenzoic acid (3,5Bd) was also derived after cupric oxidation and quantified (Goñi and Hedges, 1995). The precision for the total lignin phenols was within ±10%, while that for individual compound ranged from ±5 to ±15% based on triplicate analysis.

Lambda-8 (Λ₈),which is defined as the total weight in milligrams of the sum of C, V, and S phenols, normalized to 100 mg of organic carbon (Hedges and Parker, 1976), is used as a biomarker for terrestrial vascular plants. The acid-to-aldehyde ratios of both V and S phenols: (Ad/Al)v, (Ad/Al)s, were used as indicators of lignin degradation state prior to burial (Hedges et al., 1988). The C/V (woody/non-woody) and S/V (gymnosperm/angiosperm) ratios were plotted as indicators of the source of vascular plant (Hedges and Mann, 1979). The lignin-phenol vegetation index (LPVI) was also applied to study sources of vascular plant materials (Tareq et al., 2004; Sánchez-García et al., 2009). The 3,5Bd is used as an index of soil degradation processes, while the 3,5Bd/V indicated inputs of organic matter humification products sorbed to fine particles in soils (Houel et al., 2006).

A binary mixing model was used to resolve the non-rock-derived biospheric (OC_bio) and petrogenic (OC_petro) OC (Galy et al., 2008). Then the radiocarbon composition of the bulk OC can be expressed as follows:

where, Fm, Fm_bio, and Fm_petro are fraction modern values of bulk TOC, OC_bio, and OC_petro, respectively. OC_bio has different quantities of radioactive carbon (Fm> 0), while OC_petro does not contain radioactive carbon (Fm = 0) (Galy et al., 2008). So, equation (1) is further modified as:

Thus, Fm_bio (slope) and OC_petro (intercept/slope) can be estimated by plotting Fm ×TOC versus TOC on an X–Y plot, while OC_bio is estimated by the offset between bulk TOC and OC_petro.

A Monte Carlo simulation model was applied to track the sources of the sedimentary TOC from C₃, C₄ plants, riverbank soil, river phytoplankton, and marine algae. Assuming that the endmember parameters (δ¹³C and N/C) followed a normal distribution (mean ± standard deviation) for different OC sources in the study system ( Table S1 ), the program was run in Python 3.8.2. Briefly, 4000 out of 1,000,000 random samples from the normal distribution of each end‐member were taken in order to simultaneously optimize the following underdetermined system of linear equations:

where F_i is the fraction of the i end‐member and i = vascular C₃ plant, C₄ plant, soil OC, marine, and river phytoplankton, respectively. N/C, the inverse of C/N ratio, is used as a more sensitive end-member of terrestrial OC (Perdue and Koprivnjak, 2007; Li et al., 2017). The variation of the mean value for each end‐member was less than 0.2‰ by randomly sampling each parameter value five times, ensuring the statistical stability of the model.

The bulk carbon burial rate is calculated by the following equation:

where, SAR was determined from ²¹⁰Pb chronology; the bulk density was assumed to be 1.5 g cm^-3 in this region (Chen et al., 2006); TOC and porosity were the average values of 20-110 cm, 120-170 cm, 180-300 cm and the entire core.

The burial rate of each end-member was calculated by multiplying the bulk carbon burial rate and the fraction of each end-members from the mixing model. For example, the burial rate of lignin is calculated as:

The Origin 2021 software was used to graph the figures. Statistical differences were calculated using one-way ANOVA. Statistically significant differences were discussed within the 95% confidence interval. Principle component analysis (PCA) was performed to discriminate for any other controlling variables linked with bulk and biomarker patterns in sediment samples with all parameters being mean-normalized.

The cores displayed relatively low excess activities of ²¹⁰Pb ranging from 0.22 dpm g⁻¹ to 2.62 dpm g⁻¹ (1.33 ± 0.65 dpm g⁻¹, n=14) ( Figure 2 and Table S2 ). While grain size variation can generate an error in the downcore decay trend in relatively low excess activity samples, the core showed a supported ²¹⁰Pb level (from ²²⁶Ra activity) that were relatively invariant downcore (2.16–2.86 dpm g⁻¹) ( Figure 2 ). The surface mixed layer was around 20 cm with the depth below showing the sediment accumulation. The best fit exponential regression of excess ²¹⁰Pb activity (R² = 0.69, p<0.01) yielded a varied sedimentation rate (LSR) ranging from 1.71 cm yr^-1 (210-250 cm) to 3.18 cm yr⁻¹ (20-30 cm). The calculated geochronology dated back to 1893 for the deepest sample. Therefore, the core represents about 130 years of sediment deposition.

The TOC (%) varied from 0.43% to 1.43% (0.96 ± 0.22%, n=40) ( Figure 3A ). The TN ranged from 0.11% to 0.28% (0.24 ± 0.03%, n=40) ( Figure 3B ). A significant linear relationship between TN and TOC (R² = 0.49, p < 0.01) suggested that TN was derived predominantly from the organic origin. The C/N ratio varied from 2.27 to 9.33 (4.69 ± 1.18, n=40) ( Figure 3D ). The δ¹³C ranged from -24.54 to -22.89‰ and showed three stages with different δ¹³C features ( Figure 3C and Figure S1 ). Between 1893 and 1957, the average δ¹³C value was -23.08 ± 0.16‰. Between 1957 and 1980, the value showed a large variability (-23.59 ± 0.29‰), while after 1980, the average δ¹³C decreased to -23.86 ± 0.34‰ ( Figure S1 ), indicating an increased proportion of terrestrial OC. The Δ¹⁴C value varied from -448.40‰ to -209.55‰ (-299.26 ± 66.89‰, n=20) ( Figure 3E and Figure S1 ). The average Δ¹⁴C value was -226.84‰ (n=6), -399.04‰ (n=3), and -316.88‰ (n=9) ( Figure S1 ) in stage of 1893 to 1957, 1957 to 1980, and after 1980, respectively. Correspondingly, the average Fm in stages of 1893 to 1957, 1957 to 1980, and after 1980 were 0.78 (n=6), 0.61 (n=3), and 0.69 (n=9) ( Figure S1 ), respectively.

The clay, silt and sand content of the core ranged from 11.58 to 35.6% (26.24 ± 5.82%, n=40), 60.36 to 74.56% (66.54 ± 3.34%, n=40), 0 to 16.2% (7.22 ± 4.63%, n=40), respectively ( Figure 3F and Figure S1 ). The sediments were predominately fine-grained (<63 μm) that the sum of clay and silt ranged from 83.8% to 100%. The proportion of clay and silt were negatively correlated in all layers (R² = 0.61, p<0.01).

The median grain size (MGS) ranged from 5.83 to 26.9μm (9.88 ± 4.33μm, n=40) ( Figure 3F ). Before the 1960s, the MGS was relatively constant at a low value (6.90 ± 0.74, n=15), while they became more variable and showed relatively higher values (11.83 ± 4.61, n=25) after the 1960s. The MGS has a peak value around 2008. The correlation map between multiple correlation coefficient and numbers of end-member (EM) indicated that two EMs could fulfill the observed compositional variation required in EMA (Jiang et al., 2017) and explain the grain-size distribution pattern that the peak values were concentrated at 6.72 μm (EM 1), and 40.14 μm (EM 2), respectively. The EM1 and EM2 ranged from 36.91% to 97.42% (78.43 ± 14.78%, n=40), 2.58% to 63.09% (21.57 ± 14.78%, n=40), respectively ( Figure S2 ).

The Λ₈ values ranged from 0.43 to 2.84 mg 100mg^-1 OC (1.39 ± 0.54, n=40) ( Figure 4A ). There was a general increase in Λ₈ from the bottom to the surface sediment, that the average Λ₈ in stages of 1893 to 1957, 1957 to 1980, and after 1980 were 1.04 (n=13), 1.18 (n=6), 1.68 (n=21) mg 100mg^-1 OC, respectively. This trend indicated an increasing accumulation of vascular plant derived OC. The 3,5Bd showed no significant changes with depth (p > 0.01) ( Figure 4B ). The average (Ad/Al)v values were 0.38 (n=13), 0.35 (n=6) and 0.26 (n=21), while the average (Ad/Al)s were 0.35 (n=13), 0.36 (n=6) and 0.23 (n=21) from 1893 to 1957, 1957 to 1980, and after 1980, respectively ( Figures 4C, D ). Their insignificant depth variations (p<0.01) indicated insignificant degradation, or that the majority of bioavailable lignin has been consumed prior to delivery to the ocean (Seidel et al., 2015). The ratios of 3,5Bd/V and P/(S+V) ranged from 0.03 to 0.17 (0.08 ± 0.03, n=40), 0.11 to 0.3 (0.17 ± 0.05, n=40) with slight increase with depths ( Figures 4E, F ).

The modelling results showed historical variability in the Fm_bio, OC_petro (%), OC_bio (%), OC_bio/OC (%) in the PRE sediments ( Figure 5A and Table 1 ). The Fm_bio (0.79) values in 1957-1980 were smaller than before 1980 (0.63). Conversely, the OC_petro value (0.03) in 1957-1980 was three times higher than before 1980 (0.01), indicating more petrogenic and less biospheric OC input during 1957-1980.

The comprehensive five-end-member simulation showed that the fraction of C₃ plants (F_C3), riverbank soil (F_soil), riverine phytoplankton (F_{riverine phytoplankton}), marine algae (F_marine) and, C₄ plants (F_C4) were 17 ± 3%, 17 ± 2%, 22 ± 1%, 35 ± 10%, 10 ± 3%, respectively ( Figure 5B and Table 2 ). The fraction of total terrestrial OC (F_terr= F_C3+F_soil+F_{riverine phytoplankton}+F_C4) was 65 ± 10%. The F_C3 in 1893-1957 was 14 ± 3%, which was smaller than 1957-1980 (20 ± 2%) and 1980-2016 (19 ± 2%). The F_C4 in 1980-2016 was lower than 1957-1980(13 ± 2%) and 1897-1957(10 ± 4%). This indicated a change in plant contribution from C₄ plants to C₃ plants from the bottom to the surface sediment. The F_soil (19 ± 2%) in 1957-1980 were larger than the average value (17 ± 2%), while the F_marine (27 ± 6%) in 1957-1980 was smaller than the average value (35 ± 10%) in the entire core.

The PCA biplot included 17 normalized variables ( Figure 5C ). The first two components explained 52.5% of the sedimentary OC (PC1: 36.1% and PC2: 16.4%). The plot showed significant depth variations of the distribution and preservation of sedimentary OC across different temporal scales ( Figure 5B ). Most of the samples above 110cm (exclude 30cm) were located in Quadrant II and III, and the samples between 120 and 180cm (exclude 180cm) were located in Quadrants III and IV, while most of the deeper (190-300cm) samples were located in Quadrants I and IV.

The bulk OC burial rate increased from 1893-1957 (94 ± 22 gC m^-2 yr^-1) to 1957-1980 (138 ± 24 gC m^-2 yr^-1), but decreased after 1980 (151 ± 44 gC m^-2 yr^-1) ( Table 3 ). The trends of the end-member burial rates were similar to the bulk OC, which were increased first and then decreased except the burial rate of marine algae ( Table 3 ).

The significant correlation between δ¹³C and C/N (R² = 0.55, p< 0.01, Figure 6A ) suggested mixed source input from terrestrial and marine OC. The MGS and δ¹³C signatures were invariable in the bottom of the sediment core (> 150cm), but became coarser and more negative at discrete intervals in upper parts of the core (above 150cm) ( Figures 3C, F ) indicating a shift from marine plankton with younger ¹⁴C age to terrestrial derived OC which was older ( Figures 3E , 5A ). This was likely due to a transition (~ the 1980) from steady-state deposition with lower sediment accumulation rate ( Figure 2 and Table S2 ), to an environment influenced by both stronger anthropogenic and climate disturbance in the modern PRE (Yuan et al., 2019) with enhanced sediment accumulation rate ( Figure 2 and Table S2 ) (Owen and Lee, 2004).

More than 100 red tides have been reported in the PRE since the 1970s, and their frequency has increased in recent years due to human influences (Jia and Peng, 2003; Hu et al., 2008). However, the average marine sourced OC derived from the Monte Carlo model decreased by 20% during 1980-2016 than that from 1893-1957 ( Table 2 ). This may be related to the high degradation rates of marine algae enriched in labile compositions such as carbohydrates, sugars, amino acids, and low-molecular-weight organic acids (Hardison et al., 2013). Instead, the vascular plant biomarker of lignin did not show much degradation during the burial process ( Figures 4C, D ) that the variations were more induced by concentration changes with continuous input of relatively stable terrestrial OC, which were partly diluted by the marine-derived OC as was indicated from a lack of correlation between lignin‐phenol abundance and TOC (R² = 0.02, p<0.01).

Sources of vascular plant derived OC changed during the sediment accumulation. In general, the angiosperms (woody or non-woody tissues) contain a large amount of syringyl phenols. Therefore, the bottom depths with relatively high S/V values (0.83 ± 0.14, n=40) ( Figure 6B ) suggested an important fraction of lignin originating from angiosperm plants with higher LPVI ( Table S3 ) than that of gymnosperms (Tareq et al., 2004). The surface sediments, instead, tended to be sourced from leaves, humus, and soil ( Figure 6B ) as a result of increased soil erosion induced terrestrial plant input. The average LPVI decreased from the bottom to the surface sediment also supported a change from angiosperms to gymnosperms (Tareq et al., 2004; Li et al., 2017) that rapid urbanization has caused such loss in the angiosperm-rich farmland since the 1980s (Zhang et al., 2008) ( Table S3 ).

The sediment record showed a significant input of old carbon during 1957-1980 ( Figure 3 ; Figure 7 ), which corresponded to lower Fm and OC_bio values and higher Δ¹⁴C and OC_petro values ( Table 1 ). It happened that this period had a higher TOC value ( Figure S1 ) and terrestrial input (73%, Table 2 ). This may be caused by hydrodynamic sorting of riverine OC containing significant eroded soil with more negative Δ¹⁴C and older apparent ages (Wei et al., 2010). The significant correlation between MGS and Δ¹⁴C values (R² = 0.52, p<0.01, Figure S3 ) illustrated the importance of hydrodynamics on the transport and sorting of different grain size fractions, which should be considered in the aging processes, particularly in the estuary and continental shelf characterizing by strong hydrodynamic gradients (Bao et al., 2018a). Specifically, OC in the finer fraction was the youngest, with ¹⁴C ages and OC_petro increasing in the coarser fraction. So, the OC was preferentially associated with fine-grained, large-surface-area minerals providing stronger protection against degradation (Ausín et al., 2021). The young OC with less negative δ¹³C values ( Figure S1 ) were mainly from marine inputs during 1893-1957. Selective removal of labile marine OC was expected leading to more negative Δ¹⁴C and δ¹³C values, via preferential degradation of ¹³C-enriched, labile organic components afterwards. Correspondingly, the finer EM1 (6.72 μm) was dominant (78.4 ± 14.8%, n=40) ( Figure S2 ) in the bottom of the PRE sediments, which was beneficial to the OC preservation. Thus, the ¹⁴C age was younger with higher EM1.

The average Fm_bio value (0.71) of the PRE is the lowest in coastal China among Bohai, Yellow Sea, and East China Sea ( Table 4 ), but closer to that of the Yellow River POC (0.61), indicating older ages of biogenic OC in the PRE. The Yellow River generally has a large amount of pre-aged POC sourced from the Loess Plateau where distributes thick loess-paleosol deposits, with serious soil erosion and sparse vegetation (Eliassen, 2020). In comparison, the upper and middle reaches of the PRE drain the karst morphology areas with the carbonate rocks being the dominant lithology (Wu et al., 2020). The weathering rate of carbonate rocks is more than an order of magnitude higher than that of silicate rocks (Meybeck, 1987). As the carbonate rock weathers, the pre-aged materials including dissolved inorganic carbon are released into the river (Liu et al., 2017). Phytoplankton thus synthesize such old dissolved inorganic carbon resulting in biogenic OC with older ages and low Fm_bio values ( Table 4 ). The relatively lower Fm_bio values have also been found in some other fluvial systems such as in the Mackenzie River shelf ( Table 4 ). Warming and associated permafrost thaw exposed older biogenic OC with less OC_petro and lower Fm_bio accumulated in marine sediments (Hilton et al., 2015). Therefore, the relatively low Fm_bio composition in these aquatic ecosystems were attributed to the contribution from pre-aged OC in the drainage basin. As a result, the OC_petro contributed to 2% of TOC in the PRE. The ratio can be as low as 1% in the East China Sea sand area, and as high as 87% in the southeast Alaskan fjords ( Table 4 ). But on average, the OC_petro/TOC in the PRE was lower than the mean values of 13 ± 18% in global aquatic systems.