Deb R. Aryal · Roldan Ruiz-Corzo

Abstract Changes in forest biomass and soil organic carbon reserves have strong links to atmospheric carbon dioxide concentration. Human activities such as livestock grazing, forest fires, selective logging and firewood extraction are the common disturbances that affect the carbon dynamics of the forest ecosystems. Here, we hypothesized that such anthropogenic activities significantly reduce the carbon stocks and accumulation rates in the tropical highland forests of the Sierra Madre de Chiapas in Southern Mexico. We sampled the Pinus oocarpa Scheide dominated forests within the elevation range of 900 to 1100 m above sea level in 2010,2014 and 2017.We measured the stand structural properties and used the reported allometric equations to calculate the tree carbon stocks. Stock change approach was used to calculate carbon accumulation rates. The results showed a gradual increase in carbon storage over the 7-year period from 2010 to 2017, but the rate of increase varied significantly between the study sites.The aboveground carbon stock was 107.25 ± 11.77 Mg ha-1 for the site with lower anthropogenic intensity, compared to 74.29 ± 16.85 Mg ha-1 for the site with higher intensity. The current annual increment for the forest with lower anthropogenic intensity was 7.81 ± 0.65 Mg ha-1 a-1, compared to 3.87 ± 1.03 Mg ha-1 a-1 in the site with high anthropogenic intensity. Although at varying rates, these forests are functioning as important carbon sinks. The results on carbon accumulation rates have important implications in greenhouse gas mitigations and forest change modelling in the context of changing global climate.

Keywords Anthropogenic disturbances · Biomass ·Carbon accumulation rates · Forest carbon pools · Forest structure · Southern Mexico

Introduction

Forest ecosystems have strong control over atmospheric CO2concentration for their vital role as carbon sinks (Leet al. 2009; Ciais et al. 2013). This carbon sink removes nearly one-fourth of the annual anthropogenic CO2emissions from the atmosphere (Liu et al. 2015).During the decades of 1960 to 2000, extensive deforestation and degradation in the tropics emitted about 0.4—2.3 Pg of carbon per year, with 15—20% contribution to the total anthropogenic greenhouse emissions (Gibbs et al.2007; Houghton 2012). Recently, the large-scale deforestation has diminished and forest carbon sink has increased due to forest regrowth (Pan et al. 2011; Sitch et al. 2015). However, such estimations still have significant uncertainties related to the size and location of carbon sinks (Pan et al. 2011; Ko¨hl et al. 2015). Natural and anthropogenic disturbances can increase such uncertainties because forest ecosystems respond differently to different disturbance intensities (Bonan 2008; Kurz et al. 2008;Canadell and Raupach 2008). Along with forest type, age and climate, such disturbances determine whether forests act as carbon sinks or sources(Erb et al.2013;Aryal et al.2014; Guan et al. 2015). Some estimations show that tropical forests are still considered net carbon sources due to persistent deforestation, land use changes and other anthropogenic disturbances (Pan et al. 2011; Reich 2011).

Such disturbances make forest ecosystems undergo continuous changes in their carbon dynamics along with time (Chazdon 2003; Balvanera et al. 2005; Vargas et al.2008; Chen et al. 2016). The changes in carbon dynamics can be linked to the changes in stand structure and species composition that vary according to the type and intensity of disturbance(Aryal et al.2014;Raymond et al.2015;Mora et al. 2017). Species composition, stand structure and carbon fluxes of the forests are changed depending upon the type and intensity of the disturbance (Edwards et al.2014; Hernndez-Stefanoni et al. 2014; Aryal et al. 2015;Cohen et al. 2016). However, we still lack long-term data related to the changes in forest structure and their carbon dynamics that cover the spatial heterogeneity created by natural conditions or anthropogenic manipulations (Gunn et al.2014;2015;Bustamante et al.2016).Monitoring the temporal dynamics of forests can be used to determine ecosystem processes that drive changes in carbon fluxes and to design climate change mitigation strategies (Galindo-Jaimes et al. 2002; Birdsey and Pan 2015; Liu et al. 2015).

