East Timor as an important source of cashew (Anacardium occidentale L.) genetic diversity

Estimating relationships using genetic distances

Relationships between populations were estimated according to the methods described by Monteiro et al. (2016), specifically, the chord genetic distances introduced by Cavalli-Sforza and Edwards (DC, Cavalli-Sforza & Edwards, 1967) using the INA method computed in FreeNA (DC^INA), and Nei’s D distance (Nei, 1972) was calculated in GenALEx 6.5. The unweighted pair group method with arithmetic mean (UPGMA) and neighbor-joining (NJ) trees were produced using the ape v3.4. package (Paradis, Claude & Strimmer, 2004) for R v4.1.0 (R Core Team, 2021) based on 10,000 bootstrap values, assessed by a boot function from the poppr v2.1.0. package (Kamvar, Tabima & Grunwald, 2014). Trees were further edited in FigTree v1.4.2 (Rambaut, 2014). Genetic distance was used to measure the relationship among populations using two separate groupings: the first included the cashew populations from all three countries (East Timor, Indonesia, and Mozambique), and the second excluded the continental populations from Mozambique.

Analysis of molecular variance (AMOVA)

An analysis of molecular variance (AMOVA, Weir & Cockerham, 1984; Hill, 1996) was done using ARLEQUIN v3.5.1.3 (Excoffier & Lischer, 2010) to assess the hierarchical distribution of the genetic variation of the cashew populations analysed. Significance was assessed after 1,000 permutations. Two three-level AMOVAs were performed: one using the cashew populations in each of the three countries (MZ =2; IND =1; ET =11) as individual groups, and the second considering only the cashew populations in East Timor and Indonesia as groups. In each AMOVA, the total variance was partitioned into components to account for pairwise differences between groups (V_a, (1) MZ vs ET vs IND; (2) ET vs IND populations), differences among populations within those groups (V_b), and differences among individuals within the populations (V_c). Variance components (V_a, V_b, and V_c) were used to calculate the fixation indices (F-statistics; F_CT, F_SC, F_ST) according to Weir & Cockerham (1984).

Individual-based clustering analyses

Genetically distinct clusters were identified using two different methodologies, as previously described in Monteiro et al. (2016): a Bayesian clustering analysis using STRUCTURE (Pritchard, Stephens & Donnelly, 2000) and a multivariate analysis method called DAPC (Jombart & Balloux, 2009). These two different analyses were done because STRUCTURE uses allele frequency and LD information from the dataset directly assuming the HWE equilibrium, while DAPC does not consider a particular population genetics model outlining the genetic differentiation between and within groups (Jombart & Balloux, 2009).

Individual-based clustering analyses were performed with two datasets: one including East Timor, Indonesia, and Mozambique, and the second including just East Timor and Indonesia. The Bayesian model-based clustering algorithm in STRUCTURE v.2.3.4 was used to identify genetic clusters assuming admixture and correlated allele frequencies without using population information. In the first approach, including East Timor, Indonesia, and Mozambique, analyses were set for a burn-in period length to 100,000 followed by 1,000,000 MCMC iterations with K-values set from 1 to 14 with 10 runs computed for each K. For the second approach, including just the cashew populations in East Timor and Indonesia, the same settings were followed by configuring K-values from 1 to 12 with 10 runs in each K. Structure Harvester v0.6.94 (Earl & Bridgett, 2012) was then used to calculate ΔK ad hoc statistics from Evanno, Regnaut & Goudet (2005), which was then used to estimate the most likely K-value. CLUMPP v1.1.2 (Jakobsson & Rosenberg, 2007) was used to calculate the average replicate runs for the selected K-value to prevent problems with multimodality and label switching between iterations of STRUCTURE runs. CLUMPP results were then plotted using DISTRUCT v1.1 (Rosenberg, 2004).

The DAPC was performed in R (R Core Team, 2021) using the adegenet v1.3.1 package (Jombart, 2008) and the relative frequency of the alleles as the dataset, since presence/absence data may overlook relevant patterns in allele frequency. The find.clusters function was used to find the ideal K-value based on Bayesian information criterion (BIC) scores, maintaining default parameters and retaining all principal components (PCs). Cross validation was done using the xvalDapc function to determine the optimal number of PCs to retain in the Discriminant Analysis (DA). The generated DAPC script is available at Figshare (Barros et al., 2022).

