Recovery of chloroplast genomes from medieval millet grains excavated from the Areni-1 cave in southern Armenia

Panicum miliaceum L. was domesticated in northern China at least 7000 years ago and was subsequentially adopted in many areas throughout Eurasia. One such locale is Areni-1 an archaeological cave site in Southern Armenia, where vast quantities archaeobotanical material were well preserved via desiccation. The rich botanical material found at Areni-1 includes P. miliaceum grains that were identified morphologically and14C dated to the medieval period (873 ± 36 CE and 1118 ± 35 CE). To investigate the demographic and evolutionary history of the Areni-1 millet, we used ancient DNA extraction, hybridization capture enrichment, and high throughput sequencing to assemble three chloroplast genomes from the medieval grains and then compared these sequences to 50 modern P. miliaceum chloroplast genomes. Overall, the chloroplast genomes contained a low amount of diversity with domesticated accessions separated by a maximum of 5 SNPs and little inference on demography could be made. However, in phylogenies the chloroplast genomes separated into two clades, similar to what has been reported for nuclear DNA from P. miliaceum. The chloroplast genomes of two wild (undomesticated) accessions of P. miliaceum contained a relatively large number of variants, 11 SNPs, not found in the domesticated accessions. These results demonstrate that P. miliaceum grains from archaeological sites can preserve DNA for at least 1000 years and serve as a genetic resource to study the domestication of this cereal crop.


Results
Areni-1 millet. The Areni-1 cave located in southern Armenia produced a vast array botanical material from occupations dating to 4000 BCE ( Fig. 1 and Table 1). Exceptionally well preserved, intact Panicum miliaceum L. florets were identified from both Trench 1 and 3 based on their morphology (Table 2 and Fig. S1; methods described more fully below). Florets are the basic threshing unit of P. miliaceum, exhibiting an indurate, smooth, and shiny palea and lemma tightly enclosing an inner caryopsis, with a pointed apex and relatively blunt proximal end (Fig. 2). The vast majority of ancient millet remains across southwest Asia tend to have been preserved via charring, which frequently damages the palea and lemma, leaving only the caryopsis/seed 24 ; the desiccated remains from Areni-1 differ in that the colour and husk (palea and lemma) were maintained (Fig. 2). The lack of charring undoubtedly enhanced the preservation of DNA. The P. miliaceum florets from Areni-1 measure roughly 3.3 × 2.2 mm. These dimensions are slightly larger than charred caryopses reported from Mongolia, Afghanistan, Turkey, Iran, and elsewhere [25][26][27] , but this not surprising given that the medieval remains from Areni-1 represent a long history of selection post-dating earlier domestication events and were preserved via desiccation with the palea and lemma of the floret intact.
To confirm the contextual dating of the Areni-1 P. miliaceum, 14 C dating was performed on botanical material associated with the grains. The size of the Arerni-1 P. miliaceum did not provide sufficient mass to perform 14 C dating and extraction of aDNA from a single seed, so 14 C dating needed to be conducted indirectly using associated botanical material. The 14 C dating determined the Areni-1 grains were ≈ 1000-years-old and belonged to the Early and High medieval periods (Tables 1 and 2).
Mapping Areni-1 shotgun data to chloroplast reference genome. In the shotgun libraries, 11292a and 11295a were the only two grains that yielded a significant level of endogenous DNA and produced a large number of unique reads mapping to the chloroplast reference (> 90,000). The remaining shotgun libraries produced more than a thousand-fold fewer unique chloroplast mapped reads ( Table 3). The extraction blanks and certain grains (11293a and 11294b) contained a comparatively small number of unique mapped reads in both the shotgun and cpDNA enriched libraries, suggesting that contamination from laboratory sources or between libraries during the processing of the Areni-1 samples was not a substantial issue. In the shotgun libraries, grains 11293a and 11294b both produced fewer mapped reads than the extraction blanks, but these results are likely very stochastic as these libraries contained a very small fraction of mappable reads. www.nature.com/scientificreports/ Mapping of Areni-1 enriched cpDNA data to chloroplast reference genome. All ancient shotgun libraries were taken through two rounds of cpDNA hybridization capture regardless of endogenous DNA content in order to treat all samples in a similar manner and maximize the recovery of target sequences. Of the Areni-1 P. miliaceum enriched for cpDNA, only three grains (11292a, 11294a, and 11295a) produced libraries that contained a sufficient number of reads that mapped to the chloroplast reference genome (> 19,000) that allowed reconstruction of the Areni-1 chloroplast genomes (Table 3). No reads from the ancient libraries could be mapped to the IRs of the chloroplast genome, because the short-fragmented sequences in this data could not be accurately placed in these nearly identical loci (Fig. 3). As the IRs do not contain any mapped reads, these loci were removed from all further analysis. Libraries 11292a and 11295a covered > 99.7% of the reference (sans IRs) with at least one read, whilst the fraction of the reference covered by at least one read in library 11294a was approximately 8% less than the other two libraries.

