Vegetation–environment interactions: plant species distribution and community assembly in mixed coniferous forests of Northwestern Himalayas

One of the main goals of ecological studies is to disentangle the dynamics that underlie the spatiotemporal distribution of biodiversity and further functions of the ecosystem. However, due to many ecological and geopolitical reasons, many remote areas with high plant species diversity have not been assessed using newly based analytical approaches for vegetation characterization. Here, we classified and characterized different vegetation types (i.e., major plant communities) based on indicator species and on the influence of different environmental gradients in the Himalayan mixed coniferous forest, Pakistan. For that, we addressed the following questions: Does the vegetation composition of the Himalayan mixed coniferous forest correlate with climatic, topographic, geographic, and edaphic variables? Is it possible to identify plant communities through indicator species in relation to environmental gradients using multivariate approaches? Can this multivariate be helpful for conservation planning? During four consecutive years we assessed the vegetation composition and environmental variables (21 variables divided in geographic, climatic, topographic, and edaphic groups) of 156 50 m-trasects between an elevation of 2000–4000 m. Using newly based analytical approaches for community characterization, we found a total of 218 plant species clustered into four plant communities with the influence of environmental gradients. The highest index of similarity was recorded between Pinus-Cedrus-Viburnum (PCV) and Viburnum-Pinus-Abies (VPA) communities, and the highest index of dissimilarity was recorded between PCV and Abies-Juniperus-Picea (AJP) communities. Among these four communities, highest number of plant species (156 species) was recorded in PCV, maximum alpha diversity (H’ = 3.68) was reported in VPA, highest Simpson index (0.961) and Pielou’s evenness (0.862) were reported in VPA and AJP. The edaphic gradients (i.e., organic matter, phosphorous, pH and soil texture) and climatic factors (temperature, humidity) were the strongest environmental gradients that were responsible for structuring and hosting the diverse plant communities in mixed coniferous forest. Finally, the Himalayan mixed coniferous structure is more influenced by the spatial turnover beta-diversity process (βsim) than by the species loss (nestedness-resultant, βsne). Our analysis of the vegetation structure along the environmental gradient in the Himalayan mixed coniferous forest supported by sophisticated analytical approaches reveled indicator species groups, which are associated to specific microclimatic zones (i.e., vegetation communities). Within this focus, we side with the view that these results can support conservation planning and management for similar and different areas providing mitigating and preventive measures to reduce potential negative impacts, such as anthropic and climatic.

