Unique haplotypes of cacao trees as revealed by trnH-psbA chloroplast DNA

Introduction

The Neotropical cacao tree (Theobroma cacao L.) is widely cultivated in Central and South America, Africa, Indonesia and Malaysia. It is considered an economically important crop because its seeds are used in the chocolate industry (Wood, 2001). Trees can be traditionally classified based on agromorphological traits as Criollo, Forastero and Trinitario (Cheesman, 1944; Toxopeus, 1985). In Mesoamerica, the Criollo cacao has been widely used as food and other purposes as well for nearly 4,000 years (De la Cruz et al., 1995; Whitkus et al., 1998; Powis et al., 2011).

Motamayor et al. (2008) proposed 10 cacao genetic groups based on simple sequence repeat (SSR) analysis. Under this genetic classification, only the traditional group of Criollo has been retained as an accepted genetic group. Forastero contains members of the other nine genetic groups from South America (Motamayor et al., 2008) while Trinitarios are believed to be hybrids of various groups (Yang et al., 2013). South America has been reported to contain the highest genetic diversity of cacao trees.

Conversely, the genetic diversity of cacaos in southern Mexico was reported to be moderate to low in natural populations (Whitkus et al., 1998 using RAPD markers) and cultivated forms (Vázquez-Ovando et al., 2014 using microsatellite markers), although a wide diversity in cacao pod (fruit smooth, rugose, very rugose; apex blunt, apex point; cylindrical and oblate form; basal constriction slight/absent to very pronounced; colours green, red), and seed (elliptics and oblate forms; colour white, slightly pigment and purple intense) morphologies was observed and reported by Avendaño-Arrazate et al. (2010). In the Soconusco farms (Chiapas, Mexico), Vázquez-Ovando et al. (2014) found moderate to high allelic richness and high levels of homozygosity. The authors reported the presence of trees sharing genetic identity with those considered “Ancient Criollo” but also reported the presence of private alleles. These alleles may be associated with commercially relevant phenotypic traits that preserve their relationship with other polymorphic regions of the DNA.

The chloroplast DNA (cpDNA) and its markers have been increasingly used for studies of genetic population structure, evolution, gene flow, haplotype frequency and phylogenetic relationships. Due to its high conservation due to maternal uniparental inheritance, cpDNA is the main data source used for the construction of phylogenetic relationships in plants (Shaw & Small, 2005). In addition, the cpDNA contains variable DNA regions, which makes them useful for studies of population genetics and conservation issues (Shaw & Small, 2005; Shaw et al., 2007). These regions have been widely used to establish phylogeography patterns in alpine species (Wang et al., 2008), to gain insight into the center of origin of cultivated grape populations in Europe (Arroyo-García et al., 2006) and to explain the diversity and population structure of cultivated Chinese cherries (Chen et al., 2013).

Although cpDNA has not been commonly used in cocoa studies, the technique was employed to analyze population genetic variability and to elucidate the complex origins of cocoa varieties. Yang et al. (2011) developed cpSSRs that were subsequently used together with cpSNP markers (developed by Kane et al., 2012) to untangle the genetic origins of the Trinitario cultivar in Trinidad and Tobago (Yang et al., 2013).

The most commonly used cpDNA intergenic spacer is trnH-psbA, which has shown high variability and can be used to elucidate genetic relationships at the intraspecific level (Azuma et al., 2001; Hamilton, Braverman & Soria-Hernanz, 2003). The trnH-psbA region sequences from 10 cacao accessions deposited in the NCBI database produced only one haplotype (Kane et al., 2012), whereas Jansen et al. (2011) reported the presence of polymorphic sites, which produced a different haplotype. The main polymorphisms reported in the noncoding cpDNA region are inversions, transitions and transversions (Whitlock, Hale & Groff, 2010; Zeng et al., 2012). Few studies have reported the presence of insertions or deletions (indels), although indels are probably a common feature in the trnH-psbA spacer (Aldrich et al., 1988).

Nonetheless, the use of indels for diversity and phylogenetic analysis has been questioned (Bieniek, Mizianty & Szklarczyk, 2015; Whitlock, Hale & Groff, 2010) because the mechanism causing indels remains unclear. However, indels are informative characteristics because genetic variability detected using polymorphism due to indels or substitutions can be studied without distinction (Nei, 1987). Therefore, indels are useful markers. Moreover, the inclusion of indels in diversity and phylogenetic analyses enhances the discriminant power between species (Raymúndez et al., 2002; Hamilton, Braverman & Soria-Hernanz, 2003; Kress & Erickson, 2007; Sun et al., 2012) and even between conspecific individuals (Pérez-Jiménez et al., 2013). Therefore, the aim of this study was to evaluate the genetic diversity and describe the relationship between individuals of the Theobroma cacao L. Criollo type of the Soconusco region (Chiapas, Mexico) using the variations in chloroplast DNA revealed by the trnH-psbA spacer sequence.

