Genetic Diversity of Yuca (Manihot esculenta esculenta; Cassava, Manioc), an Indigenous Crop in the Peruvian Amazon

Abstract

Yuca (Manihot esculenta esculenta; cassava, manioc) is a native Amazonian crop represented by myriad landraces. To investigate human influences on its diversification, we conducted field observations and analyzed 13 short tandem repeat (STR) loci in 43 landraces in the Peruvian Amazon. We found a different multilocus genotype (MLG) in every landrace. However, tests for Hardy–Weinberg equilibrium found a deficit of heterozygosity at every locus (p < 0.001 for 12 of 13 loci). Further, the fraction of genetic variance due to landrace differences was greater than expected (38.84%; p = 0.001). This suggested that landrace hybridization is restricted, a finding consistent with our field observations. However, we found an excess of within-landrace heterozygosity (p < 0.001) in 39 of 43 landraces, suggesting they originated through hybridization. Mantel tests identified associations between genetic and geographic distances (p < 0.001), but their correlation coefficients were low (Mantel’s r < 0.21). In addition, AMOVA analyses revealed that differences between landraces collected from five sampled rivers accounted for just 3.05% of observed genetic variance (p < 0.001). Neighbor joining and principal components analyses also revealed little evidence of differentiation between rivers. Finally, in a comparison with a secondary sample, we found that the closest relative of 27 of 28 specimens had a landrace name different from their own, suggesting that traditional nomenclature is a poor indicator of genetic relatedness.

1. Introduction

Yuca (Manihot esculenta esculenta; also called manioc or cassava) is a field crop central to the diet of peoples throughout the New World tropics [1,2]. Genetic evidence suggests that it was domesticated ~10,000 years ago on the southern margin of the Amazon basin, where it was derived from a wild M. esculenta subspecies, M. e. flabellifolia [3]. Paleoethnobotanical evidence indicates that it then spread rapidly through the lowland humid zones of South and Central America, reaching sites in Panama as early as ~8000 years ago [4,5]. Yuca is also common in tropical Africa, Asia, and Oceania, where it arrived via trade from Brazil in the 16th century [6,7].

Yuca’s success as a crop derives from its combination of productivity, robustness, and diversity [8]. Its potato-like tubers store large quantities of starch that can be extracted using basic manual techniques, and it can be grown successfully even where soils are poor. Yuca is also naturally pest resistant, harboring a toxic defense system that deters herbivores but is easily neutralized with simple processing [9]. It is biologically diverse as well. Across its range, yuca is represented by hundreds of landraces differing in morphology, maturation time, toxicity, and other traits [10,11,12,13,14,15,16]. Their geographic distributions overlap extensively, and local areas can harbor dozens of varieties [17,18]. Further, growers readily recognize local yuca types, and usually assign traditional names to them [10,12,13,14,15,19,20]. This enables growers to maintain landraces suited for different purposes and growing conditions [18,21,22,23,24]. For instance, some varieties are used exclusively for the production of fariña (cassava meal), while others are grated and baked, and others are chopped up and roasted. Many can be used for multiple purposes. These assets have established cassava as a staple food across indigenous Amazonia, where it accounts for >50% of caloric intake in some populations [25,26,27].

Yuca’s prevalence in the New World tropics, together with evidence that it originated on the southern boundary of the Amazon Basin, raises questions about the relationships between human influences, phenotypic variation, and genetic diversity in the crop. Past studies have revealed important pressures. For instance, although yuca and its putative ancestor, M. e. flabellifolia, hybridize readily, genetic variation in yuca is a subset of that in M. e. flabellifolia, even where their geographic distributions overlap [28]. Together with the observation that yuca is rarely found away from human settlements, this suggests that hybridization between M. e. flabellifolia and M. e. esculenta has been restricted artificially, by humans [28]. Observations on local scales also reflect human impact. For instance, diversity in newly planted yuca chacras is maintained in some localities by traditional propagation practices, with growers establishing new chacras mainly using clonal cuttings from old chacras rather than seedlings, which are often inbred [18]. Further, growers selectively grow different yuca types in different landscapes, particularly differentiating between tierra firme (rarely flooded high ground) and várzea (seasonally flooded low ground) [12,29]. Thus, growers actively control both genes and environment in their management of the crop.

In a previous study, we investigated morphological, physiological, and geographic variation in yuca on five tributaries of the upper Amazon River, in Perú [17]. We identified 45 traditional landraces, and analyzed patterns of phenotypic diversity among them, including comparisons within and between rivers. We found that landraces overlapped substantially with respect to individual traits, such as leaf area and starch content, but were distinct overall. We also found that the tributaries we sampled showed little evidence of phenotypic differentiation. These patterns raised questions about genetic diversity. For instance, while patterns of phenotypic variation suggested that genetic variation was present, they did not reveal its extent. They also did not reveal the evolutionary relationships among landraces, which are potentially resolvable using genetics. Additionally, it was unclear whether the observed lack of phenotypic differentiation among rivers reflected a lack of genetic differentiation or was due to some other factor, such as environmental homogeneity. Whether the traditional nomenclature of landraces accurately reflected their relatedness, which could shape growers’ ability to manage the crop, was uncertain as well. In this study, we used genotypic data from 13 short tandem repeat (STR) loci from 43 landraces to address these issues.

2. Materials and Methods

