Does pruning affect the structural and ecological productivity of Juniper woodlands in the eastern Hindu Kush?

Study area

Chitral, previously the largest district of Khyber Pakhtunkhwa with a total area of 14,850 km², later underwent administrative division resulting in Lower and Upper Chitral districts. Chitral geographically extends between latitudes 35°15′06″ to 36°55′32″N and longitudes 71°11′32″ to 73°51′34 E, in the northern part of the province (Fig. 1). The area lies at elevation ranges of 1,100 m to 7,708 m above sea level (Khan, Ahmed & Shaukat, 2013). The climate data (2016–21) show that Chitral is characterized by arid conditions, with the majority of rainfall occurring between January and April, accounting for 73% of the total rainfall and representing the wettest months of the year. In contrast, the period from June to September experiences the lowest rainfall, accounting for only 3.68% of the total, making these months the driest of the year. June, July, and August are also the hottest months, with temperatures ranging from 34–36 °C, while December, January, and February are the coldest months in the area. The relative humidity (RH) maintains an average of 72.1% at 0300Z, with highest recorded RH in August, September, and October at 80.5, 83.7, and 84.8%, respectively (Fig. 2).

The forest of the area is dominated by conifers, particularly Cedrus deodara, Juniperus excelsa, Pinus geradiana, Pinus wallichiana, Abies pindrow, and Picea smithiana (Khan et al., 2022), followed by broadleaved species like Quercus incana, Quercus dilatata, Quercus baloot, Juglans regia, Tamarix dioica, and Betula utilis (Khan, Ahmed & Shaukat, 2013). Our present study focuses on the Mastuj Valley in Upper Chitral, located at approximately 36.3°N latitude and 72.5°E longitude. Geographically, it is bordered by Gilgit to the east, Afghanistan to the north, and District Swat to the east and south (Fig. 1). In the studied valley, Juniperus semiglobosa, Juniperus excelsa, Salix alba, Salix nigra, and scattered birch patches dominate the wild habitats.

Field data acquisition

