Effects of environmental and spatial gradients on Quercus-dominated Mountain forest communities in the Hindu-Kush ranges of Pakistan.

Quercus-dominated forests are among the most important broad-leaved evergreen forests of the Hindu Kush ranges and are currently prone to drastic anthropogenic and climatic changes. The aim of this study was to provide basic data for the development of a regional oak forest ecosystem framework for ecological restoration and management plan development to maintain local peoples' livelihoods. Hence, we analyzed distribution patterns and environmental factors that affect regional oak forests' species composition and diversity. Ward's Agglomerative clustering divided oak-dominated forest communities into three groups: i.e., Group I, dominated by Quercus baloot had an importance value index (IVI) of 89.87 ± 4.31, Group II, dominated by Quercus dilatata had an IVI of 32.16 ± 15.01, and Group III, dominated by Quercus oblongata had an IVI of 83.14 ± 4.67, respectively. The environmental factors which vary significantly within these communities were latitude, elevation, clay content and bulk density of the soil. Wilting point, saturation point, and electrical conductivity were also considered as ecosystem structural variables. Canonical correspondence analysis (CCA) indicated that community structure was affected by various environmental factors including precipitation, slope angle, elevation, clay content, and relative humidity.

Introduction

Quercus-dominated forests in Pakistan are mainly distributed in dry temperate zones ranging from 1200 to 2000 m elevation above sea level, and form pure communities or are admixed with coniferous forests in higher altitudes (Khan et al., 2010a). Forests extend to many areas across the country and Hindu-Kush ranges, covering an area of 83,840 million hectares in Khyber Pakhtunkhwa (KP) and 174,403 million hectares in the whole of Pakistan. These forests provide livelihood facilities along with ecosystem services like regulatory, provisioning and cultural services (Alamgir, 2004, Khan and Ali, 2017). Among the oak species, Q. baloot and Q. oblongata are typical elements of dry temperate areas; Q. semecarpifolia and Q. dilatata are usually distributed in moist temperate areas of Pakistan, while Q. glauca is the least common species in the moist temperate areas. The latter species is not native to Pakistan and has been planted in hilly areas of the country (Champion et al., 1965, Ahmed et al., 2006a, Ahmed et al., 2006b, Khan et al., 2010a, Khan et al., 2010b).

The forest communities and species in transition areas are potentially affected by climate change and are thus usually the first to migrate, especially when plant growth conditions are poor. This situation is abundant throughout the wet and dry temperate areas of Pakistan, particularly the Swat Highlands (Champion et al., 1965). These areas are part of the Hindu Kush ranges and the Hindu Raj extension belonging to the Sino-Japanese phyto-geographical region which spreads between Gilgit-Baltistan to Afghanistan (Ahmed et al., 2006a, Ahmed et al., 2006b). The main variables influencing vegetation patterns are known to be climate, topography and soil conditions (Siddiqui et al., 2010, Ahmed et al., 2011). The diversity of habitats, the edaphic conditions and their location in a transition between wet and dry climates make the Swat Hindukush region an excellent location to study which elements are the driving force of vegetation patterns. Steep climatic gradients and ecosystems of high elevation expansion may provide a major buffer for climatic shifts (Chapin et al., 1998, Zhang and Zhang, 2007). However, deforestation, overgrazing and land clearance for terrace farming (Ahmed et al., 2010, Shaheen et al., 2011a, Shaheen et al., 2011b) have already damaged these forests, therefore documentation of extinct plant groups and knowledge of their environmental connections is essential. Earlier plant inventories were conducted in the study region by Stewart, 1967, Beg and Mirza, 1984, Ali, 2011, Rahman, 2012, Ilyas et al., 2012, while Siddiqui et al., 2016, Ilyas et al., 2015, Rahman et al., 2017 investigated vegetation parameters. These studies concentrated on understory vegetation and either disregarded or only occasionally examined woody vegetation. In addition, some studies have also been conducted on aspects of ethno-botany and niche modelling (Ahmed et al., 2010, Sher et al., 2010, Qasim et al., 2012, Akhtar et al., 2013).

