Origin and domestication of papaya Y^h chromosome

Sequencing and analysis of a papaya MSY

The papaya HSY has been sequenced previously, using BAC libraries from the gynodioecious cultivar “SunUp,” a popular commercial transgenic cultivar from Hawaii (Wang et al. 2012). To study the MSY, the dioecious “AU9” cultivar from the USDA germplasm repository in Hilo was selected, because it was one of the few dioecious accessions available with a desirable fruit appearance. AU9 was collected from a papaya breeding program at Griffith University in Queensland, Australia, which was traced back to a papaya breeding program in Israel with plants of unknown origin. Because of the pericentric location and high repeat content of the Y-linked region, we sequenced the papaya MSY using a reiterative BAC-by-BAC approach. A physical map of the papaya MSY was constructed from the AU9 male BAC library, using probes from BACs in the HSY physical map (Gschwend et al. 2011). Candidate BACs were confirmed by a combination of PCR and BAC end sequencing, and gaps in the physical map were filled using chromosome walking (see Methods). The MSY physical map consists of a minimum tiling path of 99 BACs (Supplemental Fig. 1) with a combined nonoverlapping length of 8.1 Mb. The MSY BACs were sequenced using 454 sequencing technology (Roche) to an average depth of 50× for each BAC. Assembled contigs were anchored using the HSY as a reference, and some gaps were filled with Illumina whole-genome shotgun reads of AU9.

The HSY and MSY sequences were aligned and their nucleotide divergence was estimated to be 0.4%, with 32,517 sequence differences in total. The MSY and HSY sequences are collinear, except for an 8398-bp insertion in the MSY corresponding to a newly integrated Ty3-gypsy retrotransposon; this transposon insertion is the cause of the earlier overestimation of the 1.4% sequence divergence between MSY and HSY, which was based on only one pair of homologous BACs (Yu et al. 2008). There may be other small-scale differences between the HSY and MSY because of small sequence gaps in both assemblies. Most BACs are complete, but some highly repetitive BACs (>90% repeats) remain in multiple, ordered contigs. Overall, the gene content, exon structure, and gene order are conserved between the entire HSY and MSY regions (Fig. 1). A SVP-like gene has a copia-like retrotransposon insertion in the first intron of the HSY allele, resulting in a transcript with a partial sequence of the second intron (Ueno et al. 2014). Ninety-four transcription units were annotated on both Y chromosomes. Of these, 66 are complete, intact protein coding genes with annotated start and stop codons, and 28 (30%) are pseudogenes (truncated gene fragments, or sequences whose coding regions contain premature stop codons) (Supplemental Tables 1, 2). In HSY and MSY sequences, 39 (59%) of the intact genes have identical sequences in the coding region. The coding regions of seven genes (11%) include only synonymous differences, and 20 genes (30%) have only nonsynonymous differences (Supplemental Table 3). Below, we describe HSY and MSY sequences from more individuals, which exclude most of these genes with differences as candidates for the hermaphrodite-male difference.

The previously sequenced HSY and X revealed two major inversions, one of which (inversion 1) was inferred to have caused the initial recombination suppression between the X and Y chromosomes (Ming et al. 2011; Wang et al. 2012). The sequence of the AU9 MSY also has these inversions, which therefore distinguish both Y chromosomes from the X. In both the MSY and the HSY, inversion 1 includes 30 of the 66 intact Y-linked genes, including all 16 that have no X-linked counterparts, as well as 21 pseudogene sequences (Supplemental Table 1). Inversion 2 includes 18 genes that appear to correspond to complete transcripts, all of them with X-linked counterparts, plus three pseudogenes. Finally the MSY and HSY share a “collinear region,” adjacent to border B of inversion 2, which includes 18 genes, all with X and MSY (and HSY) copies, and two pseudogenes.

Repetitive sequences make up ∼79.2% of the MSY pseudomolecule, similar to the HSY, which has 79.3% repetitive element content. As with the HSY, LTR elements dominate the MSY repetitive element content, making up 42.7% of these sequences from the MSY (Supplemental Table 4). In contrast, the corresponding X region has a total repetitive element content of 67.2%, lower than either Y chromosome though higher than the papaya genome-wide average of 51.9% (Ming et al. 2008; Wang et al. 2012; Na et al. 2014).

Y chromosome sequence variation, population structure, and fixed Y–Y^h differences

The high sequence conservation between the MSY and HSY led us to assess their diversity in wild and cultivated papaya populations, as this can help identify both the origin of the Y^h chromosome and the differences found between all such chromosomes and the Y of males. We selected hermaphrodites from 12 gynodioceous commercial cultivars with diverse origins and agronomic traits (Supplemental Fig. 2) and males from 24 dioecious wild papaya accessions collected from 10 populations across Costa Rica for whole-genome resequencing (Supplemental Table 5). The wild papaya accessions were chosen based on the previous finding, based on four pairs of X/Y genes, that three distinct populations of Y chromosomes can be identified: MSY1, MSY2, and MSY3 (Weingartner and Moore 2012). We used individuals from the same populations.