Read length distribution and mapDamage analysis of enriched cpDNA libraries.
To support the ≈1000-year-old 14 C dating of the Areni-1 grains, sequences from the 11292a, 11294a, and 11295a cpDNA enriched libraries were examined for characteristics of authentic aDNA: short read lengths and the presence of deamination miscoding lesions 28,29 . All libraries from the three grains contained reads that were fragmented to  www.nature.com/scientificreports/ a short length, generally < 150 bp (Fig. 4a). Grain 11294a produced the shortest reads, suggesting that this sample was the most degraded of the three. In the cpDNA enriched libraries, the read length of grains 11292a and 11295a was shifted towards longer sequences, which is expected as larger DNA fragments form more stable complexes with enrichment probes and are thus recovered more efficiently in hybridization capture procedures 23 . In contrast, the enriched library for 11294a experienced a downward shift in read length, which likely stems from the stochastic effect of having only 72 unique mapped reads in the shotgun library for this sample. In mapDamage analysis of the cpDNA enriched libraries, the '5 ends of reads exhibited the increased levels of C→T misin-  www.nature.com/scientificreports/ corporations caused by deamination of cytosine typical of aDNA ( Fig. 4b and Fig. S2) 30,31 . MapDamage results are not presented for shotgun data because the library from grain 11294a contained an insufficient number of reads for this analysis. Since direct 14 C dating of the Areni-1 broomcorn millet used for sequencing was not possible due to the small mass of the grains, the fragmentation patterns and miscoding profiles are important data that support the medieval dating of the samples.
Mapping statistics for Areni-1 enriched cpDNA data. In the mapping of cpDNA enriched libraries, grains 11292a and 11295a produced similar numbers of unique reads (≈ 170,000) that mapped to the chloroplast genome reference, while sample 11294a generated approximately a tenth of the mapped reads observed in the other two grains. The maximum read depth was greater with grains 11292a and 11295a at 229 and 331 reads respectively, compared to 52 reads for 11294a (Fig. 3). Despite the differences in coverage of the chloroplast reference genome and read depth, all three samples produced chloroplast data that allowed for phylogenetic and haplotype network analysis.

Generation of chloroplast genomes from 50 modern accessions of broomcorn millet.
For comparative purposes, chloroplast genomes were generated from 50 accessions of modern P. miliaceum (Table S1). These accessions had a mainly Eurasian distribution and were weighted towards Chinese accessions with 29 of the samples originating from various locations within China. Included in these accessions were two wild (undomesticated) broomcorn millet (SampleID 31 and 167) and two accessions that have previously been used to produce nuclear reference genomes (SampleID 69 and T296) 32,33 . The sequencing of these modern accessions generated average read depths of ≥ 1356 reads across the reference chloroplast genome.

Paniceae phylogenetic analysis.
To confirm the visual identification and the taxonomic placement of the Areni-1 P. miliaceum and validate our laboratory and analytical methods, the chloroplast genomes from these medieval grains were compared to similar sequences from 11 members of the Paniceae tribe, including a modern P. miliaceum (Table S2). A phylogeny was generated of the Paniceae and Areni-1 chloroplast genomes with maximum likelihood method (Fig. 5). In the tree, the Areni-1 samples and broomcorn millet were placed together on a branch with a bootstrap support of 100, which indicated that the medieval grains were closely related to broomcorn millet and confirmed the taxonomic placement and visual identification of the medieval grains. The placement of the Panicum and Setaria species in the tree suggested the topography produced in this www.nature.com/scientificreports/ analysis was correct. All the Panicum spp. were placed on the same branch within the tree and domesticated S. italica formed a clade with the wild antecedent of this cereal, Setaria viridis 9 . www.nature.com/scientificreports/