The regional temperate woodlands supply an important natural reserve basis for the local populations and are dominated by mixed forests of oak species (Quercus spp.), which strongly vary within this forest belt between 1000 and 3500 m high sea level (MSL). The predominant forest type of Banj oak (Quercus leucotrichophora) its associated by other woody species including Rhododendron arboreum, Lyonia ovalifolia, Pyrus pashia, Alnus nepalensis, Acer cesium, Cinnamomum tamala, and Neolitsea pallens. These forests also show a great diversity of other plants including shrubs, lianas and epiphytes (Bargali et al., 2015). Due to the dense canopy, these forests produce a large quantity of litter and necromass of epiphytic mosses and lichens, affecting the soil organic matter and increasing water retention (Kumar et al., 2013). Hence, the soils have high organic matter content, high fertility, high water holding capacity and moisture content. Therefore, Banj oak woodlands contain a huge carbon pool (Gosain et al., 2015) of 268.87 t/ha (veganization + litter + ground) with a significant function for soil conservation, water preservation and a range of ecological services, regulatory and cultural supports (Bargali et al., 2015, Chakraborty et al., 2016, Dhyani et al., 2018). The forest cover in the Himalayans is influenced by changes in stand structure, pattern distribution and regeneration status (Kumar et al., 2019). About two decades ago, the regeneration of Banj oak was characterised in the Central Himalayas as scarce and insufficient (Thadani and Ashton, 1995). Recent research has shown that forests of Banj Oak have undergone changes and are rapidly deteriorating (Singh et al., 2014, Singh et al., 2016, Chakraborty et al., 2018).

Eco-climatic and edaphic factors are considered to be the dominant factors determining the pattern of vegetation, particularly in topographically complex mountain areas with steep climatic gradients (Zhang and Zhang, 2007) like the Hindu Kush Himalayan region, where many Quercus species are distributed (Ahmed et al., 2011). Besides physical factors, long-lasting anthropogenic impact including forest degradation, vegetation clearing, mountain farming, and overgrazing of the forest understory are relevant factors disturbing any forest communities in the area (Ahmed et al., 2010, Shaheen et al., 2011a, Shaheen et al., 2011b).

A preliminary study elaborating oak forests in the Swat Khyber Pakhtunkhwa (KP), Pakistan using non-numerical phytosociological techniques, provided basic information on dry oak forest communities (Beg and Mirza, 1984), but did not explore the natural requirements for forest distribution patterns, growth and dynamics. These aspects, however, are crucial to better understand the natural factors determining the distribution, ecological niche, and climatic constraints of oak forest types (Khan, 2012). So far, little is known about the extrinsic ecological and internal community factors affecting the distribution of Quercus populations in Pakistan, including the current studied area, i.e., the Swat region where these forests are highly important for the local communities in the area because of their multipurpose usage (Ahmed et al., 2010). According to previous observations (Beg and Mirza, 1984, Khan et al., 2010b), the population size of this particular taxon in the region is declining, making it urgent to identify the threatening factors to counteract Quercus forest extinction by establishing appropriate management plans. Although much research has been carried out on the impact of abiotic factors in different regions on plant diversity and species distribution, such studies in context of the oak forests needs exploration as much of the uncertainty still exists. The current study was therefore initiated i) to identify the main factors restricting the distribution and determining the composition of different Quercus populations in the study area; and ii) to identify key factors influencing the population size of forests.

Though, an inclusive investigation has not been carried out regarding the classification of forest vegetation with environmental variables, even though the area includes thickets of dry and wet forest with diverse flora. Through this study the present plant communities of Oak populations have been documented and considerable ecological and environmental levels which are responsible for their dispersal in the Swat Hindukush Mountains investigated. This study will help the growth and potential drivers of Oak plant communities to comprehend their significant ecological importance as a standard in northern Pakistan and adjacent lands.

The data obtained from this study will provide baseline information for the development of conservation strategies and for other research endeavors in the same area. Such information will also be helpful to derive the future trend lines and estimates on the stability of Quercus communities in the Swat region of the Hindukush Mountains.

Materials and methods

Study area

The present study was conducted in Swat (Fig. 1), which is bordered by Chitral, Indus Kohistan and Shangla, Bunir and FATA, and Dir in the north, east, south and west respectively (Ahmad et al., 2005). To encompass forests under different management regimes, the Oak dominated forests of both dry and moist temperate zones were sampled across the Swat Hindu-Kush range. In the temperate zone and its ecotones, oaks are intermixed with Pinus communities (Ahmed et al., 2006a, Ahmed et al., 2006b), while moist temperate Oak forest cover the lower Himalayan temperate forests of western Pakistan (Beg, 1975). All sampled plots were located in elevations ranging from 1500 to 2700 m a.s.l., and between longitude 34° 34′ to 35° 55′ N and latitude 72° 08′ to 72° 50′ E (Shinwari et al., 2003).

Fig. 1

Location of the Swat District in Khyber Pakhtunkhwa, Pakistan showing the sampling sites in black points.

Vegetation sampling