Material & Methods

Sequence alignment and data analysis

The sequence quality was checked and the electropherograms were edited using BioEdit© (Hall, 1999). Sequences were limited at the ends to avoid the presence of variable sites due to the introduction of sequencing artifacts by the polymerase (approx. 40 bp) and aligned with ClustalW 1.81 (Thompson, Higgins & Gibson, 1994). Visual inspection and manual editing of the sequences was performed to confirm the variable sites. We used two different analytical approaches based on the geographic origin and the genetic origin of the samples (Table 1). In both approaches, molecular diversity indices including the number of segregating sites (S), the number of haplotypes, the haplotype diversity (Hd) and the nucleotide diversity (πd) were estimated following the methods of Nei (1987) in DnaSP© 5.1 (Rozas et al., 2010).

To infer evolutionary relationships at the intraspecific level, we produced a network. The method used was median-joining (MD) based on parsimony criteria (Bandelt, Forster & Röhl, 1999; Polzin & Daneshmand, 2003) and was performed with the software Network© 4.6.1.3 (Bandelt et al., 1995).

Analysis of molecular variance (AMOVA), pairwise Fst values and statistical analyses of molecular variance (F_CT (test performed by permuting individuals within populations), F_ST (test performed by permuting genotypes among populations but within groups) and F_SC (test performed by permuting genotypes among groups)) were estimated using Arlequin^© version 3.0 (Excoffier, Laval & Schneider, 2005). Significance was evaluated by 99,999 random sequence permutations. To determine whether sample sites clustered on a population level, a spatial analysis of variance (SAMOVA) was conducted (Dupanloup, Schneider & Excoffier, 2002) using haplotype data and the geographic coordinates of each of the 5 sample sites. The SAMOVA was run for K = 2–5 putative populations to determine the maximum F_ST value and the highest proportion of differences between populations due to genetic variation.

The neutral evolution of chloroplast DNA was evaluated to examine whether any population had experienced historic demographic changes using Tajima’s D test (Tajima, 1989) with Arlequin^© version 3.0 (Excoffier, Laval & Schneider, 2005). It was computed for all seven geographic populations and overall without a priori populations designation; p-value were generated using 1,000 simulations under a model of selective neutrality.

Seven accessions from the NCBI database were included as references in the genetic origin approach analysis: MATINA 1/6 ( HQ336404.2), CRIOLLO-22 ( JQ228379.1) AMELONADO ( JQ228380.1) SCA 6 ( JQ228382.1), ICS 1 ( JQ228381.1), ICS 6 ( JQ228383.1) and ICS 39 ( JQ228387.1).

Table 2:

Nucleotide polymorphic sites and cpDNA haplotypes in cacao populations based on variation in the intergenic trnH-psbA spacer region.

Haplotype	Polymorphic site						Populations (see Table 1 note)
	22	134	206	309	310	487	Pop1	Pop2	Pop3	Pop4	Pop5	Pop6	Pop7	Pop8	Pop9
H1	–	T	–	A	A	A	3	2	3	2	5	1	3
H2	C	T	–	A	A	A	1		1		2			1
H3	C	T	A	–	–	–			1		1
H4	–	A	–	A	A	A				1
H5	–	T	–	A	A	–						1
H6	C	T	–	A	–	A	1
H7	–	T	–	A	–	–				1				1
H8	–	T	A	–	–	A				1				1
H9	–	T	A	–	–	–				1	1			1
H10	–	A	–	A	A	–		1
H11	–	T	–	A	–	A								1	1
H12	–	T	A	A	–	–									6

DOI: 10.7717/peerj.1855/table-2

Results

Intraspecific relationships

Figures 2 and 3 show the haplotype networks built with data from the geographic (Fig. 2) and genetic approaches (Fig. 3). The individuals belonging to each haplotype are also included. The general base has a common haplotype for the two networks (H1) that includes cacao trees from all geographic populations (Fig. 2) and four of six groups based on the genetic approach (Fig. 3). A unique set of intermediate haplotypes are derived from this common haplotype (H1) and are interconnected by one or two mutational steps in both networks. The H4–H6 haplotypes were farthest from the central clade (i.e., newly created haplotypes; Figs. 2 and 3). Haplotypes H3–H6 were singletons.

Discussion

