Tracking Montane Mediterranean grasslands: Analysis of the effects of [620201]

Tracking Montane Mediterranean grasslands: Analysis of the effects of
snow with other related hydro-meteorological variables and land-usechange on pollen emissions
Jose A. Algarraa,⁎,P a l o m aC a r i ñ a n o sb,c, Javier Herrerod, Manuel Delgado-Capelb,
María M. Ramos-Lorentee, Consuelo Díaz de la Guardiab
aCurator, Botanic Garden Detunda-Cueva de Nerja, C/Minerva, 7 edif. Zeus n°3, 18014 Granada, Spain
bDepartment of Botany, Universidad de Granada, 18071 Granada, Spain
cAndalusian Institute for Earth System Research (IISTA), Edf. CEAMA, University of Granada, Av. del Mediterráneo s/n, 18006 Granada, Spain
dFluvial Dynamics and Hydrology Research Group, Andalusian Institute for Earth System Research (IISTA), University of Córdoba, Rabanales Campus, L e o n a r d od aV i n c iB u i l d i n g ,1 4 0 7 1
Córdoba, Spain
eDepartment of Sociology, Faculty of Health Sciences, Av. de la Ilustración n°60, 18071 Granada, Spain
HIGHLIGHTS
•Temporal evolution of the high moun-
tain Mediterranean grasslands was ex-plored.
•Grass pollen emissions were used as in-
dicators of response to environmentalchanges.
•Snow-packs outside the winter season
is one of the most in fluential parameters
on the pollen index.
•Changes in land use in the preferred
habitats of grasses are also driver of
change.GRAPHICAL ABSTRACT
abstract article info
Article history:
Received 24 May 2018
Received in revised form 22 August 2018
Accepted 23 August 2018Available online 24 August 2018
Editor: Elena PAOLETTIThis paper explores the dynamics of temporal evolut ion of the high mountain Mediterranean grasslands,
(Sierra Nevada, Spain SE). The indicator used is the emission of pollen (Pollen Index, PI) with respect to
two important aspects: the incidence of the snow dynamic together with other hydro-meteorological
parameters, and the changes in land use, which can In fluence the evolution of the grasslands throughout
time. The results reveal that pollen emissions in th e last 25 years have shown a slight downward trend,
with large interannual fluctuations, which are a consequence of diverse environmental factors, both
general and speci fic to the area. One of the most in flu e n t i a lp a r a m e t e r so np o l l e nc o n c e n t r a t i o n si s
snow cover, which reinforces the importance of the p resence of snow-packs as water resource outside
t h ew i n t e rs e a s o ni nt h eH i g hM e d i t e r r a n e a nM o u n t ain environments. The changes in land use experi-
enced in the area are a driver of change, especia lly due to the losses experienced in the last
decades in the preferred habitats for many species o f grasses. It can be conclu ded that the vulnerability
of these ecosystems will be affected by an increase in winter temperatures and/or a decrease in rainfall
(climate change) and an increase in the intensity of anthropogenic activities on land use. In thisKeywords:
Alpine grasslands
Mediterranean mountains
Aerobiological monitoringLand use changeHydro-meteorological modelSnow dynamicsScience of the Total Environment 649 (2019) 889 –901
⁎Corresponding author.
E-mail address: [anonimizat] (J.A. Algarra).
https://doi.org/10.1016/j.scitotenv.2018.08.311
0048-9697/© 2018 Elsevier B.V. All rights reserved.
Contents lists available at ScienceDirect
Science of the Total Environment
journal homepage: www.elsevier.com/locate/scitotenv

context, the PI is shown as a useful indicator of global change given its sensitivity to both anthropic and
hydro-meteorological changes. In addition, it has a w ide range of spatial detection and discrimination
capacity by altitudinal dimensions.
© 2018 Elsevier B.V. All rights reserved.
1. Introduction
The grasses constitute a family of monocotyledonous plants that in-
cludes some 10,000 species of cosmopolitan distribution ( Bouchenak-
Khelladi et al., 2008 ). Among the outstanding species are those forming
vegetation units (grasslands) present in most of the terrestrial
ecozones: tropical and sub-tropical, savannas, temperate, flooded, mon-
tane, tundra, desert and xeric ( Ellenberg and Mueller Dombois, 1967 ).
In the Mediterranean Region, the grasslands, including rangeland, pas-
tures, meadows and fodder crops, occupy up to 48% of the territory
and participate in the provision of ecosystem services such as biodiver-
sity conservation, forage and food production, carbon fixation, climate
regulation, soil and water protection, pollination and nutrient storage
(Hönigová et al., 2012 ). The pastures located in the mountainous areas
of the region are also highly diversi fied and rich in endemic plants,
with endemicity rates that can exceed 30% of the total of flora taxa in
some territories ( Medáil and Quézel, 1999 ).
Sierra Nevada is one of the highest mountain ranges in the Mediter-
ranean Region, located in the southeast region of the Iberian Peninsula,
in an east-west direction. Given its pronounced altitudinal gradient, the
grasslands there present several domains. On the one hand, the
psychroxerophilic grassland resides in the most developed soils, with
a predominance of perennial grasses; on the other hand, starting at
2000 m a.s.l, the hygrophilous vegetation corresponds to humid grass-
land in summer that is covered with snow throughout most of the
year, and which is known as “borreguiles ”(Salazar et al., 2001 ). The im-
portant changes experienced in the area in the last 50 years have
highlighted the serious risk of threat that these ecosystems can suffer,
which is exacerbated by the particular environmental characteristics
of the Mediterranean area ( Beniston, 2003 ;Bravo et al., 2008 ). These
important changes in the area include both pronounced climatic varia-
tions ( Pérez-Palazón et al., 2015 ), with signi ficant reduction of the
snow cover period ( Pérez-Luque et al., 2015 ), and changes in land use
and plant cover ( Jiménez-Olivencia et al., 2015 ).
Much of the natural wealth of this territory is given by its particular to-
pographic, climatic, edaphic and geological conditions, which in turn are
the cause of great vulnerability to any change in them. The analysis and
monitoring of climatic conditions and the evolutionary dynamics of
plant communities become fundamental actions both to know the trend
of change and its effects. In the Sierra Nevada, the tracking and monitoring
of current climate conditions and future trends have been carried out for
several decades ( Zamora et al., 2016 ) in the framework of monitoring
programs such as The Sierra Nevada Global Change Observatory
(OBSNEV) ( B o n e te ta l . ,2 0 1 1 ), based on the conceptual framework and
the thematic areas proposed by the Global Change in Mountain Regions
initiative (GLOCHAMORE), through the UNESCO MaB program ( Schaaf,
2012 ) and the Global Observation Research Initiative in Alpine Environ-
ments (GLORIA) ( Pauli et al., 2007 ). The monitoring of climatic variables
is carried out with a fairly extensive network of meteorological stations,
which since 2004 has been expanding towards higher altitudes
(Herrero et al., 2011 ;Algarra and Herrero, 2014, 2016 ), where there is
almost a complete lack of data. These meteorological records are the
pillars underpinning the advance towards the most precise knowledge
of the hydrology of the region. This knowledge, in conjunction with
distributed physically-ba sed hydrological models ( Herrero et al., 2009 ),
allows for the obtaining of time series of hydrological variables, such as
the amount and duration of snow ( Pimentel et al., 2015 ), the potential
evapotranspiration ( Aguilar et al., 2010 ), or the soil moisture, with a spa-
tial resolution within the tens of meters. This information is important notonly for the analysis of the trends on the hydrometeorological variables
themselves ( Pérez-Palazón et al., 2015 ), but also because these abiotic
variables drive the changes in the plant communities that reside in this
physical environment ( Zamora et al., 2016 ).
To obtain information on the effects that environmental changes
may have on plant species and populations, it is necessary to have indi-
cators that allow predicting the intensity of the expected change and the
response to it. This is essential, especially when dealing with highly vul-
nerable ecosystems with low resilience capacity, such as high mountain
grass communities. The amount of pollen emitted into the atmosphere
by the vegetation during its reproductive process has proven to be a
valid indicator to know the factors that have greater effect on the biolog-
ical state of the vegetation and its response to the change. Thus, pollen
records have been used to know the response and environmental be-
haviour of different groups of plants to arid climate conditions
(Cariñanos et al., 2004, 2014 ), as a tool for assessing the status of endan-
gered species ( Cariñanos et al., 2014 ), and to estimate trends in evolu-
tion of forest-forming species ( Cariñanos et al., 2016 ). There are also
some speci fic studies on grasses, in which the exogenous factors and en-
dogenous processes that intervene in the release of pollen into the at-
mosphere have been analyzed ( García-Mozo et al., 2010 ;García de
León et al., 2015 ;Hernández Plaza et al., 2012 ). In addition to these fac-
tors, the impact of land-cover changes on the presence of grass pollen
has also been explored in some areas of the Mediterranean ( García-
Mozo et al., 2016 ), although given the ubiquity and large number of
existing species in some areas, this relationship is not yet clear.
In this context, it could be considered that the high dependence of the
different populations of grasses present in the Sierra Nevada to the envi-
ronmental conditions, and the intensity of climate change expected in the
area, makes it necessary to intensify the tracking of these communities
using indicators of response to these effects. Due to the Sierra Nevada's
high vulnerability to climate change, and its projections of future change,
it would be hypothesized that hydro-meteorological conditions and land
use change are some of the variables with the greatest effect on mountain
grass populations. Therefore, an analysis of pollen emissions of grasses inthe Sierra Nevada environment could be used to describe temporal and
spatial variation of the different grassland communities present in it, as
well as to highlight the main drivers participating in its dynamics. The
aim of this paper is to explore the relationship that could exist between
the evolution of the montane grasslands of the High Mediterranean
Mountain with two important groups of variables not well studied until
now. To this end, pollen emissions derived from grass communities and
their possible response to hydrometeorological and anthropogenic factors
will be analyzed. Among the first, special emphasis is placed on induced
effects by the snow dynamics. The anthropogenic approach focuses on
the possible in fluence of changes in land use that occurred in the same
area during the last decades.
2. Material and methods
2.1. Description of the area of study
The Sierra Nevada ( Fig. 1 ) is a mountain massif with a surface area of
N2000 km
2and a maximum height of 3479 m a.s.l., which extends
linearly 90 km east-west and an average width of about 35 km, in the
southeast region of the Iberian Peninsula. The Sierra Nevada presents
important natural values, which are recognized in the figures of
“Natural Park ”,“National Park ”and “Biosphere Reserve ”. Moreover, it
is one of the main centers of diversity of the Eastern Mediterranean890 J.A. Algarra et al. / Science of the Total Environment 649 (2019) 889 –901

