2.1. STUDY AREA

The study area is located in the zone culturally delimited as Maya region, roughly between 13.30 to 21.40N and from 87.10 to 93.80W (Figure 1), spanning from the lower Lempa River in El Salvador to the Isthmus of Tehuantepec in Mexico, and covering about 324000 km² of tropical lands. The region could be roughly divided into three main geographical zones: (i) the southern highlands that include the mountains of Honduras, Guatemala, and Chiapas; (ii) the southern lowlands that cover Petén in northern Guatemala, northern Honduras, Belize, and the southern portions of Yucatan Peninsula; and (iii) the northern lowlands on the central and northern portion of Yucatan Peninsula.

Figure 1. Study area. A. Elevation map indicating the main geographical zones in the Maya region, and B. Human Influence Index (HII) along the region with their respective frequencies. Yellow dots show the location of the sampled lakes (see Appendix 1 for a list of samples and their geographic coordinates). 

The Climate is warm in the lowlands and becomes cooler towards the highlands, with differences among regions mostly produced by rainfall patterns. High topographic diversity expressed through large elevations gradients and rain shadows lead to significant precipitation variability across the region. In the highlands and the mountains annual precipitation can reach up to 5000 mm, while in northwestern Yucatan it barely reaches 500 mm (^{Méndez and Magaña, 2010}). The climate is primarily controlled by the atmospheric patterns influenced by the Intertropical Convergence Zone and the Bermuda High, which together cause a markedly wet season from May to November, and a dry season from December to May (^{Wilson, 1980}; ^{Magaña et al., 2003}).

Temperature and precipitation are the main drivers of natural vegetation patterns across the Maya lands, whereas precipitation seasonality and edaphic features exert a more local control on vegetation associations. Scrubs and tropical deciduous forests naturally cover the northernmost and driest portion of the Yucatan Peninsula, whereas tropical semi-deciduous forests dominate the central and northeast areas, where conditions are still dry but less seasonal. As precipitation rises southward, tropical evergreen forest emerge, with intermissions of woody savannas around Belize and northern Guatemala, commonly associated with local edaphic conditions (^{Rzedowski, 2006}). The highland mountains are characterized by three basic types of vegetation, namely mountain mesophyllous forest, Quercus forest, and coniferous forest (^{Nixon, 2006}; ^{Rzedowski, 2006}).

Human impact has caused important changes in the natural vegetation mostly through deforestation. These changes have favored a widespread distribution of human-induced vegetation characterized manly by croplands and pastures in the entire region (^{Loveland et al., 2000}; ^{Rzedowski, 2006}). Anthropogenic vegetation in Mesoamerica is the result of both modern land uses and at least five millennia of extensive human occupation (^{Pohl et al., 1996}). The earliest cultivation of maize and the associated impact on vegetation have been reported between 5000 and 4500 years ago in Belize and Guatemala (Pohl et al., 1996; ^{Piperno, 2006}; ^{Wahl et al., 2006}). Crops of diverse sizes included maize (Zea mays), beans (Phaseolus spp.), squash (Cucurbita sp.), and papaya (Carica papaya). Also, the ancient legacy of Mayan agroforestry has been in use for 4000 years and it is widespread among present-day farmers. These activities include the Maya milpa that consists of annual rotation of crops with a series of managed arboreal species (e.g., Brosimum alicastrum, Bursera simaruba, Cecropia peltata, and Zwietenia macrophylla), which under abandonment lead to their dominance in the reestablished forest (^{Quintana-Ascencio et al., 1996}; ^{Ford and Nigh, 2009}; ^{Nigh and Diemont, 2013}). Besides the previously mentioned species, in the Maya milpa, species such as Acacia cornigera, Brosimum alicastrum, Bursera simaruba, Cecropia peltata, Zwietenia macrophylla, and Vitex gaumeri are also commonly found (^{Ford and Nigh, 2009}). Currently, new production systems such as cattle raising and other crops (e.g., sugarcane and coffee) are also widespread along the region (^{Sharer and Traxler, 2006}).

2.2. POLLEN AND HUMAN INFLUENCE DATA SET

The pollen data set was composed of modern sediment samples collected along the Maya region. Whereas part of the dataset came from previously published work (77 samples from the Yucatan Peninsula and adjacent mountains, ^{Correa-Metrio et al., 2011}), 48 new samples were collected, extending the sampling towards Guatemala, Honduras and Salvador (Appendix 1). The preparation and analysis of the new samples followed the protocols described by ^{Correa-Metrio et al. (2011)} to ensure homogeneity of the entire dataset. All pollen counts were transformed to percentages of the pollen sum to off set differences in sample size (^{Birks and Gordon, 1985}).

