Phylogenetic diversity patterns in Himalayan forests reveal evidence for environmental filtering of distinct lineages

Large-scale environmental gradients have been invaluable for unraveling the processes shaping the evolution and maintenance of biodiversity. Environmental gradients provide a natural setting to test theories about species diversity and distributions within a landscape with changing biotic and abiotic interactions. Elevational gradients are particularly useful because they often encompass a large climatic range within a small geographical extent. Here, we analyzed tree communities in plots located throughout Arunachal Pradesh, a province in northeast India located on the southern face of the Eastern Himalayas, representing one of the largest elevational gradients in the world. Using indices of species and phylogenetic diversity, we described shifts in community structure across the landscape and explored the putative biotic and abiotic forces influencing species assembly. As expected, species richness and phylogenetic diversity decreased with increasing elevation; however, contrary to predictions of environmental filtering, species relatedness did not show any clear trend. Nonetheless, patterns of beta diversity (both taxonomic and phylogenetic) strongly suggest lineage filtering along the elevational gradient. Our results may be explained if filtering is driving the assembly of species from distinct evolutionary lineages. New metrics exploring community contributions to regional taxonomic and phylogenetic beta diversity provided additional evidence for the persistence of unique communities at high elevations. We suggest that these patterns may be consistent with filtering on glacial relicts, part of once more diverse clades with convergent traits suited to climates at the last glacial maximum, resulting in random or over-dispersed community assemblages at high elevations. We propose that these high-elevation sites with evolutionarily distinct species represent possible regions for conservation priority that may provide refugia for species threatened by current warming trends.


