The Origins of Coca: Museum Genomics Reveals Multiple Independent Domestications from Progenitor Erythroxylum gracilipes

Museum Genomics: Target Capture, Assembly, and Gene Alignment

Hybridization-based exon-capture is an efficient method for sequencing select loci from genomic DNA isolated from herbarium specimens (Hart et al. 2016; Villaverde et al. 2018). We generated 270 GB of sequence data and deposited reads for all 154 samples in the NCBI Short Read Archive under BioProject PRJNA485502. Samples were collected between 1900 and 2016 with a median collection date of 1981 (Supplementary Fig. S1 and Table S2 available on Dryad). These samples presented a spectrum of DNA quantity and degradation, but these two factors had little correlation with the success of our exon capture and target assembly. A linear model fit between the genomic DNA maximum fragment size and the number of genes recovered was nearly flat (); the polynomial regression model fit the relationship with an R of 0.024 and was also mostly flat (Supplementary Fig. S1c available on Dryad). Samples with low input DNA amounts (less than 200 ng) were scattered throughout the range of observations of gene recovery. Plants collected in the moist tropics require extra care to dry, and DNA quality thus appears to depend more on initial drying and preservation quality rather than age, corroborating (Forrest et al., 2019) (Supplementary Fig. S1a, b available on Dryad).

Not surprisingly, we observed a positive correlation between the log number of reads and the number of genes with good target coverage, defined as 75% of target length (linear , polynomial  Supplementary Fig. S1d available on Dryad). This demonstrates that we recovered 90% of exons at 75% of their length with 4 coverage per site at about 10 sequencing depth, and this depth provides a good target for sequencing effort. However, our number of reads per sample is quite variable, spanning three orders of magnitude, thus highlighting the variability incurred by combining museum samples of varying qualities together in large sequencing pools. Dilution of RNA baits and smaller numbers of samples per hybridization pool prior to exon-capture will help minimize amplification bias and overall variability of sequence reads per sample.

We removed 10 of our 154 samples from the analysis because they failed to assemble at least 90 genes with good coverage (Supplementary Fig. S2 available on Dryad). Across the 424 genes, the average assembled contig length per sample ranged from 819 bp to 8279 bp with a median of 2160 bp of sequence per sample and the 424 genes had an average of 141 samples per gene. The length of concatenated exon targets was 740,646 bp (median 1746 bp). Our assembled supercontigs (exon+intron) had an average concatenated length of 1,484,570 bp, revealing that we sequenced about 740 kb of intron and flanking sequence. Gene alignments ranged from 1749–57,834 bp in length (median 7290) with 30–93% missing data (median 59%). Cleaned trimAL alignments ranged from 129 to 8700 bp in length (median 1677) with 0.05–51.73% missing data (median 9.5%).

While previous studies have included key historical collections within broader biodiversity investigations (e.g., McCormack et al. 2016; Prosser et al. 2016; Greiman et al. 2018), our study has shown modern phylogeographic analyses can rely almost entirely on historical museum collections to generate deep and useful genomic data sets for phylogeographic research. Historical museum collections can also alleviate geographic sampling shortfalls, which can lead to erroneous inferences or oversplitting geographic populations (e.g., Jackson et al. 2017; Linck et al. 2019). This applies to systems where field collection is challenging or even dangerous but also in most biological systems due to the ever-increasing difficulty to obtain permits (Prathapan et al. 2018).

Genetic Structure of Coca and Wild Relatives

For a crop of such cultural and economic significance, our scientific understanding of the origin and evolution of coca is minimal. This is not only due to the obstacles to research imposed by its international prohibition but also the botanically challenging nature of Erythroxylum, a morphologically complex and speciose clade of tropical shrubs and treelets. Genetic sampling from the world’s largest Erythroxylum collection at the Field Museum has reshaped our fundamental understanding of the systematics of this clade (White et al. 2019).

The maximum likelihood gene trees inferred from the full gene alignments, cleaned trimAL alignments, or 212 alignments with the least missing data produced nearly identical ASTRAL summary trees, the only differences were among tips within the main clades and do not affect our presented results. Of these, the ASTRAL summary tree (AST; Fig. 2a; Supplementary Fig. S3 available on Dryad) inferred from the full gene alignments had the best local posterior probability along backbone nodes. To measure support for AST versus ML-SNP topology (topological differences are described in the Supplementary Text available on Dryad), we calculated SNP-based concordance factors (sCF; Supplementary Figs. S4 and S5 available on Dryad) and gene-tree concordance factors (gCF; Supplementary Figs. S6 and S7 available on Dryad) and found AST to be the best-supported phylogeny across our sequence and SNP data sets.