Details of our study site and collection methods are provided in Wooding and Payahua [17]. In brief, our study was conducted on the Amazon River in northeastern Perú, on tributaries surrounding the city of Iquitos (Figure 1). Ecologically, the region is dominated by lowland tropical rainforest and whitewater rivers, including numerous tributaries of the Amazon’s main channel. The majority of the region is inhabited by scattered communities of smallholders, who grow yuca as a subsistence crop. Our sampling focused on five rivers, the Itaya, Nanay, Orosa, Pintuyacu, and Tahuayo. Four of the five are direct tributaries of the Amazon. The Pintuyacu is a tributary of the Nanay.

Our field observations and specimens were collected in riverfront communities, which we visited by boat. We identified collection sites opportunistically by traveling up rivers, focusing on villages with 50–100 residents. A major goal in our study was to determine the breadth of cassava diversity in the region. Therefore, we sought to identify as many landraces as possible, and conversations with growers were specifically aimed at finding novel varieties. When growers reported cultivating a landrace not yet encountered in our study, we offered to purchase a specimen. We then visited their chacras for observation, measurement, and collection. In total, our visits identified 45 named landraces (Figure 2) [17]. In addition, we opportunistically collected 17 specimens representing putative duplicates of the 45 landraces in our core sample. For instance, in addition to our first specimen of the amarilla landrace, which we collected on the Rio Tahuayo, we opportunistically collected five others referred to as amarilla by growers on the Itaya, Orosa, and Nanay Rivers.