Oak forests were identified and tagged at 30 different elevations (in total 300 sampling plots) through a participatory mapping procedure (Fig. 1) that combined satellite image analysis with traditional ecological knowledge. Consequently, a stratified sampling design was applied for oak forests associated inventories and ecological studies. A total of 30 stands were quantified using quadrate method (quadrate size: 20 m × 20 m), each stand having 10 quadrates (plots). Moreover, these 300 sampling plots were set up to investigate the tree layer of 3700 oak trees, including 1800 individuals of Q. oblongata, 1110 of Q. baloot, 415 of Q. dilatata, and 369 individuals of Q. semecarpifolia. Further, within these plots, shrub and herb vegetation was recorded using 5 m × 5 m subplots in each stand. Plant species were identified using the flora of Pakistan (Nasir and Ali, 1971). During sampling, ecological and geophysical characteristics were recorded, including site elevation, slope angle and aspect degrees by using Global Positioning System (GPS) and magnetic compass. Vegetation was recorded by listing species names, numbers, coverage, plant height and diameter. From each plot, a bulk soil sample (about 1 kg) mixed from all soil layers between 0 and 30 cm depth was collected. The soil was dried and analyzed in the Arid Agriculture Institute Swat (AAIS). The soils’ physical and chemical parameters were determined following the methods described in (Bray and Kurtz, 1945, Jackson, 2005, Tabbagh et al., 2000), including pH value, clay, silt, sand, and organic matter contents, bulk density, wilting point, field capacity, electrical conductivity, total nitrogen, available nitrogen, available phosphorus, and available potassium (Table 2).

Table 2

Average values (Mean ± stand error) of the environmental variables, i.e., topographic, edaphic, and soil parameters in the three community types (vegetation groups) isolated by Ward’s agglomerative cluster analysis.

Environmental variables	Code	Group I	Group II	Group III	F-value	P-value
Latitude	Lat	35.0 ± 0.12	34.9694 ± 0.19	34.58 ± 0.050	5.57	0.009
Longitude	Long	72.28 ± 0.031	72.29 ± 0.023	72.25 ± 0.016	1.72	0.196
Elevation (meters)	Elev	2017.7 ± 205.2	2289.5 ± 240.2	1659.4 ± 64.4	4.87	0.015
Slope (degree)	Slope	41.8 ± 1.8	41.14 ± 2.55	38.66 ± 1.14	1.15	0.33
Clay (%)	Clay	16.37 ± 1.17	17.5 ± 2.2	12.11 ± 1.02	4.65	0.018
Silt (%)	Silt	24.44 ± 2.20	26 ± 3.6	29.03 ± 2.74	0.66	0.52
Sand (%)	Sand	59.18 ± 2.66	56.38 ± 3.42	58.85 ± 3.08	0.17	0.84
pH	pH	6.4 ± 0.12	6.6 ± 0.17	6.26 ± 0.13	1.25	0.30
Organic matter (%)	Om	2.2 ± 0.78	4.6 ± 0.76	2.61 ± 0.54	2.92	0.071
CaCO3 (%)	Lime	3.16 ± 0.52	4.06 ± 0.61	3.61 ± 0.38	0.66	0.52
Nitrogen (%)	N	0.12 ± 0.04	0.18 ± 0.03	0.16 ± 0.03	0.50	0.60
Phosphorus (mg/kg)	P	8.24 ± 0.63	7.44 ± 1.12	12.94 ± 4.46	0.61	0.54
Potassium (mg/kg)	K	146.25 ± 27.07	233.42 ± 37.52	186.33 ± 26.33	1.55	0.22
Wilting point (kPa)	Wp	0.114 ± 0.4	0.11 ± 0.09	0.09 ± 0.004	4.36	0.02
Field capacity (kPa)	Fc	0.22 ± 0.6	0.23 ± 0.09	0.21 ± 0.07	1.31	0.28
Bulk density (g/cm3)	Bd	1.4 ± 0.1	1.46 ± 0.02	1.52 ± 0.01	3.43	0.046
Saturation (kPa)	Sat	0.44 ± 0.5	0.44 ± 0.008	0.42 ± 0.06	3.41	0.047
Electrical conductivity (mm/hr)	Cond	12.41 ± 1.6	12.8 ± 0.03	23.46 ± 2.80	5.44	0.01
Available water (%)	Aw	0.10 ± 0.38	0.11 ± 0.005	0.11 ± 0.04	0.46	0.63

Data analysis

Relative values of frequency, density, and dominance of individual species were expressed as percentages relative to the sum of all species in a plot. The Importance Value Index (IVI) of the tree and understory vegetation were then calculated by using the equationswhere RF is relative frequency, RD is species relative density, RB refers to relative basal area and RC is relative coverage of shrubby and herbaceous species (Cottam and Curtis, 1956, Ahmed et al., 2011). The community types were named after the first and second dominant species in the group using the IVI and its distrubtion pattern along the eleavtion gradient were dipicted in Fig. 3.

Fig. 3

Distribution of different Oak-dominated forests (Groups) along the elevation gradient in Swat District, Hindukush Range Mountain of northern Pakistan.