Open in a separate window

Figure 2.

Phylogeny and genetic clustering of coca varieties and wild relatives. a) ASTRAL-III summary (AST) tree topology with branch lengths optimized using the ML-SNP data set; scale indicates substitutions per site and branch support values denote concordance factors of SNPs; collapsed clades are marked with the number of samples; tip names indicate taxon, sample ID, and geographic provenance indicated by the two letter country code and political subdivision. Probability assignment of individuals to five or nine genetic clusters are indicated in bar plots along with the population assignment of ‘pure’ individuals for species trees and population genetic statistics. b) SVDquartets species tree inferred from SNP data with bootstrap support. c) ASTRAL-III species tree inferred from 424 gene trees with local posterior probability (LPP) and percent quartet support, respectively; branch lengths scaled by coalescent units. d) Constrained ASTRAL-III species trees indicating LPP and percent quartet support for an E. cataractarum  Trujillo/Colombian clade or a Huánuco gracilipes1 clade.

The AST phylogeny reveals the 37 E. gracilipes samples form a series of nested clades within which arise separate monophyletic lineages comprising E. cataractarum, E. novogranatense (Trujillo and Colombian cocas), Amazonian coca, and Huánuco coca (Fig. 2a). The goodness-of-fit estimators for the number of genetic clusters indicated the most significant information content in 5 and 9 clusters (Supplementary Fig. S8 available on Dryad; clustering results described in Supplementary Text available on Dryad).

The earliest diverging lineages of E. gracilipes are geographically distributed around the Amazon Basin and form a single cluster at and three clusters at (gracilipes2, gracilipes3, gracilipes4; Fig. 2a). Next, the E. novogranatense lineage is divided into two clades and two clusters () that define the Trujillo and Colombian coca varieties. The exceptions are samples trux.854 and trux.855 from the small, disjunct population on Río Chical at the Colombia/Ecuador border, revealing this population has a genealogical history closer to Colombian coca. This population could represent an early Colombian landrace derived from Trujillo coca, or it could be an admixed variety with combined Colombian and Trujillo ancestry. Additional sampling from Ecuador and Venezuela should improve our understanding of the history of this clade.

The remaining E. gracilipes samples form a single cluster, gracilipes1, and belong to well-supported clades with geographic structure (Fig. 2a). Within gracilipes1, 15 samples of Amazonian coca form a clade with perfect posterior probability (Supplementary Fig. S3 available on Dryad), most closely related to gracilipes1 samples from Ecuador and northern Peru. Finally, three gracilipes1 samples from the Ucayali and Madre de Dios Departments of Peru form a nested series leading to the well-supported Huánuco coca clade (Fig. 2a; Supplementary Fig. S3 available on Dryad). The fourth sample in this series, “cf.coca.884,” is described as cultivated on its herbarium voucher, but its exceptional height (3 m), leaf morphology (apex and venation), and clustering results indicate admixture with E. gracilipes (Image in Supplementary material available on Dryad).

The E. cataractarum clade diverges as sister to gracilipes1 and comprises a distinct genetic cluster (Fig. 2a). However, our ML-SNP (Supplementary Figs. S5 and S7 available on Dryad) and SVDquartets tree (Fig. 2b) stirred a prior suspicion that E. cataractarum could actually be sister to the E. novogranatense lineage (White et al. 2019). We tested support among our gene trees for this alternative placement and found it received lower quartet support (35% vs. 41%) and posterior probability (0.01 vs. 0.99). To further evaluate the evolutionary history generating this discordance, we used Treemix to model genomic admixture and gene flow across a tree representing the nine populations. This inferred gene flow from Colombian coca into E. cataractarum and, secondarily, from Amazonian coca into E. cataractarum. The third most significant migration edge was from E. cataractarum into gracilipes4 (Supplementary Fig. S8 available on Dryad). These grow in the vicinity of E. cataractarum (Figs. 1 and ​and2a),2a), and thus it is sensible that seeds and/or pollen from nearby Amazonian and Colombian coca farms have naturally introduced alleles into E. cataractarum. This result explains the phylogenomic discordance of the E. cataractarum clade and also indicates that it may not have a role in the evolution of coca, though the reverse appears to be true.

