The domestication of Cucurbita argyrosperma as revealed by the genome of its wild relative

Abstract

Despite their economic importance and well-characterized domestication syndrome, the genomic impact of domestication and the identification of variants underlying the domestication traits in Cucurbita species (pumpkins and squashes) is currently lacking. Cucurbita argyrosperma, also known as cushaw pumpkin or silver-seed gourd, is a Mexican crop consumed primarily for its seeds rather than fruit flesh. This makes it a good model to study Cucurbita domestication, as seeds were an essential component of early Mesoamerican diet and likely the first targets of human-guided selection in pumpkins and squashes. We obtained population-level data using tunable Genotype by Sequencing libraries for 192 individuals of the wild and domesticated subspecies of C. argyrosperma across Mexico. We also assembled the first high-quality wild Cucurbita genome. Comparative genomic analyses revealed several structural variants and presence/absence of genes related to domestication. Our results indicate a monophyletic origin of this domesticated crop in the lowlands of Jalisco. We found evidence of gene flow between the domesticated and wild subspecies, which likely alleviated the effects of the domestication bottleneck. We uncovered candidate domestication genes that are involved in the regulation of growth hormones, plant defense mechanisms, seed development, and germination. The presence of shared selected alleles with the closely related species Cucurbita moschata suggests domestication-related introgression between both taxa.

References

[1] Meyer, R. S. & Purugganan, M. D. Evolution of crop species: genetics of domestication and diversification. Nat. Rev. Genet. 14, 840-852 (2013).

[2] Zeder, M. A. Core questions in domestication research. Proc. Natl Acad. Sci. USA 112, 3191-3198 (2015).

[3] Kistler, L. et al. Gourds and squashes (Cucurbita spp.) adapted to megafaunal extinction and ecological anachronism through domestication . Proc. Natl Acad. Sci. USA 112, 15107-15112 (2015).

[4] Paris, H. S. in Genetics and Genomics of Cucurbitaceae 111-154 (Springer International Publishing, 2016). https://doi.org/10.1007/7397_2016_3.

[5] Castellanos-Morales, G. et al. Historical biogeography and phylogeny of Cucurbita: Insights from ancestral area reconstruction and niche evolution. Mol. Phylogenet Evol. 128, 38-54 (2018).

[6] Sanjur, O. I., Piperno, D. R., Andres, T. C. & Wessel-Beaver, L. Phylogenetic relationships among domesticated and wild species of Cucurbita (Cucurbitaceae) inferred from a mitochondrial gene: implications for crop plant evolution and areas of origin. Proc. Natl Acad. Sci. USA 99, 535-540 (2002).

[7] Zizumbo-Villarreal, D., Flores-Silva, A. & Marín, P. C.-G. The Archaic Diet in mesoamerica: incentive for milpa development and species domestication. Economic Bot. 66, 328-343 (2012).

[8] Chomicki, G., Schaefer, H. & Renner, S. S. Origin and domestication of Cucurbitaceae crops: insights from phylogenies, genomics and archaeology. New Phytol. 226, 1240-1255 (2019).

[9] Whitaker, T. W. & Cutler, H. C. Cucurbits and cultures in the Americas. Economic Bot. 19, 344-349 (1965).

[10] Piperno, D. R., Ranere, A. J., Holst, I., Iriarte, J. & Dickau, R. Starch grain and phytolith evidence for early ninth millennium B.P. maize from the Central Balsas River Valley, Mexico. Proc. Natl Acad. Sci. USA 106, 5019-5024 (2009).

[11] Sánchez-de la Vega, G. et al. Genetic resources in the Calabaza Pipiana Squash (Cucurbita argyrosperma) in Mexico: genetic diversity, genetic differentiation and distribution models . Front Plant Sci. 9, 400 (2018).

[12] Montes-Hernandez, S. & Eguiarte, L. E. Genetic structure and indirect estimates of gene flow in three taxa of Cucurbita (Cucurbitaceae) in western Mexico. Am. J. Bot. 89, 1156-1163 (2002).

[13] Lira, R. et al. in Ethnobotany of Mexico 389-401 (Springer New York, 2016). https://doi.org/10.1007/978-1-4614-6669-7_15.

[14] Barrera-Redondo, J. et al. The genome of Cucurbita argyrosperma (Silver-Seed Gourd) reveals faster rates of protein-coding gene and long noncoding RNA turnover and neofunctionalization within Cucurbita . Mol. Plant 12, 506-520 (2019).

[15] Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, E. M. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210-3212 (2015).

[16] Sun, H. et al. Karyotype stability and unbiased fractionation in the Paleo-Allotetraploid Cucurbita genomes. Mol. Plant 10, 1293-1306 (2017).

