Relationship between fruit phenotypes and domestication in hexaploid populations of biribá (Annona mucosa) in Brazilian Amazonia

Plant material

The fruits were purchased at fairs in central Amazonia, in the state of Amazonas (Manaus and Rio Preto da Eva), along the middle (Tefé) and upper Solimões River (Benjamin Constant, Tabatinga and Atalaia do Norte) in western Amazonia (Table 1 and Fig. 1A) (SISBIO authorization 70846-1 and SisGen registration A19391B). We classified the fruits found into five different morphotypes based on the carpel protrusions on the exocarp (Fig. 1B), as done for cherimoya (A. cherimola) (Vanhove, 2008): smooth fruits; fruits with small carpel depressions; fruits with small carpel protrusions; fruits with medium carpel protrusions and fruits with large carpel protrusions.

Chromosome counts

Chromosomal preparations were performed according to Guerra & Souza (2002), with some adaptations for A. mucosa. The seeds were germinated, and the root apices pretreated with 8-hydroxyquinoline (8 HQ) 2 mM for 24 h in a refrigerator. Roots were fixed in Carnoy (ethyl alcohol: acetic acid, 3:1) for about 6 h at room temperature (~25 °C) and then stored in a freezer (−20 °C). Slides were prepared with enzymatic digestion (1% macrozyme, 2% cellulase and 20% pectinase) for 40 min at 37 °C.

After 5 days, we used the Schweizer (1976) protocol to band the chromosomes in the slides, with the following modifications. Slides were covered with 12 µL of chromomycin A3 (CMA; 0.5 mg/ml) in a darkroom for 1 h, then washed with distilled water and covered with 12 µL of 4–6-diamidino-2-phenylindole (DAPI; 2 mg/ml) in a dark chamber for 30 min. Then, slides were mounted with 12 µL of glycerol/McIlvaine/MgCl₂ and covered with a glass coverslip. The slides were then stored for at least three days until observation under an Olympus BX51 microscope. Images were captured using the Olympus Cell View program (Olympus Corporation, Shinjuku City, Japan). Adobe Photoshop CS6 (Adobe Systems, Inc., San Jose, CA, USA) was used to assist in balancing colors, contrast and brightness, and assembling the chromosome image.

Genome size estimations by flow cytometry

Young leaves of 15 A. mucosa individuals were collected and used for genome size measurements by flow cytometry. Approximately 1 cm² of young leaf tissue was minced with fresh leaf material from an internal standard (Petroselinum crispum (Mill.), Apiaceae, 2C = 4.50 pg). Nuclear suspensions were prepared according to Doležel, Greilhuber & Suda (2007), adding 1.5 mL of WPB (woody plant buffer) (Loureiro et al., 2007) to the nucleus extract, and then the nuclear suspension was filtered through a 30 μm nylon mesh. The nuclear suspension was stained with 50 μL of 1 mg/mL propidium iodide, incubated for at least 10 min and analyzed by flow cytometry. At least 5,000 nuclei from three replicates were analyzed for each sample using a CyFlow Ploidy Analyzer (Sysmex, Kobe, Japan) cytometer. Each set of histograms from the flow cytometry analyses was analyzed using Flowing Software v2.5.1 by Perttu Terho (Turku Center for Biotechnology, University of Turku, Turku, Finland).

Statistical relationship between fruit weight, morphotypes and genome size

To test whether the samples originated from the same distribution and to analyze the relationship among genome size, morphotype and fruit weight of A. mucosa individuals (non-normal data), we performed the Kruskal-Wallis test (McKight & Najab, 2010), a nonparametric ANOVA, in R package MultNonParam (Kolassa & Jankowski, 2014).

Sequence editing, alignment and phylogenetic analysis

We sampled 50 taxa of the genus Annona available in GenBank (Benson et al., 2012), which represents 31% of the accepted species of the genus (Table S2). Among the 50 taxa of Annona, 16 were from the Rollinia clade (38% of the 42 species in Rollinia). The species Asimina incana (W. Bartram) Exell. was used as an outgroup following Li, Thomas & Saunders (2017). The plastid intergenic spacer psbA-trnH and the plastid genes matK, ndhF, rbcL and trnL were used, totaling 161 sequences (Table S3). The DNA sequences were aligned using MUSCLE (Edgar, 2004) as an extension of Geneious v.7.1.9 (Kearse et al., 2012). Subsequently, the sequences of markers for each of the species were concatenated for analysis, as done by Guo et al. (2017) in the major phylogeny for the Annonaceae.

