Species distribution modeling of red hawthorn (Crataegus monogyna Jacq.) in response to climate change

Forests have been used for basic needs, such as nutrition, shelter and protection since the beginning of human history. Non­Wood Forest Products (NWFPs) play a significant role in these areas of use. Türkiye hosts many plants with different ecological requirements owing to its geographical location and climatic characteristics (Özdemir and Özkan, 2016). While some of these species are widely distributed, others are locally distributed in limited areas. Studies indicate that species distribution and diversity are changing and may continue to change over time (Mert and Acarer, 2021; Almeida et al., 2022). It is believed that the main reason for this is climate change. This is because anomalies observed due to climate change can affect the habitat characteristics of the species. This process forces species to either adapt or relocate, which can threaten their continuity or cause significant changes in their distribution (El Gendy et al., 2023; Hosseini et al., 2024; Tekeş, 2024). Moreover, climate change is not the only factor that drives species to change. The increasing human demand for natural resources has become one of the reasons as well (Tekeş and Cürebal, 2017; Almeida et al., 2022). Considering their current use, NWFPs are one of the most prominent examples. At this point, in order to produce solutions, the most effective tool that can be used to reveal how species can react to climate change is species distribution modelling (El Gendy et al., 2023; Hao et al., 2024; Hosseini et al., 2024).

Species distribution models are inherently compatible with climate models, allowing for their integration into future climate scenarios and enabling the prediction of potential shifts in species distribution (Özkan, 2014; Zenbilci et al., 2024). In this context, the maximum entropy (MaxEnt) method is one of the most widely used species distribution modelling methods in the literature. MaxEnt is preferred due to its ability to work with non­existent data and produce highly accurate results with a small number of observation data (Özdemir et al., 2020; Hussein and Workeneh, 2021; Karakaya and Yücel, 2021; Acarer, 2024a). Additionally, integrating this method with the free MaxEnt 3.4.4 software package and open­source software like R and RStudio provides researchers with an accessible and practical solution (Phillips et al., 2004; Wisz et al., 2008; Phillips et al., 2009; Süel, 2014; Özdemir, 2022; Çıvğa et al., 2024).

In Türkiye, many species such as hawthorn, thyme, stone pine, sage, bay leaf, cyprus turpentine, Cistus, and rosehip stand out in terms of quantity (Özdemir and Özkan, 2016; OGM, 2018). Among these, hawthorn (Crataegus spp.) has a wide distribution in our country with 17 species and 34 taxa (Ulus et al., 2016). Red hawthorn (Crataegus monogyna Jacq.) is widely distributed in Türkiye (Meriçli and Melikoğlu, 2002; Özderin, 2014). Crataegus monogyna is a deciduous shrub that can grow up to 10 meters tall. This species is distributed in various habitats, from sea level to 2000 m above sea level. It is frequently found in hilly areas, maquis, oak thickets, mixed forests, and roadsides (Browicz, 1972). Crataegus monogyna flowers in April­June and its red­coloured fruits can remain on the tree until mid­winter. Therefore, it is an important food source for wild animals (Ayberk, 2004; Ünal & Arslan, 2019; Akkemik, 2020). Examination of the ecological requirements of hawthorn shows that it does not tolerate high temperatures, it can withstand temperatures as low as ­18ºC, and it is resistant to sea winds and air pollution (Genç, 2007; PFAF, 2023).

Hawthorn leaves, flowers and fruits have a wide range of uses in food and medicine. The shoots, roots, bark and fruits of hawthorn species are used in the treatment of various diseases, such as digestive system and heart diseases, gout, depression, kidney stones and hypertension. In addition, teas made from dried hawthorn fruits and flowers are widely used for the treatment of conditions such as cough, throat inflammation, and oedema (Vatansever, 2016). Crataegus monogyna is also preferred as an ornamental plant due to the colour of its fruits. Crataegus monogyna has the highest antioxidant capacity among the hawthorn species (Özyürek et al., 2012). In recent years, C. monogyna has stood out as a preferred NWFPs in arid, and semi­arid cold regions and has provided an important source of economic income (Bayar and Deligöz, 2016). Due to these characteristics, C. monogyna is the most common and frequently used species among other hawthorn species (Özderin et al.,2015; Bayar and Deligöz, 2016).