Ecological studies not only help us comprehend the interaction between vegetation and the environment 1-3 , but they are also required for monitoring global climate change responses [4][5][6] .Nonetheless, such research identifies vegetation alteration pathways, causes, and processes 7 .Species diversity variation along environmental gradients is a core concern of ecological research 8 , and it has been explained in terms of climate, productivity, biotic interaction, habitat heterogeneity, and history 9,10 .Species diversity is a measure of resilience in ecosystems since it offers resistance to environmental changes 11,12 .Mountains often have a vast altitudinal range, abrupt climatic changes along the altitudinal gradient even over short distances, and a high level of endemism, making them more relevant for such research 13 .Species diversity is reduced as altitude rises due to temperature changes, precipitation, the length of the growth season, changes in solar radiation intensity, chilly, fast winds at high altitude, and steep slope aspect 14 .
Different factors influence the distribution and composition of plant communities in mountainous forests, such as altitude and closely linked edaphic 15 and climatic factors 16 , topographic heterogeneity 17 , soil chemistry 18,19 , species competition for nutrients 20 , soil texture 21 and light availability 22 .The environments for species growth and distribution are determined by a combination of these variables.Among these, elevation played a key role in influencing the diversity, richness, and distribution of species 23 .Through run-off redistribution, altitude also affects the availability of soil nutrients and water resources 24 .The moisture regime of dip and scarp landscapes, as well as concave and convex landscapes, differs often, as does the general flora.Run-off accumulates on a range of scales, from little depressions to enormous wades (run-on).As a result, niches and habitats of all types and sizes emerge, dictating the structure and composition of vegetation 25 .
The vegetation within a forest, on the other hand, is strongly affected by the local microclimate [26][27][28] .Forests have the highest species richness due to the presence of the herb layer 29,30 .By influencing resource availability and environmental variables important to herb layer plants, tree layer diversity can impact herb layer diversity 29,31 .While there have been reports of relationships between herb and tree layer diversity 32 , most research to date have investigated herb layer diversity between forest types in different regions of the Himalaya 33,34 with only a few dominant tree species or between different monospecific stands, in particular conifer versus broad-leaved forests 29 .
Regarding the effects of environmental gradients (e.g.mountains) on vegetation structure, sophisticated computer-based analyses, such as multivariate analytical programs are very useful to unravel the vegetation structure and dynamic 35 , identifying potential environmental drivers [36][37][38] and indicator species 35,39 .For instance, classification and ordination have been used as a way to summarize the multidimensional field data in a smaller number of dimensions clustering similar habitats and stands which share common species 36 .In addition, using these statistical analyses allow us to identity ecological groups, which can be clustered based on indicator values of different environmental drivers such geographic, climatic, topographic, and edaphic.
In this context, the aim of this study was to classify and characterize different vegetation types (i.e., major plant communities) based on indicator species and on the influence of different environmental gradients in the Himalayan mixed coniferous forest, Pakistan.Due to the remote areas with difficult access, uneven terrain, and adverse geopolitical relationships, most part of Himalayan forests has not been assessed using newly based analytical approaches for vegetation characterization.In other words, this study analyzed the vegetation structure and distribution, identifying potential indicator species groups for different microclimatic conditions along the mountainous gradient.Specifically, this study addressed the following questions: Does the vegetation composition of the Himalayan mixed coniferous forest correlate with climatic, topographic, geographic, and edaphic variables?Is it possible to identify plant communities through indicator species in relation to environmental gradients using multivariate approaches?Can this multivariate be helpful for conservation planning?Overall, The study was carried out in Manoor valley, which is a mountainous valley (34.68165N to 34.83869 N latitude, and 73.57520 E to 73.73182 E longitude) with 1580 to 4677 m elevation above sea level (Fig. 1) in the Himalayan belt of northwest Pakistan (for more details see Rahman et al. [40][41][42][43]46,49,50,122 ). The Himalaya are geologically young, having formed when the northward-drifting Deccan Plateau clashed with the Eurasian continent around 50 million years ago, causing geological upheavals that gave rise to the Himalaya, which now extends over 3000 km from Pakistan to Myanmar 44 .Whole area is defined by mountain ranges on both sides of the Manoor river, which runs northeast to southwest along the valley that emerges from Malika Parbat ('Queen of Mountains, ' elevation 5279 m).
The mountains are divided into three ranges, with this ecoregion located in the middle Himalayan range, which reaches a height of roughly 5000 m.Pure fir forest, mixed oak-fir forest, and mixed coniferous forests with fir, blue pine, and spruce are among the forest types found in the ecoregion 44 .The valley has been identified as an important region of the Himalayas with Sino-Japanese vegetation 45 , which is mainly composed of interconnecting mountain ranges that sustain monsoon-driven flora in a drier, colder climate 42,46 .This Himalayan part of Pakistan possesses a rich flora due to the fact that this part of country remained undisturbed by man for a long period; this enabled many species to survive and to evolve 47 .Even across short distances, the flora varies greatly 48 , with a high degree of resource seasonality 49,50,122 and diversity in both species 42,48 and communities 51 .

