Introduction

Modern maize (Zea mays ssp. mays L.) is one of the most important and productive global food crops¹. Prior to European contact in the Americas, maize had been adapted for a wide range of environments. Its versatility is clear from its extensive geographic distribution today, ranging from 58° north latitude to 40° south latitude, and growing below sea level as well as 3,600 meters above sea level in the Andes². In Bolivia there is empirical evidence that small maize varieties can grow at altitudes of 3,800-3,900 meters above sea level (masl) in some areas of the Lake Titicaca Basin³, as well as in other sheltered areas of Oruro and Potosí, although with lower levels of production^{4, 5}Studies of archaeological samples and modern molecular data indicate that maize originally evolved in Mexico from the flowering species of grass, teosinte (Zea mays ssp. parviglumus) ∼9,000 years before the present (cal. BP) ^{6, 7}and spread into South America by ∼7,000 cal BP^{8, 9} from the lowland and Eastern North America by ∼2500 cal BP¹⁰. During this domestication process there was hybridization between teosinte and mexicana (Zea mays ssp. mexicana)¹¹. The dispersal of maize to the rest of the world followed European colonization of the Americas in the 15^th and 16^th centuries AD^{2, 12}. The evolution and diversity of maize has been shaped by both cultural and natural selection in a wide range of environmental and social contexts. For example, evidence has shown that cultural preferences were a factor in maintaining otherwise undesirable traits^{13, 14}. Maize seeds were selected and retained by farmers for use in subsequent growing years to maintain or enhance specific morphological traits, such as color and texture, based on changing cultural preferences, sometimes at the cost of reduced cultivability and/or nutrition¹⁵.

Politics and culture have had a strong influence on crop planting and food preferences in ancient (and modern) civilizations. Indeed, in the Andes, previous research showed that under the Inca empire, maize was a sacred crop¹⁶, symbolic of political power, and tightly associated with high-status or elite social groups^{6, 17}. It was also an important element in numerous cultural events and rituals¹⁸. From an agriculture-economic point-of-view, planting a crop on the most suitable land takes priority if a farmer wishes to maximize the output-to-input ratio. However, in certain civilizations, including Andean groups, maize was not a crop for which yield was considered a deciding factor. For example, in Andean highland communities, maize was not the most essential food crop during pre-Inca times¹⁹, not only because maize was ill-adapted to arid highland conditions and its cultivation was risky during the short growing season, but also because planting maize did not yield a return on investment. In the Andes during the Spanish conquest, the chronicler Cobo observed that, “Indians dropped all that they were doing in order to produce a little maize, even though it cost them more than it was worth”. It seemed to him that work in the fields was “…one of the major forms of recreation and festivals…”, and he was impressed by the “extraordinary fondness for farming”^{20, 21}. This did not, however, deter Andean farmers from planting an abundance of maize. In fact, while it was a desirable food, it was more so a symbol of empire¹⁷ and a ritual-feast crop, since its first introduction to the region around 3000 BC²².

Given the availability of archaeological maize samples, as well as the association of maize with cultural rituals of an historically powerful civilization and its antecedents, we were curious as to how – if at all – maize diversity was associated with Inca cultural influence and modern maize biogeographic patterns in South America. The current study focuses on maize samples from an offering associated with an elite stone tomb burial context from south of La Paz, Plurinational State of Bolivia, and dated ca. 1470 cal AD. We sequenced two archaeological Bolivian Maize (aBM) kernels with six libraries, using short read sequencing. In parallel, we used ¹⁴C Accelerator Mass Spectrometer (AMS) radiocarbon dating to determine the approximate age of the samples. Phylogenetic analysis of aBM, 16 archaeological and 226 modern maize samples (Zea mays ssp. mays L.), 88 modern samples of Teosinte (Zea mays ssp. parviglumis), 79 modern samples of Mexicana teosinte (Zea mays ssp. mexicana) and 1 Tripsacum genotype revealed genetic relationships among themselves (Table S1; Dataset S1). Finally, the unique single nucleotide polymorphisms (SNPs) in aBM infer that the selection imposed to enhance reproductive success in marginal highland growing environments resulted in increased maize diversity. Based on the data presented herein, we conclude that a combination of Inca culture, the social and political needs of local ethnic groups, and mitigation of certain natural environmental factors may have contributed to the specific preference for, and diversification of maize in at least this part of South America. As such, our conclusions offer insight into the possible pathways that influenced maize genetic diversity, thus providing a vital resource to investigate further the factors that drove selection for maize domestication and diversity at the cusp of the cultural influence of the Inca empire.

