Quantifying Ancient Maya Land Use Legacy Effects on Contemporary Rainforest Canopy Structure

: Human land use legacies have significant and long-lasting ecological impacts across landscapes. Investigating ancient (>400 years) legacy effects can be problematic due to the difficulty in detecting specific, historic land uses, especially those hidden beneath dense canopies. Caracol, the largest (~200 km 2 ) Maya archaeological site in Belize, was abandoned ca. A.D. 900, leaving behind myriad structures, causeways, and an extensive network of agricultural terraces that persist beneath the architecturally complex tropical forest canopy. Airborne LiDAR enables the detection of these below-canopy archaeological features while simultaneously providing a detailed record of the aboveground 3-dimensional canopy organization, which is indicative of a forest’s ecological function. Here, this remote sensing technology is used to determine the effects of ancient land use legacies on contemporary forest structure. Canopy morphology was assessed by extracting LiDAR point clouds (0.25 ha plots) from LiDAR-identified terraced ( n = 150) and more closed, more vertically diverse forests. These human land uses abandoned >1000 years ago continue to impact contemporary tropical rainforests having implications related to arboreal habitat and carbon storage.

more closed, more vertically diverse forests. These human land uses abandoned >1000 years ago continue to impact contemporary tropical rainforests having implications related to arboreal habitat and carbon storage.

Human Land Use Legacies
One interface between the disciplines of landscape archaeology and landscape ecology is the study of land use legacies. Understanding the consequences of human alteration of landscapes has become increasingly important as human populations continue to expand and natural landscapes are continually transformed [1,2]. Just as modern human societies alter landscapes, so did prior ones [3]. Land use history (e.g., forest clearing, agricultural regime, abandonment) directly influences both the biotic (e.g., presence of novel species assemblages) and abiotic (e.g., changes in soil nutrients, hydrology, topography) environments [4]. These changes may lead to long-term effects on contemporary measures of biodiversity (i.e., composition, structure, and function) that may impact current services (e.g., habitat, carbon sequestration) provided by the ecosystem (Figure 1).

Figure 1.
The interacting and cascading effects of land use legacies on modern forested systems adapted from [4]. Regions enclosed by dashed outlines include some parameters that can be measured with airborne LiDAR. Blue and red outlines represent historic-archaeological (past) influences and contemporary (present) factors, respectively.
While deforestation is still occurring at an alarming rate throughout many parts of the world, afforestation is occurring in other regions as former agricultural fields are abandoned [5][6][7]. By one estimate 50%-80% of the forests in New England are on former agricultural fields [8]. Locations that have traditionally been thought of as sparsely populated and fairly pristine prior to European colonization, such as the Americas, are now known to have been significantly altered by ancient humans [9,10]; hence, the concept of "virgin" forests is now relatively obsolete [11,12]. Evidence suggests that vast tracks of forest had been cleared for agriculture or burned by indigenous humans prior to western arrival in the New World [13][14][15]. Given the current rate of land alteration and the extensiveness of past alteration by ancient people, the persistence and implications of these land use legacies have become important questions in ecology [16][17][18].
The time since a land use has been abandoned and allowed to revert back to a natural state is a determinant of the type and extent of the legacy present [1,2,10]. In European forests, agricultural legacies have been shown to endure over a millennium after abandonment, with species composition in secondary forests that regenerated over former agricultural fields lacking native species found in ancient forests [13,19]. Another example of ancient land use legacies continuing today is found in Central America with the Maya civilization. A large percentage of forests in Central America are secondary forests that have regenerated after the collapse of the Maya civilization [10,20]. Studies have suggested these forests contain trees species that were of economic importance to the ancient Maya, which may persist as remnants of the Maya orchard-gardens [21]. In addition to species introduction by the Maya, disturbance from agriculture and settlement constructions may have altered the landscape in ways that provided optimal habitat for specific, limestone loving species, such as the Ramon tree (Brosimum alicastrum) [22,23]. Fingerprints of ancient land use may also be found in forest structure (basal area, biomass, canopy height canopy closure, etc.) independent of species composition [24,25]. Forests may take between 60-200 years to regenerate to their previous canopy height [26,27]. However, structure has been shown to be influenced by topography [28] and soil nutrients [29], which may have been irrevocably altered by past land use, making it impossible for forests to regenerate to pre-land alteration conditions [2,17].

