Species Distribution Pattern and Their Contribution in Plant Community Assembly in Response to Ecological Gradients of the Ecotonal Zone in the Himalayan Region

The ecotonal zones support populations that are acclimated to changing, fluctuating, and unstable conditions, and as a result, these populations are better equipped to adjust to expected change. In this context, a hypothesis was tested that there must be vegetation dominated by unique indicator plant species under the influence of ecological gradients in the ecotonal zone of Manoor Valley (northwestern Himalaya), Pakistan. Keeping the aforementioned hypothesis in mind, detailed field studies were conducted during different seasons in 2015-18. Line transect sampling and phytosociological characteristics (density, frequency, cover, and their relative values and Importance Value) were implemented as ecological methods. This investigation documented 97 plant species recorded from seven sampling sites. The community distribution modelling revealed that the ecological variables separate the seven sampling sites into two major plant communities (Indigofera-Parrotiopsis-Bistorta and Ziziphus-Leptopus-Quercus) recognized by TWINSPAN. The IBP communities showed a positive and significant correlation with altitude (1789.6–1896.3 m), sandy soil texture with a slightly acidic pH (6.4–6.5), and higher phosphorous (9–13 mg kg−1). In contrast with this, the ZLQ community was recognized on the southern slope under the strong influence of high electrical conductivity (2.82–5.4 dsm−1), organic matter (1.08–1.25%), calcium carbonate (5.8–7.6 mg kg−1), potassium (202–220 mg kg−1), and temperature (28.8–31.8 °C). Hence, both communities were found on opposite axes with clear differences based on the ecological gradients. NMDS clustered different species with similar habitats and different stands with common species, showing that plant species and stands were in a linear combination with ecological gradients. The IPB community has the maximum number of plant species (87 species), Shannon value (H’ = 4), Simpson value (0.98), and Pielou’s evenness value (0.96). Thus, the multivariate approaches revealed unique vegetation with sharp boundaries between communities which might be due to abrupt environmental changes.