Results

Archaeological Bolivian maize sample (aBM)

La Paz is a major city located in the highlands of west-central modern Bolivia, which is 16°30′S, 68°09′W, 3,650 masl. Dry conditions in this region and the adjacent Andes Mountain range have aided in the preservation of desiccated organic material. Our maize samples were found within a woven camelid fiber pouch that was one of many offerings associated with the mummified remains of a child reported as being found in a chullpa (stone tower tomb) south of La Paz, Plurinational State of Bolivia around 1890 (Fig. 1). Although it is impossible to know the exact location of the chullpa or the origin of the individual, the assemblage was acquired by the U.S. Consul to Chile and was sent as a gift to the Michigan State University Museum, where she was curated until 2019. Following a series of recent analyses, she has since been voluntarily repatriated to Bolivia²³ Given the purported archaeological context of the mummy and her abundant offerings recorded in the MSU Museum intake records, it has been proposed that she might represent an elite member of society. However, the nature of the coarse fabrics, simple design motifs, and restricted color palette of her textiles would suggest a more common individual, perhaps with an affiliation with a local high-status group. The vast majority of the funerary towers that exist in a radius of 60 km south of La Paz city are made of mudbricks. To our knowledge, the only exception is an Inca “cushioned” (i.e., rounded edges) stone and mudbrick tower located in Uypaca, in the nearby Achocalla Valley, 8 km southwest of the city center²⁴. Because of its characteristics, this tomb most likely belonged to a local (possibly non-Inca) high-status group and was part of a more extensive funerary site. Although it is impossible to know the exact provenience of the mummy, it would not be unreasonable to suppose that she originated from this or another similar tomb in the Achocalla Valley. Regardless, in the decades following her death, community members periodically visited the chullpa and left offerings. These included the pouch full of well-preserved maize from which we collected and analyzed samples from two kernels. A single Accelerator Mass Spectrometer (AMS) date was obtained from one kernel, yielding a radiocarbon age of 400±26 BP (D-AMS-027148; maize kernel; corrected for fractionation, no δ¹³C reported). For the date 400±26 BP the two possible calibrated age ranges are 1440-1520 cal AD (p=0.823) and 1588-1621 cal AD (p=0.177), with a median age of cal AD 1474. (Calibrated at 2σ with Calib 8.2 [IntCal 20]). Two associated maize kernels from the same pouch were subjected to molecular analysis.

Fig. 1.

Genomic base modification patterns associated with ancient DNA (aDNA), including cytosine deamination damage patterns, have been identified. Observed position-specific substitutions at the ends of the sequence reads (e.g., C→T/G→A^25-29 (Fig. S1a), as well as an increase in adenine (A) and guanine (G) residues at the 5’ end, concomitant with a reduction of C and T residues^29-31 (Fig. S1b). Both features are indicative of postmortem deamination patterns. From this, we evaluated all substitution types, error rate (Fig. S1c) and read length (Fig. S2), including read length range³² and high error rate types, which are associated with aDNA genome data³³. In total, these data are consistent with the radiocarbon results, supporting that 1) with the maize is in fact consistent with an archaeological biological sample, and 2) with based on the geographic origin of the sample, the kernels potentially descended from the period of the last Inca empire before the Spaniard Francisco Pizarro’s conquest of the Inca capital, Cuzco, in A.D. 1533³⁴.

Our study focused on identifying the location from which aBM originally came, establishing and explaining patterns of genetic variability of maize, with a specific focus on maize strains that are related to our current aBM samples. We were especially interested in exploring how these patterns of variability found in both archaeological and modern samples reflect the influence of cultural systems of behavior on both the spread of variants across the ancient landscape and selection for phenotypic traits by local populations. In particular, it appears that the Inca empire had a large effect on this process of genetic differentiation.

Genetic Identity and Genetic Relatedness