Forest Canopy Structure
The structure of a forest offers key insights into ecosystem function and biodiversity. It has been directly correlated to a number of important ecological measures, such as biodiversity and carbon stocks, i.e., aboveground biomass (AGB) [30,31]. By describing spatial patterns we can determine processes at work in a landscape [32]. Forest structure can be directly tied to species diversity, with a more heterogeneous structure often indicative of not only higher floral diversity, but also high faunal diversity as well [33][34][35]. The structure of forests is often used as a proxy for determining live AGB, an important factor in identifying sources of carbon sequestration [28,36]. Forest structure relates to the health of an ecosystem, where a higher degree of heterogeneity is suggestive of a more dynamic forest (i.e., disturbance, gaps, recruitment) [37,38]. Understanding the factors that drive forest structure can enable a better understanding of general ecosystem function [39][40][41].
A number of factors are known to affect forest structure, among them natural and human disturbance [42], climate, soil type and nutrients, and topographic position [29,43]. Topographic position, which include measures of relief (or slope) and elevation, have been shown to have a strong effect on forest structure. Forest canopy height and AGB have been directly correlated to topographic relief, with high AGB and greater canopy heights in lower slopes and valleys compared to steeper slopes which are typically associated with lower AGB and lower canopy heights [28,44]. Canopy openness has also been shown to increase on steeper slopes compared to low slopes and valleys [45].
While topographic position appears to be a strong factor in influencing forest structure, another important, if indirect, factor to consider is disturbance.
While natural disturbances such as floods, fire, and wind damage can drastically alter the structure of a forest [30,46], human disturbances can also have a significant impact [24,25,47]. However, for many human land uses, regenerating forests attain levels of vertical structure very close to original levels within 50 to 200 years [25,27]. Other forms of human land use may have longer lasting impacts, depending on the intensity of former land use [2]. Land use that alters the topography of the landscape has the potential to permanently alter the structure of the forest that regenerates, so long as the change in topography remains. One such land use can be found at Caracol, Belize, where the Maya extensively terraced across the hilly landscape. The terraces were constructed between A.D. 650-900 and are still in place today [48,49]. The terracing has essentially created fine-scale, level areas along hill slopes, effectively transforming the micro-topography of the region. Thus, Caracol provides the opportunity to investigate the extent to which human-altered topography continues to impact the ecosystem 1000 years after forest reestablishment.

Objectives
Though archaeologists typically wish to ignore or remove the vegetation that obscures the view of their ancient landscapes [50] and ecologists often wish to discount the influence of historic events on ecosystem composition, structure, and function [2], this study exploits the ability of LiDAR to simultaneously map forest canopy and below-canopy archaeological surfaces. Specifically, the goals of this research are to use LiDAR to (1) delineate the location of Maya agricultural terracing and areas where there are no detectable terraces and (2) quantify structural attributes of the forest canopy (i.e., height, openness, and vertical diversity) at Caracol. Furnished with this information, we can (3) assess the extent to which the ancient Maya agriculturally engineered (i.e., terraced) landscape impacts the LiDAR-derived forest structure of today.

Study Area
The Caracol Archaeological Reserve is located in west-central Belize along the Guatemalan border. The reserve resides across a karst landscape consisting of alternating hills and valleys, with elevations ranging between ~310 to 720 m asl. The Caracol study site encompasses an area of ~200 km 2 ( Figure 2) and is surrounded by the Chiquibil National Park. It represents a massive Maya city with extensive monumental architecture whose extent has been documented through LiDAR [51,52]. The dense subtropical moist forest of Caracol reaches average heights of 22 m and leaf off never exceeds 20%, with lowest leaf area index occurring during the end of the dry season (end of March, beginning of April) [42]. Recent human and natural disturbance of the forest has included hurricanes in 1968 and 1973 [53] and encroachment from Guatemala. These incursions entail the clearing of large swaths of forest along the border and targeted removal of desirable tree species throughout the reserve, primarily Swietenia macrophylla (Mahogany) and Cedrela odorata (Cedar) [54].