The research and climate change mitigation actions of the current decade have focused on reducing carbon emissions from deforestation and degradation. However,vast areas of the tropical forests are still disturbed by human activities like selective logging, firewood collection,understory fires,animal grazing and non-timber forest product extraction (Berenguer et al. 2014; Phalan et al.2016). Anthropogenic disturbances are the key drivers in determining the variation of forest carbon dynamics in southern Mexico where perpetual conflict between rural farming practices and conservation programs are common due to the lack of appropriate resource management and social inclusion strategies (Cayuela et al. 2006; Toledo-Aceves et al. 2011; Cortina-Villar et al. 2012). The significant loss of forest cover especially for pastureland establishment led to fragmentation and consequent formation of smaller forest islands in mountainous terrain during the decades of 70—90 s in this region(De Jong et al.1999;Cairns et al.2000;Kauffman et al.2003).In addition to the huge deforestation and selective logging in the earlier decades, the remaining forests of the region still suffer strong pressure from agricultural production, livestock ranching and other livelihood activities (Flamenco-Sandoval et al.2007;Garca-Barrios et al.2009;Cortina-Villar et al. 2012; De Jong 2013; Aryal et al. 2019). Firewood collection, resin extraction and occasional fires to foment understory grasses for livestock grazing are the common livelihood activities for the rural people living around the tropical montane forests of Chiapas, a southern state of Mexico (Ramrez-Marcial 2003; Cayuela et al. 2006).Better understanding regarding the effect of anthropogenic disturbances on stand structure, species composition and carbon dynamics can help local people and policy makers to develop alternative forest management strategies that maintain or enhance carbon sinks, thereby reducing net greenhouse gas emissions.

The forests that remove and store atmospheric CO2as of primary productivity are termed as carbon sinks and those forests that emit CO2from mortality, harvests and disturbances are called carbon sources (Myneni et al. 2001). In the stock change approach, the same forest stands are sampled and resampled at different times (T1 and T2) to observe the temporal differences in carbon stocks (Aryal et al. 2014). If carbon stock increases with time, these forests are removing atmospheric carbon (sinks). If the stock difference between T2 and T1 is negative, these forests are emitting carbon to the atmosphere (sources).Many studies quantify carbon stocks at one point of time,but studies on repeated measurements of carbon stocks and accumulation rates are not so common. Long term monitoring of forest carbon stocks provide a sound basis for forest change modelling in the context of changing global climate. Understanding carbon accumulation rates also helps in land use planning and forest management decisions.

The aim of this study was to evaluate the carbon stocks and accumulation rates of the tropical highland forests with different intensities of anthropogenic disturbances. We formulated the following research questions: (1) what are the sizes of carbon pools in the tropical highland forests of Southern Mexico; (2) are these forests acting as carbon sinks or sources; and (3) do anthropogenic disturbances affect the rates of carbon accumulation in these forests?

We evaluated the temporal changes in stand structure and carbon stocks in aboveground biomass, belowground biomass,ground litter,dead wood,and soil organic carbon pools in the forests with higher and lower disturbance intensities. Specifically, we tested the following two hypotheses: (1) the net carbon accumulation rates quantified by stock change approach shall be greater than zero making the forest ecosystems a significant carbon sink;and(2) the proximity of human settlement and possible disturbances negatively affect the carbon storage and accumulation rates in different reservoirs of the tropical highland forests.

Materials and methods

Study sites

The study was conducted at two sites,namely Niquidambar and Nambiyugua, within the municipality of Villaflores in the region Frailesca, Chiapas, Mexico (Fig. 1). The first site is submerged within the Sepultura biosphere reserve and the second site is within the Frailesca natural protection area, both being part of the Sierra Madre de Chiapas.The Sierra Madre de Chiapas is a largely forested mountain chain of southern Mexico consisting of five important natural protected areas: La Sepultura (167,310 ha), El Triunfo (119,177 ha) and Volcn Tacan(6378 ha) biosphere reserves,La Frailescana(181,350 ha)forest reserve and Pico de Loro El Paxtal (60,982 ha) state reserve(Schroth et al.2009).The climate is warm sub-humid with summer rains. The average annual temperature is 24 °C and the annual rainfall ranges from 1300 to 1900 mm.Irregular topographic features with glens and depressions are the common characteristics of the study sites with the elevations ranging from 900 to 1100 m asl.. Since Pinus and Quercus species dominate the canopy of the forests located above 900 m altitude in the Sierra Madre of Chiapas (Gonzlez Espinosa et al. 2007; Martnez-Meléndez et al. 2008), the dominant vegetation of the study area is Pinus oocarpa Scheide forest (INEGI 2013). The mid elevation forests of the Sierra Madre de Chiapas (such as Niquidambar site) had experienced more anthropogenic intensities from livestock grazing and firewood extraction(Cortina-Villar et al. 2012) compared to the forests at higher elevation and difficult access (such as Nambiyugua site).We interviewed the respective local residents to know the current and past anthropogenic disturbances on the forests. In addition, we observed the visible signs of disturbances like resin and firewood extraction as well as livestock grazing in Niquidambar site. Such disturbance signs were absent in Nambiyugua sites during repeated field observations. Easy access from nearby human settlements has increased such disturbances in Niquidambar site compared to Nambiyugua site.Table 1 presents the general biophysical characteristics of the sampling sites. We established eight circular carbon-monitoring plots of 400 m2,i.e.,four plots at each site with a minimum distance of 45 m from each other. Nested plot approach was used to sample the different carbon pools.

