ECOLOGICAL PATTERNS OF QUERCUS SPECIES ALONG ENVIRONMENTAL GRADIENTS IN THE WESTERN HIMALAYA, PAKISTAN

INTRODUCTION

The Himalaya is globally recognized as a biodiversity hotspot, supporting exceptional species richness and endemism driven by complex topography and diverse microclimatic conditions (Tali et al., 2019). However, due to anthropogenic stress, rapid land-use change, and climate impact, the Himalayan biodiversity is being lost at a faster rate (Yadav et al., 2021; Premlata et al., 2024). Excessive use of medicinal plants further contributes to this crisis, threating several endemic species, leading to their extinction (Shivali and Namrata, 2022; Shinwari et al., 2024). These impacts are particularly severe in mountainous regions, where species adapted to narrow ecological niches exhibit limited resilience to environmental change  (Agnihotri et al., 2017).The genus Quercus (oak), comprising approximately 600 species, is ecologically and economically important across the Northern Hemisphere (Aldrich and Cavender. Bares, 2011). Oaks take up various habitats, maintain biodiversity and soil stability, and provide food and shelter for wildlife. Their durable timber, acorns, and tannins have long been valued, and their longevity and versatility make them essential components of forest ecosystems. The current phenomenon in molecular genetics and breeding provides the chance of improving their management and conservation in the context of evolving environmental conditions (Tantray et al., 2017).

In Pakistan, Quercus-dominated forests are primarily distributed in dry and moist temperate zones along elevational gradients ranging from approximately 1200 to 3500 meters above sea level (MASL), forming pure and mixed stands with other forest types (Khan et al., 2010; Rahman et al., 2022). Species such as Q. baloot and Q. oblongata dominate dry temperate regions, whereas Q. semecarpifolia and Q. dilatata are more prevalent in moist temperate zones, reflecting their ecological adaptations to varying climatic conditions (Khan et al., 2010; Rahman et al., 2022).

Vegetation distribution in mountainous ecosystems is strongly regulated by environmental gradients, particularly climate, topography, and soil properties. The Hazara Division of the western Himalaya represents a critical region for studying these dynamics due to its steep elevational gradients and heterogeneous edaphic conditions (Haq et al., 2024, 2025). However, increasing anthropogenic pressures, including deforestation, overgrazing, and land-use conversion, are significantly altering forest structure and ecological processes in this region (Singh et al., 2019; Kumar et al., 2022).

Despite the ecological importance of Quercus forests, previous studies in the region have largely focused on understory vegetation, with limited emphasis on woody species and the integrated role of environmental drivers in shaping community structure (Farooq, 2019; Ahmed et al., 2021; Faisal et al., 2022). This lack of comprehensive, gradient-based analysis limits our understanding of species distribution patterns, regeneration dynamics, and ecosystem functioning in these forests. Therefore, the current study was designed to quantify the distribution and community structure of Quercus species along environmental gradients in the Hazara Division of the Western Himalaya, and identify key environmental drivers (climatic, edaphic, and topographic) influencing species composition, population structure, and regeneration potential. By integrating vegetation and environmental data, this study provides a robust ecological framework to support conservation planning and sustainable management of oak forests under increasing climatic and anthropogenic pressures.

MATERIALS AND METHODS

Study area: Hazara Division is an ecologically diverse landscape located between the 33º-44′ to 35º-35′ North latitude and from 72º-33′ to 74º-05′ East longitude and spans over Abbottabad, Haripur, Tor Ghar, Mansehra, Battagram, and Kohistan districts as detailed in Figure 1. The study was carried out in 98 different localities, comprised of 294 quadrats and 98 stands. The area is typified by high altitude meadows, fertile cultivated land, mountainous mountain ranges, snowy glaciers, picturesque tributary systems, such as rivers, streams, lakes, and ponds. Due to diverse topography, geological and habitat heterogeneity, Hazara sustains a high number of flora throughout the seasons (Rahman et al., 2023).

Figure 1: The spatial distribution of sampling sites of Himalaya of Hazara Division, KP, Pakistan.   

The climate of the area is characterized by high seasonality, hot and humid summer, cold and snowy winter. This environmental variability combined with diverse habitats, supports  vegetation cover and a high diversity of medicinal and aromatic plants of significant economic and ecological importance (Shekhar et al., 2010).