(Molero, 1994 ), has been a Long-Term Ecological Research (LTER) site
since 2008 and is a biodiversity “hotspot ”(Blanca et al., 1998 ;Myers
et al., 2000 ;Médail and Diadema, 2009 ). Of the N2100 floral taxa
cataloged in the area, almost 40% of them are exclusive to the massif
and 80 are endemics to the Sierra Nevada ( Molero and Pérez-Raya,
1987 ). Recently, it has been included in the first World Green List of
Well Managed Protected Areas, accredited by the International Union
for the Conservation of Nature ( IUCN, 2014 ).
This mountainous massif of the Mediterranean Region has a consid-
erable altitude, only surpassed in Western Europe by the Alps. It is lo-
cated in southern latitudes (37°N) and, consequently, a continuous
snow cover is likely to persist above 2500 to 3479 m a.s.l. during the
winter season, often interrupted by periods of intense melting
(Herrero and Polo, 2016 ). The expected alpine climate is modi fied by
the proximity of the Mediterranean Sea (40 km south) which signi fi-
cantly affects the snow dynamics ( Herrero et al., 2009 ).
Average temperature in the Sierra Nevada ranges from −10 °C to 10
°C above 2000 m a.s.l. during the snow season ( Pérez-Palazón et al.,
2015 ). Annual precipitation fluctuates widely (220 mm/year in the dri-
est areas up to 1000 mm/year) with a high spatial variability throughout
the area due to changes in elevation, longitude, and slope exposure
(north –south). The mean annual precipitation on the western side of
the Sierra is 550 mm at 1000 m a.s.l. and 750 mm at 2000 m a.s.l. On
the opposite side, to the east and the northeast, there is an important
rain shadow effect that diminishes this mean annual precipitation
down to 300 mm at 1000 m a.s.l. and 465 mm at 2000 m a.s.l.
(Herrero and Polo, 2016 ). The mean gradient of precipitation with ele-
vation is about 150 mm km−1. Snowfalls occur mainly from November
to April at altitudes above 2000 m a.s.l. At the Refugio Poqueira weather
station (2500 m a.s.l.), the average precipitation is 889 mm yr−1,5 9 %o f
which occurs as snow. Precipitation oscillates annually between
1426 mm for a wet year and 520 mm for a dry one. The fraction of snow-
fall with respect to the total precipitation also varies between 88% and
46%, with a general tendency to be higher with lower annual precipita-
tion values. The difference in total snowfall varies from 910 mm during
a wet year to 335 mm during a dry year ( Herrero and Polo, 2016 ). An-
other parameter of interest is the wind regime in Sierra Nevada. Thisparameter is of great relevance and it is the most extreme weather phe-
nomenon in the Sierra Nevada both by its dominant direction, condi-
tioned by the orientation of the massif, as by the intensity, which can
reach N100 km/h. This regime affects the conditions in the city of Gra-
nada, which is characterized by the high percentage of calms ( N50%
per year) and prevailing winds of the 3rd quadrant (SW), which rollto the south during the months of July and August ( Viedma Muñoz,
1998 ). It is also remarkable the daily dynamics of valley-mountain
breeze, direction W-E day and vice versa in the night hours.
All these climatic and topographical particularities are re flected in the
variety of existing habitats, distributed in 4 biogeographical sectors and
5 of the 6 bioclimatic belts described for the entire Mediterranean Region
(Rivas-Martínez, 2007 ). In our study, conditioned by the location of the
aerobiological sampler in the city of Granada, we will obtain representa-
tion of the communities located in the 4 upper belts: mesomediterranean,
supramediterranean, oromediterranean and cryoromediterranean, that
is, the altitudinal range between 700 and 3479 m a.s.l. The limits for the
thermotypes come from the latest mapping of vegetation that detail the
work area ( Molero et al., 2001 ;REDIAM, 2009 ).
2.2. Description of the target taxonomic group
The grasses are represented in the Sierra Nevada by some 200 spe-
cies, which constitutes 2% of the species of the Poaceae Family. Although
it is not the most abundant botanical family in the area, it is the most
represented in the different plant communities existing in the Sierra Ne-
vada, especially in the absence of trees ( Romero and Morales, 1996 ).
Most taxa are widely distributed, although 7% of them (14 taxa) are ex-
clusive to the Nevadense massif ( Blanca et al., 2001 ). In terms of distri-
bution, they are present both in natural and anthropized environments,
from sea level to the highest peaks of the Sierra ( Blanca and Algarra,
2011 ). In the most altered urban and peri-urban environments, ruderal
species of colonizing behavior are abundant (i.e., Bromus madritensis ,
Aegilops geniculata ,Hordeum murinum ,Lolium perenne ,Dactylis
glomerta ,Holcus lanatus ,Trisetaria panicea orSorghum halepense ,e t c . ) ;
while in the sierra the species will be distributed according to the altitu-
dinal, edaphic and water availability conditions of the substratum,
Fig. 1. Location of study area.891 J.A. Algarra et al. / Science of the Total Environment 649 (2019) 889 –901