Sampling plots and measurement of carbon pools

Aboveground carbon stocks

Within each plot, we measured all the trees having diameter at breast height (DBH) ≥7.5 cm, in 2010. The plots were further divided in to quadrants where trees were tagged and measured in clockwise direction. Consecutive measurements of the trees were made in 2014 and 2017 within the same sampling plots.Tree trunks were painted at breast height with a horizontal line to avoid the measurement error during following forest inventories. We measured tree DBH, total height, and crown diameter and compared between the study years. Tree heights were measured using clinometers. We measured the crown diameters in both cardinal directions (north—south and east—west) by vertical sighting method using meter tapes and averaged to get one diameter value. The tree DBH having more than one branch at breast height (1.3 m) was converted to one tree diameter based on the sum of the basal areas of all stems (Aryal et al. 2014).

Fig. 1 Location of the study sites in Southern Mexico

Table 1 Biophysical characteristics of the study sites located in Chiapas, Mexico

We used the available allometric equations to calculate the aboveground biomass of each tree using DBH, total height and wood density (Eqs. 1—3). Wood density values were obtained from published reports (Zanne et al. 2009;Ordón~ez Daz et al.2015).For Pinus oocarpa and Quercus peduncularis, we used the species-specific allometric equations (Vargas-Larreta et al. 2017), and for the rest of the species we used a general allometric equation (Cairns et al. 2003) with the correction by species-specific wood density(Urquiza-Haas et al.2007).We used the allometric equations developed by Vargas-Larreta et al. (2017) to calculate the biomass allocation by plant parts like stems(trunk), stem bark, branches and foliage. The sum of the parts is considered as the total biomass.

Pinus oocarpa:

Quercus peduncularis:

where AGB = Aboveground biomass of the tree (kg ind.-1); DBH = Diameter at breast height (cm); H = Total height of the tree(m);ρ = Wood density of each individual tree species (g cm-3).

ρm= Mean wood density of the trees used to generate the Eq. (0.72 g cm-3). This value was obtained from the original equation (Cairns et al. 2003) where authors calculated the mean wood density of the trees they used to develop it.

Belowground carbon stocks

The belowground (root) biomass was calculated using an allometric equation (Eq. 4) developed by Cairns et al.(1997).

where BGB = Belowground biomass; AGB = Aboveground biomass.

Both AGB and BGB were converted to carbon stocks.Once the biomass values were calculated, carbon stocks were considered as 49% of the dry biomass (Vargas-Larreta et al. 2017). Both the AGB and BGB carbon stocks were then converted to Mg ha-1by using the sampling area(400 m2) and kg to Mg conversions.

Deadwood carbon stocks

Carbon stocks in the fallen woody materials were measured by the line intersect method (Van Wagner 1982). Woody materials are the fallen branches and trunks separated from its original source (not attached to the trunk) and accumulated at ground level from 0 cm to 2 m in height. The diameter of woody material was measured along the four transects of 11.28 m oriented towards each cardinal point starting from the center of each plot. The volume of dead wood was calculated using Eq. 5.

where V is volume of dead wood(m3ha-1),L is the length of the sampling line (m), and diis the diameters of dead wood at the intersection (cm).

The dry weight of the dead wood mass was calculated by applying wood density according to the state of decomposition. Wood densities of 0.52, 0.48, 0.35 and 0.23 g cm-3were used, respectively, for fresh hard, fresh soft, slightly decomposed and highly decomposed woody material(Reyes et al.1992;Gutiérrez Vzquez et al.2010).The carbon content of the dead wood biomass was calculated using a factor of 0.47 (Fonseca et al. 2011).

Litter carbon stocks

Forest floor litter mass was sampled in eight subplots of 0.30 m × 0.30 m within each carbon monitoring plot.Fresh and decomposed litter were collected separately.Litter samples were oven-dried at 70 °C for 72 h,and then weighed to obtain dry mass. The dry litter mass was converted to carbon stocks per hectare using different fractions according to the decomposition state, 0.39 for fresh litter,and 0.27 for decomposed litter (Orihuela-Belmonte et al.2013).

Soil organic carbon stocks