We considered the tree's species importance values as the first matrix and twenty environmental variables (Table 3) including topographic, edaphic, and soil physiochemical properties as response variables (2nd matrix) in quantitative multivariate analyses (McCune and Mefford, 2005, Khan et al., 2017, Wahab et al., 2008). Cluster analysis was performed to classify the studied thirty oak-forest stands into different community types (major groups). In the clustering procedure, Euclidean (Pythagorean) and effective distance measure for ecological community analysis was combined with the Wards linkage method using PC-ORD software 6.0 (Orlóci, 1967). Generally, rare species may excessively exert an undue effect on the results and hence are often excluded before analysis (Ter Braak, 1987, Skinner et al., 1998, Huo et al., 2015, Khan et al., 2013). In this study, none of the species were excluded due to the reasons, 1) that all tree species have an importance values ≥ 5% at stand level, ii) the oak dominated forests are generally species poor. Besides, rare species are often retained in multivariate analyses because they may be better indicators of ecosystem stress than common species e.g., (Khan et al., 2013, Faith and Norris, 1989, Cao et al., 2001, Clarke and Green, 1988), given the assumption that some or all of these rare species may be more sensitive to the stressor(s). The data tabulation was made using MS Excel 2010 and each community was named by the first two species with highest mean importance values (Lepš and Šmilauer, 2003, Sarker et al., 2013). The data tabulation was made using MS Excel (2010) and analysis of variance (ANOVA) was performed to compare the environmental variables among different oak-forest assemblages. Environmental variables showing significant differences were further subjected to the Post-hoc Tukey HSD test to examine inter-group variations and determination of significance.

Table 3

Pearson’s correlations among different environmental variables associated with Quercus species from 30 different stands of Swat Khyber Pakhtunkhwa.

	Lat	Long	Ele	Slope	Clay	Silt	Sand	pH	OM
Lat	1
Long	−0.175	1
Ele	0.727***	0.009	1
Slope	0.466**	−0.18	0.648***	1
Clay	0.482**	−0.16	0.363*	0.326	1
Silt	−0.298	−0.11	−0.155	0.171	−0.116	1
Sand	0.043	0.181	−0.032	−0.318	−0.378*	−0.87***	1
pH	0.406*	0.083	0.372*	0.162	0.373*	−0.392*	0.183	1
OM	0.088	0.188	0.04	0.099	0.059	0.347	−0.352*	−0.021	1
Li	0.027	0.201	−0.119	−0.22	−0.079	−0.421*	0.431*	0.052	0.32
N	−0.185	0.105	−0.064	0.139	−0.129	0.114	−0.043	−0.024	0.686***
P	−0.199	0.022	−0.081	−0.155	−0.133	0.493**	−0.395*	−0.362*	−0.085
K	0.144	0.256	0.116	0.06	0.034	−0.176	0.147	0.372*	0.45**
Wp	0.485**	−0.16	0.36*	0.313	0.997***	−0.102	−0.39*	0.381*	0.049
Fp	0.217	−0.2	0.207	0.368*	0.768***	0.541**	−0.87***	0.07	0.256
Bd	−0.37*	0.183	−0.321	−0.42*	−0.91***	−0.267	0.694***	−0.215	−0.181
Sat	0.368*	−0.18	0.31*	0.416*	0.914***	0.269	−0.69***	0.21	0.181
Cond	−0.453	0.133	−0.38*	−0.43*	−0.90***	−0.076	0.509**	−0.3	−0.08
Aw	−0.25	−0.13	−0.114	0.2	0.003	0.991***	−0.92***	−0.354*	0.346
	Lime	%N	P	K	Wp	Fc	Bd	Sat	Con
Li	1
N	0.298	1
P	−0.229	−0.1***	1
K	0.393	0.489**	−0.317	1
Wp	−0.094	−0.133	−0.106	0.032	1
Fp	−0.348	−0.043	0.23	−0.094	0.78***	1
Bd	0.246	0.082	−0.047	0.014	−0.91***	−0.94***	1
Sat	−0.247	−0.085	0.049	−0.015	0.917***	0.941***	−0.99***	1
Cond	0.191	0.113	0.06	−0.058	−0.90***	−0.80***	0.949***	−0.95***	1
Aw	−0.446**	0.093	0.49**	−0.194	0.018	0.64***	−0.374*	0.376*	−0.172

Note: Lat: Latitude, Long: Longitude; Elev: Elevation; pH: pH (1:5) OM: % Org Matter; Li: % Lime; N: % Nitrogen; P:

Phosphorous (mg/kg); K: potassium (mg/kg). Confidence levels: *P < 0.05; **P < 0.01; *** P < 0.001.