giving rise to very diverse communities: psychroxerophilous perennial
grasses on siliceous or carbonaceous soils, pastures associated to forest
formations, and hygrophilous pastures in the highest areas of the Sierra
(Molero and Pérez-Raya, 1987 ;Molero, 1994 ). It is in this last zone
where the arti ficialization of the habitat and changes in the strict cli-
matic conditions have threatened some of these species (Suppl. 1,
Suppl. 2). Although each domain has its characteristic species, remark-able genera such as Poa(whose Poa annua species is present in all
biocimatic belts), Festuca ,Bromus and Stipa , have 20, 12 and 11 species
respectively. The genus Festuca includes some exclusive endemisms
such as F. clementei orF. pseudoeskia , and threatened species such as
F. frigida .
This diversity of habitats is re flected in the habit, with N60% of the
perennial species, compared to 40% of annual species. There is also a
great amplitude regarding the flowering period, since the species that
grow in the lower areas, closer to the sampler, have a main flowering
period between months III –VI (VII) (with the exception of Sorghum
andTrisetaria ), and those that grow in the different bioclimatic belts os-
cillate between V –VII for those of the Mesomediterranean, VI –VIII in the
Supramediterranean and (VI) VII –VIII (IX) for those of the upper belts
(Suppl. 2).
2.3. Aerobiological data
Poaceae pollen data were obtained from the Aerobiological Sam-
pling Unit of Eastern Andalusia, located in the Faculty of Sciences of
the University of Granada (37°11 ′N, 3°57 ′W, 685 m a.s.l.), from which
a continuous data series from 1991 to 2016 was obtained. The distance
to the Natural Park of Sierra Nevada is just 6.6 km in a straight line. Aero-
biological sampling was carried out according to the standardized
methodology of the Spanish Aerobiology Network, in its Quality and
Management Manual ( Galán et al., 2007 ), which recommends the use
of Hirst-type volumetric suction samplers ( Hirst, 1952 ) and the expres-
sion of results in pollen grains/m
3of air/day. Based on the results of the
series, some aerobiological parameters are established, such as the
curve of average daily values of grass pollen and the Pollen Index (PI,
sum of the daily values for each year of study), the start date of
flowering, (i.e. the date from which 1 pollen grain/m3of air/day was re-
corded for at least 5 consecutive days ( García-Mozo et al., 2009 )) and
end of the pollen season (i.e. the date in which at least
1p o l l e ng r a i n / m3of air/day followed by 5 consecutive days with nil pol-
len presence is recorded).
2.4. Hydro-meteorological data
For this study, the meteorological data were obtained from all the
weather stations located in the area of in fluence of the study area
(Suppl. 3). These stations belong to different networks at the national
level (The State Meteorological Agency (AEMET) and the Autonomous
Organization of National Parks (OAPN)), the regional level (The
Agroclimatic Information Network of Andalusia (RIA-JA)) and the local
level (The Universities of Granada and Cordoba (IISTA)). Speci fically,
we used the data for daily and hourly precipitation (27 stations), daily
temperature (30 stations), solar radiation (16 stations), wind speed
(19 stations) and relative humidity (23 stations). Data correction and
gapfilling of the meteorological data series was assessed using cross
correlation between stations. For the analysis of the general trends of
temperature and precipitation in the long term, we used the Arquilla
station (AEMET), at 1652 m a.s.l., which is the closest station to the
study area with the most complete data series and the highest altitude.
Using these meteorological inputs, the hydro-meteorological model
WiMMed ( Herrero et al., 2014 ) is used to generate distributed maps of
meteorological and hydrological variables for the whole Sierra Nevada
at a spatial resolution of 90 × 90 m from 1999 to 2016. These maps
are subsequently averaged for the different five bioclimatic levels de-
fined in REDIAM (2009) ,a n df o re a c ho ft h e five main basins in theSierra Nevada (Genil, Fardes, Guadalfeo, Adra and Andarax). The values
were obtained independently for each region and different combina-
tions of them were used. This facilitates the analysis of the direct rela-
tionships between the series of daily concentrations of pollen data and
the most representative hydro-meteorological parameters of an entire
study area, with the most homogeneous characteristics possible. The
splitting up used gives a total of 21 regions (Suppl. 3).
From WiMMed, 20 representative hydrometeorological variables
are extracted on a daily scale. These data series can be grouped into 10
meteorological (atmospheric) variables, 4 hydrological (snow-related)
variables, and 6 accumulated (both meteorological and hydrological)
variables. The meteorological variables are 1) mean daily wind speed
(Ws_m), 2) daily maximum temperature (Tmx), 3) daily minimum
temperature (Tmn), 4) mean daily temperature (T_m), 5) number of
hours per day with temperature above freezing, total daily (Tum), 6) di-
rect (Rdr), 7) global solar radiation (Rad), daily 8) rainfall (Pliq),
9) snowfall (P_n) and 10) total precipitation as the sum of both (Pre).
The hydrological variables are 1) fraction of the surface covered by
snow (SCn), 2) snow depth (h_n), 3) snowmelt (Fus), and 4) snow
water equivalent –the amount of snow accumulated on the surface in
mm (EAn). In addition, hydrometeorological variables are also com-
bined and/or temporarily treated to obtain other variables representa-
tive of the processes and cycles associated with high mountain
vegetation. For example, in order to jointly quantify all the contributions
of liquid water on the surface, rain and snowmelt are added together by
a daily simple sum (PliqMasFus). Regarding the time scale of the pro-
cesses, to quantify the fact that the pollen collected in a day depends
not only on the characteristics of that day but also on the previous con-
ditions, cumulative variables of the previous 7 days are also obtained.
The new variables obtained are 1) the total rainfall (PliqAcum7),
2) the total snowfall (P_n_Acum7), and 3) the total precipitation as
the sum of both rainfall and snowfall (Pre_Acum7), accumulated in
the last 7 days, 4) the total daily water reaching surface, rainfall plus
snowmelt (PliqMasFus), 5) the total daily water reaching surface accu-
mulated in the last 7 days (PliqMasFus_Acum7) and 6) the total snow-
melt accumulated in the last 7 days (Fus_Acum7).
2.5. The land-use
The information related to the land use variation was recovered
from the Corine Land Cover Project (CLC), which categorizes the
changes that occurred at different geographical levels, at a 1:100.000
scale. This research attended to the second and the third level of appli-
cation, which comprises 15 and 44 land use types, respectively, and in-
dicates the main variations of the surface on the total coverage in the
thermotypes considered in the zone of study, during the 1990 –2012 pe-
riod. Within the study area, 9 land use types from level 2, as well as 19
land use types from the level 3 were identi fied. Regarding the second
level of application, the following Corine classes were found: urban fab-
ric; mine, dump and construction sites; arable land; permanent crops;
heterogeneous agricultural areas; forests; shrub and/or herbaceous veg-
etation associations; and open spaces with little or no vegetation and in-
land waters ( Tables 1 and 2 ).
2.6. Statistical analysis
The variables PI, Start-date and length of the season and meteorolog-
ical trends were fitted to simple linear regression models; slopes of re-
gression equations and determination coef ficients (R2) were examined.
For the comparative analysis of the hydro-meteorological variables
with the PI, all meteorological data generated above were initially con-
sidered. Taking into account the scope of the sampler and following the
principle of parsimony, the exploratory analysis proceeds to limit as
much as possible the hydrographic basins, by altitude and by total num-
ber of variables to obtain maximum signi ficance with the least number
of variables. To clarify the weight of each variable in the pollen emission892 J.A. Algarra et al. / Science of the Total Environment 649 (2019) 889 –901

of the mountain family Poaceae, a classical Generalized Linear Model
(GLM) was applied versus the possibility of using a Principal Compo-
nent Analysis (PCA, Zuur et al., 2007 ), due to there are relatively low
correlation between hydro-meteorological variables (by a non-
parametric Spearman test). In this particular case, the response variable
is the PI, so a GLM is used, estimated by maximum likelihood, following
a Poisson distribution model with a logarithmic link function, with anon-normal distribution of errors and non-constant variance
(Faraway, 2005 ). When the data have a very high percentage of
zeros, which is the case, because there are scattered days in the series
in which there are no pollen grains, then use a variant of GLM, Zero-
inflated Poisson (ZIP) regression ( Long, 1997 ;Kleiber and Zeileis,
2008 ). For to know the better fit of the models was used Akaike In-
formation Criterion (AIC) ant to compare the result of both models
(GLM and ZIP) is compared by Vuong's non-nested test. Following
the principle of maximum parsimony and maximum deviancy (D
2,
proportion of the variance explained by the model), the analysis is
applied by testing different combinations of variables and homoge-
neous spatial areas (combination of river basins and bioclimatic
belts).
To explore the incidence of changes in land use on vegetation, we
compare the average PI data with the variation of the area obtained
from the spatial analysis, by means of a comparative analysis of para-
metric means (t-Student). All analyses were performed with the R soft-
ware ( R Core Team, 2016 ).
3. Results
3.1. Aerobiological data
Grass pollen is detected in the atmosphere of Granada from mid-
February to mid-October, showing several peaks corresponding to the
flowering of species in different habitats. The first part of the curve cor-
responds to the flowering of the ruderal species near the city and the
lower bioclimatic belt, while in May the blooming of the species that
grows in the upper bioclimatic belt of the Sierra Nevada begins
(Fig. 2 ). Between June and September, coinciding with the maximum
peak of the curve, the flowering of the species forming the pastures
takes place in the cacuminal zones of the Sierra, being able to reachvalues of 40 pollen grains/m
3of air/day.
The trend of the aerobiological parameters, despite the lack of statis-
tically signi ficant changes in 1991 –2016, shows that, in the case of the
PI, the annual value has been close to 500 in some years (1995, 2012),
while in others, the records have exceeded 3000 (2002, 2004 and
2007). Although the trend along the series is slightly downward (s =
−1.7746, R2= 0.002), there is a slight rise towards the end (Suppl. 4).
In the case of the dates of the start and the end of pollen season, opposite
tendencies are obtained: a delay in the beginning (s = 0.3608, R2=
0.0293) and an advance in the end (s = −0.1123, R2= 0.001), so the
pollen seasons are getting shorter in terms of duration (Suppl. 4).3.2. Hydro-meteorological data
Through WiMMed, 20 maps were obtained, the same number of
selected variables ( Fig. 3 ). The time range is limited to the 1999 –2016
interval in order to ensure consistency in the data. This limitation
comes from the great shortage of meteorological stations installed in
the years before 1999 in high mountain areas and the lack of hourlydata in those that did exist, which prevents having a good knowledge
of precipitation in the form of snow. The high correlation between me-
teorological variables is a frequent occurrence and usually causes mis-
use when it is overlooked ( Statheropoulos et al., 1998 ;von Storch and
Navarra, 1999 ). As a result, an exploratory correlation analysis is first
carried out (Spearman test) showing a high correlation only among 9
of the 20 variables (Spearman nonparametric test, D values N0.80,
Fig. 3 ). For the exploratory statistical treatment, the results of the 20
variables selected for the 25 established areas are averaged (by combi-
nation of river basins and thermotypes).
The analysis of the series of meteorological data from the Arquilla sta-
tion, selected to analyze genera l trends, shows that during the 1989 –2016
period (Suppl. 5), there is a downward trend in maximum and minimum
temperatures, also in the case of maximums (s = −0.08392, R
2=
0.08303, p = 0.137, s = −0.05874, R2= 0.07589, p = 0.1559, respec-
tively), although in no case is it signi ficant. A decrease in the annual
amount of rainfall was recorded (s = −3.919, R2= 0.02207, p =
0.4506) although the records have large oscillations between years
(Suppl. 5) and meaning that this decrease is insigni ficant. The pre-
cipitations in the form of snow in the Sierra Nevada present a slightly
ascending tendency for the available series of this parameter (2000
to 2014).
3.3. Relationship between daily pollen and hydro-meteorological data
Thanks to the GLM, a model is generated that shows which variables
have the greatest in fluence on the pollen records. The best fito ft h e
model selected the Genil river basin in the top three thermotypes(Fig. 3 , Suppl. 3) and only 9 parameters of the initial 20 ( Table 3 ,
Fig. 4 ): wind speed (Ws_m), maximum Tă (Tmx), minimum Tă (Tmn),
number of hours per day with temperature above 0 °C (Tum), fraction
of surface covered by snow (SCn), direct radiation (Rdr), global radia-
tion (Rad), total precipitation (Pre) and snow melt (Fus). With these pa-
rameters, a 62.1% explanation of the response variable is obtained.
There is a lower signi ficance in the contribution of water by direct pre-
cipitation (p = 0.29224) compared to the higher signi ficance due to
snow melting (p b0.001).
Although this first result is very satisfactory, the series of daily con-
centration pollen counts shows a very high percentage of zeros
(Suppl. 6), which justi fies the exploration of a variant of the GLM,
Zero-in flated Poisson (ZIP) regression. The result of both models is com-
pared by Vuong's non-nested test and a better model is obtained by ap-
plying the ZIP (AIC. ordinary-glm: 53080.13; AIC. zero-in flated:
48536.79; p b0.001). The result of the ZIP slightly modi fied the result
Table 1
Land uses-types from the Corine Land Cover 2nd level of application and variation of the area (in hectares), as a percentage of the total coverage durin g the period 1990 –2012. LEVEL 2:
level 2 of geographical application for land use categorization; Δ(YY-YY): Variation of the surface, in total hectares, attending level 2 between the revisions of 1990, 2000, 2006 and 2012;
L2 S%: Variation of the area, as a percentage of the total coverage, taking into account level 2 between the revisions of 1990 and 2012.#Anthropogenic types.
LEVEL 2 Δ(90-00) Δ(01-06) Δ(07-12) Δ(90-12)L2 L2 S%
11 Urban fabric#0.00 0.00 6.18 6.18 0.02%
13 Mine, dump and construction sites#0.00 0.67 −5.09 −4.42 −0.01%
21 Arable land#18.47 0.00 −221.30 −202.83 −0.53%
22 Permanent crops#−0.81 0.00 220.18 219.37 0.57%
24 Heterogeneous agricultural areas#72.17 0.00 867.15 939.32 2.43%
31 Forests −25.30 0.00 −1191.47 −1216.77 −3.15%
32 Shrub and/or herbaceous vegetation associations −64.49 −0.67 5797.95 5732.80 14.86%
33 Open spaces with little or no vegetation −0.05 0.00 −5476.52 −5476.57 −14.19%
51 Inland waters 0.00 0.00 2.91 2.91 0.01%893 J.A. Algarra et al. / Science of the Total Environment 649 (2019) 889 –901