After filtering and aligning the reads to the papaya SunUp draft genome and the HSY sequence (see Methods) (Ming et al. 2008), the average sequencing depth for each individual was 15.6× for autosomal sequences and 7.8× for the haploid HSY and X regions (Supplemental Table 4). In the complete data set, we identified a total of 3,301,017 polymorphisms (SNPs and small indels) genome wide and 66,579 Y-specific variants (shared by all MSY and HSY sequences and not found among the 36 X-linked ones) (see Supplemental Table 6). The vast majority of the Y polymorphisms are in intergenic regions, and only 3885 are in coding regions (3685 in introns and 180 in exons). One hundred nine of the 180 exon variants (61%) polymorphisms are nonsynonymous; 69 (38%) are synonymous; and two (1%) cause frame shifts in coding genes. Most of the Y-linked variants are private to single (or a few) individuals. For example, one of the two frame-shift variants is in CpXY26, a phospholipid-transporting ATPase that has a premature stop codon in the cultivated variety HCAR314, and the other is a premature stop codon in gene CpXY10 in the wild male Cp112; this gene is a homolog of FLOWING LOCUS T. Neither of these is a candidate for the difference between males and hermaphrodites, as they are not fixed between the sets of plants of these two sex types. We discuss the few fixed differences below.

The SNPs in Y-specific regions of the MSY were used to investigate the population structure of the Y chromosomes. Model-based analysis, maximum likelihood phylogenetic reconstruction, and principal component analysis (PCA) all confirm the previous finding of three distinct populations of Y chromosomes (Fig. 2B–D; Weingartner and Moore 2012). The MSY1 and MSY2 Y haplotypes were found in wild populations from two geographic regions on the two opposite coasts of Costa Rica (with nine and six males from six and five subpopulations, respectively), and three of these subpopulations (operational populations 8, 10, and 2 in Supplemental Table 4) include plants with both these groups. The MSY3 haplotype group was found in a total of nine wild males, all from four subpopulation samples from the north Pacific region of Costa Rica.

All the HSY haplotypes from cultivated hermaphrodites cluster with MSY3 haplotypes, strongly suggesting that the hermaphrodite Y^h chromosome evolved from this region. When the Y-specific SNPs from the same genomic regions are compared in all cases, the set of HSY sequences differ from those of the nine MSY3s by an average of 2729 SNPs, an order of magnitude fewer differences than from MSY1 (31,781) or MSY2 (22,777). To quantify the differences, we compared pairwise F_ST values, which measure variance in allele frequencies between populations in terms of the proportion of diversity between subpopulations relative to the total genetic diversity, as opposed to within them (Charlesworth 1998). High F_ST values indicate population structure, and this can arise either (1) because the populations have been isolated for long enough evolutionary time to have diverged in sequence (e.g., races or subspecies) or (2) because the sequences within populations are very similar in comparison with those from different populations, which is predicted for Y chromosomes, because nonrecombining regions have reduced effective population size (N_e) (Charlesworth 2009). As shown in Figure 3, F_ST for the Y region is low between HSY and MSY3 (F_ST = 0.02) and much higher for MSY1 (F_ST = 0.491) or MSY2 (F_ST = 0.297).

The MSY1, MSY2, and MSY3 populations are not, however, highly isolated, but rather the observed high F_ST values are likely to be a property of the Y chromosomes, which are expected to show strong population structure for reason 2 above. Y-linked sequences are expected to have low within-population diversity, due to genetic hitchhiking processes that occur in nonrecombining regions. In Silene latifolia, for example, Y chromosome population structure is also strong (Ironside and Filatov 2005). This interpretation for the minimal differentiation (F_ST = 0.02) between the HSY sequences from hermaphrodites and the MSY3 sequences is supported by estimates of diversity of Y sequences within populations. Of the total of 66,579 polymorphisms found throughout all the Y chromosomes sequenced (see above), the combined MSY3/HSY set includes only 14,058; 2729 polymorphisms were found only in the wild males, and even fewer (1097) only in the cultivated hermaphrodites. As discussed further below, many fewer variants are completely associated with the difference in sex type; four variants specific to the HSY are found in all hermaphrodites; and only two variants are fixed in the MSY sequences of all males.