Multiple Origins Hypothesis

Our phylogenomic and clustering results indicate E. gracilipes is a paraphyletic taxon with respect to the Amazonian, Huánuco, and Colombian/Trujillo coca lineages (Fig. 2). This structure elucidates a novel hypothesis of multiple origins of domestication of coca from progenitor E. gracilipes (i.e. gracilipes1 and possibly gracilipes2), refuting both Plowman’s linear-series hypothesis beginning with Huánuco coca (Plowman 1979b) and the sister-species hypotheses suggested by Johnson and colleagues (Johnson et al. 2005; Emche et al. 2011).

While there is a clear separation of the E. novogranatense lineage (Colombian and Trujillo) from the E. coca lineage (Amazonian and Huánuco) and thus at least two domestication events, the separation of the Amazonian and Huánuco varieties, which would suggest three origins of domestication, is more tenuous. The first evidence of independent origins of Amazonian and Huánuco coca comes from the phylogeny and clustering results (Fig. 2), where they are separated by samples of gracilipes1 from central and southern Peru and a clade of Ecuadorean gracilipes1. While the placement of this Ecuadorean clade is not statistically robust, three gracilipes1 samples from Peru are firmly supported as basal to the Huánuco clade (Fig. 2a; Supplementary Figs. S2–S7 available on Dryad). The second basis for three domestications is that Amazonian coca clusters with gracilipes1 at and forms its own cluster in , suggesting closer proximity to gracilipes1 than Huánuco coca.

To further evaluate this separation and the progenitor-derivative relationships under our three-domestication hypothesis, we conducted a principal components analysis and calculated population genetic statistics across individuals from the nine populations. The PCA separates individuals into populations in agreement with the clustering analysis (Fig. 3). The nine populations, including Huánuco and Amazonian, cluster independently of one another in our bivariate plots describing the first eight components (explaining 65.9% of SNP variance; Fig. 3); further supporting the evolutionary separation of these two varieties.

Open in a separate window

Figure 3.

Genetic variation of coca varieties and wild relatives. The first eight principal components are presented in four bivariate plots.

Wright’s index of fixation () showed all pairs showed moderate to high levels of differentiation (range 0.25–0.78; Table 2). By our measures of genetic distance, Amazonian and Huánuco cocas are more similar to E. gracilipes than they are to each other, and Amazonian and gracilipes1 are the least differentiated taxa across all populations (Table 2). Amazonian and Huánuco show a similar degree of differentiation as Colombian and Trujillo (0.55 vs. 0.52, respectively). Statistics for the number of private alleles (PA; not weighted by sample size), average allelic richness per SNP (AR; weighted by sample size), observed heterozygosity (Ho), and genetic diversity (Hs) are presented in Table 1. Amazonian has the highest genetic diversity among the cocas (AR, Hs). Together, these statistics all corroborate the separation of Huánuco and Amazonian coca and thus a three-origin hypothesis of coca domestication. Corroborating expectations of genetic bottlenecks during domestication events (Gross and Olsen 2010), the coca varieties have lower genetic diversity than the wild taxa (PA, AR, Hs; see Supplementary Material available on Dryad for discussion of Ho).

Testing Alternative Domestication Scenarios

We used coalescent simulations under an approximate Bayesian computation (ABC) framework to estimate statistical support for Plowman’s linear series (Scenario 1; Plowman 1979b), Johnson’s sister species (Scenario 2; Johnson et al. 2005), a two-origin hypothesis combining Huánuco and Amazonian coca (Scenario 3), and a three-origin hypothesis (Scenario 4; Fig. 4a). Simulations generated under Scenario 3 and evaluated using the logistic regression approach received the highest posterior probability of matching the observed summary statistics from our first SNP data set, but Scenario 4 was more similar to our second SNP data set (Fig. 4b). Using the direct proportions of simulated data sets closest to our observed data sets, Scenario 3 best represented both SNP data sets (Supplementary Figs. S10 and S11 available on Dryad). The ABC results consistently showed Plowman’s linear-series and Johnson’s sister-species hypotheses are unlikely historical scenarios explaining the current genetic diversity and divergence of these taxa (Fig. 4b; Supplementary Figs. S10 and S11 available on Dryad).