Figure 2.
Hillshaded DEM of the Caracol study site with 300 sample plot locations randomly distributed across terraced and non-terraced areas. Samples, which were stratified by land use type (terrace and non-terrace) and by slope (low, medium, high), were not collected near the Guatemala border, where there is ongoing deforestation (shown in pink).
Ancient Maya land use is ubiquitous throughout the study site. Thousands of residential and ceremonial structures are found throughout the ~200 km 2 area, as well as ~67 km of causeways radiating out from the epicenter (i.e., central metropolis area with a high density of temple, palace, and residential structures). The Maya built extensive limestone terraces throughout the landscape for agricultural purposes, growing maize and other crops, such as beans and squash [48,55,56]. These terraces have been reported to cover as much as 80% of the study site and can be found along hill slopes and in valleys. The ability of airborne LiDAR to geolocate accurately ground surveyed terraces hidden below the forest canopy at Caracol has been previously documented [52,57,58]. From our visual analysis and digitization of the hillshaded digital elevation model (DEM), conspicuously terraced areas represented 28.7% of the landscape while conspicuously non-terraced areas represented 35.7%. Terraces at Caracol have been shown to increase soil depth along hill slope by up to 0.6 m [48]. Studies of other Maya terraces have shown that terraces effectively prevent soil erosion along slopes, increase soil depth, and increase moisture retention [59,60].

LiDAR Data Collection, Extraction, and Sampling
In April 2009, LiDAR was flown over the Caracol study site [61]. The LiDAR survey used an Optech GEMINI Airborne Laser Terrain Mapper (ALTM) mounted on a twin-engine Cessna Skymaster. During the 9.5 h of laser-on flight time, a total of 122 flight lines were flown ~800 m above the ground surface, 66 in a North-South direction and 60 in an East-West direction, allowing optimal penetration through the dense canopy. The swath width was 520 m and flight lines were placed 260 m apart, insuring a 200% overlap. The survey yielded ~20 points per m 2 , for a total of 4.28 billion measurements, of which 295 million were classified as ground returns (1.35 ground return points per m 2 on average). The data were processed by NCALM (National Center for Airborne Laser Mapping); the final product included LiDAR point cloud files and a bare earth DEM with a 1 m horizontal resolution (see [62] for details explaining the DEM creation).
Three hundred non-overlapping 0.25 ha circular plots were randomly placed across the entire Caracol DEM using Hawth's tools [63] to the east of a ~2.5 km-wide area near the Guatemala border which has been subjected to illegal clearing and selective logging [54]. Sample plots were stratified using terraced (n = 150) and non-terraced (n = 150) land use layers in ArcGIS that had been carefully digitized off the DEM (Figure 3). Other than terraces, these areas excluded areas with other Maya structures (e.g., causeways, residences, monuments). Samples were further stratified by low (0°-10°; n = 50), medium (10°-20°; n = 50), and high (20°-40°; n = 50) slopes to insure an even representation of slopes. Using FUSION LiDAR analysis software [64], LiDAR point returns were extracted for all elevation layers for the Caracol DEM and slope, aspect, and elevation values were obtained for each of the slope classes in the sample. ArcGIS 9.3 Spatial Analyst was used to generate slope, aspect, and plots. The Density Metrics function in FUSION permitted us to sample the 0.25 ha LiDAR point clouds by extracting slices associates with height ranges. Returns were sampled in 3 m height slices and sorted into bins from the forest floor to the canopy top ( Figure 4) [65]. The final FUSION product yielded a data set of vertically binned height returns for each plot. All plots were sampled using a fixed grid with 10 × 10 m cells and all cells were summed for each height class. The proportion of total points was calculated for all height bins within each plot for the final height bin data set.