Human influence at the sampled sites was assessed using the global Human Influence Index (HII) dataset, which has a spatial resolution of 1 km² (^{Sanderson et al., 2002}). HII values range from zero for fully pristine areas, to 64 for totally urbanized regions (for details see ^{Sanderson et al., 2002}, and http://sedac.ciesin.columbia.edu/data/set/wildareas-v2-human-footprint-geographic).

2.3. DATA ANALYSIS

A Detrended Correspondence Analysis (DCA, ^{Hill and Gauch, 1980}) on the pollen dataset was used to explore the ability of pollen assemblages to reflect human-impact patterns at the community level. The DCA was meant to identify the ecological space represented by the pollen samples through the a priori interpretation of taxa ordination (^{Correa-Metrio et al., 2014}). HII values per taxon were calculated though averages of HII of the sites weighted by taxon-relative abundance.

Subsequently, HII per species and sites were compared with DCA taxa and sites’ scores, respectively, and their relationship was generalized thorough non-parametric locally weighted regression (^{Cleveland and Devlin, 1988}).

Given oversampling of some localities within the studied region, the original dataset was geographically resampled, using a 2x2 km grid, where only one sample was selected from each cell. The resampling was replicated 500 times. Additionally, the resampled dataset was submitted to a filter to avoid undesirable effects of rare species (only taxa with a minimum abundance of 1% in 5% of the samples were selected, after ^{Correa-Metrio et al., 2010}).

The magnitude, direction, and uncertainty of responses of individual taxa and community to the studied human-influence gradient were estimated though a Threshold Indicator Taxa Analysis (^{Baker and King, 2010}). The method finds values along the studied environmental gradient (HII in this case) where the largest community change occurs. An indicator value is calculated for each taxon at each candidate change-point along the human-impact gradient, and change-points with maximum indicator values are retained (for details see ^{Dufrene and Legendre, 1997}). Scores of the indicator values are subsequently standardized to obtain a Z-score for each taxon, such that positive and negative responses can be distinguished (^{Baker and King, 2010}). Community-level thresholds can be estimated as the summation of individual taxon responses represented by Z-scores. A bootstrap resampling is performed to evaluate statistical significance of Z-scores, which can be summarized as follows: (i) high purity when the Z-scores associated with the taxon are found within the same segment and show the same direction in at least 90% of the bootstrapping runs (1000 in total), and (ii) high reliability when at least 95% of bootstrapping were significantly different from a random distribution (p < 0.05). All analyses were performed using R (^{R Core Team, 2017}), especially packages TITAN 2.0 (^{Baker and King, 2013}) and Vegan 2.4 (^{Oksanen et al., 2017}).

3.1. POLLEN ASSEMBLAGES AND HUMAN INFLUENCE INDEX AT THE MAYA LANDS

With an average pollen count of 427 grains per sample, the 125 samples totaled 192 pollen taxa. The different kinds of highlands formations (mountain mesophyllous forest, Quercus forest, and coniferous forest) do not display a particular pollen signature. Therefore, the samples of these locations are grouped together as coniferous-Quercus forest (CQF). These samples are located at elevations above 1400 m a.s.l. (e.g., Atitlán, Chilangatorio, Esmeralda, and 5-Lakes) and show high percentages of Pinus, Quercus, Liquidambar, and Myrica (Figure 2). Croplands and pastures (CP) are characterized by a dominance of non-arboreal pollen (e.g., Amaranthaceae, Byrsonima, Poaceae, and Solanaceae) and high percentages of Cucurbitaceae and Zea, both taxa characteristic of croplands (Figure 2). These samples cover a broad altitudinal gradient (from 0 to 1200 m a.s.l.) and sites near to urban centers, roads, and farming areas (e.g., Amatitlán, El Espino, and Jucutuma) (Figure 2). Samples from seasonal vegetation (scrubs, deciduous, and semi-deciduous forests) do not differ in their pollen spectra, so are grouped as tropical seasonal forest (TSF). This group has the largest sample size, with 51 locations from lowlands and mountain depressions, and its pollen spectra is dominated by Acacia, Bursera, Vitex, and Mimosa. Pollen assemblages from locations classified as tropical evergreen forest (TEF), under 800 m a.s.l. (e.g., Yaxhá, Macanché, and Sacpuy), are dominated by Brosimum, Ficus, Trema, Moraceae, and Melastomataceae (Figure 2).

Figure 2. Pollen percentage diagram of mud-water interface samples from 125 lakes. Selected pollen taxa with at minimum presence of 1% in at least five samples are shown. Taxa were arranged according to their scores of Detrended Correspondence Analysis, while samples were organized according to the vegetation type of the area where the sample was recovered. 

Values of HII for the studied locations vary from 4 and 7 at the most preserved locations (Cobá and Chacan-Bata), to 43 and 48 at the most impacted localities (Calderas and El Espino). The higher values of HII are found at CP locations (26.3 in average), followed by CQF (21.8 in average), and lastly, TEF and TSF obtain similar HII (18 in average; Figure 2).

DCA scores produce a split of vegetation types along Axis 1 (eigenvalue 0.51, axis length 2.92) and Axis 2 (eigenvalue 0.24, axis length 3.1) (Figure 3C). Liquidambar, Myrica, Pinus, and Quercus show positive scores, while Brosimum, Ficus, Melastomataceae, and Moraceae display negative scores. Elements that usually characterize disturbances, such as Amaranthaceae, Ambrosia, Poaceae, Solanaceae, and Zea, are associated with scores near zero. Results from DCA ordination for the sites show a clear separation of vegetation types along the first two axes (Figure 3D). Samples from CQF are clustered in the positive end (quadrants I and IV), while CP locations show a wide distribution along the central part of Axis 2. Samples from TSF and TEF are separate in the negative side (quadrants II and III, respectively).

Figure 3. Detrended Correspondence Analysis (DCA) of modern pollen from Maya lands. A. Locally weighted non-parametric regression (loess) for human influence index as a function of DCA Axis 1 taxon scores. B. Loess regression for human influence index as a function of DCA Axis 1 sample scores. C. Ordination of taxa. D. Samples ordination with taxa selected from the Threshold Indicator Taxa analysis, negative (green) and positive (red) responses to the human influence CQF: Coniferous and Quercus forest; CP: croplands and pastures; TSF: tropical seasonal forest; TEF: tropical evergreen forest. 

A positive non-linear relationship between DCA Axis 1 taxa scores and their weighted HII is evident according to loess regression with low residual error standard (RSE, 2.3) (Figure 3A). Also, the loess analysis also suggests a breakpoint along the HII gradient between 20 and 22. The comparison between HII values and DCA Axis 1 scores per sites do not show a relationship with a loess-regression poorly adjusted (RSE = 8.3) (Figure 3B).

3.2. THRESHOLD INDICATOR TAXA ANALYSIS

The spatial resampling and the presence/persistence filter generate a reduced dataset composed of 112 samples and 76 pollen taxa. Through taxa indicators analysis, 20 of the 76 taxa are identified as pure or reliable indicators of human impact (Figure 4A). Negative responses (Z-) to HII are found for Mimosa, Vitex, Bursera, Ficus, Sapotaceae, Brosimum, Trema, and Cecropia, which show lower occurrence and abundance along a short length of human-influence gradient (average = 11). Eugenia, Mecardonia, and Acacia also decrease as HII increases but alonga wider length of gradient (average = 29). Nine taxa show positive responses to HII. Among them, Zea, Byrsonima, Pinus, and Myrica respond along short sections of the human-impact gradient (average = 13), while Solonaceae, Euphorbiaceae, Clusiaceae, Amaranthaceae, and Quercus do so along a wider gradient (average = 25).

Figure 4. Threshold Indicator Taxa Analysis. A. Inflection points (5^th, 95^th bootstrap quantile intervals) for significant pollen taxa along human influence gradient. Significant indicator pollen are taxa with IndVal p < 0.05, purity is greater than 0.90 and reliability is greater than 0.95 for 1000 bootstrap and 250 permutation replicates. Green and red circles represent change points associated with negative and positive responses, respectively. Note that Z scores are sized proportional to the magnitude of the response. B. Community changes for pollen taxa. Summation of individual taxon responses represented by Z- and Z+ scores. 

About the community responses along the human-influence gradient, the summation of Z-scores (sum(Z)) indicate a change-point at HII = 15 (5^th, 95^th quantiles of HII = 12, 20, respectively), having a sharp, well-defined peak in sum(Z-) scores (Figure 4B). Except for Mimosa, Vitex, Eugenia, and Cecropia, where significant Ztaxa have points changes outside the threshold. The sum (Z+) scores show a change-point at HII = 24 (5^th, 95^th quantiles of HII = 14, 30, respectively), but it is not considered a threshold response because of wide quantiles and poorly defined peaks in sum (Z+) (Figure 4B). Few taxa have well-defined synchronous change-points along the human-influence gradient (Figure 4A), indicating gradual changes in sum (Z+) scores with HII between 14 and 30.

4.1. MODERN POLLEN SIGNAL

Modern pollen assemblages are highly sensitive to local vegetation types, as previously reported for the Maya lands (e.g.^{Islebe et al., 2001}; ^{Domínguez-Vásquez et al., 2004}; ^{Bhattacharya et al., 2011}; ^{Correa-Metrio et al., 2011}). However, although vegetation of the studied region is represented by at least 10 different vegetation types, our pollen spectrum distinguishes four main associations: CQF, TSF, TEF, and CP (Figure 2). Human-induced vegetation such as croplands and pastures are widely distributed in the region, simplifying vegetation and diversity, and therefore blurring the ability of pollen to distinguish more-specific vegetation associations. Nevertheless, the low taxonomic resolution of pollen analysis is likely a structural problem to distinguish among vegetation types (^{Correa-Metrio et al., 2011}). In the lowlands, the pollen signal might be masked by genera shared by different forest types that are not distinguishable at the species level. In the mid-elevation and highlands, anemophilous taxa that have been widely recognized as problematic (e.g.^{Lozano-Garcia and Xelhuantzi-López, 1997}; ^{Correa-Metrio et al., 2011}; ^{Correa-Metrio et al., 2012b}) cause further masking of the pollen signal. Although a pollen-based, clear-cut distinction of forest types according to modern vegetation is not possible, the pollen signal reflects the continuous nature of the realized biological and environmental gradients.

Our results also produce a clear distinction of CP assemblages dominated by non-arboreal taxa such as Ambrosia, Amaranthaceae, Solanaceae, Zea, and Poaceae, a finding previously reported for the region (^{Bhattacharya et al., 2011}; ^{Correa-Metrio et al., 2011}). These taxa are also commonly found in large numbers within the CQF locations, which would suggest an important component of human-impacted vegetation within these areas. Nevertheless, CQF enrichment of these elements is probably also the result of the environmental hardships that vegetation experiences at high elevations through the entire year.

4.2. POLLEN ASSEMBLAGES OF THE MAYA LANDS: ECOLOGICAL AND HUMAN SIGNALS AT ASSEMBLAGE LEVEL

The DCA Axis 1 shows a clear separation of TEF and TSF on its negative side, while CQF is clustered towards the positive end (Figure 3D). Thus, DCA Axis 1 clearly reflects a temperature gradient, a finding previously reported for Central America and the Maya lowlands (e.g. ^{Bush and Colinvaux, 1990}; ^{Correa-Metrio et al., 2011}; ^{Correa-Metrio et al., 2012a}). CP samples obtained scores near zero in the DCA Axis 1 (Figure 3D), although they are not necessarily coming from intermediate temperatures. According to the loess regression, the relationship between pollen assemblages and HII values show a non-linear pattern. This is explained because taxa associated with the highest HII are arranged towards the center of Axis 1. Therefore, is probably the result of cosmopolitan taxa (e.g., Amaranthaceae, Ambrosia, Asteraceae, and Poaceae) acting as a gravity center of the ordination and therefore representing the transition between contrasting environments, in this case from lowlands to highlands (Figure 3C). Indeed, this transition zone likely reflects a simplifying effect of people on the landscape (^{Caballero-Rodríguez et al., 2017}), explaining the changes from complex tropical ecosystems to most simple temperate ecosystems in terms of its structure.

On quadrants II and III of the DCA, Axis 2 produces a clear separation of lowland systems into TEF (negative scores) and TSF (positive scores) (Figure 3D), implying a seasonality gradient previously reported for the region (^{Correa-Metrio et al., 2011}; ^{Correa-Metrio et al., 2012a}). At the center and the positive end of Axis 1, however, CP and CQF samples are widely distributed along Axis 2. This finding reflects two independent factors on these pollen spectra groups: (i) the widespread distribution of CP samples along the region in terms of the seasonality gradient, and therefore the lack of association of HII with DCA Axis 2; and (ii) the lack of a clear precipitation seasonality gradient along the CQF sampled sites (Figure 3D).

Although the temperature/elevation gradient led the DCA ordination, there is a signal of human influence in the vegetation suggested by the significant correlation between Axis1 taxa scores and their weighted HII (Figure 3A). However, loess regression on HII data as a function of Axis 1 site scores produced by the DCA does not show a pattern as clear as that of the scores of taxa (Figure 3B). Thus, pollen spectra reflect human influence on vegetation through a taxon approach, and the assemblage signal is rather weak.

4.3. HUMAN INFLUENCE INDICATION AT TAXON LEVEL

The threshold-indicator taxa analysis allows the identification of 20 pollen taxa whose presence and abundance are significantly associated with human influence (Figure 4A). Most of the elements that show a decreasing response to HII are arboreal and their distribution is confined to the lower end of the impact gradient (Figure 4A). The inflection point in the response of these elements is found below HII = 20, except for Cecropia (Figure 4A). According to the DCA, Mimosa, Vitex, Eugenia, Ficus, Brosimum, Acacia, and Trema, all negative responders to HII, are associated with TSF (Figure 3D). Alternatively, only Bursera and Sapotaceae are associated with TEF, pointing to the dominance of negative responders in the TSF and suggesting a high sensibility of these ecosystems to human impact. These forests are established mostly in the lowlands from the Yucatan Peninsula, where high seasonality probably potentiates the effects of human activities. Conversely, more stable climates through the year at the TEF locations might buffer human impact. Cecropia and Mecardonia do not show a clear grouping with a single vegetation type, but they are arranged at the TSF-TEF transition. Specifically Cecropia, a taxa widely recognized as a component of secondary forests (^{Rzedowski, 2006}; ^{Bhattacharya et al., 2011}), shows a negative response along HII gradient. This apparent discrepancy might indicate that, albeit tolerant to certain disturbance levels, this taxon is affected when the human impact gets intensified. Indeed, the inflection point of its response to HII is found at a rather high value (22, Figure 4A).

Taxa with positive responses to human influence are characterized by both non-arboreal and arboreal elements, and spread along a substantial portion of the studied HII gradient (Figure 4A). Zea, Solanaceae, Euphorbiaceae, Clusiaceae, Amaranthaceae, and Byrsonima are distributed along mid to high HII indexes, and according to the DCA are mostly associated with CP locations (Figure 3D). However, their central location on DCA Axis 1 suggests that they can be found in any vegetation association. These elements have been reported as disturbance indicators in modern and fossil pollen studies from Central America (^{Wahl et al., 2006}; ^{Bhattacharya et al., 2011}; ^{Correa-Metrio et al., 2011}). Myrica, Pinus, Quercus and Byrsonima, all arboreal to shrub elements from CQF (Figure 3D), show increases at low to moderate levels of human influence. Low intensity anthropogenic activities such as agroforestry systems, particularly in montane areas, can generate patchy open habitats, which can be rapidly colonized by opportunistic taxa (^{Cayuela et al., 2006}). However, Pinus and Quercus should be interpreted with caution because they are high pollen producers and long-distance dispersers (^{Mazier et al., 2006}; ^{Correa-Metrio et al., 2011}), implying a significant regional input that may not necessarily reflect local conditions.

Other taxa commonly identified as indicators of disturbance, such as Ambrosia, Asteraceae, Celtis, and Poaceae (^{Pohl et al., 1996}; ^{Bush, 2002}; ^{Correa-Metrio et al., 2011}), are not identified as significantly responding to human impact. However, they are clustered near the origin of DCA Axis 1 (Figure 3C), a region that our loess regressions identified as of maximum HII for both taxa and sites. The lack of significant response to the HII gradient through the threshold-indicator taxa analysis might be a result of their cosmopolitan distribution (^{Marchant et al., 2002}). Their association with human impact might be better reflected by their abundances than by their occurrence along the studied sites, and both attributes are taken into account by the method used.

4.4. THRESHOLD IDENTIFICATION

The summation of Z values to evaluate community-level responses show a substantial decline around HII = 15 (Figure 4B), suggesting almost synchronous response of all taxa that decline with the human impact. Responses to human disturbance of forests, insects, and lake biological communities showed similar thresholds for low and intermediate disturbance levels (^{Cardoso et al., 2013}; ^{Kovalenko et al., 2014}; ^{Rodrigues et al., 2016}). It is likely that once the level of human pressure exceeds a given threshold, important elements of forest cover are affected simultaneously, benefiting opportunistic species. Although we detected 7 of 11 taxa with thresholds near HII = 15, the uncertainty at the community level is rather high, pointing at the individual response of biological populations to disturbance. Thus, there is not a single threshold that explains the response of the entire community, but instead there might be a common level where the risk of taxa loss is high. On the other hand, asynchrony among positive responders is evidenced in the several change points (at 14, 18, 24, and 26 HII values) along the human-influence gradient, which produce a relatively weak aggregate signal of community response (Figure 4B). It is possible that the occupation of these taxa in human-disturbed areas depends on more complex factors; for example, the temporal structure of the disturbance.