SPATIOTEMPORAL DYNAMICS OF GIANT PANDA HABITAT: IMPLICATIONS FOR PANDA CONSERVATION UNDER A CHANGING ENVIRONMENT By Mao-Ning Tuan Mu A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Fisheries and Wildlife 2012 ABSTRACT SPATIOTEMPORAL DYNAMICS OF GIANT PANDA HABITAT: IMPLICATIONS FOR PANDA CONSERVATION UNDER A CHANGING ENVIRONMENT By Mao-Ning Tuan Mu Under the current rapidly changing environment, effective and efficient actions for biodiversity conservation rely on detailed knowledge on the spatiotemporal dynamics of species distribution and habitat. However, inadequate spatiotemporal information on species habitat has compromised conservation effectiveness, even for one of the most endangered species on Earth, the giant panda (Ailuropoda melanoleuca). To address this information gap, the objectives of this dissertation were to: (1) develop an approach for remotely detecting the distribution of understory bamboo, the panda’s staple food, across large geographic regions; (2) develop a modeling approach for monitoring panda habitat changes across space and time; (3) evaluate the effects of current conservation efforts on short-term panda habitat changes; and (4) assess the potential impacts of climate change on long-term panda habitat dynamics. Using two dominant bamboo species in Wolong Nature Reserve, China, I showed that an integration of species distribution modeling with land surface phenology obtained from high temporal resolution remotely sensed data is a promising approach for providing detailed information on understory bamboo distribution across large geographic regions. Derived from time series data acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS), eleven land surface phenology metrics successfully captured the phenological characteristics of vegetation caused by understory bamboo. In addition, a species distribution model (SDM) built using the maximum entropy modeling approach (Maxent) accurately captured the distribution of understory bamboo species across the reserve based on their phenological characteristics. I further demonstrated the usefulness of the phenology-based model for not only characterizing panda habitat across space, but also monitoring its dynamics over time. By quantitatively examining the effects of different predictor variables portraying land surface phenology on a model’s ability to be reliably applied to different time periods, I showed that a model built with phenology metrics derived from multi-year remotely sensed data had the highest temporal transferability. Based on the outputs of this model, I evaluated the effectiveness of a conservation program by investigating the spatiotemporal dynamics of panda habitat from 2001 to 2007 in Wolong Nature Reserve. Results suggested that an innovative implementation of the Natural Forest Conservation Program, which encourages active participation of local residents in forest monitoring by providing direct payments and enhancing social norms among households, is an effective instrument for panda habitat conservation. While the current conservation programs have effectively reduced the threats of land use/cover change to panda habitat, I also showed that climate change may become the next major threat to the survival of the giant panda. Focusing on the food resources of the panda population in the Qinling Mountains region, an ensemble of panda habitat projections obtained from bioclimatic envelope models indicated a substantial loss of panda habitat due to a potential shortage of food under projected climate change in the 21st century. This poses a big challenge for panda conservation in the face of climate change and suggests an urgent need for developing proactive conservation practices. This dissertation makes substantial contributions to giant panda conservation by providing useful tools and essential information for a better understanding of the spatiotemporal dynamics of panda habitat. The findings from this dissertation also have broad implications for biodiversity conservation under a changing environment. Copyright by Mao-Ning Tuan Mu 2012 ACKNOWLEDGEMENTS I would like to express my deep appreciation to all the people and organizations that made my dissertation research possible. First, I would like to thank my committee, Dr. Jack Liu, Dr. Gary Roloff, Dr. Ashton Shortridge, Dr. Andrés Viña and Dr. Julie Winkler for their guidance and advice throughout my doctoral study. I want to especially thank Andrés for always believing in me and encouraging me, and for leading me to the study field of remote sensing. In particular, I would like to offer my sincerest thanks to Jack, my advisor. I would not have achieved anything without his continuous help, support and patience, especially when I was struggling with my dissertation research in the first two years of my PhD program. I want to give thanks to the people who worked with me in remote mountains in China, including Yu Li, Wei Liu, Andrés Viña and Weihua Xu. In particular, I would not have survived the 2008 Wenchuan Earthquake without help and support from Wei and many people in Wolong Nature Reserve. I also would like to thank Gaodi Dang, Xichun Du, Junzheng Meng and Yong Wang for their help with plant species identification. I really appreciate the Forestry Department of Shaanxi Province, Forestry Department of Sichuan Province and administrations of nature reserves in the Qinling, Minshan and Qionglai Mountains for support with fieldwork logistics. Of course, it is the assistance of Zhiyun Ouyang and Hemin Zhang that made my field trips possible. In addition, I thank Scott Bearer, Xiaodong Chen, Vanessa Hull, Wei Liu, Andrés Viña and Wu Yang for sharing with me their valuable data and ideas. I also really appreciate help and support from all the members of the Center for Systems Integration and Sustainability, Department of Fisheries and Wildlife, Michigan State University. Of course, I would not have V been able to concentrate on my research without the support from my always patient and helpful wife, Chin-Mei Lee, and from our families in Taiwan. Finally, I would like to gratefully acknowledge the financial support from the U.S. National Science Foundation (Dynamics of Coupled Natural and Human Systems and Partnership for International Research and Education), the National Aeronautics and Space Administration (Land Use/Land Cover Change program, Terrestrial Ecology and Biodiversity program, and Earth and Space Science Fellowship), Michigan Agricultural Experiment Station, and the National Natural Science Foundation of China. VI TABLE OF CONTENTS LIST OF TABLES……………………………………………………………………………....IX LIST OF FIGURES………………………………………………………………………………X CHAPTER 1 BACKGROUND AND RESEARCH QUESTIONS……………………………………………1 Background………………………………………………………………………………2 Research Objectives……………………………………………………………………...11 Study Areas……………………………………………………………………………12 CHAPTER 2 MAPPING UNDERSTORY BAMBOO USING PHENOLOGICAL CHARACTERISTICS DERIVED FROM REMOTELY SENSED DATA……………………………………………...18 Abstract………………………………………………………………………………….19 Introduction……………………………………………………………………………..20 Materials and Methods……………………………………………………………......24 Field data…………………………………………………………………….…..24 Remotely sensed data and phenology metrics…………………………………..25 Phenological characteristics of forests with understory bamboo…………….....28 Overall bamboo distribution……………………………………………………..29 Individual bamboo distribution………………………………………………….31 Model validation and comparison……………………………………………….32 Results…………………………………………………………………………………..35 Phenological characteristics of forests with understory bamboo…………….....35 Overall bamboo distribution……………………………………………………36 Individual bamboo distribution…………………………………………………40 Discussion……………………………………………………………………………….44 Conservation Implications and Conclusions…………………………………………….50 CHAPTER 3 TEMPORAL TRANSFERABILITY OF WILDLIFE HABITAT MODELS: IMPLICATIONS FOR PANDA HABITAT MONITORING………………………………52 Abstract…………………………………………………………………………………53 Introduction………………………………………………………………………….....54 Materials and Methods………………………………………………………………….56 Giant panda presence data……………………………………………………56 Remotely Sensed data…………………………………………………………..57 Predictor variables……………………………………………………………….58 Analytical design……………………………………………………………......60 Model development……………………………………………………………..62 Habitat prediction…………………………………………………………….....63 VII Model evaluation………………………………………………………………..64 Model transferability…………………………………………………………….65 Model comparison……………………………………………………………….68 Results…………………………………………………………………………………..68 Multicollinearity and model complexity………………………………………..68 Model accuracy………………………………………………………………….69 Model transferability…………………………………………………………….71 Discussion……………………………………………………………………………….77 Effects of variable type on model transferability………………………….....77 Effects of time series length on model transferability………………………80 Usefulness of land surface phenology for mapping and monitoring species habitat………………………………………………………………………….…81 CHAPTER 4 INCENTIVE-BASED AND LOCAL-ENGAGED INSTRUMENT FOR SUCCESSFUL CONSERVATION OF PANDA HABITAT ……………………………………………………84 Abstract…………………………………………………………………………………..85 Introduction…………………………………………………………………………..…85 Materials and Methods……………………………………………..……………88 NFCP implementation in Wolong Nature Reserve………………………….88 Spatiotemporal dynamics of panda habitat……………………………………90 Effects of NFCP implementation…………………………………………….....91 Results………………………………………………………………………………….96 Discussion……………………………………………………………………………....98 CHAPTER 5 CASCADING EFFECTS OF CLIMATE CHANGE THROUGH TROPHIC RELATIONS: PROJECTED FOOD SHORTAGES THREATEN THE GIANT PANDA……………………102 Abstract…………………………………………………………………….…….……..103 Introduction…………………………………………………………………………....103 Materials and Methods………………………………………………………………106 Field data………………………………………………………………….…....106 Baseline climate conditions…………………………………………………....106 Future climate projections………………………………………………………107 Bioclimatic models……………………………………………………………..107 Climatically suitable areas for bamboo species………………………………..112 Extent of giant panda habitat…………………………………………………..113 Results………………………………………………………………………………......114 Discussion………………………………………………………………………………124 CHAPTER 6 SYNTHESIS AND CONCLUSIONS………………………………………………………130 Summary……………………………………………………………………………….131 Conclusions……………………………………………………………………………..134 REFERENCES…………………………………………………………………………………139 VIII LIST OF TABLES Table 2.1. Comparisons of the performance of models developed for individual bamboo species with and without elevation information, based on the results of 20 model runs………………...42 Table 3.1. Properties of habitat models developed with four different sets of land surface phenology variables in two time frames…………………………………………………………59 Table 3.2. Analyses of variance on the effects of variable type (Type) and the length of original time series data (Length), as well as their interaction effects (Type × Length) on the singledirection and overall transferability indices calculated for the habitat models developed with different land surface phenology variables……………………………………………………....74 Table 3.3. Analyses of variance on the effects of variable type (Type) and the length of original time series data (Length), as well as their interaction effects (Type × Length) on agreement coefficients (AC) and proportions of systematic disagreement (PSD), which were calculated from within- and beyond-time-frame predictions of panda habitat in the two time frames……..76 Table 4.1. Independent variables included in regression models that relate the changes in panda habitat suitability index (HSI) with the NFCP implementation and other biophysical and anthropogenic factors…………………………………………………………………………….93 Table 4.2. Summary of regression models that relate the changes in panda habitat suitability index (HSI) with a set of independent variables…………………………………………………94 Table 5.1. General circulation models (GCMs) from which downscaled future climate projections were obtained. GCMs from the IPCC Third Assessment were used for evaluating temporal dynamics of bamboo distribution and panda habitat over the 21st century. The other models are from the IPCC Fourth Assessment and were used to better capture uncertainty introduced by different GCMs……………………………………………….……………..…..109 Table 5.2. Bioclimatic variables used in the final bioclimatic models for the three bamboo species studied………………………………………………………………………………….111 IX LIST OF FIGURES Figure 1.1. Locations and topography of the Wolong Nature Reserve, Sichuan Province and the Qinling Mountains region, Shaanxi Province, China……………………………………………14 Figure 2.1. Locations of the field plots where arrow or umbrella bamboo cover was 25% or higher according to the Third National Giant Panda Survey conducted in Wolong Nature Reserve in 2001………………………………………………………………………………………...…25 Figure 2.2. Phenology metrics (A – K) derived from the smoothed values (triangles) of a time series of WDRVI values. A - Base level; B - Maximum level; C - Amplitude; D - Date of the start of the season; E - Date of the middle of the season; F - Date of the end of the season; G Length of the season; H - Large integral; I - Small integral; J - Increase rate; K - Decrease rate……………………………………………………………………………………………….27 Figure 2.3. Box plots of 11 phenology metrics calculated using 1,000 pixels randomly selected from the entire study area (All), 1,000 pixels randomly selected from forested areas (For), 356 pixels where field plots with bamboo cover as 25% or higher were located (Bam), 145 pixels where field plots with arrow bamboo cover as 25% or higher were located (Arr), and 184 pixels where field plots with umbrella bamboo cover as 25% or higher were located (Umb). The dark line inside a box indicates the median; the bottom and the top of a box show the 25th and 75th percentiles, respectively; the ends of the two whiskers indicate the lowest and the highest values within 1.5 inter-quartile ranges from a box; dots beyond whiskers show outliers whose values were outside the range indicated by the whiskers. The letters above boxes show the results of pair-wise comparisons on the phenology metric values conducted using Mann-Whitney U-tests. Two boxes share the same letter if there was no significant difference (p > 0.001) between them………………………………………………………………………………………………37 Figure 2.4. Overall bamboo distribution across Wolong Nature Reserve. The pixel-wise (a) mean values and (b) standard deviations of bamboo presence probabilities were calculated over 20 runs of the overall bamboo distribution model. The model was generated using 356 pixels with 25% or higher bamboo cover as presence locations and 11 phenology metrics as predictor variables. The 289 pixels where phenology metrics could not be determined in at least two years between 2001 and 2003 from a smoothed curve of a time series of WDRVI values by TIMESAT are represented in white………………………………………………………………………….39 Figure 2.5. Mean AUC values of the overall bamboo distribution models (a) with only one of the 11 phenology metrics and (b) with all, but one, metrics as predictor variables. A higher AUC value for a model with only one metric indicates that the metric contained more useful information for mapping bamboo distribution in the full model. A lower AUC value (larger loss of the AUC value) for a model without one metric indicates that the metric contained more information which cannot be represented by the other metrics for mapping bamboo distribution. Max - Maximum level; Amp - Amplitude; SOS - Date of the start of the season; MOS - Date of the middle of the season; EOS - Date of the end of the season; L_Int - Large integral; S_Int Small integral; I_Rate - Increase rate; D_Rate - Decrease rate………………………………….41 X Figure 2.6. Spatial distribution of arrow and umbrella bamboo across Wolong Nature Reserve. Mean presence probabilities of pixels were calculated over 20 runs of the individual bamboo distribution models containing only 11 phenology metrics (a and b) and the models containing the11 phenology metrics plus elevation (c and d). Pixels where phenology metrics could not be determined in at least two years between 2001 and 2003 from a smoothed curve of a time series of WDRVI values by TIMESAT are represented in white…………………………….………...43 Figure 3.1. Locations of panda activity signs recorded during the Third National Giant Panda Survey in 2001 and during wildlife surveys from 2006-2008…………………………………...57 Figure 3.2. The analytical design included three steps: (1) model development (arrow with double lines), (2) habitat prediction (arrow with single solid line), and (3) model evaluation (arrow with dashed line). (1) Panda habitat models were developed with four different sets of land surface phenology (LSP) variables in two time frames (t1 and t2) using panda presence data collected in 2001 and 2006-2008, respectively. (2) The models were used to predict panda habitat within (WTP) and beyond the time frame (BTP) in which the models were developed. (3) The habitat maps from WTP and BTP were evaluated using the presence data collected in the time frame in which the habitat was predicted. AUC and MPA were calculated for the habitat maps from both WTP (AUCt1t1, AUCt2t2, MPAt1t1, and MPAt2t2) and BTP (AUCt1t2, AUCt2t1, MPAt1t2, MPAt2t1). The habitat maps from WTP and BTP were compared (arrow with dasheddotted line) within each time frame, and values of agreement coefficients (ACt1 or ACt2) and the proportions of systematic disagreement (PSDt1 or PSDt2) were calculated. The AUC, AC and PSD values were then used to evaluate model transferability………………………………...…61 Figure 3.3. Box plots of the area under the receiver operating characteristics curve (AUC) and the minimal predicted area (MPA) for the panda habitat models developed with four different land surface phenology variable sets (SYVI, MYVI, SYPM and MYPM in Table 3.1) in two time frames (t1 and t2). The AUC and the MPA values were calculated from 20 variants of each habitat model when the model was developed in t1 and used to predict habitat in t1 [(a) and (e)] and t2 [(b) and (f)], as well as developed in t2 and used to predict habitat in t2 [(c) and (g)] and t1 [(d) and (h)]. Each box plot shows the maximum, 75th percentile, median, 25th percentile and minimum values. Note that y-axes for the MPA are reversed since lower MPA values indicate higher model accuracies. The letters above box plots indicate the results of pair-wise comparisons conducted using paired t-tests. The alphabetical order shows the order of model accuracy from high to low. There is no significant difference (p > 0.05 after the Holm– Bonferroni adjustment) between two models if they share the same letter……………………..70 Figure 3.4. Mean values of the single-direction and overall transferability indices for the four habitat models. These models were developed with either a time series of WDRVI (VI) or phenology metrics (PM), which were derived from either single-year (SY) or multi-year (MY) MODIS data. The indices were calculated from 20 variants of each habitat model for evaluating: (a) the ability of the model developed in Time Frame 1 (t1) to predict panda habitat in Time Frame 2 (t2), (b) the ability of the model developed in t2 to predict habitat in t1, and (c) the overall ability of the model to predict habitat beyond time frames. The error bars indicate 2 standard errors from the mean. The lines do not imply any linear relationship, but are shown just for helping visualise value differences…………………………………………………………..73 XI Figure 3.5. Mean values of agreement coefficients (AC) and the proportions of systematic disagreement (PSD) for the four habitat models. These models were developed with either a time series of WDRVI (VI) or phenology metrics (PM), which were derived from either singleyear (SY) or multi-year (MY) MODIS data. The values of AC and PSD were calculated between the maps generated from within- and beyond-time-frame predictions for the panda habitat in Time Frame 1 ((a) and (c), respectively) and in Time Frame 2 ((b) and (d), respectively). The error bars indicate 2 standard errors from the mean. The lines do not imply any linear relationship, but are shown just for helping visualise value differences……………..75 Figure 4.1. Spatial patterns of giant panda habitat change from 2001 to 2007 in Wolong Nature Reserve, China. Changes in habitat suitability index (HSI) were calculated from the outputs of the panda habitat model built with multi-year phenology metrics in Chapter 3. The locations of households and household-monitored forest parcels in the three townships (Gengda, Wolong and Sanjiang) comprising the reserve are also shown………………………………………………..89 Figure 4.2. Effects of the implementation of the Natural Forest Conservation Program (NFCP) on giant panda habitat change in Wolong Nature Reserve, China. Boxplots show the partial residuals [controlling for biophysical and anthropogenic factors (Table 4.1)] for two different NFCP monitoring types [government- (red) vs. household-monitoring (blue)]. Plotted values are the residuals of the full model plus the partial predicted values (black lines) in overall (a) and under two different NFCP payment levels (b). Results of Student’s t-tests on the difference N.S. p > 0.05)………………………….....…97 between monitoring types are shown (* p < 0.05; Figure 5.1. Location of the study area relative to current giant panda distributional range and the baseline climatically suitable areas (CSAs) for the three bamboo species studied. Blue, red, and green colors indicate the CSAs for Qinling arrow, dragon-head, and wooden bamboos, respectively, based on the bioclimatic models under the baseline climate. The mixtures of the three colors indicate overlaps of the CSAs for individual species. The brightness of colors shows the number of bioclimatic models (among the 10 models with different presence data partitions) predicting pixels as suitable for each species, with brighter colors indicating a larger number of models………………………………………………………………………………………….115 Figure 5.2. Projected elevational distributions of the climatically suitable areas (CSAs) for the three bamboo species studied. The elevation values were obtained for pixels which were projected to be climatically suitable for the wooden (a and d), dragon-head (b and e), and Qinling arrow bamboo (c and f) by at least one of 10 bioclimatic models with future climate projections for three time slices in the 21th century from four IPCC TAR GCMs (Table 5.1) under the SRES A2 (a – c) and B2 (d – f) greenhouse gas emissions scenarios. Each box plot shows the maximum, 75th percentile, median, 25th percentile, and minimum values. The grey zone on the background shows the inter-quartile range (75th percentile – 25th percentile) of elevation for the baseline CSAs (the CSAs under the baseline climate). The white line within the grey zone indicates median values and black dashed lines indicate the maximum and minimum values of elevation for the baseline CSAs………………………………………………………………..116 XII Figure 5.3. Projected future distributions of climatically suitable areas (CSAs) for the three bamboo species studied under the climate projections from four IPCC TAR GCMs (Table 5.1) for three time slices under SRES A2 and B2 greenhouse gas emissions scenarios. Blue, red, and green colors indicate the CSAs for the Qinling arrow, dragon-head, and wooden bamboos, respectively. See the legend of Figure 5.1 for the detailed information on the color representation…………………………………………………………………………………..118 Figure 5.4. Projected future distributions of climatically suitable areas (CSAs) for the three bamboo species studied under the climate projections from 15 IPCC AR4 GCMs. The climate projections were from the GCMs (Table 5.1): (a) BCCR-BCM2, (b) CGCM3.1 (T47), (c) CGCM3.1 (T63), (d) CNRM-CM3, (e) CSIRO-Mk3.0, (f) ECHO-G, (g) FGOALS-g1.0, (h) GFDL-CM2.0, (i) GFDL-CM2.1, (j) GISS-AOM, (k) IPSL-CM4, (l) MIROC3.2 (hires), (m) MIROC3.2 (medres), (n) MRI-CGCM2.3.2, and (o) PCM for the time slice of 2040 – 2069 under the SRES A2 greenhouse gas emissions scenarios. See the legend of Figure 5.1 for the detailed information on the color representation………………………………………………………...121 Figure 5.5. Temporal dynamics of the projected changes in the area of giant panda habitat over the 21st century. The changes (%) are relative to the area of panda habitat (climatically suitable areas for the three bamboo studied combined) under the baseline climate. The projected values were obtained from 1,000 combinations of the bioclimatic models for the three bamboo species (10 models for each species) under the unlimited (a and c) and no (b and d) bamboo dispersal scenarios and under multiple future climate projections for three time slices from four IPCC TAR GCMs (Supplementary Table S1) under the SRES A2 (a and b) and B2 (c and d) greenhouse gas emissions scenarios. The points indicate median values of the projections from the 1,000 projections, and the vertical bars indicate the range of projected values…………….122 Figure 5.6. GCM-related uncertainty of projected changes in giant panda habitat area for the time slice of 2040 – 2069 under the SRES A2 greenhouse gas emissions scenario. The percent changes are relative to the area of panda habitat (climatically suitable areas for the three bamboo studied combined) under the baseline climate. The projected values were obtained from 1,000 combinations of the bioclimatic models for the three bamboo species (10 models for each species) under the unlimited (a) and no (b) bamboo dispersal scenarios and under multiple future climate projections from four IPCC TAR and 15 AR4 GCMs (Table 5.1). Each boxplot shows the maximum, 75th percentile, median, 25th percentile and minimum of the projected values…..123 Figure 5.7. Percentage of projected giant panda habitat within current nature reserves in the Qinling Mountains over the 21st century. The extent of panda habitat was obtained by combining the occupied portions of the projected future climatically suitable areas for the Qinling arrow, dragon-head, and wooden bamboo, under the unlimited bamboo dispersal scenario. The percentages are associated with the panda habitat projections from 1,000 combinations of the bioclimatic models for the three bamboo species (10 models for each species) under future climate projections from four IPCC TAR GCMs (Table 5.1) under the SRES A2 (a) and B2 (b) greenhouse gas emissions scenarios. The points indicate the median values of the 1,000 projections and the vertical bars indicate the maximum and minimum values. The horizontal solid lines indicate the median values and dashed lines indicate the maximum and minimum values for the percentages under the baseline climate………………………………125 XIII CHAPTER 1 BACKGROUND AND RESEARCH QUESTIONS 1 Background Loss of biodiversity and associated ecosystem services is one of the most serious crises that humanity is facing (IUCN 2011; Millennium Ecosystem Assessment 2005). Human activities are threatening biodiversity at local to global scales through over-exploiting natural resources, damaging and degrading species habitats, introducing exotic species, polluting the environment and changing global climate (Butchart et al. 2010; Soule 1991). Evidence has shown that the current species extinction rate is higher than the background rate obtained from fossil records and the sixth mass extinction may be happening due to human activities (Barnosky et al. 2011; Pimm et al. 1995). Human-induced species extinction, population declines, species invasions and distributional range shifts and collapse not only have dramatically changed biotic structure and interactions in biological communities, but also have altered ecosystem properties and functions (Hooper et al. 2005). These changes in turn have caused declines in human wellbeing by reducing many goods and services that biodiversity and ecosystems provide (Balmford and Bond 2005). With growing human population and resource demand, the trend of biodiversity loss and associated declines in human wellbeing are expected to continue or even become worse in the future (Millennium Ecosystem Assessment 2005; Pereira et al. 2010; Thomas et al. 2004). In the face of this crisis, conservation efforts to significantly reduce the rate of biodiversity loss are critical and urgently needed (Hoffmann et al. 2010). More than 100,000 protected areas, which cover ca. 13% of the planet’s land surface, have been established for conserving terrestrial ecosystems and the biodiversity they support (Chape et al. 2005; Jenkins and Joppa 2009). However, although conservation actions have beneficial effects on global biodiversity (Hoffmann et al. 2010), current efforts are still insufficient to offset the negative 2 environmental impacts caused by human activities (Butchart et al. 2010). Since conservation resources are limited, more effective and efficient conservation practices depend on a wise allocation of resources, which in turn relies on detailed knowledge on the spatial patterns and temporal dynamics of species distribution and biodiversity patterns (Whittaker et al. 2005). For example, based on knowledge on spatial patterns of biodiversity, systematic conservation planning (Margules and Pressey 2000) can prioritize the areas with high species richness (Myers et al. 2000), with good biodiversity representation (Olson and Dinerstein 1998), or sustaining keystone or umbrella species (Branton and Richardson 2011). In addition, observed and projected spatiotemporal dynamics of species distribution and biodiversity patterns are also essential for evaluating the efficacy of current conservation efforts (Araújo et al. 2004; Hoffmann et al. 2010) and for guiding their adjustment under a changing environment (Carvalho et al. 2011; Hannah et al. 2007). However, current knowledge on species distribution and biodiversity patterns is seriously inadequate because only a fraction of the species on Earth have been described so far (Mora et al. 2011) and very limited distributional information for those described species is available (Whittaker et al. 2005). Broad-scaled information on species distribution for well-known taxa (e.g., mammals, birds and some plant groups) is mostly available in the form of checklists at country or ecoregion levels. In addition, many range maps are based on expert opinion, and thus usually have very coarse spatial and temporal resolutions (Jetz et al. in press). While this coarse information is important for guiding broad-scale conservation planning (e.g., Ceballos and Ehrlich 2002), it may sometimes lead to incorrect judgments about conservation needs and lead to problematic decision making (Jetz et al. 2008). Fine-resolution data on species distribution can be obtained from specimen collections in museums and herbaria and from field surveys and 3 species inventories. However, availability of those data is usually restricted in space and time (Jetz et al. in press; Whittaker et al. 2005), and thus may hinder effective conservation of biodiversity. Science-based conservation planning for tackling the crisis of biodiversity loss is particularly important for China, one of the most biodiverse and populous countries in the world (Liu 2010; Liu and Raven 2010). Due to its vast geographic area, which covers diverse environmental conditions and ecosystem types, China has ca. 9% of all vascular plants and terrestrial vertebrates known in the world, and is estimated to contain ca. 1 million species (Liu et al. 2003b; Liu and Raven 2010). However, continuous increase in human population and demand for natural resources, as well as the exponential growth of the economy in the last three decades (Liu and Diamond 2008) have led to serious biodiversity losses (Liu and Raven 2010). In response to this crisis, the Chinese government has established more than 2,500 protected areas, which cover more than 15% of China’s territory (Liu and Raven 2010). Nevertheless, systematic conservation planning has not been used for the design and management of many of those reserves (Liu et al. 2003b) and the representativeness of current protected areas for the conservation of the diversity of ecosystem and vegetation types was not systematically evaluated until recently (Wu et al. 2011). Assessing the effectiveness of current protected areas and other conservation efforts in conserving biodiversity, monitoring biodiversity changes and understanding underlying drivers, and developing science-based conservation planning are urgently needed for biodiversity conservation in China (Liu et al. 2003b; Xu et al. 2009a). The inadequacy of information on the spatial patterns and temporal dynamics of species distribution and their habitats may even impede effective conservation of one of the most endangered species in the world, the giant panda (Ailuropoda melanoleuca). The giant panda 4 was once distributed across most of southern and eastern China (Loucks et al. 2001; Pan et al. 2001), but currently only ca. 1,600 individuals survive within six mountain regions at the edge of the Tibetan Plateau (State Forestry Administration 2006). Habitat loss and degradation due to human activities such as cultivation and timber harvesting were the major reasons for the drastic decline of wild panda populations and their habitat during the past several decades (Liu et al. 2001; Loucks et al. 2001). As a national treasure of China and an icon of biodiversity conservation, the giant panda has received unparalleled conservation resources. For example, more than 60 nature reserves have been established specifically for giant panda conservation. However, lack of detailed knowledge on the distribution of panda habitat prior to the selection and design of reserves has led to the omission of some important habitats from the protection of current reserves (Viña et al. 2010; Xu et al. 2006) and has led to less effective reserve design (Hull et al. 2011; Liu and Li 2008). One of the major challenges for accurately mapping panda habitat is the lack of detailed information on the spatial patterns of understory bamboo across large spatial extents. Because the giant panda is an extreme dietary specialist, with more than 99% of its diet being composed of bamboo species growing under forest canopies (Schaller et al. 1985), distribution of understory bamboo is one of the most important factors determining the quantity, quality and spatial distribution of panda habitat (Bearer et al. 2008; Schaller et al. 1985). Due to the lack of essential information on understory bamboo, previous evaluations of panda habitat usually used coarse-scale bamboo maps, or used forest cover as a surrogate for bamboo distribution, resulting in considerable overestimations of the amount of panda habitat and its carrying capacity (Linderman et al. 2005). 5 Remote sensing provides spatially continuous observations of the Earth’s surface across large geographic areas, and thus constitutes a suitable tool for synoptically measuring vegetation characteristics. In some cases, remote sensing can directly detect the presence and distribution of individual species or species assemblages (Turner et al. 2003). For example, spatial distribution of some tree species can be obtained from remotely sensed imagery with very high spatial resolutions (e.g., < 1 m) through detecting individual tree crowns, or from imagery with high spectral resolutions (e.g., > 200 spectral bands) through species-specific spectral signatures (Turner et al. 2003). However, remote detection of understory vegetation is challenging due to the interference of overstory canopies (Eriksson et al. 2006; Rautiainen et al. 2007). Although many efforts have been made to overcome this challenge through advanced remote sensing techniques and image processing approaches (Linderman et al. 2004; Wang et al. 2009a; Wang et al. 2009b), applying those approaches to large areas is still limited by the availability of the required remotely sensed imagery. An effective approach for remotely detecting understory bamboo distribution across large geographic regions is needed to provide the essential information for accurately measuring panda habitat. Another information gap for effective conservation of the giant panda relates to the temporal dynamics of panda habitat. Knowledge about the changes in the quantity, quality and spatial distribution of panda habitat and the underlying drivers of the changes is essential for identifying threats, assessing the effectiveness of current conservation efforts, and guiding conservation planning under a changing environment, such as adaptive habitat management (Holling 1978; Swaisgood et al. 2011). Although repeated field surveys can provide detailed information on habitat changes, they are time consuming, expensive, and labor intensive. Therefore, repeated field surveys are usually restricted to local scales and/or to low frequencies, 6 which reduce their ability to monitor habitat dynamics across large areas or over short time periods. Although remotely sensed data have been used to investigate panda habitat changes, previous studies were focused on land cover changes, mainly forest cover changes (Jian et al. 2011; Liu et al. 2001; Viña et al. 2007; Xu et al. 2009b). Forest cover is an important factor affecting not only the habitat use of the giant panda (Bearer et al. 2008; Liu et al. 2005b; Zhang et al. 2011), but also the growth and distribution of understory bamboo (Qin et al. 1993; Taylor and Qin 1987; Wang et al. 2009c). Therefore, forest cover is a necessary component of panda habitat, and deforestation is a good indicator of loss and fragmentation of panda habitat (Liu et al. 2001; Viña et al. 2007). However, forest cover is only one of many components of panda habitat and many forests, including recently logged forests, plantations, and forests without understory bamboo, cannot provide suitable habitat to the giant panda (Bearer et al. 2008). Therefore, an evaluation based on an increase in the amount of forested areas, without consideration of forest characteristics, such as presence or absence of understory bamboo, as an indicator of panda habitat improvement may produce biased results. Better measures of panda habitat dynamics and enhanced understanding of their underlying drivers thus rely on habitat evaluation approaches that use vegetation information beyond nominal classification of land cover types. Information on projected future dynamics of species distribution and habitat is also essential for conservation planning, especially under currently rapid environmental changes (Hannah et al. 2002; Heller and Zavaleta 2009). Shifts in species’ distributional ranges due to climate change have been observed (Chen et al. 2011; Parmesan 2006; Root et al. 2003) and projected (Araújo et al. 2006; Bakkenes et al. 2002; La Sorte and Jetz 2010) for many species in diverse taxa, posing challenges to biodiversity conservation (Coetzee et al. 2009; Hannah et al. 7 2002; Heller and Zavaleta 2009). However, although repeated calls have been made for climate change impact assessments for China’s biodiversity, including the giant panda (Liu and Raven 2010; Ministry of Environmental Protection of China 2010), assessments on the potential impacts of climate change on the giant panda are still missing, both in the literature and in conservation planning. While current conservation efforts, including establishment of protected areas and several incentive-based conservation programs (e.g., the Nature Forest Conservation Program and the Grain-to-Green Program), have effectively alleviated the impacts of land use/cover change on giant panda habitat (Liu et al. 2008; Viña et al. 2011), climate change is believed to become the next major threat to the giant panda’s survival (Swaisgood et al. 2010; Yan et al. 2004). Therefore, quantitative and spatially explicit assessments of the impacts of climate change on panda habitat are urgently needed. Recent advances in remote sensing technology and species distribution modeling may provide a solution for the problem of the above-mentioned information gaps. High temporal resolution remotely sensed data, such as those collected by the Moderate Resolution Imaging Spectroradiometer (MODIS), provide detailed information on the seasonal variability in biophysical characteristics (e.g., biomass) of vegetation. The phenological characteristics of vegetation captured by remotely sensed data, i.e., land surface phenology (Friedl et al. 2006), reflect different land cover, forest, and plant functional types (DeFries et al. 1995; Reed et al. 1994; Sun et al. 2008). Since understory vegetation also affects the phenological characteristics obtained from MODIS imagery (Viña et al. 2008), which has nearly global coverage and temporally continuous availability, land surface phenology has the potential to be used for mapping understory bamboo across large geographic regions. 8 In addition, in the past two decades, species distribution models (SDMs) have developed rapidly and have become a useful tool for assessing the distribution of species and their habitat across space and time (Elith and Leathwick 2009; Guisan and Thuiller 2005; Guisan and Zimmermann 2000). SDMs relate species occurrence or abundance in geographic locations with the environmental properties (e.g., climate, soil types) and spatial characteristics (e.g., connectivity to forested patches, distance to water bodies) of those locations (Guisan and Zimmermann 2000). The established species-environment relationships can then be used to map species distribution across space (or even beyond the geographic range of species’ occurrence and abundance data), given spatially continuous information on environmental conditions. When the relationships are applied to environmental conditions for different time periods, SDMs can also be used to predict temporal dynamics of species distribution in response to the changed environment (Guisan and Thuiller 2005). In particular, SDMs are useful for assessing the impacts of climate change on the distribution of species and their habitat (Araújo et al. 2005; Pearson and Dawson 2003) and also for conservation planning in response to climate change (Araújo et al. 2011; Hannah et al. 2007; Marini et al. 2010). Although diverse statistical approaches have been used in SDMs to deal with different types of species data (e.g., presence-absence, presence-only and abundance data) and to establish species-environment relationships (Guisan and Zimmermann 2000), most SDMs are based on the concept of ecological niche, and thus SDMs are also referred to as ecological niche models (Elith and Leathwick 2009). Hutchinson (1957) defined the ecological niche of a species as an ‘ndimensional hypervolume’ within which species can maintain viable populations. Each dimension of the hypervolume is defined by each of the environmental conditions and resources that allow the species to persist indefinitely (Hutchinson 1957). Hutchinson’s niche concept 9 distinguishes environmental space from geographic space. While a species niche is defined in environmental space, the distribution of the species in geographic space is determined by the geographic distribution of the environmental conditions and resources defining the niche. The process of establishing species-environment relationships in species distribution modeling can be viewed as a process of estimating species niches, and the mapping process is to project the niches from environmental space to geographic space. While SDMs provide the approach to link species distribution between environmental and geographic space, remote sensing provides spatially-explicit environmental data. Remote sensing can measure diverse vegetation characteristics that reflect species distribution or relate to the quality of species habitat. Therefore, variables derived from remotely sensed data portraying vegetation characteristics and other environmental conditions [e.g., topographic variables from the Shuttle Radar Topography Mission (SRTM; Farr et al. 2007) and precipitation data from the Tropical Rainfall Measuring Mission (TRMM; Kummerow et al. 1998)] have been used in SDMs to map the distribution of species and their habitat (Mason et al. 2003; Morisette et al. 2006; Saatchi et al. 2008; Zimmermann et al. 2007). In addition, remote sensing has the ability to repeatedly measure land surface characteristics at the same locations, and thus to detect landscape changes that drive the dynamics of species distribution (Turner et al. 2003). Therefore, the integration of remote sensing technology and species distribution modeling can be a useful tool for not only mapping the spatial patterns of species distribution, but also monitoring their dynamics and investigating underlying drivers. This dissertation was developed with the goal to fulfill the above-mentioned information gaps that hinder effective conservation of the giant panda: lack of detailed information on understory bamboo distribution across large geographic regions, lack of detailed information on 10 short-term spatiotemporal dynamics of panda habitat and underlying drivers, and lack of information on the potential impacts of climate change on the long-term dynamics of panda habitat. It is hoped that the findings of this research and the approaches developed in the study will provide essential information and useful tools to better understand the spatial patterns and temporal dynamics of giant panda habitat, and to assist in the design of conservation practices for this endangered species, as well as many other species that occur within the panda’s distributional range, under a rapidly changing environment. It is also hoped that the findings and approaches will contribute to biodiversity conservation beyond the species and the geographic range evaluated here. Research Objectives The overall goals of this dissertation are to develop approaches for investigating spatiotemporal dynamics of giant panda habitat and to understand how panda habitat changes across space and time. Specific objectives are to (1) develop a remote sensing approach for detecting the distribution of understory bamboo, the giant panda’s staple food, across large spatial extents; (2) develop a modeling approach for monitoring panda habitat changes across space and time; (3) evaluate the effects of current conservation efforts on spatiotemporal dynamics of panda habitat; and (4) assess the potential impacts of climate change on long-term dynamics of panda habitat. The following four chapters (i.e., Chapter 2 to Chapter 5) address each of these objectives, respectively, while the final chapter synthesizes the findings of the previous chapters. In Chapter 2, I evaluate the effectiveness of land surface phenology in detecting presence of bamboo species under forest canopy. By deriving metrics portraying the phenological 11 characteristics of vegetated land surface from a time series of MODIS data, I assess the differences in those phenology metrics between forests with and without understory bamboo and develop an approach based on these differences for mapping understory bamboo distribution. In Chapter 3, the effectiveness of phenology metrics, which contain information on both forest and understory bamboo, for monitoring temporal changes in panda habitat is evaluated. I develop species distribution models with different phenology metrics, examine how those metrics affect the predictive power and temporal transferability of the models, and identify the best metric set for monitoring panda habitat changes. In Chapter 4, the best model generated in Chapter 3 is applied to detect panda habitat change in a nature reserve designated specifically for panda conservation. Using spatial regression models, I investigate the effects of the implementation of an incentive-based conservation program, as well as several biophysical and anthropogenic factors, on the spatiotemporal dynamics of panda habitat. In Chapter 5, the potential impacts of climate change on the spatial distribution of understory bamboo and giant panda habitat are investigated. By generating an ensemble of projections of future panda habitat with different bamboo dispersal scenarios under projected climates from different global climate models and greenhouse gas emission scenarios, I evaluate the temporal changes in the amount and spatial distribution of panda habitat across the 21st century. Finally in Chapter 6, I summarize the findings of previous chapters and discuss their implications for giant panda and global biodiversity conservation. Study Areas In Chapters 2, 3 and 4, I focus on Wolong Nature Reserve, Sichuan Province, China (Figure 1.1) mainly for four reasons. First, this reserve is regarded as a flagship nature reserve. 12 What we learn from Wolong Nature Reserve may not only facilitate the conservation planning in this reserve, but also provide guidance for the design and management of other reserves in China and even around the world. Second, this reserve is located within a global biodiversity hotspot (Myers et al. 2000). Efforts for conserving panda habitat in this reserve may also protect the habitats of thousands of other species living within this region. Third, field data (including the occurrence of bamboo species and the giant panda) collected during the third national giant panda survey, which covered almost all known panda habitat in the reserve, provide good on-theground information on bamboo occurrence. The data allow the development and evaluation of the remote sensing approach for mapping understory bamboo (Chapter 2). Due to the long-term study of human-panda interactions in this reserve conducted by our research team, panda occurrence data collected during different time periods are also available for developing and evaluating the modeling approach for detecting panda habitat changes (Chapter 3). Finally, knowledge on the long-term human-panda interactions and the effects of conservation policies on them provides background information for investigating the factors affecting the spatiotemporal dynamics of panda habitat (Chapter 4). In Chapter 5, I focus on the Qinling Mountains region, Shaanxi Province, China (Figure 1.1) because this region sustains a panda population which has important contributions to the genetic diversity of remaining wild pandas (Lü et al. 2001; Wan et al. 2005), but lives in a habitat patch completely isolated from other patches. Therefore, while its genetic peculiarity makes this population particularly valuable for conservation, geographic isolation makes it vulnerable to environmental changes. In addition, geographic isolation also makes it possible to assess the potential impacts of climate change on panda habitat across the entire mountain region. 13 Figure 1.1. Locations and topography of the Wolong Nature Reserve, Sichuan Province and the  Qinling Mountains region, Shaanxi Province, China.  For interpretation of the references to  color in this and all other figures, the reader is referred to the electronic version of this  dissertation.    14 Wolong Nature Reserve Wolong Nature Reserve (Figure 1.1), which lies between the Sichuan basin and the Tibetan highlands, is characterized by a wide vertical variation in topography, climate and soils, together with a diverse flora and fauna. Established in 1963 and expanded to its current extent in 1975, 2 this reserve is one of the largest nature reserves (ca. 2,000 km ) specifically designated for giant panda conservation. It protects more than 4,000 plant species and 2,200 animal species (Schaller et al. 1985), including approximately 10% of the entire wild panda population (State Forestry Administration 2006). Natural vegetation in the reserve varies along the elevation gradient (Schaller et al. 1985). Broadleaf forests are dominated by evergreen species below 1,600 m, and by a mixture of evergreen and deciduous species between 1,600 and 2,000 m. Above 2,000 m, a mixed coniferous and deciduous broadleaf forest gradually changes to a subalpine coniferous forest around 2,600 m. The forest reaches about 3,600 m where it is replaced by alpine meadows. Under forest canopies, evergreen bamboo species dominate the understory layer. While seven native bamboo species are found in the reserve, two of them, arrow bamboo (Bashania faberi) and umbrella bamboo (Fargesia robusta), are dominant and constitute the major food for giant pandas (Schaller et al. 1985). While arrow bamboo is mainly distributed between 2,500 and 3,400 m in elevation, umbrella bamboo usually occurs between 1,600 and 2,650 m. The giant panda is mainly distributed in the mixed and subalpine coniferous forests between 2,250 and 2,750 m (Liu et al. 1999). Wolong Nature Reserve entirely comprises Wolong and Gengda townships and part of Sanjiang township (Figure 1.1). All the residents of Wolong and Gengda townships (ca. 4,900 residents in ca. 1,200 households) live inside the reserve, whereas all the residents of Sanjiang (ca. 3,250 residents in ca. 900 households) live outside the reserve. Agriculture is the major 15 economic activity and fuelwood is the major energy source among the households of the three townships. Diverse human activities, including illegal timber harvesting, fuelwood collection, cultivation, livestock husbandry and infrastructure construction, are the major threats to the giant panda habitat in the reserve (Liu et al. 2001). Qinling Mountains The Qinling Mountains lie in an east-west direction across the south of Shaanxi Province, China (Figure 1.1). Because it is a divide of two major watersheds, the Yellow River and the Yangtze River, this mountain region forms a natural boundary between northern and southern China and also constitutes a climatic transition, from cold and dry in the north to warm and wet in the south (Pan et al. 2001). Due to the north-south climatic transition and its vertical zonation along elevational gradients (up to ca. 3,700 m), the Qinling Mountains exhibit high biodiversity, supporting more than 3,000 plant species, over 300 bird species, and more than 85 mammal species, including the endangered giant panda (Pan et al. 1988). This mountain region is home to ca. 275 wild giant panda individuals (about 17% of the entire wild population) and exhibits the highest population density of all the remaining panda habitat areas (State Forestry Administration 2006). Nevertheless, the panda population in this region is isolated from other major populations due to a long history of human habitation (Loucks et al. 2003). The distribution of vegetation in the Qinling Mountains follows an elevation gradient. Natural vegetation below 1,350 m is broadleaf deciduous forests, but most of them have been replaced by croplands and built-up areas (Pan et al. 2001). Mixed broadleaf/coniferous forests are mainly distributed between 1,350 and 2,400 m, coniferous forests located between 2,400 and 16 3,350 m, and alpine meadows dominate the areas above 3,350 m (Pan et al. 2001). Three understory bamboo species, Qinling arrow (Fargesia qinlingensis), dragon-head (F. dracocephala), and wooden (Bashania fargesii), account for more than 90% of bamboo cover under the forests in this region and constitute the main diet of the panda population (State Forestry Administration 2006). While Qinling arrow and wooden bamboos are distributed from 1,800 – 3,000 m and from 900 – 1,900 m, respectively, the distribution of dragon-head bamboo is patchier and overlaps with the other two species at elevations between 1,100 and 2,300 m (Pan et al. 2001; State Forestry Administration 2006). The giant pandas in this region forage for Qinling arrow bamboo at higher elevations from June to September and forage for wooden bamboo at lower elevations from October to May (Pan et al. 1988). 17 CHAPTER 2 MAPPING UNDERSTORY BAMBOO USING PHENOLOGICAL CHARACTERISTICS DERIVED FROM REMOTELY SENSED DATA In collaboration with Andrés Viña, Scott Bearer, Weihua Xu, Zhiyun Ouyang, Hemin Zhang, and Jianguo Liu This chapter is based on: Tuanmu, M.-N., Viña, A., Bearer, S., Xu, W., Ouyang, Z., Zhang, H., & Liu, J. (2010). Mapping understory vegetation using phenological characteristics derived from remotely sensed data. Remote Sensing of Environment, 114, 1833-1844. 18 Abstract Understory vegetation is an important component in forest ecosystems not only because of its contributions to forest structure, function and species composition, but also due to its essential role in supporting wildlife species and ecosystem services. Therefore, understanding the spatiotemporal dynamics of understory vegetation is essential for management and conservation. Nevertheless, detailed information on the distribution of understory vegetation across large spatial extents is usually unavailable, due to the interference of overstory canopy on the remote detection of understory vegetation. While many efforts have been made to overcome this challenge, mapping understory vegetation across large spatial extents is still limited due to a lack of generality of the developed methods and limited availability of required remotely sensed data. In this study, we used understory bamboo in Wolong Nature Reserve, China as a case study to develop and test an effective and practical remote sensing approach for mapping understory vegetation. Using phenology metrics generated from a time series of Moderate Resolution Imaging Spectroradiometer data, we characterized the phenological features of forests with understory bamboo. Using maximum entropy modeling together with these phenology metrics, we successfully mapped the spatial distribution of understory bamboo (kappa: 0.59; AUC: 0.85). In addition, by incorporating elevation information we further mapped the distribution of two individual bamboo species, Bashania faberi and Fargesia robusta (kappa: 0.68 and 0.70; AUC: 0.91 and 0.92, respectively). Due to its generality, flexibility and extensibility, this approach constitutes an improvement to the remote detection of understory vegetation, making it suitable for mapping different understory species in different geographic settings. Both biodiversity conservation and wildlife habitat management may benefit from the 19 detailed information on understory vegetation across large areas through the applications of this approach. Introduction Understory vegetation plays an important role in forest ecosystems (Gilliam 2007; Nilsson and Wardle 2005). Different structure and species compositions of understory plants, including native and exotic species, can affect regeneration of tree species, alter forest succession, and change species diversity through physical, chemical, and biological mechanisms (Royo and Carson 2006; Urgenson et al. 2009). Understory vegetation also provides essential shelter and food resources for wildlife, and thus its structure and composition is usually associated with the diversity and abundance of many wildlife species (Díaz et al. 2005; Hagar 2007). Besides its ecological functions in forest ecosystems, non-timber forest products provided by many understory plants (e.g., fibers from bamboo and rattans and medicines from medicinal herbs) are economically important for many countries (Iqbal 1993). Therefore, understanding the spatiotemporal dynamics of understory vegetation is essential not only for wildlife and biodiversity conservation (Deal 2007; Estades and Temple 1999), but also for sustainable forest management (FAO 1995). While the importance of understory vegetation is well known, detailed information on its spatial distribution across large areas and its dynamics at fine temporal resolutions is usually unavailable because conventional methods for gathering vegetation information emphasize ground-based surveys, which are time-consuming, labor-intensive, and sometimes logistically unfeasible. Although remote sensing is a useful alternative tool for gathering vegetation information across large areas and over time (Jensen 2007; Roughgarden et al. 1991), remote 20 detection of understory plants is limited due to the interference of overstory canopies. While the presence of understory vegetation influences the signals received by remote sensors, its relative contributions to the signals vary with the structure and species composition of both over- and under-story vegetation (Eriksson et al. 2006; Rautiainen et al. 2007). Because of the complex and nonlinear interactions between the reflectance of over- and under-story components, distinguishing the signals of understory vegetation from overstory canopies and characterizing understory vegetation are challenging. Many efforts have been made to overcome this challenge via advanced classification algorithms. For instance, using artificial neural networks to capture the non-linear relationship between the reflectance of over- and under-story vegetation, Linderman et al. (2004) have successfully detected understory bamboo distribution with an overall accuracy of 80% (and a kappa statistic of 0.56, which we calculated from the confusion matrix reported in Table 3 of Linderman et al., 2004). Wang et al. (2009a) classified an Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) image into three understory cover classes with a kappa statistic of 0.60 by integrating a neural network and a Geographic Information System (GIS) expert system. However, applying these methods in a larger spatial extent may be limited by the availability of cloud-free images with high spatial resolutions, and locally specific rules in the expert system. Other studies have used active sensor systems. For example, data acquired using Light Detection and Range (LiDAR) sensors have been used to characterize the three-dimensional structure of forests (Lefsky et al. 2002; Vierling et al. 2008), and map understory plants in boreal forests (Korpela 2008; Peckham et al. 2009). Combined with hyperspectral images, LiDAR data have also been used to map an understory invasive species in tropical forests (Asner et al. 2008). 21 However, the most important limitation for ecological applications of these airborne sensor data are their high acquisition costs and low availability, especially in developing countries (Vierling et al. 2008). Differences in phenology between over- and under-story vegetation have also been used to detect understory vegetation. For instance, due to the earlier senescence of tree leaves as compared to the leaves of an understory invasive shrub species, Resasco et al. (2007) found that stands with high coverage of the understory shrub could be separated from those with low/zero coverage using late-fall Landsat imagery. In addition, Chastain, Jr. and Townsend (2007) and Wang et al. (2009b) also successfully used leaf-off Landsat images to detect evergreen understory vegetation under deciduous forests (with kappa statistics of 0.755-0.806 and 0.59, respectively). However, the limited temporal windows when the phenological difference between over- and under-story can be detected further reduce the data availability. Furthermore, the inter-annual variability of vegetation phenology due to variable climatic conditions may change the optimal dates for separating over- and under-story components in different years (Resasco et al. 2007). High temporal resolution remotely sensed data, such as those collected by the Moderate Resolution Imaging Spectroradiometer (MODIS), may provide a solution for the problem of mapping understory vegetation across large areas. With high frequency of acquisitions, those data reduce the problem of cloud contamination and provide detailed information on the temporal dynamics of the land surface. As vegetation phenology causes changes in surface reflectance over time, they can be captured by multi-temporal remotely sensed data (Ahl et al. 2006; Reed et al. 1994; Schwartz and Reed 1999). The phenological patterns captured by remotely sensed data are termed land surface phenology (Friedl et al. 2006), in order to 22 distinguish them from the traditional definition of vegetation phenology. Previous studies have shown that these phenological characteristics are useful for classifying different land cover types (DeFries et al. 1995; Hansen et al. 2000; Lloyd 1990), monitoring land cover change (de Beurs and Henebry 2004; Lenney et al. 1996), and differentiating forest classes and types (Townsend and Walsh 2001). However, the phenological characteristics captured by multi-temporal remotely sensed data are affected by not only the dominant overstory vegetation, but also the understory vegetation. Viña et al. (2008) have shown that the seasonal patterns of vegetation indices derived from MODIS data were different between forests with similar overstory vegetation but with different understory coverage. This conspicuous difference in phenological patterns suggests that land surface phenology has the potential to be used for mapping understory vegetation across large areas. The main goal of this study was to develop an approach for deriving detailed information on the spatial distribution of understory vegetation using the phenological patterns detected by multi-temporal remotely sensed data. We selected the bamboo species living under the canopy of temperate forests in Wolong Nature Reserve, China as a case study. Because bamboo species dominate the understory vegetation within this region and provide essential food for several wildlife species (Schaller et al. 1985) including the endangered giant panda (Ailuropoda melanoleuca), identifying their spatial distribution has direct wildlife conservation implications. To develop and test our approach, we (1) generated phenology metrics from a time series of MODIS data and examined the phenological characteristics of forests with understory bamboo; (2) developed a spatial model for mapping understory bamboo distribution using field data, phenology metrics, and species distribution modeling; and (3) explored the potential application of this approach to mapping and differentiating individual bamboo species. 23 Materials and Methods Field data We obtained bamboo presence data from the Third National Giant Panda Survey (State Forestry Administration 2006). Survey teams collected field data in Wolong Nature Reserve following transects across different vegetation types during the summer of 2001. Along the transects, surveyors recorded, using Global Positioning System receivers, the geographic locations of 468 sites with signs of giant panda activity, including fecal droppings, feeding sites, dens, footprints, and visual sightings. In each site, surveyors identified bamboo species and assigned the coverage of understory bamboo to one of four categories (0-24%, 25-49%, 50-74% and 75-100%). Because low bamboo cover has a limited contribution to the surface reflectance registered by satellite sensors and may not provide enough food or shelter for wildlife species, we only used the locations where bamboo cover was estimated to be 25% or higher (Figure 2.1). In these locations, besides the two dominant bamboo species [i.e., arrow bamboo (Bashania faberi) and umbrella bamboo (Fargesia robusta)], F. nitida and Yushania brevipaniculata were also found. Although bamboo species may occur in shrubland, in the study area shrubland with bamboo are only found in small patches within clear-cut areas for past timber production (Reid et al. 1991; Schaller et al. 1985). In addition, only one of the 468 field sites has a vegetation type of shrubland with bamboo species in the dataset. Therefore, we excluded that site from the following analyses and focused on the bamboo under forests in this study. 24 Figure 2.1. Locations of the field plots where arrow or umbrella bamboo cover was 25% or  higher according to the Third National Giant Panda Survey conducted in Wolong Nature  Reserve in 2001.  Remotely sensed data and phenology metrics We obtained a time series of MODIS surface reflectance imagery (8-day L3 Global 250m product, MOD09Q1) acquired between May 2000 and April 2004 through the Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov/). Using the surface reflectance values in the red (RRED, 620-670 nm) and near infrared (RNIR, 841-876 nm) spectral bands, we calculated the Wide Dynamic Range Vegetation Index (WDRVI; Gitelson 2004): 25 WDRVI = (α·RNIR – RRED) / (α·RNIR + RRED) (2.1) where α is a weighting coefficient set as 0.25 following Henebry et al. (2004). This coefficient reduces the saturation problem of the Normalized Difference Vegetation Index (NDVI) under moderate-to-high biomass conditions (Gitelson 2004). The WDRVI has been proved to be linearly related with leaf area index (LAI) and sensitive to changes in LAI up to 6 (Gitelson 2004; Viña et al. 2004b). Therefore, the WDRVI appears to be more suitable than the widely-used NDVI for studying the changes in green biomass under high LAI values, such as the forests with dense understory bamboo in our study area. In order to further reduce the cloud contamination, which lowered WDRVI values, we generated a time series of 16-day composites using the maximum value between two consecutive 8-day periods. Using TIMESAT 2.3 (Jönsson and Eklundh 2004, 2006), we smoothed the time series (May 2000 – April 2004) of 16-day WDRVI composites for each pixel by means of the adaptive Savitzky-Golay filter. With these data we obtained three full phenological cycles (2001-2003) and calculated 11 phenology metrics for each cycle: (1) the base level, corresponding to the average between the starting and ending minimum values of each cycle (A in Figure 2.2); (2) the maximum level, corresponding to the highest value in each cycle (B in Figure 2.2); (3) the amplitude, calculated as the difference between the maximum and the base levels (C in Figure 2.2); (4) the date of the start of the season (SOS), determined as the date when WDRVI values increase to 20% of the difference between the maximum WDRVI value and the minimum value at the start of each cycle (D in Figure 2.2); (5) the date of the end of the season (EOS), defined as the date when WDRVI values decrease to 20% of the difference between the maximum value and the minimum value at the end of each cycle (F in Figure 2.2); (6) the date of the middle of the season (MOS), determined as the median of the two dates when WDRVI values increase 26 (decrease) to 80% from the minimum value at the start (end) of each cycle (E in Figure 2.2); (7) the length of the season, defined as the difference between SOS and EOS (G in Figure 2.2); (8) the large integral, obtained as the area under the smoothed curve between SOS and EOS (H in Figure 2.2); (9) the small integral, defined as the large integral minus the area below the base level (I in Figure 2.2); (10) the increase and (11) decrease rates, calculated as the slopes of two lines across the 20% and 80% level points on the left and right sides of the MOS, respectively (J and K in Figure 2.2). Figure 2.2. Phenology metrics (A – K) derived from the smoothed values (triangles) of a time  series of WDRVI values.  A ‐ Base level; B ‐ Maximum level; C ‐ Amplitude; D ‐ Date of the start  of the season; E ‐ Date of the middle of the season; F ‐ Date of the end of the season; G ‐ Length  of the season; H ‐ Large integral; I ‐ Small integral; J ‐ Increase rate; K ‐ Decrease rate.  27 In some pixels, phenology metrics could not be obtained on a particular cycle due to either lack of detectable phenological cycles or incomplete cycles within a year. To reduce both the effects of missing cycles and inter-annual variability in phenology metrics, we calculated the average of the annual values from 2001 to 2003 for each pixel. We treated the pixels which had less than two valid annual values as missing data and excluded them from the following analyses. Phenological characteristics of forests with understory bamboo To test whether the phenological characteristics of the forests with understory bamboo can be distinguished from other land cover types using the 11 phenology metrics, we compared them among five groups of pixels. For the first group (All), we randomly selected 1,000 pixels from the study area as a representative of the entire area. The second group (For) was a random selection of 1,000 pixels with a forest cover. Pixels with a forest cover were determined based on a binary forest cover map derived from a Landsat-5 Thematic Mapper image acquired on 13 June 2001 (Viña et al. 2007) and re-sampled to 250 × 250 m/pixel using the majority algorithm. A series of tests on the selected pixels’ representativeness of the entire study area and forest area indicated that the variation of the means of pixel values (i.e., values of each phenology metric) decreased with the increase in the number of selected pixels, but the change became negligible when more than 500 pixels were selected. On the other hand, selecting more pixels may increase the spatial autocorrelation among selected pixels. To achieve representativeness and reduce the effects of spatial autocorrelation on statistical tests (see below), we used 1,000 randomly selected pixels in this analysis. The third group (Bam) contained 356 pixels where bamboo cover was 25% or higher (including all four bamboo species, i.e., B. faberi, F. robusta, F. nitida and Y. brevipaniculata) according to the field data. The remaining two groups were composed of 145 28 pixels where arrow bamboo cover was 25% or higher (Arr), and 184 pixels where umbrella bamboo cover was 25% or higher (Umb), respectively. We then compared the pixel values between each pair of the five groups using the Mann-Whitney U-test for each of the 11 phenology metrics. We used this non-parametric test because the pixel values were not normally distributed, according to a Shapiro-Wilk normality test performed. The conclusions on significant differences drawn from this test are more conservative since the non-parametric Utest is less powerful to detect significant differences between groups, as compared to parametric methods (e.g., t-test; Sheskin 2000). For this analysis, we conducted presence vs. background comparisons (e.g., forest pixels with bamboo vs. random selection of all forest pixels), rather than presence vs. absence comparisons (e.g., forest pixels with bamboo vs. forest pixels without bamboo) mainly for two reasons. First, although possible, it is not practical to confirm absolute absence of understory species across a 6.25 ha field plot (i.e., the area of a 250 × 250m MODIS pixel). And second, the modeling algorithm used to map understory bamboo is based on the difference in the values of predictor variables obtained in presence locations and in background (random) locations (see below). Therefore, the presence vs. background comparisons examined the information content of the data directly used by the modeling algorithm to map understory vegetation. Overall bamboo distribution In order to distinguish the pixels with understory bamboo from the others based on their phenological characteristics, we used the maximum entropy modeling framework (Maxent). Maxent is an algorithm designed to make predictions based on incomplete information (Phillips et al. 2006) and has been proven to be one of the best methods for mapping species distribution 29 (Elith et al. 2006). The algorithm contrasts the environmental conditions (characterized in a multi-dimensional space defined by environmental variables) in species presence locations vs. the conditions in background locations (i.e., the entire study area). It then establishes the species-environment relationship by matching the contrasts and approaching a maximum entropy distribution (i.e., maximum uniform distribution) simultaneously (Phillips et al. 2006). The relationship can be used to estimate the probability of species occurrence across the entire study area given the spatial patterns of the environmental variables (Phillips et al. 2006; Phillips and Dudík 2008). Besides its good performance on mapping species distribution (Elith et al. 2006), Maxent also has several characteristics which make it suitable for mapping understory vegetation based on phenological characteristics obtained from MODIS data. First, it uses presence-only, rather than presence/absence data. This is important since it is not logistically feasible to confirm absence of understory bamboo in an entire 250 × 250 m area. Using the presence-only procedure, we can avoid the potential biases caused by uncertain or false absence data. In addition, like other machine learning methods (e.g., neural networks), Maxent can capture complex and nonlinear species-environment relationships, even with noise in input data (Elith et al. 2006; Phillips et al. 2006). Finally its continuous output values, i.e., species presence probabilities, make Maxent a fuzzy classifier that provides more detailed information on understory vegetation distribution than binary outputs (i.e., presence/absence). We used the software Maxent (version 3.3.1; Phillips et al. 2006) to generate a model for mapping overall bamboo distribution. For generating the model, we used the pixels whose bamboo cover was 25% or higher (including all four bamboo species) as presence data, 10,000 pixels randomly selected from the study area as background data following the suggestions of 30 Phillips & Dudík (2008), and the 11 phenology metrics derived from a time series of MODISWDRVI values as predictor variables. We used 70% of the bamboo presence data (249 pixels) as a training dataset and the remaining 30% (107 pixels) as a validation dataset. In order to reduce the potential effects of the random data partitioning, we ran the model 20 times (replicates) with different random data partitions for each run and averaged the predicted bamboo presence probabilities over the 20 runs for each pixel. This number of replicates was used because pilot tests showed the variation of model outputs, measured as the mean of standard deviations among different runs, decreased with the increase in the number of runs, but changed negligibly after 20 runs. Individual bamboo distribution Besides mapping overall bamboo distribution, we also explored the ability of our approach to map individual bamboo species. With the same method described earlier, we generated a model and obtained the average of presence probabilities over 20 model runs for each of the two dominant bamboo species, i.e., arrow and umbrella bamboo. Because the two bamboo species are distributed within different elevation ranges, we hypothesized that adding elevation information into the model would improve the model’s ability to separate the two species. To test this hypothesis, besides the model generated with 11 phenology metrics, we also generated a model using elevation as an additional predictor variable and compared the model performance (see below) between the two models. We obtained the information on elevation from a digital elevation model created by the Shuttle Radar Topography Mission and re-sampled the original data to 250 × 250 m/pixel using the nearest neighbor algorithm to keep consistent spatial resolution with the other data. 31 Model validation and comparison We used both threshold-dependent and threshold-independent methods to validate models and evaluate their performance. The Cohen’s kappa analysis, a chance-corrected measure of agreement (Cohen 1960), was selected for the threshold-dependent method because it is commonly used for evaluating classification accuracy of remote sensing imagery and also used in previous studies on mapping understory vegetation (Chastain Jr. and Townsend 2007; Wang et al. 2009a; Wang et al. 2009b). The kappa value ranges between 0 and 1 with a larger value indicating better model performance (Cohen 1960). Model performance can be judged as excellent if kappa > 0.75, good if 0.75 > kappa > 0.4, or poor if kappa < 0.4 (Araújo et al. 2005; Landis and Koch 1977). Because only presence data on bamboo distribution were available, we performed the analysis by contrasting presence pixels to those randomly selected from the study area (background pixels). To avoid the potential failure of kappa analysis with unbalanced validation data (McPherson et al. 2004), we randomly selected 100 background pixels to make the number be close to the number of presence pixels in validation datasets (i.e., 30% of bamboo presence data or 107 pixels). The threshold for cutting off continuous outputs from each model run was determined by the kappa maximization approach, which finds the threshold corresponding to the maximum kappa value (Liu et al. 2005a). Because 100 pixels (ca. 0.3% of total pixels) were not representative of the entire area, we calculated 30 kappa values for each model run by using the same presence pixels but different 100 random background pixels, and then obtained an average of the 30 values. The receiver operating characteristic (ROC) analysis, a threshold-independent method, is also a widely-used method for evaluating the accuracy of classification models (Fielding and 32 Bell 1997; Pearce and Ferrier 2000). The ROC curve is generated by plotting sensitivity values (i.e., fraction of true positive) against 1-specificity values (i.e., fraction of false positive) for every possible threshold (Hanley and Mcneil 1982). The area under the ROC curve (AUC) provides a single-value measurement of model performance. While omission errors reduce sensitivity and commission errors reduce specificity, both types of errors equally reduce the AUC value. While an AUC value of 1 indicates a perfect model, a value of 0.5 indicates a random model. A standard for judging model performance based on AUC values (Araújo et al. 2005; Swets 1988) is: excellent (AUC > 0.9), good (0.9 > AUC > 0.8), fair (0.8 > AUC > 0.7), poor (0.7 > AUC > 0.6), and failed (0.6 > AUC > 0.5). Similar to the calculation of the kappa statistic, we also used presence/background validation data for calculating AUC values. However, because AUC values are not sensitive to unbalanced validation data (McPherson et al. 2004; Zweig and Campbell 1993) and no statistical test is involved in the calculation of AUC values, the number of background pixels does not affect the calculated values if the background pixels are representative of the entire study area. Therefore, we used the default background pixels in the Maxent software (i.e., 10,000 background pixels) to calculate AUC values. By using presence/background data, the kappa and AUC values calculated in this study tend to be underestimated because some of the background pixels are actually presence pixels, which artificially increase commission errors. In addition, the degree of underestimation is determined by the proportion of actual presence pixels in the background pixels, which is, in turn, determined by the actual proportion of habitat in the entire study area. Since the actual proportion of habitat is almost always unknown, direct comparisons of the values calculated in this study with the values reported in other studies should be done with caution. Comparisons between different models using these statistics are valid only if the models are generated for the 33 same species in the same study area and are evaluated using the same presence/background data, as was done in this study (see below). In order to examine the relative importance of each phenology metric for mapping bamboo distribution, we conducted a jackknife analysis on model performance by using the Maxent software. For this, the software calculated the AUC values of models containing only one of the 11 metrics and of models containing all, but one of the metrics used as predictor variables, through a jackknife re-sampling approach. In this analysis, a higher AUC value for a model containing only one metric indicates that the specific metric is more informative for mapping bamboo distribution than other metrics. In contrast, a lower AUC value for a model without one specific metric indicates that the metric contains more information for mapping bamboo not provided by the other metrics. For comparing the performance of the individual bamboo models with and without elevation information, besides calculating kappa and AUC values, we also used the minimum predicted area (MPA) method (Engler et al. 2004). The MPA is the minimum area which is constituted by all pixels whose species presence probabilities are above a defined threshold and encompasses a specified percentage (e.g. 95%) of presence locations (Engler et al. 2004). With presence-only validation data, a model predicting species present everywhere has the best performance because it correctly predicts all presence locations, but the model is useless. Therefore, the MPA method evaluates model performance based on the parsimonious concept that a good model should predict the habitat as small as possible (i.e., with low commission errors), but it still encompasses a maximum number of presence locations (i.e., with low omission errors). In this analysis, a threshold was selected for each model run so that the pixels with probabilities above the threshold encompassed 95% of presence locations. Kappa and AUC 34 values and the proportions of MPA to the whole study area obtained for 20 runs of the models with and without elevation were compared using Mann-Whitney U-tests. Results Phenological characteristics of forests with understory bamboo Because of a lack of detectable seasonal cycles or incomplete annual cycles in at least two years, phenology metrics could not be calculated in 289 pixels (i.e., 0.76% of the total number of pixels) of the entire study area. Among these pixels, only 112 pixels (i.e., 0.79% of the total number of forest pixels) were covered by forests, according to the forest cover map, and none of them contained field plots with a 25% or higher bamboo cover. Therefore, the impact of their exclusion from the analysis is negligible due to the small proportions of these pixels in the five pixel groups. Significant differences in the values of 11 phenology metrics were found among the five groups of pixels (Figure 2.3), even with the more conservative significant tests conducted using the non-parametric method. Compared to the pixels randomly selected from the whole study area (All), the pixels with understory bamboo (Bam) had significantly (p < 0.001) higher base and maximum levels, a higher amplitude, earlier SOS and MOS, a longer length of season, and higher large and small integrals (Figure 2.3). Compared to pixels with forest cover (For), the pixels of the Bam group had a significantly higher maximum level, a higher amplitude, earlier SOS, MOS, and EOS, higher large and small integrals, and a higher increase rate (Figure 2.3). While the values of the pixels with arrow bamboo (Arr) were significantly different in eight and five of the 11 metrics from the All and For groups, respectively, they were different from the Bam group only in the maximum level and SOS (Figure 2.3). The values of pixels with umbrella 35 bamboo (Umb) were significantly different from the All and For groups in all and eight of the metrics, respectively (Figure 2.3). However, significant differences were found only in the maximum level and MOS between the Umb and Bam pixels (Figure 2.3). Overall bamboo distribution The kappa and AUC values (mean ± 2 SEM) of the 20 runs of the overall bamboo distribution models were 0.591 ± 0.018 and 0.851 ± 0.005, respectively. Both values, even though underestimated due to the use of presence/background validation data, indicated a good model according to the judgment standards (Araújo et al. 2005; Landis and Koch 1977; Swets 1988). The estimated bamboo presence probabilities for pixels ranged between 0 and 0.9 across the study area, with higher probability values occurring at low- and mid-elevations (Figure 2.4a). The highest standard deviation of the estimated presence probabilities was 0.3, but most pixels had standard deviations lower than 0.05 (Figure 2.4b). The low standard deviations indicated that the estimated probabilities did not change much with different data partitioning for training and validation datasets. Among the pixels whose phenology metrics could not be calculated, 137 were located above 3600 m, thus beyond the distribution range of forests and bamboo species in the study area, and 40 of the other pixels were not covered by forests, according to the forest cover map. Therefore, the missing data may affect the estimated bamboo presence probabilities in only 112 pixels (i.e., 0.29% of the total number of pixels). 36 d c d c b c a All For Bam Arr Umb a a b b c 100 Day of Year 200 300 200 b 100 a 0 Day of Year All For Bam Arr Umb d All For Bam Arr Umb WDRVI 0.5 1.0 1.5 a c c ab b a 0.0 a All For Bam Arr Umb All For Bam Arr Umb Day of Year 200 300 400 a WDRVI -0.5 0.0 0.5 1.0 0.5 a -1.0 -0.5 0.0 WDRVI b b a bc ab c All For Bam Arr Umb Figure 2.3. Box plots of 11 phenology metrics calculated using 1,000 pixels randomly selected from the entire study area (All), 1,000  pixels randomly selected from forested areas (For), 356 pixels where field plots with bamboo cover as 25% or higher were located  (Bam), 145 pixels where field plots with arrow bamboo cover as 25% or higher were located (Arr), and 184 pixels where field plots  with umbrella bamboo cover as 25% or higher were located (Umb).  The dark line inside a box indicates the median; the bottom and  the top of a box show the 25th and 75th percentiles, respectively; the ends of the two whiskers indicate the lowest and the highest  37 values within 1.5 inter‐quartile ranges from a box; dots beyond whiskers show outliers whose values were outside the range  indicated by the whiskers.  The letters above boxes show the results of pair‐wise comparisons on the phenology metric values  conducted using Mann‐Whitney U‐tests.  Two boxes share the same letter if there was no significant difference (p > 0.001) between  a a 100 c ab ab a b a ab a All For Bam Arr Umb b ab ab ab a 0 WDRVI/Day 0.005 0.01 0.015 bc 0 WDRVI/Day 0.005 0.01 0.015 All For Bam Arr Umb c All For Bam Arr Umb WDRVI × Day 100 200 a d c ab b a 0 a WDRVI × Day 100 300 500 b Day 250 400 them. All For Bam Arr Umb Figure 2.3. (Continued)  38 All For Bam Arr Umb Figure 2.4. Overall bamboo distribution across Wolong Nature Reserve.  The pixel‐wise (a) mean  values and (b) standard deviations of bamboo presence probabilities were calculated over 20  runs of the overall bamboo distribution model.  The model was generated using 356 pixels with  25% or higher bamboo cover as presence locations and 11 phenology metrics as predictor  variables.  The 289 pixels where phenology metrics could not be determined in at least two  years between 2001 and 2003 from a smoothed curve of a time series of WDRVI values by  TIMESAT are represented in white.  The results of the jackknife analysis on the relative importance of phenology metrics for mapping understory bamboo are shown in Figure 2.5. The models with only the maximum level, base level, or large integral had the highest AUC values (Figure 2.5a), which indicated that those metrics contained the most useful information for mapping understory bamboo. The models 39 without the SOS, small integral, or base level had the lowest AUC values (i.e., largest loss in AUC) (Figure 2.5b), and thus those metrics contained unique information for mapping bamboo distribution. Individual bamboo distribution Although the mean kappa and AUC values indicated that the performance of the arrow and umbrella bamboo models without elevation was fair to good and good to excellent, respectively (Table 2.1), the models could not effectively differentiate the distribution of the two species. While umbrella bamboo had higher presence probabilities at relatively lower elevations (Figure 2.6b), consistent with the field data and our knowledge about the general distribution pattern of this species, the estimated probabilities for arrow bamboo were also high in some low-elevation areas (Figure 2.6a). The relatively low kappa and AUC values obtained for arrow bamboo (Table 2.1) seem to reflect this overestimation at lower elevations. With the incorporation of elevation information, the mean kappa and AUC values increased and the proportion of MPA decreased significantly for both species (Table 2.1). Higher kappa and AUC values and smaller MPA suggested that model performance on mapping individual bamboo species was significantly improved with the addition of elevation as a predictor variable. In addition, the spatial patterns of estimated presence probabilities also showed that the distribution patterns of the two bamboo species can be differentiated (Figure 2.6c and d) with the incorporation of elevation information into the models. 40 Figure 2.5. Mean AUC values of the overall bamboo distribution models (a) with only one of the  11 phenology metrics and (b) with all, but one, metrics as predictor variables.  A higher AUC  value for a model with only one metric indicates that the metric contained more useful  information for mapping bamboo distribution in the full model.  A lower AUC value (larger loss  of the AUC value) for a model without one metric indicates that the metric contained more  information which cannot be represented by the other metrics for mapping bamboo  41 distribution.  Max ‐ Maximum level; Amp ‐ Amplitude; SOS ‐ Date of the start of the season;  MOS ‐ Date of the middle of the season; EOS ‐ Date of the end of the season; L_Int ‐ Large  integral; S_Int ‐ Small integral; I_Rate ‐ Increase rate; D_Rate ‐ Decrease rate.    Table 2.1. Comparisons of the performance of models developed for individual bamboo species  with and without elevation information, based on the results of 20 model runs.          *   Kappa    p‐value of Mann‐ Whitney Test    *    AUC  p‐value of Mann‐ Whitney Test     *  Threshold for MPA Proportion of MPA  *  Without  Elevation  With  Elevation  0.461 ± 0.017  0.681 ± 0.015  ‐10 < 10 0.798 ± 0.009            0.906 ± 0.004    ‐7    < 10 Umbrella Bamboo  Without  Elevation  With  Elevation  0.658 ± 0.022  0.703 ± 0.019  ‐4 < 10   0.900 ± 0.005  0.920 ± 0.005  ‐4  < 10   0.211 ± 0.015  0.229 ± 0.009    0.192 ± 0.013  0.199 ± 0.013    0.476 ± 0.009  0.223 ± 0.005    0.297 ± 0.008  0.212 ± 0.005  to Study Area  p‐value of Mann‐ Whitney Test  Arrow Bamboo    ‐10    < 10 * Values are shown as mean ± 2 SEM.  42 ‐5  < 10 Figure 2.6. Spatial distribution of arrow and umbrella bamboo across Wolong Nature Reserve.   Mean presence probabilities of pixels were calculated over 20 runs of the individual bamboo  distribution models containing only 11 phenology metrics (a and b) and the models containing  the 11 phenology metrics plus elevation (c and d).  Pixels where phenology metrics could not be  determined in at least two years between 2001 and 2003 from a smoothed curve of a time  series of WDRVI values by TIMESAT are represented in white.  43 Discussion In this study, we developed an effective approach for mapping understory vegetation across large spatial extents using remotely sensed data. By taking the advantage of MODIS data’s high temporal resolution, we captured the phenological characteristics, in terms of WDRVI values, of forests with understory bamboo using phenology metrics. We then established the relationship between bamboo presence and phenological characteristics and used it to map understory bamboo by using Maxent. With this approach, we successfully mapped the spatial distribution of bamboo species under temperate forests in Wolong Nature Reserve, China. While land surface phenology derived from time series of remotely sensed data has been used for land cover classification (DeFries et al. 1995), vegetation change detection (de Beurs and Henebry 2004, 2005), canopy phenology monitoring (Ahl et al. 2006; Viña et al. 2004a; Zhang et al. 2003), invasive plant mapping and monitoring (Huang et al. 2009; Morisette et al. 2006), and wildlife habitat characterization (Viña et al. 2008), in this study their applications have been extended to evaluate the spatial distribution of understory plant species growing under a forest canopy. By analyzing the land surface phenology characterized by phenology metrics, we found that forest pixels with understory bamboo can be distinguished from background and forest pixels in the whole study area. While higher base and maximum levels, an earlier SOS, a longer season length, and higher large and small integrals reflected the difference between forests and other land cover types (Figure 2.3), a still higher maximum level, much earlier SOS and MOS, higher yet integrals, and a higher increase rate showed the contributions of bamboo species to the land surface phenology of forest pixels with understory bamboo (Figure 2.3). The high biomass and -1 annual net primary productivity of understory bamboo (dry weights: 5-12 ton·ha and 1.2-1.9 44 -1 -1 -1 -1 -1 ton·ha ·year for arrow bamboo and 15-40 ton·ha and 1.5-3.9 ton·ha ·year for umbrella bamboo) in the study area (Taylor and Qin 1993) may account for the higher maximum level and integrals. The rapid growth in the height of bamboo shoots during the early growing season (Qin et al. 1993; Taylor and Qin 1993) may be the reason for the higher increase rate and earlier SOS of the pixels with bamboo. Since the bamboo species in the study area are evergreen, it was expected that their contribution to the green biomass was larger when overstory tree leaves undergo senescence (e.g., during winter months). However, no difference in the base level between forests with bamboo and background forests was observed. This result could be partially explained by snow cover during the winter months. Although bamboo species under evergreen coniferous forests may not cause difference in land surface phenology as much as those under deciduous or mixed forests, pure evergreen coniferous forests are rare in the study area. While firs (e.g., Abies faxoniana) are dominant in the subalpine coniferous forests above 2,600 m in elevation, birches (e.g., Betula utilis and B. albosinensis) and rhododendrons (e.g., Rhododendron oreodoxa and R. watsonii) are also abundant (Schaller et al. 1985). However, although it is not a major concern in this study, the potential effect of evergreen overstory on the detectability of phenological difference caused by understory vegetation needs further study. In addition, besides the direct contribution of understory bamboo, the difference in canopy tree species composition and density caused by different understory bamboo cover (Taylor et al. 2004; Taylor et al. 2006) may also affect the land surface phenology. Therefore, further studies on the phenology of canopy trees and understory bamboo measured on the ground are needed for understanding the phenological characteristics captured by remotely sensed data. 45 According to the jackknife analysis, the most important predictor variables containing either the most useful or the most unique information for mapping understory bamboo were the base and maximum levels, SOS, and large and small integrals (Figure 2.5). Significant differences between forest pixels with bamboo versus background pixels of the whole study area were observed in all of these five metrics (Figure 2.3). In addition, except for the base level, the other four metrics were also significantly different between forests with bamboo versus background forests (Figure 2.3). The results of these two analyses and the consistence between them indicate that the good performance of the approach developed in this study on mapping understory bamboo is due to (1) the ability of phenology metrics derived from a time series of MODIS data to capture differences in land surface phenology caused by understory bamboo, and (2) the ability of Maxent to extract and use the phenological difference for mapping bamboo distribution. Besides the good agreement between the field data and the bamboo distribution maps generated in this study, our approach has several improvements on mapping understory vegetation as compared to other methods. First, our approach is more suitable for mapping the spatial patterns and monitoring temporal dynamics of understory vegetation across large areas. While several previous approaches have been proved useful for detecting understory vegetation at local scales (Korpela 2008; Linderman et al. 2004; Resasco et al. 2007; Wang et al. 2009a; Wang et al. 2009b), their applications to broader scales may be limited due to cloud contamination (e.g., Landsat data), high acquisition costs (e.g., LiDAR data) and/or lack of images acquired during specific time periods (e.g., leaf-off seasons). In contrast, our approach solves the problem of data availability by using MODIS data, which have been acquired daily and globally since 24 Feb. 2000 and can be freely obtained. Because of the short revisiting rate (1 day), the problem of cloud contamination can be reduced by using multi-date composites. 46 With high data availability in terms of both space and time, this approach not only can be easily applied to mapping understory vegetation across large geographic areas, but also has the potential for monitoring its temporal dynamics. In addition, our approach may not be limited to specific understory species or specific areas and time periods because of its generality. Similar to some previous methods (Chastain Jr. and Townsend 2007; Resasco et al. 2007; Wang et al. 2009b), our approach is also based on the phenological difference between over- and under-story vegetation. However, our approach uses a time series of MODIS data to capture phenological differences throughout a whole year, rather than using a single image to detect differences on a specific date. Therefore, our approach does not need prior knowledge or testing on the optimal dates on which the phenological difference between over- and under-story canopy components can be detected. It also does not need retesting and adjusting the optimal dates to account for the inter-annual variability of vegetation phenology when the approach is applied to monitoring temporal dynamics (Resasco et al. 2007). In addition, as required by a GIS expert system for adjusting maps derived from remotely sensed data (Wang et al. 2009a), knowledge on the relationships between the distribution of understory vegetation and environmental variables is not necessary in our approach. Although the relationships can effectively improve the accuracy of mapping (Wang et al. 2009a), they are specific to particular vegetation types, understory species, and geographic areas. Therefore, without the requirement of specific knowledge, our approach is more general and thus is easily applicable to other vegetation types, understory species, and geographic locations. An additional advantage of our approach is its flexibility and extensibility. Although prior knowledge on the species-environment relationships is not necessary, if available, the new approach can incorporate this information easily to extend its ability to separate different 47 understory species. In this study, we showed that the approach with only the phenology metrics could not differentiate individual bamboo species effectively from the overall bamboo distribution because there was no significant difference in most phenology metrics between all bamboo pixels versus those with arrow or umbrella bamboo (Figure 2.3). However, by incorporating elevation as an additional predictor variable, we significantly improved the ability of our approach to separate the spatial distributions of the two species. Therefore, while our approach has its generality for detecting overall understory vegetation or groups of species with similar phenological characteristics across large areas, it can be applied to mapping individual species within specific areas by adding species- and/or area-specific information, such as elevation in this study. Contrasting to previous approaches which focused on either a group of similar species (Korpela 2008; Linderman et al. 2004; Wang et al. 2009a) or a single species (Resasco et al. 2007; Wang et al. 2009b), our approach provides a tool to separate individual species from a group of similar ones. This advantage would be valuable for the assessment and management of understory species biodiversity. Like any other methods, the approach developed in this study also has some limitations. First, there is always a compromise between spatial and temporal resolutions of remotely sensed data. By using WDRVI derived from MODIS data, our approach mapped understory vegetation with a spatial resolution of 250 × 250 m/pixel. Although this resolution is coarser than those of the previous approaches which use higher spatial resolution data, such as Landsat (e.g., Linderman et al. 2004), ASTER (e.g., Wang et al. 2009a) and LiDAR (e.g., Korpela 2008), a previous study has shown that a time series of MODIS data performs as well as a Landsat image on mapping wildlife habitat because the finer temporal resolution of MODIS data compensates for the disadvantage of coarser spatial resolution (Viña et al. 2008). Even though this study 48 tends to underestimate the model accuracy due to the use of presence/background, rather than presence/absence validation data, the kappa values of the overall bamboo model and the individual bamboo models with elevation generated in this study are comparable to, or even higher than, the values reported in previous studies mapping understory bamboo using higher resolution remotely sensed data (Linderman et al. 2004; Wang et al. 2009a; Wang et al. 2009b). In addition, coarse spatial resolutions may be detailed enough to reveal biologically and ecologically meaningful information in studies and applications where finer spatial resolutions are not necessary, such as those used for characterizing the habitat of wildlife species with large home ranges, as the giant pandas (Schaller et al. 1985). Second, this approach can only be applied to mapping the understory vegetation whose presence causes detectable differences in phenological characteristics. In this study, we found phenological differences between forest pixels with understory bamboo and background pixels of forests due to the high biomass and rapid growth of bamboo species. However, many understory plant species, including several non-native invasive species, have the ability to form dense understory layers and affect forest structure and function (Royo and Carson 2006; Urgenson et al. 2009). Therefore, we believe that our approach can be applied to mapping many other understory species whose presence causes differences in land surface phenology. Finally, calculating phenology metrics from a time series of WDRVI values derived from MODIS data requires more data processing time and computational resources than many previous methods. However, a global phenology product with a spatial resolution of 1 km (MOD12Q2, Zhang et al. 2003), and a 250 m product for North America (http://accweb.nascom.nasa.gov/), are being generated from MODIS data and being made freely accessible. In addition, the recent availability of software especially developed for extracting 49 phenological characteristics from remotely sensed data (e.g., TIMESAT; Jönsson and Eklundh 2004) has made the data processing easier and more efficient. With the growing interests and continuous improvements in related research on land surface phenology (Morisette et al. 2009), more data and improved tools will become available in the near future. Conservation Implications and Conclusions Understory vegetation not only has important contributions to and significant effects on the biodiversity of plant species in forest ecosystems (Gilliam 2007; Royo and Carson 2006), but also shapes the environments and provides resources for many wildlife species (Deal 2007; Díaz et al. 2005; Hagar 2007). Therefore, understanding the spatial patterns and temporal dynamics of understory vegetation is important for biodiversity conservation and habitat management. In this study, we developed an effective and practical approach for mapping understory vegetation using phenological characteristics derived from a time series of remotely sensed data. Due to the easy access, global coverage, and temporally continuous availability of MODIS data, our approach solves the problem of limited data availability that other methods may encounter when applied to larger spatial extents or finer temporal resolutions. Without the need of prior and specific information on the phenological difference between over- and under-story vegetation and on the relationships between understory vegetation and environmental variables, our approach can be easily applied to different species in different geographic areas. Due to its flexibility and extensibility, besides detecting general understory vegetation, the approach can be also used to differentiate individual species by incorporating species-specific information. The approach developed in this study could provide valuable information for ecosystem management and for biodiversity conservation. For example, while remote sensing has been 50 widely used to map and monitor the distribution of invasive plants across large spatial extents (Asner and Vitousek 2005; Huang and Asner 2009), its application for invasive understory species is limited (but see Asner et al. 2008; Resasco et al. 2007). Because of high biomass and rapid growth of many invasive species and their strong influence on species compositions, which may alter the land surface phenology, the approach developed in this study could provide a useful tool for the management of invasive understory plants at broad spatial scales. In addition, wildlife habitat management and conservation could also benefit from the new approach. Understory bamboo species, for example, are staple food for giant pandas and one of the most important factors determining the quality of giant panda habitat (Bearer et al. 2008; Liu et al. 1999; Reid et al. 1989; Schaller et al. 1985). Without the essential information on bamboo distribution, a habitat evaluation may overestimate the carrying capacity of the giant panda by more than 40% (Linderman et al. 2005). Besides providing the distribution patterns of overall understory bamboo across large areas for panda habitat evaluations, our approach can also map individual species. Because different bamboo species have unequal contributions to comprising giant pandas’ diet and determining habitat quality (Schaller et al. 1985), identifying the distribution of individual bamboo species may provide more detailed information for characterizing panda habitat. Furthermore, with the individual species information, the potential impacts on panda habitat of species-specific dynamics of understory bamboo (e.g., mass die-offs following flowering) can be incorporated into management planning. Since many other wildlife species around the world also depend on the understory vegetation whose information on spatiotemporal dynamics across large spatial extents is unavailable, habitat management and conservation of those species might benefit from the approach we developed in this study as well. 51 CHAPTER 3 TEMPORAL TRANSFERABILITY OF WILDLIFE HABITAT MODELS: IMPLICATIONS FOR PANDA HABITAT MONITORING In collaboration with Andrés Viña, Gary J. Roloff, Wei Liu, Zhiyun Ouyang, Hemin Zhang, and Jianguo Liu This chapter is based on: Tuanmu, M.-N., Viña, A., Roloff, G.J., Liu, W., Ouyang, Z., Zhang, H., & Liu, J. (2011). Temporal transferability of wildlife habitat models: implications for habitat monitoring. Journal of Biogeography, 38, 1510-1523. 52 Abstract Temporal transferability is an important issue when habitat models are used beyond the time frame corresponding to model development, but has not received enough attention, particularly in the context of habitat monitoring. While the combination of remote sensing technology and habitat modeling provides a useful tool for habitat monitoring, the effect of incorporating remotely sensed data on model transferability is unclear. Therefore, our objectives were to assess how different satellite-derived variables affect temporal transferability of habitat models and their usefulness for habitat monitoring. We modeled giant panda habitat in Wolong Nature Reserve, China with the maximum entropy algorithm using panda presence data collected in two time periods and four different sets of predictor variables representing land surface phenology. Each predictor variable set contained either a time series of smoothed wide dynamic range vegetation index (WDRVI) or eleven phenology metrics, both derived from single-year or multi-year (i.e., 3-year) remotely sensed imagery acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS). We evaluated the ability of models obtained with these four variable sets to predict giant panda habitat within and across time periods by using thresholdindependent and threshold-dependent evaluation methods and five indices of temporal transferability. Our results showed that models developed with the four variable sets were all useful for characterizing and monitoring giant panda habitat. However, the models developed using multi-year data exhibited significantly higher temporal transferability than those developed using single-year data. In addition, models developed with phenology metrics, especially when using multi-year data, exhibited significantly higher temporal transferability than those developed with the time series. This study indicates that the integration of land surface phenology, captured by high temporal resolution remotely sensed imagery, with habitat 53 modeling constitutes a suitable tool for characterizing wildlife habitat and monitoring its temporal dynamics. Using multi-year phenology metrics reduces model complexity, multicollinearity among predictor variables, and variability caused by inter-annual climatic fluctuations, thereby increasing temporal transferability of models. This study provides useful guidance for habitat monitoring through the integration of remote sensing technology and habitat modeling, which may be useful for the conservation of the giant panda and many other species. Introduction Habitat loss and degradation due to human activities and human-induced climate change have impacted and will continue to affect many animal and plant species (Sala et al. 2000). To minimize negative impacts and threats, there have been increased efforts to protect species habitats. Monitoring the spatiotemporal dynamics of species habitats is therefore essential not only for improving current conservation efforts but also for guiding future conservation strategies (Balmford et al. 2003; Lengyel et al. 2008; Pereira and Cooper 2006). Although routine field surveys can detect fine-scale changes in species habitat, they seldom provide complete spatial coverage of the areas of interest. While empirical habitat models are a useful tool for generalizing field information (Guisan and Zimmermann 2000), remote sensing technology provides synoptic information of the land surface and, in some instances, with a high temporal resolution (Turner et al. 2003). Therefore, the combination of remotely sensed data, field survey data, and habitat modeling makes it possible to map species habitats and monitor their temporal changes across large areas. Seasonal variability in biophysical characteristics (e.g., biomass) of vegetation as portrayed by multi-temporal remotely sensed data, i.e., land surface phenology (Friedl et al. 54 2006), is an important feature of the land surface for characterizing species habitat. Land surface phenology reflects different land cover types as well as different characteristics of vegetation (DeFries et al. 1995; Reed et al. 1994), and thus has been used for mapping land use change (de Beurs and Henebry 2004) and for monitoring vegetation dynamics (Beck et al. 2006; Koltunov et al. 2009). In addition, several variables representing land surface phenology have been used in habitat models for mapping plant and animal habitats at a single point in time (Morisette et al. 2006; Tuanmu et al. 2010; Viña et al. 2008). However, the usefulness of land surface phenology for monitoring the temporal dynamics of species habitat has not been assessed. When habitat models are intended to be used beyond the areas and time periods over which they were originally developed, one critical characteristic is their transferability, i.e., the ability of a model developed in one area or time period to be reliably applied in different areas or time periods. While spatial transferability has drawn increasing attention (e.g., Peterson et al. 2007; Randin et al. 2006; Zanini et al. 2009), the issue of temporal transferability has comparatively received less attention (but see Thuiller et al. 2004; Varela et al. 2009; Zharikov et al. 2009), particularly in the context of habitat monitoring. As diverse characteristics of the land surface portrayed by remotely sensed data have become increasingly available for habitat modeling (Kerr and Ostrovsky 2003; Turner et al. 2003), it is essential to assess how different use of remotely sensed data may affect model transferability. The goal of this study was to evaluate the utility of different land surface phenology variables for monitoring the temporal dynamics of wildlife habitat, particularly addressing their effects on model transferability. Using the giant panda (Ailuropoda melanoleuca) habitat as a case study, our objectives were to: (1) evaluate the predictive power and temporal transferability of habitat models derived from different land surface phenology variables; (2) identify the best 55 land surface phenology variable set for modeling, and hence monitoring giant panda habitat; (3) explore potential factors affecting model transferability; and (4) discuss implications for monitoring the temporal dynamics of wildlife habitat with the integration of habitat modeling and remotely sensed data. Materials and Methods Giant panda presence data We obtained giant panda presence data in Wolong Nature Reserve from two field datasets. The first dataset was acquired by the Third National Giant Panda Survey (State Forestry Administration 2006) during the summer of 2001. This survey covered all areas that were known or had the potential to support giant pandas (Figure 3.1). The surveyed area was 2 divided into ca. 2 km sections and each surveyor was assigned one section per day to search for and geo-reference giant panda signs (including faecal droppings, feeding sites, dens, footprints, and visual sightings) using Global Positioning System (GPS) receivers (Loucks and Wang 2004; State Forestry Administration 2006). The second dataset was obtained from wildlife surveys we conducted from August 2006 to February 2008. We followed the same procedure used in the national survey, but concentrated our survey efforts in one of the regions considered to possess the best giant panda habitat in the reserve (Liu et al. 2001; Figure 3.1). 56 Figure 3.1. Locations of panda activity signs recorded during the Third National Giant Panda  Survey in 2001 and during wildlife surveys from 2006‐2008.  Remotely sensed data We used a time series of the Moderate Resolution Imaging Spectroradiometer (MODIS) imagery (MOD09Q1) acquired between May 2000 and April 2008 for portraying phenological characteristics of vegetation. This image time series is composed of 8-day composite surface reflectance in the red (620-670 nm) and near infrared (841-876 nm) spectral bands, with a spatial resolution of ca. 250×250m/pixel. Using surface reflectance, we calculated the Wide Dynamic 57 Range Vegetation Index (WDRVI; Gitelson 2004) for each 8-day composite image. To further reduce potential effects of cloud cover on the WDRVI values, we generated a time series of 16day WDRVI composites using the maximum value between two consecutive 8-day periods. Using TIMESAT 2.3 (Jönsson and Eklundh 2004), we smoothed pixel-wise the time series of WDRVI composites by means of the adaptive Savitzky-Golay filter (Savitzky and Golay 1964). We then generated 11 phenology metrics for each of the seven full-year cycles (2001-2007) from the time series of smoothed WDRVI values to capture the shape and phenological characteristics of the smoothed curve of WDRVI values. Detailed information on the definitions and calculations of the metrics can be found in Chapter 2 and Figure 2.2. Predictor variables To assess the effects of different land surface phenology variables on temporal transferability of habitat models, we created four different variable sets and built panda habitat models for two time frames. Each variable set contained either the time series of smoothed WDRVI composites or the 11 phenology metrics, whose values were derived from single-year or multi-year (i.e., 3-year) MODIS data (Table 3.1). For the variables derived from multi-year data, the values of each WDRVI composite or phenology metric were averaged over three years. Since multi-year averages smoothed out inter-annual variability in variable values, the four variable sets allowed us to assess not only the effects of variable type (i.e., WDRVI or phenology metrics), but also those of inter-annual variability on model transferability. 58 Table 3.1. Properties of habitat models developed with four different sets of land surface phenology variables in two time frames.        a     Variables  Time Frame 1    Time Frame 2    Number of Variables (in  both time frames)    WDRVI     Single‐year  (SYVI)  Multi‐year  (MYVI)          Phenology Metrics   Single‐year  (SYPM)  Multi‐year  (MYPM)      A time series of  A time series of  Phenology metrics  Phenology metrics  smoothed WDRVI  smoothed WDRVI    in 2001  averaged over 2001‐2003 in 2001  averaged over 2001‐2003  A time series of  A time series of  Phenology metrics  Phenology metrics  smoothed WDRVI  smoothed WDRVI    in 2007  averaged over 2005‐2007 in 2007  averaged over 2005‐2007    23  23    11  11              Time Frame 1    0.87 ± 0.01  0.90 ± 0.01    0.16 ± 0.12  0.20 ± 0.12  Time Frame 2    0.87 ± 0.01  0.90 ± 0.01    0.21 ± 0.13  0.18 ± 0.13              Time Frame 1    103.5 ± 3.3  84.9 ± 2.7    39.4 ± 1.6  34.0 ± 1.4  Time Frame 2    118.8 ± 4.0  96.7 ± 3.1    41.7 ± 2.1  40.7 ± 1.4  Correlation Coefficient  b  Between Variables  Number of Terms in Final  c Models    a b  WDRVI: Wide Dynamic Range Vegetation Index;   Values were calculated between every pair of variables in the variable set and  c are shown as mean ± 2 standard errors from the mean;   Values were calculated from 20 variants of the habitat models and are  shown as mean ± 2 standard errors from the mean  59 For each land surface phenology variable set, we obtained values from MODIS data in 2001 or the averages over 2001 to 2003 to represent land surface phenology in the first time frame, and obtained values from 2007 data or averages over 2005 to 2007 to portray land surface phenology in the second time frame (Table 3.1). For single-year variable sets, phenology metrics could not be calculated for some pixels (5.7% of total pixels for 2001 and 3.4% for 2007) due to either lack of detectable seasonal cycles or incomplete cycles within a year. We excluded those pixels from further analyses. For multi-year variable sets, we excluded pixels that lacked phenology metrics in ≥ 2 years (0.8% of total pixels for both 2001-2003 and 2005-2007 periods) from further analyses. Using the values of 10,000 randomly selected pixels (ca. 26% of total pixels in the study area), we calculated Pearson’s correlation coefficients for every pair of variables, and used them as indicators of the degree of multicollinearity among variables. With this procedure, a correlation matrix was obtained for each variable set in each time frame. Analytical design The analytical design included three steps: (1) model development, (2) habitat prediction, and (3) model evaluation (Figure 3.2). Detailed methods of these steps are explained in the following sections. 60 Figure 3.2. The analytical design included three steps: (1) model development (arrow with  double lines), (2) habitat prediction (arrow with single solid line), and (3) model evaluation  (arrow with dashed line).  (1) Panda habitat models were developed with four different sets of  land surface phenology (LSP) variables in two time frames (t1 and t2) using panda presence  data collected in 2001 and 2006‐2008, respectively.  (2) The models were used to predict panda  habitat within (WTP) and beyond the time frame (BTP) in which the models were developed.  (3)  The habitat maps from WTP and BTP were evaluated using the presence data collected in the  time frame in which the habitat was predicted.  AUC and MPA were calculated for the habitat  maps from both WTP (AUCt1t1, AUCt2t2, MPAt1t1, and MPAt2t2) and BTP (AUCt1t2, AUCt2t1,  MPAt1t2, MPAt2t1).  The habitat maps from WTP and BTP were compared (arrow with dashed‐ dotted line) within each time frame, and values of agreement coefficients (ACt1 or ACt2) and  the proportions of systematic disagreement (PSDt1 or PSDt2) were calculated.  The AUC, AC and  PSD values were then used to evaluate model transferability.    61 Model development We developed panda habitat models using the maximum entropy algorithm (Maxent), a machine-learning approach for making predictions from incomplete information (Phillips et al. 2006). Maxent estimates the probability of species presence by finding the most uniform probability distribution (i.e., with the maximum entropy) as constrained by the data distribution of predictor variables associated with confirmed species locations (Phillips et al. 2006; Phillips and Dudík 2008). Maxent uses presence-only data, and thus it is especially suitable for mapping the distribution of species when confirmed absence data are difficult to obtain, as is the case for the giant panda. Maxent contrasts the values of predictor variables associated with species presence locations against the values of the same variables for all available locations (i.e., background). We randomly selected 10,000 pixels as a representation of the entire study area (Phillips and Dudík 2008), and defined the background by only using pixels where giant pandas could possibly occur following the suggestion of Phillips et al. (2009) for single species applications. For this, as giant pandas seldom occur in non-forest areas (Schaller et al. 1985), the background was defined by pixels with forest cover according to a binary forest cover map which was derived from a Landsat-5 TM image acquired on 13 June 2001 (Viña et al. 2007) and resampled to the spatial resolution of the MODIS data (i.e., 250 m). Maxent derives and uses different forms of input variables (i.e., feature types) to represent non-linear and interactive effects of predictor variables on species presence probability (Phillips et al. 2006). The contributions of these derived predictors to the model prediction are then evaluated during model development and only those having significant contributions are retained in a final model (Phillips et al. 2006). We used a combination of linear, quadratic, and 62 product feature types, which represent the means, variances, and co-variances of the predictor variables, respectively (Phillips et al. 2006). The number of terms retained in the final model was used for indicating the complexity of model structure. We developed four habitat models in two different time frames (i.e., Time Frame 1 and Time Frame 2, Figure 3.2) using our two panda presence datasets. We considered a 250×250 m pixel as a confirmed presence pixel if it contained at least one panda location according to field surveys. Field datasets contained 399 and 220 presence pixels in 2001 and 2006-2008, respectively. For each dataset, we used 70% of the presence pixels for model development and the remaining 30% for model evaluation (see below). In order to reduce the effects of data partitioning on model outputs, we randomly re-partitioned the data and created 20 different data partitions for each field dataset. Twenty variants of each model were then developed using these 20 partitions for each time frame. This number of partitions was used because a previous study showed that the variation of model outputs decreased with the increase in the number of partitions, but changed negligibly with more than 20 partitions (see Chapter 2). Habitat prediction For each time frame, 20 variants of each model were used to predict panda habitat within the time frame (i.e., within-time-frame prediction, WTP) and beyond the time frame (i.e., beyond-time-frame prediction, BTP; Figure 3.2). Therefore, in each time frame, 80 panda presence probability maps (i.e., 4 models × 20 variants) were obtained from the WTP and 80 additional maps were obtained from the BTP (Figure 3.2). 63 Model evaluation To evaluate the accuracy of model predictions, we used both threshold-independent and threshold-dependent methods. The threshold- independent method consisted in the receiver operating characteristic (ROC) curve analysis, a common method for evaluating the accuracy of classification models (Fielding and Bell 1997). The area under the ROC curve (AUC) provides a single-value measurement of model accuracy, with a value of 1 indicating a perfect prediction and 0.5 indicating a random prediction (Hanley and Mcneil 1982). In habitat modelling, models with AUC values higher than 0.7 are considered useful (Boyce et al. 2002). The ROC curve analysis is typically conducted by contrasting presence to absence data. Here, we calculated AUC values by contrasting presence pixels to those randomly selected from the study area (i.e., background pixels) as suggested by Phillips et al. (2006) and used in other studies (e.g., Marini et al. 2010). For each time frame, we used 30% of panda presence pixels and the background pixels selected during model development to calculate AUC values for both WTP and BTP (Figure 3.2). By contrasting presence to background data in the ROC curve analysis, the maximum achievable AUC value is less than 1 and is negatively correlated with the proportion of actual presence pixels in the background pixels, but the value of a random prediction is still equal to 0.5 (Phillips et al. 2006). Some concerns have been raised on the application of the ROC curve analysis for comparing the accuracy of models using different modelling approaches or among different species, because AUC values are subject to the range of model output values, reliability of species presence and absence data, and the delineation of a study area, especially when it is used to define the background in presence-only models (Lobo et al. 2008; Peterson et al. 2008). However, comparisons among different AUC values performed in this study are considered to be 64 valid since all models were generated for the same species, by the same modelling algorithm (i.e., Maxent), within the same study area, and were generated and validated with the same presence and background data (see below). The threshold-dependent evaluation method used was based on the calculation of the minimal predicted area (MPA; Engler et al. 2004). The MPA method evaluates model performance based on the parsimony concept that a good model should predict the smallest habitat area as possible (i.e., minimize commission errors as much as possible), while its omission errors are under control (Engler et al. 2004). Since the MPA depends on the actual proportion of species habitat in a study area, which is almost always unknown, it is only suitable for comparing models generated for the same species in the same area, as is the case of this study. In addition, it is relative magnitudes of MPA among models, rather than absolute values for individual models, that matter for evaluating model performance. Following Engler et al. (2004), we defined a threshold for each panda presence probability map so that 90% of presence locations in the validation dataset were encompassed (i.e., 10% omission error). Instead of using absolute area, we calculated the MPA as the ratio of the number of above-threshold pixels to the total number of pixels. Model transferability We evaluated the temporal transferability of panda habitat models based on three criteria adapted from Randin, et al. (2006). First, a model with good temporal transferability should have similar accuracy between its predictions within and beyond the time frame corresponding to its development. Therefore, for our analytical design (Figure 3.2), there should be similar accuracy between WTP in Time Frame 1 (t1) and BTP in Time Frame 2 (t2), as well as between 65 WTP in t2 and BTP in t1. Second, the model should have similar performance no matter in which time frame it was developed. That is, the accuracy of WTP (or BTP) should be similar between t1 and t2 (Figure 3.2). This criterion, together with the first criterion, implies that the model should have similar transferability in both transferring directions between the two time frames. Third, besides model accuracy, the spatial patterns of predicted habitat from within- and beyond-time-frame predictions also should be similar. Therefore, the spatial patterns of WTP and BTP in t1 (or t2) should match each other (Figure 3.2). To quantify temporal transferability based on these three criteria, we calculated five indices for each of the four habitat models. We calculated the single-direction (TIt1→t2 and TIt2→t1) and overall transferability indices (TIoverall), which were adapted from Randin et al. (2006) as: TIt1→t2 = 1 - |AUCt1t1 – AUCt1t2| / 0.5; (3.1) TIt2→t1 = 1 - | AUCt2t2 – AUCt2t1| / 0.5; (3.2) TIoverall = [0.5 × (TIt1→t2 + TIt2→t1)] / 1 + | TIt1→t2 - TIt2→t1|, (3.3) where AUCt1t1 and AUCt2t2 are AUC values for WTP in Time Frame 1 and Time Frame 2, respectively, and AUCt2t1 and AUCt1t2 are for BTP in the two time frames, respectively (Figure 3.2). TIt1→t2 and TIt2→t1 measure the ability of a model to be transferred from t1 to t2 and vice versa, respectively. They range from 0 to 1 as AUC values are typically between 0.5 and 1, and they are closer to 1 when the AUC values for WTP and BTP are similar (i.e., high transferability based on the first criterion). TIoverall measures transferability in both directions and puts a penalty on the difference between two directions. It also ranges from 0 and 1 and is closer to 1 66 when single-direction transferability indices in both directions are higher and closer to each other (i.e., high transferability based on the second criterion). We compared the spatial patterns of WTP and BTP in each time frame (Figure 3.2), by calculating an agreement coefficient (AC; Ji and Gallo 2006) for each variant of each habitat model as: AC = 1 - ∑ 2 Xi – Yi) / ∑ | - | + | Xi - |) × (| - | + Yi - |)], (3.4) where Xi and Yi are pixel values (i.e., estimated panda presence probability) of habitat maps from WTP and BTP, respectively; X and Y are the mean values of Xi and Yi, respectively; and n is the total number of pixels. AC is a standardized sum of squared difference between Xi and Yi, and its maximum value is 1, indicating a perfect agreement in pixel values between two maps (Ji and Gallo 2006). Because two maps may have a low AC value even when they show the same spatial patterns with different absolute pixel values, we also calculated the proportion of systematic disagreement (PSD), following Ji and Gallo (2006) as: PSD = 1 - ∑ |Xi - i|) × (|Yi - i|)] / ∑ 2 Xi – Yi) , (3.5) where i and i are the estimated values of Xi and Yi, respectively, from a linear regression between Xi and Yi, based on a geometric mean functional relationship model. The denominator of the main term of equation 3.5 measures the total disagreement between Xi and Yi, and the numerator is the sum of residuals from the regression line, which indicates non-systematic or random disagreement. As PSD equals one minus the ratio of non-systematic disagreement to 67 total disagreement, this index measures the proportion of systematic disagreement captured by a linear regression between the pixel values in two maps (Ji and Gallo 2006). In our case, a higher PSD indicates that a larger proportion of the disagreement between two habitat maps is due to a linear shift of pixel values, and thus the two maps show a more similar spatial pattern of habitat but just with different absolute pixel values. Model comparisons We conducted the mixed-design analysis of variance (ANOVA) for comparing AUC and MPA values among habitat models with model variants as a random factor, and used paired ttests for pairwise comparisons. With the same random factor, we used two-way ANOVA to evaluate the effects of land surface phenology variable type and the length of original time series data used for generating them, as well as their interaction effects on model transferability measured by TI (Eqs. 3.1-3.3), AC (Eq. 3.4) and PSD (Eq. 3.5). We used these parametric statistical tests after verifying the validity of the normality assumption with Shapiro-Wilk tests. All statistical tests were conducted using R 2.10.1. Results Multicollinearity and model complexity The correlation analysis among land surface phenology variables in each variable set showed that smoothed time series of WDRVI values, regardless of single- or multi-year data, were highly correlated in both time frames, but phenology metrics were less correlated with each other (Table 3.1). In both time frames, the models developed with phenology metrics or variables derived from multi-year time series data tended to have fewer terms, suggesting lower 68 model complexities (Table 3.1). In addition, each of the four models tended to contain more terms when developed in Time Frame 2 than in Time Frame 1 (Table 3.1). Model accuracy According to the threshold-independent evaluation, all habitat models had median AUC values ranging between 0.85 and 0.95 for predicting giant panda habitat within the time frame of model development (i.e., WTP; Figure 3.3a and c). In general, accuracy decreased when models were used to predict habitat beyond time frames (i.e. BTP), but all median AUC values were still higher than 0.79, indicating that they constitute useful models (Figure 3.3b and d). Significant differences in AUC values were found among the four habitat models in both time frames (F = 11.96, 267.06, 194.53 and 489.49 for AUCt1t1, AUCt1t2, AUCt2t2 and AUCt2t1, respectively; d.f. -5 = 3, 57 and p < 10 for all). In Time Frame 1, the MYVI and MYPM models had the highest predictive power within the time frame (Figure 3.3a), but the MYPM model was significantly better than the MYVI model for beyond-time-frame predicting (Figure 3.3b). In Time Frame 2, the MYVI model was significantly better for both WTP and BTP, although the ranking of the MYPM model improved when it was used for predicting habitat beyond the time frame (Figure 3.3c and d). The threshold-dependent evaluation procedure showed very similar patterns of model accuracy among the four models (Figure 3.3e-h; note the reversed y-axes). The major difference from the threshold-independent evaluation was found when the models developed in Time Frame 1 were used for predicting habitat in Time Frame 2. While the MYPM model was the best in terms of AUC values, the SYVI and MYVI models were as good as the MYPM model in terms of MPA values (Figure 3.3b and f). 69 t1t1 t2t2 t2t1 t1t2 MYPM  SYPM  MYVI  SYVI  MYPM  SYPM  t2t2 t2t1 MYPM  SYPM  MYVI  SYVI  MYPM  SYPM  MYVI  SYVI  MYPM  SYPM  MYVI  SYVI  MYPM  SYPM  MYVI  SYVI  0.4 0.3 MPA 0.2 0.1 t1t1 MYVI  SYVI  MYPM  SYPM  MYVI  SYVI  MYPM  SYPM  MYVI  SYVI  0.80 AUC 0.85 0.90 0.95 t1t2 Figure 3.3. Box plots of the area under the receiver operating characteristics curve (AUC) and  the minimal predicted area (MPA) for the panda habitat models developed with four different  land surface phenology variable sets (SYVI, MYVI, SYPM and MYPM in Table 3.1) in two time  frames (t1 and t2).  The AUC and the MPA values were calculated from 20 variants of each  habitat model when the model was developed in t1 and used to predict habitat in t1 [(a) and (e)]  70 and t2 [(b) and (f)], as well as developed in t2 and used to predict habitat in t2 [(c) and (g)] and  t1 [(d) and (h)].  Each box plot shows the maximum, 75th percentile, median, 25th percentile  and minimum values.  Note that y‐axes for the MPA are reversed since lower MPA values  indicate higher model accuracies.  The letters above box plots indicate the results of pair‐wise  comparisons conducted using paired t‐tests.  The alphabetical order shows the order of model  accuracy from high to low.  There is no significant difference (p > 0.05 after the Holm– Bonferroni adjustment) between two models if they share the same letter.    Model transferability Single-direction and overall transferability indices indicated different transferability among the four habitat models. The MYPM model had the highest values of all three transferability indices (Figure 3.4), indicating it was the most transferable among the four models in terms of model accuracy. Results of the two-way ANOVA showed that both the variable type and length of original time series data had significant main effects on model transferability (Table 3.2). The models developed with phenology metrics were more transferable than those with time series WDRVI, and the models developed with the variables derived from multi-year data had higher transferability than those with the variables from single-year data (Table 3.2). Comparing the two single-direction transferability indices indicated that model transferability was also different between transferring directions (Figure 3.4a and b). All models were less transferable from Time Frame 2 to Time Frame 1 than in the opposite direction (Figure 3.4a and b). Although the variable type and data length had the same effects on transferability in both transferring directions, their relative magnitudes were different. While the variable type 71 was more influential on TIt1→t2, the data length was more influential on TIt2→t1 (Figure 3.4a and b). Regarding the match between spatial patterns of WTP and BTP, significant differences in the agreement coefficients (AC) and the proportions of systematic disagreement (PSD) were found among models (Figure 3.5). For both time frames, the habitat maps produced by the MYPM model showed the most similar spatial patterns, indicated by both the highest agreement in predicted panda presence probabilities (i.e., highest AC) and the highest proportion of disagreement that can be captured by a linear regression (i.e., highest PSD; Figure 3.5). This indicates that the MYPM was the most transferable among the four models in terms of the spatial patterns of model predictions. Land surface phenology variable type and length of original time series both had significant effects on AC and PSD values (Table 3.3). Using the variables derived from multiyear data as predictors increased both AC and PSD values. Significant interaction effects on PSD (Table 3.3) indicated that the effect of data length was more influential on the models based on phenology metrics (Figure 3.5c and d). The variable type affected AC and PSD differently. While using phenology metrics as predictors increased the agreement in the pixel values of habitat maps (i.e., AC), it reduced the proportions of systematic disagreement in habitat predictions (i.e., PSD; Table 3.3). However, the significant interaction effects on PSD (Table 3.3) indicated that the negative effect of phenology metrics was not influential when multi-year data were used to generate the metrics (Figure 3.5c and d).       72 Figure 3.4. Mean values of the single‐direction and overall transferability indices for the four habitat models.  These models were  developed with either a time series of WDRVI (VI) or phenology metrics (PM), which were derived from either single‐year (SY) or  multi‐year (MY) MODIS data.  The indices were calculated from 20 variants of each habitat model for evaluating: (a) the ability of the  model developed in Time Frame 1 (t1) to predict panda habitat in Time Frame 2 (t2), (b) the ability of the model developed in t2 to  predict habitat in t1, and (c) the overall ability of the model to predict habitat beyond time frames.  The error bars indicate 2  standard errors from the mean.  The lines do not imply any linear relationship, but are shown just for helping visualise value  differences.  73 Table 3.2. Analyses of variance on the effects of variable type (Type) and the length of original time series data (Length), as well as  their interaction effects (Type × Length) on the single‐direction and overall transferability indices calculated for the habitat models  developed with different land surface phenology variables.          Difference *  Type    PM > VI  22.3  < 10   Length    MY > SY  113.2    2.9  Type ×Length    *  TI t1→t2  F‐value  TI t2→t1    *  p‐value    Difference  ‐3 F‐value  *  p‐value    Difference  ‐12   PM > VI  313.2  < 10 < 10     MY > SY  101.0  < 10   0.11      12.2  0.002  ‐8 TI overall    ‐8     p‐value  ‐11 PM > VI  238.7  < 10   MY > SY  42.9  < 10       14.0  0.002  PM: phenology metrics; VI: a time series of WDRVI values; MY: multi‐year data; SY: single‐year data 74 F‐value  ‐5   Figure 3.5. Mean values of agreement coefficients (AC) and the proportions of systematic  disagreement (PSD) for the four habitat models.  These models were developed with either a  time series of WDRVI (VI) or phenology metrics (PM), which were derived from either single‐ year (SY) or multi‐year (MY) MODIS data.  The values of AC and PSD were calculated between  the maps generated from within‐ and beyond‐time‐frame predictions for the panda habitat in  Time Frame 1 ((a) and (c), respectively) and in Time Frame 2 ((b) and (d), respectively).  The  error bars indicate 2 standard errors from the mean.  The lines do not imply any linear  relationship, but are shown just for helping visualise value differences.  75 Table 3.3. Analyses of variance on the effects of variable type (Type) and the length of original time series data (Length), as well as  their interaction effects (Type × Length) on agreement coefficients (AC) and proportions of systematic disagreement (PSD), which  were calculated from within‐ and beyond‐time‐frame predictions of panda habitat in the two time frames.          AC    Type  Time Frame 1  *  F‐value  p‐value              PM > VI  282.5  < 10‐   Length    MY > SY  339.6  < 10 Type × Length      6.64        Type    VI > PM  220.7  < 10 Length    MY > SY  365.1  < 10 Type × Length      552.7  < 10 PSD  *  Difference    Time Frame 2  *  F‐value  p‐value          PM > VI  192.8  < 10   MY > SY  218.4  < 10 0.018      0.14  0.717                VI > PM  19.2  < 10       MY > SY  2079.6        72.5  12 ‐12 ‐11 ‐13 ‐14   Difference  PM: phenology metrics; VI: a time series of WDRVI values; MY: multi‐year data; SY: single‐year data 76 ‐10 ‐11     ‐3 ‐16 < 10 ‐7   < 10   Discussion Model transferability is an important characteristic of empirical habitat models when they are used beyond the area and time in which they were originally developed. Diverse factors may affect spatial or temporal transferability of habitat models, such as the environmental variability of species habitat captured by field data (Bamford et al. 2009; Phillips 2008; Thuiller et al. 2004), relevance of predictor variables for describing underlying processes affecting species distribution, especially for those that vary across space or over time (Vanreusel et al. 2007; Zharikov et al. 2009), the ability of modelling methods to capture species-environment relationships (Araújo et al. 2005; Randin et al. 2006), and consistency of the relationships across space (Zanini et al. 2009) and time (Pearson and Dawson 2003). By controlling other factors, such as focusing on a single species, using the same field data for model development and validation, and using the same modelling algorithm, we showed that the methods for generating predictor variables, in particular from remotely sensed data, has significant effects on model transferability. Effects of variable type on model transferability Panda habitat models developed with phenology metrics, especially when they were generated using multi-year data, exhibited higher temporal transferability than those developed with time series of WDRVI in terms of both model accuracy and the spatial match in model predictions. The advantage of using phenology metrics as predictor variables on model transferability may be related to the reduction of (1) model complexity and (2) multicollinearity among variables. Our results showed that the models developed with time series of WDRVI had more complex structure than the models with phenology metrics. More complex models have a higher 77 chance of over-fitting the training data, and thus to lose the ability of capturing general relationships between species occurrence and predictor variables and lose their transferability (Araújo and Guisan 2006; Randin et al. 2006). Therefore, the lower risk of over-fitting caused by less complex models may be a probable reason for the higher transferability of the models developed with phenology metrics. In addition, the models developed in Time Frame 2 tended to be more complex than Time Frame 1 models. We also found that temporal transferability was lower and the effects of variable type were stronger when the models were applied from Time Frame 2 to Time Frame 1. The lower transferability of Time Frame 2 models may be due to an incomplete environmental range of panda habitat captured by field data (Thuiller et al. 2004), since the panda presence data were collected only in high quality habitat within Time Frame 2. However, the relationships among higher model complexity, lower transferability, and stronger effects of variable type for Time Frame 2 models suggest that the lower transferability could also be attributed to higher over-fitting risks caused by more complex Time Frame 2 models, and using phenology metrics in a habitat model can reduce these risks and thus increase its transferability. In regression models, multicollinearity influences the estimation of coefficients and their standard errors, affects significance tests on the coefficients, changes model structure, and thus reduces the robustness of the established species-environment relationships in habitat models (Graham 2003; Mac Nally 2000). Although no significance tests on coefficients are involved in model development of non-regression-based modelling approaches such as decision trees (Berk 2006) or Maxent, models developed with these approaches may not be completely free of multicollinearity problems (Mac Nally 2000). When models are used for predicting species habitats for different areas or time periods, spatially or temporally inconsistent correlations 78 among predictor variables may reduce the predictive ability of the established speciesenvironment relationships, even if they were determined with statistically robust approaches. Therefore, multicollinearity reduction by using phenology metrics may also be the reason for their positive effects on model transferability. Although several approaches can be used to reduce over-fitting and multicollinearity, they are not suitable for the habitat models developed with land surface phenology variables. For instance, one common approach is to remove highly correlated variables based on their biological importance for the species of interest (Graham 2003). However, for a time series of remotely sensed data representing land surface phenology, all variables are highly correlated (e.g., the lowest correlation coefficient was 0.73 in this study). Therefore, even if only two variables are selected from the time series, multicollinearity may still exist. In addition, with the selection of fewer variables, less information on land surface phenology can be represented. Principal component analyses are also commonly used for solving multicollinearity and over-fitting problems (Aguilera et al. 2006; Graham 2003), and are useful in models predicting habitat within the same area and time frame corresponding to model development (e.g., Viña et al. 2010). However, the combinations of predictor variables in principal components are datadependent and thus subject to change through time (Schowengerdt 2007). As variable values change over time, the inconstancy of principal components may limit their utility for predicting habitat beyond the time frame of model development. Fixed combinations of variables, like the Tasselled-cap transformation commonly used in digital image processing of remotely sensed imagery (Schowengerdt 2007), may be useful for predicting habitat changes over time. However, finding general and meaningful tasselled-cap components that reflect the underlying processes determining habitat quality and driving its change is quite challenging. 79 Effects of time series length on model transferability Our results showed that length of original time series data had significant effects on both model accuracy and transferability. The models developed with variables derived from multiyear remotely sensed data had higher accuracy for both within- and beyond-time-frame predictions and had higher temporal transferability than those with variables derived from singleyear data. The major advantage of using multi-year data is in smoothing inter-annual variability of the biophysical and phenological characteristics of the vegetation caused by inter-annual climatic fluctuations (Ichii et al. 2002), which may not reflect real habitat changes. For example, while higher vegetation index values in one year may correspond to higher temperatures (Ichii et al. 2002), this may not indicate long-term changes in the quality of panda habitat. Thus, a habitat model that uses land surface phenology variables derived from single-year data may be affected by the particularities of that year and thus lose its temporal transferability. In addition, different plant species may have different responses to inter-annual fluctuations of climatic conditions, which may cause more non-systematic disagreements in habitat predictions of the models with variables derived from single-year data. High sensitivity to inter-annual variability of vegetation characteristics is not specific to models developed with land surface phenology variables. Any model using remotely sensed data to reflect vegetation information may have the same problem when used for studying temporal dynamics of species habitat. While average values of climatic variables over several years are commonly used in habitat models to reflect long-term climatic conditions, variables derived from single-year remotely sensed data or even a single image are often used in habitat models (e.g., Zimmermann et al. 2007), mostly due to limited availability of remotely sensed data. While incorporating remote sensing variables into habitat models can improve model accuracy for the 80 time period corresponding to model development (Zimmermann et al. 2007), predicting species habitat across time frames and using these results to monitor habitat dynamics should be done with caution. As some high temporal resolution remotely sensed data (e.g., MODIS data) and multi-temporal data of median spatial resolution (e.g., Landsat TM, ETM, ALOS) are becoming increasingly available, variables derived from multi-season and multi-year data appear to be more appropriate for monitoring temporal dynamics of species habitat at broader spatial and longer temporal scales. Usefulness of land surface phenology for mapping and monitoring species habitats Our study showed that land surface phenology is useful for characterising giant panda habitat and also monitoring its temporal dynamics. Forest cover, understory bamboo, topography, and human disturbances are the most important documented determinants of panda habitat (Bearer et al. 2008; Liu et al. 2001). Land surface phenology not only reflects different land cover types and their dynamics (Beck et al. 2006; de Beurs and Henebry 2004), but it also reflects the characteristics of understory bamboo occurring under the forest canopy (Tuanmu et al. 2010; Viña et al. 2008). In addition, human disturbances on panda habitat are usually associated with land cover or vegetation change (Bearer et al. 2008; Liu et al. 2001). Therefore, besides capturing the characteristics of vegetation that is suitable for the giant panda (Viña et al. 2008; Viña et al. 2010), land surface phenology may also capture its temporal dynamics due to human disturbances. While previous studies have found the usefulness of land surface phenology for detecting vegetation changes due to human and natural disturbances (Eklundh et al. 2009; Koltunov et al. 2009), this study showed that changes in land surface phenology could be directly linked to changes in wildlife habitat through the use of habitat models. 81 However, the usefulness of different land surface phenology variables depends on their application. Since the habitat model developed with a time series of WDRVI derived from multi-year MODIS data (i.e., MYVI model) can produce the most accurate habitat maps in the time frame when the model was developed, it is a good tool for evaluating the habitat conditions in that particular time frame (e.g., Viña et al. 2010). Alternatively, the model developed with multi-year phenology metrics (i.e., MYPM model) reduced the problem of multicollinearity and the risk of over-fitting, and thus it appears to be the best in terms of temporal transferability. Therefore, the MYPM model constitutes a suitable tool for monitoring the temporal dynamics of giant panda habitat and providing essential information for the conservation of the species. Under changing environments, monitoring the temporal dynamics of species habitats at regional or global scales is essential for reducing biodiversity loss and maintaining sustainable ecosystem services (Balmford et al. 2003; Lengyel et al. 2008; Pereira and Cooper 2006). Combining remote sensing and habitat modelling provides a practical and efficient tool for monitoring temporal dynamics of biodiversity and species habitats at different spatial and temporal scales (Lengyel et al. 2008; Pereira and Cooper 2006). In particular, land surface phenology has been found to be sensitive to vegetation changes due to short-term human and natural disturbances (Eklundh et al. 2009; Koltunov et al. 2009) and long-term climate changes (Morisette et al. 2009; Zhang et al. 2004). Phenology-based models have also been successfully applied to predicting species habitat at different spatial scales (Morisette et al. 2006; Tuanmu et al. 2010; Viña et al. 2010). Therefore, the combination of land surface phenology and habitat modelling constitutes an excellent tool for biodiversity conservation under changing environments. 82 Although we used giant panda habitat in Wolong Nature Reserve as a case study, the approaches and conservation implications of this study can go beyond this specific species, geographical area and spatial scale. Previous studies have shown considerable variability in model transferability among species (Randin et al. 2006). The direct causes indicated in this study underlying the differences in transferability among models (i.e., model complexity, multicollinearity among variables, and relevance of variables to habitat quality and its change) have also been reported in other studies (Mac Nally 2000; Peterson et al. 2007; Randin et al. 2006; Vanreusel et al. 2007; Zharikov et al. 2009). Therefore, we believe that the suggestions provided for increasing model transferability (i.e., using phenology metrics and multi-year remotely sensed data) can be generally applied for modelling the habitat of many other species in different geographical settings. This is important, since model transferability cannot be directly evaluated for many species, particularly endangered species, due to low availability of field data collected over multiple years. In such cases, habitat models developed using remotely sensed data may still be useful for habitat monitoring if the suggestions provided in this study are taken into consideration. 83 CHAPTER 4 INCENTIVE-BASED AND LOCAL-ENGAGED INSTRUMENT FOR SUCCESSFUL CONSERVATION OF PANDA HABITAT In collaboration with Andrés Viña, Wu Yang, Xiaodong Chen, Ashton M. Shortridge, and Jianguo Liu 84 Abstract Conflicts between local people’s livelihoods and conservation goals have led to many failures of conservation practices and created a debate between establishing strict protected areas and improving local people’s livelihoods. Here we show that an incentive-based conservation instrument which encourages active participation of local people in natural resource management may provide a solution for the debate. Our empirical and spatially explicit assessment indicates the effectiveness of the instrument at conserving giant panda habitat and shows a positive contribution of local engagement to a panda habitat recovery, which has reversed a more-than30-year trend of panda habitat degradation in a world-renowned nature reserve. This study suggests that the implementation of this nationwide conservation program may achieve a greater overall benefit to ecosystems and their services by encouraging local engagement through economic incentives and social norms. It also provides implications for solving people – park conflicts not only in China, but around the world. Introduction The establishment of protected areas has long been the leading instrument for protecting biodiversity and ecosystem services worldwide (Millennium Ecosystem Assessment 2005; Naughton-Treves et al. 2005). However, by limiting or entirely excluding human access to natural resources without sufficient respect for local people’s right to use natural resources, this “fences-and-fines” strategy has caused negative social and economic impacts on the people living in and around protected areas (Adams et al. 2004; McShane et al. 2011). Since more lands and seas are being covered by protected areas in response to the crisis of global biodiversity loss (Butchart et al. 2010) and both human population and resource demand are growing (Millennium 85 Ecosystem Assessment 2005), the conflicts between natural resource exploitation and conservation are getting more serious. In response to these conflicts, along with a paradigm shift from viewing humans separately from nature to viewing humans as an important component in coupled human and natural systems (Berkes 2004; Liu et al. 2007), people-oriented conservation is rapidly becoming popular. By providing alternative livelihood options that cause less pressure on biodiversity and lead to sustainable use of natural resources (e.g., agroforestry and ecotourism), community-based conservation programs, including integrated conservation and development projects (ICDPs), have been developed to reduce resource-depleting activities (Berkes 2004; Kremen et al. 1994). Furthermore, programs of payments for ecosystem services (PES) provide direct incentives for local people to reduce resource extraction and even actively participate in conservation by compensating the costs of forgone livelihoods (Ferraro and Kiss 2002; Pattanayak et al. 2010). While these people-oriented instruments are intended to simultaneously protect biodiversity and sustain human livelihoods, their effects on natural resource conservation have been mixed (Andersson and Gibson 2007; Hughes and Flintan 2001; Pattanayak et al. 2010). The mixed effectiveness has resulted in a call for stricter management of protected areas and has provoked the heated “park vs. people” debate (Miller et al. 2011). This debate is particularly relevant to China, one of the most biodiverse and populous countries in the world (Liu 2010). In response to biodiversity loss, the number and the spatial coverage of protected areas in China have increased exponentially since the 1980s (Liu and Raven 2010). While conventional “fences-and-fines” and top-down management is prevalent among these protected areas (Liu and Diamond 2008), the livelihoods of tens of millions of poor rural people living in and around protected areas are negatively affected (An et al. 2001; Xu and 86 Melick 2007). Without adequate consideration of local people’s dependency on natural resources and potential people – park conflicts, failures in biodiversity conservation are common in China’s protected areas, even in flagship protected areas (Liu et al. 2001). At the end of the last century, the Chinese government implemented one of the largest PES-like policies in the world, the Natural Forest Conservation Program (NFCP), to protect and restore natural forests through logging bans and afforestation (Liu et al. 2008). NFCP provides payments to forest enterprises and local governments to compensate their economic losses due to a shift from timber harvesting to afforestation and forest management. With NFCP payments, ca. 0.7 million former logging and timber-processing workers have retired, obtained jobs in other economic sectors, or been hired for tree plantation and forest monitoring (Yin and Yin 2010). Since the implementation of NFCP in 1998, a considerable amount of natural forest has been protected and many new forested areas have been created (Liu et al. 2008; Yin and Yin 2010). While this is believed to be beneficial for biodiversity conservation (Liu et al. 2008; Loucks et al. 2001), it is unclear whether, and to what extent, the increase in forest cover translates into improved habitats for forest species. It is also unclear whether this PES-like program can be a solution for improving conservation in many protected areas. To address the issues raised above, it is essential to conduct rigorous empirical assessments that disentangle the effects of one instrument from the effects of other instruments and confounding factors (Andersson and Gibson 2007; Pattanayak et al. 2010). In this study, we conducted the first empirical and spatially explicit assessment of the effectiveness of NFCP in conserving the habitat of one of the most endangered species (the giant panda; Ailuropoda melanoleuca) in a world-renowned protected area (Wolong Nature Reserve). Our assessment showed an improvement of panda habitat after the NFCP implementation, which reversed a 87 more-than-30-year habitat degradation trend observed in the reserve (Liu et al. 2001). After controlling for confounding factors, our assessment also suggests that an innovative NFCP implementation, which encourages active participation of local residents in forest monitoring, is particularly effective at protecting panda habitat. These findings not only have implications for the NFCP implementation and panda conservation in China, but also for the “park vs. people” debate on biodiversity conservation around the world. Materials and Methods NFCP implementation in Wolong Nature Reserve The Natural Forest Conservation Program has been implemented in the reserve since 2000. While in most of China NFCP only involves state-owned forestry enterprises and local governments (Liu et al. 2008), in Wolong Nature Reserve it involves both the local government and local residents. Approximately one third of the total NFCP monitoring area (ca. 400 out of 2 1,200 km ) is assigned to ca. 250 household groups of various sizes (ranging from 1 to 16 households per group) for monitoring activities, and the remaining area is monitored by government officials (Figure 4.1). Each participating household in Wolong and Gengda receives an annual payment of ca. US$110 (ca. 8% of household income in 2001) for monitoring an assigned forest parcel, while the households in Sanjiang receive about half of that amount. All households within a monitoring group suffer payment reductions as punishment for any anthropogenic damage found in their co-monitored forest parcel. The differences in monitoring types (i.e., household vs. government monitoring) and payment levels make this reserve an excellent place to examine how different NFCP implementations affect conservation outcomes. 88 Gengda Households Townships Reserve Householdmonitored Wolong HIS Change 0.7 -0.7 Sanjiang Figure 4.1. Spatial patterns of giant panda habitat change from 2001 to 2007 in Wolong Nature  Reserve, China.  Changes in habitat suitability index (HSI) were calculated from the outputs of  the panda habitat model built with multi‐year phenology metrics in Chapter 3.  The locations of  households and household‐monitored forest parcels in the three townships (Gengda, Wolong  and Sanjiang) comprising the reserve are also shown.  89 Spatiotemporal dynamics of panda habitat We used a satellite-based habitat model to obtain the values of panda habitat suitability index (HSI ranging from 0 to 1 with higher values indicating higher suitability) in 2001 and 2007 for each pixel (250×250 m) across the study area. The habitat model was built based on panda feces locations obtained during a field survey in 2001 and land surface phenology metrics derived from a time series of remotely sensed imagery acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS) during multiple years (i.e., 2001 – 2003; Chapter 3). This model is particularly useful for monitoring changes in panda habitat as it captures information on the distribution of the most important determinants of panda habitat (i.e., forest cover and understory bamboo) and exhibits high accuracy at estimating HSI values across the study area [area under the receiver operating characteristic curve (AUC) is 0.853 and 0.855 for 2001 and 2007, respectively]. Details on the data, the metrics, and the modeling and validation approaches are reported in Chapter 3 and Tuanmu et al. (2011). We calculated HSI change (i.e., 2007 value minus 2001 value) for each pixel across the study area. We also estimated the areal change in suitable habitat in the entire reserve by applying a threshold to convert continuous HSI values into binary outcomes (habitat or nonhabitat). We used a threshold corresponding to a 10% omission error. While commission errors are, in general, negatively related to omission errors, we did not explicitly consider commission errors because only confirmed presence data (i.e., locations of panda feces) were available for model evaluation. To assess the effect of choosing different thresholds on results, we also calculated habitat areas using thresholds corresponding to 5 and 15% omission errors. 90 Effects of NFCP implementation To analyze the effects of NFCP on the spatiotemporal dynamics of panda habitat, we used regression models to spatially relate observed HSI changes with different NFCP implementations (i.e., government monitoring or household monitoring with high or low NFCP payments; Table 4.1) at the pixel level. To control for potential confounding effects, we included in the models several biophysical and anthropogenic factors (Table 4.1), which have been found important for determining the spatial patterns and temporal dynamics of panda habitat in the reserve (Bearer et al. 2008; He et al. 2009; Liu et al. 1999; Viña et al. 2011). We also included the HSI values in 2001 to account for the potential dependency of HSI changes on initial values. Information on the data and processing approaches to obtain these factors is provided in Table 4.1. We randomly selected 3,000 pixels (ca. 15%) from the study area below the tree line (i.e., 3,600m in elevation) for generating regression models because the area above the tree line does not constitute giant panda habitat (Schaller et al. 1985). Using the ordinary least squares (OLS) method, we generated a regression model with HSI change being the dependent variable and the above mentioned being the independent variables (standardized prior to model generation) to access the partial effect of NFCP monitoring types (i.e., government monitoring vs. household monitoring) on panda habitat change (Model 1 in Table 4.2). To account for spatial autocorrelation among pixels, we also generated spatial simultaneous autoregressive lag and error models (Model 2 and 3, respectively, in Table 4.2). We defined a spatial weights matrix for the lag and error models by considering a pixel as a neighbor of another pixel if the Euclidean distance between them is shorter than the range in the variogram of the residuals from the OLS model. Finally, to further examine the effect of monitoring types under different payment levels 91 (i.e., higher payments in Wolong and Gengda and lower payments in Sanjiang), we correlated HSI changes to four different implementation methods (2 monitoring types × 2 payment levels; Model 4 in Table 4.2). We conducted all statistical analyses using R (R Development Core Team 2011) with the packages “car”, “gstat” and “spdep”. 92 Table 4.1. Independent variables included in regression models that relate the changes in panda habitat suitability index (HSI) with  the NFCP implementation and other biophysical and anthropogenic factors.  Variable Unit Description HSI_2001 unitless value of habitat suitability index (ranging from 0 to 1) in 2001 for each pixel (250 × 250m) of the HSI map (HSI pixel) FC_2001 % percentage of forested pixels (30 × 30m) of a binary forest cover map [derived from a 2001 Landsat TM image (Viña et al. 2007)] within the surrounding eight pixels of each HSI pixel Elevation m average of elevation over the pixels (90 × 90m) of a digital elevation model from the Shuttle Radar Topography Mission (DEM pixel; Farr et al. 2007) within each HSI pixel Roughness m standard deviation of elevation over the DEM pixels within each HSI pixel Aspect_north degree deviation from north (0 – 180°) Aspect_east degree deviation from east (0 – 180°) CTI 2 m /radian Compound Topographic Index, a function of both the slope and the upstream contributing area per unit width orthigonal to the flow direction (Moore et al. 1993) Dist2Household m Euclidean distance from each HSI pixel to the nearest household Dist2Road m the nearest Euclidean distance from each HSI pixel to paved roads Monitoring dummy NFCP monitoring type (household monitoring vs. government monitoring) Monitoring × Payment dummy two monitoring types × two payment levels (high payment vs. low payment) 93 Table 4.2. Summary of regression models that relate the changes in panda habitat suitability  index (HSI) with a set of independent variables.  Variable Standardized Coefficient (standard error) Model 1 – OLS HSI_2001 FC_2001 Elevation Roughness Aspect_north Aspect_east CTI Dist2Household Dist2Road *** Model 2 – Lag *** -0.0731 (0.0027) -0.0650 (0.0029) *** *** 0.0186 (0.0024) 0.0170 (0.0023) *** *** Model 3 – Error *** -0.0947 (0.0031) *** 0.0215 (0.0029) *** Model 4 – Error *** -0.0951 (0.0032) *** 0.0216 (0.0029) *** -0.0369 (0.0030) * -0.0049 (0.0023) -0.0374 (0.0045) -0.0416 (0.0046) * -0.0475 (0.0029) -0.0027 (0.0024) -0.0032 (0.0024) -0.0043 (0.0022) *** *** *** *** 0.0088 (0.0021) 0.0092 (0.0021) 0.0129 (0.0023) 0.0131 (0.0023) 0.0038 (0.0023) 0.0019 (0.0022) -0.0020 (0.0023) -0.0022 (0.0023) *** * -0.0076 (0.0023) * * -0.0052 (0.0022) -0.0019 (0.0029) ** -0.0076 (0.0029) * -0.0045 (0.0021) -0.0183 (0.0074) -0.0199 (0.0071) -0.0085 (0.0069) 0.0008 (0.0072) * -0.0067 (0.0026) -0.0065 (0.0025) 94 * -0.0053 (0.0021) ** Table 4.2. (Continued)    Variable Standardized Coefficient (standard error) Model 1 – OLS Model 2 – Lag Model 3 – Error Model 4 – Error Monitoring Government Household *** * 0.0184 (0.0029) 0.0065 (0.0029) *** *** 0.0297 (0.0038) 0.0424 (0.0039) ** 0.0184 (0.0065) *** 0.0411 (0.0074) Monitoring × Payment *** Government – High 0.0271 (0.0076) *** Household – High 0.0539 (0.0083) Government – Low -0.0023 (0.0107) Household – Low 0.0033 (0.0151) *** Auto-regressive Term 0.4418 (0.0289) Moran’s I of residuals 0.1737 AIC *** -4583 0.0544 -4770 95 *** *** *** 0.6776 (0.0251) 0.6616 (0.0259) 0.0131 -0.0121 -4995 -5000 Results From 2001 to 2007, the average HSI values increased ca. 6.7% (from 0.210 to 0.224) across the entire reserve, and increased ca. 7.1% (from 0.366 to 0.392) under the tree line. The 2 area of suitable habitat increased ca. 3.4% (from 686 to 709 km ) during the same period. Choosing different thresholds affected the magnitude, but not the trend, of panda habitat change since the habitat area increased ca. 2.4% and 5.9% with thresholds corresponding to 5% and 15% omission errors, respectively. However, considerable spatial heterogeneity of the HSI changes was observed across the reserve. For example, larger HSI increases occurred near human settlements in Wolong and Genda townships, while large HSI decreases occurred in the north of Sanjiang (Figure 4.1). To understand what factors affected the spatial heterogeneity, HSI changes were correlated with several biophysical and anthropogenic factors, as well as different forest monitoring types, in a pixel-wise OLS regression model (Model 1 in Table 4.2) with an adjusted 2 2 R equal to 0.26. A score test for non-constant variance (χ = 1.29, df = 1, p = 0.26) and an examination of variance inflation factors (< 3 for all independent variables) conducted using the R package “car” (Fox and Weisberg 2011) indicated no heteroscedasticity or multicollinearity problems. The significant Moran’s I for model residuals (Table 4.2) indicates the need to account for spatial autocorrelation. The autoregressive terms were significant in both lag and error models, but significant Moran’s I for the residuals of the lag model and the smaller Akaike Information Criterion value for the error model (Table 4.2) suggest that the error model (Model 3) better controls for the spatial autocorrelation. 96 The error model showed that several biophysical and anthropogenic factors are significantly correlated with HSI changes (Table 4.2). Better panda habitat improvement tended to occur in areas located in south-facing slopes, at lower elevations, with lower Compound Topographic Index values, surrounded by forests, and closer to local households (Table 4.2). As expected, there was a negative relationship between HSI changes and starting HSI values (i.e., in 2001) because pixels with higher values had less room to increase. While significantly positive coefficients for both monitoring types (Table 4.2) suggest an overall improvement of panda habitat in the reserve, the significant difference between the two coefficients indicates a better habitat recovery within household-monitored areas than within government-monitored areas (Figure 4.2a). However, the significant effect of household monitoring was only observed under the higher, but not the lower payment level (Model 4 in Table 2 and Figure 4.2b). Figure 4.2. Effects of the implementation of the Natural Forest Conservation Program (NFCP) on  giant panda habitat change in Wolong Nature Reserve, China.  Boxplots show the partial  residuals [controlling for biophysical and anthropogenic factors (Table 4.1)] for two different  NFCP monitoring types [government‐ (red) vs. household‐monitoring (blue)].  Plotted values are  the residuals of the full model plus the partial predicted values (black lines) in the entire study  97 area (a) and under two different NFCP payment levels (b).  Results of Student’s t‐tests on the  N.S. difference between monitoring types are shown (* p < 0.05;   p > 0.05).    Discussion To the best of our knowledge, this study is the first empirical and spatially explicit assessment of the effect of NFCP implementation on the habitat dynamics of an endangered species. Our results showed an overall improvement of panda habitat in Wolong Nature Reserve after the NFCP implementation. This reversed a trend of continuous habitat loss and degradation since the 1960s (Liu et al. 2001; Viña et al. 2007), and thus suggests a positive contribution of NFCP implementation to panda habitat recovery. As expected, several biophysical and anthropogenic factors affected the spatial patterns of panda habitat change. The positive effect of surrounding forest cover on panda habitat recovery not only reflects the important roles of forest cover in affecting bamboo growth and distribution and in constituting panda habitat (Schaller et al. 1985), but also suggests that forest regeneration could be an important underlying process of the habitat recovery. Higher vegetation regeneration rates associated with higher temperature and solar radiation may explain the larger HSI increases observed on south-facing slopes at lower elevations. Since the Compound Topographic Index is highly correlated with several physical and chemical attributes of soils (Moore et al. 1993), its significant relationship with HSI changes suggests the influence of soil attributes on vegetation regeneration and in turn on panda habitat recovery. Our regression models also indicated better habitat recovery in areas closer to local households, where higher intensity of human activities is expected (He et al. 2009). This suggests that human 98 activities involving resource exploitation may be diminishing after the NFCP implementation, and are thus inducing a net beneficial effect on panda habitat. Comparing panda habitat suitability before vs. after NFCP implementation alone cannot totally rule out the potential effects of other temporally variable factors (e.g., changes in socioeconomic conditions of local households and other conservation instruments implemented simultaneously with the NFCP). However, by correlating spatial heterogeneity of the observed habitat dynamics with NFCP implementation, our regression models indicated a significantly positive effect of NFCP monitoring, especially household monitoring, on panda habitat. This suggests that NFCP implementation, particularly with an active engagement of local residents, is at least one of the major reasons for the panda habitat recovery. In addition, some other conservation policies may not contribute as much as NFCP to panda habitat recovery. For example, the Grain-to-Green Program (GTGP), which provides local farmers with cash or grain subsidies to encourage a conversion of cropland on steep slopes into forest or grassland (Liu et al. 2008; Yin and Yin 2010), has also been implemented in the reserve since 2000. However, its direct contribution to panda habitat recovery is negligible because GTGP-enrolled croplands only cover a small portion of the reserve [ca. 367 ha (Wolong Nature Reserve 2005) and < 0.2% of the entire reserve], and the tree seedlings and saplings planted cannot provide suitable panda habitat within the 7-year timeframe of this assessment (Bearer et al. 2008). Finally, since human population increased ca. 6% and the number of households increased ca. 23% in the three townships between 2001 and 2007 (Wenchuan Statistics Bureau 2008), an expected increase in resource consumption during the same period (Liu et al. 2003a) cannot explain the observed panda habitat recovery. 99 While further studies are warranted to understand the mechanisms behind the effectiveness of the innovative implementation of NFCP, two reasons may explain its conservation effectiveness over the previous “fences-and-fines” strategies. First, the direct payments to local residents compensate the costs of forgoing resource-depleting activities and thus may create stronger conservation incentives (Engel et al. 2008; Ferraro and Kiss 2002). The association between the stronger effect of household monitoring on panda habitat recovery and the higher payment level indicates the important role of the payments in determining NFCP’s effectiveness, and implies that NFCP benefits panda conservation when the associated payments are high enough to create conservation incentives (Engel et al. 2008). Besides creating incentives, the payments may also encourage electricity purchases and thus reduce the dependency of local residents on fuelwood as an energy source (An et al. 2002). Second, the shared responsibility for forest monitoring and sanctions among households in a monitoring group may enhance rule compliance through social norms (Chen et al. 2009; Dietz et al. 2003). Thus, people may fulfill their monitoring duties to avoid payment reductions that could harm their social relations with other members of the same monitoring group. People may also avoid causing damage on the parcels monitored by other groups if they do not want to harm their social relations with the people in those groups. Although biodiversity conservation is not a main aim of NFCP, this study empirically shows that NFCP is effective in conserving the habitat of an endangered species, in addition to protecting forests and soils (Liu et al. 2008), and even better effectiveness can be achieved by engaging local residents in forest monitoring with direct payments. Since in most of China local residents do not directly participate in NFCP implementation (Yin and Yin 2010), a greater overall benefit to ecosystems and their services could be achieved if the successful experience of 100 incentive-based and local-engaged NFCP implementation can be spread to other places. However, while our study suggests that this instrument is more effective than reserve establishment, higher conservation effectiveness and efficiency may be achieved by implementing both strategies. For instance, the costs of monitoring local people’s compliance and conservation outcomes under PES programs may be reduced if the resident monitoring is combined with regular patrols by reserve officials. The regulations applied in protected areas may also reduce the expected benefits from non-compliance and thus increase incentives to participate in PES programs (Engel et al. 2008). Therefore, the incentive-based and localengaged instrument should not be used as a substitution for, but a complement to, reserve establishment. Complementary instruments that pursue resource conservation but explicitly incorporate human needs (e.g., the combination of PES, decentralized management, logging ban and reserve regulations) may offer a potential solution to the “park vs. people” debate.  101 CHAPTER 5 CASCADING EFFECTS OF CLIMATE CHANGE THROUGH TROPHIC RELATIONS: PROJECTED FOOD SHORTAGES THREATEN THE GIANT PANDA In collaboration with Andrés Viña, Julie A. Winkler, Yu Li, Weihua Xu, Zhiyun Ouyang, and Jianguo Liu 102 Abstract Climate change is one of the most serious threats to global ecosystems because of its significant impacts not only on the survival of individual species, but also on their ecological functions. Understory plant species play important roles in forest ecosystems by regulating forest succession and structure, facilitating nutrient and energy cycling, and supplying shelter and food resources for many wildlife species. However, studies on the effects of climate change on understory plants, particularly on their role as a food resource for wildlife species, are scarce in the literature. Here we report the first quantitative and comprehensive climate change impact assessment for understory bamboo species and their trophically dependent giant pandas (Ailuropoda melanoleuca). An ensemble of projected changes in bamboo distribution associated with multiple climate change projections and bamboo dispersal scenarios indicates a substantial reduction in the distributional ranges of three dominant bamboo species in the Qinling Mountains, China during the 21st century. As these three species comprise almost the entire diet of the panda population in the region, the projected changes in bamboo distribution suggest a potentially dangerous shortage of food for this population. This study underscores the importance of incorporating inter-specific cascading effects of climate change into impact assessments and associated conservation planning. Introduction Global climate change is one of the most serious threats that human society faces in the 21st century. Not only human well-being, but also ecosystems and biodiversity have been and will be further influenced by climate change (IPCC 2007). Observations around the world have shown that the altered climate is impacting natural systems in many ways and the impacts are 103 across biological taxa, ecosystem types, and geographic regions (Parmesan 2006; Parmesan and Yohe 2003; Root et al. 2003). Among other impacts, climate change-induced shifts in species’ distributional ranges have been observed (Chen et al. 2011; Parmesan 2006; Root et al. 2003) and been projected (Araújo et al. 2006; Bakkenes et al. 2002; La Sorte and Jetz 2010) for many species. Changes in species distributional ranges may not only cause extinction of individual species and populations, but also result in changes in the structure of biological communities and properties and function of ecosystems (Parmesan 2006). Thus they pose a particular challenge to biodiversity conservation (Coetzee et al. 2009; Hannah et al. 2002; Heller and Zavaleta 2009). Understory vegetation plays an important role in forest ecosystems by regulating forest structure, maintaining ecological functions, and sustaining biodiversity (Gilliam 2007; Hagar 2007). Like many understory plants, understory bamboo species are also an essential component in many forest ecosystems (Griscom and Ashton 2003; Taylor et al. 2004). They not only influence species composition and structural complexity of forests (Griscom and Ashton 2003; Taylor et al. 2004), but also provide essential food resources for diverse wildlife species, including one of the most endangered species in the world, the giant panda (Pan et al. 2001; Schaller et al. 1985). While deforestation and forest degradation are threatening the survival of about half of all bamboo species worldwide (including many understory species) (Bystriakova and Kapos 2006), climate change may present an additional significant threat. Bamboo species are particularly vulnerable to climate change because their unusual extended sexual reproduction intervals (from 10 to 120 years) (Janzen 1976) render them less capable of adapting to the rapidly changing climate projected to occur within this century (IPCC 2007). However, knowledge of climate change-induced dynamics of bamboo distribution and cascading effects on 104 other species is limited. This leaves an important knowledge gap of particular significance for biodiversity conservation efforts. Understory bamboo is one of the most important components of giant panda habitat because pandas are extreme dietary specialists who devote most of their active time feeding on large amounts of bamboo (up to 38 kg/day) to compensate their poor ability to digest cellulose (Pan et al. 2001; Schaller et al. 1985). The dependency of giant pandas on understory bamboo not only affects current panda distribution (Bearer et al. 2008; Liu et al. 2005b; Reid et al. 1989), but is also believed to have driven historical changes in panda distribution as bamboo distributions shifted in response to climate fluctuations (Pan et al. 2001; Schaller et al. 1985). Therefore, the close relationship between bamboo species and giant pandas provides an excellent opportunity to address the knowledge gap on the impacts of climate change on understory plants and their cascading effects on trophically dependent wildlife species. The main goal of this study was to quantitatively assess the potential impacts of climate change on bamboo distribution, with an emphasis on the cascading effects on giant panda habitat. With a focus on three bamboo species dominating the forest understory in the Qinling Mountains region (i.e., Qinling arrow, dragon-head and wooden), this study investigated the potential changes in climatically suitable areas for the three bamboo species under projected climate change, and evaluated the cascading effects of these changes on the spatiotemporal dynamics of the amount and distribution of giant panda habitat. The assessment may not only provide essential information for giant panda conservation in the face of climate change, but also address the knowledge gap on the cascading effects of climate change through inter-specific interactions. 105 Materials and Methods Field data 2 To collect bamboo presence data, we established 293 field plots (ca. 300 m ) across the Qinling Mountains and geo-referenced the center of each plot using a Global Positioning System (GPS) receiver in 2005, 2007 and 2008. Field plots were established in diverse types of forest (including broadleaf deciduous, coniferous and mixed primary and secondary forests, as well as in tree plantations) and within an elevation range from 1,000 to 3,000 m, which roughly covers the elevational range of the distributions of the three bamboo species studied and the giant pandas in this region (Pan et al. 2001; State Forestry Administration 2006). We recorded the presence of Qinling arrow, dragon-head, and wooden bamboos in 86, 40, and 79 plots, respectively. Baseline climate conditions We obtained gridded climate data from the WorldClim database (http://www.worldclim.org/; Hijmans et al. 2005) as a representation of the baseline climate conditions. The WorldClim dataset was generated using a thin plate spline scheme that considers longitude, latitude and elevation [obtained from the Shuttle Radar Topography Mission (SRTM; http://www2.jpl.nasa.gov/srtm/)] to interpolate observations of monthly mean, maximum, and minimum temperature and monthly total precipitation for 1950 – 2000 from weather stations worldwide to a 30 arc-second resolution grid (Hijmans et al. 2005). The dataset additionally includes 19 biologically meaningful (i.e., bioclimatic) variables (e.g., mean diurnal temperature, mean temperature in the warmest month and precipitation seasonality), derived from the monthly temperature and precipitation gridded fields. 106 Future climate projections We obtained gridded projections of bioclimatic variables from the database of the International Center for Tropical Agriculture (CIAT; http://ccafsclimate.org/download_sres.html). The projected values are derived from monthly temperature and precipitation gridded fields downscaled from several general circulation models (GCMs; Supplementary Table 5.1) for the third (TAR) and fourth assessment reports (AR4) of the Intergovernmental Panel on Climate Change (IPCC 2001, 2007) under SRES A2 and B2 greenhouse gas emissions scenarios. Differences in monthly temperature fields and ratios of the precipitation fields between GCM simulations for a reference periods (1961 – 1990) and a future time slice (e.g., 2040 – 2069) were calculated for each GCM grid point. The differences (or ratios) were interpolated by using the thin plate spline approach to the 30 arc-second grid used by the WorldClim and then were added to (or multiplied by) the climatological fields in the WorldClim database. Bioclimatic models We used the maximum entropy modeling (Maxent; Phillips et al. 2006) to generate bioclimatic models for characterizing and mapping climatically suitable areas (CSAs) for each bamboo species. Maxent is a machine-learning and niche-based approach for establishing the relationship between geographic locations of species occurrence and corresponding environmental conditions, and for mapping the spatial distribution of suitable habitat for the species based on the species-environment relationship (Phillips et al. 2006; Phillips and Dudík 2008). This modeling approach estimates the probability of species presence (or habitat suitability) across space, given the spatially continuous values of environmental variables, by 107 finding the most evenly distributed probabilities (i.e., with the maximum entropy) as constrained by the distribution of values of environmental variables associated with species presence locations. Because it is practically impossible to confirm absence of bamboo species within an 2 area of ca. 1 km (the grid size of the climate data and the spatial unit of our bioclimatic models), the requirement of presence-only data makes maximum entropy modeling one of the most suitable approaches for this study. To generate bioclimatic models, we selected seven bioclimatic variables (Table 5.2) that were: (1) statistically important for fitting the bamboo presence data, (2) biologically important for the three bamboo species, and (3) less correlated with each other. For this, we selected statistically important variables based on the results of a jackknife analysis. In this analysis, we compared the accuracy of the model built using all 19 variables with models built with each individual variable, and with models built with all but each one of the variables. Higher accuracy of a single-variable model indicates that the variable in question contains more useful information for modeling bamboo distribution, while a larger reduction of accuracy of a leaveone-out model indicates that the omitted variable contains more unique information. Using the values of 50,000 randomly selected pixels (~14% of total pixels in the study area), we calculated pair-wise Pearson’s correlation coefficients to identify highly correlated variables. For the variables that were statistically important but highly correlated with others, our choice of which variables to retain was based on observational and experimental studies on the influence of climate conditions on bamboo species (Li 1997; Qin et al. 1993). 108 Table 5.1.  General circulation models (GCMs) from which downscaled future climate projections were obtained.  GCMs from the  IPCC Third Assessment were used for evaluating temporal dynamics of bamboo distribution and panda habitat over the 21st century.   The other models are from the IPCC Fourth Assessment and were used to better capture uncertainty introduced by different GCMs.  GCM Center SRES Scenarios Bjerknes Centre for Climate Research, Norway CCSR/NIES Center for Climate System Research, Japan IPCC # Assessment A2 2050s AR4 A2 & B2 BCCR-BCM2 Time Periods of † Downscaled Data 2020s, 2050s & 2080s TAR National Institute for Environmental Studies, Japan CGCM2 Canadian Center for Climate Modelling and Analysis, Canada A2 & B2 2020s, 2050s & 2080s TAR CGCM3.1 (T47) Canadian Center for Climate Modelling and Analysis, Canada A2 2050s AR4 CGCM3.1 (T63) Canadian Center for Climate Modelling and Analysis, Canada A2 2050s AR4 CNRM-CM3 Centre National de Recherches Meteorologiques, France A2 2050s AR4 CSIRO-Mk2 Commonwealth Scientific and Industrial Research Organisation, Australia A2 & B2 2020s, 2050s & 2080s TAR CSIRO-Mk3.0 Commonwealth Scientific and Industrial Research Organisation, Australia A2 2050s AR4 109 Table 5.1.  (Continued) GCM ECHO-G Center SRES Scenarios Time Periods of † Downscaled Data IPCC # Assessment A2 2050s AR4 Meteorological Institute, University of Bonn, Germany Meteorological Research Institute of KMA, Korea Model and Data Group at MPI-M, Germany FGOALS-g1.0 Institute of Atmospheric Physics, China A2 2050s AR4 GFDL-CM2.0 Geophysical Fluid Dynamics Laboratory, USA A2 2050s AR4 GFDL-CM2.1 Geophysical Fluid Dynamics Laboratory, USA A2 2050s AR4 GISS-AOM Goddard Institute for Space Studies, USA A2 2050s AR4 HadCM3 Hadley Centre for Climate Prediction and Research, UK A2 & B2 2020s, 2050s & 2080s TAR IPSL-CM4 Institut Pierre Simon Laplace, France A2 2050s AR4 MIROC3.2 (hires) National Institute for Environmental Studies, Japan A2 2050s AR4 MIROC3.2 (medres) National Institute for Environmental Studies, Japan A2 2050s AR4 MRI-CGCM2.3.2 Meteorological Research Institute, Japan A2 2050s AR4 PCM National Centre for Atmospheric Research, USA A2 2050s AR4 † # 2020s: 2010 – 2039; 2050s: 2040 – 2069; 2080s: 2070 – 2099 TAR: IPCC Third Assessment Report; AR4: IPCC Fourth Assessment Report 110 Table 5.2.  Bioclimatic variables used in the final bioclimatic models for the three bamboo species studied.  WorldClim Code Bio2 Name Mean Diurnal Range Description Average of the differences between monthly maximum and minimum temperatures Bio4 Temperature Seasonality Standard deviation of monthly mean temperatures Bio10 Mean Temperature of Warmest Quarter Average of the mean temperatures in the warmest three months Bio11 Mean Temperature of Coldest Quarter Average of the mean temperatures in the coldest three months Bio15 Precipitation Seasonality Coefficient of variation of monthly total precipitations × 100 Bio18 Precipitation of Warmest Quarter Total precipitation in the warmest three months Bio19 Precipitation of Coldest Quarter Total precipitation in the coldest three months 111 We randomly selected 70% of the presence data of each bamboo species for model development and used the remaining 30% for validation. To account for the uncertainty in the model structure introduced by potential biases in data partitioning, we generated 10 bioclimatic models for each species using 10 different data partitions. The variation among the results of the 10 models may also partially reflect the uncertainty associated with potentially incomplete and biased observations of bamboo presence in the field. We used the receiver operating characteristic (ROC) curve analysis (Fielding and Bell 1997) for model validation and calculated the area under the ROC curve (AUC) to measure model accuracy. While an AUC value of 1 indicates a perfect prediction, a value of 0.5 indicates a random prediction (Hanley and Mcneil 1982). As suggested for the application on presence-only models (Phillips et al. 2006), we calculated AUC values by contrasting presence pixels against those randomly selected from the study area (i.e., background pixels). We selected 10,000 background pixels within the elevational range of the distributions of the three bamboo species (i.e., 900 – 3,000 m) (Pan et al. 2001) for the ROC analysis. We then calculated average AUC values over the 10 models for each bamboo species. Climatically suitable areas for bamboo species Using the bioclimatic models and the bioclimatic variables under the baseline climate, we mapped the spatial distributions of baseline climatically suitable areas (CSAs) across the study area for each bamboo species. For the output of each bioclimatic model, we used a threshold corresponding to a 10% omission error for converting the continuous presence probabilities into binary values (i.e., suitable vs. non-suitable). We determined thresholds based on omission errors only because we were limited to presence-only validation data. Therefore, without 112 controlling for commission errors, the output binary maps may overestimate the extent of CSAs for the three bamboo species. Using the bioclimatic models and the bioclimatic variables under the projected future climate, we then projected future distributions of CSAs for the three bamboo species. We examined temporal dynamics of the CSAs over the 21th century using the climate projections from four IPCC TAR GCMs (Table 5.1) for the SRES A2 and B2 greenhouse gas emissions scenarios because fine-resolution downscaled climate projections are readily available for three time slices (2010 – 2039, 2040 – 2069 and 2070 – 2099) for these models. To better capture the uncertainty introduced by different GCMs, we also incorporated projections from 15 IPCC AR4 GCMs (Table 5.1) that had been downscaled to a fine resolution using the same methods as applied to the IPCC TAR GCMs, although available only for the time slice of 2040 – 2069 under the SRES A2 scenario. Extent of giant panda habitat We defined giant panda habitat areas as those projected to be occupied by any of the three bamboo species studied. We generated occupation maps for each of the bamboo species based on the CSA maps and on two bamboo dispersal scenarios: unlimited dispersal and no dispersal. Under the unlimited dispersal scenario, all pixels within the future CSAs projected by the bioclimatic models were considered to be occupied. Under the no dispersal scenario, only the projected future CSAs of a bamboo species that overlap with its baseline CSAs were considered to be occupied. To calculate the area of panda habitat, we counted the number of pixels that were projected to be occupied by at least one bamboo species for every possible combination of the 113 projections of occupied CSAs for each bamboo species (10×10×10 bioclimatic models × 2 dispersal scenarios). This calculation is based on the assumption that the occurrence of at least one of the three bamboo species is sufficient for providing food for the pandas. Since the presence of bamboo is a necessary but not sufficient habitat condition for the panda, this assumption may lead to an overestimation of habitat amount. We calculated the habitat area under the baseline climate and under each future climate projection for the different time slices [4 (or 15) GCMs × 2 (or 1) greenhouse gas emissions scenarios × 3 (or 1) time slices]. We then reported the projected change in the habitat area as a percentage relative to the baseline. Results Using bioclimatic models to associate bamboo presence locations recorded in the field with bioclimatic variables, we characterized and mapped current suitable climate conditions for each of the three bamboo species (Figure 5.1). The predicted spatial patterns of the climatically suitable areas (CSAs) effectively captured the observed bamboo presence locations according to model accuracy evaluations. The average AUC values (0.98, 0.96 and 0.97 for Qinling arrow, dragon-head and wooden bamboos, respectively) indicated good model performance on capturing the observed presence locations of the three bamboo species. The spatial patterns were also consistent with our understanding of the current distribution of understory bamboo (i.e., Qinling arrow and wooden bamboos are distributed from 1,800 m to 3,000 m, and from 900 m to 1,900 m, respectively, with the distribution of dragon-head bamboo overlapping between 1,100 and 2,300 m; Figure 5.2) (Pan et al. 2001; State Forestry Administration 2006). 114 Figure 5.1. Location of the study area relative to current giant panda distributional range and  the baseline climatically suitable areas (CSAs) for the three bamboo species studied.  Blue, red,  and green colors indicate the CSAs for Qinling arrow, dragon‐head, and wooden bamboos,  respectively, based on the bioclimatic models under the baseline climate.  The mixtures of the  three colors indicate overlaps of the CSAs for individual species.  The brightness of colors shows  the number of bioclimatic models (among the 10 models with different presence data  partitions) predicting pixels as suitable for each species, with brighter colors indicating a larger  number of models.    115 20102039 20402069 20702099 20102039 20402069 20702099 20102039 20402069 20702099 Figure 5.2. Projected elevational distributions of the climatically suitable areas (CSAs) for the  three bamboo species studied.  The elevation values were obtained for pixels which were  projected to be climatically suitable for the wooden (a and d), dragon‐head (b and e), and  Qinling arrow bamboo (c and f) by at least one of 10 bioclimatic models with future climate  projections for three time slices in the 21th century from four IPCC TAR GCMs (Table 5.1) under  the SRES A2 (a – c) and B2 (d – f) greenhouse gas emissions scenarios.  Each box plot shows the  maximum, 75th percentile, median, 25th percentile, and minimum values.  The grey zone on  the background shows the inter‐quartile range (75th percentile – 25th percentile) of elevation  for the baseline CSAs (the CSAs under the baseline climate).  The white line within the grey zone  indicates median values and black dashed lines indicate the maximum and minimum values of  elevation for the baseline CSAs.  116 The projections of our bioclimatic models under the climate projections from the four IPCC TAR GCMs show substantial changes in the CSAs during the 21st century with considerable differences among GCMs (Figure 5.3). Areas in the north and northwest of the Qinling Mountains become climatically suitable for the three bamboo species under some climate projections, especially for the lower-elevation species, i.e., the wooden bamboo (Figure 5.3). However, within the Qinling Mountains region, our projections suggest that the CSAs shift to, and become restricted at, higher elevations (Figure 5.2). These projected northward and upward shifts are consistent with observations and simulations for many other species worldwide (Chen et al. 2011; IPCC 2007; Parmesan 2006; Skov and Svenning 2004). The projected CSAs under the climate projections from the 15 IPCC AR4 GCMs also indicate large GCM-introduced uncertainty about the extent and spatial distribution of CSAs of the three bamboo species (Figure 5.4). However, compared to the CSA projections associated with the TAR GCMs (Figure 5.3), those associated with the newer GCMs show a higher degree of consensus since all of them indicate a reduction of CSA extent (Figure 5.4). 117 2010-2039 2040-2069 2070-2099 Figure 5.3. Projected future distributions of climatically suitable areas (CSAs) for the three  bamboo species studied under the climate projections from four IPCC TAR GCMs (Table 5.1) for  three time slices under SRES A2 and B2 greenhouse gas emissions scenarios.  Blue, red, and  green colors indicate the CSAs for the Qinling arrow, dragon‐head, and wooden bamboos,  respectively.  See the legend of Figure 5.1 for the detailed information on the color  representation.  118 2010-2039 2040-2069 2070-2099   Figure 5.3. (Continued)  The substantial changes in CSAs for bamboo species, in turn, result in drastic panda habitat changes. Our assessment suggests that almost the entire panda habitat in the region disappears by the end of the 21st century (80 – 100% projected decreases in median area) if bamboo species are not able to colonize new CSAs beyond their current distributional ranges (Figure 5.5b and d). With unlimited bamboo dispersal ability, a considerable amount of panda habitat is projected to persist over the entire century, but only if future climate is closer to the projections from two of the four GCMs (CSIRO-Mk2 and HadCM3), and then only with a lower greenhouse gas emissions scenario (SRES B2; Figure 5.5a and c). All panda habitat projections based on the downscaled projections from the 15 IPCC AR4 GCMs show substantial habitat losses (59 – 100% projected decreases in median area), with about half of the projections 119 indicating an almost complete loss, by the middle of this century under both bamboo dispersal scenarios (Figure 5.6). Compared with the habitat projections derived from the IPCC TAR GCMs, those associated with the AR4 GCMs, including those obtained from newer versions of TAR GCMs (CSIRO-Mk3.0 vs. CSIRO-Mk2 and CGCM3.1 vs. CGCM2), indicate more drastic losses of panda habitat (Figure 5.6). This comparison suggests that the projected habitat expansion associated with two of the older GCMs (CSIRO-Mk2 and HadCM3; Figure 5.5) may be overly optimistic (Figure 5.6a). 120 Figure 5.4. Projected future distributions of climatically suitable areas (CSAs) for the three  bamboo species studied under the climate projections from 15 IPCC AR4 GCMs.  The climate  projections were from the GCMs (Table 5.1): (a) BCCR‐BCM2, (b) CGCM3.1 (T47), (c) CGCM3.1  (T63), (d) CNRM‐CM3, (e) CSIRO‐Mk3.0, (f) ECHO‐G, (g) FGOALS‐g1.0, (h) GFDL‐CM2.0, (i) GFDL‐ CM2.1, (j) GISS‐AOM, (k) IPSL‐CM4, (l) MIROC3.2 (hires), (m) MIROC3.2 (medres), (n) MRI‐ CGCM2.3.2, and (o) PCM for the time slice of 2040 – 2069 under the SRES A2 greenhouse gas  emissions scenarios.  See the legend of Figure 5.1 for the detailed information on the color  representation.      121 Figure 5.5. Temporal dynamics of the projected changes in the area of giant panda habitat over  the 21st century.  The changes (%) are relative to the area of panda habitat (climatically suitable  areas for the three bamboo studied combined) under the baseline climate.  The projected  values were obtained from 1,000 combinations of the bioclimatic models for the three bamboo  species (10 models for each species) under the unlimited (a and c) and no (b and d) bamboo  dispersal scenarios and under multiple future climate projections for three time slices from four  IPCC TAR GCMs (Supplementary Table S1) under the SRES A2 (a and b) and B2 (c and d)  greenhouse gas emissions scenarios.  The points indicate median values of the projections from  the 1,000 projections, and the vertical bars indicate the range of projected values.   122 NGMH FOLB RCJ K DE I QPA S GNMF LHOBC RK J DI E QPA S Figure 5.6. GCM‐related uncertainty of projected changes in giant panda habitat area for the time slice of 2040 – 2069 under the  SRES A2 greenhouse gas emissions scenario.  The percent changes are relative to the area of panda habitat (climatically suitable  areas for the three bamboo studied combined) under the baseline climate.  The projected values were obtained from 1,000  combinations of the bioclimatic models for the three bamboo species (10 models for each species) under the unlimited (a) and no (b)  bamboo dispersal scenarios and under multiple future climate projections from four IPCC TAR and 15 AR4 GCMs (Table 5.1).  Each  boxplot shows the maximum, 75th percentile, median, 25th percentile and minimum of the projected values.  A: BCCR‐BCM2; B:  CCSR/NIES; C: CGCM2; D: CGCM3.1(T47); E: CGCM3.1(T63); F: CNRM‐CM3; G: CSIRO‐Mk2; H: CSIRO‐MK3.0; I: ECHO‐G; J: FGOALS‐ g1.0; K: GFDL‐CM2.0; L‐GFDL‐cm2.1; M: GISS‐AOM; N: HadCM3; O: IPSL‐CM4; P: MIROC3.2(hires); Q: MIROC3.2(medres); R: MRI‐ CGCM2.3.2; S: PCM. 123 Discussion While repeated calls have been made to quantify the potential impacts of climate change on China’s biodiversity (Liu and Raven 2010; Ministry of Environmental Protection of China 2010), this study provides the first quantitative and spatially explicit climate change impact assessment for understory bamboo and the trophically dependent giant panda. Our ensemble of model projections indicated potential challenges for the conservation of the three bamboo species and thus the giant panda. In the worst situations (in terms of the projected extent of bamboo distributional ranges and panda habitat) shown in our assessment, almost the entire region will become climatically unsuitable for the three bamboo species by the middle of this century (Figure 5.6). As these species currently constitute almost the entire diet of giant pandas in this region (Pan et al. 2001), the pandas may face a shortage of food, unless they can find alternative food resources, an unlikely outcome because of the panda’s specific diet needs. Although giant pandas have survived large area die-offs of single bamboo species by shifting their home ranges and foraging on non-affected bamboo species (Pan et al. 2001; Reid et al. 1989), they may face regional extinction if climate change, as projected by this study, induces simultaneous die-offs of multiple bamboo species. Other bamboo species of the region, especially those currently growing at lower elevations (e.g., Phyllostachys sulphurea), might have the potential to occupy the areas currently dominated by the three species evaluated. A concern is that low-elevation species may not be able to meet pandas’ needs mainly because of their thicker culms, which may reduce pandas’ forage efficiency (Tarou et al. 2005). Even under the most optimistic projections, the bamboos and pandas may still face substantial challenges due to climate change. The climatically suitable areas (CSAs) for the three bamboo species, although still relatively large in extent, are distant from current bamboo 124 distributional ranges and panda habitat (Figure 5.3u and x), and lie outside of current protected areas (Figure 5.7). With a long history of intense human activities surrounding the Qinling Mountains (e.g., Xi’an, at the northern foot of the Mountains, is a city with more than eight million people), keeping the projected CSAs from human disturbance, and maintaining or reestablishing ecological connectivity among those areas will be difficult. This fragmented landscape may thus hinder the range shifts of bamboos and pandas in response to climate change. Figure 5.7. Percentage of projected giant panda habitat within current nature reserves in the  Qinling Mountains over the 21st century.  The extent of panda habitat was obtained by  combining the occupied portions of the projected future climatically suitable areas for the  Qinling arrow, dragon‐head, and wooden bamboo, under the unlimited bamboo dispersal  scenario.  The percentages are associated with the panda habitat projections from 1,000  combinations of the bioclimatic models for the three bamboo species (10 models for each  species) under future climate projections from four IPCC TAR GCMs (Table 5.1) under the SRES  A2 (a) and B2 (b) greenhouse gas emissions scenarios.  The points indicate the median values of  125 the 1,000 projections and the vertical bars indicate the maximum and minimum values.  The  horizontal solid lines indicate the median values and dashed lines indicate the maximum and  minimum values for the percentages under the baseline climate.    As a global icon of biodiversity conservation, the giant panda has attracted unparallel conservation efforts and resources to reduce the loss and degradation of its habitat due to land use/cover change (Liu et al. 2001). Besides more than 60 nature reserves established specifically for panda conservation, recent implementation of national conservation programs (e.g., Natural Forest Conservation Program and Grain-to-Green Program) has halted and even reversed a trend of deforestation not only within the panda distributional range (Viña et al. 2011) but in all of China (Liu et al. 2008). While it is believed that these programs are beneficial to panda conservation (Chapter 4) because deforestation has been a major threat to both the panda (Liu et al. 2001; Pan et al. 2001) and bamboo species (Bystriakova and Kapos 2006) for decades, our assessment suggests that these benefits may be potentially compromised by climate change. Therefore, current conservation efforts may lose their effectiveness if they do not explicitly address potential impacts of climate change. Besides proactive and adaptive conservation strategies for protecting current and projected habitat, and maintaining habitat connectivity (Heller and Zavaleta 2009), more aggressive actions, such as assisted colonization (HoeghGuldberg et al. 2008) for bamboo and the panda, also need to be considered. The success of captive breeding and ongoing re-introduction projects for giant pandas provide a good foundation for their assisted colonization. However, there is an urgent need to better understand the risks and challenges associated with assisted colonization (Hoegh-Guldberg et al. 2008), and to increase their effectiveness as conservation tools in the face of climate change. 126 Like other assessments of potential climate change impacts on species distributions using bioclimatic models, the present study is based on many assumptions, which constitute sources of uncertainty in model projections (Araújo and Guisan 2006; Guisan and Thuiller 2005). However, due to the particularities of this study, some assumptions are quite reasonable. For example, bamboo species’ evolutionary adaptations to novel climatic regimes may occur at a much slower pace than the predicted climate changes because of their unusual extended sexual reproduction intervals [from several decades to more than one century (Janzen 1976)], making niche conservatism (Wiens and Graham 2005) a reasonable assumption for the species studied. In addition, uncertainties due to other assumptions, including those associated with future climate projections and bamboo dispersal ability, have been captured, at least partially, by the projection ensemble. However, like all other climate change impact assessments, our projection ensemble cannot capture all potential uncertainties and may represent a narrower range of the full uncertainty (Winkler et al. 2011). For example, our bioclimatic models did not account for many biotic and abiotic factors other than climate (e.g., soil types, canopy cover, interactions with other plant species) which also limit current bamboo distributions. Omitting these factors may lead to an underestimation of climatically suitable areas (CSAs) because the bioclimatic models may not be able to capture the areas where other limiting factors cause the absence of the bamboo species even under suitable climate. In contrast, the extent of projected future CSAs shown in this study may be overestimated without considering the influence of above-mentioned limiting factors on bamboos’ colonization. However, some of these limiting factors (e.g., soil properties) may be climate dependent (Franzluebbers et al. 2001; McKenzie and Ryan 1999), 127 and thus part of their effects may be implicitly considered by the model through climate variability. Our definition of giant panda habitat in this study is based on two assumptions which are not captured by our projection ensemble and tend to cause an overestimation of the extent of panda habitat. First, we assume that giant panda habitat is determined only by the presence of bamboo species. Because giant pandas almost exclusively depend on bamboo species as food resources (Pan et al. 2001), bamboo availability is a necessary component of panda habitat. However, without the consideration of other factors [e.g., forest cover, human activities, and pandas’ climatic tolerance and dispersal ability (Liu et al. 2001; Pan et al. 2001)] which further limit panda distribution, this assumption may cause an overestimation of future giant panda habitat or an underestimation of panda habitat loss due to climate change. The second assumption is that the presence of any one of the three bamboo species evaluated is sufficient to support the giant panda’s demand for food. Therefore, without considering the spatial coverage and the amount of biomass of each bamboo species, our assessment tends to overestimate the area of panda habitat. Furthermore, because the giant pandas in the Qinling Mountains currently forage different bamboo species in different seasons (Pan et al. 2001), a single bamboo species may not be able to sufficiently support giant pandas year round. This may also cause an overestimation of projected panda habitat. However, because the same models and assumptions were applied to predicting the baseline habitat and projecting the future habitat, the percent changes relative to the baseline, which we reported in this study, can cancel out some of the biases caused by the bioclimatic models and associated assumptions, as well as those caused by the future climate projections. 128 The implications of this study are broad and well beyond the specific case presented here. While we focused on understory bamboo as food for giant pandas, providing food resources for wildlife is just one of many ecological functions performed by understory bamboos (Griscom and Ashton 2003; Taylor et al. 2004) and by many other understory plants (Gilliam 2007; Nilsson and Wardle 2005). Although several studies have observed and/or projected that climate change drives many species to shift their distributional ranges or go extinct in diverse ecosystems (Chen et al. 2011; IPCC 2007; Parmesan 2006; Skov and Svenning 2004), the cascading effects of these changes through inter-specific interactions, especially those across trophic levels (Van der Putten et al. 2010), are usually neglected. This may result in an underestimation of the impacts of climate change on ecosystems and seriously compromise the benefits obtained from current biodiversity conservation efforts. This study underscores the importance of incorporating these cascading effects not only in impact assessments but also in conservation planning in the face of climate change. 129 CHAPTER 6 SYNTHESIS AND CONCLUSIONS 130 Summary As a global icon of biodiversity conservation and one of the most endangered and most beloved species around the world, the giant panda has attracted substantial conservation resources. However, many previous efforts for conserving wild panda populations have had limited effectiveness and efficiency, partly due to lack of detailed information on the spatial patterns of panda habitat and their temporal dynamics under a changing environment. To address these information gaps, in this dissertation I developed effective and practical approaches for obtaining detailed information on the distribution of pandas’ staple food (i.e., understory bamboo) across large geographic regions and on the dynamics of panda habitat over time. I also investigated the effects of current conservation practices on the short-term panda habitat changes and the potential impacts of climate change on long-term habitat dynamics. These approaches and information substantially contribute to giant panda conservation. Using two dominant bamboo species in Wolong Nature Reserve, I showed that a combination of species distribution modeling and land surface phenology obtained from high temporal resolution remotely sensed data is a promising approach for providing detailed information on understory bamboo distribution across large geographic regions. This approach has several improvements on mapping understory vegetation as compared to other methods. Due to high availability of MODIS data in terms of areal cover and temporal resolution, this approach solves the problem of data limitation faced by many other methods. Thus, it may be easily applied not only for mapping understory bamboo across large geographic areas, but also for monitoring its temporal changes. In addition, while previous approaches focused on either a group of similar species or a single species, this approach has the ability to separate individual species. Since different bamboo species contribute unequally to the diet of the giant panda in 131 different seasons, knowledge on the distribution of individual bamboo species is useful for better measuring panda habitat. For example, it could be helpful for examining seasonal habitat use of the giant panda and for assessing the influence of species-specific dynamics of understory bamboo on panda habitat. Since land surface phenology obtained from MODIS data contains spatially and temporally continuous information on both forests and understory bamboo, which are important components of panda habitat, it is particularly useful for investigating the spatiotemporal dynamics of panda habitat. Using Wolong Nature Reserve as a case study, I showed the usefulness of land surface phenology in a species distribution model (SDM) for monitoring temporal dynamics of panda habitat. In addition, I also examined the effects of different predictor variables portraying land surface phenology on a model’s temporal transferability, which is an important, but often overlooked, characteristic of SDMs. Taking these effects into consideration, a phenology-based habitat model constitutes a useful tool for measuring panda habitat conditions across space and time, identifying potential threats to panda habitat and populations, evaluating the effectiveness of current conservation practices, and guiding conservation planning. To demonstrate the applications of the phenology-based model to panda conservation, I evaluated the effectiveness of a conservation program by investigating the spatiotemporal dynamics of panda habitat from 2001 to 2007 in Wolong Nature Reserve. Using spatial autoregressive models, I identified a significantly positive effect of the Natural Forest Conservation Program (NFCP) on the observed habitat improvement after controlling for confounding effects introduced by spatial autocorrelation and several biophysical (e.g., topography, forest cover and soil attributes) and anthropogenic factors (e.g., distance to local 132 household as a surrogate of human activity intensity). Results of this study suggest that an innovative implementation of NFCP, which encourages active participation of local residents in forest monitoring by providing direct payments and enhancing social norms among households, is an effective instrument for panda habitat conservation. Since local residents do not directly participate in the current implementation of NFCP in other places of China, better conservation of the giant panda and many other species could be achieved if the successful experience of engaging local residents with direct payments can be spread to other places. While the current conservation programs have effectively reduced the threats of land use/cover change to panda habitat, I showed that climate change may become the next major threat to the survival of the species. Focusing on the food resources of the panda population in the Qinling Mountains region (i.e., three dominant understory bamboo species: Qinling arrow, dragon-head and wooden), an ensemble of panda habitat projections obtained from bioclimatic envelope models indicated a substantial loss of panda habitat due to a potential shortage of food under projected climate change in the 21st century. These results indicate a potentially big challenge for giant panda conservation in the face of climate change. Benefits from current conservation efforts (e.g., the improvement of panda habitat due to NFCP implementation) may be offset by climate-induced food shortages and habitat losses if the potential climate change impacts are not incorporated into conservation planning or addressed by proactive conservation practices. In response to this potential threat, it is important to maintain and increase ecological connectivity not only among current habitat patches, but also among current and projected patches. Due to unavoidable uncertainties associated with climate change impact assessments, adaptive conservation strategies should be developed so that they can be adjusted when uncertainties are resolved and/or new knowledge becomes available. Finally, since almost all 133 suitable panda habitat was projected to disappear in the Qinling Mountains under the most pessimistic climate situation modeled, more aggressive actions, such as assisted colonization for bamboo and the panda, may need to be seriously considered. While the success of captive breeding and ongoing re-introduction projects for giant pandas provide a good foundation for their assisted colonization, there is an urgent need to better understand the risks and challenges associated with this action. Conclusions The research described in this dissertation demonstrates the usefulness of species distribution modeling and remote sensing technology for providing detailed information on spatiotemporal dynamics of species habitat across large geographic regions and over long time periods. This information is usually unavailable through traditional field surveys, but is crucial for effective and efficient conservation of species and their habitats, especially under the current rapidly changing environment. Although the remote sensing and modeling approaches developed in this dissertation were applied to investigate the spatiotemporal dynamics of understory bamboo and panda habitat in China, they also have the potential to be applied to other species and geographic regions. For example, many other understory plants, including some invasive species, also have substantial effects on the species compositions, structure and function of forest ecosystems, and may form dense understory layers which can cause detectable differences in land surface phenology. Therefore, SDMs with land surface phenology variables may be applicable to investigating the spatiotemporal dynamics of those understory plants and the habitats of animal species which rely on them, and thus may provide essential information for better managing forest ecosystems and conserving forest biodiversity. 134 However, applications of land surface phenology for investigating species distribution and biodiversity patterns face some limitations. For example, there is a trade-off between spatial and temporal resolutions of remotely sensed data. To take the advantage of high temporal resolution imagery in portraying land surface phenology, a compromise of spatial resolutions is unavoidable. Although the spatial resolution of MODIS data (250×250 m) is fine enough for large-scale habitat evaluations (e.g., range-wide evaluation) for many species, it may be too coarse to reveal habitat heterogeneity suitable for local-scale applications (e.g., studies on habitat selection of animal individuals). Another limitation relating to MODIS data is the relatively short temporal depth (i.e., data are available after February, 2000). Although the time series data acquired by the Advanced Very High Resolution Radiometer (AVHRR) make it possible to apply the approach to investigating species habitat back to the 1980s, the spatial resolution of AVHRR time series is even coarser than that of MODIS data. Coarse spatial resolutions challenge the possibility of linking information on land surface phenology with biodiversity information obtained on the ground. Further research on advanced remote sensing techniques, such as multi-sensor data fusion which combines information from high spatial resolution imagery with that from high temporal resolution data, may provide a solution for this limitation. Besides these technique limitations, further research is also needed for applying land surface phenology to biodiversity studies. While land surface phenology has mostly been used for studying human-induced land use/cover changes and climate-related phenology dynamics of plant species, its use for investigating spatiotemporal dynamics of species habitat and biodiversity patterns has not received enough attention. More studies are needed to further evaluate the usefulness of land surface phenology for detecting the distribution and characterizing the habitat of different species in different geographic settings. In addition, 135 further studies are also needed to understand how the biophysical and phenological characteristics of individual plant species, as well as species assemblages, affect land surface phenology obtained from remotely sensed data, and to understand how different vegetation characteristics and plant species assemblages affect the habitat quality of animal species. This knowledge will be useful for effectively extracting more detailed information relevant to biodiversity from land surface phenology. In addition to land surface phenology, information on many other vegetation characteristics and other environmental conditions (e.g., topography and climatic conditions) can also be obtained through remotely sensed data. While most previous studies about remote detection of habitat dynamics focused on changes in land cover types, SDMs with satellitederived continuous variables have the potential for revealing more detailed habitat changes, even within a single land cover type. This dissertation provides useful guidance for the integration of remote sensing technology and species distribution modeling with a focus on model transferability and algorithms for generating variables from high temporal resolution remotely sensed data. However, further studies are needed to understand the influence of diverse satellitederived variables and modeling approaches on SDMs’ characteristics and their performance on characterizing species habitats and monitoring their temporal dynamics. Beyond the contribution in the development of tools for biodiversity research and conservation, this dissertation also has broad implications for biodiversity conservation. First, this dissertation indicates the importance of viewing humans as one important component in coupled human and natural systems (CHANS) and encouraging active participation of local people for biodiversity conservation. Conservation practices excluding human access to natural resources without sufficient respect for local people’s need for them may cause serious conflicts 136 between biodiversity conservation and local people’s livelihoods, and thus reduce the effectiveness of conservation efforts. While this dissertation evaluated the effectiveness of a payment-for-ecosystem-services scheme for panda habitat conservation, further studies about the underlying mechanism of the effectiveness are warranted. In addition, studies about effects of this conservation strategy on human systems and interactions between human and natural systems are also needed to provide comprehensive information for guiding policy making. Better understanding of the complexity of human-nature interactions under a CHANS framework is necessary for developing more effective and efficient conservation strategies. This dissertation also highlights the urgent need to consider the potential impacts of climate change on biodiversity conservation. Without adjustments for potential changes in species distributions and biodiversity patterns under changed climates, current conservation strategies may lose their effectiveness in the future and conservation benefits obtained from them may be offset. While information on the potential impacts of climate change is the foundation for developing conservation strategies in the face of climate change, SDMs have been widely used for assessing the impacts of climate change on species distribution and biodiversity patterns, and for informing conservation planning under a changed climate. However, SDMs are usually criticized for establishing correlative, not causal, relationships between species distribution and climate conditions, being based on many assumptions whose validity is questionable (e.g., consistent and equilibrant species-climate relationships), and omitting important factors affecting species-climate relationships (e.g., biotic interactions among species and dispersal ability). For example, this dissertation indicates that not only the survival of individual species (e.g., bamboo species), but also their ecological functions and their interactions with other species (e.g., providing food for giant pandas) are vulnerable to climate change. Thus it underscores the 137 importance of incorporating cascading effects of climate change through inter-specific relationships into impact assessments and associated conservation planning. Although process-based dynamic models may have a better ability to reflect mechanistic influences of climate conditions on species and to incorporate relevant biological and ecological processes, the required knowledge for establishing such models is currently unavailable for most species on Earth. Since developing conservation strategies to cope with climate change impacts cannot wait until all knowledge becomes available, SDMs provide practical alternatives for guiding conservation planning. However, given the limitations of SDMs, especial caution should be taken when using information provided by them. In particular, uncertainty about climate change impacts introduced by SDMs and associated assumptions should be taken into conservation. In conclusion, in the face of the current crisis of biodiversity loss, knowledge about the spatial patterns and temporal dynamics of species distributions are crucial for effective and efficient biodiversity conservation. This dissertation makes substantial contributions to giant panda conservation by providing effective and practical approaches and essential information on the spatiotemporal dynamics of its habitat. However, further basic and applied research is needed to foster understanding of the spatiotemporal dynamics of biodiversity and its underlying drivers, and to build links between scientific knowledge and real-world conservation under a changing environment. 138 REFERENCES 139 REFERENCES Adams, W.M., Aveling, R., Brockington, D., Dickson, B., Elliott, J., Hutton, J., Roe, D., Vira, B., & Wolmer, W. (2004). Biodiversity conservation and the eradication of poverty. Science, 306, 1146-1149 Aguilera, A.M., Escabias, M., & Valderrama, M.J. (2006). Using principal components for estimating logistic regression with high-dimensional multicollinear data. Computational Statistics & Data Analysis, 50, 1905-1924 Ahl, D.E., Gower, S.T., Burrows, S.N., Shabanov, N.V., Myneni, R.B., & Knyazikhin, Y. (2006). Monitoring spring canopy phenology of a deciduous broadleaf forest using MODIS. Remote Sensing of Environment, 104, 88-95 An, L., Liu, J., Ouyang, Z., Linderman, M., Zhou, S., & Zhang, H. (2001). Simulating demographic and socioeconomic processes on household level and implications for giant panda habitats. Ecological Modelling, 140, 31-49 An, L., Lupi, F., Liu, J., Linderman, M.A., & Huang, J. (2002). Modeling the choice to switch from fuelwood to electricity: Implications for giant panda habitat conservation. Ecological Economics, 42, 445-457 Andersson, K., & Gibson, C.C. (2007). Decentralized governance and environmental change: Local institutional moderation of deforestation in Bolivia. Journal of Policy Analysis and Management, 26, 99-123 Araújo, M.B., Alagador, D., Cabeza, M., Nogués-Bravo, D., & Thuiller, W. (2011). Climate change threatens European conservation areas. Ecology Letters, 14, 484-492 Araújo, M.B., Cabeza, M., Thuiller, W., Hannah, L., & Williams, P.H. (2004). Would climate change drive species out of reserves? An assessment of existing reserve-selection methods. Global Change Biology, 10, 1618-1626 Araújo, M.B., & Guisan, A. (2006). Five (or so) challenges for species distribution modelling. Journal of Biogeography, 33, 1677-1688 Araújo, M.B., Pearson, R.G., Thuiller, W., & Erhard, M. (2005). Validation of species-climate impact models under climate change. Global Change Biology, 11, 1504-1513 140 Araújo, M.B., Thuiller, W., & Pearson, R.G. (2006). Climate warming and the decline of amphibians and reptiles in Europe. Journal of Biogeography, 33, 1712-1728 Asner, G.P., Hughes, R.F., Vitousek, P.M., Knapp, D.E., Kennedy-Bowdoin, T., Boardman, J., Martin, R.E., Eastwood, M., & Green, R.O. (2008). Invasive plants transform the threedimensional structure of rain forests. Proceedings of the National Academy of Sciences of the United States of America, 105, 4519-4523 Asner, G.P., & Vitousek, P.M. (2005). Remote analysis of biological invasion and biogeochemical change. Proceedings of the National Academy of Sciences of the United States of America, 102, 4383-4386 Bakkenes, M., Alkemade, J.R.M., Ihle, F., Leemans, R., & Latour, J.B. (2002). Assessing effects of forecasted climate change on the diversity and distribution of European higher plants for 2050. Global Change Biology, 8, 390-407 Balmford, A., & Bond, W. (2005). Trends in the state of nature and their implications for human well-being. Ecology Letters, 8, 1218-1234 Balmford, A., Green, R.E., & Jenkins, M. (2003). Measuring the changing state of nature. Trends in Ecology & Evolution, 18, 326-330 Bamford, A.J., et al. (2009). Trade-offs between specificity and regional generality in habitat association models: a case study of two species of African vulture. Journal of Applied Ecology, 46, 852-860 Barnosky, A.D., et al. (2011). Has the Earth's sixth mass extinction already arrived? Nature, 471, 51-57 Bearer, S., Linderman, M., Huang, J., An, L., He, G., & Liu, J. (2008). Effects of fuelwood collection and timber harvesting on giant panda habitat use. Biological Conservation, 141, 385393 Beck, P.S.A., Atzberger, C., Høgda, K.A., Johansen, B., & Skidmore, A.K. (2006). Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI. Remote Sensing of Environment, 100, 321-334 Berk, R.A. (2006). An introduction to ensemble methods for data analysis. Sociological Methods & Research, 34, 263-295 141 Berkes, F. (2004). Rethinking Community-Based Conservation. Conservation Biology, 18, 621630 Boyce, M.S., Vernier, P.R., Nielsen, S.E., & Schmiegelow, F.K.A. (2002). Evaluating resource selection functions. Ecological Modelling, 157, 281-300 Branton, M., & Richardson, J.S. (2011). Assessing the value of the umbrella-species concept for conservation planning with meta-analysis. Conservation Biology, 25, 9-20 Butchart, S.H.M., et al. (2010). Global biodiversity: Indicators of recent declines. Science, 328, 1164-1168 Bystriakova, N., & Kapos, V. (2006). Bamboo diversity: the need for a Red List review. Biodiversity, 6, 12-16 Carvalho, S.B., Brito, J.C., Crespo, E.G., Watts, M.E., & Possingham, H.P. (2011). Conservation planning under climate change: Toward accounting for uncertainty in predicted species distributions to increase confidence in conservation investments in space and time. Biological Conservation, 144, 2020-2030 Ceballos, G., & Ehrlich, P.R. (2002). Mammal population losses and the extinction crisis. Science, 296, 904-907 Chape, S., Harrison, J., Spalding, M., & Lysenko, I. (2005). Measuring the extent and effectiveness of protected areas as an indicator for meeting global biodiversity targets. Philosophical Transactions of the Royal Society B: Biological Sciences, 360, 443-455 Chastain Jr., R.A., & Townsend, P.A. (2007). Use of Landsat ETM and topographic data to characterize evergreen understory communities in Appalachian deciduous forests. Photogrammetric Engineering & Remote Sensing, 73, 563-575 Chen, I.C., Hill, J.K., Ohlemüller, R., Roy, D.B., & Thomas, C.D. (2011). Rapid range shifts of species associated with high levels of climate warming. Science, 333, 1024-1026 Chen, X., Lupi, F., He, G., & Liu, J. (2009). Linking social norms to efficient conservation investment in payments for ecosystem services. Proceedings of the National Academy of Sciences of the United States of America, 106, 11812-11817 142 Coetzee, B.W.T., Robertson, M.P., Erasmus, B.F.N., van Rensburg, B.J., & Thuiller, W. (2009). Ensemble models predict Important Bird Areas in southern Africa will become less effective for conserving endemic birds under climate change. Global Ecology and Biogeography, 18, 701-710 Cohen, J. (1960). A coefficient of agreement for nominal scales. Education and Psychological Measurement, 20, 37-46 de Beurs, K.M., & Henebry, G.M. (2004). Land surface phenology, climatic variation, and institutional change: Analyzing agricultural land cover change in Kazakhstan. Remote Sensing of Environment, 89, 497-509 de Beurs, K.M., & Henebry, G.M. (2005). Land surface phenology and temperature variation in the International Geosphere-Biosphere Program high-latitude transects. Global Change Biology, 11, 779-790 Deal, R.L. (2007). Management strategies to increase stand structural diversity and enhance biodiversity in coastal rainforests of Alaska. Biological Conservation, 137, 520-532 DeFries, R., Hansen, M., & Townshend, J. (1995). Global discrimination of land cover types from metrics derived from AVHRR pathfinder data. Remote Sensing of Environment, 54, 209222 Díaz, I.A., Armesto, J.J., Reid, S., Sieving, K.E., & Willson, M.F. (2005). Linking forest structure and composition: avian diversity in successional forests of Chiloé Island, Chile. Biological Conservation, 123, 91-101 Dietz, T., Ostrom, E., & Stern, P.C. (2003). The struggle to govern the commons. Science, 302, 1907-1912 Eklundh, L., Johansson, T., & Solberg, S. (2009). Mapping insect defoliation in Scots pine with MODIS time-series data. Remote Sensing of Environment, 113, 1566-1573 Elith, J., et al. (2006). Novel methods improve prediction of species' distributions from occurrence data. Ecography, 29, 129-151 Elith, J., & Leathwick, J.R. (2009). Species distribution models: Ecological explanation and prediction across space and time. Annual Review of Ecology, Evolution, and Systematics, 40, 677-697 143 Engel, S., Pagiola, S., & Wunder, S. (2008). Designing payments for environmental services in theory and practice: An overview of the issues. Ecological Economics, 65, 663-674 Engler, R., Guisan, A., & Rechsteiner, L. (2004). An improved approach for predicting the distribution of rare and endangered species from occurrence and pseudo-absence data. Journal of Applied Ecology, 41, 263-274 Eriksson, H.M., Eklundh, L., Kuusk, A., & Nilson, T. (2006). Impact of understory vegetation on forest canopy reflectance and remotely sensed LAI estimates. Remote Sensing of Environment, 103, 408-418 Estades, C.F., & Temple, S.A. (1999). Deciduous-forest bird communities in a fragmented landscape dominated by exotic pine plantations. Ecological Applications, 9, 573-585 FAO (1995). Non-Wood Forest Products for Rural Income and Sustainable Development. NonWood Forest Products No. 7. Rome, Italy: Food and Agriculture Organization Farr, T.G., et al. (2007). The Shuttle Radar Topography Mission. Reviews of Geophysics, 45, RG2004 Ferraro, P.J., & Kiss, A. (2002). Direct payments to conserve biodiversity. Science, 298, 17181719 Fielding, A.H., & Bell, J.F. (1997). A review of methods for the assessment of prediction errors in conservation presence/absence models. Environmental Conservation, 24, 38-49 Fox, J., & Weisberg, S. (2011). An R Companion to Applied Regression. Thousand Oaks, CA: Sage Franzluebbers, A.J., Haney, R.L., Honeycutt, C.W., Arshad, M.A., Schomberg, H.H., & Hons, F.M. (2001). Climatic influences on active fractions of soil organic matter. Soil Biology and Biochemistry, 33, 1103-1111 Friedl, M., Henebry, G., Reed, B., Huete, A., White, M., Morisette, J., Nemani, R., Zhang, X., & Myneni, R. (2006). Land surface phenology NASA white paper. ftp://ftp.iluci.org/Land_ESDR/Phenology_Friedl_whitepaper.pdf 144 Gilliam, F.S. (2007). The ecological significance of the herbaceous layer in temperate forest ecosystems. Bioscience, 57, 845-858 Gitelson, A.A. (2004). Wide Dynamic Range Vegetation Index for remote quantification of biophysical characteristics of vegetation. Journal of Plant Physiology, 161, 165-173 Graham, M.H. (2003). Confronting multicollinearity in ecological multiple regression. Ecology, 84, 2809-2815 Griscom, B.W., & Ashton, P.M.S. (2003). Bamboo control of forest succession: Guadua sarcocarpa in Southeastern Peru. Forest Ecology and Management, 175, 445-454 Guisan, A., & Thuiller, W. (2005). Predicting species distribution: offering more than simple habitat models. Ecology Letters, 8, 993-1009 Guisan, A., & Zimmermann, N.E. (2000). Predictive habitat distribution models in ecology. Ecological Modelling, 135, 147-186 Hagar, J.C. (2007). Wildlife species associated with non-coniferous vegetation in Pacific Northwest conifer forests: A review. Forest Ecology and Management, 246, 108-122 Hanley, J.A., & Mcneil, B.J. (1982). The meaning and use of the area under a ROC curve. Radiology, 143, 29-36 Hannah, L., Midgley, G., Andelman, S., Araujo, M., Hughes, G., Martinez-Meyer, E., Pearson, R., & Williams, P. (2007). Protected area needs in a changing climate. Frontiers in Ecology and the Environment, 5, 131-138 Hannah, L., Midgley, G.F., & Millar, D. (2002). Climate change-integrated conservation strategies. Global Ecology and Biogeography, 11, 485-495 Hansen, M.C., Defries, R.S., Townshend, J.R.G., & Sohlberg, R. (2000). Global land cover classification at 1 km spatial resolution using a classification tree approach. International Journal of Remote Sensing, 21, 1331-1364 He, G., et al. (2009). Spatial and temporal patterns of fuelwood collection in Wolong Nature Reserve: Implications for panda conservation. Landscape and Urban Planning, 92, 1-9 145 Heller, N.E., & Zavaleta, E.S. (2009). Biodiversity management in the face of climate change: A review of 22 years of recommendations. Biological Conservation, 142, 14-32 Henebry, G.M., Viña, A., & Gitelson, A.A. (2004). The Wide Dynamic Range Vegetation Index and its potential utility for gap analysis. GAP Analysis Program Bulletin, 12, 50-56 Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., & Jarvis, A. (2005). Very high resolution interpolated climate surfaces for global land areas. International Journal of Climatology, 25, 1965-1978 Hoegh-Guldberg, O., Hughes, L., McIntyre, S., Lindenmayer, D.B., Parmesan, C., Possingham, H.P., & Thomas, C.D. (2008). Assisted colonization and rapid climate change. Science, 321, 345-346 Hoffmann, M., et al. (2010). The impact of conservation on the status of the world's vertebrates. Science, 330, 1503-1509 Holling, C.S. (1978). Adaptive Environmental Assessment and Management. New York: Wiley Hooper, D.U., et al. (2005). Effects of biodiversity on ecosystem functioning: a consensus of current knowledge. Ecological Monographs, 75, 3-35 Huang, C., & Asner, G. (2009). Applications of remote sensing to alien invasive plant studies. Sensors, 9, 4869-4889 Huang, C., Geiger, E.L., Van Leeuwen, W.J.D., & Marsh, S.E. (2009). Discrimination of invaded and native species sites in a semi-desert grassland using MODIS multi-temporal data. International Journal of Remote Sensing, 30, 897 - 917 Hughes, R., & Flintan, F. (2001). Integrating Conservation and Development Experience: A Review and Bibliography of the ICDP Literature. In, Biodiversity and Lovelihoods Issue Papers (p. 24): International Institute for Environment and Development Hull, V., et al. (2011). Evaluating the efficacy of zoning designations for protected area management. Biological Conservation, 144, 3028-3037 Hutchinson, G.E. (1957). Concluding remarks. Cold Spring Harbor Symposia on Quantitative Biology, 22, 415-427 146 Ichii, K., Kawabata, A., & Yamaguchi, Y. (2002). Global correlation analysis for NDVI and climatic variables and NDVI trends: 1982-1990. International Journal of Remote Sensing, 23, 3873-3878 IPCC (2001). Climate Change 2001: The Physical Science Basis. New York, NY, USA: Cambridge University Press IPCC (2007). Climate Change 2007: Impacts, Adaptation and Vulnerability. New York: Cambridge University Press Iqbal, M. (1993). International Trade in Non-Wood Forest Products. An Overview. Rome, Italy: Food and Agriculture Organization IUCN (2011). International Union for Conservation of Nature Red List. http://www.iucn.org/about/work/programmes/species/red_list/ Janzen, D.H. (1976). Why bamboos wait so long to flower. Annual Review of Ecology and Systematics, 7, 347-391 Jenkins, C.N., & Joppa, L. (2009). Expansion of the global terrestrial protected area system. Biological Conservation, 142, 2166-2174 Jensen, J.R. (2007). Remote Sensing of the Environment : An Earth Resource Perspective. Upper Saddle River, NJ: Pearson Prentice Hall Jetz, W., McPherson, J.M., & Guralnick, R.P. (in press). Integrating biodiversity distribution knowledge: toward a global map of life. Trends in ecology & evolution (Personal edition) Jetz, W., Sekercioglu, C.H., & Watson, J.E.M. (2008). Ecological correlates and conservation implications of overestimating species geographic ranges. Conservation Biology, 22, 110-119 Ji, L., & Gallo, K. (2006). An agreement coefficient for image comparison. Photogrammetric Engineering & Remote Sensing, 72, 823-833 Jian, J., Jiang, H., Zhou, G.M., Jiang, Z.S., Yu, S.Q., Peng, S.L., Liu, S.Y., & Wang, J.X. (2011). Mapping the vegetation changes in giant panda habitat using Landsat remotely sensed data. International Journal of Remote Sensing, 32, 1339-1356 147 Jönsson, P., & Eklundh, L. (2004). TIMESAT - a program for analysing time-series of satellite sensor data. Computers & Geosciences, 30, 833-845 Jönsson, P., & Eklundh, L. (2006). TIMESAT - A Program for Analyzing Time-Series of Satellite Sensor Data. Users Guide for TIMESAT 2.3 Kerr, J.T., & Ostrovsky, M. (2003). From space to species: ecological applications for remote sensing. Trends in Ecology & Evolution, 18, 299-305 Koltunov, A., Ustin, S.L., Asner, G.P., & Fung, I. (2009). Selective logging changes forest phenology in the Brazilian Amazon: Evidence from MODIS image time series analysis. Remote Sensing of Environment, 113, 2431-2440 Korpela, I.S. (2008). Mapping of understory lichens with airborne discrete-return LiDAR data. Remote Sensing of Environment, 112, 3891-3897 Kremen, C., Merenlender, A.M., & Murphy, D.D. (1994). Ecological monitoring: a vital need for integrated conservation and development programs in the tropics. Conservation Biology, 8, 388397 Kummerow, C., Barnes, W., Kozu, T., Shiue, J., & Simpson, J. (1998). The Tropical Rainfall Measuring Mission (TRMM) Sensor Package. Journal of Atmospheric and Oceanic Technology, 15, 809-817 La Sorte, F.A., & Jetz, W. (2010). Projected range contractions of montane biodiversity under global warming. Proceedings of the Royal Society B: Biological Sciences, 277, 3401-3410 Landis, J.R., & Koch, G.G. (1977). The measurement of observer agreement for categorical data. Biometrics, 33, 159-174 Lefsky, M.A., Cohen, W.B., Parker, G.G., & Harding, D.J. (2002). Lidar remote sensing for ecosystem studies. Bioscience, 52, 19-30 Lengyel, S., Kobler, A., Kutnar, L., Framstad, E., Henry, P.-Y., Babij, V., Gruber, B., Schmeller, D., & Henle, K. (2008). A review and a framework for the integration of biodiversity monitoring at the habitat level. Biodiversity and Conservation, 17, 3341-3356 148 Lenney, M.P., Woodcock, C.E., Collins, J.B., & Hamdi, H. (1996). The status of agricultural lands in Egypt: The use of multitemporal NDVI features derived from Landsat TM. Remote Sensing of Environment, 56, 8-20 Li, C. (1997). A Study of Staple Food Bamboo for the Giant Panda. Guiyang, China: Guizhou Scientific Publishing House Linderman, M., Bearer, S., An, L., Tan, Y., Ouyang, Z., & Liu, J. (2005). The effects of understory bamboo on broad-scale estimates of giant panda habitat. Biological Conservation, 121, 383-390 Linderman, M., Liu, J., Qi, J., An, L., Ouyang, Z., Yang, J., & Tan, Y. (2004). Using artificial neural networks to map the spatial distribution of understorey bamboo from remote sensing data. International Journal of Remote Sensing, 25, 1685-1700 Liu, C., Berry, P.M., Dawson, T.P., & Pearson, R.G. (2005a). Selecting thresholds of occurrence in the prediction of species distributions. Ecography, 28, 385-393 Liu, J. (2010). China's road to sustainability. Science, 328, 50 Liu, J., Daily, G.C., Ehrlich, P.R., & Luck, G.W. (2003a). Effects of household dynamics on resource consumption and biodiversity. Nature, 421, 530-533 Liu, J., & Diamond, J. (2008). Revolutionizing China's environmental protection. Science, 319, 37-38 Liu, J., et al. (2007). Complexity of coupled human and natural systems. Science, 317, 15131516 Liu, J., Li, S., Ouyang, Z., Tam, C., & Chen, X. (2008). Ecological and socioeconomic effects of China's policies for ecosystem services. Proceedings of the National Academy of Sciences of the United States of America, 105, 9477-9482 Liu, J., Linderman, M., Ouyang, Z., An, L., Yang, J., & Zhang, H. (2001). Ecological degradation in protected areas: the case of Wolong Nature Reserve for giant pandas. Science, 292, 98-101 149 Liu, J., Ouyang, Z., Pimm, S.L., Raven, P.H., Wang, X., Miao, H., & Han, N. (2003b). Protecting China's biodiversity. Science, 300, 1240-1241 Liu, J., Ouyang, Z., Taylor, W.W., Groop, R., Tan, Y., & Zhang, H. (1999). A framework for evaluating the effects of human factors on wildlife habitat: the case of giant pandas. Conservation Biology, 13, 1360-1370 Liu, J., & Raven, P.H. (2010). China's environmental challenges and implications for the world. Critical Reviews in Environmental Science and Technology, 40, 823-851 Liu, X., & Li, J. (2008). Scientific solutions for the functional zoning of nature reserves in China. Ecological Modelling, 215, 237-246 Liu, X.H., Toxopeus, A.G., Skidmore, A.K., Shao, X.M., Dang, G.D., Wang, T.J., & Prins, H.H.T. (2005b). Giant panda habitat selection in Foping Nature Reserve, China. Journal of Wildlife Management, 69, 1623-1632 Lloyd, D. (1990). A phenological classification of terrestrial vegetation cover using shortwave vegetation index imagery. International Journal of Remote Sensing, 11, 2269-2279 Lobo, J.M., Jiménez-Valverde, A., & Real, R. (2008). AUC: a misleading measure of the performance of predictive distribution models. Global Ecology and Biogeography, 17, 145-151 Loucks, C., & Wang, H. (2004). Assessing the habitat and distribution of the giant panda: methods and issues. In D. Lindburg & K. Baragona (Eds.), Giant Pandas: Biology and Conservation (pp. 149-154). Berkeley and Los Angeles, CA: University of California Press Loucks, C.J., Zhi, L., Dinerstein, E., Wang, D., Fu, D., & Wang, H. (2003). The giant pandas of the Qinling Mountains, China: a case study in designing conservation landscapes for elevational migrants. Conservation Biology, 17, 558-565 Loucks, C.J., Zhi, L., Dinerstein, E., Wang, H., Olson, D.M., Zhu, C., & Wang, D. (2001). Giant pandas in a changing landscape. Science, 294, 1465 Lü, Z., et al. (2001). Patterns of genetic diversity in remaining giant panda populations. Conservation Biology, 15, 1596-1607 150 Mac Nally, R. (2000). Regression and model-building in conservation biology, biogeography and ecology: The distinction between – and reconciliation of – ‘predictive’ and ‘explanatory’ models. Biodiversity and Conservation, 9, 655-671 Margules, C.R., & Pressey, R.L. (2000). Systematic conservation planning. Nature, 405, 243-253 Marini, M.Â., Barbet-Massin, M., Martinez, J., Prestes, N.P., & Jiguet, F. (2010). Applying ecological niche modelling to plan conservation actions for the Red-spectacled Amazon (Amazona pretrei). Biological Conservation, 143, 102-112 Mason, D.C., Anderson, G.Q.A., Bradbury, R.B., Cobby, D.M., Davenport, I.J., Vandepoll, M., & Wilson, J.D. (2003). Measurement of habitat predictor variables for organism-habitat models using remote sensing and image segmentation. International Journal of Remote Sensing, 24, 2515 - 2532 McKenzie, N.J., & Ryan, P.J. (1999). Spatial prediction of soil properties using environmental correlation. Geoderma, 89, 67-94 McPherson, J.M., Jetz, W., & Rogers, D.J. (2004). The effects of species' range sizes on the accuracy of distribution models: ecological phenomenon or statistical artefact? Journal of Applied Ecology, 41, 811-823 McShane, T.O., et al. (2011). Hard choices: Making trade-offs between biodiversity conservation and human well-being. Biological Conservation, 144, 966-972 Millennium Ecosystem Assessment (2005). Ecosystems and Human Well-Being. Washington, D.C.: Island Press Miller, T.R., Minteer, B.A., & Malan, L.-C. (2011). The new conservation debate: the view from practical ethics. Biological Conservation, 144, 948-957 Ministry of Environmental Protection of China (2010). China National Biodiversity Conservation Strategy and Action Plan (2011-2030): China Environmental Science Press Moore, I.D., Gessler, P.E., Nielsen, G.A., & Peterson, G.A. (1993). Soil attribute prediction using terrain analysis. Soil Sci. Soc. Am. J., 57, 443-452 151 Mora, C., Tittensor, D.P., Adl, S., Simpson, A.G.B., & Worm, B. (2011). How many species are there on Earth and in the ocean? PLoS Biology, 9, e1001127 Morisette, J.T., Jarnevich, C.S., Ullah, A., Cai, W., Pedelty, J.A., Gentle, J.E., Stohlgren, T.J., & Schnase, J.L. (2006). A tamarisk habitat suitability map for the continental United States. Frontiers in Ecology and the Environment, 4, 11-17 Morisette, J.T., et al. (2009). Tracking the rhythm of the seasons in the face of global change: phenological research in the 21st century. Frontiers in Ecology and the Environment, 7, 253-260 Myers, N., Mittermeier, R.A., Mittermeier, C.G., da Fonseca, G.A.B., & Kent, J. (2000). Biodiversity hotspots for conservation priorities. Nature, 403, 853-858 Naughton-Treves, L., Holland, M.B., & Brandon, K. (2005). The role of protected areas in conserving biodiversity and sustaining local livelihoods. Annual Review of Environment and Resources, 30, 219-252 Nilsson, M.-C., & Wardle, D.A. (2005). Understory vegetation as a forest ecosystem driver: evidence from the northern Swedish boreal forest. Frontiers in Ecology and the Environment, 3, 421-428 Olson, D.M., & Dinerstein, E. (1998). The Global 200: a representation approach to conserving the Earth’s most biologically valuable ecoregions. Conservation Biology, 12, 502-515 Pan, W., Gao, Z.S., & Lü, Z. (1988). The Giant Panda's Natural Refuge in the Qinling Mountains. Beijing: Peking University Press Pan, W., Lu, Z., Zhu, X.J., Wang, D.J., Wang, H., Long, Y., Fu, D.L., & Zhou, X. (2001). The Opportunity for the Giant Panda to Exist. Beijing, China: Peking University Press Parmesan, C. (2006). Ecological and evolutionary responses to recent climate change. Annual Review of Ecology, Evolution, and Systematics, 37, 637-669 Parmesan, C., & Yohe, G. (2003). A globally coherent fingerprint of climate change impacts across natural systems. Nature, 421, 37-42 152 Pattanayak, S.K., Wunder, S., & Ferraro, P.J. (2010). Show me the money: do payments supply environmental services in developing countries? Review of Environmental Economics and Policy, 4, 254-274 Pearce, J., & Ferrier, S. (2000). Evaluating the predictive performance of habitat models developed using logistic regression. Ecological Modelling, 133, 225-245 Pearson, R.G., & Dawson, T.P. (2003). Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful? Global Ecology and Biogeography, 12, 361-371 Peckham, S.D., Ahl, D.E., & Gower, S.T. (2009). Bryophyte cover estimation in a boreal black spruce forest using airborne lidar and multispectral sensors. Remote Sensing of Environment, 113, 1127-1132 Pereira, H.M., & Cooper, H.D. (2006). Towards the global monitoring of biodiversity change. Trends in Ecology & Evolution, 21, 123-129 Pereira, H.M., et al. (2010). Scenarios for Global Biodiversity in the 21st Century. Science, 330, 1496-1501 Peterson, A.T., Papes, M., & Eaton, M. (2007). Transferability and model evaluation in ecological niche modeling: a comparison of GARP and Maxent. Ecography, 30, 550-560 Peterson, A.T., Papes, M., & Soberón, J. (2008). Rethinking receiver operating characteristic analysis applications in ecological niche modeling. Ecological Modelling, 213, 63-72 Phillips, S.J. (2008). Transferability, sample selection bias and background data in presence-only modelling: a response to Peterson et al. (2007). Ecography, 31, 272-278 Phillips, S.J., Anderson, R.P., & Schapire, R.E. (2006). Maximum entropy modeling of species geographic distributions. Ecological Modelling, 190, 231-259 Phillips, S.J., & Dudík, M. (2008). Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography, 31, 161-175 153 Phillips, S.J., Dudík, M., Elith, J., Graham, C.H., Lehmann, A., Leathwick, J., & Ferrier, S. (2009). Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data. Ecological Applications, 19, 181-197 Pimm, S.L., Russell, G.J., Gittleman, J.L., & Brooks, T.M. (1995). The Future of Biodiversity. Science, 269, 347-350 Qin, Z., Taylor, A., & Cai, X. (1993). Dynamic Succession of the Bamboo and Forests in Panda Habitat in Wolong. Beijing: China Forestry Publishing House R Development Core Team (2011). R: A language and environment for statistical computing. In. Vienna, Austria: R Foundation for Statistical Computing Randin, C.F., Dirnbock, T., Dullinger, S., Zimmermann, N.E., Zappa, M., & Guisan, A. (2006). Are niche-based species distribution models transferable in space? Journal of Biogeography, 33, 1689-1703 Rautiainen, M., Suomalainen, J., Mottus, M., Stenberg, P., Voipio, P., Peltoniemi, J., & Manninen, T. (2007). Coupling forest canopy and understory reflectance in the Arctic latitudes of Finland. Remote Sensing of Environment, 110, 332-343 Reed, B.C., Brown, J.F., VanderZee, D., Loveland, T.R., Merchant, J.W., & Ohlen, D.O. (1994). Measuring phenological variability from satellite imagery. Journal of Vegetation Science, 5, 703714 Reid, D.G., Hu, J., Dong, S., Wang, W., & Huang, Y. (1989). Giant panda Ailuropoda melanoleuca behaviour and carrying capacity following a bamboo die-off. Biological Conservation, 49, 85-104 Reid, D.G., Taylor, A.H., Hu, J., & Qin, Z. (1991). Environmental influences on bamboo Bashania fangiana growth and implications for giant panda conservation. Journal of Applied Ecology, 28, 855-868 Resasco, J., Hale, A.N., Henry, M.C., & Gorchov, D.L. (2007). Detecting an invasive shrub in a deciduous forest understory using late-fall Landsat sensor imagery. International Journal of Remote Sensing, 28, 3739 - 3745 Root, T.L., Price, J.T., Hall, K.R., Schneider, S.H., Rosenzweig, C., & Pounds, J.A. (2003). Fingerprints of global warming on wild animals and plants. Nature, 421, 57-60 154 Roughgarden, J., Running, S.W., & Matson, P.A. (1991). What does remote sensing do for ecology? Ecology, 72, 1918-1922 Royo, A.A., & Carson, W.P. (2006). On the formation of dense understory layers in forests worldwide: consequences and implications for forest dynamics, biodiversity, and succession. Canadian Journal of Forest Research, 36, 1345-1362 Saatchi, S., Buermann, W., ter Steege, H., Mori, S., & Smith, T.B. (2008). Modeling distribution of Amazonian tree species and diversity using remote sensing measurements. Remote Sensing of Environment, 112, 2000-2017 Sala, O.E., et al. (2000). Biodiversity - Global biodiversity scenarios for the year 2100. Science, 287, 1770-1774 Savitzky, A., & Golay, M.J.E. (1964). Smoothing and differentiation of data by simplified least squares procedures. Analytical Chemistry, 36, 1627-1639 Schaller, G.B., Hu, J., Pan, W., & Zhu, J. (1985). The Giant Pandas of Wolong. Chicago, Illinois: University of Chicago Press Schowengerdt, R.A. (2007). Remote Sensing: Models and Methods for Image Processing. Burlington, MA: Academic Press Schwartz, M.D., & Reed, B.C. (1999). Surface phenology and satellite sensor-derived onset of greenness: an initial comparison. International Journal of Remote Sensing, 20, 3451 - 3457 Sheskin, D. (2000). Handbook of parametric and nonparametric statistical procedures. Boca Raton: Chapman & Hall/CRC Skov, F., & Svenning, J.C. (2004). Potential impact of climatic change on the distribution of forest herbs in Europe. Ecography, 27, 366-380 Soule, M.E. (1991). Conservation: tactics for a constant crisis. Science, 253, 744-750 State Forestry Administration (2006). The Third National Survey Report on Giant Panda in China. Beijing: Science Press 155 Sun, W., Liang, S., Xu, G., Fang, H., & Dickinson, R. (2008). Mapping plant functional types from MODIS data using multisource evidential reasoning. Remote Sensing of Environment, 112, 1010-1024 Swaisgood, R.R., Wei, F., Wildt, D.E., Kouba, A.J., & Zhang, Z. (2010). Giant panda conservation science: how far we have come. Biology Letters, 6, 143-145 Swaisgood, R.R., Wei, F.W., McShea, W.J., Wildt, D.E., Kouba, A.J., & Zhang, Z.J. (2011). Can science save the giant panda (Ailuropoda melanoleuca)? Unifying science and policy in an adaptive management paradigm. Integrative Zoology, 6, 290-296 Swets, J.A. (1988). Measuring the accuracy of diagnostic systems. Science, 240, 1285-1293 Tarou, L.R., Williams, J., Powell, D.M., Tabet, R., & Allen, M. (2005). Behavioral preferences for bamboo in a pair of captive giant pandas (Ailuropoda melanoleuca). Zoo Biology, 24, 177183 Taylor, A.H., Huang, J.Y., & Zhou, S.Q. (2004). Canopy tree development and undergrowth bamboo dynamics in old-growth Abies-Betula forests in Southwestern China: a 12-year study. Forest Ecology and Management, 200, 347-360 Taylor, A.H., Jang, S.W., Zhao, L.J., Liang, C.P., Miao, C.J., & Huang, J.Y. (2006). Regeneration patterns and tree species coexistence in old-growth Abies-Picea forests in southwestern China. Forest Ecology and Management, 223, 303-317 Taylor, A.H., & Qin, Z. (1987). Culm dynamics and dry matter production of bamboos in the Wolong and Tangjiahe giant panda reserves, Sichuan, China. Journal of Applied Ecology, 24, 419-433 Taylor, A.H., & Qin, Z. (1993). Structure and dynamics of bamboos in the Wolong Natural Reserve, China. American Journal of Botany, 80, 375-384 Thomas, C.D., et al. (2004). Extinction risk from climate change. Nature, 427, 145-148 Thuiller, W., Brotons, L., Araújo, M.B., & Lavorel, S. (2004). Effects of restricting environmental range of data to project current and future species distributions. Ecography, 27, 165-172 156 Townsend, P., & Walsh, S. (2001). Remote sensing of forested wetlands: application of multitemporal and multispectral satellite imagery to determine plant community composition and structure in southeastern USA. Plant Ecology, 157, 129-149 Tuanmu, M.-N., Viña, A., Bearer, S., Xu, W., Ouyang, Z., Zhang, H., & Liu, J. (2010). Mapping understory vegetation using phenological characteristics derived from remotely sensed data. Remote Sensing of Environment, 114, 1833-1844 Tuanmu, M.-N., Viña, A., Roloff, G.J., Liu, W., Ouyang, Z., Zhang, H., & Liu, J. (2011). Temporal transferability of wildlife habitat models: implications for habitat monitoring. Journal of Biogeography, 38, 1510-1523 Turner, W., Spector, S., Gardiner, N., Fladeland, M., Sterling, E., & Steininger, M. (2003). Remote sensing for biodiversity science and conservation. Trends in Ecology & Evolution, 18, 306-314 Urgenson, L.S., Reichard, S.H., & Halpern, C.B. (2009). Community and ecosystem consequences of giant knotweed (Polygonum sachalinense) invasion into riparian forests of western Washington, USA. Biological Conservation, 142, 1536-1541 Van der Putten, W.H., Macel, M., & Visser, M.E. (2010). Predicting species distribution and abundance responses to climate change: why it is essential to include biotic interactions across trophic levels. Philosophical Transactions of the Royal Society B: Biological Sciences, 365, 2025-2034 Vanreusel, W., Maes, D., & Van Dyck, H. (2007). Transferability of species distribution models: a functional habitat approach for two regionally threatened butterflies. Conservation Biology, 21, 201-212 Varela, S., Rodríguez, J., & Lobo, J.M. (2009). Is current climatic equilibrium a guarantee for the transferability of distribution model predictions? A case study of the spotted hyena. Journal of Biogeography, 36, 1645-1655 Vierling, K.T., Vierling, L.A., Gould, W.A., Martinuzzi, S., & Clawges, R.M. (2008). Lidar: shedding new light on habitat characterization and modeling. Frontiers in Ecology and the Environment, 6, 90-98 Viña, A., Bearer, S., Chen, X., He, G., Linderman, M., An, L., Zhang, H., Ouyang, Z., & Liu, J. (2007). Temporal changes in giant panda habitat connectivity across boundaries of Wolong Nature Reserve, China. Ecological Applications, 17, 1019-1030 157 Viña, A., Bearer, S., Zhang, H., Ouyang, Z., & Liu, J. (2008). Evaluating MODIS data for mapping wildlife habitat distribution. Remote Sensing of Environment, 112, 2160-2169 Viña, A., Chen, X., McConnell, W., Liu, W., Xu, W., Ouyang, Z., Zhang, H., & Liu, J. (2011). Effects of natural disasters on conservation policies: the case of the 2008 Wenchuan Earthquake, China. Ambio, 40, 274-284 Viña, A., Gitelson, A.A., Rundquist, D.C., Keydan, G., Leavitt, B., & Schepers, J. (2004a). Monitoring maize (Zea mays L.) phenology with remote sensing. Agronomy Journal, 96, 11391147 Viña, A., Henebry, G.M., & Gitelson, A.A. (2004b). Satellite monitoring of vegetation dynamics: Sensitivity enhancement by the Wide Dynamic Range Vegetation Index. Geophysical Research Letters, 31, L04503 Viña, A., Tuanmu, M.-N., Xu, W., Li, Y., Ouyang, Z., DeFries, R., & Liu, J. (2010). Range-wide analysis of wildlife habitat: Implications for conservation. Biological Conservation, 143, 19601969 Wan, Q.-H., Wu, H., & Fang, S.-g. (2005). A new subspecies of giant panda (Ailuropoda melanoleuca) from Shaanxi, China. Journal of Mammalogy, 86, 397-402 Wang, T.J., Skidmore, A.K., & Toxopeus, A.G. (2009a). Improved understorey bamboo cover mapping using a novel hybrid neural network and expert system. International Journal of Remote Sensing, 30, 965 - 981 Wang, T.J., Skidmore, A.K., Toxopeus, A.G., & Liu, X. (2009b). Understory bamboo discrimination using a winter image. Photogrammetric Engineering & Remote Sensing, 75, 3747 Wang, Y., Tao, J., & Zhong, Z. (2009c). Factors influencing the distribution and growth of dwarf bamboo, Fargesia nitida, in a subalpine forest in Wolong Nature Reserve, southwest China. Ecological Research, 24, 1013-1021 Wenchuan Statistics Bureau (2008). Wenchuan County Statistical Yearbook 2007 Whittaker, R.J., Araújo, M.B., Paul, J., Ladle, R.J., Watson, J.E.M., & Willis, K.J. (2005). Conservation Biogeography: assessment and prospect. Diversity and Distributions, 11, 3-23 158 Wiens, J.J., & Graham, C.H. (2005). Niche conservatism: integrating evolution, ecology, and conservation biology. Annual Review of Ecology, Evolution, and Systematics, 36, 519-539 Winkler, J.A., Guentchev, G.S., Liszewska, M., Perdinan, & Tan, P.-N. (2011). Climate scenario development and applications for local/regional climate change impact assessments: An overview for the non-climate scientist. Part II: Considerations when using climate change scenarios. Geography Compass, 5, 301-328 Wolong Nature Reserve (2005). History of the Development of Wolong Nature Reserve. Chengdu (in Chinese): Sichuan Science Publisher Wu, R., et al. (2011). Effectiveness of China's nature reserves in representing ecological diversity. Frontiers in Ecology and the Environment, 9, 383-389 Xu, H., et al. (2009a). China's progress toward the significant reduction of the rate of biodiversity loss. Bioscience, 59, 843-852 Xu, J., & Melick, D.R. (2007). Rethinking the effectiveness of public protected areas in southwestern China. Conservation Biology, 21, 318-328 Xu, W., Ouyang, Z., Viña, A., Zheng, H., Liu, J., & Xiao, Y. (2006). Designing a conservation plan for protecting the habitat for giant pandas in the Qionglai mountain range, China. Diversity and Distributions, 12, 610-619 Xu, W., Wang, X., Ouyang, Z., Zhang, J., Li, Z., Xiao, Y., & Zheng, H. (2009b). Conservation of giant panda habitat in South Minshan, China, after the May 2008 earthquake. Frontiers in Ecology and the Environment, 7, 353-358 Yan, C., Ai, X., Lin, E., & Yu, C. (2004). Research on the Impacts of Climate Change on Giant Panda Habitat Beijing: Weather Press Yin, R., & Yin, G. (2010). China's primary programs of terrestrial ecosystem restoration: Initiation, implementation, and challenges. Environmental Management, 45, 429-441 Zanini, F., Pellet, J., & Schmidt, B.R. (2009). The transferability of distribution models across regions: an amphibian case study. Diversity and Distributions, 15, 469-480 159 Zhang, X., Friedl, M.A., Schaaf, C.B., & Strahler, A.H. (2004). Climate controls on vegetation phenological patterns in northern mid- and high latitudes inferred from MODIS data. Global Change Biology, 10, 1133-1145 Zhang, X., Friedl, M.A., Schaaf, C.B., Strahler, A.H., Hodges, J.C.F., Gao, F., Reed, B.C., & Huete, A. (2003). Monitoring vegetation phenology using MODIS. Remote Sensing of Environment, 84, 471-475 Zhang, Z., Swaisgood, R.R., Zhang, S., Nordstrom, L.A., Wang, H., Gu, X., Hu, J., & Wei, F. (2011). Old-growth forest is what giant pandas really need. Biology Letters, 7, 403-406 Zharikov, Y., Elner, R., Shepherd, P., & Lank, D. (2009). Interplay between physical and predator landscapes affects transferability of shorebird distribution models. Landscape Ecology, 24, 129-144 Zimmermann, N.E., Edwards, T.C.J., Moisen, G.G., Frescino, T.S., & Blackard, J.A. (2007). Remote sensing-based predictors improve distribution models of rare, early successional and broadleaf tree species in Utah. Journal of Applied Ecology, 44, 1057-1067 Zweig, M., & Campbell, G. (1993). Receiver-operating characteristic (ROC) plots: a fundamental evaluation tool in clinical medicine [published erratum appears in Clin Chem 1993 Aug;39(8):1589]. Clinical Chemistry, 39, 561-577 160