Vegetation sampling and plant identification
During four consecutive years (2015-2018), 13 sampling sites ranging from 2000 to 4000 m above sea level were selected in the mixed coniferous forest of Manoor valley, Himalaya, Pakistan.Line transect ecological technique was used for vegetation sampling 36 .In each sampling site, three transects of 50 m and 100 m apart each other were determined.In each year of evaluation, three different transects were designed to cover the largest vegetative area of the sampling site, totaling 12 transects per sampling site and 156 transects over the four years of analyses.For each sampling site, relative values of density, frequency, cover and Importance Value Index (IVI) were calculated following the methodology of Curtis and McIntosh 53 and Buckland et al. 54 .Each sampling site's GPS coordinates were recorded.The clinometer was used to identify the mountain's aspect, which included east (E), west (W), south (S), and north (N), as well as latitude, longitude, and height (GPS).Plant specimens were collected, tagged, placed between newspapers, pressed using a plant presser, poisoned with Mercuric Chloride and Ethyl Alcohol solution, and mounted on herbarium sheets 55,56  www.nature.com/scientificreports/specimens [57][58][59] .Scientific name of plant species was cross-checked and updated with an online website (www.thepl antli st.org) of the Royal Botanic Gardens, Kew (The Plant List, 2013), accessed on 13 November 2018.Voucher numbers were given to the plant species (Supplementary data Table S1), which were then deposited in the Hazara University Herbarium in Mansehra, Pakistan.

Ethical approval
This study was permitted and approved by the "Ethical committees of the Department of Botany" as well as "Advanced Studies and Research Board, Hazara University Mansehra, Pakistan".Experimental research and field studies on plants (either cultivated or wild), including the collection of plant material, comply with relevant institutional, national, and international guidelines and legislation.

Environmental gradients
In terms of edaphology, 200-400 g soil samples were taken from three random sampling site in each transect from a depth of 0-30 cm and properly mixed to produce a composite sample 60 .The samples were packed in polythene bags and labelled with a permanent marker.Furthermore, stones and other raw materials were sieved out before the remaining samples were dried in the shade.For each soil sample, physicochemical analyses were performed, including soil texture (clay, sand, silt, loam), pH 61 , electrical conductivity (EC) 62 , organic matters (OM) 63 , nitrogen (N), phosphorus (P), potassium (K), and calcium carbonate (CaCO3) concentrations [64][65][66] .With the use of a portable weather station (Kestrel weather tracker 4000), several environmental gradients (such as barometric pressure, dew point, humidity, heat index, temperature, wet bulb, and wind speed) were also determined.Here altitude is used as proxy for all the said environmental variables 67,68 .

Statistical analyses
All the gathered data of species, topographic as well as other environmental variables data were organized to determine the association between them 69,70 .The analyses were conducted using matrices of IV data from all studied sites.
Indicator species analysis TWINSPAN (two-way indicator species analysis) was used to identify plant communities and their primary indicator species using PC-ORD version 5.0. 36.This analysis applied Sorenson Distance Measurements using Wards Linkage Method 71 and IVI to identify trends of similarity 72 .

Diversity patterns
For each stand, the following diversity patterns were calculated: species richness, Pielou's evenness, Shannon (H'), and Simpson diversity indices.The H' value reveals not only how many species there are, but also how their abundance is distributed throughout all the species in the community.Pielou's evenness reveals how plant species are evenly distributed within a recognized community.Higher H' values represent maximum diversity.
The relative abundances of the most significant species are particularly sensitive to changes in Simpson's index.
Closer to 1 indicates clustering of individuals in a few species, whereas a small number (near to 0) indicates a more equal distribution of individuals across species.
To evaluate which is the beta-diversity component that most influence the vegetation distribution and structure in the area, we used the (1) spatial turnover (Simpson pairwise dissimilarity) and (2) nestedness-resultant components (nestedness-fraction of Sorensen pairwise dissimilarity) of β-diversity applying "Sorensen" as family of dissimilarity index [73][74][75] .Dissimilarity analysis and dendrograms were conducted in the package "betapart" 76 and "dendextend" 77 respectively.

Regression models
To compare the parameters (21 parameters distributed into the four different groups: geographic, edaphic, climatic, and slope) evaluated among the communities produced by our previously analyses, we conducted a Generalized Linear Model (GLM) with Gaussian error distribution followed by Likelihood ratio test using the packages "stats" and "car" 78 , respectively.In addition, after calculating species richness, Pielou's evenness, and Shannon and Simpson diversity for each stand of each plant community, we ran a GLM followed by Likelihood ratio test.In the case of species richness, we used a Poisson error distribution, since we have a count response variable; for others we used Gaussian error.

