{'fraction_non_alphanumeric': 0.03538627768773636, 'fraction_numerical': 0.021981361426256077, 'mean_word_length': 4.828970675063964, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
HOW ENVIRONMENTAL VARIABLES CAN DETERMINE THE CHIR PINE (PINUS ROXBURGHII SARG.) DISTRIBUTION IN SWAT HINDUKUSH RANGE OF PAKISTAN: CURRENT AND FUTURE PROSPECTIVE OF THE SPECIES

Distribution patterns of plant species and their relationships with physiographic, soil and climatic variables were investigated in the subtropical pine forests across the natural limits in the Swat, Hindukush range of Pakistan. A vegetation survey identified 9 tree species belonging to 8 families and 8 genera, which were classified into 3 distinct vegetation communities with an average density ranging from 512-1231 individuals ha and basal areas of 31.6-36 m ha. Total density and basal area values for tree species located at higher elevations were lower than those in the middle and low elevation communities for all the species. We found that P. roxburghii had unimodal size distributions, suggesting that these populations are not at equilibrium and are changing over time. CCA-Ordination with the associated Monte Carlo permutation test was employed to explore the patterns of variation in vegetation distribution and identified elevation and soil organic matter as the most influential variables responsible for the changes in species composition. The MaxEnt predictive modeling results clearly indicate a significant difference in the present predicated inter-site floristic composition and distribution of the species, suggesting an overall restriction to the south, whereas, the future distribution model suggest a shift towards the neighbouring districts in the west. This study identified few significantly important environmental variables linked with P. roxburghii associations and their distribution across the study area. We recommend additional research that includes multiple sampling from varied locations and other abiotic and biotic variables throughout Pakistan and neighbouring countries for better understanding of the species niche modeling.


Introduction
Disentangling of the principal mechanisms responsible for structuring plant communities has been a central research dilemma in community ecology (McCune and Grace, 2002; Condit et al., 2011). Several factors such as biogeochemical variations and edaphic, topographic, soil, and climatic variables allow for a suite of potential limiting factors and have varied effects on the vegetation structure and function at both local and regional scales (Song et al., 2004;Miehe et al., 2009). However, the effects of such factors on vegetation distribution were not examined until the advent of multivariate statistical techniques and bioclimatic models for the forest tree species in the northern mountainous ranges of Pakistan (Khan et al., 2011). The prominent role of environmental in Swat Hindukush range of Pakistan. These forests are distributed on both sides of river Swat and in the inner valleys, classified under the subtropical dry temperate forests of the country (Champion et al., 1965) and typically occurring in the Sino-Japanese phytogeographical region of the world (Sher et al., 2014). Although, P. roxburghii is commercially and ecologically important and has high conservation value, but so far remains little explored by foresters and biologists .
In the present study, an attempt was made to explore P. roxburghii associations and their future distribution with environmental and climatic variables using advanced multivariate statistical methods and species distribution models (SDMs). Among the SDMs, MaxEnt model reliably predicts suitable habitat using presence records and pseudo-absence points Elith et al., 2006) and have been widely used for the identification of suitable habitats, potential distribution range, and plant's future distribution changes (e.g., Ashraf  Kamyo and Asanok, 2020; Khan et al., 2020). Using both these techniques our specific research objectives were to examine (i) whether the occurrence of P. roxburghii is restricted to certain compositional and environmental circumstances (ii) whether P. roxburghii is regenerating in various community types linked with certain topographic, edaphic and climatic variables by analyzing its stand structures (iii) that density and recruitment of P. roxburghii in the forests can vary the structural characteristics of the stand, which may reflect an environmental gradient from various sites strongly dominated by P. roxburghii to those rich in broadleaved species and (iv) to simulate and compare the possible present and potential future distribution of P. roxburghii associations in Swat Hindukush range of Pakistan to assist its conservation and management.

