Modelling the Impacts of Habitat Changes on the Population Density of Eurasian Skylark (Alauda arvensis) Based on Its Landscape Preferences

: The dramatic decline of the abundance of farmland bird species can be related to the level of land-use intensity or the land-cover heterogeneity of rural landscapes. Our study area in central Europe (Hungary) included 3049 skylark observation points and their 600 m buffer zones. We used a very detailed map (20 × 20 m minimum mapping unit), the Hungarian Ecosystem Basemap, as a land-cover dataset for the calculation of three landscape indices: mean patch size (MPS), mean fractal dimension (MFRACT), and Shannon diversity index (SDI) to describe the landscape structure of the study areas. Generalized linear models were used to analyze the effect of land-cover types and landscape patterns on the abundance of the Eurasian skylark ( Alauda arvensis ). According to our findings, the proportions of arable land, open sand steppes, closed grassland patches, and shape complexity and size characteristics of these land cover patches have a positive effect on skylark abundance, while the SDI was negatively associated with the skylark population. On the basis of the used statistical model, the abundance density (individuals/km*) of skylarks could be estimated with 37.77% absolute percentage error and 2.12 mean absolute error. We predicted the skylark population density inside the Natura 2000 Special Protected Area of Hungary which is 0–6 individu-als/km* and 23746 ± 8968 skylarks. The results can be implemented for the landscape management of rural landscapes, and the method used are adaptable for the density estimation of other farmland bird species in rural landscapes. According to our findings, inside the protected areas should increase the proportion, the average size and shape complexity of arable land, salt steppes and meadows, and closed grassland land cover patches.


Introduction
In the terrestrial ecosystems of the world, the dominant land-cover category is agriculture (38%), including the arable-land use type [1]. In Europe, this value is much higher, at 45% (EBCC, 2015). The agricultural land-cover category contains various land-use types with different levels of human impact. The heterogeneity and spatial structure of these land-use/land-cover (LULC) patches vary greatly across rural areas, which has strong impact on farmland-bird diversity in Europe [2,3]. Many articles have determined that the decreasing trend of farmland birds is strongly connected with the intensity of agricultural management (level of use of fertilizers etc.) [4][5][6][7]. Very few studies have investigated the dramatic decline of the abundance of farmland birds, and its connection with change in landscape structure and land-cover heterogeneity [7][8][9]. There are some regional (country)-scale studies that analyze the connection between land-cover types and farmlandbird population data [10][11][12][13][14][15]. These studies have indicated that the abundance of farmland birds is significantly connected with the intensity of agricultural cultivation, crop  identify skylark land-cover preferences on the basis of the local-scale LULC map;  analyze the impact of landscape patterns of preferred and nonpreferred land-cover classes (habitats), and estimate the impact of all LULC-related variables (proportions, shape, and size characteristics of patches, heterogeneity) on skylark abundance; and  estimate, based on our findings, the skylark population density inside the Natura 2000 Special Protection Area (SPA) of Hungary based on the HEB land cover categories. According to our hypotheses the population density of skylark is predictable based on the preferred LULC categories of skylark and landscape indices (proportion of LULC categories and shape and size related landscape metrics). The methodology is adaptable for analyzing the impact of landscape composition on other farmland-bird populations, and for predicting the population density of the skylark, in protected areas, where field observation-based datasets are not available.

Study Area
Hungary is located in the Carpathian basin (45°43′ to 48°35′N and 16°06′ to 22°53′E) in central Europe, and is part of the Pannonian biogeographical region (Figure 1). The total area is 93,033 km*, and its elevation ranges from 77 to 1014 m a.s.l. The most important land-cover type (61%) is agricultural land [31]. A further 20.7% is natural and seminatural grasslands and forest, and 5.5% is built-up area. In the 1990s, a dramatic landscape change was mainly caused by land privatization. Agricultural lands with low quality and poor agroecological conditions were abandoned [32]. The common agricultural policy of the EU (strong decline of grazing livestock) and land abandonment caused the transformation of arable lands into non-cultivated lands, and the fast and spontaneous reforestation of open grasslands [33].

