Search tips
Search criteria 


Logo of eceLink to Publisher's site
Ecol Evol. 2017 February; 7(4): 1112–1124.
Published online 2017 January 23. doi:  10.1002/ece3.2766
PMCID: PMC5306017

Determining the factors affecting the distribution of Muscari latifolium, an endemic plant of Turkey, and a mapping species distribution model


Species distribution modeling was used to determine factors among the large predictor candidate data set that affect the distribution of Muscari latifolium, an endemic bulbous plant species of Turkey, to quantify the relative importance of each factor and make a potential spatial distribution map of M. latifolium. Models were built using the Boosted Regression Trees method based on 35 presence and 70 absence records obtained through field sampling in the Gönen Dam watershed area of the Kazdağı Mountains in West Anatolia. Large candidate variables of monthly and seasonal climate, fine‐scale land surface, and geologic and biotic variables were simplified using a BRT simplifying procedure. Analyses performed on these resources, direct and indirect variables showed that there were 14 main factors that influence the species’ distribution. Five of the 14 most important variables influencing the distribution of the species are bedrock type, Quercus cerris density, precipitation during the wettest month, Pinus nigra density, and northness. These variables account for approximately 60% of the relative importance for determining the distribution of the species. Prediction performance was assessed by 10 random subsample data sets and gave a maximum the area under a receiver operating characteristic curve (AUC) value of 0.93 and an average AUC value of 0.8. This study provides a significant contribution to the knowledge of the habitat requirements and ecological characteristics of this species. The distribution of this species is explained by a combination of biotic and abiotic factors. Hence, using biotic interaction and fine‐scale land surface variables in species distribution models improved the accuracy and precision of the model. The knowledge of the relationships between distribution patterns and environmental factors and biotic interaction of M. latifolium can help develop a management and conservation strategy for this species.

Keywords: abiotic factors, biotic factors, boosted regression modeling, bulbous plant, species distribution modeling

1. Introduction

Endemic species grow naturally in restricted geographic ranges, and specific habitats and are prone to become endangered under changing environmental conditions and other threats. They also have a great tendency to become extinct if they are both rare and endemic (Işık, 2011; Lomba et al., 2010; Marcer, Sáez, Molowny‐Horas, Pons, & Pino, 2013). Sustainable management practices and the preservation of endemic and rare plants are essential for the conservation of global biodiversity because these plants are important not only for local regions but also for global biodiversity. Therefore, endemic species are important targets for global conservation efforts (Myers, Mittermeier, Mittermeier, da Fonseca, & Kent, 2000).

Muscari is a genus of 46 species, distributed across Europe, Asia, and North Africa (Govaerts, Zonneveld, & Zona, 2015). Thirty‐three of these species occur naturally in Turkey (Davis & Stuart, 1984; Demirci, Özhatay, & Koçyiğit, 2013; Güner, 2012; Pirhan, Yıldırım, & Altıoğlu, 2014; Yıldırım, 2015). Muscari latifolium J. Kirk. (Asparagaceae) (Figure 1) is an endemic bulbous plant species of Turkey with a highly local distribution in western, inner western, and southwestern Turkey, Balıkesir–Çanakkale Kazdağı Mountains, Kütahya Murat Mountain, and Antalya Akseki at altitudes between 1,100 and 1,800 m in Pinus nigra J.F. Arnold and Pinus sylvestris L. forests (Davis & Stuart, 1984). Bulbs are solitary, ovoid, and 1.5–3 cm in diameter. Leaves are usually solitary, and two are found in rare cases, erect, broadly linear, lanceolate, 7–30 cm long, and 10–30 mm wide. Flowers are carried on a scape longer than the leaves. Inflorescences are racemes 1.5–6 cm long and consist of both fertile and sterile flowers. Sterile flowers are pale violet to light blue, 4–8 mm long, and located at the top of the raceme, whereas fertile flowers are dark violet to black, 5–6 mm long, and located at the bottom of the raceme. Fruit is a capsule 7–8 mm in size (Davis & Stuart, 1984). The species, which its flowering period is in April and May, are being used as ornamental plants (Bryan & Hort, 2002) and are usually propagated from seeds (Wraga & Placek, 2009). Muscari latifolium is easy to detect even outside the flowering season because of its broad leaves. It prefers lime and slightly acidic loamy soil with potassium, high phosphorus content, and rich organic matter (Hopa, Tümen, Sevindik, & Selvi, 2013).

Figure 1

Muscari latifolium. Photographed in Gönen Dam watershed, Turkey, April 2013

It is important to know the distribution, ecological traits, and population structure of endemic plant species to manage them in a sustainable manner and to develop effective conservation strategies for them. Determining the entire distribution area of a plant species is neither feasible nor realistic merely by navigating through the area without sampling. Furthermore, despite the fact that the area is well sampled, the organism will be present outside the sampling plots. However, species distribution models (SDMs) give us the ability to predict the distribution of the species across a landscape or within a certain time frame (Elith & Leathwick, 2009; Guisan & Thuiller, 2005; Peterson, 2006). SDMs are suitable tools for understanding the realized species distribution and for estimating the species’ potential distribution for endemic and rare species in well‐surveyed areas (Marcer et al., 2013; Williams et al., 2009).

In SDMs, presence–absence, presence‐only, or abundance data are used to predict species distribution. Presence–absence data provide valuable information about the availability and prevalence of species in the research area and allow for more ecologically realistic predictions to be made (Elith & Leathwick, 2009; Phillips, Dudik, Elith, Graham, & Lehmann, 2009). Either the presence–absence or the abundance of vascular plants is affected by three main groups of factors: direct, indirect, and resource gradients (Austin, 2002; Franklin, 2009; Guisan & Zimmermann, 2000). Additionally, the occurrence of a herbaceous plant species in a forest can be affected by overstory and understory species, canopy closure, and disturbances such as human or animal activities. These biotic factors are difficult to measure and analyze, and they are often ignored in SDMs, even when they are necessary to make realistic predictions (Wisz et al., 2013). M. latifolium usually grows in forest understories. The occurrence of this species might be affected by overstory tree species, development stage of trees, canopy, and characteristics of shrub layer. Moreover, the distances of sample plots to the nearest settlement area may cause indirect human and livestock disturbance effects.