Eight soil samples were collected to a depth of 0—15 cm from each sampling plot using an AMS Auger with a cylinder of 5 cm internal diameter. Four of the parallel samples were used for soil bulk density and the other four for carbon concentration analysis. The bulk density (g cm-3) was calculated as a proportion of oven dry (105 °C for 72 h) weight to the volume of the soil sample. The fraction of organic carbon of each soil sample was analyzed by dichromate digestion method(Walkley and Black 1934; De Vos et al. 2007). The soil organic carbon stock(Mg ha-1) was calculated (Eq. 6) from the values of organic carbon content (%), the depth of sample and bulk density after correcting rock and root fragments of ＞2.0 mm size (Throop et al. 2012).

where SOC = Soil organic carbon stock (Mg ha-1);SOCC = Organic carbon content of the soil sample (%);BD = Bulk density (g cm-3); D = Depth of sampling (m);Frag = Fraction of the rock and root fragments.

Current annual increment calculations by stock change approach

The annual rates of carbon accumulation in different pools were calculated by stock change approach (Aryal et al.2014). Here, we used the term current annual increment(CAI) as the rate of increase on carbon pool size between different measurement years (Eq. 7).

where CAI is the current annual increment (Mg ha-1a-1),C1and C2are carbon stocks (Mg ha-1) in different measurement years, and T is the time interval between two measurements (a).

Stand structural properties

Tree DBH, height and crown diameters were measured during the process of forest inventories.We also calculated the abundance and dominance of the tree species. The abundance is the total count of individuals from each species within a specific area (converted to one hectare in this case). The relative abundance (percentage) was calculated as the proportion of the specific species abundance to the abundance of all the species registered. The dominance is the sum of the basal area per hectare for each tree species.We calculated the relative dominance(percentage)as the proportion of the basal area of the given species to the total basal area of all the species.

Data analysis

We used basic statistics such as mean, standard deviation and confidence interval to report the amounts of carbon in different reservoirs. Storage and carbon capture rate data were statistically analyzed using one-way analysis of variance (ANOVA) to evaluate the significant differences between forest inventories. A repeated-measures ANOVA was used to test significant differences between the two study sites considering the repeated measurements performed in the years of 2010, 2014 and 2017. To test the differences in carbon storage between study sites and forest inventories, we performed a factorial ANOVA. Tukey’s posthoc test (p = 0.05) was carried out to verify the significant differences in multiple comparisons. Linear and non-linear regression analyses were performed to show relationship between different stand structural properties.A non-linear regression analysis was performed to define the relationship between tree DBH and measurement error(ME). To calculate ME, we first selected ten trees of different sizes (7.5—53 cm DBH range) randomly. Three members of the field technicians measured those trees for ten times each counting 300 measurements in total. We used the mean and 95% confidence interval for each member and tree to quantify measurement error (Eq. 8).

where ME = measurement error (%); 0.95CI = 95%confidence interval; μ = mean from the repeated measurements.

Results

Forest carbon stocks and accumulation rates

The pine forests with low anthropogenic influence stored an average of 107.25 Mg ha-1carbon in AGB compared to 74.29 Mg ha-1in the forests with higher influence(Table 2). There was a 30.8% reduction in aboveground biomass carbon storage in high-disturbed sites compared to less disturbed sites. The temporal change in AGB was statistically significant in less disturbed site (F = 10.93,p = 0.004) but not in high disturbed site (F = 3.59,p = 0.07).Dead wood material and soil organic carbon did not vary between two sites demonstrating that they are less sensitive to the disturbances.The soil organic carbon stock to 15 cm depth was 52.56 ± 11.17 Mg ha-1at Nambiyugua site and 44.54 ± 5.33 Mg ha-1at Niquidambar site.Litter carbon stock varied significantly(p ＞0.05)and was 36% lower in highly disturbed site (5.89 Mg ha-1)compared to less disturbed sites (9.18 Mg ha-1).

The AGB current annual increment during 2010 to 2017 period was 7.81 Mg ha-1a-1for less disturbed site compared to 3.87 Mg ha-1a-1for more disturbed forests.The annual carbon accrual rates calculated by stock change approach differ significantly between two sites for AGB and litter carbon pools (Table 3). The annual increase in litter carbon was 0.46 Mg ha-1a-1in less disturbed site compared to 0.25 Mg ha-1a-1in more disturbed site.There was no significant increase in deadwood carbon pool in both sites.Wilcoxon’s matched pair test showed that all the CAI values except dead wood carbon were significantly higher (p = 0.05) than zero (0). We could not calculate accrual rates for soil organic carbon because this pool was analyzed only in 2017 (Table 3).

Distribution of AGB stock by tree size and plant parts