Using genetic distance to estimate relationships between populations

UPGMA and NJ trees were built using Nei’s D and DC^INA (FreeNA) genetic distances across accessions screened for East Timor, Indonesia, and Mozambique (Table S3), and an analysis narrowed to East Timor and Indonesia (Table S4) was also performed. Similar tree topology structures were observed with both Nei’s D (Figs. S1 and S2) and DC^INA (Figs. 2 and 3) matrices, indicating a reliable topology regardless of the algorithms used for genetic distance. Because of the similarity between the two algorithmic results, only the trees derived from DC^INA distance matrices are presented in Fig. 2 (East Timor, Indonesia, and Mozambique analyses) and Fig. 3 (East Timor and Indonesia). In the UPGMA- (Fig. 2A) and NJ- (Fig. 2B) derived trees, three clusters were identified: one cluster (I) that included the ETTR populations (ETTR1-3) from Baucau, Viqueque (ETV), Cova Lima (ETSU), and the ETK population from Manatuto, grouped with Indonesia (IND); a second cluster (II) that included the remaining East Timor populations from Manatuto (ETNA), all populations from the Bobonaro district (ETMA, ETBAT, ETSAN), and the Manufahi population (ETFA); and the third (III) cluster including just the populations from Mozambique (MZB and MZD). These results support two different genetic clusters within East Timor: one including only East Timor populations (Cluster III) and the other with a high genetic similarity with the Indonesian cashew population (Cluster II). In the NJ tree (Fig. 2B), the dendrogram derived from the DC^INA distance matrix also presented three different clusters: cluster one included the cashew populations from Mozambique, as observed in the UPGMA (Fig. 2A); the second cluster included the Bacau, Viqueque (ETV), Cova Lima (ETSU), and Indonesian populations; and the third cluster included the Bobonaro populations (ETMA, ETBAT and ETSAN) and the remaining East Timor populations (ETFA, ETNA and ETK).

When observing the cashew populations without including those from Mozambique (Fig. 3), the two clusters were similar to clusters I and II observed in the UPGMA, as shown in Fig. 2.

Conversely, the two genetic distance algorithms (Nei’s D distance, Fig. 3B, Figs. S2A and S2B) produced dissimilar NJ-generated trees to those derived from both the three-country dataset (Fig. 2A, Fig. S1A) and the East Timor/Indonesia dataset (Fig. 3A, Fig. S2B) from the UPGMA, indicating complex population structuring.

Analysis of molecular variance

AMOVA results from all three countries (MZ, IND, ET) showed that molecular variation was mainly found within individual cashew trees (76%), whereas variation among populations and among individual trees explained 13% of the total genetic differentiation, in both cases (Table 5). AMOVA results including only East Timor and Indonesian cashew populations were similar, with genetic differentiation within individual trees (73%) also higher than among individual trees (14%) or among populations (13%).

Individual-based clustering using Bayesian and a multivariate discriminant analysis to uncover population structure

Two different approaches were used to identify individual clusters in East Timor, using Indonesian cashew populations as the outgroup: (1) covering all populations from East Timor, Indonesia, and Mozambique, and (2) excluding Mozambique.

In the first approach, STRUCTURE was run considering the highest range of clusters conceivable (K = 1–15). This analysis assigned K = 5 as the optimal number of groups based on the ΔK method (Fig. S3) outlined by Evanno, Regnaut & Goudet (2005). According to the results when including populations from East Timor, Indonesia, and Mozambique, in K = 5, Mozambique was grouped in a single cluster (orange cluster; Fig. 4A); Manatuto populations were grouped into two principal clusters (brown and red); Baucau populations were mostly in one cluster (blue), except for ETTR3 that included a mix of all clusters except the cluster from Mozambique; Cova Lima was grouped in a single cluster (green cluster); Bobonaro populations were grouped into the red cluster; and the Manufahi populations were mostly grouped into the red cluster sharing genetic flow with populations from Viqueque and Indonesia, which were grouped into two clusters (brown) (Fig. 4A).