Study area
The Swat valley is part of the subtropical dry temperate areas in Hindukush range of northern Pakistan, located at 34° 34ʹ to 35° 55ʹ N and 72° 08ʹ to 72° 50ʹ E (Shinwari et al., 2003). The total area of the district is 5337 km 2 with long stretches of pines, and it is generally considered as a hub to biodiversity due to its unique hills and climate . Geographically, the area shares borders with Chitral, Indus Kohistan, Shangla, Bunir and District Dir (Fig. 1). The area has a Mediterranean climate comprising two phytoclimatic subtypes: dry and moist (Champion et al., 1965). The Swat meteorological station data revealed that the average annual temperature is 19 °C and annual precipitation averages 897 mm. Mean relative humidity remains high from January to March (75%) and dips below 40% from May to June. The most abundant substrate consists of marble sandstones, which form soils of the sandy type only, and are accompanied by clays, marls and limestone in Vertisols . P. roxburghii stands constitute the most extensive plant formation in the area between altitudinal ranges from 700 m to 1800 m above sea level. P. roxburghii occurs either in pure stands or mixed with other species (Appendix 1) in the shady zones and the valley bottoms, as well as interspersed with Olea ferruginea and Quercus baloot in drier transition zones.

Field methods
After general reconnaissance of the Swat District, 25 sites were selected for sampling that approximately covered the entire range of P. roxburghii distribution in its natural Khan et  Prior to field data collection, stand locations were selected systematically using Arc GIS v. 9. Coordinate pairs and other geophysical characteristics were downloaded as waypoints into GPS receivers so that the forest stands could be located in the field. In addition, collection of geo-referenced data about the stands was ensured robust with RedHen DX-GPS system and Nikon D300 camera to gather and save the background information with pictures as metadata. Over 2,000 photographs of forest communities were obtained and metadata was extracted with the help of BR's EXIF extractor, a freeware available online (http://www.br-softwaare.com/extracter.hotml).

Figure 1. Map showing altitudinal topographical variations in Swat District Hindukush mountain range of Pakistan
At each forest stand, 30 quadrates of 10 m × 10 m in size along a 200 meter straight transect in a suitable direction in both disturbed and undisturbed forests were phytosociologically analyzed (Uprety et al., 2014). All constituent trees of ≥ 10 cm diameter at breast height (DBH) were counted and DBH (above 1.37 m height) was measured with forestry tape to quantify species composition and structural characteristics. In addition, dead trees were identified and their density was obtained to quantify structural feature and disturbance history of the forest. For advance growth, the regeneration layer, i.e. sapling (≤10 cm) and seedling (≤5 cm) of P. roxburghii and associated tree species, were systematically sampled by laying down 5 m × 5 m quadrates in the entire stands. We extracted two cores radii from living trees of P. roxburghii at breast height parallel to the slope contour in the opposite direction using Swedish increment borers (Appendix 1) to document age and radial growth. At least 2 cores from 30 randomly chosen individual trees were obtained with an attempt to achieve the pith of tree. Samples were placed with 2 end sealed plastic straws for safety with relevant information, i.e. DBH, site name and tree number. Six saplings were harvested from three different forests at low, middle, and high altitudinal zones at ground level to help determine the mean age at coring height following the method of Rigg et al. (1998). For each forest stand, soil samples from two pits (1 kg/pit at two different plots) were extracted at a depth of 30 cm using a bucket auger, and a pooled sample of 500 g was analysed in the Swat Agriculture Research Centre (SARC).