Ordinations
To depict the floristic relationships among the key syntaxonomic units (communities), non-multidimensional scaling ordination (NMDS) and Principal Component Analysis (PCA) were used 79 using the package "vegan" 80 in the software R 4.0.0 81 .
To check how explanatory factors (geographic, climatic, edaphic, and topographic) affect plant species distribution, we used canonical correspondence analysis (CCA) and variation partitioning tests (partial CCA) 82 .First, we used the step function with permutation in the "stats" package to create the optimal model with the fewest variables, those that best explain variance 81 .Next, we used the Variance Inflation Factor (VIF) to assess multicollinearity across variables in the final model, and we excluded any variables with VIF > 10 one by one.Finally, we used CCA and partial CCA on the final model to see how much each set of variables explained in our model.

Characterization and classification of vegetation
In the mixed coniferous forest of Manoor valley, Himalaya, Pakistan, 218 plant species were documented from 13 study sites (Supplementary data Table S1).Under the effect of several observed environmental factors, TWIN-SPAN (Twoway Indicator Species Analysis) categorization identifies four plant communities (Fig. 2).TWINSPAN and indicator species analyses discovered the following communities, which are discussed below.

Variation of environmental variables among communities
Most of the gradients were substantially different (p-value 0.05; Fig. 6) in GLM analyses of the 21 examined variables arranged in distinct classes (geographic, slope, climatic, and edaphic gradients) across the four plant communities.The only variables that did not indicate a significant difference were latitude, longitude, slope angle, humidity, dew point, pH, electric conductivity, calcium carbonate, phosphorus, and soil texture (i.e., sand, silt, and clay).

NMDS and PCA
The ordination diagrams such as NMDS and PCA were used to identify the relationship between mixed coniferous forest vegetation and environmental factors (Fig. 3a-e).In both ordinations, we can see that the environmental variables separate the sampling sites in the four communities already demonstrated and confirmed by TWINSPAN.In addition, most of the communities of mixed coniferous forest was resided by 20-50% of clay, sand, and silty soil texture (Fig. S1).

