The importance of individuals of different sizes in the population maintenance of a palm species used by the Fulni-ô Indigenous People in northeast Brazil

Study area

We conducted data collection in the municipality of Águas Belas (9°07′03″S, 37°07′06″W), located 311 km from Recife, the capital of Pernambuco state in northeastern Brazil, within the Caatinga biome. The vegetation is characteristic of the Caatinga, dominated by xerophytic, deciduous, and thorny species (Araújo, Castro & Albuquerque, 2007). The regional climate features two distinct seasons: a prolonged dry season lasting 5 to 9 months and a brief rainy season occurring between May and July (Prado, 2003). Águas Belas has a semi-arid climate (BShw’ according to the Köppen, 1948), with an average annual temperature of 23.6 °C and mean annual precipitation of 621 mm (APAC, 2024). The predominant soil types are Regosols and Leptosols (IUSS Working Group WRB, 2022). The region is characterized by granitic mountain ranges and Cretaceous-origin plateaus (Andrade-Lima, 1960). Notably, these plateaus contain high-altitude swamps locally known as “hill marshes”—a unique vegetation physiognomy within the Caatinga biome. These areas support deciduous to sub-deciduous forests transitioning to sub-evergreen forests with evergreen species, exhibiting milder environmental conditions compared to surrounding vegetation types (Rodrigues et al., 2008). It is in these specific swamp environments that the Fulni-ô harvest leaves from S. coronata palms, which serve as raw material for their traditional handicraft production.

The Fulni-ô Indigenous People

The Fulni-ô are one of seven Indigenous groups in Pernambuco state (ISA, 2007). Despite close proximity to non-indigenous settlements, they have preserved their cultural traditions. The Fulni-ô maintain their native Yathê language, the musical Toré tradition, and a secret religious ritual called “Ouricuri” practiced from September to December (Silveira, Marques & Silva, 2012). The Fulni-ô inhabit a 11,500-hectare Indigenous territory located 500 m from Águas Belas city, with two primary villages spaced 4 km apart. The main village has approximately 3,430 inhabitants, while the more rural Xixiakhlá village has about 100 residents. During their religious rituals, some community members temporarily inhabit a third village—which shares the same name as the ritual, Ouricuri—situated 6 km from the main settlement.

Economically, the Fulni-ô rely on handicraft production and sales, musical performances in other municipalities, tertiary sector employment (Campos, 2011), and land leasing (Campos et al., 2018). While agriculture, hunting and fishing occur occasionally, these activities primarily serve family subsistence. Their handicrafts incorporate wood and seeds from native Caatinga species (deciduous/xerophytic vegetation) surrounding their villages, along with leaves from S. coronata palms. These palms grow in high-altitude swamps—called “serras” by the Fulni-ô—located within or near their territory. Currently, many artisans (particularly elders) obtain palm leaves through purchase from younger community members or non-Fulni-ô suppliers. The leaves are sold weekly in the village or at Águas Belas’ local market.

The Ouricuri palm, Syagrus coronata

S. coronata occurs in the states of Pernambuco, Alagoas, Sergipe, Bahia, and northern Minas Gerais, growing in Caatinga vegetation, remnants of semi-deciduous forests, and transition zones to restinga (Atlantic Forest) and Cerrado (Lorenzi, 2010). This palm species has a solitary, erect stipe reaching 3–10 m in height and 15–25 cm in diameter, covered with leaf scars. As a monoecious species, it produces both male and female flowers on the same inflorescence. The fruits are yellowish with a brownish tomentum, measuring 2.5–3.0 cm long, with a sweet-tasting mesocarp (Lorenzi, 2010). While the palm fruits throughout the year, fruit production peaks between May and August, with ripening occurring from October to December (Lorenzi, 2010).

The Ouricuri palm is widely used by Indigenous Peoples and local communities inhabiting the Caatinga region. Both fruit pulp and seeds are edible and can be consumed raw or cooked. They are also processed into ouricuri flour for making cakes, sweets, cookies, and bread (Crepaldi, Salatino & Rios, 2004). Additionally, the seeds yield an oil of high nutritional value used in cooking and soap production. The water from immature seeds, known as ouricuri milk, serves as a key ingredient in regional dishes such as “umbuzadas” (a traditional food from northeastern Brazil) and various fish preparations (de Lima Rufino et al., 2008). The leaves are used for crafting handmade items and as animal fodder (de Lima Rufino et al., 2008; Andrade et al., 2015). Medicinally, the water from immature seeds is used to treat eye infections, mycoses, and wounds (de Lima Rufino et al., 2008). Beyond its importance to local communities, S. coronata serves as a crucial food resource for the endangered Lear’s macaw (Anodorhynchus leari Bonaparte) (Andrade et al., 2015; IUCN, 2019), a bird species endemic to the “Raso da Catarina” region and Serra Branca Environmental Protection Area in Bahia.