obtained in the ordinary GLM ( Table 3 ), increasing the level of signi fi-
cance of direct precipitation (Pre), although to a lesser extent than the
contribution by melting snow (Fus).
3.4. Relationship between PI and changes in land use
The total surface that presents a change in land use is different de-
pending on the level of aggregation observed. Level 2 presents a net
change in 17.19% of the study area while level 3 reaches up to 36.97%
of the territory. In both cases, 2 of the units collect practically all the
changes ( Tables 1 and 2 ). The difference in total area is due to a change
between two units integrated in the same group. The same trend has
been identi fied in terms of Land Cover transformation from a CLC classinto another. One of the areas identi fied with the greatest increase is
“321 Natural grasslands, ”which partially stemmed from the transfor-
mation of “324 Sparsely vegetated areas ”and “323 Sclerophyllous veg-
etation ”areas. In terms of surface loss, the land use types most affected
by the land cover change over time are “324 Transitional woodland-
shrub ”and “324 Sparsely vegetated areas, ”transforming their class
mainly into “323 Sclerophyllous vegetation ”and “321 Natural grass-
lands ”(Fig. 5 ). The in fluence of human activity is observed in the
lower layers of the mountain and especially on the surface transformed
into “321 Natural grassland ”.
Evaluating the exchange rate within each period and its total com-
putation ( Table 4 ), the exchange rates in the first two periods do not ex-
ceed 1%, while in the last period the exchange rate does not exceedTable 2
Land uses-types from the Corine Land Cover 3rd level of application and variation of the area (in hectares), as a percentage of the total coverage durin g the period 1990 –2012. LEVEL 3:
level 3 of geographical application for land use categorization; 1990 to 2012: Year of actualization, surface in hectares (ha); Δ(YY-YY): variation of the surface, in total hectares, attending
level 2 between the revisions of 1990, 2000, 2006 and 2012; L3 S%: variation of the area, as a percentage of the total coverage, taking into account level 2 between the revisions of 1990 and
2012; **: major Surface increase stemmed from Land Cover Change (together they exceed 25% of the total area); *: major surface loss stemmed from Land Co ver Change (together they
exceed 25% of the total area).#Anthropogenic types.
LEVEL 3 Δ(90-00) Δ(01-06) Δ(07-12) Δ(90-12)L3 L3 S%
112 Discontinuous urban fabric#0.00 0.00 6.18 6.18 0.02%
131 Mineral extraction sites#0.00 0.67 −5.09 −4.42 −0.01%
211 Non-irrigated arable land#18.47 0.00 −221.30 −202.83 −0.53%
212 Permanently irrigated land#0.00 0.00 0.00 0.00 0.00%
222 Fruit trees and berry plantations#0.00 0.00 254.98 254.98 0.66%
223 Olive groves # −0.81 0.00 −34.79 −35.60 −0.09%
242 Complex cultivation patterns#0.01 0.00 85.84 85.85 0.22%
243 Land principally occupied by agriculture, with signi ficant areas of natural vegetation#72.16 0.00 −176.85 −104.70 −0.27%
244 Agro-forestry areas#0.00 0.00 958.17 958.17 2.48%
311 Broad-leaved forest −25.39 0.00 −609.01 −634.40 −1.64%
312 Coniferous forest 0.07 0.00 −905.04 −904.96 −2.35% *
313 Mixed forest 0.01 0.00 322.58 322.60 0.84%
321 Natural grasslands 0.00 0.00 4262.28 4262.28 11.04% **
322 Moors and heathland 0.00 0.00 2104.71 2104.71 5.45% **323 Sclerophyllous vegetation −38.04 −0.67 6309.23 6270.53 16.25% **
324 Transitional woodland-shrub −26.45 0.00 −6878.27 −6904.72 −17.89% *
332 Bare rocks 0.00 0.00 −1184.24 −1184.25 −3.07%
333 Sparsely vegetated areas −0.04 0.00 −4292.28 −4292.32 −11.12% *
512 Water bodies 0.00 0.00 2.91 2.91 0.01%
Fig. 2. Daily mean airborne Grass pollen counts for Sierra Nevada (Spain) over the period 1991 –2016.894 J.A. Algarra et al. / Science of the Total Environment 649 (2019) 889 –901

Fig. 3. Correlation graph between the 20 meteorolog ical and hydrological variables analyzed. Ws_m : wind speed in m/s; Tmx: maximum temperature; Tmn : minimum temperature; T_m: average temperature; Tum :n u m b e ro fh o u r sp e rd a yw i t h
temperature above 0 °C; SCn: fraction of snow covered surface (in decimal); Rdr: direct radiation (in MJ); Rad: global radiation (in MJ); Pre: total precipitation (mm); Pre_Acum7 : total precipitation accumulated in the last 7-days; PliqMasFus :
precipitation (rain plus snow melt) in mm; PLiqMasFus_Acum7 : precipitation (rain plus snowmelt) accumulated in the last 7-days; Pliq: precipitation (rain) in mm; PLiqAcum7 : precipitation (rain) accumulated in the last 7-days; P_n:
precipitation (snow) in mm; P_n_Acum7 : precipitation (snow) accumulated in the last 7-days; h_n: snow depth; Fus:s n o w m e l t ; Fus_Acum7 : snowmelt accumulated in the last 7-days; EAn: amount of snow accumulated on surface in mm.895 J.A. Algarra et al. / Science of the Total Environment 649 (2019) 889 –901