We used jModelTest v.2.1.6 to assess the best DNA substitution model for each individual marker (Darriba et al., 2012) evaluated with the Akaike Information Criterion (Akaike, 1974). The best fitting models were: rbcL (K80+I), trnL (HKY+I), matK (TPM1uF+G), ndhF (HKY+I), psbA-trnH (HKY+G) (Table S3). Phylogenetic relationships were inferred using Bayesian Inference with four independent runs with four Markov Chain Monte Carlo (MCMC) runs, sampling every 1,000 generations for 10,000,000 generations, implemented in MrBayes v.3.2.6 (Ronquist et al., 2012). The analysis was evaluated in TRACER v.1.6 (Rambaut & Drummond, 2014) to determine if the estimated sample size (ESS) was greater than 200, after applying a burn-in of 25%. The consensus tree was viewed and edited in FigTree v.1.4.2. (Rambaut, 2009). All phylogenetic analyses, excluding the ESS, were performed on the CIPRES Science Gateway (Miller, Pfeiffer & Schwartz, 2011).

Estimates of divergence time using a molecular clock

For the 22 taxa with chromosomal information (Table S1), we repeated the phylogenetic analyses following the same methodology described above. We estimated the divergence time between species in BEAST v.1.8.3., with input preparation by BEAUti (Drummond & Rambaut, 2007; Drummond et al., 2012). We used the “Relaxed Clock: Uncorrelated Lognormal” (Drummond & Rambaut, 2007) and the “Yule” specification model (Gernhard, 2008). Two independent runs of 200,000,000 generations were performed, sampled every 10,000 generations. To verify the effective sampling of parameters and to assess the convergence of the independent chains, we examined the posterior distributions in TRACER v.1.6. (Rambaut & Drummond, 2014). Markov Chain Monte Carlo (MCMC) sampling was considered sufficient at effective sample sizes (ESS) greater than 200. After removing 10,000 samples as burn-in, independent runs were combined, and a maximum clade credibility tree (MCC) was constructed using TreeAnnotator v.1.8.2. (Drummond et al., 2012). BEAST analyses were performed on the CIPRES Science Gateway (Miller, Pfeiffer & Schwartz, 2011). The tree calibration used secondary calibrations, based on the work of Li, Thomas & Saunders (2017); for the node of the monophyletic genus Annona (29 Mya) and the node of the group Annonae (53.3 Mya) a standard deviation of two was used for both.

Reconstruction of ancestral character state

To reconstruct the ancestral ploidy state, we obtained chromosome numbers (n) of Annona species from The Chromosome Counts Database (Rice et al., 2015) and the Index to Plant Chromosome Numbers (Goldblatt & Johnson, 1994). The genome size data (1C) were taken from the Plant DNA C-values Database (Leitch et al., 2019) and literature (Table S1).

For the 22 species with chromosome number information (Table S1), the evolution of haploid chromosome numbers was inferred using two approaches. The first analysis was performed using ChromEvol v.2.0 (Glick & Mayrose, 2014). This program determines the likelihood of a model explaining the given data within the phylogeny based on the combination of two or more parameters. Ten models were tested with different sets of program parameters. The models were fitted to the data, each with 10,000 simulations, and the best fit model was Base_Number_No_Dupl selected using the AIC (Table S4). In addition to the AIC estimates, we tested the model suitability for probabilistic models of chromosome number evolution according to Rice & Mayrose (2021). The second analysis was performed using R package phytools (Revell, 2012), considering the haploid number of species (n = 7, 14, 21 and 28). The characters were reconstructed using function “ace” of the R package Ape (Paradis, Claude & Strimmer, 2004) with the EqualRates model. In cases where multiple chromosome numbers were reported for a given taxon, the model number was used (Glick & Mayrose, 2014; Salman-Minkov, Sabath & Mayrose, 2016).

Fruit morphometrics and karyotypes

The weights of A. mucosa fruits ranged between 100–1,850 g (Table 1), with the largest fruits from the upper Solimões River region in western Amazonia. The morphotypes based on the classification of carpel protrusions were variable. The most frequent morphotype in the Manaus region had small carpel protrusions (58%), while in the upper Solimões River region the most frequent morphotype had large carpel protrusions (89%). Smooth and small depression morphotypes were rare (Table 1, Figs. 1A and 1B). In the Kruskal-Wallis analysis, there was no significant relationship between fruit size and morphotype, only a tendency for larger protrusions to occur in larger fruits (p = 0.0646).