After Tupac Inca Yupanqui became the emperor of the Inca Empire in the 15^th century, the territory of the Inca empire was greatly expanded across much of western South America, including the Qollasuyu, a broad territory located in what is now western Bolivia, southern Perú, northern Chile, and non-highland Argentina. Different ethnic groups populated this territory forming polities that were confederated at the time of Inca arrival. Several Aymara Pakaje groups inhabited the highlands and the inter-Andean valleys south of La Paz, from where the maize samples derive. Because maize was more of a sacred than subsistence crop of the Inca and, to a degree symbolic of Inca imperial control, the extent of its distribution as a crop may have been related to the extension of social and political hegemony over newly acquired territories. Given the provenience of the aBM and its age, it is very possible it was introduced into western highland Bolivia from the Inca core area – modern Peru. To test this hypothesis, we first evaluated the genetic relatedness of archaeological and contemporary modern maize samples using a genome sequence approach to assist in identifying the origin(s) of aBM and the relatedness to other maize from different locations. Our admixture f3-statistics test results suggest aBM is not admixed. Still, f3 outgroup with Tripsacum suggests that aBM is closely related to two archaeological lowland Peruvian maize (Z65 and Z61) among all other ancient maize and modern maize (Fig. 2). These results indicate that the closest relative of aBM is the archaeological maize that came from Peru. In addition, to predict the geographic location of aBM, we used spatial auto-correlations in genetic data, which allowed us to compare it to a set of archaeological maize samples of known geographic origin³⁵. Our results indicate the potential geographic location that aBM may have originally come from (locator.py –seed 549657840) and the predicted location is inside of Peru and very close to archaeological maize Z61 (Fig. 3; R2(x)=1.0, R2(y)=1.0). This suggests that aBM might not only be genetically related to the archaeological maize from ancient Peru, but also in the possible geographic location. This is interesting, as historically, evidence from recordings of capacocha rituals showed that they were able to source their food from the stations connecting high-altitude ritual sites and other regions of the empire³⁴. Indeed, data generated herein show that our aBM has a strong relationship with archaeological Peruvian maize in terms of both genetic identity and geographic location³⁶. This finding is also consistent with the above previous research revealing that one of the maize groups from Peru during Inca peak times spread into newly conquered territories, then slowly diverged to fit the local environment³⁷ following Inca cultural influence²⁰. After this, diverged groups were introduced outside of Inca conquered territories.

Fig. 2.

Fig. 3.

To investigate the genomic data from the aBM in relation to other archaeological and modern maize, principal component analysis (PCA) was used. As expected, the results showed that all the archaeological maize clustered with our aBM in the same cluster, which supports the authenticity of aBM as archaeological maize (Fig. 4). SM3, 5, and 10 were excluded due to their limited data to meet requirements for the PCA analysis. Unexpectedly, there is one modern Brazil maize in the ancient maize cluster that has overlap with archaeological Brazil maize Z6, which indicates that certain genetic traits from ancient Brazil maize might have been preserved in modern Brazil maize populations. There is also one modern tropicalis maize very close to archaeological Peruvian maize Z65. Considering the geographic origin of Z65 and its lowland characters, it is not surprising that modern tropicalis maize shares some level of genome sequence similarity. Following the establishment of an important maize culture center in Peru, maize was domesticated and diverged into different groups, which then spread along with Inca culture³⁸. We found that the Central American cluster has overlap with both North and South America, which supports the distribution route of maize^39-41. This finding also is consistent with previous research revealing that the original maize came from teosinte in Mexico, then dispersed northward into North America and southward into South America^{6, 7, 9, 42, 43}. We also found that a few types of maize from Europe have a much closer distance to the archaeological maize cluster compared to other modern maize, which indicates maize from Europe might expectedly share certain traits or evolutionary characteristics with ancient maize. It is also consistent with the historical fact that maize spread to Europe from multiple areas and time periods after Christopher Columbus’s late 15th century voyages to the Americas¹². But as shown, maize also has diversity inside the European maize cluster. It is possible that European farmers and merchants may have favored different phenotypic traits, and the subsequent spread of specific varieties followed the new global geopolitical maps of the Colonial era^{44, 45}.

Fig. 4.

Taken together, these results support our hypothesis that aBM was introduced into western highland Bolivia from the Inca core area – modern Peru, and aBM may possibly reveal a history of adaptation to tropical conditions following cultural influence of the Inca empire^{43, 46}. After this, maize, both as a primary food and as a symbol of power, was dispersed from Peru throughout the Inca empire, in the end enhancing maize diversity in South America.

Selection for aBM under Inca Culture Influence