Open in a separate window

Figure 4.

ABC tests and effective population size through time. a) Historical scenarios representing four domestication hypotheses. Scenarios 1–4 are drawn within colored boxes and the legend at top identifies the branch color for each taxon. b) Posterior probability of historical scenarios based on multinomial logistic regression procedure. c) Stairway plots reconstructing changes in effective population size through time for each of the four coca varieties, gracilipes1, and E. cataractarum. Solid lines are mean population size across 200 bootstraps and the 95% confidence interval is shown in gray.

Based on several analyses, we can be confident that ABC could discriminate these scenarios and estimate their probabilities. Principal components analyses show the posterior distributions for Scenarios 3 and 4 were able to accurately replicate the observed data (Supplementary Figs. S10c and S11c available on Dryad). The posterior predictive error rates specific to confidence in scenario choice ranged from 0.146 to 0.206 among direct and logistic estimates from the two data sets (Supplementary Table S5 available on Dryad). The global error computed over the whole prior distribution (parameters and scenarios) ranged from 0.243 to 0.307 (Supplementary Table S5 available on Dryad). We summarized calculations of scenario-based prior error rates in confusion matrices and estimates of type I and type II error (Supplementary Tables S6–S9 available on Dryad). We found the simulated scenarios received the highest probability in all cases, but the Type I error, the probability a correct scenario was not chosen, was high for Scenario 1 (0.39–0.48) and Scenario 2 (0.33–0.42; Scenario 3 0.17–0.23; Scenario 4 0.08–0.18). However, the Type I error is significantly lower if we only analyze confusions among Scenario 3, Scenario 4, and Scenario 1 or Scenario 2 (range 0.07–0.20), which are the meaningful comparisons with respect to our main results. Type II error, the probability a scenario was chosen when data were generated under an alternate scenario, ranged from 0.07 to 0.12 for Scenario 1, 0.15–0.16 for Scenario 2, 0.09–0.11 for Scenario 3, and 0.03–0.04 for Scenario 4. Model check results are presented in Supplementary Figs. S10 and S11 available on Dryad and prior and posterior distributions of the demographic parameters are presented in Supplementary Figs. S12 and S13 available on Dryad.

We chose not to evaluate support for two or three origins from additional SNP data sets because the multiple domestication hypothesis is a high-confidence conclusion given our genetic and geographic sampling. Our posterior samples converged on parameter estimations and posterior support for the domestication scenarios agrees with the other analyses. In addition, there are limitations in the ABC approach, the largest being that many demographic phenomena beyond our topological scenarios have influenced the allele frequencies and summary statistics our results are based upon (Sunnåker et al. 2013). For instance, the remarkably high observed heterozygosity of Amazonian coca could be explained by the accumulation of somatic mutations during the mostly clonal propagation of this crop, as is observed in other mostly asexual populations (Table 1; Stoeckel and Masson 2014). More sampling of E. gracilipes, delimitation of gracilipes1–4 as distinct taxa (see below), and distinct demographic models of each domestication event will more accurately bear on this history.

Timing and Locations of Domestication

We estimated the relative age of each coca variety by explicitly modeling the change in population sizes through time using site frequency spectra derived from our SNP data (Fig. 4c; Supplementary Fig. S14 available on Dryad). We did not estimate generation times or extrapolate mutation rates from other studies and therefore we cannot estimate the absolute timing of domestication events.