Studies on C. monogyna mostly focused on its chemical composition (Bahorun et al., 1994; Bahorun et al., 2003; Bernatoniene et al., 2008; Barros et al., 2011; Özyürek et al., 2012; Keser et al., 2014; Ruiz­Rodriguez et al., 2014; Tahirovic and Basic, 2014; Nabavi et al., 2015; Abuashwashi et al., 2016; Belabdelli et al., 2022), botanical characteristics (Martinelli et al., 2021) and genotype characteristics (Yılmaz et al., 2010). Potential dispersion modelling (Karataş et al., 2019) and climate scenarios integrated into these studies (Radha and Khwarahm, 2022) are limited. In Türkiye, no studies have been conducted in this direction. Therefore, the aim of this study was to model the current and future (2100) potential distribution of C. monogyna, which is distributed in the Bozdaglar Mountains region of the Aegean region in Türkiye, using the MaxEnt method. These outputs will enable more effective and informed planning of the target species in terms of sustainable forestry and NWFPs (Karataş et al., 2019). This study has a unique value in terms of the method to be used and the parameters evaluated and is one of the first studies in this field in Türkiye. In addition, this study serves as a guide for the identification and modelling of important NWFP species in each region.

MATERIAL AND METHODS

Study Area – Područje istraživanja

The study area, the Bozdaglar Mountains, is located in Manisa­Izmir provinces in the Aegean region of Türkiye and covers an area of approximately 259,000 ha (Figure 1). It is surrounded by the Gediz Plain in the north, the Küçük Menderes Plain in the south, Nif (Kemalpaşa) Mountain in the west, and the Alaşehir Plain in the east. The study area generally has a Mediterranean climate with hot and dry summers and mild and rainy winters, and it is located in the Mediterranean phytogeographic region (Atalay and Efe, 2015; Tekeş and Cürebal, 2019). Forest, maquis, Mediterranean mountain steppe, and subalpine vegetation are observed in the study area. Forest vegetation consists of red pine, larch, and chestnut trees (Günal, 1987; Bekat and Oflas, 1990).

The Bozdaglar Mountains is the largest mountain range in the Aegean region (Eken et al., 2006) and one of the richest mountains in the region in terms of vegetation. Approximately 750 taxa were recorded in the Bozdaglar Mountains, 104 of which are endemic. The Bozdaglar Mountains harbour about 40 rare taxa throughout the country, 8 of which are endemic species specific to the Bozdaglar Mountains. The Boz Dag Important Plant Area (KBA) is located in the centre of the Bozdaglar Mountains (Oflas and Bekat, 1988; Özhatay, et al. 2008). In addition to the richness of the ecosystem and plant diversity in the area, the presence of wildlife is also an important factor in determining the study area. The presence of the Bayındır­Ovacık Wildlife Development Area within the study area (Malkoç, 2017) increases the ecological value of the region. In addition, the fact that the Bozdaglar Mountains is among the areas that may be affected by climate change (GDM, 2015) is an important motivation for selecting this site. Also, mountainous ecosystems as areas where species displacement is most prominent is another reason for its selection. In mountainous ecosystems, which are characterised by high altitudes and low­temperature conditions, it is suggested that climate change may lead to more pronounced effects (Pauli et al., 2007; Pauli et al., 2012).

Figure 1. Location map of the study area Slika 1. Karta lokacije područja istraživanja

Preparation of environmental variables – Priprema varijabli okoline