After establishing the identity of aBM and its relatedness to other regional samples, we explored how aBM was shaped by human behavior, and specifically how cultural selection within the Inca empire influenced maize diversity in South America. To address this question, we examined the evolutionary relationship between aBM, its wild ancestor, and extant landraces. We constructed a Neighbour-Joining (NJ) and maximum-likelihood (ML) phylogenetic trees using the common genome-wide polymorphisms shared among aBM and modern maize samples⁴⁷. Our data included 226 samples of modern maize, 88 modern samples of Teosinte (Zea mays ssp. parviglumis), 79 modern samples of Mexicana teosinte (Zea mays ssp. mexicana), and 1 Tripsacum dactyloides as an outgroup, totaling 28,588 genome-wide SNPs. As shown in Fig. S3a and Fig. S3b, both phylogenetic trees indicate that aBM lies outside the diversity of modern maize yet clusters within a specific clade alongside one modern tropicalis maize. aBM tends to have long branches compare to tropicalis maize, which can be explained by adaption for specific local environment by time. This suggests that aBM shares a closer genetic relationship with modern tropicalis maize while being genetically distinct from other modern maize. Considering modern tropicalis maize is adapted to tropical climates exhibiting traits such as drought tolerance, heat resistance, and the ability to thrive in low-fertility soils, it is plausible that aBM might share these stress response characteristics with tropicalis maize^{48, 49}.

To investigate the selection preference on aBM under Inca cultural influence, we compared the allele frequency between aBM and other ancient maize samples. We identified 18,668 SNPs unique to aBM. Notably, two genes - Zm00001eb123120 and Zm00001eb178600-associated with traits “internode length below ear” and “nodes above ear” traits⁵⁰contain several unique SNPs in aBM (Dataset S2). This suggested that height related traits in aBM may have been selected under Inca cultural practices. Additionally, variations in the 5 or 3 untranslated regions UTR) can significantly affect gene expression. We found 83 target genes harboring unique aBM SNPs in their UTR regions (Dataset S3). Gene Ontology (GO) results showed that 3 of these genes are involved in glycosyltransferase activity, while serval of others is involved in defense response to biotic and abiotic stress and starch biosynthetic processes (Dataset S4). Combining all this information, we infer that trait associated with plant development and stress response processes might have been selected in maize under Inca cultural influence, although we could not determine whether the selection was positive or negative.

Finally, to understand positive selection in modern maize compared to ancient maize, we investigated the selective sweep using cross-population extended haplotype homozygosity (XP-EHH) method. Following the same method of previous research⁵¹, the cutoff metric for significant signatures of selection was set at a value of 6 for 5 kb windows across the genome. We only identified one target gene (Fig. 5): Zm00001eb105060, whose biological process functions include plant ovule development, chromosome segregation and asparaginyl-tRNA aminoacylation. Ovule development plays a crucial role in improving maize breeding strategies via higher seed production and hybrid vigor⁵². The positive selection of this gene in modern maize suggests that characteristics aiming for higher seed production were selected during maize diversification after aBM. This finding may relate to breeding efforts to improve seed yield in modern maize compared to ancient varieties.

Fig. 5.

In conclusion, the aBM included in the funerary accoutrements of a young elite Andean girl belongs to the family of ancient Peruvian maize groups, and likely was introduced to the Aymara region following Inca expansion in the late 15^th century. Following its dispersal as a crop across the entire Inca empire, maize was further modified by regional populations. Some genomic changes reflect selection to mitigate local environmental conditions, including response to stress. Thus, patterns of regional diversity in pre-Contact South America maize reflect complex biocultural processes. These data may also inform subsequent studies about Colonial dispersals and the origin of global maize varieties¹².

Discussion

Previous archaeological and anthropological studies have demonstrated that Inca cultural and political systems shaped access to and preferences for food. In our study we explored the evidence for the contribution of Inca culture to maize diversity in South America by analyzing 16 archaeological maize samples spanning at least 5,000 years of evolution, and 226 modern maize samples. Recognizing their shared central Andean origins, we propose that the archaeological Bolivian maize (aBM) sample came from ancient Peruvian archaeological maize. After aBM was cultivated in the Andean highland, traits adapting to climate change and environmental variability may have been selected by Inca culture. For modern maize, traits related to improved seed production through ovule development were selected compared to ancient maize. While our findings are consistent with both genetic and archaeological evidence, they raise several questions. Notably, what specific traits were selected in aBM, and how this selection occurred under Inca cultural practices is unclear, especially given our limited aBM sample size and the lack of the same genotype in archaeological maize under no Inca culture influence.