The DAPC analysis was done without assigning an a priori group. For the first dataset, the clustering analysis determined that K = 10 was the one with the best combination of mean and 95% CI of BIC (Fig. S4). However, there was no “elbow” effect observed, but a large plateau in the number of clusters with a significant overlap among confidence intervals for the nearby number of clusters. K = 5 was identified as the optimal number of clusters for later analysis to produce results comparable to the STRUCTURE results (Fig. 4B). Using the xvalDAPC function, 40 PCs (highest successful assignment: 88.2%, with the lowest mean squared error: 0.178) were retained and four discriminant functions, conserving 87.5% of the variance (Fig. S5).

Cross validation was performed using the xvalDapc function with the outcome being the number of PCA axes retained against the proportion of successful outcome predictions. The results showed 40 PCA axes were retained (considering the highest successful assignment: 88%, with the lowest mean squared error: 17%) and two discriminant functions (explaining 87.5% of cumulative variance) for inferring the five genetic clusters (Fig. S6). The loading plots from both discriminant functions identify which variables (i.e., alleles/loci) contributed the most for the five genetic clusters. Both DAPC results (K = 9 and K = 10) identified the same variables responsible for cluster assemblage (172 (mAoR6), 328 (mAoR17), 174 (mAoR35), 170 (mAoR47)), which highlights the importance of these alleles for cluster discrimination (Figs. S6 and S9).

For the Bayesian analysis of the second approach, which included populations from East Timor and Indonesia, STRUCTURE was run considering the highest range of clusters conceivable (K = 1 − 13). This analysis identified K = 2 as the optimal number of groups (Fig. S8) based on the methods outlined by Evanno, Regnaut & Goudet (2005), with no alternative ideal K. Considering the high ΔK-values displayed for K = 2 (Fig. 5), two clusters were observed: the first cluster, in green, includes cashew populations in Baucau, Cova Lima, Manatuto, and Indonesia; and the second cluster, in red, includes cashew populations in Bobonaro, Manufahi, and Viqueque (Fig. 5A). For the DAPC analysis with K = 2 (Fig. 5B), the cross-validation retained 20 PCA axes (considering the highest successful assignment: 93.5%, with the lowest mean squared error: 10%) and two discriminant functions (explaining 96.4% of the cumulative variance) for inferring the two genetic clusters. For K = 2, there were differences between the STRUCTURE and DAPC analyses. The DAPC results showed a genetic diversity in populations from Baucau that was not shared with Indonesian populations, while the Bobonaro, Manufahi, and Viqueque populations showed a common allelic diversity with Indonesia (Fig. 5B), contrary to the STRUCTURE optimal clustering results (Fig. 5A). These analytical inconsistencies may be related to different pre-requisites associated with both methods: STRUCTURE assumes that markers are not linked and that populations are panmictic (Pritchard, Stephens & Donnelly, 2000), while DAPC is better for populations that are clonal or partially clonal. In this case, STRUCTURE provides a more realistic observation of the genetic diversity of cashew populations, given the panmictic nature and absence of clonal populations in cashew varieties.

Population structuring in the East Timor germplasm

The AMOVA showed that most of the genetic diversity lies within populations with little diversity present between populations within East Timor or between populations in each country. The most diversity was observed within accessions for both groupings (all three countries and only East Timor and Indonesia) suggesting a high gene flow between populations, which is in accordance with previous studies (Freitas & Paxton, 1996). Cashew is an allogamous species favoring cross-fertilization (Freitas & Paxton, 1996), allowing intraspecific hybridizations and heightening genetic variation. Overall, outcrossing species tend to have higher genetic variation within populations, whereas selfing species, or species with a mixed mating system, have less genetic variability (Nybom, 2004). Since cashew is an outcrossing species, negative to low inbreeding coefficients (F) were observed, which agrees with former studies (Freitas & Paxton, 1996; Layek et al., 2021).

