Long-term patterns of peat accumulation and organic matter decomposition in Costa Rican peatlands

Tropical peatlands, compared to their boreal counterparts, are vastly understudied despite acting as a significant terrestrial carbon sink, sequestering 100—300 Gt of carbon. In particular, the low number of field-based studies from Latin America and the Caribbean limits our knowledge of these important wetland ecosystems. Across the tropical Panamerican region, peatland location, soil characteristics, inception ages, and carbon accumulation histories remain largely unknown. These datasets are needed to inform a mechanistic understanding of why peat develops in certain areas but not in others, both in terms of peat initiation conditions as well as the factors that enable peat to subsist over centuries and millennia. Here we present extensive, high-resolution laboratory datasets from 11 peat cores from 4 peatland types from Costa Rica (high-elevation, riverine, coastal palm swamp, and mangrove). A multi-proxy palaeoecological approach was employed to shed light on the successional pathways and past conditions that have allowed these peatlands to form, as well as to provide a first estimate of their carbon stock. The core characterization includes radiocarbon dating, loss-on-ignition, carbon and nitrogen content, and plant macrofossils. Fourier transform infrared spectroscopy (FTIR) was also used to assess changes in organic matter quality across sites and over time. The averaged peatland carbon stock in Costa Rica is estimated at 1080 MgC ha⁻¹, making these ecosystems exceptionally rich carbon stores that are comparable to values found in lowland Amazonian peatlands. Overall, this research provides a basis for understanding long-term carbon accumulation within Caribbean tropical peatlands.

1 Introduction

Tropical peatlands are estimated to occupy around 14 %, and potentially up to 30 %, of the global peatland area (Gumbricht et al., 2017; UNEP, 2022; Xu et al., 2018). Carbon stock estimates for tropical peatlands similarly vary by a factor of about two, from ∼ 120–130 Gt (Dargie et al., 2017; Leifeld and Menichetti, 2018) to 288 Gt (Ribeiro et al., 2021). The reason for these large ranges is mainly attributable to the scarcity of field-based datasets pertaining to (1) peatland distribution and surface area, (2) peat depth, and (3) peat carbon density (Gumbricht et al., 2017; Melton et al., 2022; Page et al., 2011). At times, the tropical peatland carbon pool has been estimated on the basis of peat carbon density datasets from extra-tropical regions (Yu et al., 2010), adding unquantifiable uncertainties that pertain to differences in vegetation assemblages, peat formation processes, and environmental conditions. Similarly, averaged bulk density and carbon concentration values from Southeast Asian peatlands have previously been extrapolated to the Panamerican region, highlighting the lack of local information (Page et al., 2011). These data and knowledge gaps undermine our ability to quantify the tropical peatland soil carbon pool and assess the past, present, and future roles of these ecosystems in the global carbon cycle.

In the tropical rainforest biome, the existence of peat has often been referred to as “surprising”. This idea stems from the notion that tropical terrestrial ecosystems tend to be characterized by rapid rates of microbial decomposition due to consistently hot and humid conditions (Knorr et al., 2005; Meentemeyer, 1978). The accumulation of peat in the tropics therefore requires specific conditions that either dampen the intensified decomposition process near the soil surface and/or a carbon-soil storing mechanism that bypasses the aerobic decay zone. Previous studies have proposed the following conditions and mechanisms, which are not mutually exclusive: (1) peat soils must be flooded (or waterlogged all the way up to the near surface) for long enough throughout the year to induce anoxia, which limits organic matter decay and allows for peat buildup (Chimner and Ewel, 2005; Gillman et al., 2015); (2) plant litter inputs to the peat must be very high to maintain a positive balance between litter production and peat decomposition (Chimner and Ewel, 2005); (3) the peat litter is composed of relatively recalcitrant compounds and/or lignin-rich tissues that inhibit decomposition (Couwenberg et al., 2010; Dommain et al., 2015; Hodgkins et al., 2018; Hoyos-Santillan et al., 2015; Phillips et al., 1997; Verbeke et al., 2022; Wright et al., 2013); and (4) root biomass inputs (down into the constantly anoxic zone) circumvent the intensive near-surface peat decomposition processes and lead to peat buildup (Chimner and Ewel, 2005).