For our study, we classified individuals as seedlings when they were identified in subplots during the first year of data collection with one or more entire, lanceolate, and narrow eophylls, including those that appeared in subsequent time intervals (recruits). We considered young individuals as those with fully segmented leaf blades and aerial stipes (either visible or covered by leaf sheaths) but showing no signs of reproduction. Reproductive individuals were identified by their fully segmented leaf blades, aerial stipes, and presence of reproductive structures. We further categorized all individuals into three height classes: small (≤2.5 m), intermediate (2.6–5 m), and large (>5 m). This categorization was implemented solely to illustrate the size distribution of individuals referenced in our results and discussion. During fieldwork, we recorded reproductive individuals as short as 0.9 m in height. For their handicraft production, the Fulni-ô primarily harvest leaves from intermediate and large S. coronata individuals.

Legal aspects

The research project was approved by the National Commission of Ethics in Research (CAAE 24211014.0.0000.5207), the National Indian Foundation (Authorization No. 04/AAEP/PRES/2015), the National Historical and Artistic Heritage Institute (Case No. 2000.000203/2014-35), and the System of Authorization and Information in Biodiversity (Authorization No. 41944-1). These institutions are responsible for approving research involving Indigenous Peoples and local communities in Brazil. In accordance with Fulni-ô cultural protocols, the community leaders provided verbal consent after reviewing and approving the research, which was accepted by all regulatory bodies given the Indigenous context.

Identification of harvest sites and characterization of harvest regimes of S. coronata leaves

We initially identified members of the Fulni-ô Indigenous People who utilize S. coronata leaves for handicraft production using the snowball sampling technique. This method involves identifying key informants who then recommend other experts, eventually including all community specialists in the research process (Albuquerque et al., 2014). Subsequently, we conducted a participatory workshop (Campos et al., 2018) where participants marked all leaf harvest locations on a regional map and provided information about harvesting frequency at each site. From these locations, we selected six sites representing different harvest frequencies: two with low frequency (harvested approximately once monthly per harvester), two with intermediate frequency (harvested twice monthly), and two with high frequency (harvested weekly). The absence of completely undisturbed control areas in the Águas Belas region, as confirmed by workshop participants, prevented the inclusion of a no-harvest control site.

Effects of leaf harvest on population dynamics of S. coronata

At each of the six S. coronata populations (Fig. 1), we established two permanent 50 × 50 m plots (minimum 100 m apart), totaling 1.5 ha. During initial plot establishment (July 2014, T0), we measured the total height (from ground to crown apex using graduated rods or tape measures) and tagged all non-seedling individuals with numbered plates. Over three consecutive years (June 2015/T1, June 2016/T2, July 2017/T3), we monitored survival, height growth, and recruitment of new individuals.

For seedling assessment, we established five 10 × 10 m subplots per main plot (four corners and center) at T0. All subplot seedlings were tagged and measured initially, with subsequent monitoring of survival, growth, and new recruitment during annual resampling. Plants that were missing were recorded as mortalities. Additionally, we counted infructescences on reproductive individuals quarterly during the first two study years.

Additional information about S. coronata leaf harvesting regimes

We conducted semi-structured interviews (Albuquerque et al., 2014) with harvesters to obtain additional information about harvesting regimes and the characteristics that make certain individuals preferable for harvesting. Additionally, we counted the number of new leaves on each S. coronata individual every three months during the first 12 months of the study (between July 2014 and June 2015). For this procedure, we marked the newest leaf of each individual with waterproof red ink in July 2014 and again during each quarterly visit. We also recorded the total number of leaves per individual in both July 2014 and June 2015. We summed the number of new leaves produced annually by each palm, then calculated the annual number of harvested leaves for each S. coronata individual using the following formula:

Annual harvested leaves = (Total leaves July 2014−Total leaves June 2015) + Annual number of new leaves.

Characterization of the sampling sites