The AGB carbon storage by tree size varied with time and study site. The major part of the AGB carbon stocks was found in medium size trees (15—30 cm DBH) in 2010 measurements at both sites.But in 2014,the major portion of the carbon stocks (67.58 out of 84.75 Mg ha-1) was found stored in bigger trees(＞30 cm DBH)at the site with lower anthropogenic disturbance. The contribution of bigger trees was even higher in 2017 measurements at this site. In the site with higher anthropogenic disturbance, the major portion of the AGB carbon stocks remained in medium size trees (15—30 cm DBH) in 2014 and 2017 measurements.The total AGB carbon stock values changed from 52.58 to 107.25 Mg ha-1between 2010 and 2017 at Nambiyugua (p ＜0.05, Tukey post-hoc test) and from 47.18 to 74.29 Mg ha-1at Niquidambar site with no statistically significant difference (Table 4). Although there was no significant difference in total AGB stocks between two study sites in 2010, we found that the difference was significant in 2017 measurements (Table 4).

Biomass allocations in different parts of the plant vary when the trees grow.The major part of the AGB was found in stem wood (60—69%) followed by branches (17—21%).Foliage and stem bark each contributed about 5—11%to the AGB and were higher in smaller trees.The rate of increase in biomass accumulation along with the increase in tree DBH was different among stem wood,stem bark,branches and foliage. For every increase in DBH, the biomassincrease in stem wood and branches were higher compared to that of stem bark and foliage (Fig. 2). The patterns of biomass allocation with tree size correspond to the trees of both study sites. However, the presence of bigger trees in less disturbed sites indicate that the biomass allocation in woody stems are higher in those sites compared to highdisturbed sites.

Table 2 Carbon stocks (Mg ha-1) in aboveground biomass (AGB), belowground biomass (BGB), litter, dead wood materials (DW) and soil organic carbon to 15 cm depth SOC-15 [Mean (± standard error)] in tropical highland forests of Southern Mexico

Table 3 Current annual increment (CAI) of carbon stocks (Mg ha-1 a-1) in different reservoirs during different forest inventory periods

Table 4 Carbon stocks(mean ± standard error) in aboveground biomass (Mg ha-1) by DBH size category in the tropical highland forests of Chiapas

Fig. 2 Changes in aboveground biomass allocation (kg ind.-1) in different plant parts:stem wood,stem bark,branches and foliage with the changes in tree DBH

Changes in stand structure and composition

The average DBH of the trees from all measurements was 24.7 ± 1.6 cm (mean ± 0.95 confidence interval) in Nambiyugua compared to 16.7 ± 0.9 cm in Niquidambar site.Tree height increases with the increase in diameter but it stabilizes at certain height. The maximum height of the tree was 27.0 m from the ground level.The average height of the trees in Nambiyugua site was 16.7 ± 0.8 m while that of Niquidambar was 13.7 ± 0.5 m. Tree height prediction equations based on DBH are presented for different measurement years (2010, 2014 and 2017) for both study sites (Nambiyugua and Niquidambar). Tree height increased with the increase in DBH in a log natural model at all cases (Fig. 3). We found reasonably comparable goodness of fit (R2values ranging between 0.51 and 0.81)among forest inventories and study sites (Fig. 3).

We found positive linear relationships between tree DBH (cm) and crown diameter (m). The significant relationship with a higher goodness of fit (R2= 0.71) was found in 2017 measurements compared to the earlier measurements.This demonstrated that the relationship was clearer when the forests got nearer to maturity (Fig. 4).Different to the trend of tree height, the positive linear relationship between crown diameter and DBH showed that crown diameter continue to increase with the increase in trunk girth.

Fig.3 Allometric relationships between DBH and height of the trees in Pinus oocarpa dominated tropical highland forests of Chiapas,Mexico.a Nambiyugua forest site for 2010;b Niquidambar forest site for 2010; c Nambiyugua 2014; d Niquidambar 2014; e Nambiyugua 2017; f Niquidambar 2017

Fig.4 Relationships between crown diameter and DBH of the trees in Pinus oocarpa forests.a 2010 measurement;b 2014 measurement;c 2017 measurement; d all the measurements

The total basal area of trees increased with time in both sites. The relative basal area showed that P. oocarpa was the dominant species followed by Q.peduncularis in both sites in all measurement years.However,we found that the number of tree species with ≥7.5 cm DBH increased with time in both sites.Only two species were registered during 2010 measurements, but five species were recorded during 2014 and seven species during 2017 measurements. As a result, the contribution of the dominant species to total basal area reduced with time in both sites (Table 5). The tree density ranged from 331 to 419 trees per hectare in less disturbed sites and 538 to 769 trees per hectare in more disturbed sites. Niquidambar site showed more trees per hectare compared to Nambiyugua. P. oocarpa was the most abundant species followed by Q. peduncularis in 2010 measurements but Byrsonima crassifolia occupied the second place in 2017 measurements with the relative abundance of 12.2—16.4% (Table 5). P. oocarpa contributed more than 93% of the total AGB. Although Q.peduncularis has been the second important species, the relative contribution of B. crassifolia to total AGB increased with time in both sites (Table 5).