[17] Hoff, K. J., Lange, S., Lomsadze, A., Borodovsky, M. & Stanke, M. BRAKER1: unsupervised RNA-Seq-based genome annotation with GeneMark-ET and AUGUSTUS. Bioinformatics 32, 767-769 (2016).

[18] Alonge, M. et al. RaGOO: fast and accurate reference-guided scaffolding of draft genomes. Genome Biol. 20, 224 (2019).

[19] Whitaker, T. W. & Bemis, W. P. Origin and evolution of the cultivated Cucurbita . Bull. Torrey Bot. Club 102, 362-368 (1975).

[20] Ott, A. et al. tGBS® genotyping-by-sequencing enables reliable genotyping of heterozygous loci. Nucleic Acids Res. 45, e178 (2017).

[21] Lee, T. H., Guo, H., Wang, X., Kim, C. & Paterson, A. H. SNPhylo: a pipeline to construct a phylogenetic tree from huge SNP data. BMC Genomics 15, 162 (2014).

[22] Alexander, D. H., Novembre, J. & Lange, K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655-1664 (2009).

[23] Excoffier, L. & Foll, M. Fastsimcoal: a continuous-time coalescent simulator of genomic diversity under arbitrarily complex evolutionary scenarios. Bioinformatics 27, 1332-1334 (2011).

[24] de Villemereuil, P. & Gaggiotti, O. E. A new F_ST-based method to uncover local adaptation using environmental variables . Methods Ecol. Evol. 6, 1248-1258 (2015).

[25] Luu, K., Bazin, E. & Blum, M. G. PCAdapt: an R package to perform genome scans for selection based on principal component analysis. Mol. Ecol. Resour. 17, 67-77 (2017).

[26] Caye, K., Jumentier, B., Lepeule, J. & François, O. LFMM 2: fast and accurate inference of gene-environment associations in genome-wide studies. Mol. Biol. Evol. 36, 852-860 (2019).

[27] Frichot, E., Mathieu, F., Trouillon, T., Bouchard, G. & François, O. Fast and efficient estimation of individual ancestry coefficients. Genetics 196, 973-983 (2014).

[28] Khan, A. W. et al. Super-Pangenome by integrating the wild side of a species for accelerated crop improvement. Trends Plant Sci. 25, 148-158 (2020).

[29] Wilson, H. D. Gene flow in squash species. BioScience 40, 449-455 (1990).

[30] Jarvis, D. I. et al. A global perspective of the richness and evenness of traditional crop-variety diversity maintained by farming communities. Proc. Natl Acad. Sci. USA 105, 5326-5331 (2008).

[31] Montes-Hernández, S., Merrick, L. C. & Eguiarte, L. E. Maintenance of squash (Cucurbita spp.) landrace diversity by farmers activities in Mexico . Genet. Resour. Crop Evol. 52, 697-707 (2005).

[32] Barrera-Redondo, J. et al. Landrace diversity and local selection criteria of domesticated squashes and gourds (Cucurbita) in the central Andean mountain range of Peru: Tomayquichua, Huánuco . Bot. Sci. 98, 101-116 (2020).

[33] Ranere, A. J., Piperno, D. R., Holst, I., Dickau, R. & Iriarte, J. The cultural and chronological context of early Holocene maize and squash domestication in the Central Balsas River Valley, Mexico. Proc. Natl Acad. Sci. USA 106, 5014-5018 (2009).

[34] Zizumbo-Villarreal, D. & Colunga-GarcíaMarín, P. Origin of agriculture and plant domestication in West Mesoamerica. Genet. Resour. Crop Evol. 57, 813-825 (2010).

[35] Stinnesbeck, W. et al. The earliest settlers of Mesoamerica date back to the late Pleistocene. PLoS ONE 12, e0183345 (2017).

[36] Piperno, D. R. The origins of plant cultivation and domestication in the new world tropics. Curr. Anthropol. 52, S453-S470 (2011).

[37] Lowry, D. B. et al. Responsible RAD: striving for best practices in population genomic studies of adaptation. Mol. Ecol. Resour. 17, 366-369 (2017).

[38] Zhang, J. et al. Receptor-like cytoplasmic kinases integrate signaling from multiple plant immune receptors and are targeted by a Pseudomonas syringae effector. Cell Host Microbe 7, 290-301 (2010).

[39] Zeilmaker, T. et al. DOWNY MILDEW RESISTANT 6 and DMR6-LIKE OXYGENASE 1 are partially redundant but distinct suppressors of immunity in Arabidopsis . Plant J. 81, 210-222 (2015).