All the 31 fruits that were karyotyped had the same chromosome number of 2n = 6x = 42 (Table 1 and Fig. 2). The distribution of constitutive heterochromatin (HC) was analyzed quantitatively. Five pairs of bands with a strong signal of CMA+/DAPI− were observed, three in regions close to the centromere (pericentromeric) and two in the short arms of chromosomes (with distension). There was no discernable relationship among the sizes, morphotypes and HC distributions.

The genome sizes of the 15 A. mucosa individuals analyzed revealed little variation and no significant differences between individuals (2C mean = 5.31 ± 0.12, minimum = 5.21 pg, maximum = 5.65 pg) (Table 1 and Fig. S1). The genome size was not related to fruit size (p = 0.2267) or morphotype variation (p = 0.6437).

Phylogenetic relationships, diversification rate and polyploid state reconstructions in Annona

With the 50 species of Annona examined with five molecular markers (Table S3), the genus was recovered as monophyletic (PP 100%), as was the Rollinia clade (PP 96%) (Fig. S2). The tree topology shows that A. mucosa and A. cuspidata ((Mart.) H. Rainer) shared the same most recent common ancestor (MRCA) (PP 98%) (Fig. S2).

With the 22 species of Annona for which both DNA sequences and chromosome counts were available, the tree topology was similar (Figs. 3A and 3B; Fig. S3) and A. mucosa and A. cuspidata also shared the same MRCA (PP 93%). The BEAST analysis estimated that the divergence of genus Annona occurred in the Paleocene (59.2 Mya: 95% CI [58.3–60.2]), but its diversification appears to have started in the Oligocene when the three largest groups diverged (29 Mya: 95% CI [28–30]). The BEAST analysis estimated that during the Oligocene the Rollinia clade began to diversify (23 Mya: 95% CI [17.8–28.4]), and A. mucosa diverged in the Pliocene (3.4 Mya: 95% CI [0.3–9.2]).

The reconstruction using ChromEvol revealed x = 7 as the ancestral chromosome number for the genus Annona (Fig. 3A), as did the reconstruction using the R package Phytools (Fig. 3B). With the model of Chromevol, we found one independent genome duplication event in the MRAC of A. mucosa and A. cuspidata, i.e., x = 14. With the model of Phytools, we found two independent genome duplication events in the MRAC of A. mucosa and A. cuspidata, with similar probabilities, i.e., x = 14 and x = 21. The probability of a polyploid ancestor for the Rollinia clade is low in both approaches.

Evolutionary events and their consequences for the phenotypic variation of A. mucosa

A primary potential advantage of polyploidy for domestication of fruit crops is fruit size variation (Salman-Minkov, Sabath & Mayrose, 2016; Akagi et al., 2022). Since polyploidy is associated with gigantism of reproductive parts in plants (Hancock, 2012), the ploidy level was expected to be associated with different fruit sizes, i.e., smaller fruits should present a smaller ploidy level and as the fruit size increases, the ploidy level would increase proportionally. This difference in ploidy may also represent different stages of domestication or different domestication syndromes in species of a genus (Akagi et al., 2022). This association has been observed in some fruit crops, for example, in the domestication of strawberry (Bertioli, 2019) and star fruit (Hu et al., 2021). However, this hypothesis was not confirmed for Annona mucosa. Within Annona, plants with very large fruits, such as domesticated soursop (A. muricata), are diploid, and plants with small fruits, such as undomesticated Araticum do mato (Annona sylvatica A. St.-Hil.), are octaploid (Table S1), demonstrating that wild species have variable ploidy levels also. Thus, we infer that human selection in Annona mucosa influenced fruit size regardless of the ploidy level.

The other hypothesis raised concerns the presence of variation in genome size related to different phenotypes, which was also refuted. Even in different regions of Brazilian Amazonia, the variation in genome size was not significant, suggesting stability in the genome. The genome size reported here is similar to Lorenzoni (2016): 2C = 5.42 pg (±0.12) with 15 accessions of unknown origin of hexaploid cultivated plants maintained at the Federal University of Espírito Santo (Alegre, ES, Brazil). The other two studies showed divergences in the size of the A. mucosa genome (Soares et al., 2014; Leitch et al., 2019), but methods, sample sizes and provenances are unclear.