Cassava is a diploid subspecies with 2N = 36. We extracted DNA from dried, powdered leaf samples using the EZNA Plant DNA Kit (Omega Bio-Tek, Cat:D3485-01) following the manufacturer’s instructions for dry plant material. The extracted DNA was then checked for quality using a NanoDrop ND-8000 Spectrophotometer and quantitated using the DeNovix QFX Fluorometer with the Qubit dsDNA HS Assay (Invitrogen, Cat: Q32854). These analyses successfully extracted high-quality DNA from 43 of the 45 landraces in our core sample, and failed with two (crema and motelo), which were excluded from further consideration in this paper. PCR reactions were performed for 13 STR markers previously identified by Chavarriaga-Aquirre et al. [30] and Mba et al. [31], using their published cycling profiles (

Table 1. STR markers and primers.

Marker	Forward Primer	Reverse Primer
GA12	GATTCCTCTAGCAGTTAAGC	CGATGATGCTCTTCGGAGGG
GA13	TTCCCTCGCTAGAACTTGTC	CTATTTGACCGTCTTCGCCG
GA21	GGCTTCATCATGGAAAAACC	CAATGCTTTACGGAAGAGCC
GA126	AGTGGAAATAAGCCATGTGATG	CCCATAATTGATGCCAGGTT	Chavarriaga-Aquirre et al. [30]
GA127	CTCTAGCTATGGATTAGATCT	GTAGCTTCGAGTCGTGGGAGA
GA131	TTCCAGAAAGACTTCCGTTCA	CTCAACTACTGCACTGCACTC
GA136	CGTTGATAAAGTGGAAAGAGCA	ACTCCACTCCCGATGCTCGC
GA161	TGTTCTTGATCTTCTGCTGCA	TGATTGTGGACGTGGGTAGA
SSRY32	CAAATTTGCAACAATAGAGAACA	TCCACAAAGTCGTCCATTACA
SSRY46	TCAGGAACAATACTCCATCGAA	CGCTAAAGAAGCTGTCGAGC	Mba et al. [31]
SSRY70	CGCTATTAGAATTGCCAGCAC	CGCTTGTTGTATCCATTGGC
SSRY83	TGGCTAGATGGTGATTATTGCTT	TGCTTACTCTTTGATTCCACG
SSRY169	ACAGCTCTAAAAACTGCAGCC	AACGTAGGCCCTAACTAACCC

). Amplifications were performed with a Bio-Rad C1000 Thermal Cycler (Bio-Rad), using Kapa 2G Robust HS kits (Roche, Cat#KK5515), with primer pairs purchased from Thermo Fisher Scientific. The forward primer of each pair was labelled on the 5′ end with a specific fluorophore. Genotyping reactions were carried out using Hi-Di Formamide (Life Technologies, Cat: 4311320) and GeneScan 600 Liz Size Standard (Applied Bio, Cat: 4366589), with a 3 min denaturation step at 95C followed by cooling on wet ice for 3 min. The products of these reactions were analyzed using a 3730xl Genetic Analyzer (Applied Bio) with the PowerPlex 4C matrix Standard dye set (Promega, Cat: DG4800) and POP-7 Polymer (Thermo Fisher Scientific, Cat: 4363929). The resulting raw data (.ab1) files were then imported into GeneMapper 5.0 for peak analysis.

We analyzed data using the R software package (version 4.0.3 [32]) with the geosphere [33], igraph [34], ade4 [35], poppr [36], pegas [37], and ape [38] add-on libraries.

We calculated geographic distances between collection sites (chacras, or worked fields) using two methods. Direct distances (ellipsoid distances) between sites were calculated from their latitudes and longitudes using geosphere. River distances, the distances between sites assuming travel was restricted to waterways, were calculated from a map of collection sites and their respective rivers using igraph. The pairwise direct and river distances were also tested for correlation using ade4’s mantel function.

Because our sampling strategy was aimed at identifying as many landraces as possible, it was not well suited for population genetic analysis. However, to develop a preliminary portrait of population genetic variation in our sample, we conducted two basic inquiries. First, we calculated locus-specific observed and expected heterozygosities, and performed tests for Hardy–Weinberg equilibrium using the exact test implemented in pegas. We then obtained three measures of heterozygosity within landraces using Coulon’s GENHET R program: PH_t (the number of heterozygous loci/number of genotyped loci), H_sobs (PH_t/mean observed heterozygosity across loci), and H_sexp (PHt/mean expected heterozygosity across loci) [39].

We calculated genetic distances between landraces using poppr’s bruvo.dist function, which estimates distances for STRs under a stepwise mutation model. We then used the resulting pairwise distance matrix in two further analyses. First, we used it to ascertain the neighbor joining tree relating landraces. We did this using ape’s nj function with 10,000 bootstrap replicates. We also used the matrix of pairwise distances to determine whether correlations between genetic distance and geographic distance were present. These analyses were performed using ade4’s mantel function.

To further examine genetic similarities and differences among landraces, we used principal components analysis (PCA) and discriminant analysis of principal components (DAPC) [35,36]. PCA summarizes multidimensional differences between samples in a reduced number of dimensions, capturing those accounting for the most variance (in our case two), making them easier to visualize. These analyses were performed with ade4’s dudi.pca function. DAPC is similar to PCA. However, it identifies the reduced dimensions along which populations defined a priori are most distinct, as opposed to those explaining the most overall variance. These analyses were performed with poppr’s dapc function. We then used poppr’s amova function, which determines the contributions of population structure to genetic variance, to determine the statistical significance of differences between rivers and between landraces within rivers. In addition, to ascertain the extent of differentiation among rivers, we used ade4’s dudi.pca function to obtain three measures of genetic distance, Nei’s G_ST, Hedrick’s G′_ST, and Jost’s D, which define differentiation as a function of expected heterozygosities within and between populations [40,41,42].

Finally, we examined the extent to which landrace name is an indicator of relatedness by comparing our core sample of 43 landraces with the opportunistic sample composed of 17 specimens. These were the specimens with landrace names already present in the core sample, which were thus putative duplicates. Here, we obtained STR genotypes from each duplicate, then identified its closest relative across both our core sample and the other duplicates using poppr’s bruvo.dist function. This revealed whether its closest relatives had the same landrace name or a different one, and whether or not they were from the same river.

3. Results

3.1. River Distances and Direct Distances

Consistent with their definitions, the river distances between chacras in our sample were always greater than or equal to the direct distances (Figure 3). In some cases they were much greater, with the river distance being nearly 18 times the direct distance. These were cases in which chacras were near one another by direct travel, but travel restricted to waterways entailed descending one lengthy river and ascending a second one nearby. For instance, sites on the upper Tahuayo were much farther from sites on the upper Itaya via river travel than by direct travel (Figure 1b). The minimum distance between collection sites was 0 km with respect to both direct distance and river distance, which occurred when different landraces were collected from the same chacra. The maximum distances between chacras were 167 km via direct travel and 338 km via river. Mantel tests revealed that correlations between direct and river distance were strong and statistically significant (Mantel’s r = 0.91; β = 1.7; p < 0.001).

3.2. Diversity in STRs

We found polymorphism at all 13 genotyped loci (

Table 2. Genetic diversity at individual loci. Alleles is the number of STR alleles at the given locus, Length Range gives the length of the longest and shortest STR amplification products. Obs. Hz. is observed heterozygosity, Exp. Hz. is expected heterozygosity, HWE p-value is the p-value in tests for Hardy–Weinberg equilibrium.

Marker	Alleles	Length Range	Obs. Hz.	Exp. Hz.	HWE p-Value
GA12	15	123–153	0.45	0.75	<0.001
GA13	4	133–139	0.00	0.14	<0.001
GA21	7	104–118	0.37	0.74	<0.001
GA126	10	180–214	0.66	0.79	<0.001
GA127	10	205–237	0.40	0.82	<0.001
GA131	8	99–115	0.50	0.83	<0.001
GA136	8	135–159	0.26	0.76	<0.001
GA161	10	90–128	0.84	0.80	<0.001
SSRY32	3	296–300	0.02	0.35	<0.001
SSRY46	3	264–268	0.00	0.18	<0.001
SSRY70	5	245–253	0.15	0.73	<0.001
SSRY83	2	239–241	0.00	0.04	0.013
SSRY169	2	100–102	0.00	0.17	<0.001
Mean	6.69		0.28	0.55	<0.001

). The number of alleles observed at each ranged from two (at SSRY83 and SSRY169) to 15 (GA12), with a mean of 6.69. The length range of alleles also varied across loci, ranging from two (SSRY169) to 38 (GA161), with a mean of 16.46. These patterns are consistent with prior studies of the same markers in yuca [3,30,31,43,44,45,46,47]. Observed heterozygosities at individual loci ranged from 0.00 to 0.66, with a mean of 0.28, and expected heterozygosities ranged from 0.04 to 0.83, with a mean of 0.55. Tests for Hardy–Weinberg equilibrium revealed a highly significant deficit of heterozygosity at all loci, with p < 0.001 for all but one locus, SSRY83, which was significant with p = 0.013. Within landraces, the fraction of heterozygous loci ranged from 0.00 to 0.54, with a mean of 0.29 (Table 2). Standardized within-landrace heterozygosities (H_sobs) ranged from 0.00 to 1.92 with a mean of 0.52, and standardized within individual expected heterozygosities (H_sexp) ranged from 0.0 to 0.98 with a mean of 0.52 (

Table 3. STR diversity within landraces.

Landrace	PHt	Obs. Hs	Exp. Hs	IR	HL
amarilla	0.33	1.39	0.66	0.34	0.48
andioca	0.31	1.10	0.56	0.37	0.57
añera	0.31	1.10	0.56	0.47	0.56
arpón	0.30	1.09	0.63	0.28	0.50
blanca	0.23	0.82	0.42	0.55	0.66
blanca (dark stem)	0.54	1.92	0.98	-0.06	0.23
blanca (light stem)	0.31	1.10	0.56	0.43	0.57
brava (cultivated)	0.31	1.10	0.56	0.53	0.57
brava (feral)	0.23	0.82	0.42	0.57	0.67
bufeo	0.42	1.43	0.78	0.21	0.39
cerveza	0.31	1.10	0.56	0.39	0.56
cogollo colorado	0.00	0.00	0.00	1.00	1.00
cogollo morado	0.15	0.55	0.28	0.68	0.77
colorada	0.23	0.82	0.42	0.55	0.66
gallinazo	0.00	0.00	0.00	1.00	1.00
iguano	0.38	1.37	0.70	0.30	0.45
indianino	0.31	1.10	0.56	0.41	0.55
inviernino	0.00	0.00	0.00	1.00	1.00
lobera	0.31	1.10	0.56	0.41	0.56
lobera colorada	0.46	1.64	0.84	0.12	0.33
lobera negra	0.42	1.37	0.72	0.23	0.44
lupuna	0.00	0.00	0.00	1.00	1.00
lupunillo	0.33	1.10	0.58	0.38	0.54
mano de tunche	0.15	0.55	0.28	0.77	0.78
morada	0.31	1.10	0.56	0.42	0.62
morada amarilla	0.50	1.83	0.94	0.00	0.26
motelillo	0.38	1.37	0.70	0.24	0.45
napino	0.31	1.10	0.56	0.39	0.55
palmera	0.23	0.82	0.42	0.56	0.66
palo blanco	0.31	1.10	0.56	0.42	0.57
palo negro	0.33	1.18	0.63	0.29	0.50
palo negro Antonio (narrow leaf)	0.36	1.18	0.64	0.31	0.49
palo negro Antonio (wide leaf)	0.36	1.14	0.64	0.42	0.49
palo negro Arimuya	0.17	0.59	0.31	0.62	0.75
piririca	0.38	1.37	0.70	0.23	0.44
posheco enano	0.31	1.10	0.56	0.42	0.57
posheco gigante	0.31	1.10	0.56	0.38	0.54
señorita	0.25	0.88	0.47	0.48	0.62
umishina	0.38	1.37	0.70	0.23	0.44
ungurahui	0.46	1.64	0.84	0.06	0.35
vidrio	0.23	0.82	0.42	0.55	0.66
viejillo	0.15	0.55	0.28	0.72	0.78
virginia	0.23	0.82	0.42	0.60	0.68
Minimum	0.00	0.00	0.00	−0.06	0.23
Maximum	0.54	1.92	0.98	1.00	1.00
Mean	0.29	1.01	0.52	0.45	0.59

). The highest and lowest values for the three measures were consistent across landraces, with the highest being observed in a single specimen of yuca blanca (dark stem variety) from the Rio Nanay and the lowest being observed in four landraces originating from the Rio Orosa (cogollo colorado, gallinazo, inviernino, and lupuna).

3.3. NJ Tree

Construction of the neighbor joining tree identified three clades supported across 100% of bootstrap replicates (Figure 4). We designated these A, B, and C. Each contained landraces from multiple rivers. Clade A contained six landraces, which were found on four of the five sampled rivers, the Nanay, Orosa, Pintuyacu, and Tahuayo. Clade B contained seven landraces, which were found on the Itaya, Pintuyacu, and Tahuayo. The remaining 30 landraces constituted clade C, which contained representatives from all sampled rivers. A possible exception to the lack of substructure occurred in clade C, which contained a subclade composed of eight landraces from the Orosa River. It was not well supported statistically, occurring in fewer than 50% of bootstrap replicates when constructing the overall tree. However, it contained two pairs of landraces supported by more than 50% of bootstrap replicates.

3.4. Population Structure

Mantel tests for correlation between genetic distance and geographic distance revealed statistically significant relationships (Figure 5). Genetic distance was correlated with both direct and river distance at a high level of significance, with Mantel’s r = 0.12 and p < 0.001 for river distance and r = 0.20 and p < 0.001 for direct distance, and the increase in genetic distance with geographic distance was β = 5.14 × 10⁻⁴ for direct distance, and β = 2.31 × 10⁻⁴ for river distance.

The first two principal components calculated in the PCA accounted for 11.9% and 8.7% of the variance, respectively (Figure 6a). On these axes, landraces from the Itaya, Nanay, Pintuyacu, and Nanay rivers showed no evidence of differentiation. Landraces from the Orosa River showed a different pattern. Consistent with the neighbor-joining tree and distribution of closest relatives across rivers, some landraces from the Orosa were more similar to landraces from other rivers than to other Orosa landraces. However, several Orosa landraces showed evidence of being more similar to each other than to landraces from other rivers. In addition, the broad distribution of Orosa landraces across the first two variance components of the PCA was broader than that of landraces from the other rivers, suggesting the Orosa harbors more genetic diversity.

Discriminant analysis of principal components, which identifies the axes along which populations (rivers in our study) are most distinct, clarified the patterns present in the PCA analysis (Figure 6b). It revealed that landraces from the Orosa River were separated from the other rivers, consistent with the neighbor joining tree, in which the Orosa showed some evidence of being a clade. Two of the other rivers, the Nanay and Pintuyacu, overlapped extensively, and together were somewhat differentiated from all other rivers. Similarly, the Itaya and Tahuayo rivers overlapped and were mostly distinct from the other rivers.

The AMOVA analysis revealed that levels of differentiation between rivers and between landraces within rivers were 3.05% and 38.84%, respectively. Their associated p-values were 0.016 and 0.001, thus departing significantly from expectations (

Table 4. AMOVA test of diversity within and between landraces and rivers.

Variance Source	d.f.	Sum of Squares	Variance Component	Percent Variation	p-Value
Between rivers	4	42.28	0.17	3.05	0.016
Between landraces within rivers	39	299.16	2.19	38.84	0.001
Within landraces	44	144.43	3.28	58.11	0.001
Total	87	485.88	5.65	100.00

). The AMOVA also revealed that levels of diversity within landraces were significantly greater than expected (p = 0.001), and made the largest contribution to genetic variance overall, 58.11%.

Estimates of population differentiation across the five rivers for individual loci ranged from −0.017 to 0.119 for G_ST, −0.101 to 0.496 for G′_ST, and −0.079 to 0.429 for D (

Table 5. Population differentiation by locus.

Marker	Nei’s GST	Hedrick’s G′ST	Jost’s D
GA12	0.017	0.079	0.059
GA13	0.045	0.063	0.007
GA21	0.080	0.307	0.231
GA126	−0.017	−0.101	−0.079
GA127	0.016	0.098	0.079
GA131	0.096	0.496	0.429
GA136	0.089	0.308	0.224
GA161	0.017	0.089	0.070
SSRY32	0.036	0.064	0.020
SSRY46	0.117	0.156	0.017
SSRY70	0.119	0.411	0.311
SSRY83	0.001	0.001	0.000
SSRY169	0.082	0.114	0.015
Global	0.052	0.128	0.068

). Values for combined loci were 0.052, 0.128, and 0.068 for G_ST, G′_ST, and D, respectively. Estimates of pairwise, between-river differentiation varied across measures (

Table 6. Population differentiation by river.

(a) Nei’s GST.
	Itaya	Nanay	Orosa	Pintuyacu
Nanay	0.046
Orosa	0.040	0.032
Pintuyacu	0.037	−0.004	0.035
Tahuayo	0.015	0.045	0.046	0.040
(b) Hedrick’s G’ST.
	Itaya	Nanay	Orosa	Pintuyacu
Nanay	0.168
Orosa	0.167	0.137
Pintuyacu	0.135	−0.014	0.150
Tahuayo	0.057	0.161	0.188	0.141
(c) Jost’s D.
	Itaya	Nanay	Orosa	Pintuyacu
Nanay	0.087
Orosa	0.097	0.080
Pintuyacu	0.067	−0.007	0.088
Tahuayo	0.027	0.082	0.111	0.070

). However, all three were lowest between the Nanay and Pintuyacu (G_ST = −0.004, G′_ST = −0.014, and D = −0.007), indicating that they were the least differentiated. The highest values were found between the Orosa and the Tahuayo (G_ST = 0.046, G′_ST = 0.188, and D = 0.111), indicating that they were the most differentiated. All other pairwise values were intermediate.

3.5. Landrace Names and Relatedness

Pairwise comparisons between the 60 specimens in our overall sample (the 43 specimens in our core sample and 17 putative duplicates) revealed that landrace name was a poor indicator of genetic relatedness (

Table 7. Closest relatives of named landraces observed at least twice in our sample and the tributary of origin of each.

Specimen	Closest Relative	River(s) of Origin
amarilla 1	piririca	Itaya
amarilla 2	arpón	Itaya
amarilla 3	palo blanco	Orosa
amarilla 4	palo negro Arimuya	Nanay
añera 1	lobera colorada	Itaya
añera 2	cogollo colorado	Pintuyacu
brava (cult.) 1	morada	Orosa
brava (cult.) 2	palo negro Antonio (n)	Nanay
indianino 1	inviernino	Orosa	Same river
indianino 2	palmera	Orosa	Both from same river
napino 1	palmera	Orosa
napino 2	brava (cult.) 1	Orosa
piririca 1	lupuna	Orosa
piririca 2	arpón	Itaya
señorita 1	arpón	Itaya
señorita 2	brava (cult.) 1	Nanay
umishina 1	vidrio	Tahuayo
umishina 2	palo negro Antonio (n)	Nanay
ungurahui 1	mano de tunche	Orosa
amarilla 5	lobera 1	Tahuayo, Orosa
amarilla 6	palo negro Antonio (n)	Itaya, Nanay
añera 3	morada	Itaya, Orosa
lobera 1	vidrio	Orosa, Tahuayo
lobera 2	amarilla 1	Orosa, Tahuayo	Each from different river
motelillo 1	morada	Tahuayo, Orosa	Different river
motelillo 2	motelillo 1	Itaya, Pintuyacu
motelillo 3	palo negro	Pintuyacu, Nanay
ungurahui 2	andioca	Nanay, Pintuyacu

). Of the 28 specimens with landrace names occurring more than once in our sample, all but one were most closely related to a landrace with a different name. For instance, of the six specimens named amarilla (one in our core sample and 5 putative duplicates), none were most closely related to another amarilla specimen. Depending on the individual, their closest relatives were piririca, arpón, palo blanco, palo negro Arimuya, lobera, and palo negro Antonio. Only one specimen was most closely related to a specimen with the same name. In that case, an example of motellilo collected on the Pintuyacu river was most closely related to a specimen named motelillo from the Itaya river.

4. Discussion

In a previous study, we identified 45 phenotypically distinct, traditionally named landraces across five Amazon tributaries in northeastern Perú [17]. We found that growers identified landraces using a range of traits and named them accordingly. We also observed distinctions among landraces even when they were planted in the same chacra. These observations suggested that genetic polymorphism is present, including divergence among yuca types. Our analyses of polymorphism at STR loci in this study confirmed the presence of genetic variation in yuca in the upper Amazon, and revealed previously unrecognized relationships between diversity, relatedness, traditional nomenclature, and geography.

4.1. Genetic Diversity

STR genotypes in our sample revealed extensive variation. Polymorphism was present at all 13 of the loci we examined, which combined to form 43 multilocus genotypes (MLGs) in the 43 landraces we collected. Thus, every landrace was unique. This pattern is similar to patterns observed in prior studies, which have consistently reported genetic distinctions between landraces. For instance, in an investigation of 596 wild, bitter, and sweet yuca landraces in the Brazilian Amazon, Alves-Pereira et al. [43] identified 351 MLGs defined by 14 STRs. Thus, the majority (64%) were unique, although many (36%) were identical to at least one other type. Diversity in yuca has also been reported on smaller spatial scales. In an isolated Wayāpi village in French Guiana, Duputié et al. [18] genotyped 10 STRs in 436 specimens, from 61 named landraces and found 383 MLGs. This pattern was noteworthy because it revealed that the average Wayāpi landrace harbored more than 6 MLGs. Thus, genetic diversity is found within named, Amazonian yuca landraces. Our sampling strategy, which focused on identifying different landraces, did not allow us to investigate the issue in our own sample.

Studies beyond the Americas have found genetic diversity and landrace differences as well. In an investigation of 11 villages in Uganda, Kizito et al. [46] genotyped 288 specimens representing 93 yuca varieties, and encountered patterns similar to those reported by Duputié et al. [18]. Different MLGs occurred within yuca varieties, and some MLGs were shared among them. However, on average, varieties were significantly differentiated. This was noteworthy because yuca was introduced to Africa from South America just 500 years ago, which might be expected to result in founder effects reducing diversity. Nonetheless, the two continents exhibited congruent patterns. Kizito et al.’s findings were reinforced by an analysis of eight loci in 522 bitter and sweet landraces from South America, Africa, and Oceania by Bradbury et al. [44]. Bradbury et al. found 188 MLGs, including 86 MLGs across 203 landraces in Africa. Thus, many landraces had identical MLGs. In sum, while studies of STRs in yuca performed to date are not fully comparable due to their differing sampling schemes, they uniformly characterize yuca as a genetically diverse crop that harbors polymorphism throughout its range, including diversity both within and between landraces.

We found additional evidence of differentiation between landraces in patterns of heterozygosity and the distribution of variation within and across landraces and rivers. First, our tests for Hardy–Weinberg equilibrium revealed a significant deficit of heterozygosity at every locus (p < 0.001 at 12 of 13 loci and p < 0.012 for the thirteenth) (Table 2). Departures from expectation in Hardy–Weinberg tests indicate that genotype frequencies are inconsistent with panmixia, with up- and downward departures from expectation indicating nonrandom mating. Reduced heterozygosity is a classic signature of population subdivision, with the subpopulations in our study being individual landraces [40]. It occurs because genetic drift proceeds more rapidly within subpopulations than across the population as a whole, resulting in disproportionate losses of diversity in them. Second, our AMOVA tests revealed that the fraction of overall genetic variance accounted for by between-landrace differences (38.84%) was significantly greater than expected given overall diversity in the sample. This is predicted under subdivision because STR allele frequencies in different subpopulations diverge due to genetic drift.

Our field observations of chacras and growers’ cultivation techniques were consistent with our inferences about landrace differentiation based on STRs, suggesting that landraces represent subdivided populations. We found that landraces shown to us by growers were always phenotypically distinct and were readily distinguishable (Figure 7). Phenotypic differences were evident even among landraces planted in the same chacra and thus could not be due to environmental effects. We further noted that growers are motivated to maintain separation for practical reasons. Different landraces are often preferred for different purposes, such as roasting, making casabe (yuca bread), brewing yuca beer (masato), and production of fariña (yuca meal). Admixture between landraces, if permitted, has the potential to produce plants with unpredictable phenotypes, reducing their usefulness. Growers’ tradition of propagation through cloning, rather than through cultivation of sexually produced seedings, which might be hybrids, ensures that barriers to admixture are high. Thus, observed patterns of phenotypic variation, genetic diversity, and cultivation practices all support the notion that landraces are distinct entities, and they are maintained by strong human influences.

Evidence that landraces are distinct raises questions about their origins. Our observations suggest that diversity in yuca is generated by sexual reproduction. We noted that although growers claimed to rely strictly on clonal reproduction rather than through propagation of seedlings, which might be undesirable hybrids, mature chacras usually contained some plants with flowers and fruits. Further, nearly all of the chacras we observed were interplanted, with two or more landraces growing side by side. Thus, sexual reproduction was occurring, and likely included hybridization between landraces. Seedlings produced this way were not fostered by growers. Their maturation lagged behind the main crop, so they would require extra attention, and they were treated as weeds. However, we speculate that seedlings are occasionally adopted when they show desirable traits, and may come into permanent cultivation. We also noted that abandoned chacras persist ferally, and sexual reproduction in them was likely unrestricted. This raises the possibility that growers augment or even fully establish new chacras with sexually produced feral plants.

In addition to revealing distinctions between landraces, genotypes in our sample revealed an unexpected pattern: an excess of within-landrace heterozygosity. Although heterozygosity across our sample as a whole was lower than expected, which was consistent with population subdivision, heterozygosity within landraces was higher than expected given the variation they did contain. Specifically, the proportion of loci that were heterozygous (H_sobs) was significantly greater (p < 0.001) than the proportion expected under random mating (H_sexp) in 39 of 43 of our sampled landraces. Moreover, the mean observed value was nearly twice the expectation. In addition, the results of the AMOVA indicated that the fraction of variance accounted for by within-landrace diversity, 58.11%, was the largest single contributor to overall genetic variance, and was highly significant statistically (p < 0.001). This suggested that heterozygous plants are, through some mechanism, favored in yuca chacras.

Excess heterozygosity in within landraces is consistent with the long-held hypothesis that heterozygotes have fitness advantages [48,49]. In the case of yuca, these could arise as the result of natural selection imposed by the environment, artificial selection imposed by growers, or both. Evidence of heterozygote advantage has been reported in genes mediating a range of bioprocesses, from fat storage to pesticide resistance [49]. A striking number participate in vertebrate immunity, a pattern attributed to heterozygote advantages in detecting pathogens [49]. Our STR data do not pinpoint specific genes influencing seedling or landrace survivorship. However, they do suggest that such genes are present because the STRs are unlikely to be under selective pressure themselves. Landraces’ high heterozygosity must be due to breeding regimes favoring hybridization. Given past findings, an evident possibility is that heterozygosity in genes encoding yuca’s immune system provides advantages in resisting pathogens. Another possibility is that heterozygotes exhibit phenotypes preferred by growers, such as enhanced productivity, and are thus more likely to persist as the result of human imposed pressure—artificial selection.

4.2. Relationships between Landraces

In a previous study, we found substantial phenotypic variation in our sample, even among landraces collected from the same chacra. However, we found little evidence of landrace clustering in principal components analyses, which can reflect the presence of divergent lineages. Nonetheless, we could not rule out the presence of population structure with respect to genetics, because it can occur even when phenotypic differences are not evident. This phenomenon, cryptic population structure, is common and well documented in natural populations [50].

In contrast to phenotypes, genotypes in our sample did show evidence of phylogenetic structure. However, it was weak. In the neighbor joining tree, three clades supported by 100% of bootstrap replicates were present, but the branches separating them were short (Figure 4). Thus, the differences between the three clades were consistent, but small. In addition, while the presence of additional subclades was supported with bootstrap values of 51% to 88%, the subclades always contained two lineages, as opposed to defining more inclusive groups. The tree topology also suggested that eight landraces originating from the Orosa river may form a subclade, but it was not well supported statistically (<50% of bootstrap replicates). The results of our principal components analysis were consistent with these findings (Figure 6). They revealed that the majority of landraces formed a single tight cluster, rather than separate clusters as would be expected if clades are highly differentiated. In addition, several landraces from the Orosa river showed evidence of divergence from the main cluster, suggesting that they may be separate, but they were dispersed and did not form a distinct group.

A key observation in our study was that landraces have traditional names, which often reflect their phenotypes. Moreover, some traditional names are common across our study region. For example, we observed chacras planted with amarilla at multiple locations. Señorita was similarly common. This suggested that landrace name might reflect genetic relatedness, with landraces having the same name being more closely related. If so, it would imply that landraces are not just maintained locally, they are maintained across broader geographic ranges. However, we found little evidence that landraces’ names reflected their relatedness. Across the 28 specimens with names occurring more than once in our sample, 27 were most closely related to a landrace with a different name (Table 7). For instance, of the six amarilla specimens we collected, the closest relative of each was a landrace with a different name. Our first amarilla exemplar was most closely related to a landrace named piririca, the second was most closely related to arpón, and so on. Across the entirety of our sample, the sole exception to the pattern was a specimen named motelillo collected on the Itaya river, which was most similar to a specimen named motelillo from the Pintuyacu river, more than 200 km away via river and 60 km direct distance. Thus, in our sample, landrace name was a poor indicator of genetic similarity. The explanation for this finding is not clear. Resolving it will require an investigation of growers’ naming conventions, such as whether they incorporate phenotypic information, whether they change over time, and whether names are retained when landraces are transported from location to location.

4.3. Geography and Population Structure

The geographical distribution of populations is frequently a factor shaping their genetic relationships [51]. For instance, populations often exhibit isolation by distance, in which populations near each other are more similar genetically than are populations far apart. Populations can exhibit phylogeographic structure as well, with geographic structure shaping not just genetic distance but also phylogenetic relationships. An underlying factor in both of these phenomena is migration. Isolation by distance and phylogeographic structure can only emerge if migration rates between populations are sufficiently low to allow differentiation to occur. This raises questions about the relationships between genetic variation and geography in our yuca sample. On the one hand, the presence of distinct landraces and the geographic distance between chacras, which exceeded 150 km direct distance and 325 km river distance in some cases, suggested that regional differentiation was likely. On the other, the region’s human residents travel regularly, which may result in yuca migration. Residents often transport yuca from place to place in the process of moving their households or trading, for instance (Figure 8). If sufficiently common, this might have prevented the evolution of regional differences.

The associations we found between genetic distance, physical distance, and river of origin supported the presence of population structure with respect to geography. However, the structuring was weak (Figure 5). The presence of isolation by distance was strongly supported by Mantel tests, which revealed highly significant (p < 0.001) associations between genetic distance and both measures of geographic distance. These had slopes of β = 5.14 × 10⁻⁴ and β = 2.31 × 10⁻⁴ for direct and river distance, respectively, which translated to 10% increases in genetic distance over ~190 km and ~425 km. However, their correlation coefficients (Mantel’s r) were low, 0.12 and 0.20. Thus, while they were significantly associated, geographic and genetic distance were poor indicators of one another. The weakness of population structure with respect to geography was also reflected in the results of the AMOVA. In that analysis, differences between rivers were the lowest contributor to overall variance (3.05%). Differences between landraces on the same river (38.84%) and within landraces (58.11%) were higher. These patterns suggested that the yuca population in the region is not panmictic. However, the barriers to regional admixture are low. This situation can emerge through a variety of processes. For instance, it can occur following recent, rapid population dispersals, when subpopulations have had little time to differentiate. It can also occur when migration rates between subpopulations are high but not unlimited. Our data do not indicate which factor is responsible. Yuca’s domestication within the last several thousand years suggests that recent, rapid dispersal did occur. Further, the long-term frequency of human travel by both foot and canoe suggests that regional migration rates may have been high. The increasing accessibility of motorboats suggests that recent migration rates may be particularly high today. These possibilities are not mutually exclusive. Thus, we conclude that migration rates have been high, but we have little information about the relative importance of factors affecting it.

The results of the DAPC shed further light on the relationship between geography and population structure. It supported the results of the Mantel and AMOVA tests, suggesting that that spatial structure is present but weak. In the DAPC, the clusters representing the Pintuyacu and Nanay rivers overlapped extensively. This is consistent with the Pintuyacu being a tributary of the Nanay, such that they are expected to harbor related yuca landraces. Similarly, the cluster representing the Orosa river was distant from the clusters representing the other rivers, which is consistent with its location, more than 100 km away by river (Figure 1). However, patterns of clustering were not always consistent with geography. For instance, the Itaya river cluster overlapped extensively with that of the Tahuayo, but it did not overlap with the Nanay. This is inconsistent with geography because the Nanay and Itaya enter the main channel of the Amazon within 5 km of each other, and thus would be expected to share more variation than either does with the Tahuayo, which is more than 55 km away by river.

While they did not have significance values attached to them, the pairwise genetic distances we observed between rivers recapitulated the results of the DAPC, suggesting that geography and population structure are related, but weak. Some differences between rivers were consistent with their geography. For instance, the genetic distance between the Nanay and Pintuyacu rivers was ~0.0 for all three measures (G_ST, G′_ST, and D), which is consistent with the Pintuyacu being a direct tributary of the Nanay. This was reflected in the overlap between the Nanay and Pintuyacu in the DAPC. However, other distances were inconsistent with their geography. For instance, G_ST and G′_ST between the Itaya and the Tahuayo were lower than between the Itaya and the Nanay, despite the Itaya and Nanay’s proximity (<5 km) relative to the Tahuayo, 55 km away.

Our observations of human travel likely explain the weak relationships we found between genetic variation and geography at our study site. We found that people in our study region, like people throughout Amazonía, often travel. The use of canoes is ubiquitous among populations living near rivers and other water bodies, allowing travel with heavy loads such as passengers, personal belongings and goods, and harvests. River travel certainly occurs, including the transport of yuca landraces for personal use, sharing with family, and occasionally exchange with strangers. The clonal propagation of yuca makes this straightforward because it requires only stem cuttings. We regularly observed travelers carrying yuca cuttings to establish new chacras (Figure 8). Travel overland is also common. While rainforest terrain is stereotypically viewed as impenetrable, its inhabitants are familiar with regional topography and can travel for days on established trails. Overland travel does not allow the transport of the heavy loads possible with canoes. However, overland distances between locations can be substantially shorter than travel by river, and can be covered more quickly. In addition, yuca cuttings are light (usually stem pieces roughly 5 cm in diameter and 1.5 m long) and easily carried. Therefore, we hypothesize that both river and overland travel contribute to yuca’s migration, reducing the effects of geography on genetic differentiation, but not eliminating it. Taken together, these findings suggest that while husbandry practices maintain the distinctness of landraces, human travel obscures their geographic relationships.

5. Conclusions

Our findings provide a complex portrait of genetic variation in yuca in the Upper Peruvian Amazon. We found that the many morphologically defined landraces in the region, which are favored by growers for different purposes, do differ genetically. This is consistent with the prevalence of clonal propagation in the region, which is predicted to allow the perpetuation of distinct yuca varieties even when they are grown in close proximity, because it discourages hybridization. However, hybridization cannot be completely absent. We believe it most likely occurs through two types of rare event, the survival of sexually produced offspring in active chacras, and the adoption of feral, sexually produced plants. We also discovered that there is statistically significant spatial structure in yuca populations in our study region, but it is weak. Likewise, our phylogenetic analyses confirmed that differentiation among landraces is present, but also indicated that naming conventions are a poor surrogate for genetic information in predicting relatedness among landraces We hypothesize that these patterns are due in part to human movement in the region, which is common and includes the transport of landrace cuttings.

It is important to recognize that the cultural aspects of yuca production we observed do not necessarily apply to other populations growing the crop. Observations at other sites suggest that there may be profound differences between localities. This is emphasized by a comparison of our findings with those of Duputié et al. [18]. While the yuca cultivation we observed was carried out primarily by men, Duputié et al. observed that it was exclusively carried out by women. Moreover, while we observed that seedlings were excluded from propagation, Duputié et al. found that that growers frequently incorporated novel seedlings into their crops. The landrace names we observed also showed no overlap with those found by Duputié et al. Based on this comparison, and reports such as Boster et al. [11,19,20,52], we speculate that variation in both the cultural and biological aspects of diversity in yuca is extensive throughout the Amazon, and myriad landraces and accompanying traditions remain to be uncovered.