Early evidence has shown that cultural preferences were a factor in maintaining certain aesthetically desirable traits^{13, 14}, such as color and texture, even at the cost of reduced cultivability or nutrition¹⁵ by farmers. This suggests it is possible that culture played an important role in maize evolution and diversity. To better understand how exactly Inca culture influenced maize diversity, we explored the possible origin of aBM and its closest genetic relationships to ancient maize. We identified that aBM may have originated from ancient Peru, and archaeological maize from Peru has the closest genetic relationship with aBM, indicating that aBM is a product of Incan cultural dissemination. Following the expansion of the Inca civilization, maize was introduced in new territories to demonstrate the power of the Inca empire²⁰. This might lead to the selection of specific traits in maize due to cultural preferences, although this may have conflicted with natural selection pressures.

We further investigated which traits were possibly selected under Incan culture by analyzing unique SNPs in aBM. The results show that traits related to plant development, and stress response were selected. These traits can aid in resisting drought and frost in the relatively harsh conditions of the Andean highlands. The selected traits could facilitate fermentation for the creation of chicha to meet the consumption requirement of elite lifeways¹⁷. However, we are uncertain if these traits were directly selected under Inca culture due to the limited aBM sample size, which prevented us from identifying positive selection directly between aBM and other archaeological maize groups.

Climate change is expected to increasingly impact maize seed production in agricultural areas worldwide. Traditional breeding systems have increased tolerance to biotic and abiotic stress and contributed to yield in crops^{53, 54}. In this study, we compared modern maize to archaeological maize across the entire Western Hemisphere and found that the ovule development traits were selected for breeding to improving maize seed production during its evolution.

Our results not only highlight the influence of Inca culture on maize diversity in South America but also point to the adaptive characteristics of maize in improving maize seed production during evolution. This provides a potential source for maize diversity studies and breeding programs. Most importantly, this work suggests that cultural influence plays a very important role in aiding maize diversity, alongside the traits selected by early farmers.

Materials and Methods

Plant material

The two archaeological Bolivian maize (aBM) kernels were collected prior to 1890 AD from the area south of modern La Paz, Bolivia, and later curated at the Michigan State University Museum until their repatriation in 2022⁵⁵. AMS assay dates provided a range from 401 to 582 cal. BP, with a median age of 1474 cal AD. In addition, data from 16 archaeological samples collected from previous published research were also used, including two from Chile dating from 750 to 1,100 BP, two from Brazil dating from 510 to 690 BP, three from Peru dating from 600 to 1,000 BP, two from Argentina, dating from 980 to 1040 BP and from 70 to 130 BP, three from Honduras dating from 1,740 to 2,300 BP, and three from Mexico dating to 5,310 BP (Table S1)^{36, 56}.

aDNA extraction

Archaeological maize kernels were processed in a clean laboratory environment at Michigan State University. Whole kernels were crushed in a mortar and pestle, to which 600 uL of fresh PTB (N-phenacylthiazolium bromide) buffer (0.4 mg/mL proteinase K, T2.5 mM PTB, 50 nM dithiothreitol) was added. The sample was then added to 1260 uL of lysis buffer. An additional 600 uL of PTB buffer was used to rinse the mortar and added to the lysis buffer. The individually ground maize kernels were rotated overnight at 37 °C. The remaining DNA extraction followed Swarts et al., 2017⁴¹.

Library construction, amplification, and sequence

DNA extracts were converted to Illumina-compatible paired-end DNA libraries in a clean setting using SMARTer ThurPLEX DNA-seq from Takara Bio and amplified with 7 polymerase chain reactions (PCR) cycles. The generated libraries were cleaned to remove remaining adaptors with AxyPrep Mag beads (Fisher Scientific) and analyzed on a 2100 Bioanalyzer (Agilent Technologies) to check for quality. The libraries were indexed to facilitate pooling on the NovaSeq SP at the University of Illinois (Urbana-Champaign).

Geographic location map