Quantitative and laboratory methods
Soil samples were air-dried, sieved through a 2-mm sieve, and analysed for texture (hydrometer method), pH (1:5 mixed soil-water solution using a digital pH meter model AS218), and total organic matter following the Springers-Klee method (Springer and Klee, 1954). The Kjeldahl method was used to determine total nitrogen (Bremner and Mulvaney 1982), available phosphorus was estimated following Olsen (1954), and exchangeable potassium by ammonium ions exchange using a galvanometer. Phytosociological attributes i.e. relative frequency; relative density, relative basal area, and importance value index (IVI) were calculated (Curtis and Mclntosh, 1950) for the overstory and understory species. Absolute values of density/ha and basal area m 2 /ha were obtained for tree and understory strata. The IVI values for tree species and environmental variables were subjected to PC-ORD v. 5.10 for objective classification and ordination of the forests. From various classification methods, hierarchical polythetic agglomerative cluster analysis was used and stands were merged (clustered) into groups with the results being displayed as a dendrogram (McCune and Grace, 2002). The quantitative Sorensen (Bray-Curtis) an effective distance measure for ecological community analysis (McCune and Grace, 2002) and flexible beta as a linkage method (β = -0.25) was applied (Lance and Williams, 1967) which is compatible with Sorensen distance and is space-conserving (Legendre and Legendre, 1998). We used the Kruskal-Wallis test an alternative of one-way ANOVA to compare the environmental variables among different communities.
Detrended correspondence analysis (DCA), an indirect gradient analysis, was applied to identify the major gradient that influence species distribution. Preliminary analyses were made by applying the default option of DCA (Hill and Gauch, 1980) to check the magnitude of change in species composition along the first axis (i.e. gradient length in standard deviation (SD) units). In the present study, DCA estimated the compositional gradient in the vegetation data to be larger than 4.0 SD units for the first axis; thus, canonical correspondence analysis (CCA) was the appropriate ordination method to perform direct gradient analysis (terBraak, 1986). CCA was performed using 12 environmental variables after the exclusion of calcium, magnesium and electric conductivity as these variables were highly correlated and thus showed no significant differences among community types. All the default settings were used for CCA, and a Monte Carlo permutation test (499 permutations) was used to test for significance of the eigen-values of the first conical axis. Intra-set correlations from the CCAs were used to assess the importance of environmental variables.
All the tree species were grouped into different diameter classes of 10-20, and 21-30, and so on in each community type and Weibull function was fitted following Ryniker et al. (2006). Cores obtained were mounted, sanded, and polished with sandpapers of progressively finer grit until a fine surface was obtained as detailed in Stroke and Smiley (1996). These samples were measured to a precision of 0.001 mm under a Khan et 10.6). The age of the trees was determined by counting the number of rings from the outermost ring to the pith. When the pith was not obtained in the core samples, we estimated the age of the ring closest to pith according to its shape of curvature (Xing et al., 2012). It is worth noting that the age of these trees thus obtained is the age of the stem at breast height, because the age does not include the time that the tree grow from the ground to the sampling height. Therefore, following Ogden and Ahmed (1989), rings obtained from sapling were added to the trees age in order to obtain total age of the trees. Linear regression was used to calibrate the relationship between age and size of the trees. For climate change modeling, the metadata obtained were transformed into CSV commadelimited text file format that can be used with the Maximum Entropy (MaxEnt) software . The HADCM3 A2 a climate change scenario (Collins et al., 2001) was used, which represents a grid point model that has a horizontal resolution of 3.75 × 2.5 degrees in longitude × latitude. This corresponds to a spacing between points of approximately 300 km. Bioclimatic layers (see Table 1) in GIS compatible format were downloaded from the Worldclim website and used in the analysis.