Base maps of the study area were prepared using ArcMap interface of ArcGIS programme. Digital elevation model data were downloaded from the EarthData database (https://search.earthdata.nasa.gov/ search?q=C1711961296­LPCLOUD) as ‘tif ’ files, which were then combined and clipped to match the study area’s boundaries using the ArcMap software. In the next stage, slope, aspect and shading index (Mitchell, 1999; Burrough et al., 2015) maps were produced with the ‘Spatial Analyst Tools’ tool using the elevation base. Topographic position index (TPI) and solar illumination index maps were created with the ‘Topography Tools’ plug­in prepared by Jenness (2006). Ruggedness index (Jacek, 1997; Riley et al., 1999) map was created with the ‘Terrain Tools’ plug­in. Roughness index map (Jacek, 1997; Riley et al., 1999) was produced with ‘Geomorphometry and Gradient Metrics Tools’ plug­in. Solar radiation index map (Mitchell, 1999; Thuiller et al., 2003) was produced with ‘Spatial Analyst Tools’. After these base maps were produced, Radiation Index (RI) (Equation 2.1) (Moisen and Frescino, 2002; Aertsen et al., 2010), Aspect Favourability Index (AFI) (Equation 2.2) (McCune and Keon, 2002) (Figure 3.22) and Heat Index (HI) (Equation 2.3) (Parker, 1988) maps were created with the help of different equations in the literature with the ‘Raster Calculator’ tool in the ArcMap software. The equations for these indices are as follows:

RI = (1­cos((π/180) * (Q­30))) / 2 (2.1)

In the above equation, Q stands for the aspect value (Moisen and Frescino, 2002; Aertsen et al., 2010).

AFI = cos(Q_max­Q) + 1 (2.2)

In the above equation, Q is the aspect, while Q_max value of 202.5° is the value of the slope in radians (McCune and Keon, 2002).

HI = cos(A­202.5) x tan(Slope) (2.3)

Preparation of climate data – Priprema klimatskih podataka

For today and the year 2100, 19 different bioclimate datasets created by Hijmans et al. (2005) were obtained from the WorldClim database. The data for the year 2100 were downloaded from the UKESM1­0­LL (United Kingdom Earth System Model) (Sellar et al., 2019) projection with a resolution of 30 arc seconds (~ 1 km) for the period 20812100 for SSP1 2.6, SSP2 4.5, SSP3 7.0 and SSP5 8.5 scenarios. These maps downloaded at the world scale were cut in accordance with the dimensions of the study area and prepared for statistical analyses by geometric registration (Table 1).

After all the environmental base maps of topographic and climate variables were created, a grid network of the study area with a cell size of 30×30 m was prepared. All variables were applied to the grid network with the ‘Extract Multi Values to Points’ tool in the ArcMap software. Then, all variables were saved in ASCII format with an equal number of cells.

Field studies – Terenska istraživanja

In order to carry out field studies in the areas planned within the scope of the study, reconnaissance trips were organised first. During these reconnaissance trips, the areas where growing environment factors and plant communities change were visited and the variation in these areas was taken into consideration in the selection of sample plots. In addition, care was taken to select sample plots from areas where the species would not be exposed to human impact, such as grazing and recreation. Field studies were carried out in a total of 170 sample plots (Figure 1). Consequently, ‘presence’ data for C. monogyna were recorded in 49 sample plots. In addition, the coordinate values of the sample plots were also recorded during the inventory.

Statistical analyses and modelling process – Statističke analize i proces modeliranja

Before the modelling process, Pearson correlation analysis was applied to the RStudio software to prevent the multicollinearity problem that may arise from the high correlation between 12 topographic and 19 bioclimatic variables in the areas where ‘present’ data were recorded for C. monogyna (Acarer, 2024b; 2024c). The selected variables were included in the modelling process (Table 1).

Potential distribution modelling of C. monogyna for the present and future (2100) was carried out using MaxEnt 3.4.4 software. During the modelling process, a 10­fold cross­validity test was applied. The convergence threshold value was set to 10⁵, and the iteration limit was set to 5000 (Philips, 2005). The bootstrap method was used to assess the test results of the parameter estimates. These procedures were continued until a model that best explained the relationship between dependent and independent variables was obtained.

As in other types of distribution modelling methods, the accuracy of the model obtained from the MaxEnt method should be verified. In this context, the area under curve (AUC) value was used, which includes specificity and sensitivity indices that are the most widely used validity tests in the literature (Tekin et al., 2018). Accordingly, AUC values were categorised as ‘excellent’ (AUC > 0.90), ‘very good’ (0.90 > AUC > 0.81), ‘good’ (0.80 > AUC > 0.71), ‘low’ (0.70 > AUC > 0.61) and ‘unsuccessful’ (AUC < 0.60) (Swets, 1988). The dispersion model was then generalised to the SSP1 2.6, SSP2 4.5, SSP3 7.0 and SSP5 8.5 scenarios for the period 2081­2100 based on the UKESM1­0­LL projection. Finally, by comparing the maps obtained for today and the year 2100, areas most likely to be affected by climate change in the study area were identified.

RESULTS

Variable selection – Odabir varijable

Following the correlation analysis, variables with a p­value of ≥0.80 were excluded (Gülsoy et al., 2016; Karataş et al., 2019). Consequently, 10 environmental variables were selected for the modeling process: five bioclimatic variables (BIO1, BIO3, BIO7, BIO12, and BIO13) and five topographic variables (SLOPE, HILLSHADE, RI, HI, and TPI).

To determine the variables contributing to the model, elimination was performed based on their contribution rates. Then, the best distribution model obtained was extended throughout the region, and the potential distribution map of the present day was obtained.

Model performance – Učinak modela

During the modelling process, a 10­fold cross­validity test was applied. The convergence threshold value was set to 10⁵ and the iteration value was set to 5000. From the 33 variables, those with the lowest percentage contribution to the model were removed, and this process was continued until five variables remained. The AUC value of the final model was 0.802 (Figure 2a), indicating that the model falls into the ‘good’ category according to Swets’ (1988) classification. When we look at the jackknife test graph of the model, the variables shaping the model were determined to be BIO12, BIO7, HI, TPI, and BIO13 according to their contribution rates (%) (Figure 3).

In the jackknife test graph, the light­blue colour represents the change in the stability of the model when each variable is excluded. The downward trend of the bar graph on the training gain scale indicates the rate of change in the stability of the model when the relevant variable is excluded. The dark blue colour indicates the contribution of the relevant variables alone in the model. The red colour is a descriptive graph on the cumulative contribution of all variables to the model. In this graph, as one of the model evaluation criteria, it is expected that the dark blue line should not exceed the light blue line. This indicates that there may be inconsistencies in the model and its control is important. According to the results of the jackknife test, all the variables that make up the model are consistent. The marginal response curves for each variable are shown in Figure 4.

Table 1. Codes and descriptions of bioclimate and topographic variables

Tablica 1. Šifre i opisi bioklimatskih i topografskih varijabli

Figure 2. (a) ROC curve and AUC value of C. monogyna; (b) Plot of mean deficiency and estimated area of C. monogyna Duncan's test results for the number of roots in greenhouse media and phytohormones

Slika 2. (a) ROC krivulja i vrijednost AUC za C. monogyna; (b) Prikaz srednjeg nedostatka i procijenjene površine za C. monogyna

Figure 3. Jackknife test of the variables forming the model of C. monogyna

Slika 3. Jackknife test varijabli koje čine model C. monogyna

It was observed that BIO12 (annual precipitation), which is the variable that contributes the most to the model, had a positive relationship up to 720 mm and a negative relationship above it. In addition, the potential distribution areas of the species increased in the 700­800 mm precipitation range. When we look at the second variable, BIO13 (precipitation of wettest month), which is effective in the distribution of the species, it is seen that the potential distribution areas of the species increase in the range of 135155 mm. For the third variable, BIO7 (temperature annual range), which is effective in the distribution of the species, it was determined that the potential areas increased in the range of 28­31^oC. It was observed that the BIO13 and BIO7 variables had a positive relationship up to a certain value and then a negative relationship. The fourth variable, the temperature index value, was negatively affected when it reached 1.0. In other words, there is a negative relationship with values greater than 1.0. When TPI, which is the last variable affecting the distribution of the species, is examined, it is seen that its distribution is at the level of ­140 and ­50 values. A negative relationship was observed between TPI and the variable. Other variables included in the modelling process are not included here since they do not contribute to the formation of the model. As a result of the modelling process, a potential distribution map of the target species was produced using environmental variables (Figure 5). When the obtained map is examined, it is seen that the most suitable areas are distributed throughout the area, especially in the northeastern, eastern and southeastern parts, except for the summit parts where the elevation increases in the Bozdaglar Mountains. After obtaining the potential distribution map according to today’s climatic conditions, potential distribution maps were created in terms of four different climate scenarios according to the climatic conditions of 2081­2100 for the UKESM1­0­LL projection (Figure 6). In addition to these maps, the amounts of suitable and unsuitable areas for the species under both current and future climate change conditions were calculated and are shown in Figure 7. Accordingly, the suitable distribution area of the target species, which is 182.214 ha today, is predicted to decrease to 7.311 ha according to the worst­case scenario, SSP5 8.5. In parallel, it is estimated that the current 75,490 ha of unsuitable area will increase to 250,393 ha according to the SSP5 8.5 scenario (Figure 7).

Figure 4. Marginal response curves of model variables of C. monogyna

Slika 4. Krivulje graničnog odgovora varijabli modela C. monogyna

Figure 5. Potential distribution map of C. monogyna using MaxEnt method according to current climatic conditions

Slika 5. Karta potencijalne distribucije C. monogyna korištenjem MaxEnt metode prema trenutnim klimatskim uvjetima

Figure 6. Potential distribution maps of C. monogyna using MaxEnt method according to future (2081–2100) climatic conditions (a: UKESM10-LL SSP1 2.6, b: UKESM1-0-LL SSP2 4.5, c: UKESM1-0-LL SSP3 7.0, d: UKESM1-0-LL SSP5 8.5)

Slika 6. Karte potencijalne distribucije C. monogyna korištenjem MaxEnt metode prema budućim (2081–2100) klimatskim uvjetima (a: UKESM1-0-LL SSP1 2.6, b: UKESM1-0-LL SSP2 4.5, c: UKESM1-0LL SSP3 7.0, d: UKESM1-0-LL SSP5 8.5)

DISCUSSION

Climate change poses an increasing threat to the persistence of plant species within ecosystems. Analysis of the literature suggests that the distribution ranges of these species may shift or face extinction risks due to climate change. In this context, understanding how NWFP species—important for food, medicinal, and economic purposes—will be affected by climate change and how their distribution may shift is a critical research topic. Türkiye harbors significant NWFPs species, and their utilization is gradually expanding. Therefore, to ensure the sustainability of NWFPs species, it is essential to conduct studies on their potential distribution modeling and mapping under both current climate conditions and climate change scenarios. However, in Türkiye, limited research has been conducted on the distribution modeling of NWFPs species under current climate conditions (Karataş et al., 2019; Çıvğa, 2023) and climate change scenarios (IPCC; 2013; O’Neill vd., 2016; Carbonbrief, 2018; Akyol and Örücü, 2019; Akyol et al., 2020; Arslan et al., 2020; Akyol et al., 2023).

In this context, with this study, potential distribution modelling and mapping of the present and future (2100) status of C. monogyna, one of the important NWFPs species in Türkiye, was carried out. Crataegus monogyna is a medically important species whose leaves, flowers and fruits are used in the treatment of heart failure, as well as in helping to regulate low and high blood pressure, breaking down stored fats and cholesterol in the body (Altınterim, 2012;

Figure 7. Potential suitable and unsuitable areas of C. monogyna according to present and future climate scenarios

Slika 7. Potencijalna prikladna i neprikladna područja za C. monogyna prema sadašnjem i budućem klimatskom scenariju

Yener and Ay Ak, 2021). In addition, C. monogyna is consumed by wild animals and contributes to plant and wild animal diversity (Ayberk, 2004; Ünal ve Arslan, 2019). Due to these functions, the demand for this species is high and it is subjected to anthropogenic pressure. For this reason, the sustainability of the species within its natural distribution areas should be ensured. Therefore, it is important to determine the environmental conditions affecting its distribution and reveal its status under future climatic conditions.

In our study, which aims to present the potential distribution of the target species, as a result of the modelling, the AUC value of the model was determined as 0.802 and was classified as ‘good’. The key variables influencing the model were identified as BIO12, BIO7, HI, TPI, and BIO13, ranked by their contribution rates. Similar results were found in the study on the distribution of C. azarolus L. by Jafari et al. (2022) conducted in Iran. The researchers modelled the distribution of the target species with 7 climatic variables, 3 topographic variables, land use and soil variables using the MaxEnt method. It was determined that elevation, relative humidity and average annual rainfall variables were effective in the distribution of the species. Within the scope of this study, it is seen that BIO12, one of the variables forming the model obtained for C. monogyna species, is in common with the results of Jafari et al. (2022). A similar study was conducted by Radha and Khwarahm (2022) on the distribution of C. azarolus and C. monogyna species in Iraq. The researchers modelled the distribution of both hawthorn species using the MaxEnt method and found that annual mean temperature and annual precipitation variables were the factors affecting the distribution of both species. As one of the results of the present study, BIO12 is an important variable forming the model obtained for C. monogyna species, similar to the results of Radha and Khwarahm (2022). In addition, they predicted the suitable habitat of both hawthorn species will narrow in the future compared to their current situations. As seen, their study supports our findings, showing similar outcomes across different regions.

In this study, when the potential distribution map of the target species according to the current climate conditions was analysed, the potential suitable areas were observed throughout the region, especially in the eastern part, except for the summit parts of the mountain. When the maps for future climatic conditions are analysed, the target species in the SSP1 2.6 scenario has experienced contractions in the northeastern and southern parts of the site and the potential suitable areas have been relatively displaced. In the SSP2 4.5 scenario, the distribution of the target species increases, and the suitable areas decrease. In the SSP3 7.0 and SSP5 8.5 scenarios, the suitable areas of the species narrow down and remain only in the western part of the site, while the species disappears in other places. In summary, in the potential distribution maps made according to different climate scenarios, it is predicted that the future distribution of the target species will decrease, and the majority of the potential suitable areas will disappear. These results are consistent with the results by Naghipour et al. (2021) and Radha and Khwarahm (2022). It was also predicted that the distribution of these hawthorn species would decrease significantly in the future.

CONCLUSIONS

This study is considered as an important step in understanding the impact of climate change on C. monogyna.