Panicum miliaceum phylogenetic and haplotype network analysis.
At the time the current study was started little was known about P. miliaceum chloroplast diversity or what demographic or evolutionary information the chloroplast genome from this cereal could contained. To address these deficiencies, we compared the Areni-1 chloroplast genomes to similar organelle sequences from 50 modern varieties of P. miliaceum including two wild accessions. Between the modern accessions and the Areni-1 grains a total of 63 SNPs were identified in the chloroplast genomes and these genotypes were used to generate phylogenies using maximum likelihood, Bayesian, and neighbor joining methods ( Fig. 6 and Fig. S3). Generally, across the trees there was low bootstrap support and no geographic structure, so no inference on P. miliaceum demographics could be made. However, in the trees produced with the maximum likelihood and Bayesian methodologies, the domesticated accessions were split into two large clades and in each of these phylogenies the Areni-1 chloroplast genomes grouped together as a clade. The relationship between the modern and Areni-1chloroplast genomes was further investigated in a haplotype network using the PopArt software with the median spanning network algorithm (Fig. 7). In this network, all the domesticated accessions, including the Areni-1 samples, were separated by only 1-5 SNPs. A core haplogroup, represented by the large node at the centre of the network, contained 19 accessions with a wide geographic range across Eurasia. The Areni-1 chloroplast genomes were separated from the core haplogroup by two mutations, one of which was common to all three medieval samples. The largest number of variations was observed in the wild P. miliaceum, which were separated from the domesticated accessions by [11][12][13][14] SNPs.