The geographical coordinates of each sample were obtained from previously published studies, as indicated (Table S1). Utilizing latitude and longitude information, the samples were plotted on a geographic map, resulting in the successful mapping of 17 ancient samples. Geographic data visualization was carried out in R using the ggplot2, sf, leaflet, rnaturalearth, rnaturalearthdata, ggspatial, and rnaturalearthhires packages.

Genome analysis

For the aBM sequence samples data with six libraries, we combined them into a single composite sample, referred to as aBMComb. This was achieved by mapping each individual aBM sequence sample to the B73 V5 soft-masked maize genome and subsequently removing duplicate reads before merging them. We merged the sequence sample to obtain higher coverage for the analysis was needed for this study.

Subsequently, we used MapDamage v.2.2.1⁵⁷ verified cytosine deamination profiles consistent with ancient DNA, and estimated the misincorporation frequencies. Based on these results we hard-masked 5′ thymine and 3′ adenine residues within 5nt of the two ends, where deaminations was most concentrated.

Data for 16 ancient samples with specific geographic location information (maize landraces) obtained from previously published work^{36, 56, 58} by downloading sequencing from Sequence Read Archive (SRA) (https://www.ncbi.nlm.nih.gov/sra/)⁵⁹. To ensure consistency and accuracy, we downloaded the raw sequence (SRA format) for each sample and converted it into FASTQ format using fastq-dump.

For the archaeological maize samples, adapters were removed and paired reads were merged using AdapterRemoval⁶⁰ with parameters --minquality 20 --minlength 30. All 5′ thymine and 3′ adenine residues within 5nt of the two ends were hard-masked, where deamination was most concentrated. Reads were then mapped to soft-masked B73 v5 reference genome using BWA v0.7.17 with disabled seed (−l 1024 -o 0 -E 3) and a quality control threshold (−q 20) based on the recommended parameter⁶¹ to improve ancient DNA mapping. Duplicated reads were removed by GATK tools Picard (https://gatk.broadinstitute.org/hc/en-us). The coverage of each nucleotide position was estimated by bam-readcount (https://github.com/genome/bam-readcount) (Table S2), and the overall coverage of the assembled genome sequence was measured using Qualimap v2.2.1⁶².

Variant Calling

We used SequenceTools pileupCaller (v 1.5.4.0; https://github.com/stschiff/sequenceTools) to generate pseudohaploid calls for 17 archaeological maize samples. The commands executed were: samtools mpileup -R -B - q20 -Q20 & pileupCaller -- sampleNames <ancient sample name> randomHaploid -- singleStrandMode. This process yielded 16,815,782 single nucleotide polymorphisms (SNPs).

We future filtered the SNPs using estimated information from bam-readcount. SNPs supported by at least two reads were retained for each individual archaeological maize sample. We then combined the filtered SNPs from each sample, keeping only those SNPs that were present in at least one sample. This resulted in a final dataset of 2,808,917 SNPs for the 17 archaeological maize samples.

For all subsequent analyses, the ancient genotype data was merged with the genotype dataset for 394 samples from Grzybowski et al, (“Feb 03, 2023 version files” in Dryad repository: https://doi.org/10.5061/dryad.bnzs7h4f1). The merged dataset included 2,635,047 SNPs from chromosomes 1-10 that were common in both archaeological and modern maize samples, which were used for downstream analysis.

f3 statistic and phylogenetic analyses

We created a genotype file by combining 17 archaeological maize and publicly available SNP dataset form Grzybowski et al, resulting in a total of 17-archaeological and 226-modern maize, 88-modern samples of Teosinte (Zea mays ssp. parviglumis), 79-modern samples of Mexicana teosinte (Zea mays ssp. mexicana), and 1-Tripsacum genotypes.

Next, we used plink v1.9⁶³ to prune for linkage disequilibrium with the parameters --indep-pairwise 50 5 0.1, following previous research⁴⁵, to prepare the dataset for Principal Component Analysis (PCA). The final results were plotted in R using the ggplot2 package⁶⁴.

We calculated the f3 statistic^{65, 66} using the admixr package⁶⁷, which is part of the package of ADMIXTURE software suite⁶⁸, to detect evidence of admixture involving aBM. The outgroup f3 statistics was used to measure how closely aBM and each test group are related to each other. The final results were plotted using Python.

A total of 28,588 SNPs present in both aBM and all modern maize, out of a total of 2,626,087 SNPs, was used to build the phylogenetic tree. A maximum likelihood (ML) phylogenetic tree was constructed using SNPhylo⁶⁹ with parameter: -v snps.vcf -r -o tripsacum -P tree -b 1000; A Neighbour-Joining (NJ) phylogenetic tree was constructed using VCF2PopTree⁷⁰ with the “number of differences” model. Both ML and NJ phylogenies trees were inferred and visualized using FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).

Geographical location prediction

To predict the possible geographical location of aBM, we applied a deep neural network method³⁵ using only 17 ancient maize samples, as all the archaeological maize clustered in the same cluster. A total of 2,800,190 SNPs were used to predict the possible geographical location of aBM. The command used was: locator.py --vcf SNP_file --sample_data sample_data.txt –out name.

Infer selection on aBM

Due to the lack of archaeological maize with the same genotype as aBM grows outside the influence of Inca culture, it was impossible to directly identify the selection that occurred on aBM under Inca cultural influence. To determine if any of our SNPs are related to known traits, we utilized existing information on SNP-trait associations⁵⁰. To achieve potential unique selection on aBM, we calculated the allele frequency for each SNPs between aBM and other ancient maize, resulting in allele frequency data for 49,896 SNPs. Of these,18,668 SNPs were unique to aBM. We consulted SNPVersity 2.0 from MaizeGDB^{47, 71} and found that 84 of these unique SNPs were classified as “3_prime_UTR_variant” or “5_prime_UTR_variant”. These SNPs corresponded to 83 target genes, for which we retrieved Gone Ontology (GO) terms.

Identification of selection sweeps

To detect genomic signatures of recent selective sweeps in modern maize compared to ancient maize groups, SNPs that present in at least two archaeological samples and modern samples were used. We applied the cross-population extended haplotype homozygosity (XP-EHH)⁷² method. A commonly used threshold -log10(p) = 6 was set as the cutoff threshold for identifying the signal of significant selection. The regions with -log10(p) >= 6 were designated as candidate sweeps for candidate genes located within 5 KB of these regions identified, following the approach used in previous research⁵¹. MaizeMine from MaizeGDB⁷³ was used to obtain putative gene function information.

Data availability

The authors declare that all data supporting the findings of this study are included in the manuscript and its supplementary files or are available from the corresponding author upon request. Source data are provided with this paper. The genome sequence data of ancient Bolivian Maize has been deposited in the NCBI database under BioProject accession PRJNA886637. Custom scripts and single nucleotide polymorphism (SNP) calls are available on Dryad (DOI: https://doi.org/10.5061/dryad.w6m905qtd).

Acknowledgements

We are grateful to all participants who were involved in this project. We are thankful for the support of the Michigan State University Cloud Fellowship for computational skills training. This work was supported in part through computational resources and services provided by the Institute for Cyber-Enabled Research at Michigan State University. We thank Jennifer Wai, University of Chicago, for her help with DNA extraction and library preparation. We also thank Emily Josephs, Michigan State University, and Logan Kistler, Smithsonian Institution of National Museum of Natural History, for their help with the DNA analysis. Thanks to José Capriles Flores for his insightful reading and comments. The MSU Museum provided maize samples and images for documentation and analysis prior to repatriation to Bolivia. Research in the laboratory of Brad Day was supported by the MSU Research Foundation.

Supplementary figures

Fig. S1.

Fig. S2.

Fig. S4.

Additional files

Fig. S3. Phylogenetic trees of archaeological Bolivian maize (aBM). a Neighbour-Joining tree with VCF2Pop. b Maximum likelihood (ML) tree with SNPhylo. Parv (Zea mays ssp. parviglumis), Zmex (Zea mays ssp. mexicana), Europe (modern Europe), SweetC (model sweet corn), Ztrol (modern Tropical maize), mMexico (modern Mexico), mParaguay (modern Paraguay), mBrazil (modern Brazil), mPeru (modern Peru), mEcuador (modern Ecuador), mGuatemala (modern Guatemala), and mColombia (modern Colombia).

Supplementary tables. Table S1. Geographical information of the archaeological maize samples used in this study. Table S2. Paleogenomic characterization of archaeological Bolivian maize sequence samples with six libraries. Table S3. The percentage of genomic sites covered at variable depths in the archaeological Bolivian maize (aBM) sample (aBMComb_rm5nt.rmDR.bam).

Source data