In this study, high haplotype variation was found in the chloroplast DNA from cacao trees grown in the Soconusco region. No inversions or transitions were found, although they were reported to be common in other plants (Whitlock, Hale & Groff, 2010; Zeng et al., 2012). However, we found five insertions or deletions (indels) in three poly-A regions and one A ↔ T transversion (Fig. 1). This result agreed with the findings reported by Jansen et al. (2011) in the MATINA 1/6 accession and supported the affirmation by Aldrich et al. (1988) that indels were a presumably common feature in the trnH-psbA region. In the data analysis, we included the indels as informative character states, and the high interspecific divergence of the spacer region allowed their use as a marker for DNA barcoding (Kress & Erickson, 2007). The molecular diversity indices determined in the present study were similar to the results of Zeng et al. (2012) using the same intergenic spacer, which revealed 11 haplotypes for 35 Thinopyrum intermedium samples, low nucleotide diversity (πd = 0.00473) and moderately high haplotype diversity (Hd = 0.7331) (our results for the geographic populations were πd = 0.0032 and Hd = 0.66). The results of those authors supported the use of one intergenic spacer to reveal nucleotide polymorphisms similar to our findings.

Our haplotype diversity results are contrary to those reported by Vázquez-Ovando et al. (2014). These authors reported low genetic diversity in individuals from the same region (in particular Population 4 in Mazatán) using microsatellite markers. One reason for the discrepancy may be that a larger number of individuals with Criollo ancestry was included in that study, resulting in a higher degree of homozygotes and lower population genetic diversity. Our study also included individuals from other cacao varieties that possessed greater genetic diversity, at least at the nuclear DNA level. However, the low nucleotide diversity found in this study was supported by the low genetic variability found using nuclear microsatellites. Individuals included in both studies showed great morphological pod variability that resembled the Criollo type (e.g., different degrees of roughness, color and deep grooves) reported previously by Avendaño-Arrazate et al. (2010). This finding could reveal a greater association between the morphological variability of the cacao pod with the reported allelic richness (Vázquez-Ovando et al., 2014) and the presence of polymorphic sites in several trees found in our study.

The number of haplotypes was higher than the number of polymorphic sites (Table 2). This feature is associated with ancestral species that have sufficiently diverged to accumulate mutations among different haplotypes (Roger, 1995). The haplotype number detected in the present study is unusually striking compared with other works. For example, Yang et al. (2013) found only three haplotypes based on three cpSNP markers. However, that study exclusively analyzed nucleotide substitutions, whereas in this study five indel regions were included; this difference may explain the high haplotype diversity found here. Indels have been reported to have a high mutation rate compared with other regions of the cpDNA (Ingvarsson, Ribstein & Taylor, 2003), especially when they are repeated locally (Yamane, Yano & Kawahara, 2006) such as in region 309–310 of our sequences (Fig. 1).

Several explanations are possible for the presence of more than one Criollo haplotype. First, only the maternal line gave rise to the eight Criollo haplotypes by mutation. Second, the “Criollo” phenotype had multiple provenances, indicating that the ancient haplotypes persisted over time in the Soconusco cacao farms. Third, some samples were misclassified as “Criollo” (especially MAJH02 and Carmelo, which were the most divergent “Criollo” individuals; haplotypes 4 and 6, respectively, Fig. 3). These samples possibly belong to the Admixture group rather than the Criollo. However, they are also contenders for the Modern Criollo group (i.e., individuals classified as Criollo that might have been introgressed with Forastero genes) (Motamayor et al., 2002) and preserve phenotypic traits of the ancient Criollo. Finally, heteroplasmy and haplotype polymorphisms of plastid genomes within and among individuals were documented in Malvaceae (Wolfe & Randle, 2004). These phenomena could be present in Theobroma cacao. To test those hypotheses, additional studies are needed using high-throughput sequencing of chloroplast genomes.

Population 7 (Selva Lacandona) exhibited no haplotype diversity (Hd = 0). However, haplotype H1 located in this population is considered the common ancestor because it is shared by all populations (Fig. 2). In contrast, the two individuals belonging to Population 6 (Yucatán), which exhibited different haplotypes (H1 and H5) from one another, were interrelated by only a mutational step (Fig. 2). This result shows that an individual tree belonging to a Yucatán population (as well as all other haplotypes) eventually descended from other individuals in this region where the Maya people grew cacao.

The low nucleotide polymorphism levels could be explained by rapid population expansion events in the distribution range, whereas high haplotype diversity might be due to the continuous introduction of individuals from different locations. However, these results should be interpreted with caution because, as stated above, they are limited by the sample’s size. Populations recently introduced or expanded from a small number of founders would have a common haplotype shared by most individuals and many rare haplotypes connected to the main population by a few independent mutations (Slatkin & Hudson, 1991; Avise, 2000) such as observed in the present study (Fig. 2). A similar argument was proposed based on the use of microsatellite markers (Vázquez-Ovando et al., 2014).