[40] Hayashi, S. et al. The glycerophosphoryl diester phosphodiesterase-like proteins SHV3 and its homologs play important roles in cell wall organization. Plant Cell Physiol. 49, 1522-1535 (2008).

[41] de Jesús-Pires, C. et al. Plant thaumatin-like proteins: function, evolution and biotechnological applications. Curr. Protein Pept. Sci. 21, 36-51 (2020).

[42] Batanero, E., Ledesma, A., Villalba, M. & Rodríguez, R. Purification, amino acid sequence and immunological characterization of Ole e 6, a cysteine-enriched allergen from olive tree pollen. FEBS Lett. 410, 293-296 (1997).

[43] Moreira, X., Abdala-Roberts, L., Gols, R. & Francisco, M. Plant domestication decreases both constitutive and induced chemical defenses by direct selection against defensive traits. Sci. Rep. 8, 12678 (2018).

[44] Ulm, R. et al. Distinct regulation of salinity and genotoxic stress responses by Arabidopsis MAP kinase phosphatase 1. EMBO J. 21, 6483-6493 (2002).

[45] Chen, K. et al. Abscisic acid dynamics, signaling, and functions in plants. J. Integr. Plant Biol. 62, 25-54 (2020).

[46] Martínez, A. B. et al. Differences in seed dormancy associated with the domestication of Cucurbita maxima: elucidation of some mechanisms behind this response. Seed Sci. Res. 28, 1-7 (2017).

[47] Park, J. E., Kim, Y. S., Yoon, H. K. & Park, C. M. Functional characterization of a small auxin-up RNA gene in apical hook development in Arabidopsis . Plant Sci. 172, 150-157 (2007).

[48] Yang, L., Jiang, Z., Jing, Y. & Lin, R. PIF1 and RVE1 form a transcriptional feedback loop to control light-mediated seed germination in Arabidopsis . J. Integr. plant Biol. 62, 1372-1384 (2020).

[49] Hare, P. D., Møller, S. G., Huang, L. F. & Chua, N. H. LAF3, a novel factor required for normal phytochrome A signaling. Plant Physiol. 133, 1592-1604 (2003).

[50] Liscum, E. & Reed, J. W. Genetics of Aux/IAA and ARF action in plant growth and development. Plant Mol. Biol. 49, 387-400 (2002).

[51] Burla, B. et al. Vacuolar transport of abscisic acid glucosyl ester is mediated by ATP-binding cassette and proton-antiport mechanisms in Arabidopsis . Plant Physiol. 163, 1446-1458 (2013).

[52] Seiler, C. et al. ABA biosynthesis and degradation contributing to ABA homeostasis during barley seed development under control and terminal drought-stress conditions. J. Exp. Bot. 62, 2615-2632 (2011).

[53] Wang, Y. et al. Construction of a high-density genetic map and analysis of seed-related traits using specific length amplified fragment sequencing for Cucurbita maxima . Front Plant Sci. 10, 1782 (2019).

[54] Guo, D. et al. Resequencing 200 flax cultivated accessions identifies candidate genes related to seed size and weight and reveals signatures of artificial selection. Front Plant Sci. 10, 1682 (2019).

[55] Navarro, L. et al. A plant miRNA contributes to antibacterial resistance by repressing auxin signaling. Science 312, 436-439 (2006).

[56] Holdsworth, W. L., LaPlant, K. E., Bell, D. C., Jahn, M. M. & Mazourek, M. Cultivar-based introgression mapping reveals wild species-derived Pm-0, the major powdery mildew resistance locus in squash. PloS ONE 11, e0167715 (2016).

[57] Dierckxsens, N., Mardulyn, P. & Smits, G. NOVOPlasty: de novo assembly of organelle genomes from whole genome data. Nucleic Acids Res. 45, e18 (2017).

[58] Kim, D., Paggi, J. M., Park, C., Bennett, C. & Salzberg, S. L. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37, 907-915 (2019).

[59] Alverson, A. J. et al. Insights into the evolution of mitochondrial genome size from complete sequences of Citrullus lanatus and Cucurbita pepo (Cucurbitaceae) . Mol. Biol. Evol. 27, 1436-1448 (2010).

[60] Kajitani, R. et al. Efficient de novo assembly of highly heterozygous genomes from whole-genome shotgun short reads. Genome Res. 24, 1384-1395 (2014).

[61] Ye, C., Hill, C. M., Wu, S., Ruan, J. & Ma, Z. S. DBG2OLC: efficient assembly of large genomes using long erroneous reads of the third generation sequencing technologies. Sci. Rep. 6, 31900 (2016).

[62] Li, H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34, 3094-3100 (2018).

[63] Walker, B. J. et al. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS ONE 9, e112963 (2014).

