Typology of Pure Deodar Forests Driven by Vegetation–Environment Relations in Manoor Valley, Northwestern Himalaya

: The current research was carried out to characterize the phytosociology of the forests of one of Pakistan’s most valuable tree species (Deodar) across its native range. In this context, our main hypothesis was that, along the altitudinal gradient, we would ﬁnd different plant communities that would be driven by different environmental variables (climatic, edaphic, and physiographic). Therefore, to assess the vegetation structure of the pure Deodar forests of the unexplored Manoor Valley (Northwestern Himalaya), Pakistan, frequent ﬁeld visits were carried out during different seasons of 2015–2018. Ecological methods: Line transects sampling (23 stands) and phytosociological attributes were evaluated in relation to geographical and environmental variables. Various statistical software applications (i.e., PCORD, RStudio 4.0, and R 3.6.1) were used to examine all of the gathered data of plant species and environmental variables. A total of three different plant communities (


Introduction
The study of vegetation categorization based on species co-occurrence and its relationship to environmental factors is known as phytosociology [1,2]. This research area has resulted in major vegetation appraisal techniques and methodologies that may be linked to ecosystem services and biodiversity conservation [3]. Species composition, community structure, and function are the most significant ecological features of forest ecosystems, and they all change in response to environmental changes [4,5]. Factors such as vegetation type, slope, aspect, edaphic variables, and altitude [6,7] influence the community composition, structure, and distribution pattern of diversity in mountain vegetation [8]. Elevation is an important variable in mountain ecosystems [9] since it affects the structure of the vegetation in most mountains across the world [10]. Most of the environmental factors, as well as species diversity [11,12], change concurrently along the altitudinal gradient [13]. This aspect fosters habitat heterogeneity and induces micro-environmental alteration in vegetation patterning [14,15]. Heywood and Watson [16] reported that species richness in vascular plants increases with temperate latitude. Vegetation associations with ecological diversity are considered as a proportion of the quality of the entire ecosystem [17]. For this reason, vegetation in relation to environmental variables has been an important topic in recent years [18].
In Pakistan, Deodar (Cedrus deodara (Roxb. ex D.Don) G.Don) is found in Murree, Hazara, Abbottabad, Swat, Azad Kashmir, Kaghan Valley, Kohistan Chillas, Dir, and Chitral. It is mostly a dry temperate plant, although scattered stands of this species can be found in moist temperate regions [19]. A majority of these forests are found in the transitional zone between dry and moist temperate zones, with no clear boundary between the two [20]. Deodar grows on ridge tops and moderate (17 • ) to steep (50 • ) slopes in Northern areas, ranging in elevation from 1650 m (Kaghan) to 2927 m (Chitral) in moist temperate, dry temperate, and timberline regions [21]. According to Champion et al. [19], Cedrus deodara forests gradually extend into the dry inner valleys of the Himalayas. Hussain and Illahi [20] have compiled a comprehensive report on the ecology and vegetation types of Pakistan's lesser Himalaya. They further noted that Deodar trees cannot endure extremely moist environments and are thus restricted to the Himalayan drier zones.
Following the observational vegetation surveys by Champion et al. [19] and Beg [22], plenty of quantitative phytosociological studies have been conducted across the nation by several researchers. Ahmed et al. [23] carried out an extensive quantitative sampling in the Deodar forests of the Hindu Kush and Himalayan ranges of Pakistan, and they identified Cedrus deodara, Pinus wallichiana, P. gerardiana, and Abies pindrow as the dominant species in the recorded plant communities. Phytosociological work from Ayubia National Park [24] was reported three decades ago, and it was reported about a decade ago in the Sarsawa Hills, District Kotli of Azad Kashmir [25]. The former study reported five plant communities, i.e., Acacia modesta-Cannabis sativa, Acacia modesta-Cynodon dactylon, Acacia modesta-Themeda anathera, Acacia modesta-Dodonaea viscosa, and Acacia modesta-Lantana camara. The latter identified nine communities (i.e., Pinus-Poa-Maytenus, Myrsine-Themeda-Pinus, Colebrookia-Themeda-Dodonaea, Themeda-Carissa-Adhatoda, Themeda-Dodonaea-Eriophorum, Adhatoda-Themeda, Carissa-Myrsine-Themeda, Carissa-Themeda-Dodonaea, and Dodonaea-Carissa-Pinus) at different elevations. The structure and distribution patterns of plant communities in Himalayan forests are poorly understood [26]. The Northwestern Himalaya is a unique bioregion because of its variable topographical and ecological complexity throughout a vast altitudinal range.
Nonetheless, apart from these studies, no thorough quantitative analyses have been conducted that define communities of Deodar forests, wherein the most valuable wood species reside, designating the Deodar as Pakistan's national tree. Therefore, the current research was carried out to characterize the phytosociology of the forests, owing to one of Pakistan's most valuable tree species across its native range. In this context, our main hypothesis was that, along the altitudinal gradient, we would find different plant communities that would be driven by different environmental variables (climatic, edaphic, and physiographic), due to changes in environmental factors. Overall, our study aims to move toward a better understanding of how environmental variables influence the structuring and distribution of plant species in the Deodar forest.

Study Area
The current field investigation was conducted in the pure Deodar forests of the Manoor valley, which is a hilly valley ( Figure 1) in Pakistan's Northwestern Himalayan region [27][28][29][30], with elevations ranging from 1580.8 m to 2373.8 m. In the Himalayas, monsoon winds are the principal source of precipitation as well as the primary governing factor of erosion, topography, environment, and vegetation [31]. The valley's vegetation features Sino-Japanese characteristics [32,33]. Summers are chilly and dry, and the valley has a dry-temperate climate with considerable seasonal variations. Between November and April, there is a lot of snow (an average annual snowfall of 3 m). With rising altitude, the range represents a significant spike in snow depth. In comparison to lower elevations, the winter is severe at higher altitudes [34]. From January to April, there is a distinct wet season, with June to November being the driest months of the year.

Vegetation Sampling and Plant Identification
The vegetation of pure Deodar forests was surveyed and quantified [35] along the geographic, slope, edaphic, and climatic variables [36]. A line transect ecological technique was used for the vegetation sampling [35,[37][38][39][40][41]. Sampling was carried out during 2015-2018. The surveyed sites were divided into 23 stands, and three sampling points within each stand were sampled along a 50 meter transect (total = 69 transects). The phytosociological attributes (i.e., density, frequency, cover and their relative values, and importance value) were used for the calculation of the data from each sampled point [42][43][44]. The IV was further used to rank each species, and those with the highest IV were considered as the dominant species [42,43]. Similarly, plant communities were named after three dominant species [45][46][47][48][49]. The slope angle, aspect (i.e., east, west, south, and north), and exposure of each sampled point were recorded using a clinometer, and altitude, longitude, and latitude were measured by the geographical positioning system (GPS). Plant collection, labeling, pressing, and other herbarium work methodologies were adopted [50][51][52][53][54][55][56]. Their identification was confirmed using Flora of Pakistan [57][58][59], and they were submitted to the Hazara University herbarium in Mansehra (Pakistan).

Environmental Gradients
Soil samples weighing 200-400 g were collected from three randomly selected transects (0-30 cm depth) within the studied vegetation area [60,61]. The replicated samples collected from each stand were thoroughly mixed to make a composite sample [62], stored in a sterile polythene bag, and labeled with a permanent marker. The samples were cleaned of any raw materials, such as rocks and stones, and then dried in the shade. Physicochemical tests were performed, including soil texture (clay, sand, silt, and loam), pH [63], electrical conductivity (EC) [64], organic matter (OM) [65], nitrogen (N), phosphorus (P), potassium (K), and calcium carbonate (CaCO3) contents [66]. Other climatic variables (such as temperature, barometric pressure (BP), wet bulb (WB), dew point (DP), heat index (HI), wind speed(WS), and humidity) were measured by recording data at each transect with a small remote weather station (Kestrel weather tracker 4000), and then calculating average values at the stand level [36].

Statistical Analyses
The data from all of the plant species, geographic, and other environmental gradients were compiled to see whether there was a link between them [67,68]. Various software packages, such as PC-ORD [69] and R 3.6.1 [70], were applied to analyze IV matrices (abundance of the 162 species) and values of geographic and environmental gradients obtained at 23 stands. PC-ORD software [35] was used to identify plant communities [36] by processing cluster analysis (CA), two-way cluster analysis (TWCA) [71], and two-way indicator species analysis (TWINSPAN). TWINSPAN was processed for the identification of the plant communities [35] based on species classification and sample clustering [72]. RStudio 4.0.0 [70] was used for processing non-multidimensional scaling ordination (NMDS) and principal component analysis (PCA) [73], ternary plot using the "Ternary" package [74], GLM with Gaussian error distribution and Likelihood ratio test [75], the spatial turnover and nestedness-resultant components of beta diversity, dissimilarity analysis [76,77], canonical correspondence analysis (CCA), and variation partitioning tests (partial CCA) [70]. The NMDS and PCA bi-plot ordinations were used to evaluate the correlation of environmental gradients with communities. The ternary diagram (simplex plots/Gibbs triangles) was drawn using the standard graphics functions by employing the soil texture variables data (sand, silt, and clay). GLM was conducted to compare the variables (20 variables distributed into four different groups: geographic, edaphic, climatic, and slope) evaluated among the communities produced by our previous analysis [75]. For calculating species richness by using a Poisson error distribution, we used a GLM, followed by a Likelihood ratio test. In addition, we used Gaussian error to measure each plant community's evenness, Shannon, and Simpson diversity. In the case of CCA and partial CCA, analyses fitted a complete model, and by using the step function with permutation in the package "stats", we reduced this complete model to the best model with the lowest number of variables. In addition, we used the Variance Inflation Factor (VIF) to validate the model's applicability by evaluating multicollinearity between variables in the final model.

Results
In this phytosociological investigation, 162 plant species were documented from 23 stands sampled in the pure Deodar forests of the Manoor Valley, Himalaya, Pakistan ( Table 1). The pure Deodar forests ranged from 1580.8 to 2373.8 m. Table 1. List of plant species, with their importance values (IVs) calculated based on vegetative characteristics of each stand. Plant communities were named based upon three dominant species. The indicator species' IVs are mentioned in bold, as well, with respect to plant community.

Cluster Analysis (CA), Two-Way Cluster Analysis (TWCA) and TWINSPAN
Cluster analysis (CA), using PC-ORD, grouped the vegetation into three plant communities, as indicated in the dendrograms under the impact of various environmental gradients ( Figure 2). TWCA displays the detailed distribution patterns of 162 plant species recorded at 23 sampling sites (stands). Based on influenced environmental gradients, all of these stands were clustered into three distinct groups (plant communities) (Figure 3). According to TWINSPAN classification, a total of three communities (Cedrus-Isodon-Cynodon, Cedrus-Cynodon-Dryopteris, and Sambucus-Cedrus-Desmodium) were established by clustering 162 reported plants species in 23 stands under the strong influence of environmental variables in pure Deodar forests (Table 1).

Community Distribution Modeling along the Environmental Gradients
NMDS (Figure 4a-d) and PCA are used to show the relationship between identified communities and environmental gradients (Figure 4e). The most representative environmental factors that drove the plant community structure were altitude, slope angle and aspects, NE, slope SW, K, pH, OM, loam, silt, sand, clay, temperature, HI, and WS ( Table 2). Environmental gradients distinguish the 23 sample sites into 3 communities in both ordinations, as TWINSPAN, CA (Figure 2), and TWCA have already revealed (Figure 3).  In constrained PCA ordination, the PC1 axis accounted for the most explanatory variance (21.6%), whereas the PC2 axis accounted for the least (14%). By distinguishing the vegetation of the pure Deodar forest into three communities (Figure 4a-e), the considerable effect of the altitudinal gradient was revealed. Altitude (1936-2373 m), slope angle (25-85 • ), and sandy (29-48%) and loamy soil texture all indicated a positive and substantial association in the SCD community. Potassium (216 mg kg −1 ), EC (3.3 dsm −1 ), CaCO 3 (8.45 mg kg −1 ), and organic matter (1.03%) were the edaphic variables that shaped this plant community. Furthermore, wind speed (1.45 ms −1 ) and temperature (25.8 • C) were found to have a positive and significant association. In contrast with this, the CCD community was primarily found in the lower mountainous ranges (1580.8-1982 m.a.s.l.) and showed positive and significant relationships with the northeastern slope, silty (32-58%) and sandy (15.8-55%) loamy soil texture, and barometric pressure (814.3 pa). Nonetheless, the CIC community was recognized at the middle altitudinal range (1769.8-1869.5 m). This community revealed a significant positive correlation with the northeastern to southwestern slope, pH (6.3), wet bulb (19.7), and dew point (17.7). Since silt, sand, and clay were representative in both NMDS and PCA ( Figure 4) and are expressed in percentage (Table 2), the relationship of CIC, CCD, and SCD plant communities was plotted with the said three edaphic elements in a ternary plot ( Figure 5). Likewise, a clear clustering in plant communities was observed according to soil composition based on these three elements. Based on the mean values, the CCD and CIC communities were hosted by soils having the maximum silty texture, i.e., 43.7% and 42.8%, respectively. On the other hand, the SCD community was hosted on soil that comprised a higher sandy texture (41.7%). As each plant community comprised several stands, the soil-texture values recorded for each stand were found to have no significant differences ( Figure 5). These results reveal that the tree layer determines the structure of the ground cover more than soil does.

Significance Testing of Pure Deodar Forest Communities in Relation to Studied Variables
The significance testing of recognized plant communities was performed in association with geographic and environmental gradients. The results of the GLM analyses revealed that the majority of the gradients were significantly different (p-value < 0.05; Figure 6) among the 20 examined gradients between three communities. Nonetheless, there were no significant variations in the heat index, temperature, DP, WB, OM, K, P, and CaCO 3 ( Table 2).

Beta Diversity (βsim and βsne)
Beta diversity revealed a dissimilarity lower than 50% among the three communities (Figure 7). In the βsim cluster, we observed a 47% dissimilarity between the SCD community and the other two communities and a 15% dissimilarity between CCD and CCI. In the βsne cluster, the highest dissimilarity value was 28% between CIC and the other two communities (CCD and SCD). The CCD and SCD showed 18% dissimilarity. Thus, plant community structure is strongly influenced by βsim as compared to βsne, since βsim showed almost twice the dissimilarity among communities than βsne.

Species Richness and Diversity Indices
The maximum number of associated species was recorded in the SCD community (129 species, Figure 8a) at 1789.6-1896.3 m altitude, followed by CCD (66 species), and the minimum species was reported in the CIC community (38 species). The SCD community has the highest value (H' = 3.7), followed by the CCD community (H' = 3.6) and the CIC community (H' = 3.3, Figure 8b). The SCD community has the highest Simpson's dominance value (0.98), while the CIC has the lowest value (0.96, Figure 8c). Nevertheless, the analysis revealed the CIC community with the highest value (0.97) for Pielou's evenness index, while the SCD community had the lowest value (0.93, Figure 8d).

Discussion
In this phytosociological investigation, 162 plant species were documented from 23 stands sampled in the pure Deodar forests (Manoor Valley, Himalaya, Pakistan), ranging from 1580.8 to 2373.8 m. The species richness recorded in the current study is higher than that of temperate forests in Kedernath Wildlife Sanctuary (116 species), the Central Himalaya, India [79], and the Western Himalayan moist temperate forests in Kashmir, Pakistan (122 species) [5]. Dar and Sundarapandian [80] observed that mid-elevation (2300-2800 m) coniferous forest types had higher species richness than low and highelevation broad-leaved forests. This might be because mid-elevation ranges are less disturbed than low-elevation ranges. Although altitude is the most important variable in regulating species distribution, other variables, such as topography, aspect, slope, and exposure, can significantly affect the structure and composition of species [36,81].
In the current ecological investigation, various types of analyses have been used to quantify the recorded data of all species and stands in association with geographic and environmental gradients. The structure and distribution patterns of plant communities in Himalayan forests are poorly understood [26]. Due to variable topographical and ecological complexity throughout a vast elevational range, the Northwestern Himalaya is a unique bioregion. The TWINSPAN, CA, and TWCA resulted in the recognition of plant communities in the pure Deodar forests of the studied area. Each plant community was composed of many stands, with no notable changes in the soil-texture values recorded for each stand. The considerable effect of the altitudinal gradient was revealed by distinguishing the vegetation of the pure Deodar forest into three major plant communities. Moreover, these three communities (Cedrus-Isodon-Cynodon, Cedrus-Cynodon-Dryopteris, and Sambucus-Cedrus-Desmodium) were established by clustering 162 species recorded from 23 stands under the effect of climatic, edaphic, topographic, and physiographic gradients in pure Deodar forest that ranged from 1580.8 to 2373.8 m. Ahmed et al. [21] recognized seven plant communities in 47 stands of Deodar forests in the Himalayan range of Pakistan.
The SCD community showed a strong influence with altitude (1936-2373 m), slope angle , sandy (29-48%) and loamy soil texture, wind speed (1.45 ms −1 ), and temperature (25.8 • C). In contrast with this, the CCD community showed a positively significant relationship with the northeastern slope, silty (32-58%) and sandy (15.8-55%) loamy soil texture, and barometric pressure (814.3 pa). Nonetheless, the CIC community revealed a significant positive association with the northeastern to southwestern slope, pH (6.3), wet bulb (19.7), and dew point (17.7). A similar pattern was observed in the allied area (Nandiar catchment, Battagram) of the Himalayan region by stating altitude as the governing gradient [35]. Such variables alter the species diversity as well as community structure [82,83]. Species composition, community structure, and species diversity, however, are influenced not just by environmental factors but also by local micro-abiotic filtration [84].
We found significant differences (p < 0.001) among the three communities found in the pure Deodar forests in the four indexes. The SCD community had the maximum number of species (129 species) at 1789.6-1896.3 m altitude. Variability in the number of plant species is an outcome of species interaction with a particular set of environments [85]. A number of biological and environmental factors interact to control the distribution of species richness along the altitude [86][87][88]. Nonetheless, other researchers [16] reported that species richness in vascular plants increases with temperate latitude. The SCD community has the maximum Shannon's diversity (H' = 3.7) and Simpson's dominance (0.98) values among the recorded communities. The diversity values recorded in the current study are higher than the range (1. 16-3.40) reported by various researchers from other Himalayan regions [5,89]. The Pielou's evenness index value was led by the CIC community (0.97). Many researchers reported that the physicochemical properties influence plant species richness [90] and evenness [91].

Conclusions
The current research was carried out to characterize the phytosociology of the forests, owing to one of Pakistan's most valuable tree species across its native range. In this context, this investigation documented 162 plant species recorded in 23 stands of the said forests that were further recognized into three plant communities. Each plant community was composed of many stands, with no notable changes in the soil-texture values recorded for each stand. Therefore, it could be concluded that the tree layer, rather than the soil texture, determines the structure of the ground cover. Significant differences were recorded among three communities found in the pure Deodar forests by analyzing species richness (χ 2 = 14.732, p < 0.001) and Pielou (χ 2 = 18.267, p < 0.001). Beta diversity showed a dissimilarity lower than 50% among the three communities. Simple term effects in the CCA model revealed significant (p < 0.05) differences in altitude, slope angle, slope ES, and wind speed variables. The present investigation sheds light on vegetation pattern and species contribution as a function of environmental gradients and provides a baseline for future studies.
Author Contributions: I.U.R., conceptualization, methodology, software, data curation, and writingoriginal draft preparation; A.A., supervision, conceptualization, and visualization; Z.I., supervision, conceptualization, methodology, and validation; E.S.C., data curation and helped in discussion; J.A., M.S.A. and N.A., manuscript editing, and discussion; R.K. and U.K., resources and helped in literature gathering; R.W.B., supervision, manuscript editing, and validation. All authors have read and agreed to the published version of the manuscript.

Funding:
The authors extend their appreciation to the researchers supporting project number (RSP-2021/193), King Saud University, Riyadh, Saudi Arabia.