The relatively low variability in the cultivated cacao populations was supported by the lack of neutrality revealed by the global Tajima test. Specifically, the negative Tajima’s D value (−0.93302) in Population 4 (Mazatán) could be related to a “bottleneck” event, which would indicate population expansion and not natural expansion because it was a cultivated population. The occurrence of unclear events in the past (disease, volcanic eruptions or other natural events) may have caused the almost complete disappearance of populations established by the people in the Mesoamerican region (De Sahagún, 2009, Codex Florentino). Rapid expansion due to recolonization of the populations and the probable introduction of other varieties of cacao trees not native to the region would have subjected the populations to a bottleneck events in very recent periods. However, these are presumptive weak inferences of the population history based on a single locus. The bottleneck event could also be related to the loss of alleles (haplotypes; especially rare alleles), which is much greater than the loss of genetic variance per se. Although these rare alleles contribute little to the total genetic variability, they can provide unique responses against evolutionary challenges similar to the high number of unique haplotypes found in this study (3 singletons in Population 4). The presence of both common and rare haplotypes can be the result of a directional-purifying selection process or expansion events from small populations (Hedrick, 2005). The H3–H8 haplotypes (cultivated populations) are singletons. This finding agreed with Crandall & Templeton (1993), who reported that the singletons tend to are connected to haplotypes from the same population. Population 4 (Mazatán) shows the highest haplotype diversity, which makes this population an important reservoir of genetic material at the chloroplast and possibly phenotypic levels based on the abundance of pod morphologies observed in this population.

Overall, cacao trees with high Criollo ancestry were located in the center of the haplotype network. This result was supported by the coalescence theory that predicted that the ancient haplotype should be the most common and most distributed among the populations. In concordance, derived haplotypes would be less frequent and in many cases would be private; these haplotypes would be located in regions containing the latest cultivated cacao populations. The H10 and H9 haplotypes may have been recently created because they are located at the ends of the network, possibly due to germplasm exchange with traits of interest to cacao farmers. These anthropogenic activities may have had a strong impact on the levels of variation observed in the cpDNA sequences, which explains the observed lack of differentiation. Additionally, migration over long distances and the exchange by farmers contributed to the colonization of new regions founded by a few individuals, thereby establishing different alleles via mutation and genetic drift.

Furthermore, the F_ST = 0 value determined by AMOVA revealed that all of the molecular variance occurred within populations. Indeed, the SAMOVA F_ST value (Table 4) was not sufficient to show at least moderate differentiation between populations (F_ST ≥ 0.05). This finding provides some explanations regarding the demographic history of T. cacao trees, indicating that the populations formed a priori and experienced gene flow, resulting in population homogenization. The spatial analysis revealed the highest differentiation between groups when K = 2 was tested; K = 3 (F_ST = 0.00088) grouped trees from the Yucatán, Selva and Cacahoatán in the same genetic population. This grouping is unusual because the geographic distance is longer among the three localities and may be associated with the distribution of trees in the past (i.e., the ancestral haplotype (H1) grouped individuals from Selva; one mutational step resulted in the origination of the individuals from the Yucatán, which in turn originated the individuals at Cacahoatán by the same event (Fig. 3)). Following this criterion, H4 and the non-Criollo trees belonging to H1 have a greater correspondence with the Criollo genotype, although it was previously reported to be an Admixture and non-Criollo, respectively (Vázquez-Ovando et al., 2014).

Conclusions

Indels and one transversion located in the chloroplast DNA trnH-psbA spacer region of cacao trees can distinguish individuals that are indistinguishable in other marker systems or separated by only few SSR markers, and further support use of these cpDNA markers.The molecular analysis showed low nucleotide diversity but high haplotype diversity possibly due to population bottleneck events. These results were confirmed by the negative Tajima’s D and the arrangement of the haplotype network. We identified 10 different haplotypes (from cultivated trees) of which H3–H8 resulted in singletons because they were not associated with other cacaos or with those reported in the molecular databases. The presence of these haplotypes accompanied by the low number of mutational steps might suggest a very short evolutionary history or events that led to disappearing-expanding populations in southern Mexico. These results suggest confirmation of selection of fruits from few mother trees (even as few as one) that were moved by human agents from South America into Mexico and that the Criollo complex may be homogenous based on maternal influence. Increasing sampling numbers would go a long way in establishing whether a true sub-structure of maternal origin exists. One geographic population (Pop 4, Mazatán) consisted of high frequency haplotypes, which makes this zone an important reservoir of genetic material at the chloroplast and possibly phenotypic levels because an abundance of pod morphology was also observed in this population. The genetic differentiation between populations was zero, suggesting that gene flow homogenized the populations.

Supplemental Information