Plant Collection: A total of 98 sampling sites were selected using a stratified sampling design to spatial heterogeneity along environmental gradients (temperature, elevation, slope angle, slope aspect), ensuring proportional representation of major habitat types and adequate spatial replication with sufficient inter-site distance to minimize spatial autocorrelation. At each site, three stratified quadrats were established to assess species composition and plant associations, with quadrat sizes based on plant habit: 10 m × 10 m for trees, 5 m × 5 m for shrubs, and 1 m × 1 m for herbaceous species, following standard phytosociological attributes. Field surveys were conducted over three consecutive growing seasons (2022–2024), with each site sampled during the same season (Matiullah et al., 2022; Lata et al., 2024).

 Density, relative density, cover, relative cover, frequency, and relative frequency of plant species were computed at each quadrat (Curtis, 1956). These parameters calculated the Importance Value Index (IVI) with the help of conventional formulas (Curtis and McIntosh, 1950). The computation of biodiversity indices, i.e., Species richness, the Shannon index, Simpson index, and the Evenness of Pilou was done using the standard formulas (Shannon, 1948; Simpson, 1949; Pielou, 1966). The collected plant specimens were dried, mounted, and labelled on herbarium sheets, and deposited in Herbarium of Hazara University, Mansehra. The geographic coordinates, including altitude, latitude, and longitude of each site were measured with the help of a Garmin eTrex 20x GPS device.

Frequency =  × 100

RF = × 100

Density =  

RD = × 100

RC=  

Canopy cover =

IVI = RD+RF+RC

Simpson’s Diversity Index

Shannon–Wiener Diversity Index

Species Evenness

Species Richness   

