Environmental DNA reveals arboreal cityscapes at the Ancient Maya Center of Tikal

Tikal, a major city of the ancient Maya world, has been the focus of archaeological research for over a century, yet the interactions between the Maya and the surrounding Neotropical forests remain largely enigmatic. This study aimed to help fill that void by using a powerful new technology, environmental DNA analysis, that enabled us to characterize the site core vegetation growing in association with the artificial reservoirs that provided the city water supply. Because the area has no permanent water sources, such as lakes or rivers, these reservoirs were key to the survival of the city, especially during the population expansion of the Classic period (250–850 CE). In the absence of specific evidence, the nature of the vegetation surrounding the reservoirs has been the subject of scientific hypotheses and artistic renderings for decades. To address these hypotheses we captured homologous sequences of vascular plant DNA extracted from reservoir sediments by using a targeted enrichment approach involving 120-bp genetic probes. Our samples encompassed the time before, during and after the occupation of Tikal (1000 BCE–900 CE). Results indicate that the banks of the ancient reservoirs were primarily fringed with native tropical forest vegetation rather than domesticated species during the Maya occupation.


SI: Archaeological Background
Within the Maya Lowlands lies an irregularly-shaped physiographic province known as the Elevated Interior Region (EIR; Fig.1), a karst area characterized by an acute lack of perennial surface water and almost no access to groundwater (1). Given a regional climate with a 5-month-long dry season, year-round occupation of the EIR by large numbers of Maya people and urbanization was dependent on their ability to capture and store enormous quantities of rain water during the rainy season. One of the first archaeological explorations of the northern part of the EIR by John Lloyd Stephens (and artist companion Frederick Catherwood) noted: "Among the wonders unfolded by the discovery of these ruined cities, what made the strongest impression on our minds was the fact that their immense population existed in a region so scantily supplied with water (2)." Catherwood illustrated several ancient reservoirs, including a hypothetical crosssection of one that had been cleaned out by an agricultural estate owner and included a number of chambers and wells constructed in the floor of the reservoir. Early explorers of the EIR, both scientific and economic, e.g., chicle harvesters and loggers, were utterly dependent on obtaining water during the dry season from ancient Maya reservoirs and were acutely aware of their association with ruined settlements (3)(4)(5)(6). The 19 th century reuse of ancient reservoirs was a common practice in the northern EIR (7).
The first modern archaeological excavations of ancient Maya reservoirs took place in the 1950s and 1960s as part of the University of Pennsylvania Tikal Project. Peter Harrison (8) published a short description of the Penn excavations in the Museum's journal Expedition, otherwise most of the Penn Project reservoir research remains unpublished.
Gary Gallopin and Vernon Scarborough used the University of Pennsylvania Tikal Project's detailed topographic and architectural maps to make a hydrological analysis of the reservoir system within central Tikal (9). Those analyses, along with the Lentz and Hockaday study (10) of Tikal's archaeological plant remains, formed the impetus for the creation of the University of Cincinnati Tikal Project which examined ancient Maya water, forest, and land use around Tikal. Excavations were conducted in 2009 and 2010 and laboratory analyses have been ongoing. Excavations were conducted in Palace, Temple, Hidden, Perdido, Corriental, and Pital reservoirs in and just south of the site center, as well as the Aguada de Terminos, Aguada Vaca del Monte, and Aguada Elmer in and around the Bajo de Santa Fe several kilometers east of the site center (11)(12)(13)(14). Among the findings was the discovery that the feature known as the Silting Tank or ("Spring Pool") situated topographically above the Temple Reservoir was constructed around a natural spring, likely in the Mid-to Late Preclassic period (ca. 521-216 BCE).
The Palace Reservoir included a system for periodic maintenance and waterlevel adjustment by employing stacked sluice portals in the Late Classic dam which encased a smaller Early Classic dam. Particulate analysis of sediments from several of Tikal's reservoirs indicate that the region received episodic volcanic ash fall (15,16). Geochemistry and genetic analysis of several reservoirs revealed that Temple and Palace reservoirs were highly contaminated with toxic cyanobacteria and mercury as water levels declined during Terminal Classic droughts in the 9 th century CE, whereas Corriental, Perdido, and Terminos reservoirs remained less affected (17).
The water complex that includes, in descending order, the Temple Reservoir and associated Spring Pool or Silting Tank, the Palace Reservoir, and the Hidden Reservoir, lies at the heart of Tikal surrounded by royal palaces and several major temple pyramids among other structures. The most elevated element in the complex is the so-called "silting tank," a name assigned by Scarborough and Gallopin (9) to a relatively small square tank perched above the Temple Reservoir and connected to it by a sluice. The assigned name was based on a model of water flow and sediment settling found in some other ancient Maya reservoirs within fluvial systems, e.g., at Kinal, Guatemala (18) and in other parts of the world. One problem with this interpretation at Tikal is that the "silting tank" is essentially perched at the head of a drainage that would have received very little fluvial input. Furthermore, almost all of its catchment would have been paved by the Early Classic period or earlier, hence, there was likely little sediment being deposited in this tank. Some sediment did accumulate beginning in the Late Preclassic as well as in the Classic period, with some indication that dredging may have taken place (11,12). Also arguing against the "silting tank" functioning as its name implies was the discovery of a natural spring within this feature; it began to flow when one of our excavations breached a heavy cap of clayey sediment (11,12). Thus, we suggest that this feature was a Spring Pool built to contain the water of the spring, a place that would have had tremendous symbolic value to Tikal and its rulers. Water sources were closely linked to the places of foundation for ancient Maya dynasties (19).
The Spring Pool disgorged water into the subadjacent Temple Reservoir, which also received runoff from surrounding paved surfaces. Sediments within the Temple Reservoir tank, superimposed over a thick clay liner above limestone bedrock, dated exclusively to the Late Classic period, probably indicating that the tank was dredged episodically until sometime in the Late Classic (11)(12)20). The Temple Reservoir and Spring Pool were separated from the Palace Reservoir by a limestone bedrock ridge (Fig.1). Discharge from the Temple Reservoir flowed through a channel at its northern end into the Palace Reservoir and was controlled by a coffer dam.
The Palace Reservoir was created by the damming of a natural ravine initially by a small dam in the Early Classic period, but later by the massive Palace Reservoir dam, a huge construction of rock, masonry, and clay probably anchored by bedrock projections at its north and south ends (12). Within the reservoir some steps or benches were created by a combination of quarrying bedrock, flagstone paving, and earthen embankments. Sediments within the reservoir included residual Preclassic soil overtopped by Late Classic to Terminal Classic sediments and later materials, a clear indication that the reservoir tank was periodically dredged or flushed during the Early and Late Classic periods. Well-preserved varved sediments uncovered beneath dam wall collapse in the Palace Reservoir exhibit thin bands of organic sediment, probably derived from leaves falling into the reservoir in the dry season interbedded with lenses of carbonate sediment probably derived from the weathering of plastered surfaces, representing successive dry seasons (20). Many Neotropical plants in this region are deciduous because of the dramatic seasonal fluctuations in water availability. While it is possible that some leaves may have blown into the reservoir from vegetated land lying south and downslope from the South Acropolis and Temple V, we believe it is more likely that these leaves originated from plants in the immediate vicinity of the reservoir system, especially the ridge of unbuilt land between the Palace and the Temple Reservoir and Spring Pool, as well as the south flank of the Palace Reservoir.
Excavation operations 6O, 6L and 6K in the Palace Reservoir were formed by a series of three pits linked together into a trench (Fig. 2). This trench was located some 10 m west of the present-day interior dam of the Palace Reservoir and ran south from the current topographic low onto a gradual rise in the reservoir floor. Four distinct features are visible in the excavation profile. First, at the lowest part of the trench near its north end is a narrow channel incised into bedrock and filled with a very dark organic silty clay, likely an aquic soil formed with a pooling point in the ravine channel (20). Stratigraphically above this basal clay is a zone of light gray silty clay formed by sediment deposited in the stagnant pool behind the reservoir dam; ceramic sherds embedded within this clayey sediment indicate that it was deposited during the Terminal Classic period (850-900 CE) grading upward into the Early Postclassic (900-1100 CE). Interfingering with the clayey reservoir pool sediments are a series of coarse sediments (mainly sand and gravel sized particles) in the form of thin debris flows that entered the reservoir from the south. These flows likely originated from collapses along the very steep southern embankment of the reservoir. Several larger debris flows are visible on the lidar-derived image of the Palace Reservoir (Fig. 1). Charcoal embedded within one of the lower debris flows produced a calibrated radiocarbon date range of 610-680 CE consistent with Late Classic ceramic sherds recovered from multiple flows. The material in the flows can be considered intrusive older construction fill originating in the decaying south wall of the reservoir. Notably, a section of the interior wall of the Palace dam is also known to have collapsed into the reservoir pool then subject to a buttressing repair (12,20). The uppermost stratum revealed in the trench is the modern soil formed during the Postclassic period and later sediment and organic matter, reflecting a period at which time the reservoir was no longer effectively impounding water.