[64] Li, H. & Durbin, R. Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics 26, 589-595 (2010).

[65] Humann, J. L., Lee, T., Ficklin, S. & Main, D. Structural and functional annotation of eukaryotic genomes with GenSAS. Methods Mol. Biol. 1962, 29-51 (2019).

[66] Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15-21 (2013).

[67] Stanke, M., Schöffmann, O., Morgenstern, B. & Waack, S. Gene prediction in eukaryotes with a generalized hidden Markov model that uses hints from external sources. BMC Bioinform. 7, 62 (2006).

[68] Jones, P. et al. InterProScan 5: genome-scale protein function classification. Bioinformatics 30, 1236-1240 (2014).

[69] Schneider, M. et al. The UniProtKB/Swiss-Prot knowledgebase and its Plant Proteome Annotation Program. J. Proteom. 72, 567-573 (2009).

[70] Camacho, C. et al. BLAST+: architecture and applications. BMC Bioinforma. 10, 421 (2009).

[71] Koren, S. et al. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation. Genome Res. 27, 722-736 (2017).

[72] Gong, L., Pachner, M., Kalai, K. & Lelley, T. SSR-based genetic linkage map of Cucurbita moschata and its synteny with Cucurbita pepo . Genome 51, 878-887 (2008).

[73] Haug-Baltzell, A., Stephens, S. A., Davey, S., Scheidegger, C. E. & Lyons, E. SynMap2 and SynMap3D: web-based whole-genome synteny browsers. Bioinformatics 33, 2197-2198 (2017).

[74] Goel, M., Sun, H., Jiao, W. B. & Schneeberger, K. SyRI: finding genomic rearrangements and local sequence differences from whole-genome assemblies. Genome Biol. 20, 277 (2019).

[75] Kurtz, S. et al. Versatile and open software for comparing large genomes. Genome Biol. 5, R12 (2004).

[76] Sedlazeck, F. J. et al. Accurate detection of complex structural variations using single-molecule sequencing. Nat. Methods 15, 461-468 (2018).

[77] Alexa, A., Rahnenführer, J. & Lengauer, T. Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics 22, 1600-1607 (2006).

[78] Pratas, D., Silva, R. M., Pinho, A. J. & Ferreira, P. J. An alignment-free method to find and visualize rearrangements between pairs of DNA sequences. Sci. Rep. 5, 10203 (2015).

[79] Li, S. & Chou, H. H. LUCY2: an interactive DNA sequence quality trimming and vector removal tool. Bioinformatics 20, 2865-2866 (2004).

[80] Hoffmann, S. et al. Fast mapping of short sequences with mismatches, insertions and deletions using index structures. PLoS Comput Biol. 5, e1000502 (2009).

[81] Caboche, S., Audebert, C., Lemoine, Y. & Hot, D. Comparison of mapping algorithms used in high-throughput sequencing: application to Ion Torrent data. BMC Genomics 15, 264 (2014).

[82] Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078-2079 (2009).

[83] Li, H., Ruan, J. & Durbin, R. Mapping short DNA sequencing reads and calling variants using mapping quality scores. Genome Res. 18, 1851-1858 (2008).

[84] Purcell, S. et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559-575 (2007).

[85] Wigginton, J. E., Cutler, D. J. & Abecasis, G. R. A note on exact tests of Hardy-Weinberg equilibrium. Am. J. Hum. Genet. 76, 887-893 (2005).

[86] Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, 2156-2158 (2011).

[87] Keenan, K., McGinnity, P., Cross, T. F., Crozier, W. W. & Prodöhl, P. A. diveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors. Methods Ecol. Evolution 4, 782-788 (2013).

[88] Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A. & Cresko, W. A. Stacks: an analysis tool set for population genomics. Mol. Ecol. 22, 3124-3140 (2013).

[89] Gutenkunst, R. N., Hernandez, R. D., Williamson, S. H. & Bustamante, C. D. Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. PLoS Genet. 5, e1000695 (2009).

[90] Ardelean, C. F. et al. Evidence of human occupation in Mexico around the Last Glacial Maximum. Nature 548, 87-92 (2020).

[91] Gelman, A. & Rubin, D. B. Inference from iterative simulation using multiple sequences. Stat. Sci. 7, 457-472 (1992).

[92] Malinsky, M., Matschiner, M. & Svardal, H. Dsuite-fast D-statistics and related admixture evidence from VCF files. Mol. Ecol. Resour. (2020).

[93] Cingolani, P. et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin) 6, 80-92 (2012).

[94] Zheng, Y. et al. Cucurbit Genomics Database (CuGenDB): a central portal for comparative and functional genomics of cucurbit crops. Nucleic Acids Res. 47, D1128-D1136 (2019).