Martin et al. (2019) established a DNA amount for diploid individuals in Annona and associated increased ploidy with genome size. Compared to A. neosalicifolia H. Rainer, a hexaploid species (2C = 9.64 pg), A. mucosa has approximately four pg less DNA content. Based on the genome sizes of other Annona species (Table S1), we observed that for A. mucosa cultivated in Brazilian Amazonia this characteristic is not linearly associated with ploidy, as all samples have 6x = 42 and 2C = 5.31 ± 0.12. Genome size reduction is a recurrent phenomenon in polyploid plants, which could explain this reduction in A. mucosa (Bennetzen, Ma & Devos, 2005; Šmarda & Bureš, 2010). The reduced 2C mean value in higher ploidy levels may be due to the loss of repetitive DNA, as in Psidium cattleyanum Sabine (Myrtaceae), a neotropical polyploid complex frequently managed by humans (Machado & Forni-Martins, 2022).

The species also does not follow the pattern of heterochromatic regions of the genus Annona, presenting five pairs of bands, while the expected would be three pairs, based on diploid and tetraploid species (Morawetz, 1986b). In polyploidy, the larger number of chromosomes offer a greater chance of gene gain and recombination (Grant, 1974; Guerra, 1988; Hancock, 2012), which could explain the difference in the pattern of the heterochromatic regions.

Polyploidization events and the evolutionary context of the domestication of A. mucosa

Despite the under representation of our plastidial phylogenetic hypothesis (with 50 species in general and 22 species for chromosome number reconstruction), our data corroborate the existing phylogenies published for Annona (Chatrou et al., 2009; Guo et al., 2017). This issue of lower representation is due to the limitation of available ploidy data, which is relatively common in studies involving compilation of molecular data and chromosomal counts or genome size (Sader et al., 2019; Ibiapino et al., 2022). Even though our results recovered all the main clades within Annona with high statistical support, our lower taxon representation may have favored this high internal resolution because of the increasing ratio of molecular data to the number of terminals, as already mentioned for mega diverse neotropical genera (Amorim et al., 2019). As well as the phylogenetic relationships, our divergence time estimations agree with what is known for the genus Annona and the divergence of A. mucosa and A. herzogii in the Pleistocene (Li, Thomas & Saunders, 2017).

The species of the Rollinia clade, which had previously presented polytomy (Chatrou et al., 2009; Guo et al., 2017), had good statistical support for the branches in our analysis. This clade resolution allowed us to suggest that A. mucosa and A. cuspidata share the same MRCA given our sample of 38% of the species in the Rollinia clade (Fig. S2). The ancestral chromosome number reconstructions presented in our study are divergent in relation to the polyploid ancestry of the Rollinia clade. Similar divergences have been reported in other plant groups and, depending on the methodology, different chromosomal events are suggested (Brottier et al., 2018; Moraes et al., 2020). With our results it is not possible to infer whether polyploidy is a consequence of a single or multiple ancestral polyploidization events within the Rollinia clade because we sampled only 16% of the clade’s species with chromosome data in the reconstruction (Figs. 3A and 3B).

On the other hand, our results suggest (in both approaches) that the MRCA of A. mucosa was polyploid. Because this ancestor was estimated to have diverged in the Pliocene, before human arrival in the Americas (Clement et al., 2021), the polyploidy event is inferred to have occurred well before the domestication process began. This temporality of events has already been observed in guarana (Paullinia cupana Kunth var. sorbilis Ducke), a species of vine domesticated in central Amazonia (Urdampilleta, Forni-Martins & Ferrucci, 2020), and in numerous other plant species (Salman-Minkov, Sabath & Mayrose, 2016; Zhang, Wang & Cheng, 2019).

Because of our inference about the polyploid ancestor of A. mucosa, it is expected that wild populations of A. mucosa will be polyploid. Three individuals in possibly non-anthropic environments in Peru were recorded as tetraploid (Maas et al., 1992), but only one sample was described as collected in primary forest. Additionally, an individual in an apparently wild population of A. mucosa in Veracruz, Mexico, presented hexaploidy (Serbin et al., 2022). This result reinforces the hypothesis of the presence of polyploidy before the domestication process and that polyploidy was not the main evolutionary force driving human selection of A. mucosa.