In many species distribution modeling studies, the model is established using selected variables based on the accumulated ecological literature (Porfirio et al., 2014). However, the terrain effects on plant distribution can be explained better by making use of variables derived from digital elevation models (DEMs). These are variables that may have indirect effects on the distribution and abundance of plants. Additionally, annual climate variables are usually used in plant modeling studies. However, climate data should be evaluated on a monthly and seasonal basis. Because herbaceous plants have different life cycles and many different characteristics such as root depth and stem structure, they are more affected than trees by extreme climate values, short term, and seasonal fluctuations (Brovkin, 2002). There have been limited attempts to use fine‐scale DEM‐derived variables and monthly climatic data in species distribution studies.

A variety of methods, such as BIOCLIM (Nix, 1986), MaxEnt (Phillips, Anderson, & Schapire, 2006), DOMAIN (Carpenter, Gillison, & Winter, 1993), GAM (Hastie & Tibshirani, 1990), GLM (McCullough & Nelder, 1989), and random forest (Breiman, 2001), can be used in SDMs. However, this study focuses on identifying species–environment relationships and on estimating the realistic potential distribution area of the species, not on comparing the results of different modeling methods. The aim of this study is to determine the influence of climatic, land surface, geologic, and biotic variables on the distribution of M. latifolium. The study also aims to evaluate the prediction power of models fitted with the “Boosted Regression Trees” (BRT) method based on presence/absence data and a large environmental variable data set. We also summarize the relative importance of predictor variables. The BRT method was preferred in this study because it provides highly accurate predictions of species distribution models and variable shrinkage (Elith, Leathwick, & Hastie, 2008), and it is more sensitive to local site conditions (Falk & Mellert, 2011).

2. Materials and Methods

2.1. Study area

The study was conducted in the Gönen Dam watershed area, which covers 113,700 ha and ranges from 90 to 1,400 m a.s.l. (Figure 2). According to long‐term data from the nearest meteorological station located in the Yenice Province, long‐term average of annual total precipitation is 847.3 mm, and the mean annual temperature is 12.8°C. The Gönen Dam watershed area is located in the northeast Kazdağı Mountains (formerly known as Ida Mountain) in West Anatolia (26.960‐27.540°E, 39.640‐40.100°N). The Kazdağı Mountains consist of several mountain peaks and plateaus and were classified as an IPA (important plant area) not only for Turkey but also for Europe because they contain a high numbers of endemic and rare plant species (Özhatay & Özhatay, 2005). Forests in the Kazdağı Mountains are composed of both pure and mixed conifer and broadleaf trees, such as Pinus nigra J.F. Arnold subsp. pallasiana (Lamb.) Holmboe, Pinus brutia Ten., Abies nordmanniana (Steven) Spach subsp. equi‐trojani (Asch. & Sint. ex Boiss.) Coode & Cullen, Quercus sp., Fagus orientalis Lipsky, maquies, and thickets.

Figure 2

Location of the studied area (filled blue) and distribution of Muscari latifolium incidence on a 3 × 3 km grid in Gönen Dam watershed (Turkey)

2.2. Species data

To obtain a representative sample (Araujo & Guisan, 2006) of M. latifolium occurrence in the study area, it was systematically divided into 3 km × 3 km grids. Then, a 20 m × 20 m quadrat was randomly assigned in each grid, excluding agricultural and residential areas. To avoid edge effects, the quadrats were assigned at least 50 m away from roads. A total of 105 plots in the study area were in managed forests (Figure 2). Therefore, the species incidence consists of 35 presence and 70 absence records. M. latifolium was detected at altitudes ranging from 189 to 885 m in the study area, although the reported range was between 1,100 and 1,800 m (Davis & Stuart, 1984).

This study uses M. latifolium presence–absence data as the response variable. As suggested by Lobo, Jiménez‐Valverde, and Hortal (2010) and Hijmans (2012), we paid attention to the quality of the occurrence data and collected this vegetation data in May, June, and July 2012 by carefully revisiting the study area. The presence–absence of M. latifolium was recorded in five 1 m × 1 m subplots, one in the center and four at the corners of each 20 m × 20 m quadrat. It was considered present in a sample plot even if it was only detected in one of the five subsample plots. All trees with a diameter at breast height (dbh) larger than 7 cm were measured within each sample plot. At the same time, all shrubs were identified, each shrub species was counted, and the coverage percentage of each shrub species was recorded. We collected specimens of species which could not be identified in the field and identified them later in the Forest Faculty of Istanbul University Herbarium (ISTO) using the Flora of Turkey (Davis, 1965–1985; Davis, Mill, & Tan, 1988; Güner, Özhatay, Ekim, & Başer, 2000), and these specimens were deposited in the ISTO.

2.3. Environmental data