Soil Sampling and Field Observations: Soil samples were collected in triplicates at depths of 0–15 cm from all 98 sampling sites to ensure the representative data. Collected samples were cleaned to remove debris (stones, roots, and extraneous materials) and stored in labeled polythene bags according to sampling location before laboratory analysis at the Agricultural Research Institute, Tarnab, Peshawar, Pakistan.  At each site, environmental variables were measured, including microclimatic parameters (air temperature and relative humidity using a portable hygrometer/thermometer), soil characteristics (pH, electrical conductivity, moisture contents, organic matter, nitrogen, phosphorus, potassium, calcium carbonate, and texture), and topographic factors (elevation, slope, and aspect using GPS (Garmin eTrex 20x), clinometer (MD900-T Digital Clinometer), and compass.  Soil organic matter was determined using the dichromate oxidation method, available Nitrogen (N), Phosphorous (P), and Potassium (K) were estimated using the AB-DTPA extraction method (Soltanpour and Schwab, 1977), and soil texture (sand, silt, and clay fractions) was analyzed using the sieve method, following standard protocols (Bouyoucos, 1962).

Statistical Analyses: The data on 180 plant species of 98 sampling stands were entered into MS Excel to be statistically evaluated. Species-environment relationships and community gradients were analyzed using R software (version 4.5.2) through multivariate ordination techniques, including Detrended Correspondence Analysis (DCA), Canonical Correspondence Analysis (CCA), Principal Component Analysis (PCA), and Non-Metric Multidimensional Scaling (NMDS). While Cluster Analysis (CA), Hierarchical Cluster Analysis (HCA), Two-Way Cluster Analysis (TWCA), and Indicator Species Analysis (ISA) were conducted in PC-ORD 5.0 (McCune and Mefford, 2006). Plant communities were classified using the Sørensen distance and the Ward method of linkage to measure floristic similarity among sampling sites (Sorenson, 1948; Ward, 1963). Both presence-absence and abundance data were applied to determine associations, with ISA identifying significant indicator species (P ≤ 0.05). Additionally, the Generalized Linear Model (GLM) was used to analyze the effects of abiotic factors on species abundance. The QGIS software (Version. 3.40) was used to create distribution maps of the plant communities (QGISorg 2022).

RESULTS

Quercus species distribution in the study area: A total of five Quercus species were recorded across the study area with elevational stratification.  Quercus baloot, Q. dilatata, Q. glauca, Q. incana, and Q. semecarpifolia.  Q. glauca was found in low altitude (884–1077 MASL), while Q. baloot was abundant at lower- to mid-elevation (1123–2423 MASL). Q. dilatata and Q. incana showed broad ecological amplitude, extended from lowlands (942 MASL) to high montane sites (>2500 MASL), often co-occurring in mid-elevation forests (Table 01). In contrast, Q. semecarpifolia was restricted to high elevations (2687–3306 MASL), reflecting its adaptation to subalpine habitats (Figure 2).

Figure 2: Elevational distribution of Quercus species, ranging from lowland (Q. glauca) to subalpine (Q. semecarpifolia) zones of Himalaya of Hazara Division, Pakistan.

Plants growth form in study area: The chord diagram shows the structural composition and relationships of plant growth forms observed in the study area. Plants were classified into five forms: trees, shrubs, herbs, ferns, and climbers. Herbs were most common, followed by shrubs and trees; ferns and climbers were least represented. Connecting arcs indicate correlations and distribution patterns, revealing the floristic diversity and dominance of herbaceous vegetation. This highlights herbs' ecological role in community stability, while the abundance of woody and climbing species reflects the complexity and heterogeneity of habitats within the sites (Figure 3).

 Table 1. Elevational stratification of recorded Quercus species across the study area.

Note: DMS = degree, minutes and second, MASL = meters above sea level.

Figure 3: Chord diagram representing structural composition of plant growth forms, with herbs being the most common, followed by the shrubs and trees, and then the ferns and climbers were less represented, which represents the total floristic diversity and the heterogeneity of the habitat.

Hierarchical Cluster Analysis: Hierarchical Cluster Analysis (HCA) was conducted using Ward’s minimum variance method with Bray-Curtis distance to measure floristic similarity among sampling sites. The analysis revealed clear groupings, indicating high spatial heterogeneity in species composition. The squared distance matrix showed lower dissimilarity in sites like Akhora, Arbora, and Ayubia, suggesting high ecological affinity, while more distant sites like Banda and Bala exhibited significant compositional divergence. The dendrogram illustrated hierarchical patterns of vegetation assemblages, with both small-scale clusters of high internal similarity and larger clusters representing major vegetation types. These findings underscore the role of environmental gradients in community composition and the utility of HCA for ecological interpretation and conservation planning (Figure 4).

Figure 4: Bray-Curtis hierarchical clustering of 98 sampling sites based on the composition of 180 plant species, each colour represents a distinct cluster, indicating groups of sites with similar species and plant communities

Two-Way Cluster Analysis: To classify the plant species and sampling sites at once, Two-Way Cluster Analysis (TCA) was used to find out the coherent vegetation groups in the study area. The analysis revealed separate clusters that display good ecological patterns in species distributions. A few sites, such as Akhora, Arbora, and Ayubia, were closely clustered, indicating high floristic similarity and indicative of similar environmental conditions. On the other hand, other locations like Banda and Chaiiwai were more dissimilar and divided into discrete clusters, which suggested unique species assemblages, which might have been due to local environmental gradients or anthropogenic disturbance. The species clustering also brought out diagnostic taxa that defined each cluster, with some being dominated by shade-tolerant, moisture-loving species and others are predominantly drought-tolerant taxa due to the ecological heterogeneity of the landscape (Figure 5).

Figure 5: Two-way cluster analysis of species and sampling, black dots indicate species presence, and white dots for absence across sites in the study area.

The species-area curve: The species-area curve, based on 98 stands, indicated that 180 species had been observed. The value of richness peaked toward infinity as plots were increased (Table 01). At low sample sizes, species richness increased rapidly, averaging 18.53 species per plot. With more plots added, the rate of species accumulation decreased, reaching 103.19 species at 10 plots and 173.29 species at 60 plots. Around 80 plots, the curve leveled off. A total of 180 species were collected, indicating that most species in the study area were well represented. Jackknife estimates predicted slightly higher richness (first order: 192.9; second order: 193.0). This suggests that a small number of rare species may have been missed during sampling. The average Sørensen (Bray-Curtis) distance decreased with increasing number of plots (from 0.7970 at 1 plot to 0.0000 at 98 plots). This shows that floristic similarity among samples increased with sampling effort. This pattern confirms that the survey effort was sufficient to capture the study area's floristic diversity. The derived species pool provides a robust foundation for subsequent community analyses (Figure 6).

Figure 6: Species-area curve for the study area, showing the increase in cumulative species richness with sampling sites, the solid line represents the mean species accumulation, while the dashed lines indicate the confidence intervals, demonstrating sampling completeness across sites. 

 Species Environment Relationships: Generalized linear models (GLMs) revealed distinct species-specific responses of Quercus taxa to environmental and edaphic factors. Q. baloot showed significant negative associations with altitude, slope, and temperature, while being positively influenced by relative humidity, organic matter, nitrogen, CaCO₃, and silt content (p ≤ 0.05). Q. dilatata was positively correlated with altitude, slope, temperature, and CaCO₃, but negatively with soil pH, nitrogen, potassium, phosphorus, and sand content (p ≤ 0.05). For Q. semecarpifolia, altitude, relative humidity, and organic matter had strong positive effects, whereas slope, temperature, potassium, silt, and calcium exerted significant negative effects (p ≤ 0.05). Q. incana responded positively to altitude, soil pH, nitrogen, sand, and calcium, but negatively to relative humidity, organic matter, CaCO₃, phosphorus, and silt (p ≤ 0.05). In contrast, Q. glauca showed no significant relationships with the tested variables, as all predictors were statistically non-significant (p > 0.9). These results highlight species-specific ecological preferences and environmental filtering mechanisms shaping oak distribution across the Himalayan landscape (Figure 7).

Figure 7: Generalized Linear Model (GLM) of environmental variables on Quercus species distribution in the study area. Coefficient estimates are presented for variables with statistically significant effects (P ≤ 0.05)

Plants association

Quercus incana- Asparagus officinalis- Avena fatua Community (QAA): The Quercus incana- Asparagus officinalis- Avena fatua community (QAA) was composed of 58 species. The community was dominated by Quercus incana (38.00), Asparagus officinalis (14.98), Avena fatua (12.79), Bergenia ciliata (12.20), Alnus nitida (9.49), and Sorghum tomentosum (8.91). Co-dominant species included Hypericum perforatum, Desmodium elegans, Duchesnea indica, Cannabis sativa, Quercus dilatata, and Arundinaria donax. The associated species included Pyrus pashia, Cynodon dactylon, Cedrus deodara, Ficus carica, Prunus cornuta, Galium aparine, and Medicago denticulata. The least important species included Zanthoxylum armatum, Rubus fruticosus, Morus alba, Phragmites communis, Bergenia stracheyi, and Malotus philippensis. The QAA community was found at a mean altitude of 1739 MASL with moderately steep slopes (26.6°) and a mean temperature of 23.1 °C. Relative humidity averaged 54.4%, while soils showed near-neutral pH (6.32) and moderate organic matter content (0.95%) (Figure 9: b). This community had moderate species richness (14.6 species) and a Shannon diversity index of 2.57, with high evenness (Pielou’s index = 0.96) (Figure 9: a). The combination of mid-altitude, warm climate, and well-balanced soil properties suggests that QAA supports a relatively stable assemblage with high species evenness but lower richness than higher-elevation communities.

Quercus incana- Aesculus indica- Phragmites communis Community (QAP): The Quercus incana- Aesculus indica- Phragmites communis community (QAP) supported 94 species. The community was dominated by Quercus incana (52.23), Aesculus indica (7.63), Phragmites communis (7.09), Taraxacum officinale (5.47), Pyrus pashia (5.06), and Ficus carica (4.74). Co-dominant species included Urtica dioica, Impatiens balsamina, Cynodon dactylon, Ranunculus muricatus, and Poa annua. The associated species included Pinus wallichiana, Phleum pratense, Rubus fruticosus, Solanum nigrum, and Medicago denticulata. The least represented species were Debregeasia salicifolia, Hypericum perforatum, Hedera helix, Ziziphus jujuba, Robinia pseudoacacia, and Rubus ellipticus. The QAP community occurred at slightly lower elevations (1690 MASL) with steeper slopes (29.7°). The environment was characterized by warmer mean temperatures (24.3 °C), moderate relative humidity (55.0%), and slightly acidic soils (pH 6.62) with lower electrical conductivity (Figure 9: b). Richness was higher than QAA (17.8 species), accompanied by Shannon diversity of 2.65 (Figure 9: a). Its evenness (0.93) was slightly lower than QAA, but it was a richer community with more species supported. QAP appears to be a higher nitrogen, mid-elevation assemblage with robust species representativeness under relatively warmer, mesic conditions and fertile soils.

Quercus semecarpifolia- Bergenia stracheyi- Indigofera heterantha Community (QBI): The Quercus semecarpifolia- Bergenia stracheyi- Indigofera heterantha community (QBI) consisted of 46 species. The community was strongly dominated by Quercus semecarpifolia (68.96), accompanied by Bergenia stracheyi (27.48), Indigofera heterantha (25.35), Skimmia laureola (18.40), Picea smithiana (13.97), and Cedrus deodara (13.46). Co-dominant species included Oenothera rosea, Viburnum grandiflorum, Inula helenium, Taxus wallichiana, Aquilegia nivalis, and Trillium govanianum. Associated species were represented by Salix flabellaris, Iris hookeriana, Quercus dilatata, Poa annua, Phragmites communis, and Celtis australis. The least contributing species were Rosa moschata, Cichorium intybus, Bergenia ciliata, Impatiens cylindrica, and Viburnum cotonefolium. The QBI community was distinctly montane, located at the highest mean altitude (3104 MASL) and on the steepest slopes (36.6°). It experienced the lowest mean temperatures (16.6 °C) and relatively high humidity (60.6%). Soils were more acidic (pH 5.12) with the highest organic matter content (1.92%) (Figure 9: b). Species richness was the lowest (13.8 species) with reduced Shannon diversity (2.37) and evenness (0.93) (Figure 9: a). This indicates that QBI is a stress-prone, high-altitude community in which few specialized taxa dominate, resulting in low richness but a relatively balanced representation. The adaptation of Quercus semecarpifolia and its associates to cold, moist conditions reflects the specialized nature of this community.

Quercus dilatata- Cannabis sativa-Quercus glauca Community (QCQ): The Quercus dilatata- Cannabis sativa-Quercus glauca community (QCQ) comprised 92 species. The dominant species were Quercus dilatata (29.29), Cannabis sativa (9.74), Quercus glauca (9.50), Morus alba (8.85), Fragaria vesca (6.63), and Avena fatua (6.08). Codominant included Juglans regia, Cynodon dactylon, Bauhinia variegata, Quercus incana, Desmodium elegans, and Duchesnea indica. Associated species included Olea ferruginea, Bupleurum sativus, Spiraea vaccinifolia, Andrachne cordifolia, Vitex negundo, Sorbaria tomentosa, and Hypericum perforatum. The least represented species were Impatiens balsamina, Heteropogon contortus, Inula helenium, Dalbergia sissoo, Rubus fruticosus, and Viburnum cotonefolium. The QCQ community occurred at lower altitudes (1590 MASL) with gentle to moderate slopes (29.3°). The mean temperature was the highest among all communities (25.7 °C), with relative humidity averaging 54.4%. Soils were slightly acidic (pH 6.65) and had moderate organic matter (0.68%) (Figure 9: b). Species richness was moderate (15.9 species) with Shannon diversity of 2.60 and high evenness (0.95) (Figure 9: a). This community reflects a transition zone where warm conditions and fertile soils support a balanced set of species, with diverse values suggesting relatively stable coexistence and reduced dominance effects.

Quercus incana-Prunus vulgaris- Quercus dilatata Community (QPQ): The Quercus incana-Prunus vulgaris-Quercus dilatata community (QPQ) concealed 89 species. Dominant species included Quercus incana (28.67), Prunus vulgaris (12.83), Quercus dilatata (9.85), Impatiens cylindrica (9.26), Pinus wallichiana (7.24), and Cynodon dactylon (6.23). Co-dominant species were Cannabis sativa, Chenopodium album, Pyrus pashia, and Artemisia absinthium. The associated group included Arisaema flavum, Achyranthes aspera, Pinus roxburghii, Adiantum incisum, Abies pindrow, Berberis lycium, and Morus alba. Species with the least contribution were Sonchus asper, Mentha longifolia, Malva sylvestris, Urtica dioica, Rubus ellipticus, and Dryopteris stewartii. The QPQ community was distributed at mid- to high elevations (1863 MASL), with slopes averaging 28°. The mean temperature was 22.4 °C, the relative humidity was moderately high (57.6%), and the soil pH averaged 6.54 with organic matter 0.82% (Figure 9: b). Species richness (21.3 species) and Shannon diversity (2.76) were high, although evenness (0.91) was slightly lower, indicating dominance by a few species within a diverse assemblage (Figure 9: a). The combination of intermediate climate and fertile soils contributed to high diversity, making QPQ one of the most species-rich Quercus communities.

Quercus baloot-Quercus dilatata- Cynodon dactylon Community (QQC): The Quercus baloot-Quercus dilatata- Cynodon dactylon community (QQC) was represented by 104 species. The community was dominated by Quercus baloot (IVI = 30.91) and Cynodon dactylon (24.81), followed by Quercus dilatata (9.97), Pinus wallichiana (7.88), Berberis lycium (7.51), and Impatiens bicolor (6.62). Co-dominant species included Artemisia absinthium, Olea ferruginea, Punica granatum, Sassurea lappa, Cannabis sativa, and Juglans regia. A wide range of associated species were present, such as Polygonum amplexicaulis, Morus alba, Ficus palmata, Duchesnea indica, Ajuga bracteosa, Aquilegia nivalis, and Malotus philippensis. Species with the lowest contribution to the community included Sisymbrium officinale, Trilium gavanianum, Sarcococca saligna, Emilia sonchifolia, and Salvia moorcroftiana. The QQC community occurred at elevations around 2009 MASL with steep slopes (33.8°). This community experienced the lowest mean temperatures (15.9 °C) after QBI and the highest relative humidity (60.0%). Soils were moderately acidic (pH 5.79), with organic matter averaging 0.93% (Figure 8b). It had the highest species richness (22.3 species) and Shannon diversity (2.78), although evenness was lower (0.90) than in QAA and QCQ (Figure 8a). The colder, moister conditions, combined with fertile soils, create favorable niches for diverse taxa, explaining the high richness observed. This community represents a biodiversity hotspot within the Quercus-dominated vegetation.

Figure 8: (a) Biodiversity indices (e.g., species richness, Shannon diversity, and evenness) and (b) environmental variables across six plant communities in the study area. Quercus incana-Asparagus officinalis-Avena fatua (QAA), Quercus incana-Aesculus indica-Phragmites communis (QAP), Quercus semecarpifolia-Bergenia stracheyi-Indigofera heterantha (QBI), Quercus dilatata-Cannabis sativa-Quercus glauca (QCQ), Quercus incana-Prunus vulgaris-Quercus dilatata (QPQ), and Quercus baloot-Quercus dilatata-Cynodon dactylon (QQC).

NMDS and PCA of communities: The NMDS analysis shows that vegetation distribution is mainly influenced by climatic and geographic factors. Temperature (r = 0.505, R² = 0.255) and relative humidity (r = 0.376, R² = 0.141) were the strongest climatic drivers, while latitude (r = 0.438, R² = 0.192) and altitude (r = 0.236, R² = 0.056) also had significant effects. Among edaphic variables, soil pH, organic matter, and clay showed weaker but significant relationships, whereas most other soil parameters were non-significant (Figure 9).

The principal component analysis (PCA) summarized the variation in the dataset, with the first three axes explaining 13.6%, 9.6%, and 7.6% of the total variance, respectively, accounting for nearly one-third (30.2%) of the cumulative variation. PC1 (9.9 eigenvalue) represented the strongest gradient, followed by PC2 (7.5) and PC3 (7.2), each contributing meaningful ecological separation among sites. The gradual decline in explained variance across subsequent components indicates that the major patterns of species distribution and environmental variation are well captured by the first few axes, while higher-order axes contributed only marginally (Figure 10).

Figure 9: Non-metric multidimensional scaling (NMDS) ordination showing the influence of Environmental variables on distribution and classification of plant communities in the study area. Quercus incana-Asparagus officinalis-Avena fatua(QAA), Quercus incana-Aesculus indica-Phragmites communis (QAP), Quercus semecarpifolia-Bergenia stracheyi-Indigofera heterantha (QBI), Quercus dilatata-Cannabis sativa-Quercus glauca (QCQ), Quercus incana-Prunus vulgaris-Quercus dilatata (QPQ), and Quercus baloot-Quercus dilatata-Cynodon dactylon (QQC). Distinct colour indicates the different communities

Figure 10: PCA biplot of 98 sampling stands (180 species) showing site distribution along environmental gradients; PC1 (30.2%) correlated positively with latitude, altitude, and relative humidity, and negatively with temperature and soil pH, while PC2 (12.3%) reflected variation in soil properties. The size of the circle shows the species richness while the color shows different communities. Quercus incana-Asparagus officinalis-Avena fatua(QAA), Quercus incana-Aesculus indica-Phragmites communis (QAP), Quercus semecarpifolia-Bergenia stracheyi-Indigofera heterantha (QBI), Quercus dilatata-Cannabis sativa-Quercus glauca (QCQ), Quercus incana-Prunus vulgaris-Quercus dilatata (QPQ), and Quercus baloot-Quercus dilatata-Cynodon dactylon (QQC).

Detrended correspondence analysis: The response data showed a gradient length of 4.20 SD units along axis 1, indicating a long ecological gradient. In DCA ordination, axis 1 had a gradient length of 4.20 and an eigenvalue of 0.53; axis 2, 2.89 and 0.30; and axis 3, 3.20 and 0.23. The eigenvalues of 0.53 and 0.30 on axes 1 and 2 indicated a strong species-habitat correlation, while the lower eigenvalue of 0.23 on axis 3 showed a weaker relationship. Total variance in the species data was 7.85 (Table 02). The dominant gradient along axis 1 suggested a resemblance in species composition.

DCA ordination grouped species by ecological affinities. Taxus wallichiana, Quercus semecarpifolia, Bergenia stracheyi, Aquilegia nivalis, and Viburnum grandiflorum showed strong positive correlation along axis 1. Indigofera heterantha, Quercus dilatata, Rosa webbiana, Juglans regia, and Impatiens bicolor were positively correlated on axis 2. Mallotus philippensis and Olea ferruginea were near the centroids of axes 1 and 2. Cynodon dactylon and Chenopodium album were positioned away from the centroid and showed negative correlations. Alnus nitida and Ranunculus muricatus were positively associated along axis 1, but negatively correlated with Carissa opaca and Ziziphus jujuba. Most species were centered in the ordination diagram, reflecting identical or overlapping communities and habitat preferences.

The ordination of stands also showed well-separated clusters based on species composition. Stand 27 (Dana) and stand 32 (Donga Galli) were closely related along axis 1, which represented similar vegetation structure. In contrast, stand 30 (Datal Banda) was inversely correlated with stands 26 (Dalola) and 28 (Danat). Stand 25 (Dabben) appeared at the upper end of the ordination space, indicating high species diversity, and showed a negative correlation with stands 27 and 32. Most stands, including Banbaik (stand 8), Chorkalan (stand 24), and Surbanda (stand 94), are grouped at the center of the ordination diagram. They showed similar species composition, with some overlap in stands (Figure 11).

Figure 11: Illustrating DCA analysis of plant species and sites of the study area. The triangle shape shows the sites, and colours indicate distinct clusters based on species composition.

 Table 2. Detrended Correspondence Analysis (DCA) summary showing eigenvalues and gradient lengths for the first three axes. Eigenvalues represent variance explained, while gradient lengths (SD units) indicate species turnover along ecological gradients.

Canonical Correspondence Analysis (CCA): CCA was applied to study the influence of environmental variables on plant species composition, especially Quercus species. The eigenvalues of the first three axes were 0.412, 0.298, and 0.192, corresponding to explaining 28.4%, 20.5% and 13.2% variance on the species–environment relationships (Table 03). Together, the two axes explained 48.9% of the variation and therefore should encompass most of the major ecological gradients accounting for species distribution.

The ordination biplot showed clearly differing ecological preferences of the Quercus species. Q. semecarpifolia was loaded on the positive side of axis 1 and strongly associated with higher elevation, lower temperature, and dense canopy cover. This indicates its preference to damp temperate forest environments. Quercus dilatata, in contrast, was clustered toward the positive side of axis 2, where it was strongly correlated with soil moisture and organic matter, reflecting its affinity for relatively mesic habitats. Quercus baloot was separated toward the negative side of axis 1, indicating its preference for warmer, drier conditions and lower elevation sites, consistent with its distribution in dry temperate forests. Q. incana showed positive correlations with both Axis 1 and Axis 2, indicative of its mid-elevational distribution; conversely, Q. glauca exhibited a weaker positive correlation with Axis 1 alongside a slight negative correlation with Axis 2, which positions it at the lowest elevations and reflects its divergent moisture affinities compared to Q. baloot (Figure 13, 14).

Other co-occurring species further supported these habitat associations. Taxus wallichiana and Indigofera heterantha were closely aligned with Q. semecarpifolia along axis 1, reflecting shared preferences for high-altitude moist habitats. Rosa webbiana and Berberis lycium clustered near Q. dilatata, suggesting a common response to moisture gradients. On the other hand, Cynodon dactylon and Chenopodium album were positioned opposite Q. baloot, indicating their presence in disturbed and open habitats where Q. baloot was absent.

The ordination of stands further reflected these species–environment relationships. High-altitude stands such as Dana (stand 27) and Donga Gali (stand 32) grouped strongly with Q. semecarpifolia, while middle-elevation stands (e.g., Dabben - stand 25, Dalola -stand 26) were associated with Q. dilatata. Low elevation stands, including Datal Banda (stand 30) were positioned near Q. baloot, highlighting its adaptation to relatively dry and disturbed environments. Most stands were situated near the center of the ordination diagram, showing overlapping communities and transitional zones between Quercus-dominated habitats (Figure 12).

Figure 12: CCA Ordination of Species and Environmental Variables in the study area.

Table 3. Monte Carlo permutation test results for Canonical Correspondence Analysis (CCA) axes, showing eigenvalues, canonical correlations, percentage variance explained, and cumulative variance for the first four axes.

Figure 13: The species environment relationship of Quercus incana and Quercus glauca.

Quercus incana shows a positive correlation with CCA1 (r = 0.560, τ = 0.417) and CCA2 (r = 0.278, τ = 0.294). Quercus glauca shows a positive correlation with CCA1 (r = 0.250, τ = 0.269) and negative correlation CCA2 (r = -0.068, τ = -0.129).

Figure 14: The species environment relationship of Quercus dilatata, Quercus baloot, and Quercus semecarpifolia.

Q. baloot shows strong negative correlations with both axes (r = –5.58 to –7.42), suggesting its association with the environmental gradient captured by the first two axes. Q. semecarpifolia also exhibits moderate negative correlations (r = –6.39 with Axis 1; –5.60 with Axis 2), whereas Q. dilatata displays weak correlations (r = –0.10 to –1.24). The eigenvalue for Axis 1 (0.53) and total explained variation (56–62%) indicate that the CCA model effectively captures the principal ecological gradients structuring Quercus‑dominated forests in the study area.

DISCUSSION

The present study recorded 180 species from 142 genera and 65 families across Quercus-dominated forests of the western Himalayas, confirming the ecological richness of this region. Similar levels of diversity have been reported in other Himalayan forests, where altitudinal gradients, microclimate, and habitat heterogeneity create favorable conditions for diverse taxa (Rahman et al., 2022; Singh and Negi, 2023; Yadav and Pandey, 2024).

Asteraceae, Rosaceae, Lamiaceae, and Poaceae were the dominant families, with Asteraceae alone contributing more than one-fifth of the flora. This pattern is consistent with earlier studies across Pakistan and Himalayan forests, where these families repeatedly dominate due to their broad ecological amplitude and adaptive dispersal mechanisms (Qazi et al., 2023; Gillani et al., 2024). Asteraceae and Lamiaceae have also been reported to dominate alpine vegetation, indicating that these species are adapted from forest regions to higher altitudes (Shaheen et al., 2019; Qaiser et al., 2025).

A pronounced elevational stratification of Quercus species was observed, highlighting the central role of environmental filtering in structuring forest communities. Species such as Q. glauca were confined to lower elevation, consistent with their affinity for warmer, subtropical conditions  (Jin et al., 2023), however, Q. baloot occupied the lower and mid altitudinal range as reported previously for this keystone dry temperate Himalayan forest tree species (Ahmad et al., 2022). Intermediate species, including Q. incana and Q. dilatata, exhibited broader ecological amplitudes, spanning multiple elevational zones. This pattern is consistent with niche theory, where species distributions are governed by physiological tolerances and competitive abilities along environmental gradients (Rahman et al., 2023), whereas Q. semecarpifolia was restricted to high-elevation subalpine zones, reflecting adaptations to cold, moist, and physiologically stressful environments (Rawat et al., 2022; Shekhar et al., 2022). The western Himalayas, situated at a biogeographic transition zone, further amplify these gradients, facilitating species turnover and elevational replacement (Singh et al., 2009).

The results of multivariate classification (HCA and TCA) showed that there were separate vegetation assemblages and diagnostic species, which revealed a strong spatial structuring of vegetation communities. These curves are in line with the ecological theory which posits that montane ecosystems are gradients in moisture availability, light intensities, and disturbance regimes that drive species assembly (Rahman et al., 2022). Identifying diagnostic species within these clusters is crucial, as these taxa define vegetation types and contribute to conservation priorities, as emphasized by (De Cáceres et al., 2010).