INTRODUCTION
Species diversity patterns have been extensively studied along a number of large-scale environmental gradients and have advanced our understanding of the processes shaping species assemblages. It is well established that species richness generally decreases with latitudinal distance from the equator and with increasing elevation (Stevens 1989, Rahbek 1995, Lomolino 2001, likely driven by factors including climate, energy, and potential evapotranspiration (Currie 1991, Givnish 1999, Hawkins et al. 2003, although opposing patterns have been observed (Rahbek 1995). However, a shortcoming of these observations is that analyses of species richness patterns assume that species are equivalent and independent of one another when evolutionary history might be an important additional factor shaping diversity gradients (Davies et al. 2004, Mittelbach et al. 2007, Davies and Buckley 2012, Kerkhoff et al. 2014. The potential importance of evolutionary process on diversity patterns has long been recognized. For example, hypotheses explaining the unusually high diversity in the tropics have described equatorial regions as either cradles or museums of diversity, based on concepts of speciation and long-term extinction survival, respectively (Stebbins 1974). More recently, there has been growing appreciation that evolutionary history might structure not only richness, but also the composition of species assemblages (Webb 2000, Cavender-Bares et al. 2004, Vamosi et al. 2009). For example, closely related species might share similar ecological preferences and tolerances and thus tend to be found in similar environments; however, at local scales, it is possible that competitive displacement might occur among species if resource requirements are too similar (Chesson 2000, Webb et al. 2002.
By combining information on evolutionary history with a large-scale environmental gradient, we explored diversity patterns in India's northeastern-most province, Arunachal Pradesh, home to the southern face of the Himalayas-the largest elevational gradient in the world. We used phylogenetically explicit metrics to unravel the evolutionary processes shaping the local flora, allowing us to explore shifts in community diversity and evolutionary structure with elevation. It has been suggested that communities in abiotically harsher environments (such as those found at high elevations) will tend to be composed of more closely related species because phylogenetic niche conservatism and strong environmental filtering would select for a subset of lineages adapted to these extreme environments (Webb et al. 2002, Bryant et al. 2008. However, empirical studies have sometimes shown opposite trends with phylogenetic clustering at low elevations, or in warmer climates, and over-dispersion at higher elevations, or in colder climates (Bryant et al. 2008, Gonzalez-Caro et al. 2014). In addition, simplistic interpretations inferring process from pattern have recently attracted criticism (Mayfield and Levine 2010, Gerhold et al. 2015, Kraft et al. 2015, Cadotte and Tucker 2017. Experiments provide a more powerful approach for testing mechanistic explanations, but they are not practical for exploring large-scale biogeographic patterns. A robust alternative is to first generate hypotheses based on a priori expectations and then evaluate these hypotheses with empirical data. We pursue this approach here. The apparent conflict between theory and empirical studies on the relationship between community structure and elevation highlights the need for a better understanding of the multiple processes determining community assembly. For example, it is well recognized that the convergent evolution of relevant traits in stressful environments could lead to over-dispersed or random patterns of community structure as opposed to clustering (Webb et al. 2002, Bryant et al. 2008, Read et al. 2014, Kraft et al. 2015. Similarly, the convergent evolution of traits conferring a competitive advantage could also lead to over-dispersion but might be misinterpreted as evidence for niche overlap driving competitive exclusion, particularly if the pattern is observed in a unique or filtering environment (Ghosh-Harihar 2014). Relict taxa-taxa that have shifted their ranges to refugia in montane regions during periods of warming, for example-could also have a large influence on phylogenetic community structure, but have been less well studied (Birks andWillis 2008, Stewart et al. 2010). Previous work has explored the genetic imprint of refugia on intraspecific genetic variation and local population dynamics, often through phylogeographic approaches (Hooghiemstra and van der Hammen 1998, Tribsch and Sch€ onswetter 2003, Mayle 2004, Vargas 2007, Provan and Bennett 2008; but to our knowledge, refugia have not been examined within the context of community phylogenetic structure. The relict taxa that have survived successive climate-driven extinction cycles might often be range restricted (Habel and Assmann 2009) and phylogenetically distinct from the regional community (Fryxell 1962, Provan andBennett 2008). The presence of relict species might thus tend to increase local phylogenetic diversity and contribute to patterns of community phylogenetic dispersion.
The Eastern Himalayan region offers a heterogeneous landscape along one of the largest elevational and climatic gradients on Earth; its highest peak reaches 7060 m above sea level (Mount Everest reaches 8848 m). The vegetation in the region varies from tropical forests to subtropical, temperate, and gymnosperm-dominant alpine forests (Roy and Behera 2005). The rapid transition between forest and climatic zones makes the Himalayan region especially interesting for phylogenetic diversity studies. However, it remains largely underexplored. Previous work in the Eastern Himalayas has suggested that the highest species richness occurs in forest transition zones, between tropical semi-evergreen to subtropical evergreen, and subtropical evergreen to broadleaf forests (Behera and Kushwaha 2007), and it is estimated that 30-40% of the~6000 plant species in Arunachal Pradesh are endemic (Myers 1988, Baishya 1999, Roy and Behera 2005. While there are conflicting reports on the relative frequency and distribution of endemics in the region, with evidence for high endemism in both the low-elevation tropics and the highelevation alpine regions (Behera et al. 2002, Roy andBehera 2005), a study of endemism and species diversity in nearby Nepal reports that the highest proportion of endemic vascular plant species is found between 3800 and 4200 m (Vetaas and Grytnes 2002).
Here, we analyzed data on forest plots established between 2007 and 2010, which are distributed throughout Arunachal Pradesh. The plots were initially designed to investigate the species richness in the region, but because they encompass a vast, 4000-m elevational span, these data provide a unique opportunity to explore shifts in community structure and richness. We used a combination of phylogenetic and taxonomic measures of diversity, indices of phylogenetic dispersion, and a regional phylogeny to evaluate elevational trends in richness and phylogenetic diversity and to identify phylogenetically distinct sites that might represent climatic refugia.

Study site
The study sites are located in the northeastern state of Arunachal Pradesh, India (27.06°N, 93.37°E; Fig. 1). Arunachal Pradesh borders the states of Assam and Nagaland and shares international borders with Bhutan, Myanmar, and China. Geographically, it is the largest among the northeastern Indian states (83,743 km 2 ) and has a population density of 17 individuals per square kilometer (Census of India 2011). The climate of Arunachal Pradesh ranges from tropical to alpine and varies with the elevation.
The site data include species-level identifications of the trees and shrubs found within 352 belt transects (referred to as plots herein). Prior to establishing sampling sites, major landscape elements were identified on 166 topographic sheets encompassing the region, and plots were selected based on accessibility and to maximize environmental diversity. Plots ranged in elevation from 87 to 4090 m above sea level, representing four distinct forest ecosystems: tropical evergreen/semi-evergreen, subtropical broadleaf/ pine, temperate broadleaf/coniferous, and alpine.
Following preliminary examination of the data, several plots were excluded from the study, including those located in plantations or in fields without any tree or shrub species present, as noted by the field researchers. We derived a disturbance index for the remaining 291 plots based on qualitative observations made in the field. The types of disturbance noted were as follows: humus presence, lopping, cut stumps, litter, soil removal, and grazing. We assigned a value of 1 for each type of disturbance if the value was noted as either recurrent or occasional and a value of 0 if the disturbance was noted as absent. We summed these values to generate an index value for each plot, which we then regressed against elevation. We found no relationship between disturbance and elevation (R 2 adj = À0.003, P = 0.9715). Because of sampling practicalities, there was some variation in plot size (which ranged from 500 to 5000 m 2 ); we therefore divided the number of individuals of each species by total plot area, yielding the number of individuals per m 2 . Not all individuals could be identified to species level and the fraction of identified individuals varied among plots. However, we chose to analyze all 291 plots after determining that elevational trends remained similar among groups of plots with differing percent identification (Appendix S1: Table S1).

Phylogeny reconstruction
A molecular phylogeny was reconstructed for the tree species in the study plots using nucleotide sequence information from GenBank. We used three plant DNA barcodes: rbcL, matK, ITS1 and ITS2. However, not all three barcodes were available for all species (Kress et al. 2005, Hollingsworth et al. 2009). When barcodes were not publicly available on GenBank but were available for a sister taxon sampled in our study, we included the sister taxon and added the missing species post hoc as tip polytomies. If gene information was missing and a given species did not have a representative sister species in the phylogeny, we looked for sequences for another regionally occurring species from that genus (occurrence was based on the Flora of Arunachal Pradesh, Hajra et al. 1996). We were unable to locate information on congenerics for three species (Balakata baccata, Khasiaclunea oligocephala, and Oxyspora paniculata) and thus included these taxa as polytomies to their closest relatives present in the phylogeny (Ostodes paniculata, Breonia oligocephala, and Melastoma malabathricum, respectively) based on the APG3 phylogeny (Bremer et al. 2009). This iterative process allowed us to generate a DNA matrix for 206 of the 279 species in the regional pool. The final DNA sequences used for constructing the phylogeny are listed in Appendix S2: Table S1.
Sequences were aligned using MAFFT ver. 7 (Katoh and Standley 2013) and trimmed using BioEdit (Hall 1999). We then concatenated the sequences for the different gene regions using SequenceMatrix ver. 1.7.8 (Vaidya et al. 2011) yielding a combined matrix with a length of 4365 bases. We inferred the phylogeny in MrBayes ver. 3.2.2 (Ronquist and Huelsenbeck 2003) by partitioning the data for each sequence and assigning the appropriate evolutionary model, as determined by modelTest in the phangorn R library (Schliep 2011). The genes rbcL, ITS1, and ITS2 were assigned the GTR+G+I model, while matK was assigned the GTR+G model. The phylogeny was constrained at the order or family level by assigning species to their known clades following the APG3 phylogeny (Bremer et al. 2009). We ran 25 million generations and excluded the first 25% as burn-in. One hundred phylogenic trees were randomly selected from the posterior distribution and rooted on the fern, Angiopteris evecta. Each tree was then made proportional to time, with the function chronopl from the R package ape (Paradis et al. 2004 . Missing taxa were included at this stage as polytomies, using taxonomy as a guide, with the function add.species.to. genus from the R package phytools (Revell 2012). The resulting phylogenic trees thus included all 279 taxa from the region. All subsequent phylogenetic analyses were conducted on these 100 phylogenic topologies. A sample phylogeny is included in Appendix S3: Fig. S1.

Analysis of plant communities
For each plot, we calculated species richness, phylogenetic diversity (Faith 1992), and the netrelatedness index (Webb et al. 2002) using the R library picante (Kembel et al. 2010). We calculated both Faith's PD (Faith 1992), which is equal to the sum of branch lengths represented in a community, and the standardized effect size of phylogenetic diversity which corrects for species richness with a tip-swap algorithm assuming the regional phylogeny (279 species) as the species pool. The net-relatedness index (NRI) incorporates evolutionary information from the phylogeny to calculate the average relatedness of species within a community relative to a null expectation of random community assembly and is independent of species richness. We calculated both abundance weighted and non-weighted NRI using the same null model as for standardized effect size of phylogenetic diversity. We used linear regression to explore how these metrics varied with elevation, which was normalized with a log-transformation.
Next, we calculated two pairwise measures of beta diversity. First, we used Sorensen's index to contrast species composition between plots using the vegdist function from the vegan R library (Oksanen et al. 2007). Second, we used a phylogenetic equivalent of the Sorensen's index to calculate phylogenetic beta diversity between plots using the phyloSor function in picante, which quantifies the proportion of shared branch lengths (Bryant et al. 2008, Graham andFine 2008). The phyloSor function returns a pairwise similarity matrix which can easily be converted into a dissimilarity matrix when appropriate (1-similarity). We also calculated standardized effect sizes of beta diversity for each of the 100 phylogenetic topologies using the function phylosor.rnd, in picante, to control for the effect of taxonomy in our analysis (Bryant et al. 2008, Kembel et al. 2010. A more detailed description of this calculation and the results can be found in Appendix S4. We then determined the local contributions to regional beta diversity (LCBD) using the method of Legendre and De C aceres (2013). R-code for implementing this function is available from http://adn.biol.umontreal.ca/~numericalecology/ Rcode/. This metric identifies plots with unique or unusual composition. To identify individual species that contribute disproportionately to beta diversity, we used a feature of LCBD, which reports species contributions to beta diversity (SCBD) by identifying species with high abundances in relatively few sites (Legendre and De C aceres 2013). Because we were also interested in phylogenetic patterns, we used a simple extension of this metric to estimate phylogenetically informed LCBD (herein referred to as PLCBD) by using the phylogenetic beta diversity distance matrix in place of the Euclidean distance matrix of species compositional dissimilarities (phyloSor outputs converted to represent dissimilarities for this analysis). We did not calculate significance values for PLCBD due to the extensive computational requirements associated with iterating across 100 separate phylogenetic trees, and we were more interested here in the overall patterns of PLCBD across the landscape rather than the statistical significance of any particular plot. We note that other useful methods also allow the exploration of the cladewise contribution to phylogenetic structure (Borregaard et al. 2014), although our focus here was on the effects of a few, phylogenetically disparate, species on community structure.
We explored patterns in beta diversity by first contrasting Sorensen's index with the phylogenetic ❖ www.esajournals.org equivalent using a mantel test. We next used partial mantel tests to separately explore the relationship between both taxonomic and phylogenetic beta diversity and distance (space) or elevation (environment). We then explored the relationship between plot contributions to beta diversity (LCBD and PLCBD) and elevation using simple linear regression. Finally, we assessed how the relationship between local contributions to beta diversity (LCBD and PLCBD) and environment (elevation) changed across space using geographically weighted regressions (GWR). This method can reveal spatial structure in regression coefficients by weighting data such that data points that are closer in space have a higher weight in the regression. Geographically weighted regression models were fit assuming adaptive weights in the spgwr package (Bivand et al. 2017). Additional data on altitude were extracted from the World-Clim database (Hijmans et al. 2005).
All analyses were performed using R ver. 3.0.2. (R Core Team 2015).

Overall patterns of diversity
Both species richness and phylogenetic diversity decreased with increasing elevation (SR: R 2 adj = 0.285, P < 0.001; PD: R 2 adj = 0.215, P < 0.001). Standardized effect size of PD increased with elevation, but the relationship was weaker, and plots with significantly higher PD than expected were located throughout the landscape (Appendix S5). In addition, we observed a significant, but weak, negative relationship between phylogenetic dispersion (indexed by NRI) and elevation (R 2 adj = 0.133, P < 0.001 and R 2 adj = 0.132, P < 0.001 for unweighted and weighted NRI, respectively). Thus, phylogenetically clustered communities were marginally more often found at low elevations and communities became increasingly over-dispersed at higher elevations, contrary to our initial hypotheses. We found some evidence for the opposite relationship when gymnosperms and ferns were removed from the analysis, but the correlations were weaker (non-weighted: R 2 adj = 0.047, P < 0.001; weighted: R 2 adj = 0.039, P < 0.001; Appendix S6: Fig. S1). A mantel test of the Sorensen's dissimilarity matrix and the phylosor similarity matrix (transformed to dissimilarity) revealed a strong relationship between the pairwise beta diversity metrics (mantel r = 0.703, P = 0.001). Plotting the two distance matrices revealed that species turnover occurred at a faster rate than phylogenetic branch turnover (Appendix S7: Fig. S1). We also observed a strong relationship between taxonomic beta diversity and elevation, which remained significant when we corrected for differences in geographical distance among plots (Table 1). However, the relationship between phylogenetic beta diversity and elevation was stronger, both with and without correcting for distance among plots (Table 1); results were quantitatively similar using standard effect sizes (Appendix S4). We found similar patterns when gymnosperms and ferns were removed from the analyses (Appendix S8: Table S1).
Sensitivity analyses demonstrated that the relationships between PD, NRI, and phylogenetic beta diversity with elevation were not sensitive to the proportion of individuals identified per plot (Appendix S1: Table S1).

Local contributions to beta diversity
The site with the lowest LCBD (least distinct community) was located near the geographical center of the study site ( Fig. 2a; 94.704°E, 27.727°N ). The site with the lowest PLCBD was also fairly central, but slightly to the east ( Fig. 2b; 95.648°E, 28.249°N). Both sites were at relatively low elevations (LCBD min = 528 m, PLCBD min = 295 m). In general, plots located on the periphery of study area tended to have higher contributions to beta diversity than plots in central, low-elevational sites (Fig. 2a, b). We found a significant, relationship between the strength of contribution to beta diversity and elevation, with plots contributing more at higher elevations (R 2 adj = 0.279, P < 0.011; Fig. 2c). A similar trend was found for the phylogenetic Note: Distance represents the pairwise geographical distance between each plot, while elevation represents the pairwise elevational differences between plots. local contribution to beta diversity (PLCBD; R 2 adj = 0.163, P < 0.001; Fig. 2d). The relatively low r-squared can, in part, be explained by the triangular relationship in the data, with plots at lower elevations having higher variance in their contribution. The distribution of both LCBD and PLCBD was qualitatively similar when gymnosperms and ferns were removed from the analyses (Spearman rank correlations; LCBD = 0.939, PLCBD = 0.853). The taxa which contributed the most to beta diversity, as indexed by the SCBD (Legendre and De C aceres 2013), were Castanopsis indica, Duabanga grandiflora, Quercus sp., Pinus roxburghii, and Musa sp.; these taxa are restricted in their distribution but have high local abundances.
The GWR models significantly improved the fit of the relationship between LCBD and elevation (Quasi-global R 2 = 0.7912, AIC GWR = À4142.45, AIC LM = À3861.72), with high-elevation plots showing a stronger positive relationship with elevation and low-elevation plots showing weaker and sometimes negative relationships with elevation (Fig. 3a, b). In contrast, the correlation between PLCBD and elevation was weaker (Quasi-global R 2 = 0.3059), and the improvement in fit of the GWR model compared to the unweighted linear model was less pronounced (AIC GWR = À3722.05 vs. AIC LM = À3679.64), suggesting less spatial variation in model fits (Fig. 3c).

DISCUSSION
We explored shifts in tree community structure and richness across one of the largest elevation gradients in the world, the Himalayas of Arunachal Pradesh, India. We found that species richness and phylogenetic diversity declined with Fig. 2. Maps show the spatial distribution of local contribution to beta diversity (LCBD) (a) and local contribution to phylogenetic beta diversity (PLCBD) (b). For (a) and (b), symbols are shaded by contribution, where red indicates higher contributions to beta diversity. The plots with the lowest contributions are colored white and identified with arrows. These plots represent the least unique sites for LCBD (a) and PLCBD (b), respectively. We also show the change in LCBD (c) and PLCBD (d) of each plot with increasing elevation. elevation, a result that is consistent with expectations and existing ecological theory. In general, elevational declines in richness are hypothesized to be due to factors similar to those driving the decline in richness observed along the latitudinal gradient, such as the reduced availability of resources, colder temperatures, and increased extinction rates at regional scales (Lomolino 2001, McCain andGrytnes 2010). A reduction of resources (lush soils and nutrients, for example) and colder temperatures at high elevations can limit the number of individuals and select for species with specific niche attributes (McCain and Grytnes 2010), with only those species possessing the appropriate traits and adaptations able to establish and thrive in these environments.
Several lines of evidence in our study suggest that environmental filtering is contributing to Fig. 3. Maps of the spatial relationships between elevation (shown for the region in a), local contribution to beta diversity (LCBD), and local contribution to phylogenetic beta diversity (PLCBD), calculated with geographically weighted regressions (GWR). LCBD is plotted at each study site as a solid symbol shaded by the strength of contribution (b), where red symbols represent stronger contributions to beta diversity (as with Fig. 2). PLCBD is similarly plotted for each site in (c). The map colors in (b) and (c) indicate the spatial variation in GWR coefficients for elevation modeled against LCBD and PLCBD, respectively. Red background colors represent positive coefficients, while blue background colors represent negative coefficients (refer to legend for b and c).
shifts in community structure with elevation, including high local endemism and rapid turnover with changes in elevation. Although we found that species turnover was more rapid than turnover of phylogenetic branch lengths, which has previously been interpreted as evidence for niche conservatism (Jin et al. 2015), we did not find any evidence for phylogenetic clustering (sensu Webb et al. 2002), which would also be expected with high conservatism and strong filtering. The absence of strong evidence for phylogenetic clustering indicates that these filtered taxa do not represent recent radiations within one or a few clades, but instead are dispersed across the phylogeny.
There are two broad, but not mutually exclusive, explanations for the lack of pronounced phylogenetic structuring along the elevational gradient despite evidence of strong species filtering. First, important traits could demonstrate convergent evolution such that distant relatives share similar ecological habitats. Second, in high montane regions, filtering may operate on evolutionarily distinct lineages with shared traits that are remnants of once more diverse cold-adapted clades. Much previous work has focused on the former (Jobb agy and Jackson 2000, Kraft et al. 2007 (simulations), Losos 2011, Read et al. 2014); here, we explore the latter and consider the phylogenetic evidence for refugia, perhaps maintaining relict lineages that were more widespread at the last glaciation, and structuring communities in the high elevations of the Himalayas.
Areas of high topographic relief like mountain ranges have been linked to the presence of climate refugia (Weber et al. 2014) because they provide cooler, more stable climates during warming periods (Stewart et al. 2010). These refugia provide suitable conditions for species that have retreated to microclimates resembling those of the last glacial maxima (Vetaas andGrytnes 2002, Ohlem€ uller et al. 2008). Such refugia might be increasingly important for many species given current warming trends (Opgenoorth et al. 2010). However, identifying relict communities or species is challenging, usually requiring detailed population genetics on a regional scale to reconstruct patterns of migration (Hampe et al. 2003, Petit et al. 2005, Vargas 2007). Here, we draw inference from the phylogenetic structure of high-elevation communities.
Our approach combines knowledge of the evolutionary relationships among species with information on shifts in community composition with elevation, allowing us to identify diversity patterns that might reflect the distribution of relict lineages and climate refugia. For example, if glacial relicts represent survivors from once more diverse clades, perhaps a result of higher extinction rates of related species (Cain 1944, Fryxell 1962, Brooks and Bandoni 1988, we might expect them to have a large impact on the phylogenetic structure of local communities. The presence of glacial relicts would not only disrupt patterns of phylogenetic clustering predicted at high elevations in strongly filtered communities, but also contribute to the uniqueness or beta diversity of those communities. Although we did not find strong evidence for higher phylogenetic diversity within higher elevation plots in Arunachal Pradesh, we show that high-elevation plots do indeed contribute disproportionately to regional beta diversity. Because highly contributing plots are those that contain communities with relatively greater species uniqueness (Legendre and De C aceres 2013), this would be consistent with the presence of narrow ranged and evolutionarily distinct endemics at higher elevations.
While we found a general trend for plots with high contribution to beta diversity to be at higher elevations, the relationship between elevation and contribution to beta diversity varies across the landscape. Notably, we observe large variation at lower elevations, and within low elevations, there is no consistent relationship between local contribution and elevation. At these elevations, we suggest that there may be a greater incidence of landscape modification and anthropogenic influence (Menon et al. 2001, Bhuyan et al. 2003, Roy and Behera 2005, which may have reshaped community structure. For example, species with high individual contributions to beta diversity more commonly found at low elevations include Castanopsis indica and Musa sp. and may represent naturalized species cultivated for their useful properties (timber, food). In contrast, as we move to higher elevations, there is a more predictable relationship between local contribution and elevation, consistent with the process of environmental filtering, such that higher elevation plots represent a more unique subset of regional diversity. Species with high individual contributions to beta diversity found at mid-to high elevations include Pinus roxburghii (Puri et al. 2011, IUCN RedList 2015 and Pinus wallichiana (Saqib et al. 2013, IUCN 2015, which are endemic to the region. The relationship between plot level contributions to local phylogenetic beta diversity (PLCBD) and elevation was weaker and showed little spatial variation, suggesting that the presence of non-indigenous naturalized species at lower elevations does not skew the phylogenetic structure of plant communities as much as they influence their taxonomic composition. The trend for more phylogenetically unique communities at higher elevations is again consistent with increasing filtering of evolutionarily distinct species as we go up in elevation. We suggest these evolutionarily distinct species might represent glacial relicts: Species from once more widespread communities adapted to the colder climates that characterized the last glacial cycle, but which have subsequently contracted their ranges toward high elevations.

SUMMARY
Despite observing a general decrease in species richness and phylogenetic diversity with elevation, we did not find any evidence for increased phylogenetic clustering along the elevational gradient, which is typically expected if communities are being filtered by environment. Communities at higher elevations tended to show greater dispersion. However, we showed that species and phylogenetic turnover increased with elevation and that high-elevation communities were more distinct. We therefore suggest that environmental filtering is an important process structuring communities in this region, but that the species filtered into high-elevation communities are not close relatives. We propose these communities might represent assemblages of evolutionary distinct lineages that were adapted to more widespread conditions during the last glaciation. These high-elevation plots may thus act as "inter-glacial refugia." Our analysis is first to look at community structure along this extensive elevational gradient, and although the strong effect of environmental filtering may not be surprising, we show additional evidence to suggest that high-elevational communities in this region are particularly unique, and argue that they should be targets for conservation in the future.
A better understanding of richness patterns ultimately requires researchers to collect data in isolated, overlooked, and hard to access regions around the world. We suggest that regions with unique species, high endemicity, and distinct geography should become priorities for research and conservation. By understanding the historical factors that have shaped them, these communities might provide insights into responses to future environmental change, not just at the individual species level, but also at the level of the ecological assemblage. High-altitude refugia, which we identify here, may be important conservation targets because they can provide an escape from generally increasing temperatures globally by matching to the cooler climates resembling the conditions under which many taxa may have evolved.

ACKNOWLEDGMENTS
We thank our colleagues at NERIST for collecting, compiling, and sharing these data. We also thank the Davies Lab Group for their valuable input. This study was supported through an Excellence Award from the Quebec Center for Biodiversity Science, as well as a Merit Award from Concordia University to Stephanie Shooner. Additional support was provided by NSERC Discovery Grants to both Jonathan Davies and Selvadurai Dayanandan and funding was provided by the Department of Biotechnology (Government of India) for fieldwork. Baishya, A. 1999