We selected environmental predictor variables used in previous SDM studies (Beaumont, Hughes, & Poulsen, 2005; Elith et al., 2006; Lobo et al., 2010; Warren & Seifert, 2011) and added fine‐scale topographic variables and monthly climatic data. Monthly climatic variables and bioclimatic variables were obtained from WorldClim database ( These data are a set of global climate layers with a spatial resolution of approximately 1 km2 (Hijmans, Cameron, Parra, & Albert, 2005).

In addition to climate data, this study used fine‐scale topographic variables obtained from terrain analysis that affect microclimate and other ecological processes. A total of 60 topographic variables such as slope, aspect, and curvature were derived from the ASTER DEM with a 30‐m resolution using the SAGA GIS terrain analysis functions (Conrad et al., 2015).

Solar radiation affects vegetation pattern, plant distribution, and growth by influencing near‐surface air temperature, soil temperature, and soil moisture within a region (Bennie, Huntleya, Wiltshirea, Hill, & Baxtera, 2008; Coblentz & Riitters, 2004). Continuous surface solar radiation data could be obtained from interpolation of weather station data, meteorological satellite data, and modeling solar radiation with GIS, and we preferred to use the latter method to calculate spatial solar radiation considering practical and widespread usage in natural studies. The “potential incoming solar radiation” module of SAGA GIS can be computed solar radiation for an instant time or a given day/week/month/year. Monthly solar radiation (direct solar radiation, diffuse solar radiation, total solar radiation, direct‐to‐diffuse solar radiation ratio, and the duration of solar radiation) was calculated taking the terrain shade effect into account using SAGA GIS (Conrad et al., 2015) under clear‐sky conditions.

There is a strong connection between bedrock composition and vegetation (Hahm, Riebe, Lukens, & Araki, 2014). A bedrock map was obtained from a 1/25.000 scale geological map prepared by the General Directorate of Mineral Research and Exploration (MTA). Bedrock type is the only categorical variable that was used in the study.

According to the literature (Davis & Stuart, 1984) and our observations in the field, M. latifolium requires specific habitat conditions and plant associations to survive and maintain its population. Therefore, some properties of trees and the shrub layer were used to determine the habitat of the plant and to estimate the potential distribution of the plant. The number of tree species per diameter class (8‐ 10.9, 11–19.9, 20–35.9, 36–51.9, 52–79,9 larger than 80 cm) of the 24 tree species existing in the sampling plots was calculated by the cumulative number of trees using the R package “vegclust” (De Cáceres, Font, & Olivia, 2010). The abundance‐cover value, richness, Shannon, Simpson, inverse Simpson, evenness, j evenness, and Berger indices of 73 species in shrub layer were also used.

To handle the effect of humans and livestock, we used proximity to the nearest residential area and the population of the area. The Euclidian distances of sample plots to the nearest residential areas were calculated using the “r.grow.distance” function on GRASS GIS (GRASS Development Team, 2014), and a raster output map was obtained. This variable was taken as it is indicating the impact of indirect human and domestic livestock grazing.

These direct, indirect, and resource variables obtained from GIS data layers used in the study were uploaded to the spatial point vector layer of sample plots using SAGA GIS software. Thus, a data matrix consisting of 416 aforementioned environmental variables (Table 1) and one response variable was prepared for further operations. Preprocessing was performed to achieve better model results before analyses were performed. First, zero‐variance predictors were removed for computing cost even though tree‐based models are impervious to this type of predictors (Kuhn & Johnson, 2013). Because we have more predictors than samples, we handled multicollinearity of DEM‐derived data by the simple five steps way suggested by Kuhn and Johnson (2013) instead of using a variance inflation factor. We did not do multicollinearity analysis for climatic variables because determining the true month of influential climatic variables and BRT is less sensitive than other methods for collinearity (Dormann et al., 2013). After preprocessing, 247 predictor variables remained for use in analysis. Figure 3 shows the study analyses process.

Figure 3

Schematic representation of the analysis steps used in the study

Table 1

Environmental variables used to model Muscari latifolium distribution in the study area (numbers of variable given in the parenthesis)

2.4. Statistical methods

2.4.1. BRTs

To specify the factors affecting the species’ distribution, we used BRTs (aka gradient boosting tree). BRT is a machine learning technique and has important advantages for tree‐based methods. Not only can it fit complex nonlinear relationships, but it can also handle interaction effects between predictors automatically (Elith et al., 2008). Detection of important relationships from large sets of predictor variables can be achieved (Barker, Cumming, & Darveau, 2014). Relatively poor predictive performance drawbacks of single tree models are tackled by BRT (Elith et al., 2008). Wisz et al. (2008) evaluated 12 algorithms for 46 species at three sample sizes (10, 30, and 100 records) and found that gbm was the best performing prediction algorithm at sample sizes 30 and 100.

2.4.2. Model building

We used the dismo (Hijmans, Phillips, Leathwick, & Elith, 2015), gbm (Ridgeway, 2013), and raster (Hijmans & Etten, 2013) packages from the R statistical environment (R Development Core Team, 2014) for fit models, assessing relative contributions, making predictions, and mapping distribution. To prevent overfitting and determining user‐defined parameters used in BRTs, we evaluated tree complexity (1, 3, 5, 7), learning rate (0.1, 0.05, 0.01, 0.005, 0.001, 0.0005), and bag fraction (0.5, 0.75). Based on tenfold cross‐validation results, we selected 7 for tree complexity, 0.5 for bag fraction, and 0.005 for learning rate to achieve more than 1,000 trees suggested by Elith et al. (2008). Using these parameters, we built models with 105 M. latifolium incidences and a 247 environmental variable matrix. To reduce environmental noninformative variables, we simplified this model with the “gbm.simplfy” function (Elith & Leathwick, 2014) and removed 233 environmental variables. Simplification builds many models and drops unimportant variables using methods similar to backward selection in regression (Elith et al., 2008). Thus, 14 environmental variables (Table 2) remain to be used in the further steps.

Table 2

Most important variables selected according to final model performance

2.4.3. Model evaluation

We assessed the predictive performance of models using repeated subsampling processes. Ten random subdata sets were created from the entire data set. Each partition was created randomly selecting 70% (n = 74) presence/absence localities as training data, and the other 30% (n = 31) were selected as testing data. We used the area under a receiver operating characteristic curve (AUC) to evaluate the performance of each model. This metric is calculated from the receiver operating characteristic (ROC) plot that gives the false‐positive error rate (1‐specificity) on the x axis and the true positive rate (sensitivity) on the y axis (Franklin, 2009). The AUC is determined through summing the area under the ROC curve and taking the value between 0.5 and 1.0. Although Harrell (2001) states a threshold of 0.8 AUC value for models is necessary, Franklin (2009) states that a threshold of 0.5–0.7 AUC is considered poor, 0.7–0.9 AUC is considered moderate, and >0.9 AUC is considered high model performance. We created 10 models with tenfold cross‐validated train data sets using 14 environmental and one response variables. Then, predictive performance of these models was calculated on 10 replicate test data sets.

2.4.4. Variable contributions and response curves

While assessing predictive performance for environmental variables, contribution to the model was also calculated over the 10 BRT model replicates. The most influential variables according to the sum of the relative influences of environmental variables in all models were selected and evaluated to determine the ecological requirements of the species.

2.4.5. Spatial prediction

A final spatial prediction map was created from 13 of the 14 most important variables except Sorbus torminalis (L.) Crantz cover value. Potential spatial distribution of the M. latifolium prediction map was produced using a raster layer of these most important variables, and 1 of 10 models has the best prediction power. This map was produced with only part of the study area because not all of the forest survey data were up to date. These field survey data were interpolated with the regularized spline with tension method (Mitasova et al., 1995) which gives good prediction results for forest tree size attributes (Destan, Yılmaz, & Şahin, 2013).

3. Results

3.1. Model performance

The relationship between M. latifolium distribution and environmental variables was analyzed using 10 repeated BRTs models. These models’ accuracy was determined compared to test data sets. The overall average accuracy AUC value is 0.8. In total, 2 of the 10 models (m1, m2) were the most successful with an AUC value 0.93 (Table 3). Three models (Model 3, 4, and 9) gave AUC values that can be considered successful in the 0.80–0.9 range. While four models (m5, m6, m8, and m10) had AUC values between 0.70 and 0.8, only one model (m7) had an AUC value lower than 0.70 (0.68).

Table 3

Performance of 10 repeated boosted regression tree models

3.2. Variable contributions and response curves

According to their relative contributions from 10 repeated BRT models, the seven most influential variables (the relative contribution average is greater than five) account for about the 70% of relative importance. Fourteen variables included in the final model in decreasing order of relative importance are ranked as follows: bedrock type (Bedrock), number of Quercus cerris L. (Qc1), precipitation of wettest month (Bio13), number of P. nigra (diameter >36 cm—Pn4), Northness, sunset September (Sunsetsep), S. torminalis cover value in shrub layer (Sortorm), proximity to residential areas (Growdist), temperature seasonality (Bio4), minimum curvature (Mincur), direct‐to‐diffuse insolation ratio in July (Dir2difJul), duration of insolation in November (Durinsnov), direct‐to‐diffuse insolation ratio in March (Dir2difMar), and direct‐to‐diffuse insolation ratio in November (Dir2difnov) (Table 4).

Table 4

Minimum, maximum, and average relative contributions (%) of the most influential environmental predictors calculated using tenfold cross‐validated BRT models of 10 random subsampled train data sets

Among those fourteen variables, bedrock type was the most influential variable on the distribution of M. latifolium. Six bedrock types are contained in 89 of 105 sample plots (85%) (Table 5). The numbers of sample plots where the species was absent on granodiorite, sandstone, Miocene‐aged andesitic tuff, Oligocene‐aged andesitic tuff, schist, and gneiss–mica‐schist bedrock types were 12, 12, 13, 9, 6, and 5, respectively, while the numbers of sample plots in which the species existed were 3, 2, 1, 8, 5, and 9, respectively (Table 5). Muscari latifolium was present more often than it was absent in plots containing only the gneiss–mica‐schist bedrock type (five absent, nine present).

Table 5

Presence/absence of Muscari latifolium on the six main bedrock types

Occurrences were closely associated with overstory trees. Qc1 was the second most important variable, and Pn4 was the fourth most important variable. The presence of the species in the field is closely associated with Qc1 and Pn4. Qc1 has a negative effect if the number of trees is less than five, and Pn4 also has a negative effect if the number of trees at this diameter class is less than three. We investigated these associations from the data set and found that according to the data set, P. nigra did not occur in six of 35 sample plots where M. latifolium was present while Q. cerris was not detected in 13 of 35 sample plots. Additionally, Q. cerris did not occur in two of six sample plots where M. latifolium was present, but P. nigra was absent. P. nigra did not occur in 23 of the 70 sample plots where M. latifolium was absent, and Q. cerris also did not occur in 45 of these plots. Quercus cerris did not occur in 14 of 23 sample plots in which both M. latifolium and P. nigra were absent.

The third most important variable was Bio13 (December is the wettest month in the study area). A minimum of 135 mm precipitation in December precipitation is associated for M. latifolium (Figure 4). The responses of M. latifolium to northness indicate that the species mostly occurs in the northwest and northeast. The Sortorm cover value is more than 1 in shrub layer which is positively associated with distribution of M. latifolium. Muscari latifolium is also positively affected when the distance to residential areas is between 2,000 and 6,000 m and temperature seasonality (standard deviation *100) (bio4) is >66°C. Mincur is another influential variable that has a positive effect when curvature increases. The occurrence of M. latifolium was also associated with the solar radiation variables. The distribution of M. latifolium is negatively affected when the average monthly duration of insolation in November exceeds 5 hr, the direct‐to‐diffuse insolation ratio of July is >7, the direct‐to‐diffuse insolation ratio of November is >1.5, and the direct‐to‐diffuse insolation ratio of March is >2.5, but influenced positively if the sunset of September is later than 17:00 local time.

Figure 4

Partial dependence plots for the 14 most influential variables

3.3. Spatial prediction map

We also assessed the probability of presence/absence points of M. latifolium from a spatial prediction map (Figure 5). The spatial prediction map covered part of study area containing seven presence and 25 absence records of M. latifolium. The average and maximum probability value of presence was 0.64 and 0.97, respectively, and absence was 0.18 and 0.51, respectively.

Figure 5

Potential spatial distribution map of Muscari latifolium obtained using the most influential variables upper left: green area shows the spatially predicted area within the whole study area

4. Discussion

Plant distributions are limited not only when one environmental factor is less than the minimum required but also more than the maximum tolerance for a particular species (Billings, 1952). In this study, we used SDM to improve our understanding of the relationship between M. latifolium distribution and environmental factors. Species distribution modeling provides us valuable information that is useful in management and effective conservation strategies, particularly for rare and endemic plant species. Additionally, assessing the potential impact of climate change on species distribution (Thuiller, Brotons, Araújo, & Lavorel, 2004) that can be projected by SDM allows the development of strategies for sustainable management. The BRTs modeling approach applied here gives a realistic picture of a potential distribution of M. latifolium in the Gönen Dam watershed that can be used for these aims.

Our results showed that the fine‐scale distribution of M. latifolium is controlled mainly by geological, climatic, topographic, solar radiation, and biotic variables at the study area. Analysis performed on these biotic and abiotic variables showed that there were 14 factors that mostly influenced the species’ distribution (Table 4). These variables create the most favorable growth environment for this species.

Bedrock type is proved to be the most influential variable on the distribution of M. latifolium. This is because bedrock is the main factor affecting soil properties such as climate, relief, altitude, and living organisms (Beieler, 1975; Hartmann & Moosdorf, 2012). Moreover, bedrock has an important role explaining differences in vegetation (Hahm et al., 2014). This is mainly related to the fact that soil is developed from different bedrocks in different textures, which may affect the species’ distribution. Sandy soils, where M. latifolium is present, were formed mostly from granodiorite and sandstone. On the other hand, clay soils, where M. latifolium is absent, were derived from schist and mica schist.

Several climatic variables are also proved important for the distribution of M. latifolium. The increase in temperature seasonality had a positive effect on the habitat suitability of M. latifolium, while the species is unable to tolerate lower temperature seasonality. This is likely related to seasonal thermoperiodicity which is the most important factor controlling growth, development, and flowering in geophytes most of which need warm–cold–warm period to their annual life cycle (Khodorova & Boitel‐Conti, 2013). Tolerances of individual species for extreme seasonality are generally conserved across phylogeny. Therefore, temperature seasonality can be used to accurately predict the range limits of species in SDMs (Wiens, Graham, Moen, Smith, & Reeder, 2006). Precipitation during the wettest month (December in the study area) is thought to be a limiting factor of M. latifolium to survive and maintain its population when it is <135 mm. According to Doussi and Thanos (2002), Muscari seeds need a rainy season in early winter to germinate in the Mediterranean climate. December precipitation may affect the distribution of M. latifolium by influencing its germination.

Solar radiation appears to be an important factor on M. latifolium distribution particularly in March, July, September, and November. The sunset in September later than 17:00 has a positive effect on the distribution of M. latifolium. Although bulbous plants seem to be dormant in autumn and winter, active developmental processes continue in this period using reserves which are in the underground organ and are also affected by temperature conditions (Khodorova & Boitel‐Conti, 2013). Therefore, M. latifolium might require more exposure to sunlight in September, whereas it exists only as bulb and seed below soil in this period. On the other hand, the increment of duration of insolation in November and the increment of direct‐to‐diffuse insolation ratio in March and November have a negative effect on the distribution of the species. Doussi and Thanos (2002) indicated that exposure to daylight caused a decrease in germination rate in some Muscari species even if it led to the emergence of secondary seed dormancy. November solar radiation, March solar radiation, and the wettest month (December) precipitation variables directly affect plant germination; therefore, they are noteworthy factors affecting the distribution of the plant in the area. The direct‐to‐diffuse insolation ratio in July might also be associated with maturation process and spreading of seeds.

Northness is also an influential topographic variable on the distribution of M. latifolium; it mostly occurs in the northwest and northeast aspects in the study area. Northness is an important explanatory variable on a fine‐scale because it refers to the solar radiation contrast between north and south faces and it is a limiting factor on the growth period along north faces related to snow cover duration (Lasseur, Joost, & Randin, 2006). Minimum curvature, another topographic influential variable, has a positive effect when curvature increases. To investigate this relationship, we visually interpreted the M. latifolium occurrence map draped over the minimum curvature map and most of the presence was detected at the steep slope convergence areas, mainly on spurs and ridges that have relatively higher minimum curvature values. The minimum curvature is likely to affect soil properties in such a way that it is favorable for the establishment of M. latifolium (Shary, Larisa, Sharaya, & Mitusov, 2002).

The occurrence of the species in the study area was closely associated with biotic variables characterized by overstory tree species, particularly P. nigra and Q. cerris, the coverage value of S. torminalis in the shrub layer, and proximity to residential areas. Muscari latifolium occurs in, pure P. nigra forests, mixed P. nigra and Quercus sp. forest and oak‐dominant mixed deciduous forest. This might be explained by the influence of forest overstory on the herb layer through modifications of resource availability (light, water, and soil nutrients) (Barbier, Gosselin, & Balandier, 2008). Additionally, López, Larrea‐Alcázar, and Ortuño (2009) found that several herbaceous species are associated exclusively with the shrub undercanopy and he suspected that this is caused by facilitation. The proximity of the sampling plots to settlement areas has positive effects only when the plots are between 2,000 and 6,000 m away. The relationship between the distribution of M. latifolium and settlement areas is very complex and hard to explain. The negative effect observed when the sampling plots are closer than 2,000 m to a settlement may be explained as a result of grazing. Ruminants grazing was observed in some areas during the field surveys. In the same way, Louhaichi, Salkini, and Petersen (2009) determined that the number of geophyte species and the percentage of geophytes in a grazed area were dramatically lower than in ungrazed areas in semiarid Mediterranean Ecosystems. Also Chaideftou, Thanos, Bergmeier, Kallimanis, and Dimopoulos (2009) stated that the seeds of many species (such as Muscari neglectum) that exist in the vegetation of grazed areas could not be found in the seed soil bank of the Mediterranean oak forest. Grazing affects species distribution and composition adversely, and more pressure might contribute decline of the species. Although “r.growdist” function of GRASS GIS software (GRASS Development Team, 2014) gives the proximity to the settlement area, calculating it using a function that takes the terrain into account may give better results. Biotic variables are important to understand the fine‐scale distribution and abundance of species (Meineri, Skarpaas, & Vandvik, 2012) and improve both the fit and the predictive power of distribution models (Pellissier et al., 2010).

Our model establishes the importance of geologic, climatic, topographic, solar radiation, and biotic variables to the occurrence of M. latifolium. Due to a lack of regional information on S. torminalis cover and the lack of a raster map, this variable was removed from the spatial prediction map model of M. latifolium. Biotic variables’ data related to vegetation can be obtained from remote sensing images and can increase the accuracy of models (Swatantran et al., 2012; Wilson, Sexton, Jobe, & Haddad, 2013). However, it is not easy to obtain the cover value of S. torminalis with high accuracy from remote sensing images. Nevertheless, this variable provides an important contribution to the knowledge of the habitat requirements and ecological characteristics of the species. Moreover, the potential distribution map of M. latifolium obtained in this study provides a good basis for the management, conservation, and climate change strategies of this species in the study area, although it did not include S. torminalis cover values.

Ideally, ecologically most relevant variables for a species should be used within SDMs. However, when studied species is endemic and priori information is unavailable, the number of variables that could potentially be used to predict species distribution is almost infinite and has collinearity. Hence, variable selection becomes an important issue that BRTs can bring a solution to. Other important issue that we paid attention in the current study is the true absences that provide potentially relevant information on species ecology (Thuiller et al., 2004).

In conclusion, this study provides significant contribution to the knowledge of the habitat requirements and ecological characteristics of M. latifolium. The distribution of this species is explained by a combination of biotic and abiotic factors. The information obtained in this study can be used to support management, conservation, and, if needed, restoration programs for this species. M. latifolium was listed in the low critical (LC) category by Ekim et al. (2000) and endangered (EN) category by Özhatay and Özhatay (2005). The red list category of this fragmented and limited distributed species needs to be revised, and this potential distribution map may help in this effort. The absence of some biotic factors, such as dispersal limitations, and overstory trees, prevent the model from being more robust. However, bioclimatic variables and solar radiation variables were detected as influential factors that affect the distribution of M. latifolium and can provide valuable information about ecological characteristics of M. latifolium. The BRT model used in the study has reasonable model performance and simplifying mechanism that reduce uninformative variables easily.

Conflict of Interest

None declared.


This work was supported by the Scientific Research Projects Coordination Unit of Istanbul University (Project Number 25242).


Yilmaz H, Yilmaz OY, Akyüz YF. Determining the factors affecting the distribution of Muscari latifolium, an endemic plant of Turkey, and a mapping species distribution model. Ecol Evol. 2017;7:1112–1124. doi:10.1002/ece3.2766.


  • Araujo M. B., & Guisan A. (2006). Five (or so) challenges for species distribution modelling. Journal of Biogeography, 33, 1677–1688.
  • Austin M. P. (2002). Spatial prediction of species distribution: An interface between ecological theory and statistical modelling. Ecological Modelling, 157, 101–118.
  • Barbier S., Gosselin F., & Balandier P. (2008). Influence of tree species on understory vegetation diversity and mechanisms involved—A critical review for temperate and boreal forests. Forest Ecology Management, 254, 1–15.
  • Barker N. K., Cumming S. G., & Darveau M. (2014). Models to predict the distribution and abundance of breeding ducks in Canada. Avian Conservation and Ecology, 9, 7.
  • Beaumont L. J., Hughes L., & Poulsen M. (2005). Predicting species distributions: Use of climatic parameters in BIOCLIM and its impact on predictions of species’ current and future distributions. Ecological Modelling, 186, 251–270.
  • Beieler V. E. (1975). Soil survey of Chelan Area, Washington parts of Chelan and Kittitas Counties, US Natural Resources. Washington, DC: Soil Conservation Service.
  • Bennie J., Huntleya B., Wiltshirea A., Hill M. O., & Baxtera R. (2008). Slope, aspect and climate: Spatially explicit and implicit models of topographic microclimate in chalk grassland. Ecological Modelling, 216, 47–59.
  • Billings D. (1952). The environmental complex in relation to plant growth and distribution. The Quarterly Review of Biology, 27, 251–265. [PubMed]
  • Breiman L. (2001). Random forests. Machine Learning, 45, 5–32.
  • Brovkin V. (2002). Climate‐vegetation interaction. Journal of Physics IV, 12, 52–57.
  • Bryan J. E., & Hort F. I. (2002). Bulbs. Portland, OR: Revised Edition Timber Press.
  • Carpenter G., Gillison A. N., & Winter J. (1993). DOMAIN: A flexible modelling procedure for mapping potential distributions of plants and animals. Biodiversity Conservation, 2, 667–680.
  • Chaideftou E., Thanos C. A., Bergmeier E., Kallimanis A., & Dimopoulos P. (2009). Seed bank composition and above‐ground vegetation in response to grazing in sub‐Mediterranean oak forests (NW Greece). Plant Ecology, 201, 255–265.
  • Coblentz D., & Riitters K. H. (2004). Topographic controls on the regional‐scale biodiversity of the south‐western USA. Journal of Biogeography, 31, 1125–1138.
  • Conrad O., Bechtel B., Bock M., Dietrich H., Fischer E., Gerlitz L., … Böhner J. (2015). System for automated geoscientific analyses (SAGA) v. 2.1.4. Geoscientific Model Development, 8, 1991–2007.
  • Davis P. H. (ed.) (1965. –1985). Flora of Turkey and the East Aegean Islands, Vols. 1–9. Edinburgh: Edinburgh University Press.
  • Davis P. H., Mill R. R., & Tan K. (eds.) (1988). Flora of Turkey and the East Aegean Islands, Vol. 10. Edinburgh: Edinburgh University Press.
  • Davis P. H., & Stuart D. C. (1984). Muscari Miller In Davis P. H., editor. (Ed.), Flora of Turkey and the East Aegean Islands, Vol. 8 (pp. 245–263). Edinburg: Edinburgh University Press.
  • De Cáceres M., Font X., & Olivia F. (2010). The management of vegetation classifications with fuzzy clustering. Journal of Vegetation Science, 21, 1138–1151.
  • Demirci S., Özhatay N., & Koçyiğit M. (2013). Muscari erdalii (Asparagaceae, Scilloideae), a new species from Southern Turkey. Phytotaxa, 154, 38–46.
  • Destan S., Yılmaz O. Y., & Şahin A. (2013). Making Objective Forest Stand Maps of Mixed Managed Forest with Spatial Interpolation and Multi‐Criteria Decision Analysis. Iforest, 6, 268–277.
  • Dormann C. F., Elith J., Bacher S., Buchmann C., Carl G., Carré G., … Lautenbach S. (2013). Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. Ecography, 36, 27–46.
  • Doussi M. A., & Thanos C. A. (2002). Ecophysiology of seed germination in Mediterranean geophytes. 1. Muscari spp. Seed Science Research, 12, 193–201.
  • Ekim T., Koyuncu M., Vural M., Duman H., Aytaç Z., & Adıgüzel N. (2000). Türkiye Bitkileri Kırmızı Kitabı (Red Data Book of Turkish Plants). Ankara: TTKD.
  • Elith J., Graham C. H., Anderson R. P., Dudik M., Ferrier S., Guisan A., Hijmans R. J., Huettmann F., Leathwick J. R., Lehmann A., et al. (2006). Novel methods improve prediction of species’ distributions from occurrence data. Ecography, 29, 129–151.
  • Elith J., & Leathwick J. (2009). Species distribution models: Ecological explanation and prediction across space and time. Annual Review of Ecology and Systematics, 40, 677–697.
  • Elith J., & Leathwick J. (2014). Boosted Regression Trees for ecological modeling. Retrieved from [accessed 15 January 2015].
  • Elith J., Leathwick J. R., & Hastie T. (2008). A working guide to boosted regression trees. Journal of Animal Ecology, 77, 802–813. [PubMed]
  • Falk W., & Mellert K. H. (2011). Species distribution models as a tool for forest management planning under climate change: Risk evaluation of Abies alba in Bavaria. ‐. Journal of Vegetation Science, 22, 621–634.
  • Franklin J. (2009). Mapping species distributions: Spatial inference and prediction. Cambridge: Cambridge University Press.
  • Govaerts R., Zonneveld B. J. M., & Zona S. A. (2015). World checklist of Asparagaceae. Retrieved from [accessed 15 October 2015].
  • GRASS Development Team (2014). Geographic resources analysis support system (GRASS) software, version 6.4.1. Open Source Geospatial Foundation. Retrieved from [accessed 15 January 2015].
  • Guisan A., & Thuiller W. (2005). Predicting species distribution: Offering more than simple habitat models. Ecology Letters, 8, 993–1009.
  • Guisan A., & Zimmermann N. (2000). Predictive habitat distribution models in ecology. Ecological Modelling, 135, 147–186.
  • Güner A. (2012). Muscari Mill In Güner A., editor; , Aslan S., editor; , Ekim T., editor; , Vural M., editor; & Babaç M. T., editor. (Eds.), Türkiye Bitkileri Listesi (Damarlı Bitkiler) (pp. 98–100). İstanbul: NGBB & Flora Araştırmaları Derneği.
  • Güner A., Özhatay N., Ekim T., & Başer K. H. C. (ed.) (2000). Flora of Turkey and the East Aegean Islands, Vol. 11. Edinburgh: Edinburgh University Press.
  • Hahm W. J., Riebe C. S., Lukens C. E., & Araki S. (2014). Bedrock composition regulates mountain ecosystems and landscape evolution. Proceedings of the National Academy of Sciences USA, 111, 3338–3343. [PubMed]
  • Harrell F. E. Jr (2001). Regression modeling strategies: With applications to linear models, logistic regression, and survival analysis. New York, NY: Springer.
  • Hartmann J., & Moosdorf N. (2012). The new global lithological map database GLiM: A representation of rock properties at the Earth surface. Geochemistry, Geophysics, 580 Geosystems, 13, 1–37.
  • Hastie T. J., & Tibshirani R. (1990). Generalized additive models. London: Chapman and Hall.
  • Hijmans R. (2012). Cross‐validation of species distribution models: Removing spatial sorting bias and calibration with a null model. Ecology, 93(3), 679–688. [PubMed]
  • Hijmans R. J., Cameron S. E., Parra J. L., & Albert D. L. (2005). Very High resolution interpolated climate surfaces for global land areas. International Journal of Climatology, 25, 1965–1978.
  • Hijmans R. J., & Etten J. V. (2013). raster: geographic data analysis and modeling. R package version 2.1‐25. Retrieved from [accessed 10 January 2015].
  • Hijmans R. J., Phillips S., Leathwick J., & Elith J. (2015). Package ‘dismo’. Species distribution modeling. R package version 1.0‐12. Retrieved from [accessed 10 January 2015].
  • Hopa E., Tümen G., Sevindik E., & Selvi S. (2013). Kaz dağları’nda Yetişen (Balıkesir) Endemik Muscari Mill. (Liliaceae) Taksonları Üzerinde Karşılaştırmalı Morfolojik ve Ekolojik Araştırmalar. BİBAD, 6, 1–5.
  • Işık K. (2011). Rare and endemic species: Why are they prone to extinction? Turkish Journal of Botany, 35, 411–417.
  • Khodorova N. V., & Boitel‐Conti M. (2013). The role of temperature in the growth and flowering of geophytes. Plants, 2, 699–711. [PubMed]
  • Kuhn M., & Johnson K. (2013). Applied predictive modeling. New York, NY: Springer.
  • Lasseur T., Joost S., & Randin C. F. (2006). Very high resolution digital elevation models: Do they improve models of plant species distribution? Ecological Modelling, 198, 139–153.
  • Lobo J. M., Jiménez‐Valverde A., & Hortal J. (2010). The uncertain nature of absences and their importance in species distribution modelling. Ecography, 33, 103–114.
  • Lomba A., Pellissier L., Randin C., Vicente J., Moreira D., Honrado J., & Guisan A. (2010). Overcoming the rare species modelling paradox: A novel hierarchical framework applied to an Iberian endemic plant. Biological Conservation, 143, 2647–2657.
  • López R., Larrea‐Alcázar D., & Ortuño T. (2009). Positive effects of shrubs on herbaceous species richness across several spatial scales: Evidence from the semiarid Andean subtropics. Journal of Vegetation Science, 20(4), 728–734.
  • Louhaichi M., Salkini A. K., & Petersen S. L. (2009). Effect of small ruminant grazing on the plant community characteristics of semiarid Mediterranean ecosystems. International Journal of Agriculture and Biology, 11, 681–689.
  • Marcer A., Sáez L., Molowny‐Horas R., Pons X., & Pino J. (2013). Using species distribution modelling to disentangle realised versus potential distributions for rare species conservation. Biological Conservation, 166, 221–230.
  • McCullough P., & Nelder J. A. (1989). Generalised linear models. London: Chapman and Hall.
  • Meineri E., Skarpaas O., & Vandvik V. (2012). Modeling alpine plant distributions at the landscape scale: Do biotic interactions matter? Ecological Modeling, 231, 1–10.
  • Mitasova H., Mitas L., Brown W. M., Gerdes D. P., Kosinovsky I., & Baker T. (1995). Modeling spatially and temporally distributed phenomena, new methods and tools for GRASS GIS. International Journal of Geographical Information Systems, 9, 433–446.
  • 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. [PubMed]
  • Nix H. A. (1986). A biogeographic analysis of Australian elapid snakes In Longmore R., editor. (Ed.), Atlas of elapid snakes of Australia, Australian flora and fauna series, number 7 (pp. 4–15). Canberra, ACT: Australian Government Publishing Service.
  • Özhatay N., & Özhatay E. (2005). Kazdağı In Özhatay N., editor; , Bayfield A., editor; & Atay S., editor. (Eds.), Türkiye'nin 122 Önemli Bitki Alanı (pp. 37–40). İstanbul: WWF – TDHKV.
  • Pellissier L., Anne Bråthen K., Pottier J., Randin C. F., Vittoz P., Dubuis A., … Guisan A. (2010). Species distribution models reveal apparent competitive and facilitative effects of a dominant species on the distribution of tundra plants. Ecography, 33, 1004–1014. doi:10.1111/j.1600‐0587.2010.06386.x
  • Peterson A. T. (2006). Uses and requirements of ecological niche models and related distributional models. Bioinformatics, 3, 59–72.
  • 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., Dudik M., Elith J., Graham C., & Lehmann A. (2009). Sample selection bias and presence‐only distribution models: Implications for background and pseudo‐absence data. Ecological Applications, 19, 181–197. [PubMed]
  • Pirhan A. F., Yıldırım H., & Altıoğlu Y. (2014). Muscari serpentinicum sp. Nova (Asparagaceae): A new species from western Anatolia, Turkey. Ot Sistematik Botanik Dergisi, 21, 1–14.
  • Porfirio L. L., Harris R. M. B., Lefroy E. C., Hugh S., Gould S. F., Lee G., … Mackey B. (2014). Improving the use of species distribution models in conservation planning and management under climate change. PLoS ONE, 9, e113749. doi:10.1371/journal.pone.0113749 [PubMed]
  • R Development Core Team . (2014). A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; Retrieved from [accessed 10 November 2014].
  • Ridgeway G. (2013). gbm: generalized boosted regression models. Version 2.0. Retrieved from [accessed 10 January 2015].
  • Shary P. A., Larisa S., Sharaya L. S., & Mitusov A. V. (2002). Fundamental quantitative methods of land surface analysis. ‐. Geoderma, 107(1–2), 1–32.,00136-7
  • Swatantran A., Dubayah R., Goetz S., Hofton M., Betts M. G., Sun M., … Holmes R. (2012). Mapping migratory bird prevalence using remote sensing data fusion. PLoS ONE, 7(1), e28922. [PubMed]
  • 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. doi:10.1111/j.0906‐7590.2004.03673.x
  • Van der Maarel E. (1979). Transformation of cover‐abundance values in phytosociology and its effects on community similarity. Vegetatio, 39, 97–114.
  • Warren D. L., & Seifert S. N. (2011). Ecological niche modeling in Maxent: The importance of model complexity and the performance of model selection criteria. Ecological Applications, 21, 335–342. [PubMed]
  • Wiens J. J., Graham C. H., Moen D. S., Smith S. A., & Reeder T. W. (2006). Evolutionary and ecological causes of the latitudinal diversity gradient in hylid frogs: Treefrog trees unearth the roots of high tropical diversity. American Naturalist, 168, 579–596. [PubMed]
  • Williams J. N., Seo C., Thorne J., Nelson J. K., Erwin J., O'Brien J. M., & Schwartz M. W. (2009). Using species distribution models to predict new occurrences for rare plants. Diversity and Distributions, 15, 565–576.
  • Wilson J. W., Sexton J. O., Jobe R. T., & Haddad R. M. (2013). The relative contribution of terrain, land cover, and vegetation structure indices to species distribution models. Biological Conservation, 166, 170–176.
  • Wisz M. S., Hijmans R. J., Li J., Peterson A. T., Graham C. H., & Guisan A. (2008). Effects of sample size on the performance of species distribution models. Diversity Distribution, 14, 763–773.
  • Wisz M. S., Pottier J., Kissling W. D., Pellissier L., Lenoir J., Damgaard C. F., Carsten F., Dorman C. F., Forchhammer M. C., Grytnes J., et al. (2013). The role of biotic interactions in shaping distributions and realised assemblages of species: Implications for species distribution modelling. Biological Reviews of the Cambridge Philosophical Society, 88, 15–30. [PubMed]
  • WorldClim database . Retrieved from [accessed 09 November 2014].
  • Wraga K., & Placek M. (2009). Review of taxons from genus Muscari cultivated in department of ornamental plants in Szczecin. Herba Pol, 55, 348–353.
  • Yıldırım H. (2015). Muscari atillae (Asparagaceae): A new species from Eastern Anatolia, Turkey. Phytotaxa, 213, 291–295.

Articles from Ecology and Evolution are provided here courtesy of Wiley-Blackwell