37.07%. The exchange rate in this period is slightly higher than the rate
observed throughout the period (1990 –2012). It is noteworthy that the
highest areas (cryoromediterranean) have a 0% change in the first two
periods, while the latter reaches 6.11% of total thermotype despite
being highly stable habitats.
When the PI is compared to the net surface change in the 3 intervals
(1990 –2000, 2001 –2006, 2007 –2012) of the CLC that includes the pre-
ferred grass habitats in the Sierra Nevada, there are no signi ficant differ-
ences between the means of the intervals with the mean of the
complete PI series (t-student; p = 0.1979, p = 0.2229; p = 0.7013,
respectively). This is consistent with the scarce variation observed in
thefirst two intervals; however, it is not consistent with the variation
of the last period (33.15%, Table 4 ).
4. Discussion
4.1. Relationship between pollen and hydro-meteorological data
The analysis of pollen emissions derived from the different grassland
formations present in the Sierra Nevada highlights the diversity of fac-
tors involved in them, since in addition to the high number of existing
taxa, N200, the altitudinal gradient favors the existence of numerous
ecosystems, with particular environmental conditions of each biocli-
matic belt. This makes a general analysis more complex. In this study,
the meteorological parameters indicated as the main exogenous factors
with a direct effect on flowering attributes, temperatures and rainfall
were considered. Both parameters showed an altered trend with
respect to historical values, more pronounced at higher altitudes
(B o n e te ta l . ,2 0 1 5 ), and with negative values in the case of rainfall
(Ruiz Sinoga et al., 2011 ). This could explain the dynamics observed in
the aerobiological parameters of the start and end of flowering, which
also have very oscillating values throughout the series, in direct relation
with the conditions of temperature and water availability. Nevertheless,
the in fluence of these parameters will be different depending on the
environmental and topographic characteristics in which the grasslands
are located. Thus, in grassland communities that extend through the
lower and middle zones of the mountain range, the rainfall regime di-
rectly affects the optimal development of flowering, especially if they
occur in the period of 1 to 2 weeks prior to the onset ( Fernández-
González et al., 1999 ;Cariñanos et al., 2004 ;García-Mozo et al., 2010 ).
Ascending in height, other speci fic parameters are incorporated, as it
has been revealed by the GLM. The variables showing most signi ficant
and positive relationship with PI were the snow covered surface fraction
(SCn), the global radiation (Rad), the direct radiation (Rdr) or the num-
ber of hours per day in which the temperatures are higher than 0 °C
(Tum) stand out, which allows for knowing the optimal conditions for
the phenological development of grassland forming species in the
highest bioclimatic belts.The most consistent relationships were found with the indicators
averaged in the Genil basin. In itself, this fact is a positive indicator of
the method since the sampler ( Fig. 1 ) is located in the lower area of
the Genil basin. And at this point it is exposed to the majority in fluences
of the basin not only when the easterly winds prevail, but also with the
Katabatic nocturnal winds coming from the high mountain levels thatrun down the valley in the same direction ( Herrero and Polo, 2016 ;
Montávez et al., 2000 ).
The statistical analysis revealed that the variable with greater weight
in the PI was the fraction of surface covered with snow (SCn), even
ahead of variables related to temperature or water supply that could
be expected to have greater weight. Firstly, this result supports our
starting hypothesis regarding the importance of snow in the dynamics
of mountain grasses. On the other hand, the correlation between vari-
ables is positive, which could be indicative of the adaptation of this
group of grasses to the presence of snow on the peaks, and the harsh cli-
matic conditions that these species can withstand. SCn could also be-
have as a complex variable that represents and encompasses other
factors and variables that allow its presence and permanence in the
peaks: temperature, water supply, vernalization period and even
radiation.
Temperatures have proven to be one of the parameters with the
most in fluence on the phenological behavior of grasses, as it has been
shown in several studies ( Myszkowska, 2014 ;García-Mozo et al.,
2009 ). While some models have used the weighted sum of tempera-
tures above a certain value as a parameter for predicting pollen concen-
trations ( Cannell and Smith, 1983 ), in other studies the chilling units,
defined as the number of cold hours that accumulate below a certain
threshold, have been used. In the studied area, the effective temperature
threshold to promote flowering would be closer to that of alpine grass-
lands, in which the vernalization process and short-day conditions
(short photoperiod) would place the accumulated temperature thresh-
old at around 9 °C ( King and Heide, 2009 ). This could explain that the
peak of the daily average curve values throughout the series is recorded
between June and July, since as it is characteristic in alpine ecotypes.
Plant species in high mountain areas usually present a shorter vernaliza-
tion period ( Colasanti and Coneva, 2009 ) and a secondary induction
that requires a transition to long days enhanced by moderate-high tem-
peratures and a photoperiod of N10 h, permitting the development of
inflorescences and anthesis ( Heide, 1990, 1994 ;Evans, 1969 ). In this
process, the global radiation (Rad), variable with which it has shown a
positive relationship would also had an effect. This variable is very
high in the Sierra Nevada on the frequent sunny days during the
flowering season ( Aguilar et al., 2010 ;Herrero et al., 2011 ).
The water supply, both directly in the form of rain and that from
snow melting, is another parameter that has resulted in a signi ficant
positive relationship with the PI. Grasses, as a majority group of herba-
ceous plants, have a rapid response to precipitation both in the periods
immediately before flowering ( Recio et al., 2010 ), as during the growing
Table 3
Results of standard GLM and Zero-in flated Poisson, coef ficients and signi ficance levels. Ws_m : wind speed in m/s; Tmx: maximum temperature; Tmn : minimum temperature; Tum :n u m –
ber of hours per day with temperature above 0 °C; SCn: fraction of snow covered surface (in decimal); Rdr: direct radiation (in MJ); Rad: global radiation (in MJ); Pre: total precipitation
(mm); Fus: snowmelt.
Standard generalized linear model Zero-in flated Poisson
Factor Estimate Std. error z value Pr( N|z|) Signif. Estimate Std. error z value Pr( N|z|) Signif.
(Intercept) −8,292,894 0,108,903 −76,149 b2e−16 *** −5,861,971 0,120,724 −48,557 b2e−16 ***
Ws_m 0,005012 0,001657 3026 0,002 ** 0,006049 0,001690 3580 0,0003 ***Tmx −0,015467 0,001221 −12,663 b2e−16 *** −0,016094 0,001238 −13,002 b2e−16 ***
Tmn −0,034269 0,001332 −25,729 b2e−16 *** −0,027621 0,001349 −20,478 3,40E −93 ***
Tum 0,058038 0,001039 55,835 b2e−16 *** 0,051883 0,001094 47,426 b2e−16 ***
SCn 0,278,409 0,017868 15,581 b2e−16 *** 0,161,051 0,018705 8610 b2e
−16 ***
Rdr −0,097326 0,001046 −93,029 b2e−16 *** −0,073502 0,001189 −61,819 b2e−16 ***
Rad 0,172,709 0,001440 119,970 b2e−16 *** 0,131,258 0,001711 76,712 b2e−16 ***
Pre 0,000761 0,000723 1053 0,292 0,002458 0,000743 3310 0,0009 ***Fus 0,034213 0,000593 57,744 b2e−16 *** 0,033494 0,000610 54,871 b2e−16 ***896 J.A. Algarra et al. / Science of the Total Environment 649 (2019) 889 –901

season ( Craine et al., 2010 ), which can even favor the successive
flowering in different periods of the year ( González-Minero et al.,
1998 ). The water coming from the melting of snow becomes a funda-
mental resource in the area since the oromediterranean hygrophilic
grasslands (which grow on flooded soils after thawing) depend on it
(Salazar et al., 2001 ). This con firms the great importance of the presence
of accumulations of snow outside the winter season (the so-called
“neveros ”), which occasionally persist until the following winter. On
the one hand, the analysis yields the highest positive values such as
those related to the presence of snow and contribution of water by
fusion (SCn and Fus). And on the other hand, the direct precipitation
factor (Pre) appears with a much lower value and signi ficance. In the
previous analysis, this variable lacked any signi ficance (standard
GLM). This is consistent with the usual characteristics in the high Med-
iterranean mountain, where the precipitation in summer is scarce or to-
tally absent. So it reinforces the hypothesis about the high dependence
of the vegetation on the contribution of water by melting snow in these
environments ( Giménez-Benavides et al., 2007 ). As it has already been
observed in the other high mountain areas of the world, increasingly
shorter and faster periods of snow melt can generate pronounced
changes in pitlands and alpine grasslands ( Rühland et al., 2006 ;
Burrows, 1977 ;Zimmermann and Kienast, 1999 ). In the Sierra Nevada,
although in the last 25 years no signi ficant differences in the composi-
tion and in the abundance of flowering of the borreguiles-forming spe-
cies have been detected, signi ficant changes have been observed in the
flowering attributes ( Pérez-Luque et al., 2015 ), and in the diversity and
increase of the endemicity rate ( Fernández-Calzado et al., 2012 ).
The wind speed (Ws_m) is another in fluential parameter with
which correlation is obtained. It is a prominent agent in the process of
snow accumulation and in the dynamics of dispersion and accumulation
of pollen emitted at high altitude ( Navares and Aznarte, 2016 ;
Cariñanos et al., 2013 ), permitting the transport of pollen grains madeat medium and even long distance ( Oteros et al., 2015 ;Rojo et al.,
2015 ). Although some authors consider pollen from Poaceae to have a
limited dispersion capability, and they are mainly recorded near the
source of the emission ( Peel et al., 2014 ), other studies consider that
an aerobiological sampler may collect pollen grains of emission sources
located within a radius of 30 km, and even greater depending on thecharacteristics of the pollen and the sampling site ( Katelaris et al.,
2004 ;Rojo et al., 2016 ). In the studied area, the wind regime is strongly
conditioned by the orientation of the massif, which generates a daily
dynamic of valley-mountain breeze, direction W-E, which facilitates
the transport of pollen emissions originated in the mountains to the
city during night hours. This process would also be facilitated by the
smaller size of the pollen grains of grassland-forming species in relation
to the grasses that participate in other mixed vegetation formations (i.e.
forests), which have a larger size as an adaptation to lower dynamics of
wind flow ( Radaeski et al., 2016
). Although the entire Poaceae family
shares the same pollen morphology, analysis by optical microscopy
allows differentiating size ranges. In the samples obtained in the aerobi-
ological sampler located in the city of Granada, Poaceae pollen grains of
size range between 24 and 31 μm were frequently observed, compatible
with those of Festuca and Poaspecies ( Geisler, 1945 ).
4.2. Relationship between pollen and changes in land use
The changes in land use experienced in the area highlight the trans-
formations suffered by the territorial surface, and in direct relation, by
the plant cover. These changes are considered a factor of in fluence on
pollen emissions ( Rojo et al., 2016 ;Maya-Manzano et al., 2017 ). The se-
ries of maps analyzed con firm that the typologies of soil that have had
the greatest loss have been the “324 Transitional woodland-shrub ”
(6904.72 ha) and the “333 Sparsely vegetated areas ”(4292.32 ha). To-
gether, the total number of cases includes some of the preferred habitats
Fig. 4. GLM Poisson graphics results for 9 variables selected (red lines: mean of the data; dot lines: Cook's distance). (For interpretation of the reference s to colour in this figure legend, the
reader is referred to the web version of this article.)897 J.A. Algarra et al. / Science of the Total Environment 649 (2019) 889 –901

