{'fraction_non_alphanumeric': 0.033617569474265384, 'fraction_numerical': 0.04435310717133918, 'mean_word_length': 4.803488786044856, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Archaeological Central American maize genomes suggest ancient gene flow from South America

Significance Maize is a global food staple with great economic and cultural importance. Archaeogenomic studies have revealed a process of protracted maize domestication and multiple waves of human-mediated dispersal in the Americas. Maize first arrived in South America as a partial domesticate, where the domestication syndrome became independently fixed and improved varieties developed away from the influence of wild gene flow. We demonstrate that hybrids of some of these improved varieties were likely reintroduced back to Central America. We hypothesize that this backflow of South American genetic material may have contributed to the development of a more productive staple, which was related to the growth and aggregation of human populations, and the formation of more complex social and political structures regionally.


Morphometric Permutation Analysis
We performed 100,000 random permutations to estimate the probability of the observed 3 difference in each morphometric variable (Fig. S1). Each iteration of the test randomly 4 reassigned values with replacement to maximize the number of independent 5 simulations. The absolute difference between the simulated groups were calculated and 6 each value was compiled to generate a sampling distribution. Permutation p-values 7 were calculated using the proportion of the 100,000 simulations that were larger than 8 the observed difference. All statistics and permutation test were performed using R (1). 9 10 Radiocarbon Dating 11 12 Each maize cob was subsampled for AMS 14 C dating at in the Human Palaeoecology 13 and Isotope Geochemistry Laboratory at the The Pennsylvania State University. Each 14 Samples were then pretreated with repeated baths in1MHCl and NaOH at 70 °C for 15 30min on a heater block. A final acid wash removed secondary carbonates formed 16 during the base treatment. Samples were then returned to neutral pH with two 15-min 17 baths in deionized water at 70 °C to remove chlorides, and dried on a heater block. 18 Sample CO2 was produced by combustion at 900 °C for 3 h in evacuated sealed quartz 19 tubes using a CuO oxygen source and Ag wire to remove chloride compounds. Primary 20 (OX-1) and secondary (FIRI-D/F, FIRI-H) standards and a Queets Wood background 21 were selected to match the sample age and underwent the same chemical steps for 22 quality assurance. Samples were graphitized at the Keck Carbon Cycle Accelerator 23 Mass Spectrometer facility at the University of California, Irvine and AMS 14 C 24 measurements were made using a modified NEC 1.5SDH-1 instrument; National 25 Electrostatics Corporation. All 14 C ages were δ 13 C-corrected for mass-dependent 26 fractionation with measured 13 C/ 12 C values (2) and calibrated using OxCal version 4.3 27 (3) using the IntCal13 northern hemisphere curve (4). 28 29 DNA extraction and genomic data collection 30 26 maize samples were prepared in the dedicated ancient DNA clean lab facilities at the 31 Penn State Anthropology Department (n=23), and the Smithsonian Institution's Museum 32 Support Center (n=3). Standard protocols to prevent and detect contamination were 33 utilized (5), including strict workflow procedures, frequent cleaning with bleach and 34 ethanol, use of complete personal protective equipment, and the preparation and 35 sequencing of negative control reactions. Both labs are equipped with filtered air 36 handling and are regularly decontaminated. 37 We extracted DNA, prepared sequencing libraries, and screened samples following 38 established protocols for highly degraded ancient DNA (Supplemental Methods), 39 identifying EG84, EG85, and EG90 as suitable for genomic sequencing. 40 DNA was extracted from archaeological maize cob (n=24), stem (n=1), and leaf (n=1) 41 tissue exactly following the protocol described by Wales and Kistler (6). At Penn State, 42 we prepared DNA sequencing libraries exactly as described in Kistler et al (7), and at 43 the Smithsonian we employed the Blunt End Single Tube procedure (8) with 44 modifications described in (9), dual indexing with primers and primer sequences 1 described in (10), and Platinum Taq High-Fidelity (Invitrogen) for library amplification. 2 Samples were pooled in roughly equimolar ratios, and screened on a NextSeq 550 3 High-output flow cell with 75bp single-end reads, and a HiSeq X10 lane with 150bp 4 paired-end reads. Samples with sufficient endogenous DNA for genome-scale analyses 5 were sequenced completely on a HiSeq X10. All sequencing was carried out at Admera 6 Health, South Plainfield, NJ. 7 Sample reads were adapter-trimmed and paired reads were merged using 8 AdapterRemoval 2 (11), and mapped to the maize reference genome (Zea mays B73 9 RefGen_v4; (12)) using the Burrowes-Wheeler Aligner aln function (13) with seed 10 disabled to improve ancient DNA mapping (11) and a minimum mapping quality of 20.

11
We used mapDamage 2.0 (14) to verify cytosine deamination profiles consistent with 12 authentically ancient DNA, and all 5' thymine and 3' adenine residues were hard-13 masked within 5nt of sequence ends where deamination was most concentrated. All 14 analyses were restricted to the strictly mappable fraction of the maize genome, as 15 previously described (15); mappability mask previously published in (15). 16 Using the set of 17,672,809 SNPs described in (15), we generated pseudohaplotype 17 SNP calls at all sites with a minimum 2x consensus, exactly as described previously 18 (15). Using this approach we recovered 1,786,417, 3,312,860, and 2,243,175 SNPs 19 from EG84, EG85, and EG90 respectively. In addition, we followed the approach of (16) 20 and combined the three samples for most analyses, treating them as a single population 21 sample. We merged the alignment files for the three samples, and re-called the SNP 22 panel as a single set of pseudohaplotypes in this case, yielding 7,666,836 SNPs for 23 analysis. We combined new SNP calls from El Gigante with the previously reported set 24 of SNP calls available at (15), consisting of 109 modern genomes and 11 ancient 25 genomes before culling for missingness during analyses. For analyses assuming SNPs 26 in linkage equilibrium (e.g. model-based clustering), we pruned for linkage using the 27 plink "--indep-pairwise" function. The complete SNP dataset including separate and 28 combined El Gigante maize is available on Dryad.

30
Genomic analyses 31 32 Model-based clustering 33 34 We used ADMIXTURE (17) to estimate ancestry proportions under model-based 35 clustering in maize genomes, excluding teosintes and the partially domesticated mid-36 Holocene genomes from Mexico (18,19). We included all genomes with at least 25% of 37 sites called, enforced a minimum 50% of samples called to retain a site, minimum minor 38 allele frequency of 0.02, and using the LD-pruned SNP set. This analysis included 39 4,252,422 SNPs and 98 genomes, using the combined El Gigante dataset and setting 40 k=5 as previously established (7). We ran 100 independent analyses, and compared the 41 results by final log-likelihood (lnL). lnL values were bimodal. The majority of runs, 72%, 42 clustered tightly together with a higher lnL range, with the remaining 28% represented a 43 more diffuse lower tier. Among the better supported upper tier, the El Gigante genome 44 contains an estimated 97.73%-98.04% Pan-American ancestry, and 1.95-2.26% South 45 American ancestry. No other ancestry cluster contributes significantly (max 0.0016%) in 1 any run. The small proportion of South American ancestry is attributed to the 2 Andean/Pacific group. Given genetic and historical ties between this and the Lowland 3 lineage (7), we interpret this as a generic signature for an ancestral South American 4 gene pool. 5 6 f-statistics 7 8 We used previously released scripts for f4 and outgroup-f3 calculations ((15) 9 plink2freq.pl, f3.pl, f4.pl), and included the complete, unpruned dataset at all sites where 10 the outgroup Tripsacum dactyloides was present. We used a block jackknife resampling 11 procedure with 5Mb blocks to estimate standard error and calculate a Z-score to assess 12 fit to the null hypothesis. Statistical significance for rejecting the null hypothesis was 13 concluded where |Z| > 3. 14 15 Ancestry informative marker (AIM) domestication analysis discovery 16 17 We used a perl script (Dryad: AIM.pl) to calculate ancestry informativeness (In; (20)  18 between pairs of populations at all SNPs called in at least half the samples from each 19 group. Following previous research into maize domestication status (18), we designated 20 all SNPs with In ≥ 0.1 as ancestry informative markers. For domestication analysis, we 21 compared 1) all modern domesticated maize with ≤ 20% ancestry in the highland 22 Central American cluster with 2) all parviglumis and mexicana genomes. Highland 23 Central American maize carries previously documented admixture from highland 24 mexicana teosinte (21), and thus all sites are not reliably maize-like for In determination. 25 We assessed a panel of previously identified genes associated with domestication (22)  26 containing at least 10 AIMs in the region containing the gene, plus 10kbp upstream and 27 downstream, yielding 199 total genes. For each genome, we calculated the proportion 28 of teosinte-like alleles in each domestication gene region with at least 5 called alleles at 29 AIMs, as described above. The El Gigante set of proportions could then be compared 30 against the set of comparable values for maize and teosinte as described in the maize 31 text.