Composition and structural patterns
Among the 500 plots measured in 25 forest stands, 9 woody plants belonging to 8 families of 8 genera were identified and classified into 3 communities by cluster analysis (Fig. 2). These communities were clearly isolated in the nonmetric multidimensional scaling (NMS) ordination (results not shown). The summary of associated physiographic, and soil physical and chemical characteristics of these vegetation types were shown in Table 2. Group I was located at high elevations and slopes and low pH, nitrogen (%), K + , enriched with herbaceous species and grasses i.e. Khan    Significant differences between the species composition (A = 0.4037, P < 0.001), structural parameters (i.e. density; A = 0.3520, P < 0.001; basal area; A = 0.3260, P < 0.001) and environmental matrix (A = 0.1865; P < 0.001) were obtained using MRPP. The result of cluster analysis was further clarified by a pair-wise comparison of the communities with high A-values recorded for Group I and II (A = 0.4508, P ≥ 0.001) and Group I and III (A = 0.4106, P ≥ 0.001), whereas substantial similarity was found between Group II and III (A = 0.1727, P < 0.001).
The average Importance values, density and basal areas for the tree species are presented in Tables 3 and 4 according to the Phytosociological groups. The plots from 8 forest sites contained a single species declared as mono-specific community (Group-I) of P. roxburghii with an average density of 634 individuals ha -1 and 47 ± 7.66 basal area m 2 /ha. The proportion of saplings was higher (25%) as compared to seedlings (15%) and dead logs (18%) which significantly contributed to the overall density in this community ( Table 4). Group-II had eleven sites and 6 species dominated by P. roxburghii (IV = 60.6 ± 1.86) and Q. incana (IV = 16.8 ± 3.18), with a total density of 1231 individual's ha -1 and 35.31 basal area m 2 /ha. The density of P. roxburghii and juveniles, i.e. seedlings and saplings was significantly lower as compared to Group I and III, respectively. However, dead trees shared 17% of the total tree density, which is comparatively higher than that of Group III, formed by six sites with seven species led by P. roxburghii (IV = 58.33%) and Acacia modesta (IV = 15%) in the arboreal forest ( Table 2). The overall, density in this group is higher than the prior group (Group II) with an average density of 502 individuals/ha and 25.6 basal area m 2 /ha of the dominant species. The main companions in these communities are Q. baloot with an average IVI that ranged from 5 -9.0% followed by Punica granatum and Monotheca buxifolia. All these species contributed < 2% density and basal area in these communities. However, Olea ferruginea, Ficus palmata, P. dutii and P. granatum were minor associates with < 5% of IVI and 1% density ha -1 and basal area m 2 /ha (Tables 3 and 4).

Vegetation-environment relationship
Twelve environmental factors measured were used in CCA to explore the patterns of species distribution. The correlation among the environmental variables showed Khan et (Fig. 3). This means that the positive axis species are extending their population while the negative axis species showed a similar underlying gradient but in an opposite direction to P. roxburghii along different environmental regimes in the study area.

Diameter patterns
The diameters of P. roxburghii distributed in different stands were pooled based on cluster analysis and presented in the form of histograms that show a unimodal pattern (Fig. 4a-c). The mean diameter of the species was 8.1, 7.8 and 8.7 cm in Group I, II and III, respectively ( Table 5). The standard deviation around the means was approximately Khan et Table 5 indicating that the 2 parameters of Weibull function suggest the lowest mean values for the diameter distributions of P. roxburghii in all groups.

Age and growth rate patterns
Age structure and growth rate pattern of P. roxburghii in different vegetation types (Group I-III) indicated a pattern similar to that of diameter. Bell-shaped age distribution was found in Group I, where the juvenile stages, i.e. 1-10 cm and 10-20 cm classes, accounted for 3 to 4% and 6 to 8% of the total individuals. Substantial numbers of individuals were young (58%) while old trees (>100 years) were rarer (13%) in this group. The mean annual increment was 5.0 (SD = ± 2.92 cm) which is higher than that of the other forest groups located in the area. Comparing Group II with prior (Group I) and proceeding group (Group III), the age structure was generally different. In Group II the youngest individuals were comparable (14%) while young were predominant, becoming 64% of the total. Only 10% of the individuals were between the ages ranged from 100-140 years (Fig. 5), and old trees above this range were entirely absent. Mean annual increment was 3.8 ± 1.92 cm which is higher than the individuals in the prior group. Group III, shared a similar proportion of individuals in the youngest category (9%), while the majority (48%) of them were between 40 and 100 years old. Old trees (>100) were rare and no individuals were found > 180 years old, followed an almost bellshaped pattern. Mean annual increment of P. roxburghii was 2.7 ± 0.98 which is highest among all the groups. Statistically, significant relationships were obtained between age and diameter (y = 0.053x -58.001; R = 0.852 P < 0.001) and age and height (y = 0.0453x -50.002; R = 0.810; P < 0.001) of P. roxburghii in the entire groups using regression equations.