for many species of grasses such as high mountain grasslands or agricul-
tural lands with important natural vegetation, in which some of the
most contributive species to the atmospheric pollen spectrum, such as
Dactylis glomerata ,Lolium rigidum ,Trisetaria panicea and Vulpia genicu-
late, are abundant ( León-Ruiz et al., 2011 ). This process could explain
the general downward trend of the PI of grasses along the series, by
adding to the remoteness of the ruderal populations, the impact of an-
thropogenic activities ( Chen et al., 2014 ;García-Mozo et al., 2016 ).
Something similar has been detected in the pollen levels of other herba-
ceous species in which changes in land use have been drivers of a signif-
icant decline in the PI, ( Cariñanos et al., 2014 ;Tormo-Molina et al.,2001 ). In the communities of grasses that live in the Sierra Nevada,
the situation has been similar, since up to 40% of the formations of
scrubland and grassland have been replaced by extensive plantations
of conifers, in particular of Pinus sylvestris (Jiménez-Olivencia et al.,
2015 ), in which the high density of individuals and the low light inten-
sity from under the canopy do not favor an adequate development of
herbaceous.
An increase of nearly 10% of the total area is observed in one of the
CLC types that can contribute the most to pollen records ( “321 Natural
grasslands ”) and, in general, the types that contribute the most to
grass (almost in the third period). However, the recorded pollen data
Fig. 5. Change in land use during the period 1990 –2012. CLC in the framework of grass population, scale 1: 100,000.898 J.A. Algarra et al. / Science of the Total Environment 649 (2019) 889 –901

does not seem to re flect an expected increase consistent with this data.
It is also seen that changes in the upper layers, as expected, are very
light in terms of total area and regional. In the cryoromediterranean
belt, changes are registered only in the last period, with most of them
being a decrease of “333 Sparsely vegetated areas ”and a similar
increase of “321 Natural grasslands ”.
However, the two types with greater surface gain (up to 27.29% of
the total) turn out to be two types that, a priori, would provide a large
amount of grasses with their corresponding pollen emissions
(Table 3 ):“323 Sclerophyllous vegetation ”(6309.23 ha) and “321 Nat-
ural grasslands ”(4262.28 ha). To this apparent contradiction are
added two observations of interest: 1) the fact that a large part of the
surface of the cryoromediterranean belt has varied, which is inconsis-
tent with the real change of vegetation in these heights. And 2) the
value of the area lost in “333 Sparsely vegetated areas ”is practically
the same as that gained in “321 Natural grasslands ”(Table 3 ). These
inconsistencies are explained when a methodological change is discov-
ered in the cartography generated by the CLC ( García-Álvarez, 2018 ).
The last update of 2012 was developed from SIOSE (Land occupation in-
formation system in Spain), a database of land occupation at a national
level with a scale higher than that used in the CLC (equivalent
1:25,000). This increase in the degree of detail, together with the differ-
ent cartographic origin, seems to explain the sudden appearance of new
land uses in 2012 (i.e. fruit trees and berry plantations or Agro-forestry
areas) and certain substitutions between types (i.e. “321 Natural grass-
lands ”and “333 Sparsely vegetated areas ”). Therefore, it can be as-
sumed that there is a slight change in land use that creates a slight
decrease in pollen records, while the changes re flected in the last period
are of little help. Likewise, this situation would have a lesser effect at
higher elevations, where changes in land use have been less intense,and therefore climatic conditions continue to be the main drivers of
changes detected in plant communities.
5. Conclusions
The results obtained in this work highlight the important role that
different factors, both natural and derived from changes in land use,
have in the development and maintenance of the high mountain Med-
iterranean grasslands, in particular, in the case of the Sierra Nevada
(Spain SE).
Of the most in fluential variables, snow coverage stands out, followed
by radiation, temperature and indirect water supply (snow melting).
Conversely, direct water input (precipitation) and wind speed are
revealed as the least relevant variables. In relation to snow, the impor-
tance of the presence of snow accumulations outside the winter season
(snow cover) is pointed out, since they become reservoirs of water
available for the grassland-forming species in the periods of greatest re-
productive activity. This point to the enormous sensitivity of the high
mountain Mediterranean grasslands in the face of climate change. A de-
crease in these resources can have direct effects on the diversity, stabil-
ity and response of the most vulnerable and threatened plant
communities. Changes in land use made in recent years, in which
there have been signi ficant losses in some of the preferred habitats of
high mountain grasses, have also contributed to the adequate develop-
ment of the Alpine Grasslands communities, being of greater intensity
in areas of anthropogenic in fluence. As a follow-up measure of thechanges that take place, accurate cartographic information must be
counted, although the spatial analysis must rely on the collection of spe-
cific data in the field. In this context, the PI is shown as a useful indicator
of global change given its sensitivity to both anthropogenic and hydro-
meteorological changes. In addition, it has a wide range of spatial detec-
tion and discrimination capacity by altitudinal dimensions.
Acknowledgements
We want to thank REDIAM, AEMET and the National Parks Network
for the data provided, both meteorological and cartographic. We also
want to show our gratitude to the Spanish Ministry of Economy and
Competitiveness (MINECO) for their support through the project
FENOMED CGL2014-54731-R. We greatly appreciate all the comments
andfine suggestions of the reviewers and Editor Elena Paoletti which
have signi ficantly improved this paper.
Appendix A. Supplementary data
Supplementary data to this article can be found online at https://doi.
org/10.1016/j.scitotenv.2018.08.311 .
References
Aguilar, C., Herrero, J., Polo, M.J., 2010. Topographic effects on solar radiation distribution
in mountainous watersheds and their in fluence on reference evapotranspiration esti-
mates at watershed scale. Hydrol. Earth Syst. Sci. 14, 2479 –2494.
Algarra, J.A., Herrero, J., 2014. 3.3 monitoring the physical characteristics of the snow
layer. In: Aspizua, R., Barea-Azcón, J.M., Bonet, F.J., Pérez-Luque, A.J., Zamora, R.
(Eds.), Sierra Nevada Global-Change Observatory: Monitoring Methodologies.
Consejería de Medio Ambiente, Junta de Andalucía, Granada, pp. 36 –37.
Algarra, J.A., Herrero, J., 2016. Monitoring the physical characteristics of the snow layer.
In: Zamora, R., Pérez-Luque, A.J., Bonet, F.J., Barea-Azcón, J.M., Aspizua, R. (Eds.),
Global Change Impacts in Sierra Nevada: Challenges for Conservation. Consejería de
Medio Ambiente y Ordenación del Territorio, Junta de Andalucía, Granada, pp. 36 –38.
Beniston, M., 2003. Climatic change in mountain regions: a review of possible impacts.
Climate Variability and Change in High Elevation Regions: Past, Present & Future.Springer, Netherlands.
Blanca, G., Algarra, J.A., 2011. Flora del Espacio natural de Sierra Nevada. In: Serra, L. (Ed.),
Jornadas Estatales de estudio y divulgación de la flora de los Parques Nacionales y
Naturales. CAM, Alcoy, pp. 21 –36.
Blanca, G., Cueto, M., Martínez-Lirola, J.M., Molero-Mesa, J., 1998. Threatened vascular
flora of Sierra Nevada (Southern Spain). Biol. Conserv. 8, 269 –285.
Blanca, G., López Onieva, M.R., Lorite, J., Martínez Lirola, M.J., Molero, J., Quintas, S., Ruiz
Girela, M., Varo, M.A., Vidal, S., 2001. Flora amenazada y endémica de Sierra Nevada.
Universidad de Granada, Consejería de Medio Ambiente, Junta de Andalucía.
Bonet, F.J., Aspizua-Cantón, R., Zamora, R., Sánchez, F.J., Cano-Manuel, F.J., Henares, I.,
2011. Sierra Nevada Observatory for monitoring global change: towards the adaptive
management of natural resources. In: Committee, Austrian MaB (Ed.), Biosphere Re-serves in the Mountains of the World. Excellence in the Clouds?Austrian Academy of
Sciences Press, Vienna, pp. 48 –52.
Bonet, F.J., Pérez-Luque, A.J., Pérez-Pérez, R., 2015. Análisis de tendencias (2000 –2014) en
la cubierta de nieve mediante satélite (sensor MODIS). In: Zamora, R., Pérez-Luque,
A.J., Bonet, F.J., Barea-Azcón, J.M., Aspizua, R. (Eds.), La huella del cambio global en Si-erra Nevada. Retos para la conservación. Consejería de Medio Ambiente y Ordenación
del Territorio, Junta de Andalucía, Granada.
Bouchenak-Khelladi, Y., Salamin, N., Savolainen, V., Forest, F., van der Bank, M., Chase,
M.W., Hodkinson, T.R., 2008. Large multi-gene phylogenetic trees of the grasses
(Poaceae ): progress towards complete tribal and generic level sampling. Mol.
Phylogenet. Evol. 47, 488 –505.
Bravo, D.N., Araújo, M.B., Lasanta, T., Moreno, J.I.L., 2008. Climate change in Mediterranean
mountains during the 21st century. Ambio 37, 280 –285.
Burrows, C.J., 1977. Alpine grasslands and snow in the Arthur's Pass and Lewis Pass re-
gions, South Island, New Zealand. N. Z. J. Bot. 15, 665 –686.
Cannell, M.G.R., Smith, R.I., 1983. Thermal time, chill days and prediction of budburst in
Picea sitchensis .J .A p p l .E c o l .2 0( 3 ) ,9 5 1 –963.Table 4
Changes in land uses compared to pollen records by intervals and by thermotypes, with total and partial net change. ΔSurface: surface net variation between intervals in hectares and
percentage of total study area. PI Mean: PI are calculated by means of each interval (mean ± standard deviation).
1990 –2000 2001 –2006 2007 –2012 Total net change
ΔSurface 181.64 (0.47%) 1.33 (0.00%) 28,729.01 (74.44%) 14,268.19 (36.97%)
Cryoromediterranean 0.00 0.00 79.46 79.46 (6.11%)Oromediterranean 4.39 0.00 11,056.01 11,056.01 (89.44%)
Supramediterranean 177.25 1.33 17,593.53 17,593.53 (70.57%)
PI Mean 1611.00 ± 661.75 2658.83 ± 1079.14 1871.67 ± 1000.71 1920.71 ± 934.03899 J.A. Algarra et al. / Science of the Total Environment 649 (2019) 889 –901