Discussion

Forest carbon stocks and accumulation rates

In this study,we evaluated the carbon storage capacity and the variations caused by anthropogenic disturbances in the tropical highland forests. The gradual increase in AGB carbon stocks with time in both study sites showed that these forests are in a continuous process of carbon recovery.The forests of both sites were strongly influenced by selective logging during the decade of 1970s,which was demonstrated by the absence of larger trees of more than 50 cm DBH in 2010 measurements.Forest carbon recovery has been effective since the Mexican government launched conservation strategies like the protected area programs and the reduction of deforestation and degradation actions.

Carbon stock values in our study are within the reported range of coniferous and other highland tropical forests (a list of previous studies presented in supplementary information). AGB carbon stocks in our study were higher compared to the P. oocarpa forests of Opataro, Honduras(Alberto and Elvir 2008) but was lower compared to the similar forests of Nuevo Refugio, Chiapas (159 Mg ha-1)reported by (Rodrguez-Larramendi et al. 2016). Although the anthropogenic disturbances are not reported,dense and sparse forest gradient indicate the presence of either natural of anthropogenic disturbances in Opataro sites. There was no reported disturbance in the Nuevo Regio sites of Chiapas. The AGB carbon stock in the forest site with higher anthropogenic disturbance in our study was lower compared to many other highland and humid tropical forests of Mexico (De Jong et al.1999;Hughes et al.1999;Ordón~ez et al. 2008; Aguirre-Calderón and Jiménez-Pérez 2011;lvarez and Rubio 2013). However, there are studies reporting even smaller amount carbon in aboveground biomass (Pimienta de la Torre et al. 2007; Acosta-Mireles et al. 2009; Aryal et al. 2018). Aryal et al. (2018) reported 46.71 to 65.07 Mg ha-1carbon in aboveground biomass of the tropical forests disturbed by frequent livestock grazing.The belowground biomass carbon stocks in our study were similar to that obtained by(Rodrguez-Laguna et al.2009),who report 16.51 Mg ha-1in the Cielo Biosphere Reserve,Tamaulipas, Mexico. Litter carbon stock was different between the two study sites in 2017,but not in 2010,which explain that litter production was higher in less disturbed site than in more disturbed site. The litter carbon stock values obtained in this research (5.9—9.2 Mg ha-1) were higher than that of successional semi-evergreen tropical forests (2.2—3.6 Mg ha-1) of southern Yucatan Peninsula(Aryal et al. 2015). We assume that the Pine forests in our study site produce greater amount of needles, but with a lower rate of decomposition, so that higher amount of carbon is accumulated in this reservoir compared to broadleaved tropical forests. However, litter production and decomposition experiments should be carried out to confirm this hypothesis. Carbon stocks of the fallen woody materials fluctuate particularly because of the history of disturbance such as firewood extraction, forest fires,hurricanes and variation in wood decomposition rates.The lower amount of carbon in the dead wood pool in our study sites was probably because of firewood extraction by local communities.

Tree density per hectare was different between the study sites but the AGB carbon stocks were higher in the site with lower number of trees per hectare. Niquidambar site was found to be denser compared to Nambiyugua site but the tree basal area was higher in Nambiyugua. Higher carbon stocks in Nambiyugua sites were due the presence of larger trees.Thus,the differences in tree size in terms of DBH and tree height determined the carbon stock differences between study sites.DBH growth,not the increase in tree height, principally explained the temporal increase in carbon stocks. The relative contribution of the bigger trees to AGB carbon stocks increased gradually over time.