The bulk of the literature on tropical peatland carbon storage dynamics stems from Southeast Asia (Cobb et al., 2024; Dadap et al., 2022; Dommain et al., 2015; Hergoualc'h and Verchot, 2011; Hoyt et al., 2020; Ruwaimana et al., 2020; Sasmito et al., 2025) and, to a lesser extent, the Congo Basin (Crezee et al., 2022; Dargie et al., 2017; Garcin et al., 2022; Young et al., 2023). In the lowlands of Latin America and the Caribbean (LAC), a region that may encompass the largest tropical peatland complexes in the world (Gumbricht et al., 2017), regional knowledge on peatland distribution, extent, and dynamics remains sparse. Across this region, the most extensive work has been performed in Western and Central Amazonia (e.g., Dargie et al., 2024; Draper et al., 2014; Hastie et al., 2022; Lähteenoja and Page, 2011; Lähteenoja et al., 2009, 2012, 2013; Lawson et al., 2014, 2026; Roucoux et al., 2013; Winton et al., 2025), leaving Central America and the Caribbean particularly understudied. One exception is the Changuinola peat deposit in Bocas del Toro, Panama, where detailed paleoecological information reveals the historical peatland developmental pathways (Cohen et al., 1989; Phillips, 1995; Phillips and Bustin, 1996; Phillips et al., 1997; Swindles et al., 2024). A few studies have also looked at carbon and water dynamics in the Changuinola peatlands (Baird et al., 2017; Girkin et al., 2025; Hedgpeth et al., 2025; Hoyos-Santillan et al., 2016; Sjögersten et al., 2020). But beyond Changuinola, there is little knowledge on where peat is found, how deep, old, and carbon-rich it might be, and what factors control long-term peat accumulation processes in the Caribbean region (for a review, refer to Rabel and Loisel, 2024).

In this study, we present high-resolution, multi-proxy analyses from a series of peat cores that were collected across nine sites that span four distinct hydro-ecoclimatic and geomorphic settings. The main goal of the study is to characterize peat accumulation dynamics across these diverse sites and decipher spatial and temporal differences across these peatlands. First, we provide the paleoecological history from four of our study sites to document changes in the major peat-forming plant macrofossil types and humification to identify successional pathways. Second, we present peat basal ages for our nine study sites to assess the timing of peat inception in Costa Rica. Third, we provide a suite of geochemical measurements (organic matter content, dry bulk density, carbon and nitrogen content, $\mathrm{C}/\mathrm{N}$ ratio) to estimate carbon storage. Lastly, through the use of infrared spectroscopy, we endeavor to improve our understanding of the peat formation process itself. By correlating our plant macrofossil and von Post humification analyses with Fourier-Transform infrared spectroscopy (FTIR) data, we identify shifts in organic matter quality that coincide with successional changes and/or humification.

2 Study area and general approach

While Costa Rica is rarely mentioned as a peat-bearing country, a recent probabilistic map suggests that it may harbor 1456 km² of peat, equivalent to 3 % of its land area (Rabel et al., 2025; Fig. 1). It is important to note that Costa Rica does not yet possess an official peatland inventory. In July and August 2023, our team traveled across Costa Rica and visited over 15 sites (Loisel et al., 2024) that had been identified on the basis of published wetland and vegetation maps as well as few existing peat depth point data (SINAC, 2018b; Villegas-Mejía, 2018). We corroborated those data with global-scale peatland probability maps (Gumbricht et al., 2017; Melton et al., 2022) and a preliminary map that was produced by Peters and Tegetmeyer (2019). In most places that were visited, we found peaty soils (i.e., organic-rich soil layers in excess of 30 cm in thickness, sensu Lourenco et al., 2023). Our site selection was ultimately based on site accessibility as well as regional representativity of Costa Rica's peatland cover, such that we made sure to sample the following peatland types: coastal peat swamps, mangroves, riverine complexes, and montane peatbogs (Jiménez, 2016).

Figure 1 Locations and photos of sites from this study. (a) The country-wide Costa Rica map; the dark green area represents the peatland area (Rabel et al., 2025). Photos of the sites include: (b) the coastal palm swamp Gandoca (GAN), (c) the Red Mangrove site (RM), (d) the riverine wetland complex Medio Queso (MQ), and (e) the montane Bog 70 peatland. (f) A climate envelope for Costa Rican peatlands, where each circle represents a pixel that has been classified as peat (modified from Rabel and Loisel, 2024); AIR: Airport site; CN: Caño Negro site. Source for panel a: Esri, TomTom, FAO, USGS | Powered by Esri. Photo credits for panels (b)–(e): Hannah Mitchell and Patrick Campbell.

Costa Rican peatlands are found under a broad range of mean annual precipitation and precipitation seasonality (Rabel et al., 2025; Fig. 1), and different peatland types are delineated within this climate space. As such, the main peatland types mentioned above are also distinct in terms of their climatology: the Caribbean coastal sites (palm swamps and mangroves) receive high annual precipitation amounts with low precipitation seasonality (and annual precipitation increasing along a north-to-south latitudinal gradient along the coast), whereas the (inland) riverine sites are intermediate in both annual precipitation and its seasonality, which contrast with the montane sites that are characterized by lower annual precipitation with higher precipitation seasonality.

2.1 Study sites

The work herein stems from an exploratory trip; we opted for visiting a larger number of sites while obtaining fewer data points for each site. The general approach that was used to sample sites included: (1) trace a transect line from the edge of a potential peatland towards its center, (2) walk (or boat) along said transects, and (3) probe the soil along the transects to confirm the presence (or absence) of peat. The length and orientation of these transects varied for each site and mainly depended on site conditions and time constraints. Peat depth and coordinates were recorded for each probing.

The sections below describe the four main peatland types that were sampled and provide information for each coring site. Table 1 offers a summary of the most relevant information.

Table 1 Summary information for the nine Costa Rican peatland study sites. Mean annual temperature and annual precipitation represent 30-year averages (1970–2000) (data from Fick and Hijmans, 2017). Full cores consist of the entire peat column, from the peatland surface to the peat-to-mineral interface. Bottom cores consist of the deepest 50 cm drive that was collected at the site; they include the peat-to-mineral interface. These bottom sediment cores are primarily used for basal age determination.

2.1.1 Coastal palm swamps

These three sites are located along the southern Caribbean coast of Costa Rica, between the towns of Gandoca and Limón. All three sites are dominated by palm swamp species Raphia taedigera and harbor other plants, shrubs, and tree species, such as Campnosperma panamensis (Cortés, 1998; Peters and Tegetmeyer, 2019). These coastal sites have developed in back-barrier topographic lows, in close proximity to the coastline (Fig. 1). The Gandoca site (GAN) is located in the southern part of the country, close to the Panama border; two peat cores were collected along a 1 km-long transect that crosses the site (GAN-4; GAN-SF). The Manzanillo site (MAN) is located approximately 10 km north of GAN and it is part of the Gandoca-Manzanillo Wildlife Refuge and Ramsar site #783 (Ramsar, 2026). Lastly, the Airport site (AIR) is located just south of Limón, in Moín. These three sites were accessed by foot.

To our knowledge, the only studies of coastal palm swamps in Costa Rica stem from prospective work in the 1980s and 1990s that provided preliminary peat depth and calorific values and was aimed at assessing the potential economic development of these peat deposits as fuel resources (Obando et al., 1995; Thayer et al., 1995). Some botanical surveys have also been performed (e.g., Jiménez, 2016).

2.1.2 Mangrove peatlands

The Red Mangrove site (RM) is located ∼ 100 m away from the Gandoca site (Fig. 1). These two ecosystems are separated by Laguna Gandoca, a meandering stream that reaches the Caribbean Sea. The northern side of the stream supports peat-forming Raphia taedigera palm swamps, whereas the southern shore is instead characterized by brackish water, tides, and Rhizophora mangle communities.

As only 1 % of Costa Rica's mangrove forests are located along the Caribbean coast, our RM site presents a rare opportunity to study one of the largest (and preserved) mangrove stands along the Caribbean coast (Cortés, 1998, 2016). There are local studies on the floristic composition of these mangroves (Álvarez-Sánchez and Piedra-Castro, 2020) and a physico-chemical characterization of the Gangoca Lagoon (Coll et al., 2004), but otherwise there is little information on mangrove soil properties from this region with the exception of a study pointing to high soil organic matter (Rovai et al., 2018).

2.1.3 Riverine peatlands

These three sites are located in the north-central portion of Costa Rica, near the Nicaraguan border (Fig. 1). The riverine sites all constitute peatland/wetland complexes that experience relatively large water table level fluctuations on an annual basis, with surface flooding that can reach 2 m during the wet season (Moreno-Casasola and Warner, 2009; Pérez-Castillo et al., 2024). The “Medio Queso” site (MQ) is located approximately 1 km away from the town of Los Chiles and is found within the watershed of the Rio Medio Queso. This wetland complex is used for pasture and it is frequently burned to maintain land clearing. Recent studies at the site evaluated the impacts of fires on vegetation and surface peat chemistry (Camacho-Navarro et al., 2017; Pérez-Castillo et al., 2024). The MQ study site is primarily colonized by rushes (Eleocharis interstincta), though many sedge and small bush species are also found (Obando and Malavassi, 1993), along with small “islands” of Cuba palm (Acoelorraphe wrightii). Two peat cores (MQ-9; MQ-10) were collected along a 500 m-long transect that runs from the western edge of the site towards Río Medio Queso.

Our other two riverine study sites are located within a lagoon complex that is primarily fed by Ríos Frío and Mónico. The lagoon is part of the Caño Negro Wildlife Refuge and Ramsar site #541 (Ramsar, 2026), which is classified as a mixed-use protected wild area that is home to local communities who make sustainable use of the ecosystem through subsistence and cultural connections. The two study sites (San Jeronimo (SJ) and Laguna Martin (LM) are located within 1 km of one another; the vegetation communities are dominated by Eleocharis equisetoides and Scleria microcarpa (SINAC, 2018a). The SJ and LM sites were accessed by boat and the peat cores were extracted from the boat platform, as the coring sites were flooded by approximately 1m of water at the time of core collection (August 2023). We did not find any published soil information about Caño Negro, though there exists a botanical guide that was published by the Government of Costa Rica (SINAC, 2018a).

2.1.4 Montane peatbogs

These two sites are located in the Cordillera de Talamanca, Ramsar site #1286 (Ramsar, 2026), in central Costa Rica. In this high-elevation area (∼ 3000 m a.s.l.), peatlands tend to be very small in size (often <0.05 km²; Rabel et al., 2025). Their existence has been linked to the mountain range's Late Pleistocene glacial history, with montane peatlands forming in small, poorly drained depressions following the last deglaciation around 10 000 years ago (Horn, 1990). The small peatland-bearing depressions are prone to seasonal flooding (Jiménez, 2016). Today, the area is considered part of a subalpine paramo ecosystem, itself located within a cloud forest characterized by a cool and moist climate. The two sites that we sampled, Bog 70 and Bog 68, were named after the km markers on the side of the road (Rodgers and Horn, 1996). These sites are dominated by a unique mixture of peat moss (Sphagnum sp.), reindeer lichen (Cladonia rangiferina), and tropical plants such as tree ferns (Blechnum auratum), bromeliads (Puya dasylirioides), and numerous herbaceous species (see Rodgers and Horn (1996) for details). Sediment cores have been retrieved from these sites in the 1990s and used for pollen-based paleoenvironmental reconstructions (Hooghiemstra et al., 1992; Islebe et al., 1995, 1996), but studies of the peat characteristics themselves were not conducted, to our knowledge.

3 Methods

3.1 Field sampling strategy

In the field, a “full core” was retrieved from the deepest portion of four peatlands using a Russian-style peat borer equipped with a 50 cm long, 5 cm diameter chamber (Table 1). A full core consists of the entire peat column, from the peatland surface to the peat-to-mineral interface. The multi-proxy analyses presented in this study were performed on these four cores. In addition to these full cores, several “bottom cores” (n=7) were collected along transects within all 9 peatlands (Table 1). These bottom cores consist of the deepest 50 cm long drive that was collected at the site and they include the peat-to-mineral interface. These bottom sediment cores are primarily used to analyze patterns in peat basal ages along transects, which may help understand the history of peatland expansion across individual basins. In the field, all freshly collected cores were placed in PVC pipes and wrapped in plastic film and aluminum foil. The samples were kept at 4 °C until laboratory analysis.

3.2 Laboratory preparation and analysis

In the laboratory, all peat cores were photographed and sliced into contiguous, 1 cm-thick slices, which were then placed into labeled plastic bags and stored at 4 °C. The following sections describe the methodological approach that was followed to prepare our samples for different proxy analyses. Although not every analysis was conducted at the same resolution, most measurements were made on the same peat slices to allow for direct comparison across methods (see Table A1 for a list of proxies and analytical resolution).

3.2.1 Peat geochemistry (loss-on-ignition and elemental analysis)

To determine peat water content, dry bulk density, and organic matter content, the loss-on-ignition method (LOI) was followed (Dean, 1974). Every 2 cm, a sub-sample (1 cm³) was extracted from the peat slice, weighed, and sequentially oven-dried to determine water content and dry bulk density (40 °C, 48 h) as well as organic matter content (550 °C, 4 h). These measurements were also used to calculate organic matter density.

To quantify carbon and nitrogen content, 1 cm³ sub-samples were collected every 8 cm along the cores. Samples were oven-dried (40 °C) and ground into a fine powder using a Retsch MM 300 Vibration Mill Mixer (Verder Company, Germany). The homogenized samples were then placed into tin capsules. Carbon and nitrogen content was measured using a Costech ECS 4010 Elemental Analyzer (Costech Analytical Technologies, Inc., Valencia, CA USA). Different sample weights were used on the basis of organic matter content, with smaller sample weights for organic-rich samples and vice versa. The carbon-to-nitrogen ratio ($\mathrm{C}/\mathrm{N}$) is presented as a mass ratio. The $\mathrm{C}/\mathrm{N}$ ratio of peat can be used as a proxy for decomposition and humification (Biester et al., 2014; Kuhry and Vitt, 1996; Malmer and Holm, 1984). Organic carbon density was subsequently calculated for each sample by multiplying the organic matter density (from LOI data) with the measured carbon content (Loisel et al., 2014).

3.2.2 Peat organic matter quality (FTIR spectroscopy)

FTIR spectroscopy was used to analyze peat organic matter quality (Artz et al., 2008; Holmgren and Nordén, 1988; Krumins et al., 2012; Uhelski et al., 2022). FTIR spectroscopy has been shown to be a useful proxy for humification index in a peat column (Biester et al., 2014; Broder et al., 2012; Tfaily et al., 2014; Upton et al., 2018). Sub-samples (1 cm³) were collected every 8 cm along the four full cores. Samples were oven-dried (40 °C) and ground into a fine, homogenized powder using a Retsch MM 300 Vibration Mill Mixer. Spectra were acquired on a Thermo Nicolet 380 FTIR spectrometer (Thermo Electron Scientific Instruments Corporation, Madison, WI USA) by averaging 100 scans at 8 m⁻¹ resolution over a range of 4000–550 cm⁻¹. Samples were placed directly on an attenuated total reflectance (ATR) crystal, and light force was applied so the samples came into direct contact with the crystal. ATR correction and atmospheric suppression took place to account for variability in background CO₂ and water vapor. A background scan was collected before every sample. The spectra were then manually baseline-corrected in the regions between 4000–3750 and 2600–1800 cm⁻¹. A previously published R script was used to determine the relative abundances of various organic functional groups in predetermined regions of the spectra (Hodgkins et al., 2018). With the baseline correction performed prior to peak height calculation, raw peak heights normalized by total spectral area were used as estimates of relative abundance of each moiety in the sample, facilitating the comparison of peak intensities by measuring them all from the same baseline within a given spectrum.

Many FTIR spectra from our sample set showed evidence of silicate minerals in addition to peat organic matter, as evidenced by prominent features at ∼ 3620, and 3700 cm⁻¹ (see Fig. 12). Because the Si-O feature at ∼ 1005 cm⁻¹ overlaps substantially with the C-O peak at ∼ 1030 cm⁻¹, the relative abundance of carbohydrates cannot be distinguished from silicate minerals in our samples. The peak intensity at ∼ 1005 cm⁻¹ is therefore defined as “Carbohydrates + Silicates”. Aromatic relative abundance was determined from the sum of the intensities of the peaks at ∼ 1510 and ∼ 1615 cm⁻¹, and aliphatic relative abundance from the peaks at 2850 and 2920 cm⁻¹. These relative abundances are not absolute concentrations but instead provide semi-quantitative estimates of contributions from each moiety that can be used to compare relative quantities and changes across the sample set.

3.2.3 Peat composition (plant macrofossils and von Post humification index)

To characterize peat composition, plant macrofossils were analyzed following previously described procedures (Loisel et al., 2009; Mauquoy et al., 2002). Sub-samples (1 cm³) were analyzed every 8 cm (or at higher resolution near stratigraphic boundaries) along each full core. Samples were mixed with deionized water and suspended in a petri dish, then scanned under a stereomicroscope. The first step involved the visual estimation of the percentage of mosses, herbaceous, ligneous, mineral, and unidentified plant material (UOM) (Mauquoy et al., 2010). A peat type category (bryophyte, herbaceous, ligneous) was assigned to each sample on the basis of the plant macrofossil analysis, where the dominant plant type found in each sample was used to infer the peat type (Loisel et al., 2014). We note that the “ligneous” peat type was split into two sub-categories: ligneous-general vs. ligneous-mangrove to better represent our peat material. Lastly, we specified the presence of “root material” when possible. This differentiation is often challenging, given structural similarities between aerial vs. subterranean plant tissue. But in the case of mangrove roots, the distinction was safe to make.

To estimate the degree of peat decomposition, the von Post humification scale (ASTM, 1987) was used at 8 cm increments. This widely-used approach estimates peat humification by pressing and manipulating a small sample of peat (∼ 1 cm³) in one's hands (Kaila, 1956; Rydin and Jeglum, 2013; Stanek and Silc, 1977).

3.2.4 Peat basal ages and chronologies (radiocarbon dates and age-depth models)

Peat inception ages were determined using radiocarbon (¹⁴C) dating (Table 2). Along each core, the deepest (and presumably oldest) peat sample was identified using LOI results. A 30 % organic matter threshold was used to identify the peat-to-mineral boundary along our peat cores (Joosten and Clark, 2002; Lourenco et al., 2023). Given the humic nature of the peat, dated samples consisted of bulk, root-free samples (63–125 µm). Samples were wet-sieved using deionized water, oven-dried at 40 °C, and packaged in aluminum foil. Additional samples were extracted along the four full cores and similarly prepared for ¹⁴C dating, to generate peat accumulation histories. Samples were shipped to Lawrence Livermore National Laboratory's Center for Accelerator Mass Spectrometry (CAMS), where they underwent acid-base-acid (ABA) treatment (Gillespie and Hedges, 1983; Hedgpeth et al., 2025; Norris et al., 2020). Dates were first individually calibrated to calendar age (cal. yr BP) using the Intcal20 calibration curve (Reimer et al., 2020) in CALIB v8.2 (Stuiver and Reimer, 1993). As for constraining peat chronologies of the full cores using multiple ¹⁴C ages, a Bayesian age-depth modeling approach was used (package rbacon) v.3.5.2 in R (Blaauw and Christen, 2011).

Table 2 AMS radiocarbon dating results from the new peat cores in Costa Rica.

* The modern age for sample CAMS-195269 returned a fraction modern of 1.0520 ± 0.0041.

3.2.5 Statistical Analysis

First, an analysis of variance (ANOVA) was used to test the effect of “peat type” on the peat geochemical properties (Loisel et al., 2014). Second, a correlation matrix was completed in R with the packages energy, corrplot, dplyr, and ggplot2 to compare the relationship between peat proxies (Rizzo and Szekely, 2024; Wei and Simko, 2024; Wickham, 2016; Wickham et al., 2023). To account for possible non-linear trends in the data, we opted for a Spearman correlation (Zar, 2005). Third, a Principal Coordinates Analysis (PCoA) based on the Bray-Curtis dissimilarity was used to capture variations in peat properties across peat types and sites, and over time (Tian et al., 2019). The analysis was performed in R with the vegan package (Oksanen et al., 2025). Bray-Curtis distances were chosen to understand differences in ecological compositions while also handling the differences in relative abundances gathered from the FTIR analysis (Kellerman et al., 2015; Wang et al., 2022). Lastly, PERMANOVA and Distance-based redundancy (dbRDA) analyses were performed in conjunction with the PCoA to statistically test differences between sites (Anderson, 2001; Legendre and Anderson, 1999; Tian et al., 2019).

4 Results and Discussion

4.1 Geochemical characterization of Costa Rican peatlands

4.1.1 Peat properties

Combining all the loss-on-ignition measurements (n=279) and elemental analysis data (n=70) across our four cores provides an overview of peat water content, dry bulk density, organic matter content, organic matter density, carbon and nitrogen content, and $\mathrm{C}/\mathrm{N}$ ratio for Costa Rican peatlands. Table 3 and Fig. 2 present the data ranges as well as the mean values and standard deviations for all these variables. One-way analyses of variance (ANOVA) revealed that all peat types were statistically different from each other in the case of water content ($F_{\mathrm{2},\mathrm{276}}=\mathrm{11.94}$, p<0.0001) and organic matter content ($F_{\mathrm{2},\mathrm{276}}=\mathrm{44.66}$, p<0.0001). In the case of dry bulk density, the ANOVA revealed an overall effect of peat type ($F_{\mathrm{2},\mathrm{276}}=\mathrm{4.95}$, p<0.01), with the bryophyte type significantly different (Tukey's LSD: p<0.001) from herbaceous and ligneous peats. Lastly, for organic matter density, we also report an overall effect of peat type ($F_{\mathrm{2},\mathrm{276}}=\mathrm{18.38}$, p<0.0001), with the herbaceous peat type significantly different (Tukey's LSD: p<0.001) from the bryophyte and ligneous peat types. Statistical tests were not applied to %C, %N, and $\mathrm{C}/\mathrm{N}$, given the lower number of samples available. We also note that the correlation between %C and %OM in the Costa Rican samples was strong and statistically significant (Spearman correlation = 0.847; p<0.0001, Fig. 3). The slope (0.51) of the linear regression indicates that %C is about half the OM content, similar as the relationship found in the northern peatland database (0.49; Loisel et al., 2014), supporting the idea that OM content and %C are closely related in peatlands. With that said, the intercepts differ between this new dataset and the one from the northern peatland database (9 vs. 25, respectively), pointing to a generally higher %C in northern peatlands.

Table 3 Descriptive statistics for Costa Rican peat geochemical properties. The dataset is compared with values from a northern peatland compilation (shown in parentheses) (Loisel et al., 2014).

Figure 2 Frequency distribution histograms and downcore profiles of peat properties in Costa Rican peatlands. (a) Water content (%), (b) Bulk density (g cm⁻³), (c) Organic matter content (%), (d) Organic matter density (g cm⁻³), (e) Nitrogen content (%), (f) Carbon content (%), (g) Carbon/nitrogen ratio, and (h) Organic carbon density (gC cm⁻³).

Figure 3 Relationship between organic matter content (%OM) and carbon content (%C) in Costa Rican peat.

Combining the geochemical data together by peat type (bryophyte, ligneous, herbaceous) provides the first dataset of its kind for Costa Rica (and the Caribbean more broadly), with direct usability for regional carbon analysis and modeling (Dehaen et al., 2026; Swails et al., 2026; Fig. 4). Bryophyte peat (n=10, samples are only from the montane site) is characterized by high organic matter content (86 ± 14 %), the highest of all peat types. Herbaceous peat (n=84, with data coming from our montane and riverine sites) yielded moderately high organic matter content (75 ± 17 %). Costa Rican herbaceous peats are about 50 % more dense than their northern counterparts, in terms of both bulk density and organic matter density (BD: 0.23 ± 12 and OMD: 0.15 ± 0.05 g cm⁻³ in Costa Rica vs. 0.12 ± 0.08 and 0.09 ± 0.04 northern peatlands). Carbon content in herbaceous Costa Rican peats were lower (38.5 ± 10.3 %) than in northern peats (50.5 ± 4.9 %). Ligneous peat (n=185, with data from the riverine, coastal swamp, and mangrove) yielded lower organic matter content (56 ± 17 %). Bulk density and organic matter density from ligneous peats was higher in Costa Rica than in northern peats (BD: 0.25 ± 12 and OMD: 0.13 ± 0.05 g cm⁻³ in Costa Rica vs. 0.11 ± 0.05 and 0.1 ± 0.03). Carbon content was almost double in northern peatlands (50.9 ± 4.0 %) than in Costa Rica (26.7 ± 9.8 %).

Figure 4 Box plots of geochemical data by peat type. (a) Organic matter content (%), (b) Bulk density (g cm⁻³), (c) Organic matter density (g cm⁻³), (d) Carbon content (%), and (e) Nitrogen content (%).

4.1.2 Costa Rican peatland carbon stock

To estimate peat carbon storage, previous studies have often reverted to a formula that combines peat extent, average peat depth, average bulk density, and an assumed carbon content of 50 % or 51 % (Gorham, 1991; Loisel and Yu, 2013; Yu et al., 2010). Here, we applied this formula and used a peatland extent of 1456 km² (Rabel et al., 2025), an average peat depth of 150 cm (which is the average depth for the 9 sites presented in this study), an average bulk density of 0.24 gC cm⁻³, and an average carbon content of 30 %. These values yielded a peat carbon store of 1.573 × 10¹⁴ gC (or 0.1573 GtC) for Costa Rica.

While the approach described above is simple, we argue that it provides a good first-order estimate of the Costa Rica carbon stock. Of course, every one of the values used above carries uncertainties, starting with the peatland extent, which was performed using a probabilistic approach (Rabel et al., 2025) that may under- or over-estimate the extent. Illustrating this uncertainty is the range of peat extent for Costa Rica that was provided in previous (global) peat mapping efforts, from 577 to 2670 km² (Melton et al., 2022; UNEP, 2022). Using these peatland extents results in carbon stores ranging between 0.0623 and 0.2883 GtC for Costa Rican peatlands. An important improvement brought about by our new dataset is that previous peat-carbon store estimates often use 0.1 g cm⁻³ as their bulk density in combination with an assumed 50 % carbon content, yielding an organic carbon density of 0.05 gC cm⁻³. The results for Costa Rica suggest an organic carbon density of 0.072 gC cm⁻³ (0.24 g cm⁻³ × 30 %C), about 30 % greater than the value previously used in carbon assessment (e.g., Yu et al., 2010). While this result should be further confirmed through additional measurements, it cautions against the broad application of well-established knowledge from extra-tropical peatlands to discuss peat-carbon characteristics in tropical peatlands. Estimates created out of these values would greatly underestimate the tropical carbon store. We also note that other authors have documented the higher density of tropical peatland soils before, including palm swamps in Colombia (0.19 g cm⁻³; Winton et al., 2025).

At the plot scale, Costa Rican peatlands were found to store 1080 MgC ha⁻¹ on average (150 cm × 0.24 g cm⁻³ × 30 %C), which is a massive concentration of belowground carbon. The riverine site (MQ) was characterized by the greatest carbon stock (1652 MgC ha⁻¹), followed by the peat-rich mangrove (RM; 1164 MgC ha⁻¹), the (young) peat swamp (GAN; 706 MgC ha⁻¹) and the montane peat bog (BOG70; 646 MgC ha⁻¹). The average value of 1080 MgC ha⁻¹ is 5–10 times greater than the typical carbon stock of a tropical rainforest (e.g., Baker et al., 2024; Malhi et al., 2009). When compared to other peatland carbon stores, the Costa Rican peatlands rival with the Brazilian veredas (∼ 1200 MgC ha⁻¹; Verona et al., 2026), the Peruvian Amazon lowland peat swamp forests (∼ 234–1067 MgC ha⁻¹; Bhomia et al., 2019), peatland pole forests (∼ 1133 MgC ha⁻¹; Honorio Coronado et al., 2021), and open peatlands (∼ 655 MgC ha⁻¹; Draper et al., 2014). Similar results have been reported from the Colombian lowland peatlands (490–1230 MgC ha⁻¹; Winton et al., 2025).

4.2 Peat basal ages and local patterns of peatland expansion in Costa Rica

Peat inception ages range from 507 to 10 967 cal. yr BP (Fig. 5) and tend to cluster on the basis of peatland type/geomorphological setting. The montane peatlands exhibited the oldest ages, with BOG68 dating back to 10,967 ± 248 cal. yr BP. Conversely, the youngest peatlands were observed along the Caribbean coast, with the coastal palm swamp near Limón (AIR site) only dating back to 507 ± 24 cal. yr BP. The inland riverine sites encompass a wider age range, between 988 ± 64 and 8361 ± 96 cal. yr BP. Combining all peat basal ages from this study suggests that most peatlands in Costa Rica are very young (less than 2000 years old), but this finding is partly due to a sampling artefact, as most of the sites are located in the coastal area and/or in proximity to river channels.

Figure 5 Peat inception ages from Costa Rican peatlands. Peat basal age vs. depth is displayed, with peatland types shown in different symbols and colors.

The timing of peatland inception has been used in many studies worldwide to help constrain changes in hydroclimatic conditions (Dommain et al., 2014; Thomas et al., 2025; Yu et al., 2016). In the case of Costa Rica, the dataset shows a wide range of peat inception ages, with no particular time period during which peats primarily initiated. We have not found peats of pre-Holocene age in Costa Rica, which is a common theme across Latin America and the Caribbean (e.g., Lawson et al., 2026). Overall, the spatiotemporal pattern of peat inception primarily relates to landscape-scale conditions, with high-elevation sites being the oldest and coastal sites being the youngest.

While the age of the montane peatlands was likely constrained by deglacial processes (Horn, 1990), the relatively recent inception and expansion of peatlands along the Caribbean coast (∼ 1000 cal. yr BP) was probably linked to increased stability in base level, which combines stable sea level conditions as well as lowered variability in freshwater and sediment inputs from the inland rivers that feed the coastal lowlands, providing conditions that support sediment deposition and peat formation (e.g., Dommain et al., 2014). Regional sea-level stabilization is thought to have taken place around 6000 cal. yr BP (Castañeda-Posadas et al., 2022), suggesting that other local factors might have prevented peatland establishment or development along the Costa Rican coast for many thousands of years, such as dynamic floodplains, earthquakes, erosion caused by storms, etc. Across the Caribbean region, and albeit a relatively low number of inception ages available in the literature (n=31), coastal peatland initiation has taken place from ∼ 9000 to 1000 cal. yr BP, with no discernable spatial pattern (Malpica-Piñeros et al., 2026; Rabel and Loisel, 2024). As for the riverine sites, our limited sampling does not allow us to generalize, but our three sites (MQ, SJ, LM) lay within a lowland area that was flooded by Lago Nicaragua (also known as Lago Cocibolca) during the late Pleistocene (Fig. 1; Bergoeing and Protti, 2006). Holocene reduction in lake size and fluctuations in the lake level, which could have been caused by a changing climate as well as regional tectonic activity (Kutterolf et al., 2023; Slate et al., 2013), have likely caused large shifts in the regional hydrological regime and associated changes in the size of the adjacent wetland complexes that span hundreds of km². Alas, the general lack of paleodata from the region makes it difficult to further address the cause and timing for peatland inception, though a 5700-year-old record from Lago Nicaragua reveals a potential impact of pre-Colombian agriculture (5400 cal. yr BP) on the region as well as a general drying trend throughout the mid- to late-Holocene bodies (Slate et al., 2013) that may have favored wetland formation to the detriment of freshwater lake.

Two of our study sites combine a full core and a bottom core, allowing us to provide some preliminary information on the spatio-temporal pattern of peatland expansion in Costa Rican peat basins. In the coastal site Gandoca (GAN), two cores were retrieved towards opposite ends of a 1 km-long transect that runs perpendicular to the shore. The full core (GAN) is located close to the beach and is characterized by 130 cm of peat that is aged at 668 cal. yr BP. Near the other end of the transect is core GAN-SF, with 170 cm of peat that dates back to 863 cal. yr BP. These results point to a somewhat homogeneous back-barrier depression that could have filled up more or less synchronously. This preliminary finding is supported by our probe peat depth data points along the GAN transect (n=10), which are all between 130 and 170 cm. Interestingly, we found a similar depth pattern at the AIR site, with peat depth values around 140 cm (n=7).

In riverine wetland complex Medio Queso (MQ), 2 cores were also collected along a 600 m-long, west-east transect from the edge of the site towards Rio Medio Queso. The full core (MQ; depth = 180 cm, age = 7683 cal. yr BP) and the bottom core (MQ-10; depth = 230 cm; age = 8361 cal. yr BP) are located approximately 100 m apart, with MQ-10 being closer to the river. At this site, deeper (and likely older) peats are found along the edge-to-center transect, as corroborated by the probing depths (n=10), which start at 40 cm near the edge of the site and reach 230 cm near the river, suggesting a progressive infilling of a local depression found in a fluvial setting.

4.3 Peat accumulation histories across four Costa Rican sites

4.3.1 Coastal palm swamp site (Gandoca – GAN)

The GAN core is characterized by 130 cm of peat that is underlain by fine mineral particles. Peatland initiation began at 668 cal. yr BP and the plant macrofossil record, which is dominated by ligneous palm remnants, suggests that a palm swamp community has likely occupied this site since peat inception (Fig. 6). The peat humification index is consistently high, with slightly more humified peat towards the base of the peat column. The apparent rate of peat accumulation is steady and relatively rapid (0.2 cm yr⁻¹), and the age-depth relationship follows a linear trend. The peat stratigraphy is high in organic matter throughout the profile (mean = 70.2 %). Expectedly, bulk density tends to be relatively low (mean = 0.164 g cm⁻³), with a sharp increase that corresponds to the mineral interval found underneath the peat (not shown). As for peat quality, the GAN core is characterized by generally increasing relative abundances of aromatic and aliphatic carbon with depth and, conversely, generally decreasing abundance of carbohydrates (+ silicates), with some variability in near-surface samples (16 and 24 cm) and towards the mineral transition at the end of the peat (136 cm) (Figs. 6 and A1).

Figure 6 Accumulation history, macrofossil composition, geochemical, and FTIR results for the coastal palm swamp (GAN) site. (a) Age-depth model, (b) Plant macrofossil panels include: herbaceous macrofossil composition (%), ligneous macrofossil composition (%), unidentified organic matter (UOM) composition (%), and von Post humification. (c) Loss-on-ignition panels include: organic matter content (%), bulk density (g cm⁻³), (h) organic matter density (g cm⁻³). (d) Elemental analysis panels include: carbon content (%), nitrogen content (%), $\mathrm{carbon}/\mathrm{nitrogen}$ ratio. (e) Fourier Transform Infrared Spectroscopy data include: carbohydrate + silicate abundance, aromatic abundance, aliphatic abundance.

4.3.2 Mangrove site (Red mangrove – RM)

The RM core is 187 cm long and primarily composed of organic-rich material. The record begins as a palm swamp that initiated around 1020 cal. yr BP and persisted until 465 cal. yr BP (Fig. 7). The underlying sediment consists of sandy material, and the stratigraphy is interrupted by one mineral-rich layer at ∼ 740–620 cal. yr BP (110–134 cm). Around 465 cal. yr BP (83 cm), the stratigraphy abruptly shifts to a dominance of Rhizophora plant remnants (mangrove roots). We caution that, while the shift from palm swamp to mangrove seems to take place at 465 cal. yr BP according to our Bayesian age-depth model, said model also rejected a ¹⁴C date from a sample that was picked at the ecological transition (83.5 cm) but that did not fit the rest of the curve per the principle of superposition and the assumption of a constant rate of accumulation used in the Bayesian approach (¹⁴C age: 220 ± 30, which yields a median age of 189 cal. yr BP with a 2-sigma range of 0–310 cal. yr BP; Fig. 7). This age reversal, if representative of the “true” sediment history, would suggest that the “mangrove root deposit” is in fact much more recent than indicated by the age-depth relationship, possibly as young as a few decades only.

Figure 7 Accumulation history, macrofossil composition, geochemical, and FTIR results for the mangrove (RM) site. (a) Age-depth model, (b) Plant macrofossil panels include: herbaceous macrofossil composition (%), ligneous macrofossil composition (%), mangrove root macrofossil composition (%), unidentified organic matter (UOM) composition (%), and von Post humification. (c) Loss-on-ignition panels include: organic matter content (%), bulk density (g cm⁻³), (h) organic matter density (g cm⁻³). (d) Elemental analysis panels include: carbon content (%), nitrogen content (%), $\mathrm{carbon}/\mathrm{nitrogen}$ ratio. (e) Fourier Transform Infrared Spectroscopy data include: carbohydrate + silicate abundance, aromatic abundance, aliphatic abundance.

The peat humification index is generally mesic to humic, with the exception of relatively undecomposed peat following the transition from palm swamp to mangrove. The apparent rate of peat/sediment accumulation is steady and relatively rapid (0.2 cm yr⁻¹; similar to the nearby GAN core), and the age-depth relationship follows a linear trend –with the caveat of a possible hiatus at the transition from palm swamp to mangrove (83 cm), as mentioned above. The peat stratigraphy is characterized by moderate organic matter content (mean = 43.83 %), the lowest average of all our sites. The RM results are characterized by wide changes between OM-rich (mean = 47.45 %) vs. mineral-rich (mean = 18.80 %) intervals, generally in a reverse relationship with bulk density. These changes in organic matter content and bulk density follow a ∼ 200-year cycle (Fig. 7). As for peat quality, the increases in aromatic and aliphatic relative abundances were not as evident as they were in the GAN core due to the mineral layer affecting samples from 112–128 cm; only at the bottom of the core did the aliphatic and aromatic abundances appreciably increase (Figs. 7 and A2). Aromatic and aliphatic abundances tend to correlate with carbon and organic matter contents, with greater abundances corresponding to OM- and C-rich material.

4.3.3 Riverine site (Medio Queso – MQ)

The MQ core presents 180 cm of peat material with an inception age of 7683 cal. yr BP. Initially, the site was likely a palm swamp (Fig. 8), as shown by the presence of ligneous palm fragments akin to those found at our coastal (GAN) and mangrove (RM) sites. The palm swamp ecosystem persisted for at least two millennia (until ∼ 5500 cal. yr BP), but perhaps for a longer time interval. A hiatus from ∼ 5500 to ∼ 2500 cal. yr BP (75–55 cm) challenges our ability to determine the timing of the transition into the modern-day herbaceous peatland. The apparent long-term rate of peat accumulation for the palm swamp (7683–5500 cal. yr BP) is 0.05 cm yr⁻¹, which is an order of magnitude slower than our observations at GAN and RM. The peat swamp portion of the profile is characterized by relatively high OM % (>50 %) and low BD (<0.3 g cm⁻³), with the exception of two layers with ∼ 30 % OM and densities of ∼ 0.5 g cm⁻³. (Fig. 8).

Figure 8 Accumulation history, macrofossil composition, geochemical, and FTIR results for the riverine (MQ) site. (a) Age-depth model, (b) Plant macrofossil panels include: herbaceous macrofossil composition (%), ligneous macrofossil composition (%), unidentified organic matter (UOM) composition (%), mineral composition (%), and von Post humification. (c) Loss-on-ignition panels include: organic matter content (%), bulk density (g cm⁻³), (h) organic matter density (g cm⁻³). (d) Elemental analysis panels include: carbon content (%), nitrogen content (%), $\mathrm{carbon}/\mathrm{nitrogen}$ ratio. (e) Fourier Transform Infrared Spectroscopy data include: carbohydrate + silicate abundance, aromatic abundance, aliphatic abundance.

There is a long time interval from 5500 to ∼ 2500 cal. yr BP (75 to 55 cm) during which a two-order-of-magnitude slowdown in apparent rate of peat accumulation is recorded (0.007 cm yr⁻¹). This interval coincides with high bulk density values (>0.4 g cm⁻³) and extremely decomposed peat (10 on the von Post scale). These conditions may indicate secondary decomposition processes (i.e., reactivation of diagenetic processes following a disturbance) and/or a drastic change in the net carbon balance of the peatland during that time interval, which resulted in little to no net peat accumulation for two millennia. Our field notes also indicate the presence of thin mineral layers throughout this part of the profile, which could pertain to flooding. The uppermost peat section (top 55 cm) corresponds to the past 2500 years and is composed of herbaceous material with an averaged apparent rate of peat accumulation of 0.02 cm yr⁻¹. This herbaceous section is characterized by noticeably lower $\mathrm{C}/\mathrm{N}$ values and the presence of herbaceous plant remnants (up to 20 %). We note that the uppermost section of the record (0–18 cm) was compressed during the coring process due to the presence of surface water, making it difficult to assess changes in peat properties in the near-surface samples. Nevertheless, OM % is highest in the uppermost peat layers (60 %–95 %) and corresponds to low bulk density values (<0.2 g cm⁻³), despite the peat being already very humified near the top of the core (von Post = 8).

4.3.4 Montane site (BOG70)

The BOG70 core is our shortest but oldest record, at 100 cm-long but with an initiation age of 9759 cal. yr BP. The ecosystem likely started as an herbaceous-dominated peatland (Fig. 9) atop clayey sediment. Ligneous fragments were found starting around 1823 cal. yr BP (42 cm). Sphagnum and Amblystegiaceae remnants appear very recently in the record, at −65 cal. yr BP which corresponds to 2015 calendar years (18.5 cm); these mosses form the modern-day moss carpet. The apparent long-term rate of peat accumulation was the slowest of all 4 study sites at 0.01 cm yr⁻¹, which is 15–20 times lower than the coastal palm sites. Peat humification is high along the entire record though it trends towards lower humification towards the top of the profile (top ∼ 20 cm). Organic matter content is high throughout the majority of the core (mean = 86.25 %); this montane site presents the highest organic matter content of all 4 study sites. Our results are corroborated by a prior study that reported high organic matter content (60 %–96 %) for these montane sites (Rodgers and Horn, 1996). The FTIR data show that aromatic and aliphatic relative abundances increase with age down the core (Figs. 9 and A4). Notably, high aliphatic abundance appears contemporaneous to the herbaceous peat. Similarly, our previous work in Patagonian peat demonstrated positive correlations between herbaceous macrofossils and aromatic abundance (Leri et al., 2025). In the spectra from samples at the surface of the core, the two aliphatic C-H bands at 2850 and 2920 cm⁻¹ are less resolved than in all the other samples (Fig. A4, top sample), likely indicating comparatively higher content of carbohydrates and other labile organic components that become degraded further down the core. The attenuated aliphatic bands are a hallmark of fresh Sphagnum (Heller et al., 2015). At the base of the peat column (105 cm), an increase in clay content is inferred from spectral features at 3620 and 3695 cm⁻¹, representing contributions from silicate minerals, and a sharp decrease in the relative abundance of aromatic and aliphatic components.

Figure 9 Accumulation history, macrofossil composition, geochemical, and FTIR results for the montane (BOG70) site. (a) Age-depth model; the grey vertical band between 5–20 cm signifies a section of greater uncertainty. In this case, the reason is that the modern age (at 18 cm) forces the age-depth curve to substantially change its accumulation rate and thus the slope between said modern age and the surface, when compared to the slope between the basal peat and the age at 18 cm. (b) Plant macrofossil panels include: herbaceous macrofossil composition (%), ligneous macrofossil composition (%), Sphagnum macrofossil composition (%), Other moss macrofossil composition (%), mineral composition (%), and von Post humification. (c) Loss-on-ignition panels include: organic matter content (%), bulk density (g cm⁻³), (h) organic matter density (g cm⁻³). (d) Elemental analysis panels include: carbon content (%), nitrogen content (%), $\mathrm{carbon}/\mathrm{nitrogen}$ ratio. (e) Fourier Transform Infrared Spectroscopy data include: carbohydrate + silicate abundance, aromatic abundance, aliphatic abundance.

4.4 Spatiotemporal patterns in peat properties

4.4.1 Peat organic matter quality (FTIR data)

Correlations across all sites (Figs. 10, 11) show carbohydrates + silicates negatively correlated with aromatics (Spearman = −0.881, p<0.0001), aliphatics (Spearman = −0.794, p<0.0001), and acids (Spearman = −0.679, p<0.0001). Aromatic and aliphatic relative abundances positively correlate with %C and %OM (Figs. 10, 11). All these associations are stronger at the site level (Fig. 10), and downcore trends are complex. For instance, the coastal swamp site (GAN) and the montane site (BOG70) both show decreases in carbohydrates + silicates with depth inverse to the increases in aromatic and aliphatic abundances. Assuming that carbohydrates dominate the signal at ∼ 1030 cm⁻¹, this pattern is expected, as carbohydrates represent comparatively labile organic matter, whereas aromatic and aliphatic moieties tend to show greater long-term stability and are thought to “accumulate” in older, more refractory peat (Leri et al., 2025; Verbeke et al., 2022). Aliphatic abundances are greater at the montane, herbaceous-dominated site (BOG70) than in the lowland, ligneous-dominated swamp (GAN). As for the riverine (MQ) and mangrove (RM) sites, large variations along the core do not follow the expected increase in recalcitrant fractions with depth.

Figure 10 FTIR data from our four Costa Rican peatlands. (a) The relative abundances of carbohydrates + silicates are plotted against aromatics, aliphatics, and acids, (b) Changes in carbohydrates + silicates, aromatics, aliphatics, and acids are shown as time series (i.e., along peat core depth).

Figure 11 Spearman Correlation Matrix to establish statistical relationships across sample variables. Carb + sil: carbohydrate + silicate abundance; Aro: aromatic abundance; Acid: acid abundance; Aliph: aliphatic abundance; Humif: von Post humification index; Water: water content; Bd: bulk density; Omc: organic matter content; N: nitrogen; C: carbon; CN: $\mathrm{Carbon}/\mathrm{Nitrogen}$ ratio.

Average spectra from all of the depths in each core were compared (Fig. 12) in an attempt to detect site-specific signatures. As mentioned above, the coastal peat swamp (GAN) and the montane site (BOG70) exhibit similarities in terms of their averaged carbohydrate + silicate and aromatic contents, though the herbaceous site (BOG70) is characterized by greater aliphatic content. Silicates are evident in all four cores and appear particularly intense at the mangrove (RM) and riverine (MQ) sites, as evidenced by the kaolinite bands (3620 and 3695 cm⁻¹; Kloprogge, 2018; Müller et al., 2014; Seaton et al., 2024; Sengyang et al., 2015). Given that the mangrove and riverine sites are prone to flooding, it is perhaps not surprising to see this abundance of mineral compounds. The mangrove site also displayed comparatively less intense aromatic (1510 and 1630 cm⁻¹) and aliphatic (2850 and 2920 cm⁻¹) bands, which could indicate a lower degree of humification compared with the coastal swamp and the montane peatbog.

Figure 12 Fourier Transform Infrared Spectroscopy in peatlands. (a) FTIR spectra from our four full cores from Costa Rica. (b) Example of an FTIR spectrum and peak finding for a peat sample. The major functional compounds are identified after baseline correcting and normalization within the custom R script (Hodgkins et al., 2018). Wavenumbers relate to carbon compounds discussed in the text.

4.4.2 Differences in peat composition and properties across sites

Peat composition and properties vary between and across the four sites (Figs. 13 and A5; Table A2). The first PCoA axis, explaining 79.8 % of the observed variation, places carbon content as well as the aromatics, aliphatics, and acids on the left side of the ordination vs. the carbohydrates + silicates and bulk density on the right side of the ordination. The second axis (12 %) was not strongly linked to any environmental variables. The PCoA plot shows that the peat properties from most samples from the montane peatbog (BOG70) and the coastal peat swamp (GAN) sites can be distinguished from those of the other two sites. In contrast, the mangrove (RM) and riverine (MQ) sites, regardless of peat type, depth, or age (Fig. A5), were scattered across the horizontal axis. Lastly, the PCoA plot that has the samples color-coded by depth (Fig. A5) shows a cluster of the shallower samples towards the left-hand side of the statistical space, suggesting the role of age on peat properties.

Figure 13 Principal Coordinates Analysis (PCoA) using Bray-Curtis distances overlaid with envfit vectors. The data in the ordination are visualized by peat type and core. Additional plots can be found in Fig. A5. Envfit vector abbreviations: carb + sil: carbohydrate + silicate abundance; bd: bulk density; cn: $\mathrm{carbon}/\mathrm{nitrogen}$ ratio; c: carbon content; acid: acid abundance; aliph: aliphatic abundance; aro: aromatic abundance; omc: organic matter content; water: water content; n: nitrogen content.

Notably, bulk density and carbohydrates + silicates covary (Fig. 13, on the right-hand side of axis 1). While carbohydrate relative abundance alone would be expected to negatively correlate with bulk density, the ordination instead suggests a strong positive relationship between bulk density and the intensity of the signal at ∼ 1030 cm⁻¹. We attribute this result to contributions from Si-O bands in mineral-rich samples, giving rise to an association between bulk density and the peak intensity at ∼ 1030 cm⁻¹. This finding demonstrates that caution must be applied in the interpretation of FTIR spectra from peat samples with clay mineral content to ensure that the band at ∼ 1030 cm⁻¹ is not taken to represent carbohydrates alone when there are clay minerals present.

Subsequent PERMANOVA analysis (Table A3) revealed differences between sites (R²=0.434, $F_{\mathrm{3},\mathrm{65}}=\mathrm{16.64}$, p≤0.001) and among peat types (R²=0.282, $F_{\mathrm{3},\mathrm{65}}=\mathrm{8.52}$, p≤0.001), explaining 43.4 % and 28.2 % of the variation in the Bray-Curtis ordination, respectively. When tested together, site and peat type explained 50.1 % of the variation (R²=0.501, $F_{\mathrm{6},\mathrm{62}}=\mathrm{10.38}$, p≤0.001). The homogeneity of dispersion (betadisper), which shows the internal variation within each group, was determined for both sites ($F_{\mathrm{3},\mathrm{65}}=\mathrm{5.21}$, p=0.003) and peat types ($F_{\mathrm{3},\mathrm{65}}=\mathrm{1.63}$, p=0.19, Fig. A6). Dispersion was shown to be significantly different across sites, with the mangrove (RM) and the riverine (MQ) sites displaying greater heterogeneity, but not significantly different among peat types. With that said, there could be covariance with uncharacterized parameters within those sites.

Lastly, a distance-based redundancy analysis (dbRDA) constrained by site and peat type was run to constrain the community variables (Table A3). Site and peat type together ($F_{\mathrm{6},\mathrm{62}}=\mathrm{8.253}$, p≤0.001) explained 44.4 % of total dispersion in the Bray-Curtis ordination. The majority of the variation was explained by the sites themselves ($F_{\mathrm{3},\mathrm{62}}=\mathrm{14.02}$, p≤0.001), with “peat types” adding explained variance ($F_{\mathrm{3},\mathrm{62}}=\mathrm{2.49}$, p=0.024). A final pairwise PERMANOVA was run to determine which sites (Table A4) and peat types (Table A5) differed specifically from one another. Among the sites, all pairwise comparisons were statistically different from one another, except for MQ and RM (R²=0.071, $F_{\mathrm{1},\mathrm{38}}=\mathrm{2.912}$, p=0.066). The strongest difference was observed between BOG70 and RM (R²=0.667, $F_{\mathrm{1},\mathrm{29}}=\mathrm{58.764}$, p=0.0015). Among peat types, all comparisons were statistically significant, except for bryophyte and herbaceous peats (R²=0.103, $F_{\mathrm{1},\mathrm{21}}=\mathrm{2.402}$, p=0.092). Ligneous and ligneous (mangrove) showed strong differentiation from bryophyte (Ligneous: R²=0.153, $F_{\mathrm{1},\mathrm{38}}=\mathrm{6.908}$, p=0.004; Ligneous (mangrove): R²=0.687, $F_{\mathrm{1},\mathrm{10}}=\mathrm{21.976}$, p=0.009) and herbaceous peats (Ligneous: R²=0.144, $F_{\mathrm{1},\mathrm{55}}=\mathrm{9.275}$, p=0.004; Ligneous (mangrove): R²=0.412, $F_{\mathrm{1},\mathrm{27}}=\mathrm{18.924}$, p=0.004). Ligneous and ligneous (mangrove) also displayed significant differences from one another (R²=0.090, $F_{\mathrm{1},\mathrm{44}}=\mathrm{4.343}$, p=0.025).

4.4.3 Peatland successional pathways in Costa Rica

Ecological succession can be inferred from the plant macrofossil results. It is noteworthy that three out of our four study sites started as palm swamps (GAN, RM, MQ). While GAN has remained a palm swamp throughout its development (800 years), the mangrove site (RM) switched from palm to mangrove within the past 500 years. While ecological succession in wetlands and peatlands tends to go from wetter to drier environments, there exist many examples of opposite trends (from drier to wetter) that can be caused by changes in local hydrology that might have been induced by a lowering of the base level (land subsidence, sometimes caused by earthquakes or large erosional storms such as hurricanes) or an increase in sea level (causing flooding and saltwater intrusion) (McCloskey and Liu, 2012; Phillips and Bustin, 1996; Urquhart, 2009; Urrego et al., 2019). The coastal area where our mangrove site is located is both geomorphologically and tectonically active, with an estuary in close proximity and a history of paleo-earthquakes (Denyer, 1998). As for the riverine site (MQ), the switch from palm swamp to herbaceous peatland took place around 5500 years ago and suggests a transition into a drier (though still wet) environment. The river that transects this peatland complex (Rio Medio Queso) is a tributary of Rio San Juan, itself a major outlet of Lago Nicaragua (Fig. 1). Mid-Holocene fluctuations in the lake level would have likely caused large shifts in the hydrological regime and botanical nature of the regional wetland complexes that span hundreds of km² in the area.

The BOG70 core (montane peatland) indicates a recent switch from herbaceous to bryophyte peat that might have taken place within the past few decades (Fig. 9). It is important to note that the modern ecological assemblage of this peatland is unusual for Costa Rica; Sphagnum and lichen species are typically found in cool environments of the mid- and high-latitude regions of the world. A transition from herbaceous-dominated fen to moss-dominated bog (or poor fen) only requires the peatland surface to become isolated from the groundwater (e.g., Hughes and Barber, 2003). It is possible that regional changes in hydroclimate such as an increase in rainfall or fog have allowed for Sphagnum peat to colonize this site where it quickly developed dense moss carpets that promote their own sustenance (van Breemen, 1995). Alternatively, we may be witnessing the autogenic development of a rain-fed or fog-fed system (dominated by Sphagnum moss) that is transitioning away from a groundwater-fed herbaceous through vertical peat accretion (Hughes, 2002). Atop the Cordillera de Talamanca, where the site is found, temperatures today are cool and much of the moisture is brought in the form of fog (cloud forest). While those conditions may be suitable for Sphagnum and lichen growth, the question about their provenance remains.

4.4.4 Peat type as a driver of organic matter quality in Costa Rica

It has been shown that tropical lowland peats typically contain high levels of aromatics and low levels of carbohydrate compounds, indicative of humified, but stable, conditions that allow for long-term carbon accumulation and storage (Hodgkins et al., 2018; Verbeke et al., 2022). Our results generally corroborate the finding that tropical peat deposits contain substantial relative abundance of aromatic (and aliphatic) components that become enriched with core depth and age. Our carbohydrate data are more difficult to interpret definitively, due to aforementioned interference from silicates in the FTIR signal.

Peat was generally humified throughout all the profiles (von Post and $\mathrm{C}/\mathrm{N}$). As for peat types, our analysis suggests greater abundance of aliphatic compounds in herbaceous peat when compared to the other types, suggesting a potential signature for herbaceous plants. In addition, the relationships between OM content, aliphatics, and aromatics were positive, suggesting that stable components become enriched as organic matter accumulates (e.g., Hodgkins et al., 2018; Leri et al., 2025; Verbeke et al., 2022). Future work will further explore the relationships between peat organic matter quality and botanical composition in more detail.

The relative abundance of carbohydrates was expected to decrease downcore (representing a progressive decay of comparatively labile organic compounds) and in turn induce relative increases in the more stable components of organic matter (aromatics and aliphatics) (Broder et al., 2012; Tfaily et al., 2014; Upton et al., 2018). We identified these trends at the montane (BOG70) and the coastal peat swamp (GAN) sites, notwithstanding interference in the carbohydrate signal from some silicates in the samples. However, in the case of the riverine (MQ) and mangrove (RM) sites, the higher clay mineral content likely precluded such trends in the carbohydrate data, as the C-O band was heavily overprinted by the Si-O signal. Future analyses may identify a way to extract reliable measures of carbohydrate content from FTIR data in peat sediments with high silicate mineral content.

5 Conclusions and future work

This study offers the first detailed analysis of peat sediments for Costa Rica. The lab-based work presents important new datasets on peat geochemical properties (bulk density, organic matter content, carbon and nitrogen content). A key finding is that bulk density values tend to be greater in Costa Rican peatlands (average = 0.24 g cm⁻³), and perhaps across tropical peatland lowlands, than what is typically reported for extra-tropical peatlands (average = 0.1 g cm⁻³), possibly due to elevated clay mineral content as revealed in the FTIR spectra. The enrichment in aromatic and aliphatic carbon shown by the FTIR data confirms the humified, stable nature of the tropical peats, both in the lowlands and at high-elevation sites. Importantly, this observation applies across our four (diverse) sites, regardless of peatland type, age, or geomorphological setting, suggesting that rapid and intensive decomposition within the uppermost portion of the peat profile can lead to long-term carbon storage; we note however that the recalcitrance of the parent plant material remains to be tested. In terms of peatland dynamics, a broad range of peatland inception ages were found across Costa Rica (∼ 500 to 11 000 cal. yr BP). Spatiotemporal patterns are thus far related to geomorphological setting, with older peats found in the montane region, intermediate and variable peat ages found along inland river floodplains, and young peats found along the coast. Averaged apparent rates of peat accumulation in the lowlands have been rapid and steady over time, with age-depth models generally displaying linear trends. Lastly, the total peat carbon store in Costa Rica is estimated at 0.1573 Gt, on the basis of our new analyses combined with a recent peatland probability map. Overall, this study provides critical, first of its kind information on the characteristics of Costa Rican peats, deepening our overall knowledge on tropical peatlands.

Appendix A

Figure A1 Stacked overlay of the FTIR spectra, Gandoca (GAN) core.

Figure A2 Stacked overlay of the FTIR spectra, Red Mangrove (RM) core.

Figure A3 Stacked overlay of the FTIR spectra, Medio Queso (MQ) core.

Figure A4 Stacked overlay of the FTIR spectra, BOG70 core.

Figure A5 Principal Coordinates Analysis (PCoA) using Bray-Curtis Distances overlaid by envfit vectors separated by core. The data in the ordination is additionally visualized by (a) depth and peat type and (b) age and peat type.

Figure A6 (a) Betadisper box plot among cores, and (b) betadisper box plot among peat types.

Table A1 Summary of the proxies and analytical resolution used to analyze the Costa Rican peat cores.

Table A2 Statistical summary for the Environmental Fit (envfit) variables. Carb + sil: carbohydrate + silicates abundance; Aro: aromatic abundance; Acid: acid abundance; Aliph: aliphatic abundance; Humif: von Post humification index; Water: water content; Bd: bulk density; Omc: organic matter content; N: nitrogen content; C: carbon content; CN: $\mathrm{Carbon}/\mathrm{Nitrogen}$ ratio.

Table A3 Statistical summary for PERMANOVA, Betadisper, and dbRDA analyses.

The significance levels are reported as follows: * = 0.05, ** = 0.01, *** = 0.001.

Table A4 Statistical summary for pairwise PERMANOVA among cores.

Table A5 Statistical summary for pairwise PERMANOVA among peat types.