Introduction
The term "ecotone" was first coined and used for the transition zone by Frederic E. Clements in 1905 [1] and is considered as the basic unit of landscape ecology [2], having extremely unique vegetation [3], sensitive ecosystems, and sharp plant community boundaries [4], due to abrupt environmental changes [5]. The natural layering of ecosystems that occurs at different elevations owing to differing environmental conditions is referred to as altitudinal zonation in mountainous locations [6]. Nonetheless, a wide range of environmental variables, ranging from direct effects of temperature and precipitation to indirect characteristics of the mountain itself, as well as biological interactions among species, determine the boundaries of the altitudinal zones found on mountains, resulting in discrete plant communities [7][8][9]. They are natural assemblages of co-evolved individual populations, which form visible units [10,11]. Ecotones often appear along ecological gradients [12]. Such gradients are created because of spatial shifts in elevation, climate, soil, and many other ecological factors [13]. These regions are sometimes regarded as dynamic zones of community interaction that are inherently unstable over time [14]. Ecotones, as indicated by Odum [15], are more than just a boundary or an edge; the idea of an ecotone presupposes the presence of active interaction between two or more ecosystems, with features that are not present in either of the neighboring ecosystems.
Ecotones have been investigated from an ecological perspective for the past four decades [13,[16][17][18] and have recently received considerable attention in the context of biodiversity protection. Ecotones are "natural laboratories" where researchers may study a range of evolutionary processes, such as speciation and the emergence of new species. It is greatly influenced by various environmental gradients, i.e., edaphic, climatic, and physiographic variables [19]. In those regions, environmental conditions that fluctuate over time and space might be useful indicators of environmental change and ecological responses to climate change. Some researchers claim that ecotonal zones support populations that are acclimated to changing, fluctuating, and unstable conditions, and that as a result, these populations are better equipped to adjust to expected change. Ecotones' efficacy as early predictors of climatic changes and the ways that ecological groups and systems adapt to change might be a promising future endeavor. Ecotones can support distinct species or indicators that are less common or do not exist elsewhere. Therefore, a hypothesis was tested that there must be vegetation dominated by unique indicator plant species in the ecotonal zone. The current study was planned to determine the species distribution pattern, species contribution, and major plant community assembly in association with the ecological gradients of the said zone.

Study Area
The current project lies in the ecotonal zone of the Manoor Valley. Manoor Valley is reached from the main Kaghan Valley road at the junction 'Mahandri' [20] and is about 50 km north of Balakot [21,22]. Geographically, Manoor Valley ( Figure 1) is situated in northwestern Pakistan (34.68165 N to 34.83869 N latitude and 73.57520 E to 73.73182 E longitude) as a part of the mountainous series collectively known as the Himalayas [23]. The entire area is formed by crosswise ridges of mountains on either side of the Manoor river, which flows in a northeast to southwest direction down the valley emerged from Malika Parbat ('Queen of Mountains', elevation 5279 m).

Vegetation Sampling
Detailed field visits were conducted to the ecotonal/transition zone of Manoor Valley, Himalaya, Pakistan. Line transect ecological technique was used for vegetation sampling [24][25][26][27][28][29]. The study area was subdivided into seven sampling sites, and three transects each of 50 m at each sampling site were randomly selected and sampled [30]. The interval distance between the transects were 100 m and 200 m interval between each sampling site. The individuals of plant species falling precisely on the line were recorded. The phytosociological parameters such as density, frequency, cover, and their relative values and Importance Value (IV) were calculated following the methodology of Curtis and McIntosh [31] and Buckland et al. [25]. Plant specimens were collected, serially tagged with a field number, pressed between the blotting papers for drying, poisoned, and mounted on standard herbarium sheets following the recommended procedures [32,33]. All the

Vegetation Sampling
Detailed field visits were conducted to the ecotonal/transition zone of Manoor Valley, Himalaya, Pakistan. Line transect ecological technique was used for vegetation sampling [24][25][26][27][28][29]. The study area was subdivided into seven sampling sites, and three transects each of 50 m at each sampling site were randomly selected and sampled [30]. The interval distance between the transects were 100 m and 200 m interval between each sampling site. The individuals of plant species falling precisely on the line were recorded. The phytosociological parameters such as density, frequency, cover, and their relative values and Importance Value (IV) were calculated following the methodology of Curtis and McIntosh [31] and Buckland et al. [25]. Plant specimens were collected, serially tagged with a field number, pressed between the blotting papers for drying, poisoned, and mounted on standard herbarium sheets following the recommended procedures [32,33]. All the specimens were identified with the help of Flora of Pakistan and other available litera-  [34][35][36]. Afterwards, the plant specimens were deposited in the Herbarium of Hazara University, Mansehra, Pakistan for accession numbers.

Ecological Gradients
GPS locations were recorded for each sampling site. Aspect of the mountain, i.e., east (E), west (W), south (S), and north (N), was determined with the help of clinometer, and latitude, longitude, and altitude using geographical positioning system (GPS). Regarding edaphology, soil samples of 200-400 g were collected from three random points within each sampling site from depth of 0-30 cm [37] and mixed thoroughly to make a composite sample [38]. Samples were placed in polythene bags and properly tagged with permanent marker. Moreover, rocks and other raw materials were removed by sieving and then the remaining samples were shade dried. Physicochemical analyses, i.e., soil texture (clay, sand, silt, loam), pH [39], electrical conductivity (EC) [40], organic matters (OM) [41], nitrogen (N), phosphorus (P), potassium (K), and calcium carbonate (CaCO 3 ) concentrations were determined for each soil sample [42][43][44]. Moreover, other climatic gradients (i.e., barometric pressure, dew point, humidity, heat index, temperature, wet bulb, and wind speed) were determined with aid of handheld weather station (Kestrel weather tracker 4000).

Data Analyses
All the data of plant species and ecological gradients recorded in the field were transferred to the MS-Excel spreadsheets and arranged to determine the relationship among them [45,46]. The analyses were conducted using matrices of IV data from all of the recorded plant species (97 species) and ecological gradients from seven sampling sites (Tables S1 and S2) using two separate spreadsheets like IV of species x sampling sites and ecological gradients data x sampling sites [19,47]. A georeferenced map was prepared using ArcGIS version 10.1 to show the study area. For multivariate ordination analyses, CANOCO 5 [45], PC-ORD [48], and R 3.6.1 [49] were used. Species area curve (SAC) was drawn to evaluate the adequacy level of area sampled in relation to the vegetation [19]. Cluster analysis (CA), two-way cluster analysis (TWCA) [50] and two way indicator species analysis (TWINSPAN) were carried out using PC-ORD software [51] for the identification of major groups [19]. Moreover, based on species categorization and sample clustering [52], TWINSPAN was processed for the recognition of the plant communities/indicator species [51].
Non-multidimensional scaling ordination (NMDS) and Principal Component Analysis (PCA) were done to visualize the floristic associations among the main taxonomic units (communities) [53] using the package "vegan" [54]. Analysis of Similarities (ANOSIM p < 0.001, 5039 permutations) was carried out to determine the significant variations in recorded plant species composition between communities. The ANOSIM was processed in R 3.5.1 software using the package "vegan" [54]. To compare the parameters (diversity and ecological) evaluated among the two communities, we conducted a GLM with Gaussian error distribution (except for species richness, in which we used Poisson distribution) followed by Likelihood ratio test using the 'stats' and 'car' [55] packages, respectively. In addition, we measured Pielou's evenness, Shannon, and Simpson diversity indices for each stand of each plant community using Gaussian error [19,56].

Results
A very narrow range of Himalayan Subtropical-Temperate Ecotone was observed mostly on the southern slopes and foothills of the Manoor river. In total, 97 species were recorded from seven stands.

Vegetation Characterization
The species area curve (SAC) analysis showed that the maximum number of plant species appeared in the sixth stand and the species curve became parallel after that, as no new species were recorded later, which clearly shows that the sampled area was enough for the targeted vegetation ( Figure 2a). According to TWINSPAN (Figure 2b), CA (Figure 3), and TWCA ( Figure 4), a total of two different major plant communities (Indigofera-Parrotiopsis-Bistorta and Ziziphus-Leptopus-Quercus) were recognized by clustering all the recorded plant species (97 species) in seven stands by the influence of ecological gradients. They ranged from 1780.8 m to 1896.3 m. Two different clusters were observed by the TWINSPAN method, which showed a high cluster heterogeneity value (Lambda = 0.5801). Each community was comprised of different indicator species (i.e., Indigofera heterantha, Parrotiopsis jacquemontiana, Bistorta amplexicaulis and Ziziphus sp., Leptopus chinensis, Quercus incana) and was recorded at different altitudes ( Figure 2b).

Vegetation Characterization
The species area curve (SAC) analysis showed that the maximum number of plant species appeared in the sixth stand and the species curve became parallel after that, as no new species were recorded later, which clearly shows that the sampled area was enough for the targeted vegetation ( Figure 2a). According to TWINSPAN (Figure 2b), CA ( Figure  3), and TWCA ( Figure 4), a total of two different major plant communities (Indigofera-Parrotiopsis-Bistorta and Ziziphus-Leptopus-Quercus) were recognized by clustering all the recorded plant species (97 species) in seven stands by the influence of ecological gradients. They ranged from 1780.8 m to 1896.3 m. Two different clusters were observed by the TWINSPAN method, which showed a high cluster heterogeneity value (Lambda = 0.5801). Each community was comprised of different indicator species (i.e., Indigofera heterantha, Parrotiopsis jacquemontiana, Bistorta amplexicaulis and Ziziphus sp., Leptopus chinensis, Quercus incana) and was recorded at different altitudes ( Figure 2b).

Vegetation Characterization
The species area curve (SAC) analysis showed that the maximum number of plant species appeared in the sixth stand and the species curve became parallel after that, as no new species were recorded later, which clearly shows that the sampled area was enough for the targeted vegetation ( Figure

Vegetation along the Ecological Gradients
Correlations between Himalayan Subtropical-Temperate Ecotone plant communities and ecological gradients were illustrated through NMDS (Figure 5a-d) and PCA ( Figure  5e). Communities are related to specific ecological variables such as geographic, slope, edaphic, and climatic variables. Altitude, slope NE, slope N, K, pH, OM, silt, sand, clay, temperature, heat index, wind speed, and altitudinal density were the most representative

Vegetation along the Ecological Gradients
Correlations between Himalayan Subtropical-Temperate Ecotone plant communities and ecological gradients were illustrated through NMDS (Figure 5a-d) and PCA (Figure 5e). Communities are related to specific ecological variables such as geographic, slope, edaphic, and climatic variables. Altitude, slope NE, slope N, K, pH, OM, silt, sand, clay, temperature, heat index, wind speed, and altitudinal density were the most representative ecological variables which drive the plant community structure and diversity along the altitudinal gradient. In both ordinations, the ecological variables separate the seven sampling sites into two major plant communities already demonstrated by TWINSPAN (Figure 2b). In constrained PCA ordination, the maximum explanatory variation was accounted for PC1 axis (34.3%) and lower variation on the PC2 axis (33.7%). Maximum strength was recorded for the most important ecological gradients, i.e., the profound influence of the altitudinal gradient was revealed by dividing the vegetation of the Himalayan Subtropical-Temperate Ecotone into two communities (Figure 2b-d). The IBP communities showed a positive and significant correlation with altitude (1789.6-1896.3 m), sandy soil texture with a slightly acidic pH (6.4-6.5), and higher phosphorous (9-13 mg kg −1 ). In contrast to this, the ZLQ community was recognized on the southern slope. The edaphic gradients structuring the indicators and associated plant species of this community were high electrical conductivity (2.82-5.4 dsm −1 ), organic matter (1.08-1.25%), calcium carbonate (5.8-7.6 mg kg −1 ), and potassium (202-220 mg kg −1 ). Temperature (28.8-31.8 • C), barometric pressure (815.6-816.4 pa), and heat index (29.8-34.3) were found to have a positive and significant influence when compared to the IPB community. As a result, both communities were discovered on opposing axes, with distinct differences based on ecological gradients.  Table 1. IPB: Indigofera-Parrotiopsis-Bistorta and ZLQ: Ziziphus-Letopus-Quercus.   Table 1. IPB: Indigofera-Parrotiopsis-Bistorta and ZLQ: Ziziphus-Letopus-Quercus.

Analysis of Similarity (ANOSIM)
The ANOSIM was conducted to measure the species contribution for the site distribution pattern. Significant variation in plant species composition between communities was found (ANOSIM p < 0.001, 5039 permutations). Out of 97 species, 16 were significant and greatly contributed to the variation in plant community ordination in NMDS (Table 1).

Significance Testing of Plant Communities in Relation to Studied Variables
GLM analyses of the 20 studied variables distributed in different groups (geographic, slope, climatic and edaphic gradients) between two plant communities showed that most of the gradients were significantly different (p-value < 0.05; Figure 6). Nonetheless, slope angle, humidity, dew point, wind speed, wet bulb, organic matter, soil texture, potassium, and calcium carbonate did not show significant differences ( Figure 6).

Species Richness and Diversity Indices
Species richness values ranged from 30 and 87 species in two plant communities (Figure 7a). The highest number of plant species was reported in the IPB community (87 species) at 1789.6-1896.3 m altitude, followed by ZLQ (30 species) at the altitudinal range of 1780.8-1789.3m. The Shannon Diversity index values ranged between 2.4 and 4 ( Figure  7b). The IPB community had the maximum value (H' = 4), followed by the ZLQ community (2.4). The Simpson's dominance index values ranged between 0.86-0.98. The IPB community had the maximum value (0.98) between the recorded communities. Moreover, the lowest Simpson dominance value has been calculated for the ZLQ community (Figure 7c). The Pielou's evenness index values ranged between 0.75 and 0.96 (Figure 7d). The IPB community had the maximum value (0.96); on the other hand, the ZLQ community had the lowest Pielous's evenness value.

Species Richness and Diversity Indices
Species richness values ranged from 30 and 87 species in two plant communities (Figure 7a). The highest number of plant species was reported in the IPB community

Discussion
Mountains are the most remarkable land forms on the earth's surface with major vegetation zones based upon environmental variations [57]. Ecotones commonly coincide with areas of abrupt climatic transition along ecological gradients [16]. They can be found at a variety of spatial dimensions [58], ranging from continental-scale biome transitions to small-scale ecotones where plant communities and microhabitats coincide [12]. A very narrow range of Himalayan Subtropical-Temperate Ecotone was observed mostly on the southern slopes and foothills of the Manoor river. Aspect encourages habitat heterogeneity and induces micro-environmental variation in vegetation patterns [59,60]. The community distribution pattern (Figure 5a-e) revealed that ecological variables divided the seven sampling sites and their 97 native species into two major plant communities, identified by CA, TWCA, and TWINSPAN classifications. Ordination techniques have been widely used to reveal plant species diversity, distribution, and community recognition along ecological gradients [61,62]. Furthermore, multivariate statistical approaches can assist ecologists effectively in describing and elaborating the structure of vegetation, as well as quantifying the effects of ecological gradients on vegetation and/or a group of species [63,64]. Statistical tools can provide an efficient means of reducing the complexity inherent in natural vegetation and of detecting important environmental factors that explain this complexity [48,65].
Indigofera-Parrotiopsis-Bistorta and Ziziphus-Leptopus-Quercus were identified as the major plant communities (ranges from 1780.8 m to 1896.3m) by grouping all the recorded

Discussion
Mountains are the most remarkable land forms on the earth's surface with major vegetation zones based upon environmental variations [57]. Ecotones commonly coincide with areas of abrupt climatic transition along ecological gradients [16]. They can be found at a variety of spatial dimensions [58], ranging from continental-scale biome transitions to small-scale ecotones where plant communities and microhabitats coincide [12]. A very narrow range of Himalayan Subtropical-Temperate Ecotone was observed mostly on the southern slopes and foothills of the Manoor river. Aspect encourages habitat heterogeneity and induces micro-environmental variation in vegetation patterns [59,60]. The community distribution pattern (Figure 5a-e) revealed that ecological variables divided the seven sampling sites and their 97 native species into two major plant communities, identified by CA, TWCA, and TWINSPAN classifications. Ordination techniques have been widely used to reveal plant species diversity, distribution, and community recognition along ecological gradients [61,62]. Furthermore, multivariate statistical approaches can assist ecologists effectively in describing and elaborating the structure of vegetation, as well as quantifying the effects of ecological gradients on vegetation and/or a group of species [63,64]. Statistical tools can provide an efficient means of reducing the complexity inherent in natural vegetation and of detecting important environmental factors that explain this complexity [48,65].
Indigofera-Parrotiopsis-Bistorta and Ziziphus-Leptopus-Quercus were identified as the major plant communities (ranges from 1780.8 m to 1896.3m) by grouping all the recorded plant species (97 species) and the seven stands were strongly influenced by ecological gradients. Plant communities that can be described both physiognomically and floristically have a well-defined structure in relation to abiotic and biotic factors [31,66,67]. Using more or less the same techniques, a similar pattern of species distribution was reported from another mountainous area [68]. The evergreen Quercus incana has also been reported as an indicator species in the subtropical-temperate ecotonal zone of the Nandiar catchment, Batagram [69], as have Bisorta amplexicaule and Indigofera heterantha at a lower xeric elevation in the ecotonal range of District Khurram, Pakistan [18]. Such similarities clearly reveals that some abiotic factors have the influence on the distribution of vegetation [19]. Our results revealed that the IBP community had been strongly influenced by altitude (1789.6-1896.3 m), sandy soil texture with slightly acidic pH (6.4-6.5) and higher phosphorous (9-13 mg kg −1 ). In contrast with this, the ZLQ community was recognized on the southern slope under the strong influence of high electrical conductivity (2.82-5.4 dsm −1 ), organic matter (1.08-1.25%), calcium carbonate (5.8-7.6 mg kg −1 ), potassium (202-220 mg kg −1 ), and temperature (28.8-31.8 • C). Thus, both the communities were clustered with different leading indicator species, and similar habitats and different stands with common species, showing that plant species and stands were in a linear combination with ecological gradients. Similarly, many other researchers [70] stated that aspect, altitude, and soil depth are the most influential ecological variables for determining the community structure in the Naran Valley.
The highest number of plant species were reported in the IPB community (87 species) at 1789.6-1896.3 m altitude. The IPB community had the maximum Shannon value (H' = 4), Simpson value (0.98), and Pielou's evenness value (0.96). The distribution of species richness along the altitudinal gradient is governed by series of interacting biological and climatic factors [71][72][73]. The gradual decrease in species richness along with increasing altitude is considered as a general pattern [74]. In such studies, species richness [75,76] and abundances tend to peak in ecotonal regions, according to several investigations, however outliers exist. Plant species are distributed in different type of habitats; but in their own territory, they usually show an abundance that represents their special ecological optimum [57]. Due to this reason, the composition of distinct units is regarded as a function of changing habitat conditions along ecological gradients [77].

Conclusions
The current study provides the baseline and first insights into spatial distribution and vegetation mapping in response to ecological gradients of the ecotonal zone in Manoor Valley, Himalaya, Pakistan. The community distribution pattern revealed that the ecological variables separate the seven sampling sites into two major plant communities. Both the communities were clustered with different leading indicator species, sharing a similar habitat and different stands with common species, showing that plant species and stands were in a linear combination with ecological gradients. Significant variation in plant species composition between communities was recorded (ANOSIM p < 0.001, 5039 permutations). Out of 97 species, 16 were found significant that greatly contributed to the variation in plant community ordination. Hence, the multivariate approaches revealed unique vegetation with sharp boundaries between communities, which might be due to abrupt environmental changes. In the context of global change, including land-use and climate change, research on the influence of ecotones on biodiversity is an essential future direction. Since climate change is projected to be quick and intense in ecosystem border zones [78] and these regions might possibly act as "early warning" indicators of global changes by studying variations in ecotone locations over time [79,80]. The response, however, is dependent on the geographical and temporal scales studied, and it may be a helpful indication primarily at global spatial scales and on fairly coarse durations.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/plants10112372/s1, Table S1. List of plant species with their Importance values (IV) calculated based on vegetative characteristics of each sampling sites. Table S2. Means of environmental variables measured at three different transects of each sampling site were recorded (wind speed averages are given as integers).