Luminosity, air temperature and humidity in each plot were recorded quarterly between June 2015 and June 2016 using a thermohygrometer. These environmental variables were measured at the four corners and the center of each permanent plot between 10 am and 12 pm, totaling four measures. We also sample soil at a depth of 20 cm at three points in each plot, and subsequently mixed them and sent to the Laboratory of Environmental Soil Chemistry of the Federal Rural University of Pernambuco for chemical and physical composition analysis. We assessed soil pH, Ca+, Mg+, Al+, Na+, K+, P+, H+Al, organic carbon (OC) and organic matter (OM) content. These variables were included due to differences in slopes and other landscape characteristics among the study areas, which led us to question whether they could also influence the demographic responses of the population. The characteristics of sampling sites are depicted in Table 1 and in Fig. 2.

Data analysis

We calculated the average number of leaves harvested annually from each S. coronata individual per harvest frequency by summing the annual harvested leaves of all individuals within the same harvest frequency category and dividing by the total number of individuals in that frequency category (Table S1).

We performed a principal component analysis (PCA) (Legendre & Legendre, 2012) to assess differences among sites based on environmental and anthropogenic predictors. The variables included in the PCA were: soil pH, Ca²⁺, Mg²⁺, Al³⁺, Na⁺, K⁺, P⁺, H+Al, OC, OM, air humidity, air temperature, and luminosity.

Demographic data were used to construct an integral projection model (IPM) to calculate population vital rates using continuous variables rather than size classes, as in matrix models (Easterling, Ellner & Dixon, 2000). The IPM describes changes during a discrete period in a population whose structure is characterized by a continuous size variable, usually diameter or height. The initial population size is described by a probability density function n(x, t), which represents the proportion of individuals of size x at time t. Thus, the model for calculating the proportion of individuals of size y at time t + 1 is defined as follows: ${\displaystyle n\left(y,t+1\right)={\int }_{\Omega }\left[p\left(x,y\right)+f\left(x,y\right)\right]n\left(x,t\right)dx={\int }_{\Omega }k\left(y,x\right)\left.\right]n\left(x,t\right)dx}$

where k(y, x) represents all the probabilities of transition from an individual size x at time t to an individual size y at time t + 1, including new recruits. The kernel function is integrated over a set of all size possibilities (Ω) and has two components: the survival and growth function p(x, y) and the fertility function f(x, y). The fertility function is positive for reproductive individuals at time t (parents, x) and small individuals at time t + 1 (offspring, y), obtaining the value 0 for all the other individuals. The function p(x, y) incorporates the growth and survival of all individuals (Easterling, Ellner & Dixon, 2000). The survival-growth function is composed of two components: ${\displaystyle p\left(x,y\right)=s\left(x\right)g\left(x,y\right).}$

The s(x) component represents the probability of survival of an individual of size x, and the g(x, y) component represents the probability of growth from size x to size y. The fertility function is also composed of two components: ${\displaystyle f\left(x,y\right)=f1\left(x\right)f2\left(x,y\right).}$

The component f1(x) represents the average number of offspring produced by an adult of size x, which was estimated by the ratio of the number of reproductive individuals per plot in one year to the number of seedlings recruited in the next year (Table S4). In the second component, f2(x, y) represents the probability that an adult of size x will produce an offspring of size y.

We calculated alternative statistical relationships for growth and survivorship (vital rates) as functions of plant size (independent, size, size² size³). Additionally, we subsequently applied model selection methods based on the Akaike information criterion (AIC) to determine which model provided the best fit to the data. We tested the following functions: vital rate∼1, vital rate∼size, vital rate∼size + size², and vital rate∼size+ size²+ size³. To represent the final results of the probabilities for each vital rate calculated for the kernel K, graphs were constructed for each interval and sampled population. Furthermore, we selected the model with the lowest AIC score or the simplest model, using ΔAIC ≤ 2 (Tables S2 and S3).

As some populations had very few individuals, the two populations with the same leaf harvest frequency were grouped together. In this sense, a model was built by time interval for each frequency of leaf harvest, totaling nine models (three different harvest frequencies and three time intervals). The IPM was constructed based on the total height, as all individuals (seedlings, juveniles and adults) had this measurement. The deterministic population growth rate (λ) was obtained through the kernel for each of the harvest frequency and time intervals to establish trends in population growth. Thus, λ < 1 indicates that the population is declining, λ = 1 indicates stability, and λ > 1 indicates that the population is growing (Caswell, 2000; Caswell & Fujiwara, 2004).