Cariñanos, P., Galán, C., Alcázar, P., Dom ınguez, E., 2004. Airborne pollen records response
to climatic conditions in arid areas of the Iberian Peninsula. Environ. Exp. Bot. 52 (1),
11–22.
Cariñanos, P., de la Guardia, C.D., Algarra, J.A., De Linares, C., Irurita, J.M., 2013. The pollen
counts as bioindicator of meteorological trends and tool for assessing the status of en-
dangered species: the case of Artemisia in Sierra Nevada (Spain). Clim. Chang. 119
(3–4), 799 –813.
Cariñanos, P., Alcázar, P., Galán, C., Domínguez, E., 2014. Environmental behaviour of air-
borne Amaranthaceae pollen in the southern part of the Iberian Peninsula, and its
role in future climate scenarios. Sci. Total Environ. 470, 480 –487.
Cariñanos, P., Casares-Porcel, M., de la Guardia, A.V.D., De la Cruz-Márquez, R., de la
Guardia, C.D., 2016. Charting trends in the evolution of the La Alhambra forest (Gra-
nada, Spain) through analysis of pollen-emission dynamics over time. Clim. Chang.135 (3 –4), 453 –466.
Chen, B., Zhang, X., Tao, J., Wu, J., Wang, J., Shi, P., Zhang, Y.Y.C., 2014. The impact of cli-
mate change and anthropogenic activities on alpine grasslands over the Qinghai-
Tibet Plateau. Agric. For. Meteorol. 189 –190, 11 –18.
Colasanti, J., Coneva, V., 2009. Mechanisms of floral induction in grasses: something
borrowed, something new. Plant Physiol. 149, 56 –62.
Craine, J.M., Gene Towne, E., Nippet, J.B., 2010. Climate controls on grass culm production
over a quarter century in a tallgrass prairie. Ecology 91 (7), 2132 –2140.
Ellenberg, H., Mueller Dombois, D., 1967. Tentative physiognomic-ecological classi fication
of plant formations of the earth. Ber. geobot. Inst. ETH. Stiftg. Riibel, Zurich. 37,
pp. 21 –55.
Evans, L.T., 1969. The nature of flower induction. In: Evans, L.T. (Ed.), The Induction of
Flowering, Vol. III. 306. CRC Press, Boca Raton, FL, pp. 306 –323.
Faraway, J.J., 2005. Extending the linear model with r (texts in statistical science). Gener-
alized Linear, Mixed Effects and Nonparametric Regression Models. Chapman & Hall/
CRC, New York.
Fernández-Calzado, M.R., Molero-Mesa, J., Merzouki, A., Casares-Porcel, M., 2012. Vascular
plant diversity and climate change in the upper zone of Sierra Nevada, Spain. Plant
Biosyst. 146, 1044 –1053.
Fernández-González, D., Valencia-Barrera, R., Vega, A., Díaz de la Guardia, C., Trigo, M.M.,
Cariñanos, P., et al., 1999. Analysis of grass pollen concentrations in the atmosphere of
several Spanish sites. Polen 10, 123 –132.
Galán, C., Cariñanos, P., Alcázar, P., Domínguez, E., 2007. Quality and Management Manual
of the Spanish Aerobiology Network. Servicio de Publicaciones de la Universidad deCórdoba, Spain.
García de León, D., García-Mozo, H., Galán, C., Alcázar, P., Lima, M., González-Andújar, J.L.,
2015. Disentangling the effects of feedback structure and climate on Poaceae annual
airborne pollen fluctuations and the possible consequences of climate change. Sci.
Total Environ. 530 –531, 103 –109.
García-Álvarez, D., 2018. The in fluence of scale in LULC modeling. A comparison between
two different LULC maps (SIOSE and CORINE). In: Camacho Olmedo, M., Paegelow,
M., Mas, J.F., Escobar, F. (Eds.), Geomatic Approaches for Modeling Land Change Sce-narios. Lecture Notes in Geoinformation and Cartography. Springer, Cambridge,
pp. 187 –213.
García-Mozo, H., Galán, C., Belmonte, J., Bermejo, D., Candau, P., Díaz de la Guardia, C.,
Elvira, B., Gutiérrez, M., Jato, V., Silva, I., Trigo, M.M., Valencia, R., Chuine, I., 2009.Predicting the start and peak dates of the Poaceae pollen season in Spain using
process-based models. Agric. For. Meteorol. 149, 256 –262.
García-Mozo, H., Galán, C., Alcázar, P., Díaz de la Guardia, C., Nieto-Lugilde, D., Recio, M.,
Hidalgo, P., González-Minero, F., Ruiz, L., Dominguez-Vilches, E., 2010. Trends in
grass pollen season in southern Spain. Aerobiologia 26, 157 –169.
García-Mozo, H., Oteros, J.A., Galán, C., 2016. Impact of land cover changes and climate on
the main airborne pollen types in Southern Spain. Sci. Total Environ. 548 –549,
221 –228.
Geisler, F., 1945. A pollen study of thirty two species of grasses. Butler University Botanical
Studies. vol. 7 (Article 6).
Giménez-Benavides, L., Escudero, A., Iriondo, J.M., 2007. Reproductive limits of a late
flowering high-mountain Mediterranean plant along an elevated climate gradient.
New Phytol. 173, 367 –382.
González-Minero, F.J., Iglesias, I., Jato, M.V., Aira, M.J., Candau, P., Morales, J., Tomas, C.,
1998. Study of the pollen emissions of Urticaceae, Plantaginaceae and Poaceae at
five sites in western Spain. Aerobiologia 14, 117 –125.
Heide, O.M., 1990. Dual floral induction requirements in Phleum alpinum . Ann. Bot. 66,
687 –694.
Heide, O.M., 1994. Control of flowering and reproduction in temperate grasses. New
Phytol. 128, 347 –362.
Hernández Plaza, E., Navarrete, L., Lacasta, C., González-Andújar, J.L., 2012. Fluctuations in
plant populations: role of exogenous and endogenous factors. J. Veg. Sci. 23, 640 –646.
Herrero, J., Polo, M.J., 2016. Evaposublimation from the snow in the Mediterranean moun-
tains of Sierra Nevada (Spain). Cryosphere 10, 2981 –2998. https://doi.org/10.5194/
tc-10-2981-2016 .
Herrero, J., Polo, M.J., Moñino, A., Losada, M.A., 2009. An energy balance snowmelt model
in a Mediterranean site. J. Hydrol. 371, 98 –107. https://doi.org/10.1016/j.
jhydrol.2009.03.021 .
Herrero, J., Aguilar, C., Millares, A., Moñino, A., Polo, M.J., Losada, M.A., 2011. Mediterra-
nean high mountain meteorology from continuous data obtained by a permanent
meteorological station at Sierra Nevada, Spain. European Geosciences Union General
Assembly 2011. Vienna, Austria , pp. 3 –8 April 2011. http://meetingorganizer.coper-
nicus.org/EGU2011/EGU2011-12893.pdf .
Herrero, J., Millares, A., Aguilar, C., Egüen, M., Losada, M.A., Polo, M.J., 2014. Coupling spa-
tial and time scales in the hydrological modelling of Mediterranean regions:
WiMMed. CUNY Academic Works http://academicworks.cuny.edu/cc_conf_hic/315 .
Hirst, J., 1952. An automatic volumetric spore trap. Ann. Appl. Biol. 39, 257 –265.Hönigová, I., Va čkář, D., Lorencová, E., Melichar, J., Götzl, M., Sonderegger, G., Ou šková, V.,
Ho
šek, M., Chobot, K., 2012. Survey on grassland ecosystem services. Report to the
EEA –European Topic Centre on Biological Diversity. Nature Conservation Agency
of the Czech Republic, Prague, p. 78.
IUCN, 2014. The green list for protected areas global standard. IUCN GLPA Steering Group
2014-05-15.
Jiménez-Olivencia, Y., Porcel-Rodríguez, L., Caballero-Calvo, A., Bonet, F.J., 2015. Evolución
de los usos del suelo en Sierra Nevada en los últimos 50 años y cambio de paisaje. In:Zamora, R., Pérez-Luque, A.J., Bonet, F.J., Barea-Azcón, J.M., Aspizua, R. (Eds.), La huella
del cambio global en Sierra Nevada. Retos para la conservación. Consejería de Medio
Ambiente y Ordenación del Territorio, Junta de Andalucía.
Katelaris, C., Burke, T., Byth, K., 2004. Spatial variability in the pollen count in Sydney,
Australia: can one sampling site accurately re flect the pollen count for a region?
Ann Allergy Asthma Immunol 93, 131 –136.
King, R.W., Heide, O.M., 2009. Seasonal flowering and evolution: the heritage from
Charles Darwin. Funct. Plant Biol. 36, 1027 –1036.
Kleiber, C., Zeileis, A., 2008. Applied Econometrics With R. Springer-Verlag, New York.
León-Ruiz, E., Alcázar, P., Domínguez-Vilches, E., Galán, C., 2011. Study of Poaceae phenol-
ogy in a Mediterranean climate, which species contribute most to airborne pollen
counts? Aerobiologia 27, 37 –50.
Long, J.S., 1997. Regression models for categorical and limited dependent variables. Ad-
vanced Quantitative Techniques in the Social Sciences. 7. Sage Publications, Thou-
sand Oaks, CA.
Maya-Manzano, J.M., Sady ś, M., Tormo-Molina, R., Fernández-Rodríguez, S., Oteros, J.,
Silva-Palacios, I., Gonzalo-Garijo, A., 2017. Relationships between airborne pollen
grains, wind direction and land cover using GIS and circular statistics. Sci. Total Envi-
ron. 584 –585, 603 –613.
Médail, F., Diadema, K., 2009. Glacial refugia in fluence plant diversity patterns in the Med-
iterranean Basin. J. Biogeogr. 36 (7), 1333 –1345.
Medáil, F., Quézel, P., 1999. Biodiversity hotspots in the Mediterranean Basin: setting
global conservation priorities. Conserv. Biol. 13, 1510 –1513.
Molero, J., 1994. Baetic and sub-Baetic mountains, Spain. In: Davis, S.D., Heywood, V.H.,
Hamilton, A.C. (Eds.), Centres of Plant Diversity. A Guide and Strategy for Their Con-
servation, Europe, Africa, South West Asia and the Middle East. vol. 1. WWF-UICN
Publication Unit, Cambridge, pp. 55 –58.
Molero, J., Pérez-Raya, F., 1987. La Flora de Sierra Nevada. Avance sobre el catálogo
florístico nevadense. Universidad de Granada.
Molero, J., Pérez-Raya, F., López-Nieto, J.M., El Aallali, A., 2001. Cartograf¡a y evaluación de
la vegetación del Parque Natural de Sierra Nevada. Segunda fase. Universidad de Gra-
nada y Consejería de Medio Ambiente, Junta Andalucía, Granada, p. 239.
Montávez, J.P., Rodríguez, A.J., Jiménez, J.I., 2000. A study of the urban heat island of Gra-
nada. I. J. Climatol. 20, 899 –911.
Myers, N., Mittermeier, R.A., Mittermeier, G.C., da Fonseca, G.A.B., Kent, J., 2000. Biodiver-
sity hotspots for conservation priorities. Nature 403, 853 –858.
Myszkowska, D., 2014. Poaceae pollen in the air depending on the thermal conditions. Int.
J. Biometeorol. 58, 975 –986.
Navares, R., Aznarte, J.L., 2016. Forecasting the start and end of pollen season in Madrid.
ITISE 2016: Advances in Time Series Analysis and Forecasting, pp. 387 –399.
Oteros, J., García-Mozo, H., Alcázar, P., Belmonte, J., Bermejo, D., Boi, M., Cariñanos, P., Díaz
de la Guardia, C., Fernández-González, D., González-Minero, F., Gutiérrez-Bustillo,
A.M., Moreno-Grau, S., Pérez-Badía, R., Rodríguez-Rajo, F.J., Ruíz-Valenzuela, L.,
Suárez-Pérez, J., Trigo, M.M., Domínguez-Vilches, E., Galán, C., 2015. A new method
for determining the sources of airborne particles. J. Environ. Manag. 155, 212 –218.
Pauli, H., Gottfried, M., Reiter, K., Klettner, C., Grabherr, G., 2007. Signals of range expan-
sions and contractions of vascular plants in the high Alps: observations
(1994 –2004) at the GLORIA master site Schrankogel, Tyrol, Austria. Glob. Chang.
Biol. 13, 147 –156.
Peel, R.G., Orby, P.V., Skjoth, C.A., Kennedy, R., Schlünssen, V., Smith, M., et al., 2014. Sea-
sonal variation in diurnal atmospheric grass pollen concentration pro files. Biogeosci-
ences 11, 821 –832.
Pérez-Luque, A.J., Pérez-Pérez, R., Bonet, F.J., 2015. Evolución del clima en los últimos 50
años en Sierra Nevada. In: Zamora, R., Pérez-Luque, A.J., Bonet, F.J., Barea-Azcón,
J.M., Aspizua, R. (Eds.), La huella del cambio global en Sierra Nevada. Retos para la
conservación. Consejería de Medio Ambiente y Ordenación del Territorio, Junta de
Andalucía.
Pérez-Palazón, M.J., Pimentel, R., Herrero, J., Aguilar, C., Perales, J.M., Polo, M.J., 2015. Ex-
treme values of snow-related variables in Mediterranean regions: trends and long-
term forecasting in Sierra Nevada (Spain). Proc. IAHS. 369, pp. 157 –162. https://
doi.org/10.5194/piahs-369-157-2015 .
Pimentel, R., Herrero, J., Zeng, Y., Su, Z., Polo, M.J., 2015. Study of snow dynamics at
subgrid scale in semiarid environments combining terrestrial photography and dataassimilation techniques. J. Hydrometeorol. 16, 563 –578.
R Core Team, 2016. R: A Language and Environment for Statistical Computing. R Founda-
tion for Statistical Computing, Vienna, Austria https://www.R-project.org/ .
Radaeski, J.N., Bauermann, S.G., Pereira, A.B., 2016. Poaceae pollen from Southern Brazil:
distinguishing grasslands (campos) from forests by analyzing a diverse range of
Poaceae species. Front. Plant Sci. 7, 1833.
Recio, M., Docampo, S., García-Sánchez, J., Trigo, M.M., Melgar, M., Cabezudo, B., 2010. In-
fluence of temperature, rainfall and wind trends on grass pollination in Malaga
(western Mediterranean coast). Agric. For. Meteorol. 150 (7 –8), 931 –940.
REDIAM, 2009. Mapa de pisos bioclimáticos de la masa forestal de Andalucía. Consejería
de Medio Ambiente y Ordenación del Territorio, Junta de Andalucia, Sevilla http://
www.juntadeandalucia.es/medioambiente/site/rediam .
Rivas-Martínez, S., 2007. Mapa de series, geoseries y geopermaseries de vegetación de
España. Memoria del mapa de vegetación potencial de España. Parte 1. Itinera
Geobot. 17, 1 –222.900 J.A. Algarra et al. / Science of the Total Environment 649 (2019) 889 –901