The stairway plots of gracilipes1 and E. cataractarum show these taxa have experienced a gradual population size increase followed by a very recent and drastic decrease (Fig. 4c). All coca varieties except Trujillo are marked by a stable historical population size that is smaller than gracilipes1 through much of its history. Colombian and Amazonian have about the same size, but Huánuco is smaller. Colombian coca shows the earliest departure from this stability and experienced a population decrease followed by a more recent increase, consistent with a domestication bottleneck. For Trujillo coca, population sizes have increased and then leveled off. Huánuco and Amazonian coca depart from the background population size more recently and at about the same time, although confidence intervals show changes in Amazonian could be more recent. Amazonian shows a pattern of recent population decrease whereas Huánuco coca appears to have experienced a bottleneck—a small decrease in population size followed by a large increase, then leveling off (Fig. 4c). In accordance with expectations based on historical coca farming practices, Huánuco coca is inferred to have the largest effective population size at present and Amazonian has the smallest. If we remove singletons, the model reconstructs generally similar population size histories, but without the recent population decreases in gracilipes1 and E. cataractarum, and without a rebound after the bottleneck for Colombian coca (Supplementary Fig. S14 available on Dryad). The relative coalescence times from our ABC parameter estimations also show the Trujillo/Colombian domestication was oldest and, in the case of three origins, that Amazonian coca coalesced with E. gracilipes the fewest generations ago (Supplementary Figs. S12b and S13b available on Dryad).

In addition to the Colombian coca stairway plot, the long branch subtending the E. novogranatense clade (Colombian and Trujillo; Fig. 2a) and the statistics (Table 2) also corroborate the hypothesis that this is the oldest coca crop. Archaeological leaf fragments and chewing paraphernalia reveal coca culture was well established 8000 years BP (Plowman 1984; Dillehay et al. 2010). Evidence from the archaeological record, mating systems, and hybrid crosses (see Bohm et al. 1982) led Plowman to believe that Colombian coca was derived from Trujillo. Under neutral evolutionary models, we would expect higher genetic diversity in progenitor taxa (Gross and Olsen 2010; Feng et al. 2020), but we see Trujillo has fewer private alleles and lower allelic richness, genetic diversity, and observed heterozygosity than Colombian coca (Table 1). However, these statistics are also influenced by different demographic histories (see Mortimer 1901; Plowman 1984; Gootenberg 2008), so we maintain the working hypothesis established by Bohm et al. (1982) that Trujillo coca or a common ancestor was the progenitor of Colombian coca. Thus, a likely region of this domestication event was in Ecuador or northern Peru, near the current range of cultivation.

Our phylogenetic and clustering analyses indicate Huánuco coca was domesticated in the eastern Andean foothills of southern Peru. The two gracilipes1 samples from Madre de Dios, Peru (grac.298, grac.654) cluster with and are placed at the base of the Huánuco coca clade (Fig. 2a). Neither of these samples was collected in proximity (100 km) to coca farms, but this does not rule out the possibility of postdomestication gene flow from coca farms in this region into gracilipes1. A remarkable sample (cf.coca.884) from central Peru is described as cultivated but is morphologically intermediate and admixed between with gracilipes1 (Image in Supplementary material available on Dryad). A morphologically wild-type E. gracilipes from the montane forest in Ecuador (grac.100) was also inferred to be slightly admixed with Huánuco coca, but this is likely due to postdomestication gene flow because the Ecuadorean gracilipes1 forms a clade separate from Huánuco coca (Fig. 2a).

The Amazonian clade is most closely related to gracilipes1 individuals from Loreto, Peru, and Amazonian Ecuador (Fig. 2; Supplementary Figs. S5 and S7 available on Dryad). At , Amazonian coca is not distinct from gracilipes1, but it does emerge as a separate genetic cluster when nine clusters are defined (Fig. 2a). Lastly, the clustering and results, as well as the fact that Amazonian coca has the highest genetic diversity of any coca variety (Table 1), support the hypothesis that Amazonian is the most recently domesticated coca crop. This is in line with linguistic and ethnographic observations that indigenous agricultural practices, preparation, consumption, and terminology surrounding coca are remarkably consistent across the Amazon basin, and thus thought to have been more recently dispersed throughout its current range (Plowman 1986). Thus, Amazonian coca was either recently derived from Huánuco coca or it was independently domesticated from gracilipes1 in Amazonian Ecuador or northern Peru.

We thank the following herbaria for use of material: F, HUEFS, MOL. Thanks to E. Gardner and M. Johnson for guidance with library preparation and Hyb-Seq analyses, J. Walsh for assistance with DNA isolation, and the following people for help in the field and lab: A. Daza, M. Huinga, J. Janovec, F. Parra, E. Vosburgh, and J. Wells. Thanks to Proyecto Khoka in Medellín for collaboration. Specimen collection in Colombia was performed under the general collecting permit Resoluciones ANLA 1177 9 octubre 2014 and 1386 29 octubre 2015 awarded to Universidad de los Andes.