A population structuring analysis revealed that the distribution of genetic diversity did not follow a clear geographic trend, despite a well-defined clustering observed between cashew populations in Mozambique with the remaining populations in East Timore and Indonesia. The unique clustering attributed to the cashew populations in Mozambique may be because the cashew populations chosen in Mozambique may be different varieties than the ones used by farmers in both Indonesia and East Timor. There may also be differences because of the geographical isolation of the island countries, Indonesia and East Timor, compared to continental Mozambique. When narrowing the analysis to East Timor, a complex population structuring was observed, linked to district-associated genetic diversity. Within the country, the Viqueque, Manufahi, and Bobonaro districts displayed unique allelic diversity, while Baucau, Cova Lima, and Manatuto had a high genetic similarity with the Indonesian cashew population. These dissimilar results when including vs excluding the cashew populations in Mozambique highlights the complex genetic diversity of cashew trees in the context of a continental (Mozambique) country where orchards are older and have incorporated improved cashew varieties. These results are in accordance with previous results in India (Archak et al., 2009), where cashew genetic diversity lies within geographical populations and where the sharing of allele frequencies among populations does not translate into an in-country population structuring.

The clustering of cashew populations in this study had an uncommon, yet existing relationship with geographical region under a district-wise distribution, which is contrary to previous reports in India (Archak et al., 2009), where no relationship with the geographic region was observed. One of the reasons for a district-specific genetic diversity pattern may be related to seed availability and farmer preference by district, with varieties selected based on local performance and yield. Considering that few cashew varieties have been introduced in East Timor, this genetic diversity distribution may be associated with an inter-exchange of seed material adapted to similar ecological conditions. Despite the high genetic diversity observed and attributed to the high heterozygosity, allogamous nature, and high gene flow found in cashew (Mitchell & Mori, 1987; Borges, 2018), in certain East Timor districts, such as Bobonaro, Viqueque, and Manufahi, the dissimilar genetic clustering from the remaining districts suggests a diverse genetic build-up, which could be attributed to different cashew varieties being planted.

The fixation index F (also called the inbreeding coefficient) exhibits values from −1 to +1. Values close to zero are expected under random mating, while substantial positive values indicate inbreeding or undetected null alleles (Monteiro et al., 2016). Negative values denote an excess of heterozygosity due to negative assortative mating or heterozygote selection. Overall, positive F-values were observed across all populations (Table 4), revealing that populations are at or near the Hardy-Weinberg equilibrium, further supported by the lower observed heterozygosity values under the HWE compared to the expected values (Table 4). Among the cashew populations collected (associated with different geographic regions), the inbreeding coefficient was lowest in the northern Bobonaro district (ETSAN, F =  − 0.04; ETMA, F = 0.14; ETBAT, F = 0.01, Table 4) and in the southern Viqueque district (ETV, F =  − 0.09), possibly because growers from both regions exchange seeds with the other regions. In contrast, the Cova Lima district (south) showed the highest inbreeding coefficient (F = 0.26). Both the Bobonaro and Viqueque districts displayed a significant share of planting material exchange with Bobonaro getting seeds from Indonesia and Viqueque getting seeds from Australia. In the studied populations, positive F values were observed, which may indicate some level of inbreeding. The difference between expected heterozygosity and observed heterozygosity might be due to evolutionary factors (such as inbreeding) or internal genetic factors (such as gene incompatibility). The introduction of new cashew seeds for new orchards might be another explanation; when establishing new cashew orchards, some producers used seeds from a single tree with good traits (high yield and large nuts). Moreover, according to Waples (2015), the genetic makeup of a given population can vary over time in response to evolutionary forces that affect the heterozygosity of the population relative to the Hardy-Weinberg equilibrium.

Considering the 11 populations studied from East Timor, population differentiation (F_ST = 0.129, Table 5) was relatively moderate in comparison to a previous study of cashew trees in Côte d’Ivoire (Kouakou et al., 2020) where a low differentiation was indicative of a common origin. However, moderate differentiation could also be associated with the lower number of molecular markers used in this study (12 SSRs vs 18 SSRs in Kouakou et al., 2020). In India, a similarly low genetic diversity among cashew trees was associated with cashew trees being a relatively recent introduction to the country (Archak et al., 2009), which is also reported in other countries (e.g., Benin, Kouami et al., 2020). In East Timor, there are two possible explanations for the cashew genetic diversity observed in this study: (i) multiple cashew varieties were introduced from the Bobonaro and Viqueque districts, these regions being regarded as the cashew introduction hotspots in East Timor; or (ii) new cashew varieties were introduced in the other districts in East Timor.

