Maps of plant hardiness zones are useful tools for determining the extreme limits for the survival of plants. Exploration of projected climate change effects on hardiness zones can help identify areas most affected by climate change. Such studies are important in areas with high risks related to climate change, such as the Mediterranean Sea region. The objectives of this study were to (i) map plant hardiness zones for Albania and (ii) assess the projected effects of climate scenarios on the distribution of hardiness zones. Hardiness zones were affected by IPCC AR5 climate scenarios. The most extreme hardiness zone (6a) disappeared while a new, warmer zone (10b) appeared, reflecting rising temperature trends during the cold season. The shifts in spatial distribution of hardiness zones may represent opportunities for introducing new species to Albanian agriculture and forestry; however, the introduction of new species would require further studies on the variability of plant hardiness zones at local scales.
Many environmental factors control the survival and distribution of plant species. The most important are often weather-related ones, such as minimum and maximum temperatures, frost frequency and duration, precipitation amount and distribution, cloud cover, and relative humidity (Ouellet and Sherk 1967; McKenney et al. 2007; Pliscoff et al. 2014). There are several ways to assess the survival and adaptability of plants to climatic fluctuations, including the use of plant hardiness zones (severity classes) based on absolute annual minimum temperatures (USDA model) and heat zones (American Horticultural Society 1997) based on the average number of days each year with temperatures above 30°C as well as broader climatic zone classification systems (Daly et al. 2008, 2012; Widrlechner et al. 2012; McKenney et al. 2014). However, the method using plant hardiness zones based on minimum temperatures is one of the most widely adopted in temperate regions because it reflects the long-term survival of perennial plants and the degree to which they can endure low winter temperatures (Widrlechner et al. 2012). Because of the simplicity of this method, which relies on just one climatic factor, the maps of plant hardiness zones are useful and practical (McKenney et al. 2007; Widrlechner et al. 2012). Hardiness zone maps based on the USDA model have been developed for many parts of the world, including Australia (Dawson 1991), China (Widrlechner 1997), Europe (Heinze and Schreiber 1984), Japan (Hayashi 1990), South Africa (Pienaar 1996), Ukraine (Widrlechner et al. 2001), and Turkey (Yalçın and Tolunay 2012).
Climatic constraints represent an important link between prevailing meteorological conditions, as measured at weather stations, and their specific relevance for land-use activities (Stone and Meinke 2006). A change in climatic constraints implies that new opportunities for and/or risks to land-use management could develop. Therefore, an exploration of climate change scenarios on plant hardiness zones could identify geographic areas where the range of land-use options might be expected to change and help farmers and land managers select the best options for adapting to the expected changes. For many years, most studies related to increased concentrations of atmospheric greenhouse gases have focused on impacts under average climate scenarios. During the last decade, however, changes in climatic variability and extreme weather events have received increased attention (Brown et al. 2008). Such studies are especially important in areas with a high risk of climate change, such as the Mediterranean region (Resco de Dios et al. 2007).
Albania is a small country, approximately 28 000 km2, but because of its mountainous, high-relief terrain, proximity to the Mediterranean Sea, and overall atmospheric circulations, it encompasses diverse climate conditions (Paxian et al. 2015). This climatic diversity provides unique opportunities for the development of diverse horticultural production systems. According to the climate classification of Albania developed by the Hydrometeorological Institute (Hydrometeorological Institute 1988; Jaho et al. 1988), there are 4 major climate zones and 13 subclimate zones. Although this is one of the most widely used climate classifications in Albania, its practical application is difficult for two major reasons. First, the interpretation of this classification for economic purposes, such as production agriculture and forestry, is difficult. Second, this classification is not compatible with other international classifications and so cannot facilitate transferability of horticultural species across borders. In addition, the existing classification is based on outdated climate data and does not reflect current and future trends under different projected climate scenarios (Moss et al. 2008, 2010). Several studies have been undertaken to better understand the Albanian climate for agricultural purposes (Teqja 1994; Kopali et al. 2012). Teqja (1994) elaborated a method to assess the variability of crop yields caused by climatic factors and made a classification of Albanian territory based on this indicator. Kopali et al. (2012) identified 26 homogenous climatic zones of Albanian territory based on historic data of monthly precipitation, temperatures, and relative humidity. Both authors used 1988 historic climatic data compilation from the Hydrometeorological Institute of Albania (Hydrometeorological Institute 1988; Jaho et al. 1988).
The objectives of this study were to (i) develop a classification of plant hardiness zones (severity classes) for Albania that is based on absolute annual minimum temperatures compatible with international classifications and (ii) assess the impact of future climate scenarios on the distribution and extent of plant hardiness zones and the distribution of some horticultural species.
2. Materials and methods
a. Study area characteristics
Albania is located in southeastern Europe and has a total area of 28 748 km2 (Fig. 1).
It has a diverse climate that varies within short distances because of its mountainous, high-relief terrain. The mean elevation is about 750 m but varies from 0 to 2700 m above sea level. The western coastal area is relatively flat and is dominated by a Mediterranean climate, with mild winters and hot summers. In this area, the mean winter temperature is about 7°C and the mean summer temperature is about 24°C but can range to more than 30°C, especially in the southwestern part. The mountainous areas, which make up close to two-thirds of the territory, are influenced by a more continental climate. Precipitation varies greatly between regions. The mean annual precipitation for the coastal areas ranges from 100 to 150 cm and falls mostly during fall and winter. The mountainous areas, especially the Albanian Alps, can receive as much as 250 cm of precipitation, with 20%–30% of it in the form of snow.
b. The analysis of climate impact on plant hardiness indicators
The absolute monthly minimum temperature during the 1970–85 period was recorded in February (−13.7°C). Also, the month of February had the highest number of days with temperatures below −12° and −13°C as compared with January, which was the second coldest month on record. Initially, an assessment was conducted to assess differences between calculated minimum mean February temperature and actual absolute minimum February temperature for 167 climate stations during the period 1970–85 (described in detail in section 2e). This was necessary to allow us to adjust WorldClim (http://www.worldclim.org/; Hijmans et al. 2005) minimum temperature data predictions for the 1950–2000 period and climate scenarios to the absolute minimum temperature. The adjustment allowed for a comparative analysis of the influence of climate scenarios on the plant hardiness indicators. Ideally, a direct comparison between actual annual minimum absolute temperatures and predicted absolute minimum temperatures from climate scenarios would have been used. However, only mean monthly minimum temperatures from the WorldClim (Hijmans et al. 2005) data were available. Also, the approach based on the current climate data does not imply that the same exact relationships between mean minimum and absolute minimum temperatures will be applicable for the future climate scenarios. Increased climate fluctuation and uncertainty (Deser et al. 2012a,b) will most likely affect the strength and variability of this relationship.
This analysis was followed by comparisons between different forcing scenarios from downscaled IPCC Fifth Assessment Report (AR5) phase 5 of the Coupled Model Intercomparison Project (CMIP5) data at 30-s resolution from the WorldClim data (IPCC 2013; Moss et al. 2008, 2010). The data provided four different scenarios of representative concentration pathways (RCPs) of greenhouse gas (GHG) (rcp26, rcp45, rcp60, and rcp85) for two time periods (Table 1): 2050 (average for 2041–60) and 2070 (average for 2061–80) based on forecasts (Fig. 2). The RCPs (rcp2.6, rcp4.5, rcp6, and rcp8.5) described four possible ranges of radiative forcing values in 2100 relative to preindustrial values [+2.6, +4.5, +6.0, and +8.5 W m−2 (Weyant et al. 2009)]. The dataset provided four climate parameters: (i) tn—monthly average minimum temperature (°C × 10), (ii) tx—monthly average maximum temperature (°C × 10), (iii) pr—monthly total precipitation (mm), and (iv) bc—bioclimatic variables (not used in this study).
c. Comparisons between the 1950–2000 period and future climate scenarios
Based on the elevation gradient and topography, the study area was divided into five major physiographic regions: (i) western low plains (0–125 m MSL), (ii) western midplains (126–250 m MSL), (iii) hills and foothills (251–500 m MSL), (iv) mountain foothills (501–1000 m MSL), and (v) mountains (>1001 m MSL). Based on the influence of low temperatures on plant community development and distribution, the mean annual temperature values (Widrlechner et al. 2012) were divided into severity classes or plant hardiness zones (Table 2).
Changes in spatial extent of various plant hardiness regions and trends over time under various climate change scenarios were assessed by using ArcMap software (ESRI 2010). Gridded maps of February monthly minimum temperatures from the WorldClim (Hijmans et al. 2005) data were downloaded and clipped to the extent of the study area. The WorldClim (Hijmans et al. 2005) data grids were adjusted by subtracting a predicted temperature value derived from the regression equation of the difference between the measured minimum mean February temperatures and absolute minimum February temperatures for 167 stations between 1970 and 1985, as described in the section 2e. The grids were reprojected to the European Albers Equal Area projection and classified based on the USDA plant hardiness zone classification system (Daly et al. 2012). The areal extent of plant hardiness zones was calculated for the 1950–2000 period and future climate scenarios (rcp26, rcp45, rcp60, and rcp70) for the 2041–60 and 2061–80 periods based on the predicted minimum absolute temperatures for February rather than mean temperature.
d. Future climate scenarios influence on areal extent of selected plant species
We selected seven horticultural plants (Pistacia vera, Acca sellowiana, Actinidia deliciosa, Musa acuminata, Persea americana, Ananas comosus, and Mangifera indica) that have been already introduced to the Albanian market. Their hardiness zones were determined based on work conducted by other researchers (Ferguson 1999; Herrera 1997; Gilman 2014; Gilman and Watson 1994, 2014a,b,c; Reiger 2016). In this analysis, we considered long-term benefits of introducing such species for attaining reasonable and consistent yield levels. We therefore excluded the extreme cold areas and limited our analysis to the potential plant hardiness zones above the lowest theoretical hardiness zone for these plant species. This allowed for a direct comparison of the lowest theoretical hardiness zone for these plant species with the hardiness zones based on current and projected future climate scenarios.
e. Statistical analysis
A simple regression analysis was conducted to assess the relationship between mean absolute minimum annual temperature and February mean temperature for 167 climate stations for the 1970–85 period. The analysis was also conducted by physiographic regions to assess their influence on the relationship. An additional simple regression analysis was conducted to assess the relationship between measured elevations at the climate stations and elevations derived from Shuttle Radar Topography Mission (SRTM) data (Jarvis et al. 2008). A “robust fit” analysis was conducted to test for outliers. The analysis revealed the presence of 9 outliers. However, the outliers were not excluded from the analysis to raise awareness about the variability in elevation data. Also, we could not make any assumptions about elevation values from stations as being more accurate than SRTM elevations values and vice versa. This facilitated the temperature predictions for plant hardiness zones generated from SRTM elevation grids (Jarvis et al. 2008) and WorldClim (Hijmans et al. 2005) temperature data grids. Daly et al. (2008) similarly predicted temperature and precipitation on 30-arc-s (800 m) grids for the conterminous United States using the Parameter-Elevation Regressions on Independent Slopes Model (PRISM). An analysis of variance (ANOVA) was conducted to compare absolute mean minimum temperature changes among climate stations, WorldClim (Hijmans et al. 2005), and IPCC AR5 climate scenarios within each physiographic region. Mean comparisons were made by using the LSMeans Tukey honestly significant difference (HSD) test. A pairwise t test was used to compare mean absolute minimum annual temperatures between WorldClim (Hijmans et al. 2005) and climate stations and between WorldClim (Hijmans et al. 2005) and IPCC AR5 climate scenarios. The null hypotheses were rejected at a significance level of 0.05.
3. Results and discussion
a. Measured mean absolute Tmin from climate stations versus mean Tmin from the WorldClim data
The WorldClim (Hijmans et al. 2005) data including climate scenarios provided only predicted mean minimum temperatures (Tmin). However, the severity classes of plant hardiness zones are based upon absolute annual minimum temperatures. To adjust the predicted mean Tmin from the WorldClim (Hijmans et al. 2005) data to the absolute Tmin, it was necessary to assess the relationship and the difference between measured mean Tmin and measured absolute Tmin from climate stations. The regression analysis showed a high degree of correlation (Fig. 3a) between measured mean Tmin and measured absolute Tmin from climate stations [R2 = 0.92; root-mean-square error (RMSE) = 1.3]. The degree of correlation varied by physiographic region. The mountain foothills region had the lowest R2 (0.65) followed by western low plains (R2 = 0.77). The RMSE varied from 0.9°C (western low plains and hills and foothills) to 1.7°C (mountain foothills). The regression equation showed an average difference of −7.3°C between mean Tmin and predicted absolute Tmin; and the difference varied between −6.7° and −7.6°C. On average, the differences between regions and average difference (−7.3°C) was only −0.3°C, thus the predicted Tmin from WoldClim (Hijmans et al. 2005) were adjusted to absolute Tmin by subtracting the average difference of −7.3°C. The use of predicted Tmin, from the WorldClim (Hijmans et al. 2005) data adjusted to absolute Tmin to determine the plant hardiness zones (severity classes), necessitated the assessment of the spatial trends in the difference. There were no observable trends in the difference with stations as shown by the overlap between measured absolute Tmin and mean Tmin (Fig. 3b). Also, the predicted absolute Tmin was highly correlated with measured absolute Tmin (R2 = 0.92; RMSE = 1.3) (Fig. 4a).
The WorldClim (Hijmans et al. 2005) data provide a grid that extends beyond the point data stations and use underlying assumptions on relationships between elevation and temperature. In this context, we first assessed the relationship between predicted elevation data from SRTM (Jarvis et al. 2008) and measured elevation data from the climate stations. The SRTM elevations were highly correlated with measured elevational data from the climate stations (R2 = 0.92; RMSE = 120 m) (Fig. 4b). However, the robust fit analysis indicated that there were nine outlier stations. The removal of the outliers improved the relationship (R2 = 0.96; RMSE = 81 m). The purpose of this analysis was to support the use of the SRTM elevation grid for extrapolating beyond climatic stations. We cannot prove nor assume that station elevation data are more accurate than SRTM, and the presence of the outliers highlights the variability of elevation as another source of uncertainty in our predictions. A second assessment examined the relationship between predicted absolute Tmin and SRTM elevation. The SRTM elevation explained 65% of the variability in predicted absolute Tmin (R2 = 0.65; RMSE = 1.9). Similar trends and relationships between elevation and Tmin have been found by other researchers. For example, Lookingbill and Urban (2003) found that Tmin is related to elevation (R2 = 0.58) in the Oregon western Cascades. Relationships among measured Tmin, absolute Tmin, and elevation allowed for comparisons of climate scenario impact on the spatial distribution and extent of plant hardiness zones.
As expected, there was uncertainty associated with these relationships. Although there was no clear evidence of biases or trends in Tmin errors related to spatial distribution and/or elevation, it should be noted that local variability due to physiographic characteristics does influence these relationships. Researchers have found that local- to landscape-scale physiographic characteristics do influence regional spatial patterns of temperature. For example, Dobrowski et al. (2009) found that between 70% and 80% of the variability in temperature patterns in the Lake Tahoe region in the United States was explained by regional patterns and the remaining percentage by landscape-scale physiographic characteristics. Interestingly, the RMSE associated with the difference between regional- and local-landscape-driven temperature patterns found by Dobrowski et al. (2009) was between 1.2° and 2.0°C, which is very similar to the values found in our study. The magnitude of RMSE is especially important in the context of severity classes based on the range values of plant hardiness zones used to classify the study area. The threshold value for separating each plant hardiness zone is 2.8°C, which is more than twice the overall 1.3°C RMSE in our study area. The actual application of plant hardiness zones needs to account for these uncertainties, especially for differences in the long-term variance of extreme low-temperature events (Widrlechner et al. 2012). Local applications especially may need to be adjusted based on landscape physiographic characteristics. Also, the grid resolution of the WorldClim data is approximately 350 m × 350 m, which warrants caution in interpreting the exact boundary of plant hardiness zones. This is particularly important in high-relief areas where elevation changes greatly within a short distance, such as the narrow valleys in the mountainous and southern coastal areas of Albania. Because of the uncertainties associated with absolute Tmin predictions and the relatively coarse resolution of predicted Tmin from the WorldClim data, the impact of climate scenarios on the plant hardiness zones should be assessed in terms of general trends over large areas.
b. Mean absolute Tmin between climate stations, WorldClim, and IPCC AR5 climate scenarios based on physiographic regions
The mean difference between climate stations and the WorldClim data was only 0.3°C and was not significant (p value = 0.06) (Table 3; Fig. 5). The differences between the WorldClim data and IPCC AR5 climate scenarios were all significant (p value < 0.0001), showing an increasing trend overall with increasing projected CO2 emissions (Figs. 6a,b). The mean absolute Tmin values difference between climate stations and the WorldClim data was 0.9°C but varied according to physiographic regions, from 0.6°C (mountain foothills and mountains) to 1.4°C (western low plains) (Table 4). The mean absolute Tmin were significantly different among regions and decreased with IPCC AR5 climate scenarios (Tables 5 and 6). While values for the western low plains and western midplains were not significantly different for all absolute Tmin data sources, the value differences were significant for the hills and foothills. As expected, the hills and foothills physiographic region had lower Tmin, followed by the mountain foothills and mountains, which had the lowest absolute Tmin values (Figs. 6a,b).
There was a wide range of variability in mean absolute Tmin, 1.7° to −10.2°C, between physiographic regions within each data source. However, within each physiographic region the differences between different data sources were smaller and varied between 0.6° and 4.3°C. The smallest difference, 0.6°C, was observed between rcp2650 and rcp8550, while the largest (a mean of 4.3°C) was between the WorldClim data (1950–2000 period) and the rcp8570 climate scenario (2061–80 period).
c. Severity classes based on the plant hardiness zones using predicted mean absolute Tmin from the WorldClim 1950–2000 data
As expected, the spatial distribution of plant hardiness zones was largely influenced by the relationship between temperature and elevation (Fig. 7). Daly et al. (2012) observed similar trends when developing a new plant hardiness zone map for the United States. The distribution of plant hardiness zones corresponded, to a great extent, with physiographic regions (Fig. 1) derived from the classification of elevations. Zones 9b and 10a were generally restricted to western Albania and coincided mostly with the western low plains and western midplains.
The mildest zones (9b and 10a) together covered approximately 22% of the country (Tables 7 and 8). Plant hardiness zones 9a, 8a, and 8b were the most dominant (15 600 km2) and covered approximately 56% of the country. Their distribution coincided mostly with the hills and foothills and mountain foothills. Plant hardiness zones 7a and 7b covered approximately 21% of the country, an area of 5 500 km2. These zones, along with plant hardiness zones 6a and 6b, which had a very limited extent (2.2% or 620 km2), coincided with the mountains physiographic region. A visual comparison of the plant hardiness map of Albania with the plant hardiness zone map of Europe developed by Heinze and Schreiber (1984) showed reasonable similarities with regard to the presence of severity classes and their spatial distribution. The same severity classes (6, 7, 8, 9, and 10) in the Heinze and Schreiber (1984) map were also present in the map developed from this study. However, as expected, the map from this study was more detailed spatially and with regard to subzones a and b within each severity zone, which are not present in Heinze and Schreiber (1984). The increased level of detail for plant hardiness zone and subzones would be more suitable for horticulturists developing strategies and practices for managing current and future plant species.
d. Potential climate change impacts on Albanian plant hardiness zones and the distribution of some horticultural crops
The spatial distribution and areal extent of plant hardiness zones were clearly affected by the various IPCC AR5 climate scenarios (Figs. 8 and 9) in diverse ways (Figs. 10a,b and 11a,b). Plant hardiness zone 10a experienced the greatest percent areal increase, from 4.8% to 15% on average, with a maximum of 18.4% (Tables 7 and 8). In addition, a plant hardiness zone (10b) new to the historic climate scenario, appeared and expanded to 12.6% of national coverage for the most aggressive climate scenario (rcp8570) (Table 8). In contrast, plant hardiness zones 6b, 7a, and 7b experienced the greatest percent areal decreases, from 2% to 0.1%, 6.3% to 0.6%, and 13.6% to 3.6%, respectively. Also, plant hardiness zone 6a disappeared after the rcp2650/2670 IPCC AR5 climate scenarios, indicative of rising temperatures during the cold season. The least fluctuating plant hardiness zones, in terms of areal extent, were zones 8b and 9b. Zone 8b only increased from 17% to 21.5%, while zone 9b varied slightly from 16.3% (rcp8570) to 18% (rcp2670).
Shifts in the spatial pattern of plant hardiness zones may present new risks for the long-term survival of some plant species but also provide opportunities for introducing new horticultural plants. In this study, we analyzed the potential area for the cultivation of six horticultural crops in Albania (Table 9). As expected, the future climate scenarios predicted increases in area for all six horticultural crops. For example, the area extent of kiwi (Actinidia deliciosa), which is already cultivated to a limited degree in Albania, could increase from 77.9 km2 to as much as 92 km2.
Also, the potential cultivation area of Pistacia vera and Acca sellowiana, which are currently cultivated in Albania, would increase in the future according to all climate change scenarios analyzed in this study (Table 9). The remaining three species, Musa acuminata, Persea americana, and Ananas comosus, could potentially be cultivated, especially under the most aggressive climate change scenarios. Although these projected shifts in climate represent opportunities for introducing new crops, studies will be needed to address the interactions of climate and management practices to mitigate future climate fluctuations and uncertainty as predicted by other studies (Tous and Ferguson 1996; Paxian et al. 2015), especially at local levels. This is critical for the introduction of species to areas that would be considered on the fringe of their range of adaptation. The areal extend of these species is based only on biological limits as dictated by the plant hardiness zones; however, other factors need to be considered. The potential cultivation area for Pistacia vera depends also on elevation, development stage, etc., (Herrera 1997). For example, kernel development is affected at elevations greater than 1 500 m, because of cool summers, while temperatures below −12.2°C could kill the trees, especially young ones (Herrera 1997). Economic factors may also play a role, thus the areal extend provided in this study are only projections that may vary depending on a combination of climate and economic factors.
We studied the climate impact on the changing geography of cold-season extreme temperatures and its implications for plant hardiness zones in Albania by comparing historical climatic data from WorldClim (http://www.worldclim.org/; Hijmans et al. 2005) for the 1950–2000 period with future climate scenarios based on downscaled IPCC AR5 phase 5 of the Coupled Model Intercomparison Project (CMIP5) data. The measured mean Tmin and absolute Tmin from climate stations were highly correlated, with a mean difference of −7.3°C. The distribution of the mean differences did not show any spatial bias and allowed for the assessment of the influence of future climate scenarios on the distribution of plant hardiness zones. The IPCC AR5 climate scenarios affected the distribution and extent of plant hardiness zones. The disappearance of the most severe plant hardiness zones and the appearance of new, warmer plant hardiness zones under rcp2650/2670 IPCC AR5 climate scenarios reflect rising temperatures during the cold season. However, the uncertainty associated with the use of predicted absolute Tmin from the relationship between mean Tmin and absolute Tmin from climate stations was about 1.3°C, or almost one-half of the 2.8°C threshold value used to separate plant severity zones. The uncertainty may further increase because of increased climatic fluctuations. Thus, the shifts in spatial distribution of plant hardiness zones need to be viewed critically. Plant hardiness zones related to climate change present opportunities for growing new horticultural crops; however, more studies are needed to assess the interactions between climate change and management practices for these new species, especially at local levels.
The authors thank Purdue University for funding the publication of this work (Grant 21030000). A portion of this research was funded by USDA Current Research Information System Project 6020-21310-009-00D Sustainable Small Farm and Organic Production Systems for Livestock and Agroforestry. The authors also thank Jenny Sutherland at the USDA–Natural Resources Conservation Service National Soil Survey Center for her meticulous editing and revisions. We are also grateful for the comments and suggestions made by many colleagues at the USDA–NRCS National Soil Survey Center. We thank the reviewers in particular for their detailed comments and suggestions that substantially improved the quality of this article.