Discussion
In this study, DNA extraction and hybridization capture enrichment techniques were used to recover aDNA from medieval broomcorn millet grains found in Areni-1, a cave site located in the Armenian Highlands a region that has an important history in the domestication and dispersal of many crops. Despite variable quality of aDNA found in the Areni-1 grains, we were able to recover three chloroplast genomes from five grains assayed.
Despite a similar well-preserved physical appearance, the endogenous DNA content varied among the Areni-1 grains used in this study as did the quality of the aDNA. For example, in the millet enriched for cpDNA, grain 11295a contained the highest level of endogenous DNA and the longest read length, while grain 11294a contained the least amount of endogenous DNA with the shortest read length ( Table 2, Fig. 3a). The environment of the cave appears to be suited to prevent decomposition of biological material and macromolecules as excavations at Areni-1 have unearthed vertebrate soft tissues including a mummified goat 34 and brain tissue in a human skull 35 , demonstrating that the conditions in the cave can preserve ephemeral biological material. It is not clear why the Areni-1 millet exhibited different levels of DNA preservation. Grains with poor preservation may have undergone different harvesting, threshing, or other anthropogenic treatments that influenced DNA stability. Post excavation factors, such as storage or handling conditions, may have also influenced the different levels of DNA preservation. Further, the micro-environment has been shown to impact the survival of DNA in mammalian bones 36 and similar processes may have influenced the decay of aDNA in the Areni-1 millet grains. Local variations in pH, Variable preservation of aDNA from grains is not unique to the Areni-1 millet as a study of 6000 year old barley found large differences in recoverable endogenous DNA in grains excavated from the same occupation layer of a cave 38 . Regardless of the varying integrity of the aDNA, the recovery of chloroplast genomes from the Areni-1 samples demonstrates that broomcorn millet grains can, under the right conditions, preserve DNA for at least 1000 years. Even with a poorly preserved sample (11294a), hybridization capture enrichment was able to recover sufficient cpDNA to perform a chloroplast phylogenetic and haplotype analysis.
In the Paniceae grass chloroplast phylogeny (Fig. 5), the Areni-1 and P. miliaceum reference clustered with several species native to North America (Panicum virgatum, switchgrass; Panicum capillare, witchgrass) and Whiteochloa capillipes, which is from Australia, New Guinea, and Indonesia. Panicum sumatrense, a cereal crop originally from India also clustered with the Areni-1/ P. miliaceum. Panicum is the largest genera of the grass family Paniceae and has a cosmopolitan distribution so the clustering of the Areni-1/ P. miliaceum with widely distributed species is not unusual 39 . While it is unlikely that the specimens from Areni-1 are directly related to these species, it is possible that some of them share a common ancestor from Southeast Asia, further supporting the east to west introduction of P. miliaceum in Eurasia. It should be noted that the chloroplast phylogeny indicates a close relationship between broomcorn millet and P. capillare (bootstrap support = 100) an observation that supports a previous study of Panicum nuclear DNA, which found that diploid P. capillare or a close relative is the likely maternal genome donor to tetraploid broomcorn millet 40 .
To investigate the demographic history of P. miliaceum, chloroplast phylogenies were generated for the Areni-1 samples and modern accessions using maximum likelihood, Bayesian, neighbour joining methodologies. In these analyses only 63 SNPs were identified in a group that included 48 modern domesticated accessions, two wild accessions, and the three ≈1000-year-old grains from the Areni-1 cave. Accordingly, there was little differences between these chloroplast genomes which was reflected in the phylogenies with low bootstrap support values through much of the trees. No conclusions to be drawn on the demographics of broomcorn millet form these analyses ( Fig. 6 and Fig. S3). The inability of chloroplast genomes to resolve intraspecies relationships is not unique to the current study. Research on the genera Amaranthus and Euphrasia found that whole www.nature.com/scientificreports/ chloroplast genomes were highly conserved and lacked the power to discriminate between some intraspecies accessions 41,42 . Yet, in all these trees the Areni-1 chloroplast genomes fell within the diversity of domesticated accession, indicating the Areni-1 grains are a fully domesticated variety of P. miliaceum and not a wild or partially domesticated form of the cereal. However, in the maximum likelihood and Bayesian phylogenies, the domesticated accessions were split into two clades with no geographical structure and each clade containing accessions from locations throughout Eurasia ( Fig. 6 and Fig. S3a). These results are similar those found in a study by Hunt et al. 43  www.nature.com/scientificreports/ processes may be involved with the discrepancies between the P. miliaceum nuclear and chloroplast genomes, however the low levels of diversity in the chloroplast genomes are likely the main contributing factor for these differences. The P. miliaceum chloroplast genomes contained sufficient information for gross scale separation of the different accessions but was insufficient to accurately resolve fine scale relationships. In the haplotype network (Fig. 7), the majority of the haplotypes were separated by 1-5 SNPs, which did not allow for any conclusions to be drawn about the demographic history of P. miliaceum. Wild accessions contained a relatively high number of variants, 11 SNPs, that were not present in the domesticated accessions. Several scenarios that are not mutually exclusive could explain this observation. One possibility is that the wild accessions sampled may not be directly related to the broomcorn lineages that were domesticated. Another possible scenario is the SNPs observed in the wild accessions may represent diversity lost in the chloroplast genome during a domestication bottle neck. Such a loss of diversity has been reported in a study by Leigh et al. 46 , which described a reduction in both alleles and haplotypes in the chloroplast genomes of different wheat species after domestication. However, a much larger sampling of chloroplast genomes from wild accessions will be needed to resolve the relationship of undomesticated and domesticated P. miliaceum.
Domestication of crops involve the selection of plant traits that better meet the needs of humans. Chloroplasts are not only responsible for photosynthetic energy production, but are also involved with other cellular functions including synthesis of amino acids, nucleotides, fatty acids, phytohormones, vitamins, and the production of metabolites that are important for plant interaction with biotic and abiotic factors in the environment 16 . Since the chloroplast regulates such a broad range of cellular functions that are essential to P. miliaceum for survival and success as a cereal, traits under the control of this organelle will come under selection by humans. The similarity of the Areni-1 chloroplast genomes to those from modern domestications may represent a limited range of genotypes that human find beneficial to P. miliaceum.
The chloroplast inheritance pattern will aid in maintaining chloroplast genotypes that are beneficial as a crop. Although unproven in P. miliaceum, as far as the authors can determine, maternal inheritance of chloroplast in P. miliaceum is the mostly likely scenario because maternal inheritance occurs in most higher plants and maternal inheritance has been shown to be the primary inheritance pattern in Panicum virgatum, a relative of P. miliaceum 47,48 . With maternal inheritance there is no meiosis or recombination so chloroplast genotypes beneficial to humans will not be broken up and can be more readily maintained in a cultivated plant 49 .
It is possible that mapping of promiscuous DNA that has broken off the chloroplast genome and become incorporated in either the nuclear or mitochondrial genome has influenced the Areni-1 chloroplast results 50,51 . Effort was made to minimize the impact of promiscuous cpDNA on the Areni-1 results by requiring a relatively high read depth (≥ 10) and at least 95% of the reads to call a mutation as an alternative allele. However, as it will be difficult to distinguish between promiscuous cpDNA and DNA from the chloroplast, proving that promiscuous cpDNA has not influenced the current results is extremely challenging. Regardless, the agreement of the morphological identification of the Areni-1 millet as broomcorn millet and the placement of the Areni-1 millet with P. miliaceum in the Paniceae chloroplast phylogeny (Fig. 5) indicate that promiscuous cpDNA did not have a significant impact on the interspecies analysis.
In conclusion, here we generated 53 P. miliaceum chloroplast genomes, which includes three that are ≈1000-yeaers-old that can be used as resources to study this cereal crop and other C4 plants. A phylogeny generated from these chloroplast genomes divided the domesticated forms of broomcorn millet into two clades, an observation that parallels a report using nuclear DNA. The chloroplast genomes from domesticated accessions displayed low diversity and wild accessions sampled in this study contained SNPs that were outside the diversity of the domesticated forms of P. miliaceum. Lastly, this study demonstrates that broomcorn millet grains can preserve significant levels of DNA for hundreds if not thousands of years and can serve as an important genetic resource to investigate the evolution and domestication of this cereal.