Studies on forest carbon stocks are increasing in recent years, but most of the studies focus only the carbon stock values. Studies reporting the temporal changes in carbon stocks and accumulation rates are not so common in this type of forests.The lack of sufficient information related to the rate of carbon accumulation according to forest type,age and disturbance gradients increase uncertainty on the carbon offset capacity of the forests(Luyssaert et al.2008;De Jong et al. 2010; Nabuurs et al. 2013). The annual carbon accumulation rates in our study sites(3.9—7.8 Mg ha-1a-1) are nearly similar to that of young tropical secondary forests (4.0—7.3 Mg ha-1a-1) of Chiapas (Orihuela-Belmonte et al. 2013) and slightly higher than Pinus oocarpa forest (1.7—5.6 Mg ha-1a-1) of Santa Ana Honduras(Alberto and Elvir 2008).Nevertheless,it is important to note that both of these studies reported the mean annual increments, calculated as carbon stocks divided by the estimated age of the forest. Carbon accumulation rate in our study was calculated by the stock change approach, where we measured the stocks at certain time intervals and the net increase was referred to current annual increment(Aryal et al.2014).The stand level AGB carbon increment was 2.7 to 3.5 Mg ha-1a-1in the tropical forests of Panama including tree growth and recruitment (Chave et al. 2003). In a semi-evergreen tropical secondary forest of Yucatan Peninsula, the AGB current annual increment rates (tree growth + recruitments - mortality) varied between 2.4 to 4.7 Mg ha-1(Aryal et al. 2014). Carbon accumulation rate in more disturbed site was 50% lower compared to less disturbed site in this study. Some earlier studies also reported that carbon accumulation decreased with increasing past land use intensities (Hughes et al. 1999; Silver et al. 2000;Orihuela-Belmonte et al. 2013; Aryal et al. 2018). Our results showed that carbon accumulation rates of tropical highland forests are relatively high demonstrating the strong potential for local greenhouse gas mitigation(Silver et al. 2000; Fehse et al. 2002;varez and Rubio 2013;Poorter et al. 2016).

The variation in wood density related to the age and size of the tree can explain the higher rate of carbon accumulation when trees grow with time.The allometric equations used to quantify biomass for P. oocarpa and Q. peduncularis (Vargas-Larreta et al. 2017) do not consider the variation in wood density. Wood density analysis by tree size can help understand its contribution to the changes in biomass growth.Site specific variations of wood density in P.oocarpa forest have been reported in Chiapas(Gutiérrez Vzquez et al. 2010). They mention that the P. oocarpa forests from lower elevation showed higher wood density compared to higher elevation forests but there are other studies reporting the reverse patterns of wood density distribution(Campos et al.2007).Since our study sites are found within the elevational range of 900 to 1100 m asl.,wood density differences between two sites are considered insignificant.

Disturbance intensity and carbon accumulation

Reconciling ecosystem services and food production has been the greater challenge of this century (Phalan et al.2011;Collas et al.2017).The higher disturbance intensity,especially from livestock grazing, firewood and resin extraction explains the lower annual carbon accumulation rates in Niquidambar site. The carbon monitoring plots were located nearer to the adjacent human settlements,which increased the intensity of such disturbances.Although within the natural protected areas, the rural people living in Niquidambar site live in the midst of conservation activities and the necessity of fulfilling their basic living requirements.A recent study demonstrated that pantropical forest degradation and fragmentation reduce the carbon stocks at forest edges up to 25% compared to forest interiors (Chaplin-Kramer et al. 2015). Human interventions and physical factors like the edge related desiccations were responsible for such reduction in biomass carbon stocks. In addition to lower biomass accumulation rates, Niquidambar site showed lower carbon stocks in forest floor litter carbon pool.Nambiyugua site is located in relatively inaccessible location where the anthropogenic disturbances like cattle grazing, resin extraction, and firewood collection are not reported.However, both sites have suffered from occasional forest fires in the years of prolonged drought. The reduction in forest carbon stocks related to past disturbance is also reported in other parts of the tropics (Silver et al. 2000;Chazdon 2003; Aryal et al. 2014).

Both sites are located within the landscape matrices with abundant availability of seeds and propagules of native species facilitating the forest recovery process (Chazdon 2003). As shown in the temporal analysis of species composition(Table 5),we noted that new species appeared with time and their relative contribution to abundance and dominance increased. We attribute such changes in tree species composition and structure to the gradual change in environmental conditions favoring the growth and development of species more adapted to the warmer climates.Some authors warn that the unfavorable ecological conditions created by climate change can lead to the significant reduction of the montane forests by altering the species composition and ecological interactions (Ponce-Reyes et al.2012).Some of the non-endemic tree species growing in the forests are possibly dispersed by birds and bats from the adjacent cultivated areas or home gardens. This is because that the regional species pool and the extent of seed dispersal can affect species composition of the changing forests (Chazdon et al. 2007). Such dispersion can be attributed to the establishment of human settlements in the forest frontiers, increasing cropped areas and finally the increasing human intensities on natural ecosystems.High accumulation of pine needles on the forest floor may prevent the pine seeds from penetrating to the soil,limiting the appearance of new individuals and thus lower number of trees per hectare but with higher DBH in the Nambiyugua site. Cattle movement may facilitate seed to soil contact and germination of new individuals leading to higher tree density in the Niquidambar site. However, this effect needs further research because the major focus of this study was to differentiate carbon stocks and accumulation rates across disturbance gradient. The 2017 tree density data obtained in our study (419—769 ind. ha-1)were higher compared to that reported in Santa Ana,Honduras, where Pinus oocarpa forest tree density was 258—320 ind. ha-1(Alberto and Elvir 2008).