Skylark-Abundance Data
In Hungary, a countrywide bird-monitoring survey has been conducted every year (like in 2015) by approximately 800 field surveyors who add their field-observation datasets into the Hungarian Common Bird Monitoring Database (MMM) [34][35][36]. The volunteers were not randomly distributed across Hungary. The survey allowed that the observers choose their area of observation. Each observation point received two spring visits, and the abundance of birds was observed (by hearing and visually) within a 100-meter radius of each point. There is a minimum 500 m distance between the observation points. The surveyors left a minimum of two weeks between visits in mid-April and mid-June. The count was accomplished between 5:00 and 10:00, when wind speed was less than 5 m/s and there was no rain. Each observation point contains the average number of observed birds which were counted at the point in the two spring visits [34,35]. In 2015, surveyors counted 6763 skylark individuals across 3049 field observation points (mean value: 2.22, maximum: 34, standard deviation: 4.38.). We used MMM survey points from 2015 in the study area because the HUB land-cover map was also available from that time scale. We analyzed the proportion and spatial configuration of the landscape in the 600 m radius surrounding of the MMM observation points. 600 m buffer zone was chosen, because many author found that landscape composition and land cover types have the highest impact on the abundance of this species within this radius [10,37]. Land use types also have an effect on abundance of skylark population within 600 m buffer radius [38]. We used a very detailed (20 × 20 m mapping units) country scale LULC HEB maps [39] for analyses of the LULC characteristics inside these buffer zones. Unfortunately, the more detailed country scale statistical datasets about the crop structure surroundings of the observation points were not available. Most of the MMM observation points (43%) is situated inside the NATURA 2000 SPA Protected areas, where the grasslands are mowed onetime every year after 15th of June.