Statistical Analysis
The following forest structural characteristics were calculated using the LiDAR point height bins for each sample plot: average canopy height, canopy openness, and a vertical diversity index. Average canopy height was calculated by taking the average of the maximum heights for each 10 × 10 m window in the 0.25 ha plot. The measure for canopy openness was calculated by taking the ratio of ground returns divided by the sum of all other returns (ground returns/all non-ground returns) [66]. Jost diversity index was calculated for the proportion of returns that had been binned into 3-m height classes for each plot and used as an approximation of vertical diversity [67]. Permutational multivariate analysis of variance (PerMANOVA) was used to test for significant differences in forest structure between non-terraced and terraced areas. PerMANOVA allows for a non-parametric analysis of variance, which is ideal for non-normal data sets. PerMANOVA first calculates a distance matrix using a user selected distance measure. The test statistic (F-ratio) is calculated directly from the distance matrix. Permutations are then used to generate p-values to determine significance [68]. The Adonis function in the vegan package in R was used to perform the PerMANOVA analyses in our study [69].
For the first analysis, land use (terraced vs. non-terraced) was used as our factor, while slope, aspect, and elevation were held as covariates. A distance matrix (response variable) was calculated from the measured forest structural characteristics (i.e., average canopy height, canopy openness, and vertical diversity index) using the vegdist routine in the vegan package [69] for the R statistical software program. Vegdist uses the Bray-Curtis (Sorenson) distance measure, which is appropriate for community ecology datasets [70]. In a second PerMANOVA, we used a distance matrix calculated directly from the point return height bin data set for each plot. For the second analysis, land use was retained as our factor, while slope, aspect, and elevation were held as covariates. A total of 9999 permutations were used for each analysis.

Results and Discussion
PerMANOVA results from the LiDAR derived forest measurements indicated significant differences between terraced and non-terraced land uses (p = 0.016). Slope was also highly significant in explaining the variation between samples (p = 0.001), as was elevation (p = 0.001). The interaction of slope and aspect (p = 0.026), as well as the interaction of slope and elevation (p = 0.001), were also considered significant in explaining the variation seen between samples (Table 1). Using the LiDAR height bin distance matrix, PerMANOVA revealed a significant proportion (p = 0.005) of the variation was explained by land use categories, i.e., terraced vs. not-terraced. Slope was again found to be significant in explaining the variation in the proportion of points found in each height bin (p = 0.001). In this analysis, aspect did not contribute significantly to the explanation of variation. In addition, elevation did not explain a significant proportion of the variation, but the interaction effects of slope: elevation and land use: slope were significant in explaining the variation of height bin distributions (Table 2). When analyzed together as a distance matrix, the forest structure variables (canopy height, canopy openness, and vertical diversity index) were significantly impacted by the presence of terraces. When the individual trend lines of our three forest structural variables are graphed ( Figure 5), a few patterns emerged. On non-terraced plots there is a trend of vertical diversity decreasing as the slope increases; however, this decrease occurs at a slower rate for terraced land use. A similar pattern exists with average canopy height, with canopy height decreasing as slope increases, but this decrease occurs at a slower rate on terraced land. Canopy openness increases as slope increase, with canopy openness on terraces increasing at a slower rate compared to non-terraced land. At high slope conditions average canopy height and vertical canopy diversity are the most pronounced.
Histograms generated from the LiDAR height bin data showed a similar trend, that is, non-terraced bins showed a distinct gradient among low, medium, and high slopes, with medium slope point values falling between the low and high slope. However, on terraced land there is a trend of reduced variation among slope categories. Within a height bin, different slope categories do not show as steep "stair step" pattern in terms of LiDAR point returns ( Figure 6). This difference was most pronounced in the 0-3 m, 21-24 m, 24-27 m, and 27-30 m height bins. Another effect of terraces that can be deduced from this figure is an upward shift in the height of median energy (i.e., HOME) [66] from non-terraced to terraced plots which is indicative of higher AGB.
The shaping of forest structure by topographic variables may be explained by differences in tree species and the effect of topographic position on tree growth. Observations from previous studies have noted that terraces contain flora, such as palms, that are typically found in forest valleys [71]. A sister study [72] showed how tree species composition varies across terraced and non-terraced areas, with terraces acting as a type of environmental "bridge" between slope and non-sloped areas. Prior studies in Mesoamerican forests have shown topography to strongly influence the composition of tree species, with forests in valleys, along slopes, and along ridges forming distinct tree communities [73,74]. Tree species with specific topographic requirements have a direct impact on forest structure. Furthermore, though tree morphology is fairly plastic, different species with their specific architectural constraints, to a certain degree, may contribute uniquely to the general canopy structure [31].