We first performed DCA-ordination (Jongman et al., 1995) to decide whether a unimodal (Ter Braak, 1987) or linear (Lepš and Šmilauer, 2003) response curve should be used in ordination analysis to explore the principal patterns of compositional variation along the first three axes, and also to elucidate the relevance of oak assemblages identified by the cluster analysis. Since the DCA gradient length was greater than 4.1 for axis 1 to 3.52 for axis 2, we preferred to use RDA and CCA respectively, which may give accurate results in such cases (Jongman et al., 1995, Zuur et al., 2007). To proceed with further analysis, we then tested the intercorrelations among different environmental variables measured in this study. The results for Pearson correlation coefficients showed that some variables were strongly inter-correlated, which were thus removed from further analysis. The inclusion of such strongly inter-correlated variables may likely yield unreliable ordination results (Sarker et al., 2013). For instance, wilting point (r = 0.997), field capacity (r = 0.768), bulk density (r = 0.915), saturation (r = 0.914) and conductivity (r = 0.901) were strongly correlated with soil clay content. Available water was highly correlated with silt (r = 0.991) and sand content (r = 0.925), respectively, and strong interrelationships existed among them (Table 4). Finally, fifteen variables were retained in ordination (Table 5).

Table 4

Results of canonical correspondence analysis of species-environmental variables operating on oak-dominated forests in Swat District, Pakistan.

Total variance (“inertia”) in the species data: 3.8440
	Axis 1	Axis 2	Axis 3
Eigenvalue	0.756	0.446	0.277
species variance %age	19.7	11.6	7.2
Cumulative % variance	19.7	31.3	38.5
Pearson Correlation	0.88	0.72	0.64
Kendall (Rank) Correlation	0.52	0.60	0.40

Table 5

Correlations and biplot scores for 15 factors operating on oak-dominated forest in Swat District, Pakistan.

		Correlations*			Biplot Scores
S.No	Variables	Axis 1	Axis 2	Axis 3	Axis 1	Axis 2	Axis 3
1	Elevation	0.88	0.043	−0.241	0.344	0.022	−0.108
2	Slope	0.51	0.057	−0.614	0.219	0.028	−0.275
3	Clay	0.569	−0.29	−0.065	0.245	−0.144	−0.029
4	Silt	−0.068	0.374	−0.224	−0.029	0.186	−0.1
5	Sand	−0.215	−0.206	0.24	−0.092	−0.103	0.108
6	P.H(1:5)	0.357	−0.119	0.14	0.154	−0.059	0.063
7	Org Matt	0.364	0.225	0.266	0.157	0.112	0.119
8	Lime	0.126	0.156	0.5	0.054	0.078	0.224
9	N	0.01	0.202	0.057	0.004	0.101	0.026
10	P	−0.237	0.154	−0.017	−0.102	0.077	−0.007
11	K	0.174	0.362	0.263	0.075	0.18	0.118
12	Tpa	0.157	0.302	0.005	0.067	0.15	0.002
13	Tmi	0.146	−0.306	−0.472	0.063	−0.152	−0.211
14	Precip	−0.184	0.49	−0.338	−0.079	0.244	−0.151
15	RH	−0.529	0.049	0.203	−0.227	0.024	0.091

The multivariate ordination methods “Redundancy Analysis (RDA)” and ”Canonical Correspondence Analysis (CCA)“ were performed to elucidate the relationships between oak-forest assemblages recorded in 30 stands and fifteen environmental variables. However, results showed that the total variance explained by CCA (38.5%) was higher than for RDA (24.2%), apart from the ”arch effect“ which might be a potential problem (Zuur et al., 2007) justifying the use of CCA-ordination to be more reliable in exposing the underlying structure in the vegetation composition and establishing main links between environmental variables and oak-forest assemblages. All default settings were used for CCA, and a Monte Carlo permutation test (499 permutations) was used to test for the significance of eigenvalues of the first conical axis. Intra-set correlations from the CCAs were used to assess the importance of each environmental variable used in this study.

Results

Classification of Oak-dominated communities

Ward’s Agglomerative Clustering technique segregated a total of 18 tree species, recorded in the 30 studied oak-dominated forest stands, into three major groups at 65% information in the abundance of data retained (Fig. 2). The composition of different species recorded in each community, with associated environmental variables is presented in Table 1, Table 2, respectively. The classification procedure identified three distinct oak forest communities, which are described in detail below together with their corresponding environmental variables.

Fig. 2

Classification of 30 sample stands dominated by Quercus species into three groups (vegetation types) using Ward’s Agglomerative clustering procedure. Note: St. (Stand number); I (Group-I); II (Group-II); III (Group-III).

Table 1

Importance value index of tree species classified after applying a Ward’s agglomerative clustering procedure (x represents absence of the species in particular group).

Species Name	Species Code	Group I	Group II	Group III
		Mean ± SE	Mean ± SE	Mean ± SE