33
To assess the domestication status of individual genes, we used a likelihood-based 34 method to test whether alleles at AIMs for a given gene were more likely drawn from a 35 population resembling modern maize or modern teosinte. We considered all AIMs in 36 each domestication gene regions as above, and computed a gene's likelihood of 37 originating from a reference population (maize or teosinte) on the basis of the sample 38 allele's frequency in the reference population. Log-likelihood was therefore calculated 39 as the sum of the natural log of the frequency of the test sample's alleles in a reference 40 population: 41 42 ln 43 Where the test sample's allele was not present in the reference population, we set the 44 allele frequency at a nominal 0.01 to preempt the log of 0, assuming that the test allele 45 could easily be unsampled or lost to drift. We then computed Bayes Factors (BF) 1 following (23) as the ratio of lnL values for maize and teosinte. We concluded 2 "substantial" evidence for one of the competing hypothesis when BF ≥ 3 or BF ≤ 1/3, 3 and "strong" evidence when BF ≥ 10 or BF ≤ 1/10, following (23).

5
For visualization (Fig. 2), we normalized the log-likelihood ratios on a scale of -1 to 1 as: 6 7 8 9 10 AIM analysis for South American affinity 11 12 For South American affinity AIM analysis, we compared 1) modern domesticated maize 13 with ≥95% combined Andean/Pacific and Lowland South American ancestry with 2) 14 modern domesticated maize having ≤5% combined Andean/Pacific and Lowland South 15 American ancestry, and computed In to identify AIMs as above. We then assessed 16 affinity to South American ancestry by computing the proportion of these geographic 17 AIMs carrying South American alleles in each individual sample. 18 19 20 Admixture graph fitting 21 22 We included all samples with Pan-American or South American ancestry (≥99% on the 23 basis of model-based clustering), plus El Gigante maize, all parviglumis, and Tripsacum 24 dactyloides for an outgroup. We divided the Pan-American lineage into northern and 25 southern geographic sets across the Isthmus of Panama, yielding 6 total populations. 26 We first used AdmixTools (24) to calculate all permutations of f4-statistics as input for 27 AdmixtureGraph (25), which we used to explore the permutation space of graphs. We 28 first exhaustively enumerated all 3885 possible graphs relating these 6 populations with 29 up to one admixture event, and fitted each of these to the f4-statistics. Each graph was 30 fit five times, retaining the best scoring fit (as evaluated using the "best_error" score). 31 None of the graphs without any admixture events provided good fits. Among those with 32 one admixture event, two graphs provided decent fits to the data, each with four minor 33 outlier f4-statistics (after these, the next best graph had nine outliers). The first of these 34 was the topology (Tripsacum,(parviglumis,(South America,(El Gigante,(Pan-Am 35 North,Pan-Am South))))) as shown in Figure 3a, with an admixture event from the 36 ancestor of the South America lineage into Pan-Am South lineage. The second had the 37 same structure, but with the admixture instead from the Pan-Am South lineage into the 38 South American lineage. These two topologies are equivalent with respect to the f4-39 statistics and thus achieved identical fits.

