This study analyzed how the environmental conditions constrained the species richness and composition in the four river valleys of Central Nepal i.e. two from Manaslu Conservation Area (MCA) and two from Sagarmatha region. Topographical, bioclimatic and measured variables were used to analyze their effects on the vascular plant diversity along elevation and land use gradients. Altogether, 148 plots were established at five elevation levels between 2,200 m and 3,800 m above the mean sea level. Four land use types namely crop field, meadow, exploited forest and natural forest were sampled at each elevation level. Altogether, 790 species of vascular plants belonging to 114 families were recorded; Asteraceae had the highest number of species (84) followed by Rosaceae (52) and Poaceae (50). Explorative data analysis of species composition by canonical correspondence analysis (CCA) showed that the topographical variables explained the composition better than both the bioclimatic set of variables and the logger data. However, all groups of variables revealed significant effects on species composition. Generalized Linear Model (GLM) also revealed significant effects of elevation, land-use types, slope angle, aspect, temperature and precipitation on species richness. K e y w or d s : Canonical correspondence analysis, elevation, generalized linear model, land use types, multivariate analysis, species richness Effects of the environment on species richness and composition of vascular plants in Manaslu Conservation Area and Sagarmatha region of Nepalese Himalaya S. K. Rai1,2,3, S. Sharma2, K. K. Shrestha2, J. P. Gajurel2,3, S. Devkota3, M. P. Nobis3 and C. Scheidegger3 Species diversity patterns are governed by a varied set of biotic and abiotic factors. Keeping biotic interactions at one end, the abiotic environmental drivers of species distribution has gained much attention in recent studies (Guisan and Zimmermann, 2000). There are several environmental relationships that can be used to describe patterns of species distributions as well as species richness. Changes of species distributions along the latitudinal and elevation gradients are well known since the advent of modern biogeography (Lomolino, 2001; Colwell et al., 2004). The effect of latitude on species richness has been known for a long time (Pianka, 1966; Stevens, 1989). Stevens (1989) has compiled the published literatures showing the effect of latitudinal gradients in the species richness at regional as well as local scales. Species richness and their distribution are also affected by the elevation gradients (Stevens, 1992; McCain and Grytnes, 2010), for example, in mammals (McCain, 2007), birds (Island, 2012) and vascular plants (Trigas et al., 2013). However, both latitude and elevation alone cannot elucidate all the causal biological factors, instead they are proxy for numerous variables such as temperature, moisture energy and so on that change along the elevation (Körner, 2007), topography (Hofer et al., 2008) and latitude (Carpenter, 2005). Land use and geographic factors such as aspect and slope also play important roles in distribution of species in any area (Sanders and Rahbek, 2012). In the Himalaya of Nepal and adjoining countries, the species richness along the elevation gradients have shown the mid-elevation peaks for vascular plant species (Vetaas and Grytnes, 2002; Bhattarai and Vetaas, 2003), ferns (Bhattarai et al., 2004), bryophytes (Grau et al., 2007), lichens (Baniya et al., 2010) and reptiles (Chettri et al., 2010). Those studies have often focused on elevation 1 Department of Plant Resources, Kathmandu, Nepal. E-mail: sanjeevkrai4@gmail.com 2 Central Department of Botany, Tribhuvan University, Kathmandu, Nepal 3 Swiss Federal Institute for Forest, Snow and Landscape Research, Switzerland 3 Banko Janakari, Vol. 26, No. 1 4 pattern in the species richness taken as proxies of changes in temperature, energy and water availability (Bhattarai et al., 2004). In contrast, the topographical variables such as slope angle, aspect or regional differences were rarely analyzed in the Himalayan region (Paudel and Vetaas, 2014).The same hold for microclimates such as point temperature and water availability which might affect upon the species distribution (Geiger et al.,1995). In addition, different land use types also indicate different species communities with varying species richness and pattern. The settlements in the mountains of the Himalaya chiefly rely on agro-pastoral system. The shifting and open grazing system is practiced in the mountain areas. Besides crop farming, the mountain people keep herds of cattle for the supply of food and economic needs. Their energy source is mainly the firewood collected from the nearby forests. All the above activities can lead to the degradation of the natural habitats which affect upon the species diversity in different ways (Cousins, 2009; Honnay et al., 2005). In most of the cases, the species diversity declines in the degraded area due to the fragmentation of the natural forests (Tilman et al., 1997; Maitima et al., 2009). These losses are linked with the disturbances and changes in the nutrient cycling processes such as organic carbon in the soil (Maitima et al., 2009), and available nitrogen (Li et al., 2006). This study aims to find out the effects of most widely used environmental variables such as temperature, precipitation and topography at local as well as regional scales. We have also selected four land use types with an aim to show that species distribution pattern are also the function of land use types. The principle research questions are: (i) How the species richness and composition vary along the altitudinal, precipitational and other topographical indicators?, (ii) How the species are distributed in the different land use types? and (iii) Which types of environmental variables are most suitable to explain the species richness and composition in the Himalaya? Materials and methods Study area The study was conducted in the four river valleys of the two regions of Nepal: Manaslu Conservation Area in Gorkha district and Sagarmatha region in Solukhumbu district (Fig. 1). The study was conducted during 2011 to 2013. Fig. 1: Map of Nepal showing the study districts Gorkha at the center and Solukhumbu in the east In Manaslu Conservation Area (MCA), two river valleys viz. the Nubri and the Tsum (Fig. 2a) were studied. The Nubri valley starts from the confluence of Budhi Gandaki river and Siyar khola (river) near Lokpa. This valley runs along the Budhi Gandaki river upwards in north-west direction. Our study area started from Gap (2,200 m) to Samagaun (3,700 m) located between 28˚31’48.9” N and 28˚35’22.5” N latitude and between 84˚38’29.6” E and 84˚49’51.9” E longitude. The vegetation on the bank of the river near 2,200 m is broad leaved consisting of species such as Benthamidia capitata, Michelia kisopa, Pinus wallichiana and Quercus semicarpifolia. Above 2,500 m altitude, there is a dense forest of Tsuga dumosa, and above 3,000 m altitude, the forest is changed into larch forest (Larix himalaica). At 3,400 m altitude near Shyala, exists a dense forest of Abies spectabilis associated with Hippophae salicifolia and Cotoneaster spp. Similarly, the Tsum valley is oriented towards the north-east along the Siyar khola after the confluence with the Budhi Gandaki river. Our study area is located between 28˚26’19.3” N and 28˚36’56.2” N latitude and between 84˚54’44.3” E and 85˚06’40.4” E longitude. The lower elevation consists of alder (Alnus nepalensis) and pine ( Pinus wallichiana) forests. They are replaced by hemlock (Tsuga dumosa) and Himalayan fir ( Abies spectabilis) at around 3,000 m altitude. The north facing slope of the valley harbors dense vegetation. Larix himalaica forest is dominant at around 3,400 m altitude near Rachen Gumba. The north facing slopes possess more vegetation cover than the south facing slopes. Betula utilis is found upto 3,800 m altitude near Kalung. Most of the south-facing slopes consist of open meadows Rai et al. Banko Janakari, Vol. 26, No. 1 5 intersected by small human settlements such as Chumling, Gho, Chhekampar and Nile. Fig. 2a: Map of Gorkha district with plots overlaid on Nubri river valley on the left and Tsum river valley on the right Fig. 2b: Map of Solukhumbu district with plots overlaid on Dudhkunda river valley on the left and Dudhkoshi river valley on the right In Sagarmatha region, we studied the Dudhkoshi and the Dudhkunda (Fig. 2b) river valleys. The region is famous for the world’s highest mountain, Sagarmatha (Mt. Everest, 8,848 m) and the Sagarmatha National Park. The Dudhkoshi river valley runs northwards along the bank of Dudhkoshi river. The studied plots are located between 27˚40’18.1” N and 27˚49’48.3” N latitude and between 86˚42’3.2” E and 86˚44’25.2” E longitude. The plots located at 2,200 m at Surke and Nakchung and those at Muse and Sengma at 2,600 m elevation are outside the Sagarmatha National Park whereas the rest of the plots are within the boundaries of the Park. The vegetation of the site starts from Schima- Castanopsis and alder (Alnus nepalensis) at 2,200 m and is replaced by Pinus-Rhododendron at mid elevation (3,000 m) and is further replaced by Silver fir-birch-rhododendron at Khumjung (3,800 m). The study area at Dudhkund valley is located between 27˚30’39.9” N and 27˚39’49.1” N latitude and between 86˚34’34.5” and 86˚37’01.6” E longitude, and lies towards the west of Dudhkoshi river valley; the two valleys are separated by a chain of mountains. The Dudhkund valley does not fall inside the Sagarmatha National Park area. The plots, laid at 2,200 m and 2,600 m elevation, are near the settlements and the forests are managed by the local Community Forest User Groups (CFUGs). The forests above 3,000 m elevation are managed by the Government as national forest. The crop fields are not found at and above 3,000 m altitude except one at Taksindu. The study started at Boldok-Kholaghari (2,200 m). Going upwards from Phera (2,600 m), Taksindu (3,000 m) and Sarkaripati (3,400 m), our highest plot was located near Sasarbeni (3,800 m). The vegetation at 2,200 m is Schima-Castanopsis-Alnus, Pinus and then followed by Pinus-Quercus-Rhododendron at mid-elevation. Abies spectabilis forest can be noticed at Sasarbeni (3,800 m). Study design Five elevation levels were investigated with a regular elevation interval of 400 m starting from 2,200 m to 3,800 m. At each elevation level, four land use types were considered viz. (i) natural forest, (ii) exploited forest, (iii) meadow and (iv) crop field (Scheidegger et al., 2010). The category of the land use types were based on the visual observation in accordance with the methods of FAO (Gregorio and Jansen, 2000). The crop fields are cultivated areas where the vegetative cover is created by anthropogenic activities, and so become bare during off-crop season. The meadows are isolated patch or wide area of grazing land where the tree species are less than 20%, and they are also affected by anthropogenic activities such as livestock grazing and grass Rai et al. Banko Janakari, Vol. 26, No. 1 6 collection. The natural forests are far from the human settlements which are rarely intervened by anthropogenic activities. The exploited forests comprise the vegetation not planted by humans but influenced by their actions. This does not require human activities to be maintained in the long-term as compared to the crop fields. All the four land use types were assessed for species records on both sides of the river. Two sample plots (25 m x 2.5 m) were selected randomly per land use type at each elevation level (e.g. 2,200±50 m) on the one side of the river, and the same number were replicated on the another side of the river (Scheidegger et al., 2010). Each plot was divided into 5 m x 2.5 m sub-plots for species record. Thus, each elevation level consisted of eight sample plots (Fig. 3). Crop fields were not found at the elevations of 3,400 m and 3,800 m except a few in some valleys. A total of 148 plots were sampled during the study period of 2011 – 2013. Fig. 3: Schematic diagram of sampling plot design of the study per land use type in all elevation (C= crop field, E = exploited forest, M = meadow, F = forest and the straight line at the center represents the river) All the species within each plot were recorded. If the same species occurred in the next plot, it was recorded as “1” (in the presence of the species). The species recorded in the first plot but not in the second plot were recorded as “0” (in the absence of the species). Two replicate plots of the same land use type were later merged into one. Each plot was visited twice in order to record as many species as possible. To reduce the sampling bias caused by spatial auto-correlation, the replicate of each plot was established at least 50 m away from the first plot (Magurran, 2004). Data source Plant species records as response variable Most of the flowering plant species were identified in the field by using the books written by Polunin and Stainton (1984) and Stainton (1988). The specimens unidentified in the field were identified at the National Herbarium and Plant Laboratories (KATH), Godawari, Lalitpur. The voucher specimens were submitted to the KATH Herbarium. For nomenclature of the species, we followed the Angiospermic Phylogenetic Group (APG III) system (Chase et al., 2009). In the case of the unresolved names (according to APG III), the nomenclature of Press et al. (2000) was adopted. On the other hand, the nomenclature of Iwatsuki (1998) and Fraser-Jenkins (2008, 2011) were used in the case of pteridophytes. The individual species’ presence/absence data in each studied plot were used as the response variable in the current study. Environmental variables as predictor variables The following sets of environmental variables were selected as predictor variables (Table 1). I. The first set of predictor variables included the microclimate (temperature and humidity) data recorded by the logger installed in the field, from 2011 to 2013. The HOBOs (Onset Computer Corporation, Bourne, MA 02532, USA) were used to record air humidity and air temperature 2 m above the ground level in each plot. The HOBOs recorded data in every 30 minutes interval. The soil temperatures were recorded at 10 cm below the ground level using Button (Maxim Integrated, San Jose, CA 95134, USA) in each plot. The soil temperature data were recorded after every 3 hours. The mean, minimum and maximum values of the year-round data were derived using the recorded data afterwards (Table 1). The non-available (NA) values of the data were replaced by the mean of the respective variables so that there would be no loss of data rows in the data frame. II. The second set of predictor variables included the bioclimatic variables extracted from the Worldclim-Global climate data (Hijmans et al., 2005). The data were obtained in 30 arc Rai et al. Banko Janakari, Vol. 26, No. 1 7 seconds (0.93 km x 0.93 km= 0.86 sq. km.) resolution. The latitude and longitude of each plot recorded with the help of Garmin 60S GPS were supplied in the DIVA GIS ver. 7.5.0. The software extracted the interpolated values of the bioclimatic variables from the WORLD CLIM database for each plot. Out of the 19 bioclimatic variables as defined by the USGS Data Series 691 (O’Donnell and Ignizio, 2012), only 10 less correlated variables were chosen for further analysis (Table 1). III. The third set of data contained the information of the topography of the studied area, and were directly recorded in the field. Garmin GPS 60S was used to record the elevation of the plots. Brunton Compass was used to record the aspect while Clinometer was used to record the slope angle of the sample plots. The land use types, the regions and the valleys were considered as the categorical variables and all the others were taken as the ratio variables (Table 1). The above set of variables contained large number of variables. The Hmisc (Harrell et al., 2016) Package was used to check the collinearity among the environmental variables. The Pearson correlation coefficient was used to describe the relationships between the variables. The highly correlated variables (r> 0.7) were taken for analysis (Dormann et al., 2013). Data analysis Initial data recording and management were done using MS Excel and MS ACCESS. The further analyses were performed on the R ver. 3.1.2 (R Core Team, 2015). R-package vegan (Oksanen et al., 2015) was used for the multi-variate ordination analysis. Detrended Correspondence Analysis (DCA) was performed for the species data (Hill and Gauch, Table 1: The list of environmental variables selected from three sets Set Variable acronym Contained information (1) Loggers’ data MaxT.H maximum air temperature recorded by HOBO MaxT.iB maximum soil temperature recorded by iButton MeanT.iB mean soil temperature recorded by iButton MinT.iB minimum soil temperature recorded by iButton MaxH.H maximum air humidity recorded by HOBO meanH.H mean air humidity recorded by HOBO MinH.H minimum air humidity recorded by HOBO (2) Bioclimatic BIO1 annual mean temperature BIO3 Isothermality of temperature BIO5 maximum temperature of warmest month BIO6 minimum temperature of coldest month BIO8 mean temperature of wettest quarter BIO10 mean temperature of warmest quarter BIO14 precipitation of driest month BIO15 precipitation seasonality BIO17 precipitation of the driest quarter BIO19 precipitation of the coldest quarter (3) Spatial REG two regions (Manaslu and Sagarmatha) VAL four valleys HABIE exploited forest HABIF natural forest HABIM meadow ALTG recorded elevation ASP aspect SLOP slope angle Rai et al. Banko Janakari, Vol. 26, No. 1 8 1980) showing the gradient length of the first ordination axis higher than 2.5 standard units. Therefore, we used the Constrained Ordination Method, the unimodal model of the Canonical Correspondence Analysis (CCA) (Ter Braak, 1986). The inertias of all the predictors were compared among each other in order to find out the amount of variances explained by them. The diversity indices like Shannon-Wiener, Simpson and Inverse Simpson indices were calculated using “vegan” R Package (Oksanen et al., 2015). Generalized Linear Model (McCullagh and Nelder, 1989) with quasi-poisson distribution for counts were used to evaluate the relationships between the species richness as response variable and different environmental predictors. The model was fitted against the null model to check for its robustness and performance. The second order polynomial function was also tested, but Fisher’s alpha was not significant. Thus, we proceeded with the first order linear model. Results and discussion The study revealed 790 vascular plant species of 337 genera within 114 plant families. The highest number of species were recorded for Asteraceae (84 spp.) followed by Rosaceae (52 spp.), Poaceae (50 spp.) and Fabaceae (38 spp., Fig. 4). Fig. 4: Bar diagram showing the representative families, number of species on the Y-axis and families on the X-axis (families representing more than 10 spp. are included) Species composition The Detrended Correspondence Analysis (DCA) of the species values against the plots studied was performed. All of the DCA axes were more than 2.5 standard units; therefore, the data were further analyzed using CCA. The species data were further constrained separately by the logger data, bioclimatic data and topographical variables for CCA analysis. The performances of the variables are presented in Table 2. The CCA plots show the effect of the environmental variables on the species composition (Fig. 5, 6 and 7). The distribution of the species were found to be affected by the temperature along the CCA axis-1 and the humidity along the CCA axis-2 (Fig. 5) This clearly showed that the temperature and humidity were controlling environmental factors for the distribution of the species (Table 2). In terms of percentage, the variation explained by the CCA axis-1 and the CCA axis-2 were ~41.8% and ~27.7%, respectively; thus, 69% of the variation were explained by the two CCA axes (Table 3). Fig. 5: CCA plot showing the species composition constrained by humidity and soil temperature; the crosses indicating the species, the circles indicating the plots and the arrows showing the predictors The precipitation seasonality (BIO15) possesses the longest gradient length to shape the species distribution. Isothermality (BIO3) refers to the percentage of the mean diurnal range divided by the annual temperature. Thus, the growing days and length of the days which shape the temperature pattern has also significant contribution for species distribution. Precipitation of the driest month (BIO14) is another contributor for species distribution. Annual mean temperature (BIO1), mean temperature of coldest month (BIO6), mean temperature of the warmest quarter (BIO10) and mean temperature of the wettest quarter (temperature combined with the precipitation, BIO8) were found to have the significant effect on the species composition in the study areas (Table 2 and Fig. 6). Rai et al. Banko Janakari, Vol. 26, No. 1 9 The CCA axes of the bioclimatic variable were found to have performed less than the CCA axes obtained from the logger data. The CCA axis- 1 was found to have explained 24.39% of the variation followed by the CCA axis-2 (20.19%), the CCA axis-3 (14.82%) and the CCA axis-4 (12.24%). Thus, a total of 72% of the variation was found to be explained by these four axes (Table 3). Fig. 6: CCA plot showing the species composition constrained by the bioclimatic variables; the crosses indicating the species, the circles indicating the plots and the arrows showing the predictors The results obtained by constraining species with the annual temperature and precipitation mean and their derivatives show that not only the mean, minima and maxima of the temperature and precipitation are important but also their combined effect are equally important to shape the distribution of the species in the given environmental hyper-volume (Hutchinson, 1957). The predictor variables constructed with the derivatives of temperature and precipitation alone and combined have the physiological role in the germination, growth and proliferation (Wright et al., 2006). Soil temperatures are important for the physiology of the cell, water availability and nutrient uptake from the soil (Korner, 2003). Temperature is related with the energy balance as well (Scherrer et al., 2011). Topographical variables also show significant effect upon the species composition (Fig. 7, Table 2). Table 2: The test statistics expressed by the environmental variables while constraining the species composition (by “margin” i.e. each marginal term analyzed in a model with all other variables) Variable Set Code Df Chi Square F Pr(>F) Significance codes Loggers MeanT.iB 1 0.1385 1.7665 0.001 * * * MinT.iB 1 0.1450 1.8498 0.001 * * * MaxH.H 1 0.1187 1.5135 0.002 * * MinH.H 1 0.1303 1.6617 0.001 * * * Residual 143 11.2110 Bioclimatic BIO1 1 0.1259 1.6474 0.001 * * * BIO3 1 0.1664 2.1761 0.001 * * * BIO6 1 0.1078 1.4100 0.001 * * * BIO8 1 0.1132 1.4808 0.001 * * * BIO10 1 0.1185 1.5506 0.001 * * * BIO14 1 0.1617 2.1150 0.001 * * * BIO15 1 0.1458 1.9071 0.001 * * * Residual 140 10.7033 Spatial REG 1 0.2408 3.2757 0.001 * * * VAL 1 0.2368 3.2217 0.001 * * * HABI 3 0.5225 2.3695 0.001 * * * ALTG 1 0.3770 5.1288 0.001 * * * ASP 1 0.1242 1.6903 0.001 * * * SLOP 1 0.1109 1.5089 0.001 * * * Residual 139 10.2174 Significance codes: ‘***’ for P=0.001, ‘**’ for P=0.002 Rai et al. Banko Janakari, Vol. 26, No. 1 10 Fig. 7: CCA plot showing the species composition as shaped by the topographical variables; the crosses indicating the species, the circles indicating the plots and the arrows showing the predictors Elevation (ALTG) was found to be one of the significant variables for the species composition in our study (Table 2). It is a surrogate of a number of environmental factors,e.g. temperature, which in turn stands for energy, water etc. Slope angle (SLOP) and aspect of the plots were also found to be significant contributors for the species composition (Table 2). More than 60% of the variation was found to be explained by the three CCA axes produced by constraining species composition with the topographical variables. The CCA axis-1 explains ~25% followed by the CCA axis-2 (~21.8%) and the CCA axis-3 (~15.2%) (Table 3). The valleys (VAL) were also found to be significant for species composition (Table 2). The Sagarmatha region receives more annual precipitation (average 1,640.95 mm) as compared to the MCA region (average 545.36 mm, Hijmans et al., 2005). The valleys in the Sagarmatha region are geographically nearer to the Bay of Bengal, the origin of the monsoon rain system, and are less rain-shadowed by the high mountains. In contrast, the MCA valleys are geographically farther from the Bay of Bengal and rain-shadowed by Mt. Ganesh (7,422 m). Species richness For each environmental variable, annual model was first created and was tested with the first order Generalized Linear Model (GLM). Transect-wise species richness was taken as response variable which regressed against different environmental variables as predictor. These included land use types (LUT), elevation (ALTG), precipitation seasonality (BIO15), annual precipitation (BIO12), slope angle (SLOP) and aspect (ASP) of the plots. These developed models were tested Table 3: Percentage of variation explained by the CCA axes when species richness were constrained with the predictor variables Data set Constrained inertia CCA axes Eigenvalues Percentage variation explained Cumulative variation % Loggers' Set 0.750 CCA1 0.3132 41.77 CCA2 0.2076 27.68 69 CCA3 0.1223 16.31 86 CCA4 0.1072 14.30 100 Bioclimatic Set 1.258 CCA1 0.3068 24.39 CCA2 0.2539 20.19 45 CCA3 0.1864 14.82 59 CCA4 0.1539 12.24 72 CCA5 0.1512 12.02 84 CCA6 0.1081 8.60 92 CCA7 0.0975 7.75 100 Spatial Set 1.744 CCA1 0.4344 24.97 CCA2 0.3780 21.72 47 CCA3 0.2643 15.19 62 CCA4 0.2277 13.09 75 CCA5 0.1937 11.13 86 CCA6 0.1073 6.17 92 CCA7 0.0846 4.86 97 CCA8 0.0542 3.11 100 Rai et al. Banko Janakari, Vol. 26, No. 1 11 among each other by using “F” statistics. Over- dispersed residual of errors were standardized after application of “quasipoisson” family of distribution of error. The significant environmental variables with deviance and “F” values are indicated in Table 4. The graphics of some more interpretable and statistically significant variables are shown in Fig. 8a–8d. Fig. 8a: Species richness versus elevation of the plots Fig. 8b: Species richness versus land use types Fig. 8c: Species richness versus precipitation seasonality Fig. 8d: Species richness versus annual precipitation Note: In Fig. 8b, C = crop field, E = exploited forest, F = natural forest and M = meadow; in Fig. 8c, units are precipitation coefficients and in Fig. 8d, precipitation is in mm. The species richness increased with the increase in the elevation of the plots studied. The previous studies in Nepal showed the unimodal richness pattern with elevation (Baniya et al., 2010; Grau et al., 2007; Vetaas and Grytnes, 2002). Those studies analyzed long elevation gradients whereas this study considered relatively short elevation gradient between 2,200 m and 3,800 m above the Table 4: Test statistics of the generalized linear model (GLM) of species richness against the individual environmental variables Code Predictors Resid. df Resid. dev. Deviance F Pr(>F) Significance codes LUT Land Use Types 144 1460 75445 2608.9 < 2.20E-16 * * * ALTG Elevation 146 1357 75547 8574.1 < 2.20E-16 * * * BIO15 Precipitation Seasonality 146 1373 75532 8460.7 < 2.20E-16 * * * BIO12 Annual Precipitation 146 1462 75443 7934.5 < 2.20E-16 * * * SLOP Slope Angle 146 1464 75441 7928.0 < 2.20E-16 * * * ASP Aspect 146 1467 75437 7895.5 < 2.20E-16 * * * Significance codes: ‘***’ for P< 0.001 Rai et al. Banko Janakari, Vol. 26, No. 1 12 mean sea level. The short gradient in our study was not sufficient to test the species richness humps. However, there are studies which show the plateau of species richness of birds at high elevation (Patterson et al., 1998). An elevation limit of species occurrence is expected for high mountains e.g., the Himalaya, always covered with snow and the permafrost. The hump shaped unimodal distribution of species richness are expected for such restriction in the absence of any environmental gradients (Colwell and Lees, 2000; Colwell et al., 2004) or isolation from other zonal communities (Lomolino, 2001). However, hump is a union of linear segments at local scale. The result was obtained from only 1,600 m elevation. Thus, the result from this study could be a local phenomenon rather than the large-scaled unimodal pattern found by the earlier researchers. This interpretation resembles quite similar to that of Baniya et al. (2012). Four land use types namely (i) crop field, (ii) meadow, (iii) exploited forest, and (iv) natural forest were studied. The exploited forests were more species-rich, followed by the meadow, the natural forest and the crop field. The soil use intensity and fragmentation are thought to be loss of biodiversity (Cousins, 2009; Honnay et al., 2005; Maitima et al., 2005). This explains the less richness in the crop field. The species richness in the exploited forest is described by the intermediate disturbance hypothesis (Connell, 1978) and some empirical studies (Townsend and Scarsbrook, 1997). In our study, the species richness was found to have increased significantly with the increase in the annual precipitation and seasonality (Fig. 8c and 8d). Precipitation seasonality is the coefficient of variation of the monthly precipitation. The four valleys studied have different precipitation seasonality, which is explained by this study. The different valleys receive varying degree of precipitation shaping different scale of species richness and their pattern (O’Brien, 1993; Pauses and Austin, 2001). The species richness and composition pattern are also affected by the slope and aspect of the sampling plots (Nuzzo, 1996). The south-facing and steeper slopes are drier than the north-facing slopes, and more number of species is expected towards the wet areas (Kassas and Zahran, 1971; Pook and Moore, 1966). The temperature is also significantly affected by the aspects in the mountain environments at point-scale (Kroner, 2003; Parker, 1991). The variation in the slope and aspect, thus, result in the variation of the soil moisture, nutrient cycling and availability of energy dissipation (Mohammad, 2008) resulting in different composition and richness (Carmel and Kadmon, 1999). Conclusion Altogether 790 vascular plant species belonging to 114 families were recorded from six river valleys studied. Asteraceae (84 spp.) was the most dominant family among them. The three sets of environmental variables were used to study their effect on the species composition and species richness of vascular plants. The loggers recorded the microclimate data of each plot. Soil temperature and humidity of the plots affected the plants composition significantly. Out of 19 bioclimatic variables only seven showed significant effect on the plant composition. Annual mean temperature (BIO1), isothermality of the temperature (BIO3), minimum temperature of the coldest month (BIO6), mean temperature of the wettest quarter (BIO8) were the temperature related variables. Precipitation of the driest month (BIO14) and precipitation seasonality (BIO15) also were significant variables. The topography of the plots (elevation, aspect and slope) affected the vascular plant composition significantly. Nearly 50 percents of the variations were explained by two axes of the CCA in all three sets of environmental variables. Four land use types were considered during the study. These land use types also affected the species richness and composition significantly. The results of the study are in accordance with the previous studies. However, the unimodal hump of the species richness distribution was not revealed due to shorter elevation gradient in this study. Acknowledgements We are grateful to the CDB-WSL Project run by the Central Department of Botany, Tribhuvan University, Nepal for providing the fund and logistics for the field trip. We are also thankful to the Swiss Federal Institute for Forest, Snow and Landscape Research, WSL, Switzerland for providing us scholarships for data management and processing. This study was funded by the Swiss National Science Foundation (JRP Rai et al. Banko Janakari, Vol. 26, No. 1 13 IZ70Z0_131338/1 to CS). The first author is indebted to the Ministry of Forests and Soil Conservation, the Government of Nepal for granting him a three-year study leave to accomplish the study. Our sincere thanks go to Dr. Keshab Raj Rajbhandary, Dr. Khem Raj Bhattarai and Rita Chhetri for identification of plants at the National Herbarium and Plant Laboratories (KATH), Godawari. The employees of KATH and Wilder Places Treks, Kathmandu are duly acknowledged. We appreciate the contribution of Ms. Laxmi Sankhi, Ms. Srijana Shah and Mr. Hem B. Katuwal in data collection during the field-work. We acknowledge Dr. Chitra Baniya for sharing his ideas and analysis methods while developing this paper. References Baniya, C. B., Solhøy, T., Gauslaa, Y. and Palmer, M. W. 2010. The elevation gradient of lichen species richness in Nepal. The Lichenologist 42 (1): 83–96. Baniya, C. B., Solhøy, T., Gauslaa, Y. and Palmer, M. W. 2012. Richness and composition of vascular plants and cryptogams along a high elevation gradient on Buddha Mountain, Central Tibet. Folia Geobot 47: 135–151. Bhattarai, K. R. K. and Vetaas, O. R. 2003. Variation in plant species richness of different life forms along a subtropical elevation gradient in the Himalayas, east Nepal. Global Ecology and Biogeography 12 (4): 327–340. Bhattarai, K. R., Vetaas, O. R. and Grytnes, J. A. 2004. Fern species richness along a central Himalayan elevation gradient, Nepal. Journal of Biogeography 31 (3): 389–400. Carmel, Y. and Kadmon, R. 1999. Effects of grazing and topography on long term vegetation changes in a Mediterranean Ecosystem in Israel. Plant Ecology 145: 239–250. Carpenter, C. 2005. The environmental control of plant species density on a Himalayan elevation gradient. Journal of Biogeography 32 (6): 999–1018. Chase, M. W., Fay, M. F., Reveal, J. L., Soltis, D. E., Soltis, P. S., Peter, F. and Kenneth, J. 2009. An update of the Angiosperm Phylogeny Group classification for the orders and families of flowering plants : APG III, 105–121. Chettri, B., Bhupathy, S. and Acharya, B. K. 2010. Distribution pattern of reptiles along an eastern Himalayan elevation gradient, India. Acta Oecologica 36 (1): 16–22. Colwell, R. K. and Lees, D. C. 2000. The mid- domain effect: geometric species richness, 15 (2): 70–76. Colwell, R. K., Rahbek, C. and Gotelli, N. J. 2004. The mid-domain effect and species richness patterns: what have we learned so far? The American Naturalist 163 (3): E1–23. Connell, J. H. 1978. Diversity in tropical rain forests and coral reefs. Science 199: 1302– 1310. Cousins, S. A. O. 2009. Extinction debt in fragmented grasslands: paid or not? Journal of Vegetation Science 20: 3–7. Dormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré, G. and Lautenbach, S. 2013. Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography 36 (1): 27–46. Fraser-Jenkins, C. 2011. Nepal’s little known pteridophytes, the hidden work of David Don, and the geography and distribution of Indo-Himalayan ferns. www.groups.yahoo. com/group/Indian-Ferns accessed on 13 June, 2012. Fraser-Jenkins, C. 2008. Taxonomic Revision of Three Hundred Indian Subcontinental Pteridophytes with a Revised Census-list. Bhisen Singh Mahendra Pal Singh, Dehra Dun, India. Geiger R., Aron, R. H. and Todhunter, P. 1995. The Climate near the Ground. Friedr. Viewegund. Sohn Verlagsges. mbH, Braunschweig, Wiesbaden, 327–406. Grau, O., Grytnes, J. and Birks, H. J. B. 2007. A comparison of altitudinal species richness patterns of bryophytes with other plant groups in Nepal, Central Himalaya. Journal of Biogeography 34 (11): 1907–1915. Rai et al. Banko Janakari, Vol. 26, No. 1 14 Gregorio, A. D. and Jansen, L. J. M. 2000. Land Cover Classification System (LCCS): Version 1.0. FAO. www.fao.org. accessed on 01 March, 2011. Guisan, A. and Zimmermann, N. E. 2000. Predictive habitat distribution models in ecology. Ecological Modelling 135 (2–3): 147–186. Harrell, F. E. Jr, with contributions from Charles Dupont and many others. 2014. Hmisc: Harrell Miscellaneous. R package version 3.17–2. https://CRAN.R-roject.org/package=Hmisc accessed on 10 July, 2014. Herzog, S. K., Kessler, M. and Bach, K. 2005. The elevation gradient in Andean bird species richness at the local scale: a foot hill peak and a high elevation plateau. Ecography 28: 209–222. Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. and Jarvis, A. 2005. Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology 25: 1965–1978. Hill, M. O. and Gauch, H. G. 1980. Detrended correspondance analysis: an improved ordination technique. Vegetation 42: 47–58. Hofer, G., Wagner, H. H., Herzog, F. and Edwards, P. J. 2008. Effects of topographic variability on the scaling of plant species richness in gradient dominated landscapes. Ecography 31: 131–139. Honnay, O., Jacquemyn, H., Bossuyt, B. and Hermy, M. 2005. Forest fragmentation effects on patch occupancy and population viability of herbaceous plant species. New Phytologist 166 (3): 723–736. Hutchinson, G. E. 1957. Concluding remarks. Cold Spring Harbor Symposium on Quantitative Biology 22: 415–427. Island, H. 2012. Bird species richness along an elevation gradient in a forest at Jianfengling, Hainan Island, China. Zoological Studies 51 (3): 362–371. Iwatsuki, K. 1998. An enumeration of the Pteridophytes of Nepal. http://www.um.u- tokyo.ac.jp/publish_db/Bulletin/no31/ no31018.html accessed on 14 June, 2013. Kassas, M. and Zahran, M. A. 1971. Plant life on the costal mountains of the Red Sea, Egypt. Journal of Indian Botanical Society 50a: 571–586. Körner, C. 2003. Alpine Plant Life. Springer- Verlag, Berlin Heidelberg GmbH, Germany. Körner, C. 2007. The use of “altitude” in ecological research. Trends in Ecology and Evolution 22: 569–574. Li, Y., Xie, D. and Wang, S. 2006. Impact of land cover types on the soil characteristics in karst area of Chongqing. Journal of Geographical Sciences 16 (2): 143–154. Lomolino, M. V. 2001. Elevation gradients of species-density: historical and prospective views. Global Ecology and Biogeography 10: 3–13. Magurran, A. E. 2004. Measuring Biological Diversity. Blackwell Science Ltd, USA. Maitima, J. M., Mugatha, S. M., Redi, S. R., Gachimbi, L. N., Majule, A., Lyaruu, H., Pomery, D., Mathai, S. and Mugisha, S. 2009. The linkages between land use changes, land degradation and biodiversity across East Africa. African Journal of Environmental Science and Technology 3 (10): 310–325. McCain, C. M. and Grytnes, J. A. 2010. Elevation Gradients in Species Richness. John Wiley and Sons Ltd, Chichester, UK. McCain, C. M. 2007. Area and mammalian elevation diversity. Ecology 88 (1): 76–86. McCullagh, P. and Nelder, J. 1989. Generalised Linear Models. 2nd edition. Chapman and Hall, London, UK. Mohammad, A. 2008. The effect of slope aspect on soil and vegetation characteristics in Southern west Bank. Bethlehem University Journal 27: 11–27. Nuzzo, V. A. 1996. Structure of cliff vegetation on exposed cliffs and the effect of rock climbing. Canadian Journal of Botany 74: 607–617. Rai et al. Banko Janakari, Vol. 26, No. 1 15 O’Brien E.M. 1993. Climatic gradients in woody plant species richness: towards and explanation based on an analysis of southern Africa’s woody flora. Journal of Biogeography 20: 181–198. O’Donnell, M. S. and Ignizio, D. A. 2005. Bioclimatic predictors for supporting ecological applications in the Conterminous United States. US Geological Survey Data Series 691, 10. Oksanen, J., Blanchet, F. G., Kindt, R., Legendre, P., Minchin, P. R., O’Hara, R. B., Simpson, G. L., Solymos, P., Stevens, M. H. H. and Wagner, H. 2015. vegan: Community Ecology package. R package version 2.3–0. https://CRAN.R-project.org/package=vegan. accessed on 13 March, 2016. Parker K. C. 1991. Topography, substrate and vegetation patterns in the northern Sonoran desert. Journal of Biogeography 18: 151–163 Patterson, B. D., Stotz, D. F., Solari, S., Fitzpatrick, J. W. and Pacheco, V. 1998. Contrasting patterns of elevation zonation for birds and mammals in the Andes of south-eastern Peru. Journal of Biogeography 25: 593–607. Paudel, S. and Vetaas, O. R. 2014. Effects of topography and land use on woody plant species composition and beta diversity in an arid Trans-Himalayan landscape, Nepal. Journal of Mountain Science 11 (5): 1112– 1122. Pauses, J. D. and Austin, M. P. 2001. Patterns of plant species richness in relation to different environments: an appraisal. Journal of Vegetation Science 12: 153–166. Pianka, E. R. 1966. Latitudinal gradients in species diversity: a review of concepts. American Naturalist 100: 33–34. Polunin O. and Stainton A. 1984. Flowers of the Himalaya. Oxford University Press. New Delhi, India. Pook, E. W. and Moore, W. E. 1966. The influence of aspect on the composition and structure of dry sclerophyll forest on Black Mountain, Canberra. Australian Journal of Botany 14 (2): 223–242. Press, J. R., Shrestha, K. K. and Sutton, D. A. 2000. Annotated Checklist of the Flowering Plants of Nepal. The Natural History Museum, London, UK. R Core Team. 2015. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project. org. Sanders, N. J. and Rahbek, C. 2012. The patterns and causes of elevation diversity gradients. Ecography 35 (1): 1–3. Scheidegger, C., Nobis, M. P. and Shrestha, K. K. 2010. Biodiversity and livelihood in land-use gradients in an era of climate change – outline of a Nepal-Swiss research project, 7–17. Scherrer, D., Körner, C., Korner, C., Körner, C., Ko, C. and Korner, C. 2011. Topographically controlled thermal-habitat differentiation buffers alpine plant diversity against climate warming. Journal of Biogeography 38 (2): 406–416. Stainton, A. 1988. Flowers of the Himalaya: A Supplement. Oxford University Press, New Delhi, India. Stevens, G. C. 1989. The latitudinal gradient in geographical range: how so many species coexist in the tropics. The American Naturalist 133 (2): 240–256. Stevens, G. C. 1992. The elevation gradient in altitudinal range: an extension of rapoport’s latitudinal rule to altitude. The American Naturalist 140 (6): 893–911. Ter Braak, C. J. F. 1986. Canonical correspondance analysis: a new eigenvector t technique for multivariate direct gradient analysis. Ecology 67 (5): 1167–1179. Tilman, D., Knops, J., Wedin, D., Reich, P., Ritchie, M. and Sieman E. 1997. The influence of functional diversity and composition on ecosystem processes. Science 277: 1300- 1302. Townsend, C. R., and Scarsbrook, M. R. 1997. The intermediate disturbance hypothesis, refugia and biodiversity in streams. Limnology and Oceanography 42 (5): 938–949. Rai et al. Banko Janakari, Vol. 26, No. 1 16 Trigas, P., Panitsa, M. and Tsiftsis, S. 2013. Elevation gradient of vascular plant species richness and endemism in Crete-the effect of post-isolation mountain uplift on a continental island system. PloS One 8 (3): e59425. accessed on 10 June, 2013. Vetaas, O. R., and Grytnes, J. A. 2002. Distribution of vascular plant species richness and endemic richness along the Himalayan elevation gradient in Nepal. Global Ecology and Biogeography 11 (4): 291–301. Wright, I. J., Reich, P. B., Atkin, O. K., Lusk, C. H., Tjoelker, M. G. and Westoby, M. 2006. Irradiance, temperature and rainfall influence leaf dark respiration in woody plants: evidence from comparisons across 20 sites. The New Phytologist 169 (2): 309–19. Rai et al.