Quercus oblongata D.Don (Syn: Quercus incana Roxb.)	Qo	X	x	83.14 ± 4.67
Quercus dilatata A.Kern.	Qd	1.37 ± 0.90	32.16 ± 15.01	0.57 ± 0.57
Quercus baloot Griff	Qb	89.87 ± 4.31	7.16 ± 5.41	x
Quercus semecarpifolia Sm.	Qs	0.37 ± 0.37	41.16 ± 18.84	x
Pinus roxburghii Sarg.	Pr	2.62 ± 1.84	3 ± 3	3.14 ± 2.02
Pinus wallichiana A.B.Jacks.	Pw	X	2.66 ± 1.97	x
Olea ferruginea Wall. ex Aitch.	Of	2.25 ± 2.25	6.33 ± 6.33	6.14 ± 3.67
Diospyros lotus L.	DI	X	x	3.07 ± 1.58
Diospyros kaki L.f.	Dk	X	x	0.21 ± 0.21
Melia azedarach L.	Ma	X	x	1.78 ± 1.78
Cedrus deodara (Roxb. ex D.Don) G.Don	Cd	2.62 ± 2.62	x	1.21 ± 1.21
Taxus baccata L.	Tf	X	1 ± 1	x
Abies pindrow (Royle ex D.Don) Royle	Ap	X	0.83 ± 0.83	x
Picea smithiana (Wall.) Boiss.	Ps	X	3.83 ± 2.80	x
Pinus gerardiana Wall. ex D.Don	Pg	0.75 ± 0.75	x	x
Juglans regia L.	Jr	0.25 ± 0.25	x	x
Parrotiopsis jacquemontiana (Decne.)	Pj	X	1.33 ± 1.33	x
Aesculus indica (Wall. ex Cambess.) Hook.	Ai	X	0.5 ± 0.5	x

Group I consist of Q. baloot, Pinus roxburghii and Cedrus deodara dominated stands. The group is represented by eight stands in the study area. This community was the most diverse of Quercus-dominated forests of the Swat region and consisted of eight tree species. Q. baloot is the dominant species with an importance value of 89.87 ± 4.31, followed by P. roxburghii with importance values of 2.62 ± 1.84 and C. deodara with IVI of 2.62 ± 2.62. Other co-dominant species were Q. dilatata, Q. semecarpifolia, Olea ferruginea, Pinus gerardiana, and Juglans regia (Table 1).

Group II consists of Q. semecarpifolia and Q. dilatata. This forest type is represented by seven stands in the study area and is dominated by Q. semecarpifolia, with an average importance index value of 41.16 ± 18.84, and by Q. dilatata with an average IVI value of 32.16 ± 15.01. Co-dominant species include Q. baloot, P. roxburghii, P. wallichiana, O. ferruginea, Abies pindrow, Picea smithiana, Parrotiopsis jacquemontiana, and Aesculus indica. Group III consists of Q. oblongata and O. ferruginea. Q. oblongata is the dominant tree species and is distributed in 15 of the studied stands. Q. oblongata and Olea ferruginea show average information values of 83.14 ± 4.67 and 6.14 ± 3.67, respectively. Co-dominant species of this group include Q. dilatata, P. roxburghii, Diospyros lotus, Diospyros kaki, Melia azedarach, and C. deodara. Moreover, the analysis of variance (ANOVA) test results reveals that latitude, elevation, clay content of the soil, wilting point and electrical conductivity were the most significant variables based on their influence with respect to major groups (Table 2).

Understory vegetation of the oak forest

The understory plants associated with Oak dominated forest are presented in Supplementary data Table S1, which was structured on the basis of tree dominated communities. Group III was found to be the most diverse and is represented by 42 different understory species, which are dominated by Berberis lycium (7.82 ± 3.54) followed by lower plants i.e., Aspelenium and Adiantum venustum (6.6 ± 2.6 and 6.05 ± 3.29), respectively. Group I and II are of lower diversity and are represented by 23 and 17 species, respectively. Group II is dominated by Aspelenium, B. lycium and Indigofera heterantha (9.86 ± 5.49, 9.14 ± 4.71 and 8.93 ± 8.93), while cluster III has A. venustum and B. lycium as dominant species (9.34 ± 6.15 and 7.66 ± 4.9), indicating the same understory as recorded in the completely oak-dominated forests. The understory vegetation belongs to 28 families, of which the dominant family is Fabaceae represented by 10 species, followed by Lamiaceae and Asteraceae represented by 8 species each, while the remaining minor families are represented by either 2 or 1 species. The understory vegetation has 11% herbs, while the remaining are shrubs (Supplementary data Table S1).

Correlations with environmental variables

