NET PRIMARY PRODUCTIVITY OF A MOUNTAIN FOREST ECOSYSTEM AS AFFECTED BY CLIMATE AND TOPOGRAPHY

Response of terrestrial ecosystems to changing climate has become an issue of central importance for land managers and policymakers. Climate extremes and trends have a strong control on productivity of semi-arid mountain ecosystems. Located in a transition zone from continental type climate to relatively mild Black Sea type climate in Turkey, the Ilgaz Mountains with their rich biodiversity, provide a unique opportunity to evaluate vulnerability of a typical semi-arid mountain ecosystem to climate change. Therefore, we evaluated spatio-temporal variation of annual net primary productivity (NPP) of Ilgaz Mountains predicted by Moderate Resolution Imaging Spectroradiometer (MODIS), as affected by topography and climate between 2000 and 2010. The annual MODIS NPP ranged from 500 to 912 g .m -2 .y -1 . Elevation, slope aspect, and vegetation type were significantly correlated with MODIS NPP The MODIS NPP was highly sensitive to droughts, and the mean MODIS NPP generally decreased across the study pixels in the study period. The response of MODIS


INTRODUCTION
Net primary productivity (NPP) is one of the key indicators of the ecosystem productivity (Zhou et al., 2007;Zhengchao et al., 2011;Liu et al., 2013;Reeves et al., 2014;Reeves et al., 2014;Li and Pan, 2018), and it has been widely used in evaluating vegetation response to climate (Tang et al., 2010;Chen et al., 2012;Fang et al., 2013;Fang et al., 2016) (Zhang et al., 2014;Hao et al., 2016;Erşahin et al., 2016;Berner et al., 2017;Kumar et al., 2018;Xiao et al., 2019;Li and Pan, 2018). NPP is the net carbon sequestered by terrestrial vegetation (Zhou et al., 2007;Fang et al., 2016;; it is the difference between photosynthesis production and respiration consumption (Xiao et al., 2018), and it quantifies atmospheric carbon fixed and accumulated as biomass by plants (Zhao and Running, 2010). NPP represents the growth status of an ecosystem and it is an index for terrestrial ecosystems (Zhou et al., 2007;. Therefore, NPP is commonly used as an indicator to evaluate ecosystem response to climate change (Kumar et al., 2018;. The effect of climate factors on NPP has become a hot research topic (Hao et al., 2016).
Similarly to many other terrestrial ecosystems, forest productivity is vulnerable to drought and temperature extremes as reported in many studies (Zhou et al., 2018 and literatures therein). The forest vulnerability to climate change is mainly adjudged by changes in the phonological characteristics, tree line shifts, distribution of forest types and productivity (Kumar et al., 2018 and references therein). The NPP has been used to evaluate response of forest ecosystem to climate in many studies (e.g., Yong et al., 2007;Gao et al., 2013;Zhang et al., 2014;Fang et al., 2013;Erşahin et al., 2016;Liu et al., 2018;Gang et al., 2019;Zhang et al., 2019). Zhou et al. (2018) reported a substantial decrease in forest NPP in southeastern China in 2009droughts (from August 2009to April 2010. Topography may alter climate change impacts on vegetation (Daly et al., 2010). Evidences suggest that topography has an important control on climate change impact on local vegetation (Daly et al., 2010). Drought negatively impacts water and carbon utilization of vegetation particularly in water-limited localities on a landscape (Gang et al., 2019). Also, local climate regimes lead to differences in species phenology, distribution and diversity (Rodrigo, 2000;Chung et al., 2006). Studies (Lundquist and Cayan 2007;Ashcroft et al., 2009) showed that climate change resulted in complex spatial variations on the ground due to differences in topography.
Ilgaz Mountains, in which our study area is located, are located in a transition zone from Dry-subhumid/ Semiarid Continental Central Anatolian to Mid-latitude Humid Temperate Coastal Black Sea climate. Ilgaz Mountains National Park with its complex topography provides a unique advantage to study the response of vegetation productivity to climate and topography. The study area is highly rich in vegetation; it has 170 plant species, 18 of which are endemic. It is one of the remarkable national parks of Turkey (Öner & Abay, 2005). Quantum and spatial and temporal variation of NPP and its relation to topography, climate, and vegetation may provide a thorough understanding of vegetation vulnerability and its adaptability under changing climate.
The objective of this study was to evaluate the response of vegetation of the study area to climate between 2000 and 2010. The novelty of this study relies on the consideration of both topography and climate variables in evaluating the response of NPP to changing climate and on the use of quality control charts (QCCs) to evaluate the negative impact of summer droughts on vegetation productivity in the study area. The specific objectives were to 1) evaluate temporal trend and pattern in MODIS NPP between 2000 and 2010, 2) evaluate topography-plant type-NPP spatial relations, and 3) analyze climate biological productivity feedbacks as affected by topography and vegetation type in the Ilgaz Mountains.

Study Site
The study area is located in Ilgaz Mountains (Figure1), which are located in a transitional climate from dry sub-humid semi-arid continental Anatolian type to Mid-latitude humid temperate coastal Black Sea type according to Iyigun et al. (2013). We used MODISpredicted NPP (MODIS NPP) at 24 contiguous pixels (a typical representative area for the mountains), each with 1 km 2 , from 2000 to 2010 (11 years). Therefore, the study area was 24 km 2 in the surface area ( Figure 1). The MODIS NPP data were available for only 24 contiguous pixels and for only period of 2000-2010. Therefore, we had to use those 24 pixels and study period of 11 years based on the NPP data availability. The study area is a complex terrain with varying altitudes and orientation of hillslopes that may result in diverse impacts of climate change on ecosystem productivity. The study pixels highly differ in distribution of elevation, slope aspect (Table 1) and vegetation type and density (Table 2)   by north-facing slopes, while no north-facing slope is found at the P16, P21, P23, and P24. Elevation in the study area ranges from 1434 and 2065 above sea level and precipitation increases, and temperature decreases orographically with increasing elevation. This nonuniform distribution of topography and climate variables results in differences in distribution of vegetation type and cover density in the study area (Table 2). In general, elevations > 1500 m may be characterized as subalpine and those < 1500 montane forests according to Nakawatase and Peterson (2006). Elevations > 1800 m are mainly dominated by mountain pastures and Pinus Nigra and those <1800 m are dominated by pure stands of Abies nordmanniana subsp. equi-trajani, Pinus nigra, and/or their mixtures. Soil parent material comprises sandstone, siltstone, clay stone, sandy limestone, clay limestone, pudding stone, and mudstone and are relatively well-developed to a depth of 70 cm or deeper at slightly and moderately sloping localities where approximate percent slope is between 3 and 10, while weakly developed at depths of 50 cm or shallower at moderate to steep slopes where approximate percent slope is between 11 and 25. In some steeply sloping localities where approximate percent slope is > 30, the soil profile becomes highly shallow (30 cm or shallower) with a thin A horizon overlying C horizon or sprolite. Soil texture is mainly fine, and it changes between clay and sandy clay loam.

NPP Data and Climate Variables
We used remotely sensed NPP (g . m -2 . y -1 ) data of U.S. National Aeronautics and Space Administration (NASA) Earth Observing System (EOS) of Moderate Resolution Imaging Spectroradiometer (MODIS) from 2000 to 2010. MOD17A3 is a vegetation NPP of the terrestrial ecosystems calculated by MODIS Terrestrial Research Group with the BIOME-BGC model (H. Liu et al., 2018). NPP data were converted to projected coordinate system by MRT tool. Mean, maximum and minimum elevation (above sea level) were obtained from the Shuttle Radar Topography Mission (SRTM). Crown closure map is used to determine plant type and coverage in the pixels.
The monthly values for temperature and precipitation were obtained from the two nearest climate stations, one is located at 950 m, approximately 20 km to the study area, and another is located at 1980 m, in the study area. The monthly precipitation and temperature values were interpolated across the area by normal distance weighting.

Evaluation of Climate extremes-MODIS NPP relations by Quality Control Charts
Impact of drought on MODIS NPP was evaluated by quality control charts (QCC). Typically, a QCC includes three lines: a center line (representing grand mean of yearly MODIS NPPs between 2000 and 2010), an upper limit control line, and a lower limit control line. If a sample mean falls outside either the upper or lower control line, the process is judged to be out of control, suggesting that ecosystem productivity has shifted (Ott, 1993). This method has been used successfully by Erşahin et al. (2016) to evaluate response of some Anatolian forests to climate extremes between 2000 and 2010.

Evaluation of MODIS NPP Trend
Significance of change trend of MODIS NPP was evaluated using the following equation 1 (Wang et CERNE BILGILI, et al al., 2019). Where, x i and y i are the years from 2000 to 2010 (from 1 to n, n = 11) and the NPP for each year, respectively. A positive value of r indicates an increasing NPP against time increase, while a negative one indicates decreasing NPP within the study period. A r was deemed significant in the case corresponding P-value was ≤ 0.05. We used the ordinary slope (S) of least square regression to evaluate significance of the trend of MODIS NPP by Equation 2 (Wang et al., 2019). Where, x i represents years from 2000 to 2010 (1 to n, n = 11) and y i represents the NPP for the year i. A positive value of S indicates an increasing trend of NPP, while a negative one indicates a decreasing trend (Wang et al., 2019). We tested significance of the calculated slope. The null hypothesis was rejected when P-value ≤ 0.05.

Evaluation of Growth-Climate Relations
Spearman's correlation coefficients were calculated between NPP and climate variables of mean, maximum and minimum monthly temperatures and monthly total precipitation to evaluate their associations to NPP for the length of the study period. Significant correlations were attributed to P < 0.05. Differences of means across pixel-averaged and year-averaged values were evaluated by analysis of variance (ANOVA) and means were grouped using multiple comparison test of Duncan (P < 0.05). Log-transformed values of NPP were used in testing differences of means in pixel-averaged NPP values due to that the values were strongly negatively skewed according to Webster (2001) who suggests log 10 -normal transformation of data with an absolute value < 1.0.

Descriptive Statistics of MODIS NPP
The 11-year pixel-averaged values for MODIS NPP ranged from 500 g . m -2 . yr -1 at pixel 19 (P19) in 2007 to 912 g . m -2 . yr -1 at P24 in 2001 (Tables 3 and 4). The greatest variability of MODIS NPP occurred at P13 and the lowest at P5 and P6 (Table 3). The MODIS NPPvalues, except for the year of 2010, were positively [ 1] [2] 11-year trend of annual MODIS NPP was declining in majority of pixels, as slope (S) and correlation coefficients (r) for those pixels evidenced. The slope of the temporal trend is significantly negative in approximately 35% of the study area, indicating a significant continuous degradation of vegetation productivity in the study area during the study period.
Relationship between MODIS NPP, Topography and Vegetation Type Table 6 shows correlation coefficients among MODIS NPP, topography, and vegetation. The MODIS NPP was significantly negatively correlated with mean altitude above sea level, southeast, and northwest aspects; and significantly positively correlated with southwest and west aspects ( Table 6). The MODIS NPP was correlated with vegetation composition (Table 6). For example, Pinus nigra (Pn) and grasslands were negatively and Abies nordmanniana subsp. equi-trajani + Pinus sylvesters (A+Ps) was positively correlated with MODIS NPP. The vegetation composition was correlated with topographic variables; Pn tended to exist on south-facing slopes, while A tended to exist on southwest-and west-facing slopes, and A+Ps at north-and northeast-facing slopes (Table  6). Table 6 further depicts that Pn and grasslands were tended to coexist contrary to A+Pn and grasslands, and similarly to Pn and grasslands, Ps and A were tended to coexist, too.

DISCUSSION
Our pixel-specific means of annual MODIS NPP between 2000 and 2010 ranged from 500 to 912 g m -2 . yr -1 . Pixel 6 (P6) had the greatest and P19 had the lowest mean MODIS NPP (Table 3). The CV of MODIS NPP was below 10% for all pixels, which is considered low according to Mulla and McBratney (2001). Values for skewness indicated that only one pixel (P11) has a strongly right skewed and 4 pixels have moderately skewed MODIS NPP distribution, while rest of the 24 pixels possessed slightly skewed distribuiton of MODIS NPP according to Webster (2001).
However, the spatial variability of the mean annual MODIS NPP was high as significantly different mean values of MODIS NPP for pixels indicated (Table 3). Our values for annual MODIS NPP were similar to most of those reported in previous studies. For example, Chhabra and Dadhwal (2004) reported mean NPP of 666 g . m -2 . yr -1 for India for the period from June 1998 to May 1999 and Bala et al. (2013) reported mean NPP of 800 g . m -2 . yr -1 for Indian Himalayan forests. Li and Pan (2018) reported kurtotic and strongly negatively skewed according to Webster (2001), who suggested that a distribution with a skewness > absolute 1.0 is strongly skewed. These strong negative skewness-values indicate the presence of lowvalued outliers at some localities. Pixel-averaged values (for 24 pixels) of MODIS NPP ranged from 535.9 g . m -2 . yr -1 at P19 to 841.3 g . m -2 . yr -1 at P6 (Table 4). Also, the P10, P13, and P14 were statistically uniform, while they were significantly different from the rest of the pixels. The MODIS NPP was less skewed and kurtotic spatially than temporally (Tables 3 and 4). Each of P17, P19, and P20 were significantly different from the rest of the 23 pixels.

Climate MODIS NPP relations
May, April, and June mean monthly temperatures (MMT) had the most drastic influence on the MODIS NPP in the study area (Table 5). MODIS NPP was negatively correlated with May MMT at 75% of the study area (18 of 24 pixels) and positively with April MMT at 55% study area. Monthly total precipitation (MTP) of January negatively affected MODIS NPP (Table 5) in the entire study area although most of the correlations were not significant ( Table 5). The study area was highly spatially variable in the response of MODIS NPP to monthly variables of climate. For example, MODIS NPP was correlated negatively with January total precipitation (TP) at P5, P9, and P17; with October TP at P14 and P17; and with August TP at P13. In general, correlation between MODIS NPP and monthly climate variables were consistent in direction, but not in extent. For example, January TP was correlated negatively with MODIS NPP at 21 of 24 pixels, while the correlation coefficient highly variable, ranging from -0.01 to -0.72 (in absolute value).
Our results of quality control charts (QCC) suggested that MODIS NPP was strongly affected by 2007 drought (Figure 2) as shown by that the annual NPP for 2007 fell below lower control limit at approximately 92% of the study area. Figure 2 also shows that annual MODIS NPP for many of the pixels fell blow lower control limit in 2010. In addition, some localities experienced a production shift in 2005 and some others in 2002 and 2003, which may be attributed to differences in topography and plant cover. For example, the MODIS NPP fell below lower control limit at the P14 in 2002, 2007, and 2010. The P14 is dominated by northwest-and north-facing slopes, and the northwest-facing slopes are significantly negatively correlated with MODIS NPP. Also, the northwest aspects are positively correlated with Pinus sylvestris (Ps). Similarly, vegetation at P18 is highly vulnerable to climate extremes, contrary to that at P5. The Figure 2 also depicts that the that the overall mean annual NPP ranged from 789.14 to 829.15 g . m -2 . y -1 in a large area comprising forest steppe, typical steppe, desert steppe, and desert. Liu et al. (2013) reported that their simulated annual NPP showed distinct spatial patterns from the south to the north due to combined effect of land cover and climate conditions and that the simulated NPP was highly vegetation type-specific.

Topography, Vegetation, Climate, and MODIS NPP Interrelations
Annual MODIS NPP was positively correlated (P < 0.05) with April mean temperature in 45% (11 of 24 pixels) and with April maximum temperature in 12% (3 out of 24 pixels) of the study area (Table 4). Therefore, April temperature was highly important, positively impacting the NPP in the study area (Table 5), while reason behind this significant impact on MODIS NPP is not clear. Fang et al. (2016) reported significant controls of April minimum temperature and June and July precipitations on the NPP of Pinus koraiensis in Changbai Mountains of China, and they concluded that effect of climate change on NPP was complex and that elevated monthly low temperatures increased NPP, while monthly high temperatures had the negative effect on NPP. Similarly, Tang et al. (2010) reported positive relationships between forest NPP and variation of April and May monthly temperatures, which they attributed to lengthening of the growing season in spring. On the other hand, the Table 5 depicts that May temperature and June temperatures and total precipitation had a considerable negative impact on MODIS NPP in the study area.
Slope aspect was a significant variable, mediating the vegetation-climate relation. For example, P12 responded to no monthly climate variables, while P13 responded significantly negatively to May maximum and minimum temperatures, June minimum temperature, and August minimum and maximum temperatures; and positively to August total precipitation. Plant composition and vegetation density were similar, while distribution of aspects was highly different between those two pixels. Therefore, this dissimilarity in vegetation response to monthly climate variables between the two pixels may be attributed to difference in their topography. Spatial distribution of vegetation type was significantly controlled by slope aspect and elevation (Table 6). Table 6 shows that Pn tended to exist on the southeast-, south-, and southwest-facing slopes and Abies nordmanniana subsp. equi-trajani (A) tended to exist on southwest-and west-facing slopes.
Similarly to slope aspect, elevation above sea level had a significant control on MODIS NPP in the study area (Table 6). Mountain pastures and Pn were positively correlated with elevation; while both of them negatively correlated with MODIS NPP, which indicated a significant mediating effect of elevation on MODIS NPP via its interrelation between slope aspect, vegetation type, and climate in the study area. Table 6 shows that southeast, south, and northwest aspects tended to increase, and northeast, southwest, and west aspects tended to decrease with increasing elevation. This multiple effect of slope aspect, elevation, climate, and vegetation type on MODIS NPP would be the key factor controlling spatio-temporal variation of MODIS NPP in the study area. Precipitation typically increases and air temperature decreases with elevation in the study area. However, the net result of temperature-precipitationbiological productivity interactions was negative as correlation between mean elevation and MODIS NPP shows (Table 6).
NPP showed a decreasing temporal trend (Figure 2), suggesting that the vegetation is vulnerable to climate in the study area. However, the pixels differ in extent of vegetation vulnerability to climate as differences in slopes of the pixel-specific MODIS NPP temporal trends depict. The term "vulnerability" has been considered as the degree to which the ecosystem is susceptible to and is unable to cope with adverse effects (Kumar et al., 2018).
Mixed results have been reported by others on response of NPP temporal trends to climate variability across scales. For example, Zhang et al. (2014) reported that global NPP showed a decreasing trend from 2000 to 2009, which was strongly controlled by temperature and precipitation, and that the trend was locationspecific. Similarly, Fang et al. (2013) found that trends of NPP closely followed the trend of precipitation from 2000 to 2010 and that NPP was correlated negatively with air temperature. Zhang et al. (2014) showed that climate change had a negative impact on NPP at 31.9% and positive impact at 40% of the Yangtze River Basin of China from 2001 to 2010 and that the negative effect was mainly caused by water stress. Similarly, Fang et al. (2013) showed significantly different trends of NPP between different locations in Xijiang region of China between 2000 and 2010. Yong et al. (2007) found that the correlations between vegetation growth and climate variables were highly location-specific over the regions of Yunnan, Hainan, Taiwan and southeast coastal areas in China. Bala et al. (2013) reported increasing trend of remotely sensed NPP over India for the period between S is slope for temporal trend of MODIS NPP and r is the correlation coefficient between MODIS NPP and time (years). * and **: Significant at 0.05 and 0.01 level of significance, respectively.   (Table 1), and except P18, A and A+Ps are common vegetation cover to those pixels ( Table 2).
The mean annual MODIS NPP fell below lower control limit in several years during the study period indicating a production shift, while the production shift is highly pixel-specific (Figure 2). Severest decrease occurred in the worst year of 2007, in which a serious summer drought was reported (FAO, 2017). The year 2007 stood out from the years within 2000 to 2010 with MODIS NPP fell below the lower quality limit in the 96% Our QCCs showed that vegetation productivity shifted in some pixels in the study area in some years others than 2007. For example, MODIS NPP fell below lower control limit at 7 pixels in 2010, 5 pixels in 2005, 3 pixels in 2002, and 1 pixel in each of 2003 and 2009. There could be multiple reasons for spatially variable sensitivity of MODIS NPP to variation climate attributes. Differences in site-specific interactions among factors of topography, vegetation type, soil, and microclimate can be the principal factors behind the spatially varied response of vegetation productivity to climate variation in the study area. For example, P8 was significantly hit in 2002,2005,2007, and 2010, while P6 was hit only in 2007. P8 and P6 are highly similar in elevation and vegetation type; despite they are different in distributions of southeast-, south-, and southwestfacing slopes, suggesting that large differences in drought vulnerability between these two pixels may be attributed to dissimilarities in slope aspect-vegetation drought sensitivity relations between those two pixels. Similarly, P19 and P20 stood out from the rest of the study area by their distinct temporal pattern of MODIS NPP. P19 and P20 were different from the rest of the pixels in that the slopes for the temporal trend of MODIS NPP for those pixels ( Figure 2) were considerably lower than the rest of the 24 pixels. Those two pixels comprise considerable acreage of mountain pastures (Table 2). Elevation (mean, maximum, and minimum elevation) is greatest for P19 and P20 (Table 1). Therefore, again, topography induced differences in vegetation could be the major cause of the highly different temporal response of the MODIS NPP of those two pixels from the rest of the study area.
We believe that pixel-specific differences in plant type and composition, topographic features, and micro-climate topography relation would be important determinants of the pixel-specific response of MODIS NPP to the climate extremes. Plants at different elevation may differ in their growth responses, when they are under stress (Wei et al., 2018). Although not considered in this study, pixel-specific differences in soil type pattern may have a tremendous influence on NPP response to extreme climate events as suggested elsewhere (Wei et al., 2018), and this has a merit of a comprehensive evaluation.
Our comparisons suggest that MODIS NPP can provide useful predcits of NPP which enabled examination of spatiotemporal variation at landscape scale. Numerous of regional and global scale studies of NPP have shown the reliability of the MODIS NPP as reported by . Turner et al. (2006) reported a relatively strong agreement between MODIS NPP and ground measured values of NPP across nine sites covering a wide range of vegetation types and they concluded that MODIS NPP products were responsive to general trends in the magnitude of the NPP associated with local climate and land use, while tended to overestimate at low productivity sites and underestimate at high productivity sites. Similarly, Chen et al. (2012) reported a significantly strong relationship between MODIS NPP and field measured NPP although MODIS systematically overestimated the NPP. Tang et al. (2010) reported a correlation coefficient of 0.86 between their modelled NPP and MODIS NPP and they concluded that the two data sets were highly similar in forests of New England. Although MODIS NPP algorithm might overestimate forest NPP, it is a powerful tool for monitoring the terrestrial NPP due to its spatial coverage and temporal continuity (Tang et al., 2010). Studies have shown that MODIS17A3 data can be applied with accuracy in China (Li and Pan, 2018).
In this study, we deemed that vegetation growth is principally controlled by climate, simply ignoring the human influence, since most of the study area is a national park. However, many different factors, such as animal activities, pests, and diseases would still influence vegetation growth in the study area as noted by Gang et al. (2019). Also, the accuracy of the results was constrained by resolution of the remotely sensed data. As the MODIS NPP algoritm excludes the effect of soil water and nutrient stress on carbon fixation , it may overestimate the NPP as stated above. Other factors, such as forest type, stand age, and regeneration scheme may play important roles in ecosystem productivity-water use relations as noted elsewhere (W. .

CONCLUSIONS
MODIS NPP decreased across the study area during the period considered in our analyses, and the temporal trend of MODIS NPP was significantly negative in 38% of the study area. The vegetation productivity response to climate was highly site-and time-specific in the study period. The NPP was tidily related to climate variables such as monthly precipitation and temperature, as well as to variables of elevation, slope aspect, and vegetation type. Drought response of MODIS NPP was highly site and time-specific, too, and this was attributed to spatio-temporal variations in vegetation response to climate, which was mediated by topography to a large extent. Our results contributed to a comprehensive understanding of the NPP response mechanisms to climate change in a reserved mountainous area between 2000 and 2010. The information gained in this study can be used to develop adaptive management strategies to prepare for the impact of future climate on those forests. However, given the discrepancies in forest NPP between different approaches, continued efforts at inter-comparison and cross-validation are needed for reducing the uncertainties in the evaluation of plant-topography-climate feedbacks and interactions. This study has some limitations including no measured NPP data to validate NPP MODIS data in the study area, which deserves further efforts for improving the accuracy of the analysis. In this regard, controlled experiments and long-term observations can help in-depth understanding of vegetation responses to the changing climate.