Prospective perturbation analysis of elasticity was performed for each sampling interval to verify which individuals in the population and which vital rates had the greatest influence on the population growth rate (De Kroon et al., 1986). ${\displaystyle e\left(y_0,x_0\right)=\frac{K\left(y_0,x_0\right)s\left(y_0,x_0\right)}{\lambda }}$

where e (y, x) is equal to the integral of 1 and can be interpreted as the proportional contribution of K (y, x) to population growth (Ellner & Rees, 2006). In the equation above, s(y, x) represents the sensitivity formula ${\displaystyle s\left(y,x\right)=\frac{v\left(y\right)w\left(x\right)}{\left(v,w\right)}}$

where 〈v, w〉; is the product 〈v, w〉 = ∫xv(x)w(x)dx (Ellner & Rees, 2006).

Confidence intervals for population growth rates with 95% confidence intervals were calculated using bootstraps, with 2,000 replicates. In each resampling, 10% of the individuals from each population and time interval were randomly removed.

Life table response experiments (LTRE) were performed to compare the population dynamics under different harvest frequencies (low vs. intermediate and low vs. high) using the following equation: ${\displaystyle \lambda t-\lambda c=\Sigma \left(Kt-Kc\right)\left(\partial \lambda /\partial |km\right)}$

where λt corresponds to the population growth rate under exploitation pressure (intermediate or high) and λc to the population growth rate without exploitation pressure (low). Kt and Kc are the kernels corresponding to the populations under exploitation pressure and without exploitation pressure, respectively. Finally, (∂λ/∂K) is the sensitivity of λ to the perturbation of the mean kernel element (km) (Merow et al., 2014).

The statistical analyses were performed in the R 4.2.3 environment (R Development Core Team, 2023) using the packages IPMpack (Metcalf et al., 2013), fields (Nychka et al., 2015), gamm4 (Wood, Scheipl & Wood, 2017) and popbio (Stubben, Milligan & Nantel, 2008).

Harvesting regimes

Based on interviews with 27 harvesters, the preferred individuals for leaf harvesting were those with the largest and widest leaves. Harvesters typically remove all leaves except for two or three of the youngest leaves. Respondents reported an average interval between harvesting events on the same palm ranging from at least 6 months up to a maximum of 3 years. The mean number of leaves harvested annually per individual was 2.8 in low-frequency harvest sites, 3.61 in intermediate-frequency sites, and 5.14 in high-frequency sites (Table S1).

Population structure

A total of 549 S. coronata individuals were sampled between 2014 and 2017, with an average density of 366 individuals per hectare. Population abundance varied substantially among study sites. The highest abundance occurred in intermediate-harvest-frequency populations (n = 273, 546/ha), followed by high-frequency (n = 190, 380/ha) and low-frequency (n = 86, 172/ha) populations. The tallest individuals were recorded in intermediate-harvest-frequency populations, followed by high and low-frequency populations (Fig. 3). Populations with intermediate harvest frequency showed height class distributions closest to a reverse-J shape in nearly all censuses except the final one (Fig. 3).

Characterization of leaf harvest sites

The first two principal component axes (eigenvalues = 5.397 and 2.978) collectively explained 86% of the environmental variation in the study area (Fig. 4). Temperature, soil pH, and H+Al showed the strongest loadings on principal component 1 (PC1), while potassium was most influential on principal component 2 (PC2) (Fig. 4). These variables were highlighted as they exhibited eigenvector values ≥0.70. Luminosity and Mg⁺ were excluded from interpretation due to eigenvector values <0.6 on both PC axes.

The two sites with low-harvest-frequency populations showed substantial environmental and soil differentiation. Soil pH and air temperature were particularly important in sites with the lowest harvesting frequencies, whereas soil nutrients played a greater role in sites experiencing intermediate and high harvest frequencies.

Influence of leaf harvest on the population dynamics of Syagrus coronata

All populations decreased in the last year of harvest but exhibited different trends between the first and third time intervals. The populations subjected to low leaf-harvesting frequency grew at the second sampling, decreasing sharply at the third one (Table 2). Populations with intermediate frequency of leaf harvest decreased in all sampling intervals, while highly harvested populations increased in the second time interval, then decreased in the next one (Table 2). Moreover, they had the lowest population growth rates (λ).