SI Materials and Methods
Because eDNA is typically highly degraded in archaeological contexts (21), obtaining long sequence reads generally is not possible. The 18S rRNA gene has been used elsewhere to study eukaryotic biodiversity (22). Highly conserved flanking regions allow for development of universal primers to amplify the informative variable sequences in rRNA. An added advantage is that ribosomal sequences are repeated in the genome, which increases the probability of amplification. Differences in 18S rRNA genes, however, have not proved to be sufficient to resolve plant taxa. In our previous studies (17), fungal species were abundant in reservoir sediments. Fungi can metabolize decaying matter and can proliferate in the sediments of the reservoirs, increasing their biomass compared to the plants. Extremely deep sequencing would have been needed to detect the rare plant genes under these circumstances. In this study, we primarily targeted plastid genes, unique to plants, from multiple regions and protein coding genes (rbcL, rp110, ycf1 and matK), non-coding spacers (trnT-trnL, trnL-trnF, and trnH-psbA), and the internal transcribed spacer regions (nrITS) of nuclear ribosomal DNA. Similar sets of target genetic markers have been successful in previous ancient vegetation studies (23)(24)(25).
Sediment samples from stratigraphic layers within four Tikal reservoirs and several aguadas (ponds) were collected during excavations in 2010. To avoid contamination from modern soils, and eliminate cross contamination from other strata, we collected samples in a column (26) from a selected characteristic stratigraphic profile in each pit or excavation unit. This was accomplished by first shaving each profile with a sharp, clean trowel where the samples would be collected. Samples were taken in 10 cm increments starting at the bottom of the pit and working toward the top. Sediment samples destined for molecular and pollen analyses were extracted from the freshly shaved wall surface, placed in sterile plastic bags (Whirl-Packs) and labeled. In addition, flotation samples (2 liters each) were collected from the same profile after we finished collecting the pollen and molecular samples. This approach gave us the benefit of microremain and macro-remain plant data from each context. The macro-remain samples were processed by water flotation (27) and the pollen samples were sent to a laboratory in the US for pollen extraction. At the pollen lab, 2 g of sediment were removed from each sample in a laminar flow hood with a sterile spatula then the bags were immediately resealed. Results of the pollen and macrobotanical analyses have been published previously (27)(28)(29). Although we did have success extracting and identifying significant amounts of pollen from the ponds (aguadas) in the seasonal swamps (bajos) near Tikal, we were unsuccessful at doing so from the site core reservoir samples.
After arrival at the University of Cincinnati Paleoethnobotanical Laboratory, 30 samples designated for molecular genetic analysis were placed in a -80° C freezer until processing could begin. Just prior to DNA extraction, samples were thawed to 4° C then inserted, under sterile conditions, into tubes with glass beads then sealed prior to homogenization in a bead-beating machine.
DNA was extracted from sediment samples using DNeasy PowerSoil Kits (Qiagen) as described previously (17,(30)(31). In brief, archaeological samples were extracted in a DNA clean lab facility dedicated to the study of cyanobacteria at the University of Cincinnati Department of Biological Sciences. Because of the general focus of this lab on cyanobacteria, it was highly unlikely that our samples would be contaminated with plant DNA from other experimental efforts. Nevertheless, rigorous protocols to avoid contamination were employed, including standardized workflow procedures, the wearing of personal protective equipment, and the preparation and processing of negative controls. Specimens underwent DNA extraction in batches of four each. The Uaxactun garden sample was extracted separately from the other samples discussed in this study and the hood was sterilized with bleach and 70% ethanol between each extraction. As always, established lab protocols (32) designed to avoid cross-contamination were followed scrupulously. This method was employed to successfully extract DNA from Tikal archaeological samples. The DNA content of samples was initially evaluated by Qubit Fluorometer Quantitation.
In order to process the fragmented genetic data recovered from the DNA extractions, we employed the services of RAPiD Genomics LLC (Gainesville, FL). Upon arrival of our 29 archaeologial extracts at the RAPiD Genomics facility, our samples were analyzed, under sterile conditions, using PicoGreen technology (Molecular Probes, Eugene OR), a commonly used assay for fluorescence enhancement coupled with spectrophotometry for accurate dsDNA quantitation (33). Control DNAs used in this process were from commercial sources. In the spectrophotometer each sample was read against a blank 100 µl quartz microcuvette containing only TE buffer. From these assays we learned that 8 archaeological samples and our control sample from the Uaxactun garden contained sufficient DNA for whole genome amplification and library preparation. The results of this assay can be found in Table S10.
In addition, RAPiD Genomics was enlisted to design genetic probes that could capture common plant genes across a variety of taxa. To help create the probes, we provided sequences from GenBank for the regions nrITS, matK, ndhF, psbA-trnH, rbcL, rpl10, trnL-trnF, trnT-trnL, and ycf1 for 68 species that are known, or suspected, to have been used by the ancient Maya (see Tables S1 and S9). These regions were selected because they have been well documented for variability at the genus and/or species level, and because they are well represented in GenBank. From these sequences, RAPiD Genomics designed probesets that they used to enrich DNA from our sediment extractions using their Capture Seq protocol prior to high-throughput Illumina sequencing. Following are the details of the probe design, the library preparation and the bioinformatics analysis performed.
The probes, which focused on plastid and nuclear genes, were based on the sequences we provided and were arranged in a 3x tiling array. All possible probes were designed in silico on the candidate genes (common plant genes across a variety of taxa), with start-end coordinates provided as target subsequences. From the total of all possible probes within the candidate sequences, no filters were applied to select a set of 11,391 probes that were used for hybridization. Briefly, DNA was sheared initially to a mean fragment length of 400 bp; fragments were end-repaired, followed by incorporation of unique dual-indexed Illumina adapters and PCR enrichment. Sequence capture and library preparation were performed utilizing their high-throughput workflow with proprietary chemistry. Following the probe capture process, only 6 of our 8 remaining archaeological specimens and the control sample were shown to contain adequate eukaryotic DNA to merit sequencing. Of the 6 archaeological samples with eukaryotic DNA, 4 were chosen because of their high archaeological significance (the other 2 were from A horizons) leaving us with five samples for sequencing including the Uaxactun control sample.
Regarding negative controls, these were carried out through the full sample processing pipeline alongside commercial samples as positive controls. These controls were evaluated for potential contamination at routine QC checks, which included DNA quantification and other proprietary measures. In general controls, both positive and negative, were included in library preparation to ensure successful library construction, these control samples typically are evaluated following library preparation and are typically not sequenced. No evidence of contamination or other challenges were observed via the positive or negative controls for this project.
Parenthetically, the five remaining samples provide some indication about the possible effects of degradation of DNA. The garden sample had 251 identifiable DNA sequences while the four archaeological samples had an average of 32.5 sequences per sample so obviously there were degradation factors at work, but fortunately there was still enough DNA preserved in the archaeological sediments to generate useful data.
The five remaining samples with adequate amounts of eukaryotic DNA and high archaeological interest were pooled equimolar and sequenced using HiSeq 2x150. Raw sequencing reads were trimmed using Trimmomatic software to remove any sequencing adapters. (http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomati cManual_V0.32.pdf) Each sequencing run was performed in triplicate and each run included positive and negative controls resulting in ~500 total runs. Following quality control, merging of contigs, and removal of chimeras and tiny fragments, we recovered DNA sequences with an average length of 253 bp. At a sequence identity of 97%, we recovered 380 Operational Taxonomic Units (OTUs). The taxa-specific gene abundance determined from these OTUs aided in discerning the vascular plant vegetation surrounding both reservoirs.
A sliding window with a minimum quality score of 30 and a minimum accepted length of 50 bp were the parameters used. Trimmed reads were aligned to the target sequences with the Burrows-Wheeler alignment algorithm (34) using default settings. Reads that mapped to the target sequences were then assembled using SPAdes http://cab.spbu.ru/files/release3.12.0/manual.html). The resulting assemblies from SPAdes were presented to the authors for the identification of sequences (nodes).
To identify the plant DNA sequences generated by Rapid Genomics, we used the Basic Local Alignment Search Tool (BLAST) algorithm and the GenBank nucleotide database of the National Center for Biotechnology Information (NCBI). BLAST hits were initially filtered by Bitscore and Evalue to include only the best hits for each of the assembled sequences for every sample. BLAST hits were further filtered (35,36) to eliminate hits involving plants that were not native or likely domesticates (37) in Precolumbian Guatemala. To make systematic use of the BLAST algorithm, we closely followed the decision tree presented in Fig. S1 to interpret results consistently and conservatively.

Fig. S1
. Decision tree for interpreting National Center for Biotechnology Information (NCBI) BLAST searches. Our eDNA gene sequences were compared to the NCBI database and their BLAST program evaluated the closest fit for our gene sequences by generating a Bitscore and Evalue. Generally, the higher the Bitscore the closer the match of our unknown sequences to sequences from plant species stored in the database. In numerous cases, only one plant species would be listed with the highest Bitscore. If the plant identified was native to Guatemala (35,36) or was a known New World cultigen (37), then our identification was clear and definitive. Because gene sequences can be conserved among even distantly related plants, however, it was not uncommon to find more than one species with a top ranked Bitscore. When this occurred, we followed the decision tree outlined below. Note that in all cases, we took a conservative approach and assigned a broader taxon if there was any ambiguity.

One top ranked hit
More than one top ranked hit         Table S6. Accelerator mass spectrometry (AMS) radiocarbon ( 14 C) dates for Tikal reservoirs and related contexts discussed in this text. AMS dates were obtained from Beta Analytic (Miami, FL, USA) and the National Ocean Sciences Accelerator Mass Spectrometry Facility (Woods Hole, MA, USA). Data provided include AMS radiocarbon sample composition (SOM = soil organic matter), provenience, depth in cm, measured radiocarbon years before present (BP), and calibrated age at two sigma margins for error. Samples were collected from wet and dry cores and excavation profiles. These data were extracted from larger tables published previously (17,40).        Table S9. Plant species and marker sequences (Genbank) used to create the genetic probe.