Discussion
It has long been the goal of ecological studies is to disentangle the dynamics that underlie the spatiotemporal distribution of biodiversity 83 , and further functions of the ecosystem 84,85 .Both biodiversity and ecosystem functions are driven by specific drivers of contemporary environments, i.e., biotic 86 and abiotic variables [87][88][89] .www.nature.com/scientificreports/Understanding how community composition varies in response to environmental variability is important in order to understand biodiversity 90 , productivity, and ecological stability 91,92 .With the advent of modern methods and techniques 93 , a number of contributions to modern plant ecology have addressed the question of complex vegetation patterns 70,[94][95][96] .TWINSPAN categorized 218 plant species from 13 Manoor Valley mixed coniferous forest study sites into four primary plant communities.When such formation of ecological linkages is based on indicator values 97,98 , each plant community often contains one or more indicator species 99,100 .As the study area lies in the Himalayan region, the vegetation predominantly exhibited characteristics of Sino-Japanese nature.The communities were classified within this region were categorized based on a range of environmental gradients i.e., soil pH, organic matter, electrical conductivity, nitrogen, potassium, phosphorous contents, soil texture, slope angle and aspect, and altitude etc.This allows our results to be compared with the studies already undertaken in other adjacent Himalayan regions 101 .Our findings are consistent with those of 102 , who identified four plant associations during an eco-floristic investigation of Beer Hills along the Indus River in Pakistan.Similarly 103 , documented seven plant communities and five major forest types while studying phytosociological analysis of Western Himalayan forests of Muzaffarabad, Azad and Jammu Kashmir, Pakistan.Nonetheless, the current study plant communities were structured by the influence of various ecological variables like geographic, slope, edaphic and climatic gradients.It also revealed that the indicator species of each plant community were linked to the particular set of environmental gradients.A region's forest communities evolve through time, but altitude, slope, latitude, aspect, precipitation, and humidity all have a part in their development and composition [104][105][106] .Ecosystems respond to several simultaneous changes in the environment 107 , which affect community diversity and distribution 94,108 .In each microhabitat type, specialist plant communities thrive, and are composed of specific taxa that have adapted to the unique environmental conditions of that particular microhabitat 109 .Aboveground 28 and underground communities work together to control whole-ecosystem processes and reactions to changes in the environment 110,111 .
The study area is richly diverse with Shannon Diversity index (H') values ranged between 2.89 and 3.68.Among recognized communities along the altitude and other gradients, VPA community has the maximum value (H' = 3.68) among all recorded communities followed by PCV community (3.62) and AJP (2.96).This bimodal response of diversity measures might be due to the anthropogenic pressure 112 .The higher anthropogenic activities probably be linked with the easy access to the such sites 113 .Furthermore, most of the communities of mixed coniferous forest were recorded on the sloppy surfaces and edges of the mountains that might be the reason of their diversity.Gehlhausen et al. 114 observed the forest edges with maximum diversity as compared to the forest interiors.Similar pattern was reported by Khan 51 from the neighbor valley (Naran Valley, Northwestern Himalaya).Such diverse species also reveals their wider ecological amplitude.
In the topographic class, we showed that altitudinal gradients offer an important range of different environmental variables, highlighting the existence of micro-climates that drive the structure and composition of plant species in each micro-region.We compared our results with other vegetation surveys carried out in allied 51 and neighboring Himalayan regions [115][116][117] which similarly described the influential role of altitude.As a result, the northern slopes sustained the denser and homogenous growth of conifers in Manoor Valley, Himalayan forests, Pakistan.In addition, the scores of variation partitioning indicated that topography was the main driver in mixed coniferous forests, whereas those results were found according to the findings of 118 .
The ordination study demonstrated that varied edaphic gradients had a substantial impact on the mixed coniferous forests, which were categorized into four primary plant communities.The structure of communities and plant growth, ground cover, and natural regeneration capabilities are all affected by edaphic gradients 119 .Likewise, the most influential contributing gradients were soil texture (loam, silty loam, clay and clay loamy), organic matter and potassium.In the present study, the loamy soil with high organic matter significantly contributed in structuring VPA community.Similarly, Dvorský et al. 116 reported the contribution of soil moisture as the influential variable for shaping the community structure.This study depicted that the vegetation of mixed coniferous forests was strongly influenced by soil texture gradient and their contribution assessment was visualized by the direct ordination approaches.For instance, the PCV community was recorded under the influence of silty loam soil texture.Nonetheless, studies of 115,116 have reported similar findings, they also concluded that soil texture was a strong contributing gradient in structuring the vegetation from the same terrain of Himalayas.
The results showed non-significant differences among communities in relation to electric conductivity, pH and phosphorous, but slight variations were noticed in their average means.At the higher altitudinal sites, maximum humification was recorded and that might be due to the higher cover which escorted to a drop-in soil pH.The lowest soil pH at higher elevations may further led to an increase in phosphorus content through mineralization.Current results were compared with the previous reports on Himachal Pradesh (Northwestern Himalaya), India 117 , which revealed similar findings by stating non-significant impact of soil pH and phosphorus on the vegetation due to high elevation.Lastly, the possible reasons of all these resemblances may be just because of the matched environmental conditions hosted by the mountains of Himalayan region that in turn governs the biotic and abiotic variables.
Finally, beta-diversity based on the turnover of species is the trait that most influence the distribution of plant species in the Himalayan mixed coniferous forest of Pakistan.In other words, instead of decreasing the number of species along the altitudinal gradient (used as a proxy for climatic and other environmental gradients 120,121 ) and under different climatic conditions, there is a species turnover, i.e., plants that live in high elevations might have different traits than those living in low elevations 122,123 .This variation might have allowed them to survive in this harsh or different environments, leading to a variation in plant community along the altitude.Changes in the environment alter the diversity and organization of plant communities by altering the spectrum of species features that may be effective in new environments 124 .
In addition to geographic, topographic, and edaphic gradients, climatic gradients also represented a vital role in hosting the major plant communities of mixed coniferous forests.The most significant contributing gradients were temperature, heat index, wind speed, wet bulb, and barometric pressure.The analytical approaches revealed the positive and significant correlation of PCV community with temperature, heat index and barometric pressure.This might be due to the region hosting this community were located at the lower altitudinal ranges as compared to other major groups.As we all know, the response of vegetation structure to changes in environmental gradients has a significant impact on its development.Among all the recorded plant communities, PCV community was the dominant one based on the number of associated plant species (156 species).This variation in the number of associated plant species within communities might be due to variability in the values of edaphic and other environmental gradients 125,126 which are responsible in sustaining the growth of various associated species 127 .Plant communities can be described in a way that assists management decisions for a variety of ecological communities 128 .