Data were collected in August 2022 from three representative juniper forest stands located in close proximity based on pruning intensity as the key criterion. The first sampling site represented forested land exhibiting a moderate level of pruning. Inhabitants in this area predominantly engaged in mixed forest management practices, particularly to facilitate grazing activities by utilizing the forest floor. Our second sampling site was situated alongside the agricultural and residential lands and adjacent to orchards. Notably, pruning practices in this setting were more pronounced, driven by the necessity for sunlight exposure essential for cereal crops and understory vegetation. In contrast, the third sampling site, characterized by minimal anthropogenic disturbance, resulting in limited pruning activities, therefore, served as the control group (proxy) (Fig. 3). Vegetation and soil attributes in pruned and non-pruned sites were sampled using 20 × 20 m (400 m²) plots for tree stratum and 5 × 5 m subplots for understory vegetation along 200 m line transects across each stand (Tesfaye, Gardi & Blaser, 2019). The entire plants both in the overstory and understory strata were identified with the help of the flora of Pakistan (Stewart, Nasir & Ali, 1972), pictorial guides, and were confirmed from Plantworld.com (https://plantworld.com/; Kew Botanical Garden). Furthermore, the density (individuals ha⁻¹) of both trees and understory species was calculated, and the Shannon–Weiner diversity (H’) index was determined (Wu et al., 2015).

Key dendrometric parameters, including tree height, DBH, DRC, and crown base height (CBH) were measured using standard forestry methods (Bezos et al., 2012). During the sampling phase, these three sites were classified into three categories; non-pruned stand (hereafter NPS), moderately pruned stand (30–40%; hereafter MPS) and intensive pruned stand (≥ 60%; hereafter IPS), based on the ratio of the crown base height (CBH) to the total tree height (TH) (Huang et al., 2023). For categorization of the pruning intensities we used the given equation CBH/TH × 100, for details see in Fig. 4.

For growth rate assessment, more than 40 healthy trees of the targeted species in each forest stand were bored at breast height in damage-free parts through Swedish increment borer (five mm diameter) following Fritts (1976) and Speer (2010). The extracted core samples were then stored in plastic straws along with the site description, following standard methodology (Kahveci, Alan & Köse, 2018). In the laboratory, wood samples were air-dried and subsequently sanded with different girth-size sandpapers until the visibility of the growth rings (Fritts, 2001). To assess age and growth rates, annual growth rings were carefully counted with the help of magnifying lenses and a binocular microscope, following Ahmed et al. (2011). To evaluate reproductive efforts, the numbers of female cones were randomly collected from a total of five trees at 4–5 branches at each stand. Likewise, we selected 100 subsamples of cones for counting the filled seeds (living embryos) by cutting it transversely (Ortiz, Montserrat & Salvador, 2002). Generally, 3–5 and rarely six seeds are present per female cone of J. excelsa and 2–3 seeds in J. semiglobosa (Din, 2018).

Environmental and soil variables measurements

The geospatial coordinates including latitude, longitude, and elevation of each site were recorded using a geographical positioning system (GPS). Soil samples (weighing one kg each) were collected in triplicate within each stand with one sample taken from each of the two corners and the center, at various depths (ranging from 0–15 cm and 15–30 cm), followed by Weil & Brady (2017). These samples were combined to create a composite sample which was then air-dried and passed through a two mm sieve to remove any extraneous materials (Ullah et al., 2022). Finally we analysed the soil samples for various parameters, such as CaCO₃ percentage, soil texture properties (clay, sand, and silt %), pH, and organic matter (OM %), following Bartels, Bigham & eds (1996) and Dane, Topp & eds (2002). Additionally, electrical conductivity (EC in dS m⁻¹), total phosphorus (mg kg⁻¹), and potassium content (mg kg⁻¹) were also analysed in the Agricultural Research Station (ARI) Chitral, Khyber Pakhtunkhwa. By using clay and sand percentages of soil, we further calculated soil hydraulic properties (SHP) such as conductivity (mm hr⁻¹), saturation (%), bulk density (g cm⁻³), wilting point (kPa), water availability (%), and field capacity (kPa), through online calculator (http://www.dynsystem.com/netstorm/soilwater.html) following Saxton et al. (1986).

Calculation of biomass and carbon stocks

We employed a previously established allometric equation (Eq. (1)), to calculate trees biomass and carbon stock (Chave et al., 2014). Predictors such as DBH, height, and species-species wood density (obtained from the Pakistan Forest Institute, Peshawar) were used for biomass calculations (Gebeyehu et al., 2019). (1)${\displaystyle \mathrm{AGB}\left(\mathrm{Kg}\right)=0.0673{\left(p{\mathrm{D}}^2\mathrm{H}\right)}^{0.976}.}$ In the equation, AGB represents aboveground biomass, p is the wood density of the species (g cm⁻³), D is the DBH (cm), and H is the height (m). Carbon estimation was performed by assuming 50% of the total biomass following IPCC (2006) guidelines, whereas below biomass estimation was derived using a simplified root-to-shoot ratio, i.e., 1:5 (20%) of AGB (Pearson, 2007).

Soil carbon stocks

Soil organic carbon was determined using soil depth (cm), bulk density (g cm⁻³), and SOC %, by using Eq. (2) (Pearson, 2007). (2)${\displaystyle \text{SOC stocks}=\mathrm{BD}\times \mathrm{d}\times \mathrm{SOC}\text{\%}\times \mathrm{CF}}$ where SOC represents soil organic carbon, BD is the bulk density (g cm⁻³), d is soil depth (cm), and CF is the conversion factor for rock fragments. Since CF was not applicable in our case due to the absence of significant rock fragments, therefore, we assumed CF = 1 for its omission

Calculation of CO₂ sequestration

CO₂ sequestration was calculated by multiplying the conversion factor (3.67), which represents the ratio of the atomic weights of carbon (12) to carbon dioxide (44), accounting for the additional oxygen atoms in CO₂, with TLC (AGBCS, BBCS) using Eq. (3) (Afzal & Akhtar, 2013; Toochi, 2018). (3)${\displaystyle {\mathrm{CO}}_2\left({\mathrm{ton CO}}_2\right)=3.67\times \mathrm{TLC}}$ where CO₂ is CO₂ sequestration and TLC is total living carbon (Mg C ha⁻¹, where 1 Mg = 10⁶ g).

Statistical analysis

Descriptive statistics such as mean and standard error (mean ± SE) were used to illustrate the mean differences among these three forest stands (Fadairo et al., 2023). Differences in the means among the sites were evaluated using a one-way analysis of variance (ANOVA) at a 0.05 significance level (Chaiya et al., 2023), whereas to elucidate the linear associations between distinct pruning treatments and selected variables, we used Pearson correlation (r) by following Bezos et al. (2012). Data normality (ANOVA) was assessed using the Shapiro–Wilk test, with a significance level of 0.05, if the p-value from the Shapiro–Wilk test was greater than 0.05, the data were considered normally distributed (Shapiro & Wilk, 1965). For statistically significant variables, pair wise comparisons were conducted using Tukey’s HSD test (Pearson, 2007). A 95% confidence interval for the correlation coefficients was computed using Fisher’s z- transformation, where the correlation coefficients were transformed to z-scores, and confidence intervals were constructed in the z-scale before back-transforming to the original correlation scale (Steur et al., 2020). All statistical analyses were performed using IBM SPSS Statistics version 30.0.0.0 (IBM Corp., Armonk, NY, USA) and OriginPro version 2024b software (OriginLab Corporation, Northampton, MA, USA).

Forest characteristics and physical environment

The forest comprises six native tree species from four families, including two species of Juniperus (J. semiglobosa and J. excelsa) in the overall overstory composition. The sites are geographically located at elevations between 2,607 and 2,668 m (a.s.l.), facing northwest (126–153°) and southwest (260°) aspects (Table 1). J. semiglobosa is the dominant species in all forest stands, with importance values of 31–36%, followed by J. excelsa (IV = 23–29%), except in the intensively pruned site (IPS), where Elaeagnus angustifolia co-dominates with 18.5% IV. Salix alba contributes 16.45% in the moderately pruned site (MPS), while two other species contribute less than 5%. The highest IV for Salix alba (15.21%) was observed in the IPS, with Hippophae rhamnoides having the lowest IV across all sites. Populus alba, considered a minor species, ranged from 3 to 8.6% IV in the MPS and IPS but was absent in the non-pruned site (NPS). Tree density was highest in the MPS (301 stems ha⁻¹), while basal area was greatest in the NPS (26.84 m² ha⁻¹). The IPS had the lowest density (263 stems ha⁻¹) and basal area (11.36 m² ha⁻¹) (Table S1). Generally, these forest stands had non-saline soils (EC = 0.13 to 0.99 dS m⁻¹) with slightly acidic (IPS = 6.84) to alkaline nature (NPS & MPS ≤ 7.21). Organic matter ranged from 0.31 to 0.98% in the three forests, whereas CaCO₃ was high in MPS followed by NPS. Additionally, soil texture was predominantly sandy loam in all stands, with high K⁺ content (133–171 mg kg⁻¹). Soil phosphorus was generally low, ranging from 10.2 to 12.53 mg kg⁻¹. Hydraulic properties such as saturation, bulk density, wilting point, water availability, and field capacity also varied across the three stands but did not differ significantly (Table 1).

Impact of pruning on dendrometric and reproductive traits

In our study, we analyzed the effects of pruning on both Juniperus semiglobosa and Juniperus excelsa to identify the potential species-specific responses. However, there was a narrow margin in variation, and statistically non-significant differences between these two species in their response to pruning intensities were found. As a result, we collectively analyzed the data of both Juniperus species (junipers) for more concise and generalized findings. The density of juniper trees ranged from 137 to 156 (individuals ha⁻¹) with a basal area of 5.91 to 21.44 (m² ha⁻¹) and did not differ significantly between the experimental and control sites (Table 2). Similarly, the mean DBH and DRC ranged from 17 ± 0.75 to 20 ± 0.76 cm, and 23 ± 0.8 to 25 ± 1 cm respectively, with the highest values recorded at the moderately pruned site (MPS). Among the key structural attributes, the average tree volume (0.101 ± 0.012 m³) was higher in moderately pruned sites (MPS) compared to intensively pruned sites (0.08 ± 0.008 m³); however, the difference was not statistically significant across the sites. The youngest recorded juniper tree was 17 years old, and the oldest was 55 years old at both moderately and intensively pruned sites. In contrast, older trees, ranging up to 169 years with an average age of 66 years were found in the non-pruned sites. Subsequently, a statistically significant difference in tree age (f = 98; p = 1.50E-27) was observed among the three experimental sites. Dendrometric parameters, including height and growth rate, along with reproductive traits, showed significant differences (p ≤ 0.001) among the sites (Table 2). Specifically, the moderately pruned site had the tallest trees (µ= 5.98 ± 0.26 m), followed by the intensively pruned site, while the non-pruned sites contained shorter trees (µ= 4.8 ± 0.24 m). In contrast, seed viability was higher (2.1 ± 0.2 seeds per cone) in both moderately and intensively pruned sites, although cone production varied significantly (f = 7.96; p = 0.006), being negatively affected by pruning, with the highest number of cones recorded in the non-pruned site (100 ± 4 cones tree⁻¹). Additionally, our results revealed that both moderately and intensively pruned junipers had higher growth increments (0.53 ± 0.02 and 0.34 ± 0.01 cm year⁻¹ respectively), compared to non-pruned site (Table 2).

Impact of pruning on understory species composition and diversity

Pruning significantly impacts understory vegetation composition, density, and diversity. A total of 12 species were recorded in the understory stratum, predominated by shrubs (nine species) and a few herbs (three species) from 10 different families. The species composition in moderately and intensively pruned sites was largely similar, except for certain species, such as Rosa webbiana, Berberis lycium, Sophora mollis, Lepidium sativum, and Taraxacum officinale. Both the moderately pruned site (MPS) and intensively pruned site (IPS) hosted the highest number of species (nine species each). Understory densities across the three sites ranged from 5,411 to 7,714 individuals ha⁻¹, with the highest density in the MPS and the lowest in the NPS. Compared to the NPS, moderate and intensive pruning increased understory vegetation density by 11.45% and 7.77%, respectively. Taraxacum officinale, a member of the Asteraceae family, was the most abundant herbaceous species, while Berberis lycium and Rosa webbiana were the most prominent shrubs across all sites (Table 3). Pruning also affected species diversity, with the Shannon–Weiner diversity index (H’) being highest in the IPS (H’ = 2.07), followed by the MPS (H’ = 1.93). Species richness was similar in both pruned sites (nine species) but significantly lower in the NPS (five species). Species evenness (J’) was comparable between the MPS (0.88) and NPS (0.89), but highest in the IPS (0.94) (Table 4).

Impact of pruning on tree carbon biomass

Silvicultural treatments and their intensities had a significant impact on carbon biomass in juniper woodlands (Table 5). Moderate pruning proved to be more effective, increasing total living carbon (TLC) storage (4.95 Mg ha⁻¹), while higher pruning intensities led to a reduction in TLC stock (3.1 Mg ha⁻¹). Additionally, soil organic carbon (SOC) stock was adversely affected by more intense pruning. The highest SOC stock (nine Mg ha⁻¹) was observed in moderately pruned stands, whereas intensive pruning resulted in the lowest SOC stock (four Mg ha⁻¹). As a result, moderate pruning demonstrated the highest carbon sequestration potential (18.17 tons of CO₂ sequestered), while intensive pruning significantly reduced this capability to 11.38 tons of CO₂ sequestered.

Inter-correlation among dendrometric and reproductive characters

The linear relationships (Pearson, at 95% CI) among various dendrometric and reproductive traits under different silvicultural treatments are illustrated in Figs. 5 to 7. In addition, summary of Pearson correlation among the variables are presented in Table 6. In the non-pruned stands (NPS), tree DBH exhibited a strong positive correlation with DRC, and volume (r > 0.92, CI [0.85–0.91]), and moderate correlations with height (r = 0.40, CI [0.16–0.60]), age (r = 0.33, CI [0.07–0.54]), and tree growth rate (r = 0.46, CI [0.22, 0.46]). Similarly, DRC followed a comparable correlation pattern, except for height (r = 0.53, CI [0.31–0.70]). However, no significant correlation was found between reproductive traits and dendrometric variables. Height was significantly correlated with age (r = 0.47, CI [0.24–0.65]) and volume (r = 0.63, CI [0.44–0.77]), suggesting that age may serve as a good predictor for both these traits. In the moderately pruned stand (MPS), the pruning effects became evident as DBH showed slightly lower correlations with DRC (r = 0.80, CI [0.67–0.74]) and volume (r = 0.88, CI [0.71–0.74]) compared to NPS, but a higher correlation with height (r = 0.68, CI [0.54–0.59]), indicating enhanced height increment due to moderate pruning. Interestingly, in MPS, DBH had a non-significant weak correlation with growth rate, while DRC maintained a similar relationship with height. A notable shift was observed in growth rate, which was negatively impacted by age (r =  − 0.53 CI [−0.48–0.15]), a pattern absent in NPS. In the intensely pruned stand (IPS), DBH correlations weakened further with DRC (r = 0.57, CI [0.52–0.58]) and height (r = 0.48, CI [0.45–0.53]), reflecting the impact of intense pruning on these parameters. The reproductive traits, including seed and cone production, displayed a similar correlation trend as in NPS and MPS. Moreover, the growth rate in IPS was negatively correlated with age (r =  − 0.30, CI [−0.45 to −0.33]) and height (r =  − 0.23, CI [−0.23 to −0.02]), mirroring the pattern in MPS. The negative correlation between age and growth rate became more pronounced, while age and height maintained a positive association.