41
We then took the first of these graphs, studied the four outlier statistics that it did not 42 perfectly predict, and hypothesized a second admixture event from the parviglumis 43 lineage into the lineage ancestral to El Gigante, Pan-Am North and Pan-Am South. This 44 improved the fit and left only one minor outlier statistic (|Z|=3.3), though the inferred 45 admixture proportions were not stable across repeated fits. We then refitted this graph 1 using qpGraph (24), which uses f2 and f3-statistics in addition to f4-statistics, obtaining 2 stable admixture proportions and achieving a good fit without any outlier f-statistics 3 (largest |Z| = 2.7).

5
Genome Size Estimation 6 7 Heterochromatic knob content is the primary determinant of genome size differences in 8 maize, and the proportion of sequence reads mappable to heterochromatic knobs has 9 been demonstrated to reflect genome size estimates from flow cytometry (26,27). We 10 therefore used the proportion of reads mapping to the 180bp knob fraction of the maize 11 genome as a proxy for genome size, following (26). We used a custom mapping 12 strategy modeled after (28) to independently map sequence reads using a method 13 capable of assigning reads to highly repetitive transposable element and 14 heterochromatic knob fractions of the genome. We first generated a unique 15 transposable element (UTE) and heterochromatic knob (knobC) reference set exactly 16 as described in (28,29) (reference fasta files curated on Dryad), and obtained the 17 maize filtered gene set version ZmB73_5b_FGS_genes.fasta from (ftp.gramene.org). 18 We created SMALT (https://www.sanger.ac.uk/tool/smalt-0/) indexes from these three 19 reference targets with a step size of 3 and a word length of 12. We then used SMALT to 20 first attempt mapping to the knobC set, then passed remaining unmapped reads to the 21 UTE. Elements of the FGS also occur in the UTE, and therefore reads were treated as 22 transposable rather than genic in origin if they could be assigned to the UTE, regardless 23 of gene occupancy. Finally, only reads failing to map to both repetitive databases were 24 handed down to the FGS for genic read alignment.

26
The proportion of all mapped reads assigned to the 180bp knob elements of the knobC 27 fraction was summarized in terms of RPKM-reads per kilobase (following (29) and 28 used for genome size analysis. In this case, RPKM180bp = R/(K × M × 10 −6 ), where R is 29 the number of reads mapped to 180bp knob elements in the knobC database, K is the 30 combined length of the 180bp knob elements in the knobC database, and M is the total 31 number of reads mapped to the combined knobC, UTE, and FGS genomic fractions-32 the complete mappable set of reads. Because of sequence-based and genomic biases 33 in ancient DNA degradation (30), including specifically in maize (18), we did not attempt 34 genome size estimation in archaeological maize.  The absolute difference between the simulated groups were calculated and each value was compiled to 6 generate a sampling distribution. Permutation p-values were calculated using the proportion of the 7 100,000 simulations that were larger than the observed difference. Three maize samples with dates after 8 1900 BP were exclude in the tests.