Rojo, J., Rapp, A., Lara, B., Fernández-González, F., Pérez-Badia, R., 2015. Effect of land use
and wind direction on the contribution of local sources to airborne pollen. Sci. Total
Environ. 505, 860 –869.
Rojo, J., Orlandi, F., Pérez-Badia, R., Aguilera, F., Ben Dhiab, A., Bouziane, H., et al., 2016.
Modeling olive pollen intensity in the Mediterranean región through analysis of
emission sources. Sci. Total Environ. 551 –552, 73 –82.
Romero, A.T., Morales, C., 1996. El papel de las gramíneas en la flora de Sierra Nevada. In:
Chacón, J., Rosúa, J.L. (Eds.), Primera Conferencia Internacional Sierra Nevada.Conservación y Desarrollo Sostenible. volumen 2, pp. 335 –352 Granada.
Rühland, K., Phadtare, N.R., Pant, R.K., Sangode, S.J., Smol, J.P., 2006. Accelerated melting of
Himalayan snow and ice triggers pronounced changes in a valley peatland from
northern India. Geophys. Res. Lett. 33, L15709.
Ruiz Sinoga, J.D., García Martí, R., Martínez Murillo, J.F., Gabarrón Galeote, M.A., 2011. Pre-
cipitation dynamics in Southern Spain: trends and cycles. Int. J. Climatol. 31,
2281 –2289.
Salazar, C., Lorite, J., García-Fuentes, A., Torres, J.A., Cano, E., Valle, F., 2001. Ap h y t o s o c i o –
logical study of the hygrophilous vegetation of Sierra Nevada (Southern Spain).
Studia Geobotanica 20, 17 –32.
Schaaf, T., 2012. Haciendo frente al cambio global y climático en las montañas. La utilidad
de las reservas de la Biosfera como enclaves para la observación y el seguimiento.Jornadas de Investigación del Parque Nacional y Parque Natural de Sierra Nevada,
Granada.
Statheropoulos, M., Vassiliadis, N., Pappa, A., 1998. Principal component and canonical
correlation analysis for examining air pollution and meteorological data. Atmos. En-
viron. 32, 1087 –1095.
Tormo-Molina, R., Silva-Palacios, I., Muñoz-Rodríguez, A.F., Muñoz, J.T., Corchero, A.M.,
2001. Environmental factors affecting airborne pollen concentration in anemophilous
species of Plantago. Ann. Bot. 87, 1 –8.
Viedma Muñoz, M., 1998. Análisis de las Direcciones de los vientos en Andalucía. 1. Nim-
bus, pp. 153 –168.
Von Storch, H., Navarra, A., 1999. Analysis of Climate Variability. Applications of Statistical
Techniques. 2 edition. Springer, Berlin.
Zamora, R., Pérez-Luque, A.J., Bonet, F.J., Barea-Azcón, J.M., Aspizua, R. (Eds.), 2016. Global
Change Impacts in Sierra Nevada: Challenges for Conservation. Consejería de Medio
Ambiente y Ordenación del Territorio, Junta de Andalucía, Granada.
Zimmermann, N.E., Kienast, F., 1999. Predictive mapping of alpine grasslands in
Switzerland: species versus community approach. J. Veg. Sci. 10, 469 –482.
Zuur, A.F., Ieno, E.N., Smith, G.M., 2007. Analysing Ecological Data. Springer, New York.901 J.A. Algarra et al. / Science of the Total Environment 649 (2019) 889 –901