Table 3 shows the correlations among different environmental parameters. We found significant positive correlations of slope with latitude and elevation (r = 0.46; p < 0.01 and r = 0.64, p < 0.001 respectively). Soil texture and the physiological properties wilting point and field capacity show significant positive correlations with clay particle percentage (r = 0.99 and r = 0.76, both p < 0.001), while bulk density and electrical conductivity show significant negative correlations (r = −0.91 and r = −0.90, p < 0.001. Organic matter percentage shows significant negative correlations with sand particle percentage (r = −0.35, p < 0.05), lime content shows positive correlation with sand particle percentage (r = 0.43, p < 0.05), and nitrogen percentage shows significant positive correlations with organic matter percentage (r = 0.68, p < 0.001).

Fig. 4 displays the influence of some significant geographic factors (latitude and elevation) and environmental factors (clay content, wilting point, electrical conductivity) on different oak-dominated plant major groups that were further subjected to Post hoc analysis of variance “ANOVA” (Tukey HSD) for inter-group variations. Latitude (Fig. 4a) revealed clear differences between the three forest groups (ANOVA, F-value = 4.05; P = 0.021; df = 3). Forest sites of group III (Q. oblongata and O. ferruginea) were clearly separated from the other forest types by elevation (ANOVA, F-value = 4.25; P = 0.015, Fig. 4b). Besides, group III was also different concerning the soil parameter clay content (ANOVA, F-value = 3.84; p = 0.03; df = 3) (Fig. 4c) and wilting point (F-value = 3.71; P = 0.037, Fig. 4d). Electrical conductivity (Fig. 4e) strongly differed among all three Quercus-dominated forest communities (ANOVA, F-value = 3.97; P = 0.02; df = 3).

Fig. 4

Box plots of environmental variables that significantly affect the distribution of oak-dominated forest in the study region. Namely, these factors include latitude (a), elevation (b), clay content of the soil (c), wilting point (d), and electrical conductivity (e), respectively.

Canonical correspondence analysis

This analysis was used to analyze the relationship between species composition and the environmental variables as shown in Fig. 5 and Table 3, Table 4, Table 5.

Fig. 5

CCA-ordination plot for vegetation and environmental variables operating on oak-dominated communities in the studied area. For explanation of species codes, please refer to Table 1.

The percentage variance of the species environment relationship for the CCA axes and their Pearson‘s correlations were r = 0.88 and 0.72 for axis 1 and 2, respectively. The species-environment correlation coefficients for the two axes of CCA were r = 0.52 and 0.60. The highest coefficients with the first CCA axis were found for slope (r = 0.61; P < 0.001), clay content (r = 0.56; P < 0.001), pH (r = 0.35; P < 0.05), % organic matter r = 0.36; P < 0.05), and RH (r = 0.52; P < 0.001). The CCA analysis showed that slope, clay content, pH, elevation, precipitation, and relative humidity were the most important factors affecting oak-dominated forests in the Swat region.

Discussion

Quercus dominated forests in the Swat region were classified by using cluster analysis, which was found to be a suitable technique for classification of plant communities (Beg and Mirza, 1984). Vegetation data were classified into three communities dominated by different Quercus species, with co-dominant associated species Q. dilatata, O. ferruginea, and P. roxburghii. The results are in agreement with Jackson (2005), who analyzed forests dominated by P. wallichiana. A similar structure of vegetation on the basis of communities’ composition and environmental variables were reported as defined by Khan et al. (2010b). The current study focuses on the classification and ordination of four Quercus-dominated forest types along environmental gradients and soil variables, while Khan et al. (2010b) studied the phytosociology, structure, and soil physiochemical characteristics of Q. baloot dominated forests in Chitral district. While this study is similar concerning the use of soil variables, it differs in using multivariate analysis to identify the environmental variables that affect community structure and composition in the Swat region. A similar study conducted on Q. baloot in the upper Dir district (Khan, 2012) reported the same associated species and DCA ordination. That previously published study had only little emphasis on the community structure, while the present approach is more comprehensive, covering maximum community parameters in numerical analysis, which were lacking in earlier studies (Beg and Mirza, 1984). The forest community at highest elevations ranging from 2800 to 3000 m was dominated by Q. semecarpifolia and Q. dilatata. These results confirm the findings of (Champion et al., 1965) on the distribution of different Quercus species, revealing that Q. baloot and Q. oblongata are found in dry temperate regions, while Q. semecarpifolia and Q. dilatata occur in dry as well as in moist temperate regions. The community occupying the middle elevation belt within an altitudinal range of 1000–2200 m is dominated by mixed stands of Q. baloot and P. roxburghii. These findings correspond to results by Khan, 2012, Wahab et al., 2008, who studied the community structure of Q. baloot and pine forests in Dir district of Khyber Pakhtunkhwa and Afghanistan, respectively.