Sources of uncertainty

Identifying the sources and degree of uncertainty are the key aspects of carbon monitoring, reporting and verification. Such assessments help researchers find possible strategies to reduce error in their estimates of carbon stocks and change rates. The sampling plot design and measurement in this study was based on the national forest and soil inventory method to reduce the variations caused by differences in plot design and size. Since the repeated measurements of trees were made by different persons at different measurement years, we generated a measurement error curve from the repeated measurement of the different sized trees by different persons(Fig. 5).Measurement error decreased with the increase in tree DBH and was less than one percent for the trees of ≥7.5 cm DBH. Since our carbon stocks and change rates were estimated only for the trees of ≥7.5 cm DBH,the maximum error related to tree measurement was less than one percent.

Fig. 5 Measurement error curve developed from the percentage of 0.95 confidence interval obtained from the repeated measurements of trees ranging from 7.5 to 54 cm DBH

Choice of allometric equations for biomass estimation can be another source of uncertainty (Rojas-Garca et al.2015; Henry et al. 2015). Recently published species specific equations were used for P. oocarpa and Q.peduncularis, two dominant species (Vargas-Larreta et al.2017), and a general equation with wood density adjustment was used for the rest of the species (Cairns et al.2003; Aryal et al. 2014). For the sake of comparison, we fitted our data with four published allometric equations and noted that the differences in estimated AGB stocks between models increased with increasing tree size(Fig. 6). Although initially similar, the biomass estimates of Navar(2010)and Ayala(1998)equations were lower for bigger trees and fitted poorly with our data. Cairns et al.(2003) equation fitted well after wood density adjustment but the equation was developed from tropical broadleaf forest species. The Vargas-Larreta et al. (2017) equation,that we used, fitted well with the data but the estimated biomass was slightly higher in bigger trees compared to other equations. The results of increasing carbon accumulation rates with the increase in tree size in our study is consistent with Stephenson et al.(2014),who demonstrated that large, old trees do not act simply as senescent carbon stores but actively accumulate carbon at an increasing rate when tress grow. Tree size related increase in carbon accumulation rate is also supported by the metabolic scaling theory, which indicates that mass growth rate should increase continuously with tree size (Enquist et al.1999, 2007). We chose to use the later equations (Vargas-Larreta et al. 2017) because it was developed to estimate biomass by plant parts, and was developed specifically for the two most dominant species found in the forest ecosystems of our study site.

Use of the stand or species specific allometric models for biomass estimations can improve the regional and national greenhouse gas inventories as proposed by IPCC to ensure consistent and accurate estimations of carbon sinks and sources (Henry et al. 2015). Repeated measurements of stand structure and changes in carbon stocks reduce the error of temporal extrapolation. Moreover, our approach of simultaneous measurements of five carbon pools (aboveground biomass, belowground biomass,ground litter, dead wood material and soil organic carbon)reduces the error caused by the time lag in evaluating different carbon reservoirs of the forest ecosystems.

Fig. 6 Application of four published allometric equations to all the data from this study in order to quantify aboveground biomass (Mg ind.-1)

Conclusion

Pine forests of Chiapas, part of the tropical highland forests,are the important carbon reserves since most of the flat land forests of the region have already been converted to pasture and agricultural lands. This study shows how ecosystem properties like stand structure, species composition, carbon stocks and accrual rates that emerge across anthropogenic gradients change over time. Persistent human disturbances like firewood extraction, occasional animal grazing and resin extraction reduce the carbon capture capacity of these frontier forests. Carbon stocks in biomass and soil were apparently similar but the current annual increment rate was lower in the forests with higher anthropogenic disturbance.In the long term,we can expect a clear differential vegetation response to the type and intensity of human disturbances in terms of carbon dynamics. Nevertheless, the current annual increments were higher than zero in both sites. The positive current annual increment rates in both sites indicate that the forests are acting as active carbon sinks.The results can be used as important tools to the dynamic vegetation models and greenhouse gas mitigation strategies. Future research is necessary to explain how other processes like litter production, fine root turnover and organic matter decomposition affect the size and variability of soil organic carbon pool to different depths in these forests where anthropogenic disturbance intensities continue to increase.

AcknowledgementsWe thank BIOMASA A.C. and Mexico REDD + program for supporting part of the fieldwork. We are thankful to Carrie Mitchell for English revision of the manuscript.We acknowledge the constructive comments from the reviewers on the earlier version of the article.