Genetic Diversity among Cowpea (Vigna unguiculata (L.) Walp.) Landraces Suggests Central Mozambique as an Important Hotspot of Variation

Agronomy 2020, 10, x; doi: FOR PEER REVIEW www.mdpi.com/journal/agronomy Article 1 Genetic diversity of cowpea (Vigna unguiculata L. 2 Walp) landraces suggests Central Mozambique as an 3 important hotspot of variation 4 Ana Maria Figueira Gomes 1,2,3, David Draper 4,5, Pedro Talhinhas 1, Paula Batista Santos 1, 5 Fernanda Simões 6, Nascimento Nhantumbo 2, Rafael Massinga 2, José C. Ramalho 1,7,*, Isabel 6 Marques 1,*, Ana I. Ribeiro-Barros 1,7,* 7 1 Plant Stress and Biodiversity Lab, Forest Research Center (AMFG, JCR, IM, AIRB) and Linking Landscape 8 Environment, Agriculture and Food (PBS, PT), Instituto Superior de Agronomia (ISA), Universidade de 9 Lisboa, Portugal. E-mail: [email protected] (AMFG), [email protected] (PT), 10 [email protected] (PBS), [email protected] or [email protected] (JCR), 11 [email protected] (IM), [email protected] (AIRB) 12 2 Divisão de Agricultura, Instituto Superior Politécnico de Manica (DivAG-ISPM), Mozambique. E-mail: 13 [email protected] (NN), rafael.m[email protected] (RM) 14 3 TropiKMan Doctorate Programme, Nova School of Business and Economics, Universidade NOVA de 15 Lisboa, Portugal. 16 4 National Museum of Natural History and Science and Centre for Ecology, Evolution and Environmental 17 Change. Universidade de Lisboa, Portugal. E-mail: [email protected] 18 5 UBC Botanical Garden & Centre for Plant Research, and Department of Botany, University of British 19 Columbia, Canada. 20 6 Unidade .de Investigação em Biotecnologia e Recursos Genéticos (UIBRG), Instituto Nacional de 21 Investigação Agrária e Veterinária, I.P. (INIAV), Oeiras, Portugal. E-mail: [email protected] 22 7 GeoBioSciences, GeoTechnologies and GeoEngineering (GeoBioTec), Faculdade de Ciências e Tecnologia, 23 Universidade NOVA de Lisboa, Monte de Caparica, Portugal. 24 * Correspondence: [email protected] or [email protected] (JCR), 25 [email protected] (IM), ar[email protected] (AIRB) 26 Received: date; Accepted: date; Published: date 27 Abstract: Cowpea is a multiple purpose drought-tolerant legume crop grown in several dry tropical 28 areas. Its domestication center is thought to be East or West Africa where a high level of genetic 29 diversity is apparently still found in many landraces. However, detailed genetic information is 30 lacking in many African countries limiting the success of breeding programs. In this work, we have 31 assessed the genetic variation and gene flow in 59 Vigna unguiculata (cowpea) landraces spanned 32 across six agro-ecological zones from Mozambique, based on nuclear microsatellite markers. The 33 results revealed the existence of high genetic diversity between the landraces, even in comparison 34 to other world regions. Four genetic groups were found, with no specific geographic pattern, 35 suggesting the presence of gene flow between landraces. In comparison, the two commercial 36 varieties had lower values of genetic diversity, although still close from the ones found in local 37 landraces. The high genetic diversity found in Mozambique sustains the importance of local 38 landraces and on farm protection in order to enhance genetic diversity in modern varieties of 39 cowpea worldwide. 40 Keywords: Africa; cowpea; genetic diversity; landraces; microsatellites 41 42 1. Introduction 43 Agronomy 2020, 10, x FOR PEER REVIEW 2 of 13 Cowpea (Vigna unguiculata L. Walp), also known as black eye pea, is a major annual grain 44 legume mostly grown in dry tropical areas of Latin America, South Asia and Africa [1]. It is cultivated 45 mainly for its grains, which have a high content of proteins (20-32%) and carbohydrates (50-60%). 46 Both grains and leaves, are also rich in the amino acids lysine and tryptophan, vitamin C, iron and 47 zinc [2]. Cowpea has therefore an essential role in the human diet in many developing countries being 48 referred as the “poor man’s meat” [3]. As a legume, it is also an important component of traditional 49 cropping systems since it fixes atmospheric nitrogen and contributes to soil fertility improvement 50 particularly in smallholder farming systems where little or no fertilizer is used [4]. The bulk of 51 cowpea production and consumption is sub-Saharan Africa, namely West and Central Africa [1], 52 where its nutritional value and tolerance to drought place this crop in an unique position to the 53 continent’s efforts to establish nutrition sensitive food systems that are more likely to help curb 54 malnutrition, particularly among the most vulnerable – pregnant or lactant women and children 55 under five [5]. Although cowpea is known to be drought tolerant when compared to other crops, the 56 productivity of cowpea varieties is hampered by erratic rainfall and many are sensible to heat [1]. 57 Thus, appropriate agronomic practices could improve the performance of new varieties, under 58 different agro-ecological zones. Indeed, physiological, and metabolic studies show a progressive 59 acclimation of cowpea plants to stress [6] and differential drought responses of landraces with 60 contrasting tolerance levels [7]. 61 Despite being native to Africa [8], the domestication center of cowpea is unclear but thought to 62 be either in East or West Africa where a high morphological and genetic diversity is found, followed 63 by a sub-domestication region in India [8-10]. European accessions usually cluster together with those 64 from West Africa and were likely imported from this region [10]. Breeding lines in America also show 65 a high genetic similarity with African accessions [11] although local American landraces show a high 66 genetic divergence [10]. In addition, regions like East Africa and Oceania show the lowest genetic 67 diversity suggesting the presence of bottlenecks or founder effects during cowpea migration to these 68 areas [10]. 69 Because of this domestication history linked to a center of origin in Africa, cowpea research has 70 been underway in several African countries for many years. Breeding activities in sub-Saharan Africa 71 involving germplasm collection, evaluation and screening for the identification of lines with high 72 yield potential resulted in a diverse cowpea germplasm collection constituted by more than 15000 73 cultivated cowpeas from 89 different countries [1]. Additionally, a core collection of more than 2000 74 accessions based on geographical, agronomical and botanical descriptors has been established in The 75 International Institute for Tropical Agriculture (IITA) genebank with the aim of discovering new 76 traits related with stress tolerance for the development of new breeding lines [12]. On the other hand, 77 cowpea has several features of a classical model plant for genomic studies, such as a relatively small 78 diploid (2n=2x=22 chromosomes) genome of ∼613Mbp, a short annual life-cycle and a highly selfing 79 nature [13]. 80 The limited number of cowpea breeding programs in Mozambique has contributed to the 81 country ineffectiveness in taking the advantage of the continent’s high genetic potential. A significant 82 pool of cowpea landraces is thought to be available, but the limited detailed information about their 83 diversity and agronomic potential makes it difficult for breeding programs to thrive. Thus, the 84 characterization of cowpea genetic resources available in Mozambique is of extreme importance for 85 conservation and breeding, since it is the second most cultivated legume crop in the country, 86 occupying an extension of ca. 380 000 ha, with an average yield of 0.275 t ha-1 [1]. Unlike commercial 87 varieties, landraces maintained by farmers usually have high levels of genetic variability as they have 88 evolved from years of uncontrolled cross-regional and infield genetic exchange, even between 89 previously released and discontinued open pollinated varieties [14], not being subjected to selection 90 over a long period of time. However, knowledge about their variability is usually limited [15]. 91 Therefore, the aim of this study was to assess the genetic diversity of cowpea landraces from five 92 agro-ecological regions across three provinces of Mozambique, using Single Sequence Repeat (SSR) 93 markers. 94 2. Materials and Methods 95 Agronomy 2020, 10, x FOR PEER REVIEW 3 of 13 2.1. Plant material 96 Fifty nine cowpea landraces corresponding to 10 populations were sampled in six agro97 ecological zones (AEZ) in the provinces of Manica, Sofala and Zambezia, where cowpea is grown as 98 an integral component of local cereal-legume cropping systems (Fig. 1): R3 (North and Central Gaza 99 and Western Inhambane), R4 (Medium altitude areas of Central Mozambique), R5 (Low altitude 100 areas of Sofala and Zambezia), R6 (Dry areas of Zambezia and Southern Tete), R7 (Mid-altitude areas 101 of Zambezia, Nampula, Tete, Niassa and Cabo Delgado) and R10 (High altitude areas of Zambezia, 102 Niassa, AngoniaMaravia and Manica). Additionally, two widely used commercial cultivars (IT16 103 and IT18) released by the Mozambican Institute of Agricultural Research (IIAM) and bred through a 104 partnership with the International Institute of Tropical Agriculture (IITA) in Nigeria were also used 105 in this study. 106 107 Figure 1. Left: Location of Mozambique in East Africa. Right: Studied landraces of Vigna unguiculata. 108 Population codes follow Table 2. Colors indicate the different eco-geographical zones (AEZs) of 109 Mozambique based on [16]. 110 2.2. DNA extraction and nSSR amplification 111 The 61 samples used in this study were genotyped based on nine polymorphic nuclear simple 112 sequence repeats (SSR’s) previously developed by [17]: VuUGM05, VuUGM22, VuUGM31, 113 VuUGM33, VuUGM39, VuUGM40, VuUGM68, VuUGM71 and VuUGM74. Based on an initial 114 survey, we selected these nSSR markers since they produced robust, highly polymorphic amplified 115 bands among the entire collection of cowpea samples. Total genomic DNA was extracted from 50 mg 116 of ground leaves using the InnuSPEED Plant DNA Kit (Analytik Jena Innuscreen GmbH, Germany) 117 according to the manufacturer’s protocol. The average yield and purity were assessed 118 spectrophotometrically by OD230, OD260 and OD280 readings (Nanodrop 2000, Thermo Fisher 119 Scientific, Waltham, MA, USA) and visualized by electrophoresis in 1% agarose gels under UV light. 120 Amplifications were performed in 15 μl reactions containing: 1.25U TaKaRa Hot startTaq 121 polymerase, 1X Buffer I, 1mM dNTPs, 5 μM Primer F and R and 100 ng DNA under the following 122 Agronomy 2020, 10, x FOR PEER REVIEW 4 of 13 PCR conditions: an initial denaturation at 95 ℃ for 5min followed by 35 cycles of denaturation at 65 123 °C (20 sec), annealing at 56 ℃ for 30 sec and a final extension at 60°C for 30min. Allele sizes were 124 determined using GeneMapper 3.2 (Applied Biosystems; UK). 125 2.3. Genetic diversity and population structure 126 For each nSSR locus and landrace, genetic diversity was assessed by calculating the total number 127 of alleles (Na), mean expected heterozygosity (He), mean observed heterozygosity (Ho), allelic richness 128 (AR), and inbreeding coefficient (FIS) using FSTAT 2.9.3.2 [18]). GenAlEx 6 software was used to 129 estimate the mean expected heterozygosity (He) and mean observed heterozygosity (Ho) for each 130 population, as well as the number of private alleles [19]. The selfing rate (s) was estimated as s = 2FIS/(1 131 + FIS) [20]. An analysis of variance was used to detect significant differences between sites for the 132 measured genetic values. Grids for all significant genetic parameters were generated in R and are 133 based on a grid with a cell size of 30 seconds (which corresponds to approximate 1 km in the study 134 area) applying a 1.5-degree circular neighbourhood diameter. The circular neighbourhood is used to 135 re-sample the genetic composition of a single sample to all surrounding grid cells, with a size of 30 136 seconds, within a diameter of 1.5 degree around its location. In this way, the genetic composition of 137 each sample is representative for the area within the defined buffer zone. 138 2.4. Population structure and differentiation 139 The Bayesian program STRUCTURE v.2.3.4 [21] was used to test whether any discrete genetic 140 structure exists among the landraces and regions sampled. The analysis was performed assuming a 141 number of clusters from K =1 to K = 8, with 10 repetitions per K. Models were run assuming ancestral 142 admixture and correlated allele frequencies with 50,000 burn-in steps, followed by run lengths of 143 300,000 interactions for each K. The optimum K was determined using STRUCTURE HARVESTER 144 [22], which identifies the optimal K based both on the posterior probability of the data for a given K 145 and the ∆K [23]. To correctly assess the membership proportions (q values) for clusters identified in 146 STRUCTURE, the results of the replicates at the best-fit K were post-processed using CLUMPP 1.1.2 147 [24]. POPULATION 1.2 [25] was used to calculate the Nei’s genetic distance [26] among individuals 148 and to construct an unrooted neighbour-joining tree with 1000 bootstrap replicates. A Principal 149 Component Analysis (PCoA) was also constructed in GenAlEx6 [27] to detect the genetic relatedness 150 among individuals based on Nei’s genetic distance. We estimated genetic differentiation among 151 locations using an analysis of molecular variance (AMOVA) with ARLEQUIN 3.11 [28]. Molecular 152 variance was quantified among populations and within populations considering AERs and wild 153 cowpea versus cultivars, using an AMOVA using 10,000 permutations at 0.95 significance levels in 154 ARLEQUIN 3.11 [28]. 155 2.5. Spatial analysis and genetic diversity rarefaction 156 Grids for genetic parameters were generated in DIVA-GIS (www.diva-gis.org), based on a grid 157 with a cell size of 2.5 minutes (which corresponds to approximatly 4.5 km in the study area) and 158 applying a circular neighborhood with a diameter buffer of one degree (corresponding to 159 approximate 111 km). The circular neighborhood was used to illustrate the allelic composition of each 160 sampled site representative for the area within the defined buffer zone. Genetic diversity rarefaction 161 considered the spatial average of several population parameters such number of alleles (NA), 162 observed heterozygosity (Ho), inbreeding coefficient (FIS) and % selfing rate (s). 163 3. Results 164 3.1. Genetic diversity 165 The total number of alleles varied between 49 in VuUGM74 and 145 in VuUGM40 (Table 1). For 166 each locus, observed heterozygosity values (Ho) ranged from 0.014 in VuUGM74 to 1 in VuUGM40 167 and expected heterozygosity (He) ranged from 0.016 in VuUGM74 to 0.806 to VuUGM33. FIS values 168 Agronomy 2020, 10, x FOR PEER REVIEW 5 of 13 varied between -0.008 and 0.857 (respectively for loci VuUGM68 and VuUGM31; Table 1) across the 169 loci studied. 170 Table 1. Characteristics and genetic diversity statistics of the nuclear microsatellite (nSSR) primers 171 used in the genetic study of Vigna unguiculata. For each locus, the total number of alleles (Na), mean 172 expected heterozygosity (He), mean observed heterozygosity (Ho), and the fixation index (FIS) 173 obtained from the 61 studied samples are shown. 174 Primer name Primer sequence 5′-3′ Gene Bank ID Na Ho He FIS VuUGM33 F: AAAGGTGGGGGATTATGAGG FG853417 83 0.907 0.806 -0.091 R: TGTCCAATCCTGATGGATGA VuUGM71 F: TTCACAACCTGTCCACCTCA FG819327 125 0.143 0.548 0.783 R: GGCGTCCCAACAGATAAGAA VuUGM05 F: GCGGGATTCTATTCCAGTGA FC459955 82 0.174 0.617 0.767 R: TCCATTGGGTTTCTCAACCT VuUGM39 F: CGAAAAAGCATGATCAACCA FG863845 97 0.149 0.749 0.851 R: CCCCTTTCGCTAAAATTTCC VuUGM22 F: CAATCACCATTCACCAAACA FG908248 112 0.181 0.629 0.749 R: TATTGGGACTCAGGTCTTGG VuUGM31 F: TGGTTCACTTCCCATATTGTC FG932695 122 0.136 0.711 0.857 R: AGGCAGAGACGAAGGAGTGA VuUGM40 F: TTCTACATGGTTTTGGGGTCA FG864565 145 1.003 0.671 -0.426 R: GAGCTTGCCCTCAAGAATTG VuUGM68 F: TGATTGATGGTGGTGTAGCC FG807949 59 0.415 0.397 -0.008 R: GCACTTCACTCATCGTTGCT VuUGM74 F: GCCTCCTCTCACAAACTTGC FF547768 49 0.014 0.016 0.018 175 A total of 327 alleles were found among the set of V. unguiculata landraces, varying significantly 176 between sites (P<0.001; Table 2). The number of alleles varied geographically from 14 in the coastal 177 area of Muchela to 71 in the dry western area of Tambara (Fig. 2). Allelic richness varied between 178 1.250 in Muchela and 1.751 in Gurué with no statistical differences being found between areas 179 (P=0.452; Table 1). However, the number of private alleles varied significantly across areas (P<0.001; 180 Table 2) with the highest number being found in Gurué, Tambara and Machaze (Fig. 3). 181 Table 2. Genetic diversity within the cowpea genotypes studied. The number of samples analysed 182 (N), total number of alleles (NA), mean allelic richness (AR), mean observed heterozygosity (Ho) and 183 expected heterozygosity (He), inbreeding coefficient (FIS) and % selfing rate (s) are shown for each 184 population. 185 Populations Province AEZ N NA AR Ho He FIS s Gurué (GUR) North Zambezia R10 6 46 1.751 0.389 0.688 0.506 60% Namarroi (NAM) North Zambezia R7 4 23 1.534 0.379 0.454 0.250 25% Muchela (MUC) Central Zambezia R7 4 14 1.250 0.222 0.535 -0.412 74% Lucas Branco (LUC) South Zambezia R7 4 22 1.432 0.426 0.577 -0.032 41% Nhamatanda (NHA) Central Sofala R4 4 31 1.682 0.278 0.479 0.707 59% Maringué (MAR) Central Sofala R5 3 22 1.503 0.407 0.494 0.310 30% Tambara (TAM) North Manica R6 23 71 1.612 0.320 0.654 0.592 68% Sede nova (SED) North Manica R6 3 23 1.562 0.222 0.451 0.323 69% Matsinho (MAT) Central Manica R4 3 23 1.577 0.221 0.451 0.156 67% Machaze (MAC) South Manica R3 5 29 1.555 0.267 0.500 0.549 64% IT-16 Commercial cultivar R4 1 12 1.333 0.333 0.167 - Agronomy 2020, 10, x FOR PEER REVIEW 6 of 13 IT-18 Commercial cultivar R6 1 11 1.222 0.222 0.111 - 186 187 Figure 2. Map of the number of alleles (A) and observed heterozygosity (B) in 30 seconds (1km) grid 188 cells applying a 1-degree circular neighborhood. Dashed lines indicate the agro-ecological zones [16]. 189 190 Figure 3. Population structure of Vigna unguiculata based on 9 SSRs and using the best assignment 191 result retrieved by STRUCTURE (K = 4). Each individual sample is represented by a thin vertical line 192 divided into K coloured segments that represent the individual’s estimated membership fractions in 193 K clusters. Landraces and province are indicated below. AEZs are indicated in individual labels with 194 different colours for better visualization. The two cultivars are also indicated. 195 The mean observed heterozygosity varied significantly between 0.222 (Muchela, Sede Nova and 196 Matsinho) and 0.426 (Lucas Blanco) (P<0.001; Fig. 2), and the mean expected heterozygosity varied 197 between 0.451 (Matsinho) and 0.654 (Tambara) without statistical differences (P=0.481; Table 2). FIS 198 A B Agronomy 2020, 10, x FOR PEER REVIEW 7 of 13 values varied significantly between sites (P<0.001; Table 2), ranging from negative values of -0.412 in 199 the coastal area of Muchela to positive values of 0.707 in the central area of Nhamatanda (Fig. 4). The 200 rate of self-fertilization in V. unguiculata also varied significantly between sites (P<0.001; Table 2) with 201 the lowest values found in the northern region of Namarroi (25%) and the highest in the coastal area 202 of Muchela (74%) (Fig. 4). 203 The two cultivars had a low number of alleles (IT-16: 11 and IT-18: 2) and allelic richness (IT-16: 204 1.333 and IT-18: 1.222) constrained by the small sampling size. However, although the observed 205 heterozygosity (IT-16: 0.333 and IT-18: 0.222) was higher than the expected one in both cultivars (IT206 16: 0.167 and IT-18: 0.111; P<0.001 in both cases), it was also lower than the ones found in most local 207 landraces (Table 2; P<0.001). 208 209 Figure 4. Map of the fixation index (A) and selfing rate (B) in 30 seconds (1km) grid cells applying a 210 1-degree circular neighborhood. Dashed lines indicate the agro-ecological zones [16]. 211 3.2. Genetic structure of V. unguiculata 212 The Bayesian clustering program STRUCTURE found the highest LnP(D) and ΔK values for K = 213 4 (Fig. S1). Results showed a high degree of admixture between populations without any specific 214 geographic pattern or clustering considering the different AEZs (Fig. 4). One cluster was 215 predominant and grouped all landraces from North Zambezia, and most landraces from Sofala and 216 Central Manica; the second cluster characterized Central and South Zambezia landraces; the third 217 clustered landraces from North Manica as well as Central Sofala; the fourth cluster was exclusively 218 composed by landraces from South Manica (Fig. 2). The two cultivars clustered with one the 219 predominant group found in several populations, although both cultivars showed signs of admixture 220 with the other clusters. 221 In accordance with these results, the NJ tree separated all groups assigned by STRUCTURE 222 revealing again no general correlation with the geographical distribution of landraces (Fig. 5). All 223 individuals from R3 and R7 were clustered into two different clades, one with 65% and the other with 224 34% bootstrap support (BS) value (Fig. 5). Most individuals from R6 clustered in the same group (57% 225 BS) while R4, R5 and R10 were clustered into two different groups. The two cultivars were nested 226 within the wild populations, although in two different separated groups. 227 A B Agronomy 2020, 10, x FOR PEER REVIEW 8 of 13 The PCoA spatially separated the landraces analysed into three main groups (Fig. 6). In 228 accordance to the NJ tree, the landraces from R3 and R7 were separated from the main group: the 229 first being on the up-left of axis 2 that accumulated 21.14% of variance while the second on the down230 left of axis 2. All remaining landraces were clustered in a heterogeneous group containing also the 231 two cultivars. 232 233 Figure 5. Population structure of Vigna unguiculata based on 9 SSRs and using the best assignment 234 result retrieved by STRUCTURE (K = 4). Each individual sample is represented by a thin vertical line 235 divided into K coloured segments that represent the individual’s estimated membership fractions in 236 K clusters. Landraces and province are indicated below. AEZs are indicated in individual labels with 237 different colours for better visualization. The two cultivars are also indicated. 238 239 Figure 6. Unrooted neighbour-joining tree of the studied Vigna unguiculata landraces including the 240 two cultivars, based on Nei’s Da genetic distance. Numbers associated with branches indicate 241 bootstrap values (BS) based on 1000 replications. Only BS above 30 are shown. Colours of branches 242 Agronomy 2020, 10, x FOR PEER REVIEW 9 of 13 indicate the four genetic groups found in STRUCTURE. AEZs are indicated in branch labels with 243 different colours following Fig. 3. 244 3.3. Genetic differentiation between populations 245 Overall, genetic differentiation was significantly low (AMOVA FST = 0.199, P < 0.001). The 246 analysis performed over the landraces sampled indicated that only 19.92% of the genetic variation 247 was attributed among AEZs (Table 3). The highest molecular variance was found among genotypes 248 within landraces (47.39%), followed by the one found within genotypes (32.69%; P<0.001; Table 3). 249 Remarkably, a very low molecular variance was found between wild cowpea versus the cultivars 250 (0.12%) being most of the variance found among individuals within samples (65.58%; Table 3). 251 Table 3. Analysis of molecular variance (AMOVA) for the sampled populations of Vigna 252 unguiculata. 253 Source of variance d.f. Sum of squares % of variance Among landraces Among AEZs 6 77.612 19.92 Among genotypes within landraces 54 207.109 47.39 Within genotypes 61 60.001 32.69 Among cowpea landraces vs. cultivars Among samples 1 4.772 0.12 Among individuals within samples 58 279.949 65.58 Within individuals 61 60.000 34.30 4. Discussion 254 Landraces harbor a genepool of unexplored alleles that constitute an unique set of genetic 255 resources for breeding to improve productivity, nutritional value, adaptation and resilience to 256 climate change [29-32]. Given their evolutionary history and adaptation to local conditions, landraces 257 usually have higher genetic diversity and environmental resilience than modern varieties [33-36]. 258 However, such richness tends to be lost because most of the current intensive agricultural systems is 259 based on few high-input and high-yielding cultivars [37]. Thus, a comprehensive characterization of 260 landraces towards the development of conservation and breeding strategies, is among the main clues 261 to face the major agricultural challenges related to population growth and environmental risks. 262 Despite the ongoing agricultural changes in Africa, according to our data, the nine 263 microsatellites employed in this study were highly polymorphic and revealed the existence of high 264 genetic diversity between landraces of V. unguiculata landraces from Mozambique (Table 1). A total 265 of 327 alleles were found among the 59 cowpea landraces, which can be attributed to high genetic 266 heterogeneity (Table 2). Indeed, the genetic diversity values found within the studied landraces (Ho: 267 0.2220.426; He: 0.4510.654) were much higher than the ones reported for cultivated cowpeas. For 268 instance, high-density single nucleotide polymorphism (SNP) genotyping using the Cowpea iSelect 269 Consortium Array studied population structure and genetic diversity in a set of 91 worldwide 270 cowpea accessions and found an average PIC and He of 0.25 and 0.31, respectively [8]. Similar results 271 were obtained by Huynh et al. [10] and Xiong et al. [9] using respectively, 422 cowpea landraces and 272 768 cowpea genotypes, collected in 56 countries. 273 In comparison, the two commercial cultivars (IT-16 and IT-18) had a very low number of alleles 274 and heterozygosity values, and cluster analyses (PcoA or NJ tree) showed no clear differentiation 275 between these modern varieties and landraces. Pairwise genetic distances reported in other studies 276 have also shown that African landraces were close to wild cowpea samples [10]. This suggests that 277 genetic diversity of these two commercial varieties is still close from the ones found in landraces 278 although more individuals are needed to accurately determine if genetic erosion is occurring. 279