The IPM analyses revealed that individuals from populations subjected to low leaf harvest frequency exhibited equal growth probabilities during the first and second time intervals. In the third time interval, however, intermediate and large-sized individuals showed greater growth compared to small individuals (Fig. 5). Small individuals from populations subjected to intermediate leaf-harvesting frequency showed greater height growth in the first and second time intervals, but experienced higher mortality rates in the second and third intervals. In these populations, adult size did not change significantly. Among the populations subjected to the highest harvest frequencies, the younger individuals experienced greater height growth during the first time interval, but their mortality rate was also high. Taller individuals had a lower height growth rate and a lower mortality rate at all time intervals (Fig. 5).

The elasticity analyses showed that, in populations subjected to lower harvest frequencies, the growth of smaller individuals had a greater influence on the population growth rate (Fig. 6). Conversely, in populations subjected to intermediate and high harvest frequencies, the presence of larger individuals had a greater influence on the population growth rate (Fig. 6).

The LTRE analyses showed that the stasis of smaller individuals contributed positively to the increase in the population growth rate in the populations subjected to intermediate harvest frequencies. However, in these same areas, the increase in height of individuals of intermediate size contributed negatively to the population growth rate. Among the populations subjected to the highest harvest frequencies, the presence of larger individuals contributed positively to the first time interval, while the growth of small and intermediate individuals contributed negatively. In the second interval, the stasis of smaller individuals and the growth of intermediate individuals negatively contributed to the population growth rate. In the third interval, the growth and regression of the intermediate individuals negatively contributed to the population growth rate (Fig. 7).

Limitations of the study

We acknowledge that the lack of a control area (with no harvesting) limited our conclusions regarding the factors that drive the demographic response of S. coronata populations. Without a control area, it is difficult to determine whether harvesting alone or environmental factors are responsible for the observed population patterns. Despite this limitation, studying multiple sites for each of the three harvesting frequencies helped us better understand the influence of this factor on the species’ demography. Another factor that somewhat limited our conclusions was the small size of the studied populations, though we believe this was mitigated by conducting data collection over four years (across three time intervals), which provided more IPMs for interpretation.

We also recognize that the lack of information about past land management conditions, such as grazing and burning practices, makes it difficult to draw definitive conclusions about the factors that best explain the current demographic characteristics of the populations. We recommend that future studies incorporate historical land-use data to clarify how these practices may have shaped population structure over time.

While we cannot state with complete certainty which factors are primarily responsible for the demographic responses of S. coronata, our study makes important contributions by demonstrating that maintaining individuals of different sizes is crucial for population viability. Furthermore, this study helped identify the most suitable areas for leaf collection by the Fulni-ô Indigenous People, supporting both their extractive activities and handicraft production—practices that represent traditional knowledge and an essential way of life. By including environmental variables in our study, we have shown that these factors are important drivers and should be considered in studies examining harvesting effects on plant demography.

To address these limitations, we recommend that future studies in this field include control areas in their datasets. Another suggestion is to implement long-term population monitoring.

Conclusion

Our study showed the importance of individuals of different sizes for maintaining S. coronata populations, whose leaves are used by the Fulni-ô Indigenous People in northeast Brazil. In populations with lower harvest frequencies, the growth of smaller individuals most influenced population growth rates, while in populations with intermediate and high harvest frequencies, larger individuals were most important for population growth. Despite exhibiting the highest population growth rates, populations subjected to low harvesting frequencies showed limited seedling recruitment and concerning demographic structure. Based on our results, we recommend the following management practices for S. coronata populations in collaboration with the Fulni-ô Indigenous People:

(i) Maintain current leaf harvest in areas with low harvest frequencies, as these populations showed high growth rates.

(ii) Protect smaller individuals (under 2.5 m) in low harvest frequency areas. Although their growth contributed significantly to population growth rates, the low recruitment observed highlights the need to protect seedlings and understand what limits new individual establishment in these areas.

(iii) Preserve intermediate-sized (2.6–5 m) and larger individuals (over 5 m) in areas with intermediate and high harvest frequencies. The survival of these size classes, which produce better quality leaves for Fulni-ô handicrafts, was most important for population growth in these declining populations. We suggest alternating harvest years in these areas (e.g., harvesting every other year).

(iv) Develop agreements between the Fulni-ô and tenant farmers on Indigenous lands. Farming practices can make environments less suitable for palm establishment. Vegetation clearing and livestock trampling may have contributed to low recruitment in Site 2 (Cigana), where harvest frequency is lower.