Conclusions
Multivariate analyses categorized 218 plant species found at 13 study sites in the mixed coniferous forest into four distinct plant communities, each with its own indicator species.Among all the recorded plant communities, PCV community was the dominant one based on the number of associated plant species (156 species).The impact of numerous ecological factors shaped these plant assemblages.In the topographic class, altitude was shown to be the most important gradient, followed by latitude and longitude.The northern slopes fostered the denser and more uniform development of conifers in Manoor Valley, Himalayan forests, Pakistan.In addition, the scores of variation partitioning indicated that topography was the main driver in mixed coniferous forests.The most influential contributing edaphic gradients were soil texture (loam, silty loam, clay and clay loamy), organic matter and potassium.Likewise, the most significant contributing environmental gradients in structuring and hosting the plant communities of mixed coniferous forest were temperature, heat index, wind speed, wet bulb, and barometric pressure.The topmost indicator species and other associates of VPA community was hosted by loamy soil with higher sandy texture, calcium carbonate and OM as compared to other communities.These indicator species could be used to observe changes in plant communities as a result of changes in the environment or management.As a result, recognizing such an indication might be used to manage species in a range of microhabitats with varying soil types and climatic conditions.
Assessing the abiotic and biotic variables that drive the ecosystem dynamics is one of the main goals nowadays, mainly due to the continuous process of climate change and anthropogenic impacts.Studies like this can help in understanding the structure of plant communities, observing how each community responds to a certain environmental change.In this context, it is possible to identify (1) indicator species, (2) anthropic impacts, and (3) climate and soil changes in certain environments, or even provide mitigating and preventive measures to reduce these impacts.Finally, the Himalayas, where this study was developed, is a highly biodiverse region with high endemism and degradation, being classified as a biodiversity hotspot 129 .Studies evaluating the structure of plant communities residing in these locations are therefore essential for maintaining biodiversity.

Figure 3 .
Figure 3. Non-Multidimensional Scaling (NMDS) among plant communities and environmental gradients: (a) geographic, (b) slope, (c) edaphic and (d) climatic.(e)The association between several measurable environmental factors and the communities depicted in coloured circles using Principle Component Analysis (PCA).Based on a 95% confidence level, each circle with a distinctive colour represents a distinct community.PCV, Pinus-Cedrus-Viburnum; VPA, Viburnum-Pinus-Abies; AJP, Abies-Juniperus-Picea; PAJ, Picea-Abies-Juniperus.The length of the arrows indicates the effect and strength of each environmental gradient, while the direction of the arrows indicates the correlation of each environmental gradient.Positive correlation was found for gradients on the same axis, whereas negative correlation was found for gradients on opposing axes.Big circle in (e) demonstrates the centroid of each plant community.

Figure 6 .
Figure 6.Boxplots showing the variations of the studied variables of the four plant communities evaluated in this study (GLM result, and its associated p-values are displayed at each whisker boxplot).Plant communities in x-axis were plotted in ascending order according to elevation gradient (low to high elevation).PCV, Pinus-Cedrus-Viburnum; VPA, Viburnum-Pinus-Abies; AJP, Abies-Juniperus-Picea; PAJ, Picea-Abies-Juniperus.

Figure 7 .
Figure 7. Variation partitioning findings (partial CCA) and contribution (percent) of the four variable groups analysed are shown in a Venn diagram.Negative values signify zeros, indicating that the explanatory factors explain less variance than random normal variables 82 .

Table 1 .
The contribution and ranking of the environmental variables used in our CCA model.

Table 2 .
Results of variation partitioning (partial CCA) of four variable groups studied (see Fig.7for individual fraction letters code).