Predictive modeling
MaxEnt predicted the existing distribution of P. roxburghii as restricted to the south of the valley due to its significantly higher population density in these areas (Fig. 6a) whereas, the future distribution model shows even worse situation for the species, i.e., the entire shift of population is to the neighbouring district located in the west. Only a few patches at the western border may still survive because of the favourable conditions available at the end of the century (Fig. 6b). Ground-truth surveying (qualitative) was carried out to check the validity of the present predictive model. The Jackknife analysis (JA) of the present probability distribution of the area under cover (AUC) for the species indicates that all environmental variables are contributing to the AUC (over 0.80) except bio-7 (Temperature annual range), which has a very small share in the gain of AUC. The highest contribution was recorded for climatic variable bio-19 which is the precipitation of the coldest quarter (Fig. 7). The sensitivity and 1-specificity graph shows the best fit of the model with gain of 0.977 and 0.970 gain for AUC of training and test data, respectively (Fig. 8). The figure also clearly indicates the best fit of the model; both training and test omission that are very close with each other and close to the predicted omission (threshold of 0.5). In the future prediction model of P. roxburghii, the trend remains the same; the most important variable was found to be bio-19 and the least contributor was bio-7 (Figs. 9 and 10).

Discussion
Variability of environmental conditions could lead to spatial segregation of the flora and structural alterations in plant communities (Dufour et al., 2006). In the present subtropical pine forest, we have identified three distinct plant communities that vary considerably in species composition and occupied different ecological niches in this region. Contrary to monospecific community of P. roxburghii (Group I), P. roxburghii and Q. incana association (Group II) appeared to prefer medium elevation, degree of slopes, sand (%), pH and potassium, while P. roxburghii and A. modesta preferred the sites with low elevations, pH, lime (%) and high physical and chemical properties ( Table 2). P. roxburghii, typically as a pure community, was mainly restricted to relatively high elevation slope sites, with low amount of clay, silt, nitrogen (%) lime (%) and potassium content. It has been reported by Siddiqui et al. (2009) that P. roxburghii mostly occurs as pure populations in most of the forest patches in northern Pakistan. However, on some drier hill sites it is associated with broadleaved and deciduous tree species like O. ferruginea, M. buxifolia, P. granatum and F. palmata in low and middle elevation zones . Several, other studies of vegetation classification were also conducted in Hindukush and Himalayan mountain ranges in Pakistan, showing that P. roxburghii occurs as a pure community and also a dominant associate with broadleaved and deciduous trees almost in the same elevation ranges Nafeesa et al., 2007). These results correspond to our findings owing to similar ecogeographical regions. However, in neighboring countries like India, Nepal, and Bhutan this species is restricted to the Monsoon belt with summer rain forming associations with P. ponderosa, Q. leucotricophora and Q. semecarpifolia etc., which may be due to the different climate regimes and eco-geographical regions. These differences indicate that vegetation structure, community pattern and differential species distribution are affected by a broad array of biotic, environmental and climatic interactions that overlap and govern community structure in a complex manner (Sarker et al., 2014).  In subtropical forests, natural vegetation often responds to several gradients simultaneously and different combinations of gradients produce divergent responses to the set of gradients (Khan et al., 2013). Our results of ordination support such a continuum and indicate that spatial distribution patterns of P. roxburghii associations do not follow a single environmental gradient, rather, an assortment of gradients account for its compositional variation. This variation could be attributed to physiographic (e.g. elevation, aspect), soil physical (e.g. clay, silt and sand particles) and chemical properties (lime and P) as found previously (Siddiqui et al., 2009;Ahmed et al., 2011;Khan et al., 2014) or climatic variables . These results show that species composition was affected not only by physiographic but also by soil properties, probably due to high variability of the sub-alpine and alpine environments (Champion et al., 1965;Marini et al., 2007). The ordination bi-plot species data shows that species with greater loading on the axes are Q. baloot, O. ferruginea, F. palmata and Q. incana which occupied the negative end of Axis 1, whereas, A. modesta, P. dutii, and M. buxifolia the upper positive end along with P. roxburghii. Such a distribution pattern shows that positive axis species are extending their populations while the negative axis species show a similar underlying gradient but in the opposite direction to P. roxburghii along different environmental regimes. These compositional changes in the current forests seem to be governed by altitudinal gradient acting as an elevation driver among the physiographic factors responsible for compositional variation in the high mountains range (Hong et al., 2015). Several factors like temperature, humidity, snowfall, solar radiation, etc. are associated with elevation, as demonstrated by several authors (e.g., Rana et al., 2011).
Few forest stands with trees in excess of > 110 cm DBH were observed in subtropical P. roxburghii forests and the majority are represented by trees < 50 cm DBH. Ahmed et al. (2006)  The basal area results were confirmed by diameter size analysis, which indicted that all of the communities examined had a unimodal (Bell shaped) distribution, most likely due to high mortality or growth suppression of individuals in the smaller and greater size classes. It has been observed that P. roxburghii grow rapidly in height and diameter at early stage of life, with growth rate decreasing with age (Personal observations). In the current study the mean annual increment in the diameter ≥ 10 and ≤ 110 was 5.0 ± 2.92 in the highlands, 3.8 (±1.92), at middle elevations and 2.7 (±0.98) years/cm in the lowlands. These values indicated that P. roxburghii growth is also sensitive to climate as climatic parameters vary considerably on high altitude (Wang et al., 2004; Huang and Zhang, 2007). Strong relationships were obtained between diameter and age and growth rates, analogous to previous studies (e.g., Khan et al., 2014). Generally, the predication of diameter distribution of stands is of great interest to forest managers for the evaluation of forest resources and predicting future silvicultural treatments (Nano and Montero, 2002). Hence, the use of appropriate statistical models play an important role (Sheykholeslami et al., 2011) by indicating whether the density of smaller trees in a stand is sufficient to replace the current population of larger trees and to help evaluate potential forest sustainability (Rubin and Manion, 2006). We used the three-parameter Weibull function in the present study, which proved effective for fitting the diameter distributions of pine forests. This work will provide baseline information and will substantially increase our knowledge of diameter distributions of pine forests in northern Pakistan. The proposed model diameter distribution will be exceedingly useful for further inventories and the management of these forests as the present study was the first attempt in the area. Other species found in the present study including Q. incana, O. ferruginea, A. modesta, P. granatum and M. buxifolia, were all present as overstory and understory except in Group I, but at low densities. Q. incana and P. granatum have been shown to establish themselves in open canopy on northern aspects. However, their ability to replace P. roxburghii cannot be predicted from the present study, although it seems likely. Due to low sample size Weibull function for these species was not interpreted.
It is now a well-established fact that the climate change is real and the average global temperature in on the rise affecting organisms in one way or another (IPCC, 2013; Root et al., 2003). Generally, organisms show a well pronounced response to the change though the fact that some responses are less pronounced or remain unknown requires dedicated scientific investigations. The current study concludes that most of the tree species in the District will respond to the changing climate in the area and will either shift their habitats, reduce/increase their distribution or in some cases go extinct in the area. In corroboration with the findings of Song et al. (2004), altitude had significant effect on the distribution of species. They have also reported the effect of climate change on the northward movements of the tree species, including Abies spp. and Picea spp. The variable Bio-19 (precipitation of the coldest quarter) was found to be the most important variable in the present predictive distribution model of P. roxburghii. This environmental variable remains equally important for the future predictive model of P. roxburghii distribution. Another personal observation related to P. roxburghii is that it is the most preferred plant by the Forest Department of Pakistan for the purpose of reforestation, but the future distribution model tells a different story about the future of the species i.e., extreme reduction in distribution. It is evident from the results of the study that the associated species with P. roxburghii will have considerable impact on its distribution and density in the future but this warrants further studies to evaluate the impact of reduction and loss of this important plant on the other plant species and communities.

Conclusions
The use of numerical methods and species distribution model (SDM) were exceedingly useful in the exposition of the current classification, compositional variation with relation of environmental factors and in simulating the future distribution trend of P. roxburghii communities in Swat Hindukush range of Pakistan. In this paper, we described the first comprehensive investigation into environmental factors that significantly affect the distribution P. roxburghii communities in a priority conservation area of northern Pakistan renowned for eco-tourism. Despite the low variance in the data explained, the studied variables provide useful insight on plant distribution. Thus, the research approach demonstrated here can help in conserving the remaining natural patches of the forest by providing a basis for vegetation monitoring, mapping and assessing site qualities a priori.