The moderate genetic variability observed among the populations screened could also be due to the relatively recent introduction of cashew trees into East Timor on an evolutionary timescale and the allogamous nature of cashew, resulting in high gene flow and exchange of genetic material. The unexplored, yet genetically diverse accessions from East Timor could be used in future cashew breeding programs to improve yield, quality, and other traits.

East Timor as an unexplored, yet important source of genetic cashew diversity

The wide distribution of cashew in its primary center of diversity, Brazil, has been attributed to mature fruit floating down water currents as well as bats dispersing cashew seeds (Johnson, 1973). However, outside Brazil, cashew distribution has been attributed primarily to anthropogenic efforts, rather than through natural means (Archak et al., 2009). The extent and distribution of genetic diversity, as revealed by the present study, provide some clues to cashew introduction and expansion in East Timor. The existence of substantial overall genetic diversity and two distinct genetic population groupings in East Timor point to multiple introductions comprising different founder populations. These results are in accordance with a recent report in Burkina Faso that suggests several historical cashew introductions with well-performing cashew cultivars disseminated through the same route across producing areas (Moumouni et al., 2022). If the introduction occurred as a one-time event, the founder effect would have been reflected with a low genetic diversity. Cashew expansion in East Timor has been promoted mainly through the introduction of different varieties from Brazil, Indonesia, Australia, as well as a so-called “native” cashew. According to information obtained from the National Directorate of Industrial Crops and Agribusiness (Ministry of Agriculture, Forestry and Fisheries (MAFF) of East Timor), the Portuguese introduced the cashew to East Timor (ET) in the 18th century, where it was planted ornamentally in various districts. Since the establishment of cashew as an industrial crop in the 1990s (Buss & Ferreira, 2010; Do Céu Guterres, 2017), there has been no standard implemented for which varieties to include in orchards. In the Cova Lima district, cashew orchards were implemented in the 1990s and few new cashew varieties have been introduced. In the Manatuto and Baucau districts, cashew orchards were planted with accessions known for higher yield and that were better-adapted to local agro-ecological conditions. During fieldwork in the Bobonaro and Viqueque districts, land farmers reported that several cashew accessions were being introduced, namely from Indonesia and Australia; the different genetic clusters in these districts may be associated with these introductions based on the preferences of the farmers, thus shaping the current East Timor cashew genetic diversity panorama.

Considering India’s importance as the place cashew was first introduced as a commercial crop in Asia (Singh, 2018), future studies on genetic cashew diversity in East Timor should include Indian cashew populations to help explain the historical introduction of cashew trees in the country. In East Timor, a relatively significant genetic diversity for an introduced species was observed, supporting the possibility of cashew being introduced repeatedly over time. Similar observations have been reported in India (Archak et al., 2009). Contrary to the results of this study, a significant level of redundancy (homogenous group) was reported within the Nigerian cashew germplasm (Aliyu, 2012), highlighting a narrow genetic diversity. Results gathered from East Timor are of major importance, since the country aims to invest in cashew nut exports. Assessing the genetic diversity of cashew trees in East Timor will help define a valorized genetic signature of cashew nuts from East Timor. The results of this study are, thus, also relevant to the local government of East Timor, since the country would compete directly with major world cashew nut exporters, such as India and Vietnam.

This study provides useful information on the genetic diversity of cashew populations in a largely understudied country, where cashew nuts are becoming an important export. The genetic diversity results from this study identify a genetic cashew diversity hotspot. Cashew has been implemented under a monoculture system and cashew orchards have been increasing to meet global market needs (Monteiro et al., 2017). With the rising potential of pests and diseases (Monteiro et al., 2015; Monteiro et al., 2017), the current narrow genetic diversity of cashew orchards is a major concern to the future sustainable production of cashews. To ensure a solid future, cashew varieties with improved agronomic traits, such as higher yield as well as biotic and abiotic stress tolerance should be introduced. Also, the increase of the genetic diversity of on-farm orchards at the short run and the identification of genetic resources for the development of cashew genetic management should not be neglected.