We also examined diversity for genes in the recombining regions of Chromosome 1 (the same chromosome whose fully sex-linked region carries the sex-determining genes). Under our interpretation, these recombining regions, the so-called pseudoautosomal region (PAR), should not show strong population structure. Population structure was indeed much less in the PAR. First, F_ST values for PAR sequences are very similar between the HSY sequences and MSY sequences from populations with predominantly MSY1, MSY2, and MSY3 haplotypes (F_ST = 0.11, 0.12, and 0.10, respectively) (Fig. 3), consistent with little population structure for this genome region. The F_ST values for the PAR are statistically indistinguishable between the wild populations (F_ST MSY1 vs. MSY2 = 0.08, MSY1 vs. MSY3 = 0.09, MSY2 vs. MSY3 = 0.09) (Supplemental Table 8), contrasting with the distinct structure observed for the fully Y-linked region. Finally, STRUCTURE, phylogenetic analysis, and PCA analysis of the PAR gene sequences are consistent with the F_ST analysis and show no evidence of population structure correlated with that of the Y-linked region (Fig. 4A–C). Instead, all the analyses of PAR sequences show one set of populations containing all the sequences from males and separate from the hermaphrodites, suggesting that gene flow occurs between the wild papaya populations.

Gene flow is supported by the observation mentioned above of multiple Y sequences within some wild populations. Coexistence of very different Ys is not unexpected. First, even if gene flow brings different Y chromosome haplotypes into the same population, they will remain distinct because recombination between them is suppressed. Second, fixation of one Y haplotype across the species is unlikely given the large geographic region in which wild papaya populations occur, which means that mutations arising in one region will not reach distant populations for many generations.

The HSY is a product of human domestication

The scarcity of fixed differences between MSY3 and HSY sequences suggests very recent divergence. In contrast, the fully sequenced AU9 MSY does not cluster with any of the Costa Rican Y populations and is highly diverged from the other Y-linked sequences (Fig. 5). This divergence led to the previous conclusion that the hermaphrodite–male divergence time is long and greatly predates domestication. With our new sequence results, this is clearly incorrect. Rather, this accession is unusual and is not closely related to the ancestor of cultivated hermaphrodite papayas.

We estimated the MSY3–HSY divergence time using a Bayesian approach in BEAST (Heled and Drummond 2010), using all 105,953 synonymous sites in coding sequences of the HSY and MSY regions. We used a molecular clock of 7 × 10⁻⁹ substitutions per site per year estimated for Arabidopsis, corrected by a factor of 0.672 to take account of the slower molecular evolution in papaya (see Methods). The estimated synonymous site divergence between the HSY and MSY3 yields a time of ∼4000 yr ago, with a 95% highest posterior density (HPD) of 1400–6700 yr ago. This is within the time of origin of agriculture in Mesoamerica (Fig. 5; Pope et al. 2001). The estimate for the HSY from the MSY1 group is 29,000 yr (95% HPD = 24,600–34,510) and from MSY2 is 36,000 yr (95% HPD = 27,000–42,840), both well before the origin of agriculture even in the Fertile Crescent in the Middle East ∼10,000 yr ago (Gupta 2004). We also estimated the divergence time between HSY and the MSY from AU9 using the same molecular clock. This yielded 66,000 yr (95% HPD = 54,600–78,540), similar to the 73,000 yr reported previously (Yu et al. 2008).

Patterns of divergence and diversification between the Y chromosomes and the PAR

If hermaphrodites indeed resulted from human domestication, the HSY should (unlike the wild MSY and the PAR) show evidence for a bottleneck, either due to (1) an initial selection of just one or a few hermaphrodite individuals or to (2) subsequent selection for desirable characteristics. We therefore estimated nucleotide diversity (π). Our new estimates for both Y-linked and other regions are similar to previously published estimates based on Sanger sequencing, indicating that the high-throughput sequencing approach (see Methods) yields reliable estimates (Weingartner and Moore 2012). The π value estimated for our pseudoautosomal sequences is slightly lower in cultivated than wild papaya (0.0017 vs. 0.0020), perhaps due to an initial bottleneck, but the difference is not significant (P = 0.15, based on a Wilcoxon test using all the values estimated for all genes sequenced). In contrast, nucleotide diversity is much lower among the HSY sequences (0.0004) than the wild MSY ones (0.0021, similar to the value for the PAR diversity estimated from the wild populations); the HSY–MSY difference is significant and is also significantly lower than the PAR value (in both cases Wilcoxon test, P < 1 × 10⁻⁵, Fig. 6A). The low diversity specific to the HSY suggests a reduction in diversity during domestication.

Finally, we used Tajima's D-test to assess whether the Y chromosomes show evidence for a severe bottleneck, which should lead to negative Tajima's D-values (Tajima 1989). The HSY has the expected strongly negative mean Tajima's D (mean −1.26, P < 0.05, Wilcoxon test), consistent with a bottleneck or selective sweep during domestication. In contrast, the mean Tajima's D-values across the PAR region of Chromosome 1 are slightly negative, but close to zero in both males and hermaphrodites (−0.017 and −0.154, respectively), and that for the MSY is 0.13, also not significantly different from zero (Fig. 6B).