Land-Cover Database-Hungarian Ecosystem Basemap
The digital LULC HEB was created by the Hungarian Ministry of Agriculture. The basis year of this database is 2015. This very high resolution LULC dataset was based on other LULC maps of the European Copernicus Program, such as Urban Atlas, Corine Land Cover and High-Resolution Layers, and Sentinel-2 images. The dataset has a 20 × 20 m resolution (minimal mapping unit) and three category levels. Six classes in Level 1, 22 classes in Level 2, and 56 classes in Level 3 (see Table A1). The database also contained three additional LULC categories in Level 4. We used the second level for analysis, and regrouped the LULC classes to reduce the number and the likelihood of autocorrelation between them. Our dataset for statistical analyses contained the following main LULC categories inside the buffer zones (  In our investigations, we aggregated the LULC categories of the HEB database, such as "forest", "wetlands and water surfaces" LULC categories (Table A1). The HEB web map and its documentation are freely available (downloadable) on this website: http://alapterkep.termeszetem.hu/ (accessed on Feb 15, 2021). [39]. 46% of the country is arable land and cereals take the 62% of the arable lands. According to the country scale statistical datasets, the proportion of the crop structure in Hungary is 23% wheat, 26% grain maize, 14% sunflower, 7% barley, 5% rape and 7% fodder crops inside the arable lands (Hungarian Central Statistical Office [40]).

Landscape Metrics
The HEB database was applied to calculate size-and shape-related landscape metric parameters. Patch-level landscape indices were calculated for each LULC patch of the HEB database with the V-LATE 2 extension of Arc GIS 10.3 software [41]. Patch level metrics, created for individual land cover patches, characterize the spatial character and context of patches. These patch metrics serve primarily as the computational basis for developing a landscape metric. During our landscape metrics analyses, we calculated the following patch-level landscape metrics, which represent size and shape characteristics of land-cover patches ( Table 1). The mean patch size (MPS) has been widely applied in landscape ecology, since it is commonly agreed that the occurrence and abundance of different species and species richness strongly correlates with the mean patch size. The shape complexity of individual LULC types was quantified by using landscape metrics (MFRACT). We applied the Shannon Diversity Index (SDI) to determine the landscape heterogenety [25]. We calculated these landscape indices (MPS, MFRACT, SDI) inside the 600 m radius buffer zones. where aij represents the area of the j** patch in the i** class, ni represents the number of patches in the i** class, n represents the number of patches (>0).

MFRACT
Mean fractal dimension index equals 2 times the logarithm of the patch perimeter (m) divided by the logarithm of patch area (m*).

= ∑ 2ln ln
where pij represents the perimeter of the j***patch in class i**, aij represents the area of the j***patch in class i**, ni represents the number of patches in the i** class, n represents the number of patches (1-2).
Landscape Heterogeneity SDI The Shannon diversity index (SDI) provides more information about area composition than simply area richness (i.e., the number of landcover types present).

= − ( * ln( ))
where (m) represents the number of different land-cover types, Pi = the relative abundance of different land-cover types in each BMMU quadrant or LUCAS transect.

Statistical Analyses
To understand the relationship between LULC types and skylark abundance, first we had to identify those LULC categories which are selected (used as habitat) by skylark or are avoided. We applied a preliminary test to identify the group of correlated land-cover and landscape index variables using variance inflation factors (VIFs), and the explanatory variables were not linearly related. VIF values were between 0 and 1.9, which shows that the multicollinearity is low between the variables (LULC types and indices). The arableland category was ignored from statistical analyses (model) because in Hungary and other European countries, the agricultural land is the matrix (dominant LULC type) in the landscape, so the proportion of this category shows strong autocorrelations with other LULC types. We used generalized linear models (GLM) to determine the impact of land cover and landscape structure (composition) on skylark abundance. We applied negative binomial models (link = log) to account the overdispersion of skylark-abundance data (tested by overdispersiontest function of AER package in R). Models with all possible combinations of explanatory variables were generated, and we established Akaike's information criterion to rank them with the "dredge" function from the MuMin package in R [44]. We used model averaging for competitive models (delta AICc < 2) to include uncertainty arising from the high number of candidate models (Table A3) [45]. The significance of the variables was estimated by the LmerTest package [46]. We constructed two groups from the LULC categories of the HEB database based on GLM results, namely, preferred (significant positive relation) and nonpreferred (significant negative relation) land-cover types. We analyzed the relationship between the landscape metrics of the preferred (as habitats) and nonpreferred land-cover types, and the skylark abundance data with negative binominal GLM and model averaging. In the next step in our investigation, we analyzed the shape and size characteristics of those LULC types which showed significant positive relation with skylark abundance. These land-cover types were separated into arable lands and grasslands because we wanted to analyze the effect of arable land and grassland landscape metrics on the skylark population. In this model, the arable land category has been used. The distribution of landscape metric variables was not normal, so logarithmic transformation was used to normalize the data. These variables were in different dimensions, so we created a range function in R that transformed the variable values into a number between 0 and 1: where Range function is a number that describes the given number between 0 and 1, na.rm = T means that NA values were removed, min is the minimal value of the list, and max is the maximal value of the list. On the basis of the output of the statistical model, we could describe the optimal landscape configurations for this species.

Model Validation
We calculated the predicted marginal effects (ggeffects package in R) of the preferred land-cover types and their landscape metrics on the skylark population [47]. To validate our model, we set up a training and a testing group (66.6% and 33.3% proportion, respectively) with random sampling (sample.split function from caTools 1.17 package) on the basis of our dataset in R statistics software. We used the predict function from the car package to calculate the estimated skylark-abundance data. Model accuracy was measured by three indices: Spearman's rank correlation to show the relationship between observed and predicted values, mean absolute error to show the distance of the predicted values from the observed values [48], and mean absolute percentage error to show the percentage of error between observed and predicted values [49].

Prediction of Skylark Population in Natura 2000 SPAs
We could estimate skylark population density using the 600 m buffer areas and the HEB dataset. The centers of the buffer zones were in a regular grid (1200 × 1200 m) inside the Natura 2000 SPA dataset. We used the Natura 2000 SPA areas as the basis of our prediction site, because of the Eurasian skylark is a very common indicator species of agrarian landscapes (Natura 2000 Annex I. list). In Hungary the Natura 2000 SPA areas are typical agrarian landscapes which contain Urban areas (1.5%), Croplands (31.7%), Grasslands and other herbaceous vegetation (21.7%), Forest and woodlands (27.8%) and wetland and water surfaces (17.2%). Mowing of the grasslands inside the Natura 2000 sites is regulated by the law. The mowing machine should cut the grass 10 cm above the soil surface. Mowing should not begin before 1 of July, to protect the ground nesting birds. The number of the animals and the method (it is different based on the grassland type) are also regulated by the law. Prediction was performed based on the model results that analyzed the connection between the preferred area and the landscape metrics. Landscape indices were calculated inside the Natura 2000 SPAs. The estimated skylark population was calculated by the predict function in R software. Figure 3 shows the spatial distribution of the 600 m buffer zones.

Relationship between Land-Cover Proportions and Skylark Abundance
Based on GLM results, we identified two main groups (classes) of the LULC categories of the HEB database. Preferred LULC categories that were considered the habitats of the Eurasian skylark because they showed significant positive relation with skylark abundance were those such as salt steppes and meadows, and closed grasslands in hills and mountains. The closed-grasslands LULC category showed the highest significant relation, thereby having the most important effect on skylark abundance. The arable-land LULC category is also a preferred category according to the literature [11,18,29,50]. The nonpreferred group (class) of LULC categories contains land-cover types with significant negative relations with skylark abundance: built-up land, green urban areas, complex cultivation patterns, forests, and wetlands and water surfaces. The complex-cultivation-pattern LULC category had the strongest negative association with the skylark population, followed by wetland and water surfaces, and green urban areas. The relative importance of the significant variables was 100% in all cases (Table 2).

Relationship between Landscape Structure (Compositon) and Skylark Abundance
The landscape metrics that describe the shape and size characteristics of the preferred and nonpreferred LULC classes showed different directions of significant relation with skylark abundance (Table 3). The metrics that describe the shape complexity and size of the LULC patches of preferred LULC categories of the HEB database showed significant positive relations with skylark abundance. The shape complexity (MFRACT index) of the preferred LULC patches has stronger influence on the skylark abundance than the mean patch size (MPS). The shape complexity and size of the nonpreferred LULC categories had significant negative relation with skylark abundance in this case, MPS had higher association with skylark abundance. (Table 3). Land-cover heterogeneity, described with SDI, had a significant negative effect on skylark abundance, which showed that this species prefers a homogeneous landscape.

Impact of Preferred Land-Cover Categories and Their Landscape Metrics
Total grassland proportion had the highest association with skylark abundance, as shown in Table 2; the average size of arable-land patches (MPS) was more important from an abundance point view of this species than the mean patch size (MPS) of grassland patches. The complexity of grassland patches (MFRACT) had a significant positive association with skylark abundance, while the shape characteristics of arable land had no significant relationship with skylark abundance. The predicted marginal-effect graphs visualize the above-described connections between proportions of LULC categories, size-and shape-related landscape indices, and the estimated population density changes of the skylarks (Figure 4). According to the modeled population density changes, in the case of 100% grassland coverage of a hypothetical landscape, we could find about 4-6 skylark individuals/km*. While the connection between the change in proportions of different land-cover types showed a near exponential curve, landscape metrics showed almost flat linear connections with estimated skylark abundance.
On the basis of our results (Table A2), we could create an equation that describes and estimates the skylark population in a given landscape: where is the skylark number density (individual/km*), MPSarable land is the mean patch size of arable land, MPSgrassland is the mean patch size of grasslands, MFRACTgrassland is the mean fractal dimension of grasslands, Areaarable land is the proportion of arable land, and Areagrassland is the proportion of grasslands.

Model Validation
According the validation of our results, there was a significant Spearman's correlation between the observed and predicted skylark abundance values. Mean absolute error shows the distance between the predicted and observed abundance values of this species, which is +-2.12. Mean absolute percentage error (MAPE) shows the prediction accuracy of the model in percentage; in this case, it was 37.77%. The accuracy of this model based on the MAPE was 62.23% (Table 4). If the model contains just the land cover types, the MEA is 2.95; MAPE is 46.56% and the Spearman correlation coefficient is 0.493.

Prediction of Skylark Population of Natura 2000 Special Protection Areas of Hungary
The spatial distribution of the predicted skylark population in each 600-meter zone of the Natura 2000 SPAs of Hungary was very diverse ( Figure 5). The total investigated Natura 2000 SPA was 13,514 km** which cover the most valuable agroecosystems and rural landscapes of Hungary. Based on model prediction (predict function in R) inside these protected areas, approximately 23,746 skylark individuals were predicted. The density of this species is the highest in the agricultural-landscape-dominated areas of the great Hungarian plain.

Discussion
There are several publications analyzing the relationship between skylark and LULC [4,[51][52][53][54] in local small study areas, but our very detailed LULC dataset (HEB) offers a unique opportunity to obtain regional (country)-scale information about this relationship In our study, we considered both datasets describing proportions of LULC categories and landscape indices that describe the shape and size characteristics of preferred (habitat) and nonpreferred LULC categories. Based on our research findings, population density (individuals/km*) could be estimated because there was a significant statistical relationship between proportions, the shape and size characteristics of different LULC types, and the abundance of this farmland bird. One new finding from our research is that, for the estimation of skylark population density, it is necessary to consider landscape indices together with the proportions of different LULC categories because shape (mean fractal dimension) and size (mean shape size) characteristics of these LULC categories also have significant association with skylark abundance. Based on our finding we have predicted the number skylarks inside the Natura 2000 SPA areas in Hungary.

Impact of Proportions of LULC Categories on Skylark Abundance
We could select two LULC groups (classes) from the land-cover types of a very detailed (20 × 20 m resolution) LULC map. Nonpreferred types had negative significant relation with skylark abundance. These were built-up and green urban areas, which negatively affected the population because of the lack of openness and the high proportion of constructed surfaces. Our findings are confirmed by other international publications [4,6,10]. The complex cultivation pattern land-cover type has negative significant relation with skylark data. Other authors underlines that the skylarks do not prefer heterogeneous agricultural lands because this rural landscape contains many different LULC patches, including also those that are not preferable to the skylark, like vineyards, fruit and berry plantations (because of its height, they obscure the view) [10,11,16,17]. Small parcels of, annual crops, city gardens pastures, fallow lands and/or permanent crops somewhere with scattered houses. Forest and wetland LULC categories are well-known nonpreferred land-cover types of the skylark. The skylark is a typical farmland bird; therefore, it is not a surprise that wetland areas, water bodies, and water courses are not suitable habitat types for this species. The main reason of the negative significant relation of the forest is the lack of openness, which is very important for the skylark [10,11,16,55]. In our research we were not take difference between the type of forests, because according to previous studies every types of forest areas are not habitats of this species.
In the estimation of skylark population density, the preferred land-cover types had higher weights (were more important) than those of the nonpreferred LULC categories. Arable land is a well-known habitat type of this farmland bird species according to the international literature [11,18,19,56,57]. Unfortunately, in Hungary is no available detailed country scale spatial statistical data about the cultivated crop types inside the arable lands (cropland) areas. According to the available most detailed Hungarian LULC dataset, the HEB dataset the 57% of Hungary is covered by agricultural fields and its 81% is arable land (Cropland). Grassland and pasture areas are also preferred LULC categories for skylark, [7,8,10,11,[58][59][60][61]. The HEB dataset allow us to analyze the impact of different types of grassland on skylark abundance. We did not find significant statistical relations with open sand steppes and open rocky grasslands because the number of 600 m circle radius observation points of LULC categories have been low, and these landscape conditions (too-fragmented grassland areas with very short and very sparse vegetation) are not suitable for breeding skylarks [57,62]. There was a significant positive relation between skylark abundance, and the LULC categories of salt steppes and meadows, and closed grasslands. Each LULC category is suitable for nesting breeding skylarks because of the medium vegetation height and optimal proportion inside the 600 m radius circles. Our results are similar with those of others, who described strong relationship between closed grasslands and meadows and skylark abundance, the reason of this relation could be the larger amount of food [63][64][65][66]. According to our findings for the prevention of the farmland bird habitats, the EU agri-environmental policy should pay more attention to the management of salt steppes and meadows, and closed grasslands. To increase the population density of skylark, the mean patch size and the proportion of these land cover types (compare to all) in the landscape should increase. In case of the protected grassland areas, one of the biggest ecological problem is the spontaneous spreading of the bush vegetation, which can reduce the skylark habitats. If we want to stop this process, and keep the openness of the landscapes, we should reduce the size and the shape of the bush and forest patches inside these grassland areas. Therefore, we must eradicate the spontaneously spread bush vegetation (which often full of invasive species) by the proper way of grazing or haymaking, the grasslands can keep its size, shape, and openness characteristics in the protected landscapes. This kind of management of protected areas can preserve not only the vegetation diversity of grasslands but it has also important key factor in the skylark habitat protection.

Impact of Land-Cover Categories and Their Landscape Metrics
The landscape metrics of the preferred LULC classes showed positive significant relation with skylark abundance, meaning that, if arable-land and grassland proportion and shape complexity was higher, then the skylark population would also be higher. The landscape metrics of the nonpreferred LULC classes showed negative significant relation with the skylark population, meaning that, in landscapes with small size and in compact-shape nonpreferred LULC categories, skylark population density (abundance) would be higher.
LULC landscape heterogeneity has a negative effect on the skylark in this scale, where one land cover patch can contain more parcels. If landscape heterogeneity increases, the skylark population declines. This species prefer the homogenous LULC structures, which is in accordance with the results of other authors [10,11,16,17,62].
The grassland proportion had the highest association with the skylark population. This species usually nests and feeds in grasslands. The proportion of arable land has a high association with skylark abundance, but the level of its significance is lower. In the case of the MPS, the opposite phenomenon was observed: the MPS of arable lands (arable land patches of HEB) had a higher effect on skylark abundance than that of grassland. The skylark does not prefer small size arable lands (parcels) and grassland fields in that scale, where one arable land patch can contain more parcels [7,51,60,64]. According to Uuemaa et al. 2009 most bird species react more strongly to the composition land cover than to the configuration of landscapes [25]. Our results also show that the LULC proportions and mean patch sizes have stronger impacts on the abundance of this species, than the shape (fractal dimension index) characteristics of the habitat patches. The mean-absolute-percentage-error value (37.77%) was acceptable since, for a more precise prediction, we would have to use more variables (e.g., species and quantity of insects, used pesticides, parcel management) that are not accessible in country-scale analysis. We can determine that the landscape indices improved the model accuracy, based on the Table 4.

Predicted Population Inside the Natura 2000 SPAs
In Hungary, the latest estimated country-wide Eurasian skylark population is from 1999-2002. There is no spatially detailed population estimate. This study is the first estimate for Natura 2000 protected areas in Hungary. There some early 2000s studies about the skylark densities in Europe (Table 5). The studies listed above do not use the shape and size related landscape indices for estimation of the skylark abundance (density). With the combination of the detailed pointbased bird census data, detailed country-wide LULC dataset and landscape indices we can get a more precise prediction of skylark population. Our results are comparable with these previous estimations and the density values are similar [63,67-69].

Conclusions
Landscape composition (proportions, and shape and size characteristics of LULC categories) has significant association with the skylark population. The salt steppes and meadows, and closed grassland serve as habitat for the Eurasian skylark. This study provides new information about the relationship between landscape metrics of the habitat types (shape and size characteristics of patches) and skylark abundance. Fractal dimension index, which describes the shape complexity of grassland patches has a positive impact on the skylark abundance, while the shape complexity of non-habitat types shows opposite relationships with the skylark density. We analyzed them together and could estimate the association of these landscape composition variables (proportions, shape and size characteristics of LULC classes) with skylark abundance. We could estimate skylark population density inside Natura 2000 SPAs in Hungary.
The outcomes of this study can be used for further land use planning, and the habitat design of Natura 2000 SPAs and other protected areas of the rural landscapes. According to our findings, inside the protected areas should increase the proportion, the average size and shape complexity of those LULC types (arable land, salt steppes and meadows, and closed grassland), which shows positive relations with the abundance data of skylark. It is feasible by stopping the spontaneous reforestation and eradicating the spontaneously spread vegetation (especially invasive bush species). The grazing or mowing, the protected grasslands can preserve the size, shape and openness characteristics of these skylark habitats. This kind of environmental management forms help to conserve the habitat types of skylark. The skylark is an area sensitive species and it is an indicator species of farmlands, so the shown methodology is adaptable for analyzing the impact of landscape composition on other farmland bird populations [56,[70][71][72]. The skylark is considered as indicator for monitoring of agricultural landscapes, because its abundance shows strong relationships with other farmland bird species [73].