Materials and methods
Areni-1 archaeology of botanical material. To identify the plant material at Areni-1, a comprehensive archaeobotanical sampling strategy was adopted for the cave sediment. Sediment samples with a volume of 5 L were collected throughout the excavation and separated using dry sieving instead of flotation. The majority of plant remains found at Areni-1 were desiccated and contact with water was likely to damage this material. After passing each sediment sample through a 1 mm sieve, material > 1 mm was kept for study and < 1 mm was discarded, which allowed for some small weed seeds and plant parts to be lost. Similar approaches to sieving have been used at other sites that contained copious amounts of desiccated plant material 52 . Sieved material was hand sorted into categories of plant remains, bones, and pottery fragments. Millet grains were among the plant material recovered from Areni-1 excavations, which also included other domesticated cereals including wheat and barley. To date, all millet found at Areni-1 came from Trench 1 and 3 within the cave site and in stratigraphic horizons attributed to medieval occupations 11 .
Areni-1 millet grain morphology and identification. Identification of the Areni-1 millet grains were made by comparing morphological features of the ancient florets with modern Panicum, Echinochloa, and Setaria specimens accessioned from Turkey and Armenia (curated within the Archaeobotany reference collection at the University of Connecticut). To help further hone the identification, a number of images, keys, and textual descriptions of morphological features and species distributions were also consulted, including reference manuals 27,[53][54][55][56] , the Flora of Armenia 57 , Flora of Turkey 58 , Flora of Iraq 59 , and archaeobotanical publications detailing various additional Paniceae species found across Asia and Europe 25,60-64 . Textual descriptions within flora typically provide size ranges for spikelets rather than florets or caryopses. These descriptions remained helpful, however, given that the spikelet lengths of many taxa (particularly within the Setaria genus), were shorter than www.nature.com/scientificreports/ the florets observed here, and could easily be eliminated. Length to breadth ratios were also considered along with the surface texture of the florets. The Areni-1 millet grains were photographed and measured using a Nikon AZ microscope with associated NIS Elements software.
Carbon dating of botanical material from Areni-1. A single Vitis sp. pedicel from Trench 1 and a pool of three millet grains from Trench 3 (Fig. S1) were sent to the NSF Arizona AMS Laboratory (University of Arizona) for radiocarbon dating. Fraction Modern Carbon and Radiocarbon Age were calculated as weighted averages of combined machine runs to reduce overall error.
Areni-1 millet library construction and cpDNA enrichment. For genetic studies, five millet grains were chosen at random for the extraction and enrichment of cpDNA ( Table 2). Extraction of aDNA and construction of Illumina libraries from the Areni-1 millet was performed in a dedicated low-DNA laboratory where no work had previously been done with millet. Isolation of aDNA from the Areni-1 grains and negative blank controls was performed with a two-step protocol, which has been previously described [65][66][67] . Shotgun and cpDNA enriched libraries were generated using standard aDNA methods and a detailed description of these protocols is given in the Supplemental Methods.
Mapping of Areni-1 millet sequencing data. Raw 73 and SAMTools and plotted with R using the ggplot2 package. For all libraries, except the shotgun library for 11294a, read length analysis was performed using 19,000 reads, a number based on the smallest number of mapped reads produced by a library in this dataset. The shotgun library for 11294a produced 72 mapped reads and this number was used for the read length analysis.
Sequencing of modern broomcorn millet accessions. DNA was extracted from seedling leaves of 50 modern broomcorn millet accessions (Table S1)  Mapping of modern broomcorn millet sequencing data. Fastq reads from the modern broomcorn millet accessions were firstly mapped to a broomcorn millet nuclear reference genome (NCBI Assembly: GCA_003046395.2) using BWA mem (v0.5.11-foss-2016b: https:// arxiv. org/ abs/ 1303. 3997v2). From the resulting bam file, reads that did not map to the nuclear reference were extracted using SAMTools (v1.12) and remapped to the broomcorn millet chloroplast genome (Accession: KU343177.1) using BWA mem. Duplicates were marked with Picard Tools to prevent processing of these reads in downstream analysis. Mapping to the chloroplast reference genome produced average read depths of between 1356 and 4548 across the chloroplast genome (Table S2).
Variant calling in the Areni-1 and modern millet sequencing data. Variants were called in parallel in the Areni-1 and modern millet using freebayes (v 1.0.2-GCC-4.9.3-binutils-2.25: Garrison and Marth, 2012, arXiv:1207.3907) with the following parameters: --min-base-quality 30 --min-coverage 10 --hwe-priors-off --ploidy 1 --pooled-continuous --genotype-qualities --min-alternate-fraction 0.95 to generate vcf files. The value for the --min-alternate-fraction was chosen to mask polymorphic sites caused by promiscuous cpDNA (chloroplast DNA that has migrated to the nuclear or mitochondrial genomes) 74 . Freebayes automatically disregards reads marked as duplicate, so only unique reads from the modern millet were analysed. Chloroplast genome consensus sequences were generated for the Areni-1 millet using the FastaAlternateReferenceMaker function of GATK (v 3.5-Java-1.8.0_71) 75 using the freebayes vcf files and the KU343177.1 reference. Two sections of the Areni-1 consensus sequences, bp 81851 to 104489 and 117101 to 139826, were removed with Geneious Prime (Build 2019-04-26) using 11295a as a guide 76 , because these regions represent the IRs of the chloroplast genome where ancient reads could not be mapped (Fig. 3).
Paniceae phylogenetic analysis. Chloroplast reference genomes for 10 members of the Paniceae tribe of grasses (which includes members of the Panicum genus) were downloaded from the NCBI website (Table S1). The reference sequence for the Zea mays chloroplast genome was also downloaded for use as an out-group. Although the Grass Phylogeny Working Group II, which used three genetic markers, placed Whiteochloa capillipes in the Cenchrinae subtribe 77 , this species was included in the analysis because a more recent study of com-  78 . The Areni-1 consensus sequences and the grass references were aligned using MAFFT (v 7.130-foss-2016b-with-extensions) and the default settings 79 . P. miliaceum variety 69 from the modern accessions (Table S2) was included the MAFFT alignment to serve as the broomcorn millet chloroplast genome. After alignment, the Paniceae genomes were trimmed of the IRs with Geneious Prime using sample 11295a as before. Poorly aligned regions were removed from the alignment using the trimAl program with the default parameters on the Phylemon 2.0 website (http:// phyle mon2. bioin fo. cipf. es) 80 . A maximum likelihood tree was built for the alignment with RAxML 81 (v. 8.2.10) with 1,000 rapid bootstraps and the "GTRGAMMA" model. Tree visualization was performed with FigTree (v 1.4.3; http:// tree. bio. ed. ac. uk/ softw are/ figtr ee/) with rooting on the branch containing the maize chloroplast genome.
Broomcorn millet phylogenetic and haplotype analysis. The translated genotypes of the Areni-1 and modern millet were extracted from the freebayes vcf file and exported as a multiple sequence file using bcftools (v1.12) 71 . An MAFFT alignment with the IRs and poorly aligned regions removed was produced as before. Phylogenies were then produced for the chloroplast genomes using three methodologies: (1) maximum likelihood using RAxML as describe above; (2) Bayesian using MrBayes (v3.2.7a, https:// nbisw eden. github. io/ MrBay es/ downl oad. html with 5 million generations 82 ; (3) neighbor joining MEGA11 (https:// megas oftwa re. net/) with 1000 bootstrap replicates 83 . Trees were generated for the phylogenies using FigTree as above. A PHYLIP format file generated from the MAFFT alignment using Geneious Prime was used to produce a haplotype network in the PopART program (http:// popart. otago. ac. nz/ index. shtml) with the median joining network algorithm 84 . Use of plant material. The collection and use of the medieval and modern millet grains complied with all relevant institutional, national, and international guidelines and legislation on the study of plant material. No permission was needed for any of the millet grains used in this study.  www.nature.com/scientificreports/