The verifiability among Quercus-dominated communities was mostly due to aspect and elevation (Coop and Givnish, 2007). While high altitudes are dominated by the Q. semecarpifolia and Q. dilatata community, the middle elevation belt is dominated by the community of Q. baloot and P. roxburghii. In the lower elevation range, a community of Q. oblongata and O. ferruginea occurs. Species richness is mainly affected by altitude and aspect (Khan, 2012). Species richness was moderate at low altitude (Q. oblongata – O. ferruginea community), low at high elevation (Q. semecarpifolia – Q. dilatata community), and highest at the middle forest belt in the Q. baloot – P. roxburghii community (Table 1). When analyzing broad-leaved evergreen forest of Monotheca buxifolia, (Khan et al., 2010a), it was revealed that communities in the middle altitude belt showed higher species richness compared to low and high elevations. The environmental factors having significant effects on the vegetation and community structure of Quercus-dominated forests in the northern Swat area of Pakistan were identified in earlier studies, i.e., (Khan et al., 2010b, Beg and Mirza, 1984, Khan, 2012). Coop and Givnish (2007) classified different forest communities of Nepal and elaborated on the environmental factors affecting these forest types for better conservation. As confirmed by the present study, the community distribution was, in many cases, determined by soil moisture, soil nutrient concentration, precipitation, and elevation (Timilsina et al., 2007, Miehe et al., 2009).

The over story of oak-dominated forests was found to be associated with 60 understory plant species, which shows that oak forests host a considerable number of associated species. Interestingly, a study conducted in 1984 reported only 36 species (Beg and Mirza, 1984), which points to an increase in accompanying species in recent decades. This increase in understory species may indicate a decrease in forest cover and area due to anthropogenic activities, which allows new understory species to occupy space in disturbed and opened forests. The oak-pine mixed forests in temperate areas hosts 92 different species (Song et al., 2009), which share strong similarities with the understory of oak-dominated forest reported in this study. The dominance of families for understory species slightly differed from that found in other studies Khan et al. (2017) reported Asteraceae as dominant family followed by Lamiaceae, while in oak forests Fabaceae are dominant, followed by Lamiaceae and Asteraceae. The plant composition reported in this study and by Beg and Mirza (1984), is identical, with herbs dominated communities, as both studies focused on oak-dominated forests.

In this study, most of the environmental factors showed insignificant correlations with plant community occurrence indicating that these factors contribute little or nothing to the vegetation structure and community dynamics of Quercus-dominated forests. A study on forest communities in the Yellow River Delta of Eastern China (McCune and Mefford, 2005) revealed that structure and diversity overlapped between evergreen broad-leaved and coniferous forests, and that 40.8% of the undetermined variation in species distribution patterns was attributed to non-environmental, i.e., biological and stochastic factors. The regional oak forests were exposed to many threats in different situations, i.e., by military operations during 2005 to 2010, by floods in 2010 and 2012, and after that by intensive cutting by local and military people for different purposes.

Conclusions

Oak-dominated forest types in the Swat region of Pakistan were classified by considering the dominant and associated species. The most important environmental variables determining their distribution, namely elevation, precipitation, relative humidity, and slope aspect were identified. Vegetation is an important transient and undoubtedly changed ecosystem component, whereas influential environmental factors, particularly soil parent material are stable at temporal scales (Abella and Covington, 2006). This continuance suggests that the current regional Oak-forest ecosystems are probably products of human impact and have existed in this form since humans settled in the area and started to use the forests. Therefore, a precise understanding of the environmental parameters for sustaining oak-dominated forest communities using different numerical methods can provide the most comprehensive site information database improving our capability to establish appropriate protection measures and restoration guidelines for conservation strategies. The oak-dominated forests in the study region contain constituent sites widely differing in stand properties, understory vegetation and disturbance regimes. Stands of Q. dilatata and Q. semecarpifolia show indications of extreme anthropogenic disturbances and overgrazing. Stands in close neighbourhood within the oak forest ecosystems, however, showed a high degree of similarity, suggesting analogous reference conditions. The present study delivers an initial oak forest ecosystem research framework for ecological restoration in this region which should be expanded, and further developed, by more detailed research focusing on i) the temporal trends in species composition, anthropogenic regimes, and stand structure conditions, ii) Quantifying specific differences between the existing site conditions and undisturbed reference conditions, iii) Recognizing specific oak- forest communities and ecosystem-specific vegetation successional sequences for restoration goals, and iv) Conduct restoration trails with replication across the different community types to understand habitat-specific response to human disturbance and on-going climate change.

Funding

The research work was conducted on self-support basis with no involvement of any funding agency or organization.

CRediT Author Contribution Statement

N. Khan made the conceptualization and supervised the research work. A. Rahman performed the data curation and formal analysis. R. Ullah writes up the manuscript and software work. A. Bräuning thoroughly reviewed and edited the manuscript for clarity. I.U Rahman provides expertise in visualization and validation of the results.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.