Figure 5.
Relationships of LiDAR-derived values of (a,b) average canopy height (c,d) canopy openness, and (e,f) vertical canopy diversity from terraced and not terraced plots as a function of topographic slope. Significant differences between terraced and non-terraced areas in average values (detected with a t-test) for low (0°-10°), medium (10°-20°), and high (>20°) slopes are designated with p-values.
Variation in species composition alone may not explain the difference observed in forest structure. Edaphic factors related to topography, i.e., nutrient, water, and light availability, directly influence tree growth [28,29]. Differences in forest structure over terraces can be explained by altered edaphic conditions [55,56,59]. Steeper hillsides typically have thinner soils and reduced water availability; terraces increase soil depth and water availability, which impacts forest structure [44,75]. Terraces also retain nutrients that would otherwise be leached out of the soil [59,71]. This study identified a significant interaction between slope and elevation on canopy structure (Tables 1 and 2). As position upslope increases, nitrogen decreases, which can have an effect on the growth of trees, resulting in reduced basal area and canopy cover [29]. Beyond the implications for human land use legacies, this study illustrates how relatively small changes in topography can result in long-lasting changes to rainforest structure [44]. These results indicate that the terraces constructed by the Maya over 1000 years ago continue to influence the forest structure at Caracol. The terraces have significantly altered the micro-topography of the terrain, resulting in a corresponding echo in forest structure. While topography (i.e., slope, elevation, and aspect) appears to act as the driving factor in forest structure variation, the addition of terraces dilutes the topographic gradient, reducing the variation in forest structure from low-lying valleys to highly-sloped hills.

Conclusions
Agricultural terracing transforms hillslope topography; these oft-hidden archaeological features can last for millennia and are readily delineated with LiDAR [52,57,58]. While likely engineered for soil and water conservation, here we used airborne LiDAR to reveal the persistent, unintentional impacts that ancient Maya terracing has had on the structure of the tropical rainforest that regenerated after agricultural abandonment [60]. This study showed that the presence of terraces dampened the effect of topographic slope on structural variables (i.e., height, openness, and vertical diversity) that define a forest canopy. With steeper slopes, terraces yielded ~8% taller, ~20% more closed-canopy, and ~7% more vertically-diverse forests. As canopy height decreased with an increase in slope on non-terraced lands, height differences were muted with terracing (Figure 5a,b). Canopy openness, a measure related to understory light penetration, increased as slope increased; however, on terraced areas, openness varied less as a function of slope (Figure 5c,d). Vertical canopy diversity, an estimate of vertical spatial heterogeneity that has been correlated to avian diversity [33], followed a similar pattern as average canopy height; its variance decreased along the topographic gradient with terracing (Figure 5e,f). Terraces mediated the effect of slope on the proportion of LiDAR returns in the vertical, 3-m height bins, i.e., when terraces were present, forest structure did not differ as much with slope as when terraces were absent ( Figure 6). Though not directly measured, these canopy changes most likely correspond to higher levels of aboveground biomass or stored carbon [66].
The ability of LiDAR to concurrently measure ground and canopy surfaces elucidated this heretofore unknown, long-term anthropogenic impact on these forests. After Caracol was abandoned ca. AD 900, ancient Maya agricultural practices continue to influence contemporary forest structure. These effects should have consequences on ecosystem services ( Figure 1) related to arboreal habitat and carbon sequestration abilities of this biodiverse, high-biomass region, which is part of the Mesoamerican Biological Corridor [76]. Though this study only examined terraces at Caracol, these anthropogenic features are found throughout the Maya region [59,77], including the Petén in Guatemala [78] and the southern Yucatán in Mexico [20]. With the prevalence of relic terraces throughout Central American landscapes and elsewhere [79][80][81] comes the understanding that past agriculture practices have left lasting legacies on ecosystem structure and function. LiDAR demonstrated the synergy between landscape archaeology and landscape ecology. This approach can be extended across other large expanses of temperate and tropical forests where humans have left significant fingerprints from past land use [3,8,12,81] that currently lay obscured below vegetation. Moreover, this remote sensing study gives us pause over our current land use practices [82] and reminds us that today is tomorrow's past.