The aim of any breeding program is to generate high-yielding genotypes, and in this study an additional goal was to identify genotypes that were also high-N fixing. It is particularly useful to examine genotype performance at Yorito, where growing conditions are representative of the marginal production regions in Honduras. At Yorito, 14 genotypes had above-average yields coupled with above-average SNF performance (Figure 5). These included four Honduran conventional varieties (HON56, HON53, HON55, and HON77), four landraces (HON2, HON49, HON65, and HON66) and six PPB varieties (HON5, HON26, HON35, HON57, HON72, and HON78). Of the six PPB varieties, five were developed through participatory varietal selection, and one was developed through participatory plant breeding. These genotypes are dispersed throughout the phylogenetic tree (Figure 2A), suggesting a lack of close genetic relatedness; however, three of the high-yielding high-fixing PPB varieties share common genotypes in their pedigrees: HON5 and HON26 have a

common landrace parent, Cincuenteño (HON48); HON26 shares a conventional parent, Tio Canela 75 (HON55), with HON78; and, HON5 has Amadeus 77 as a parent, which is a daughter of Tio Canela 75. There was no apparent relationship between release date and higher yields, nor was there a temporal trend for SNF performance.

**Figure 5.** Genotype performance for nitrogen fixation (x axis) and yield (y axis) at Yorito. For this trial, average yield was 799 kg ha−1 and average nitrogen fixation was 43.3 %Ndfa, and these values are shown on the plot in dashed lines, dividing the plot into quadrants. Genotypes in the upper right quadrant are higher yielding and higher nitrogen-fixing.

#### **3. Discussion**

#### *3.1. Genotype Origins and Pedigree Explain Honduran Panel Structure*

The patterns observed in the phylogenetic tree and STRUCTURE diagrams derived from the SNP compositions of the lines in the Honduran panel largely agree with what is known about their geographic origins and their pedigrees, but there are also a few exceptions. The two large groupings in the dendrogram based on SNP similarities (groups I–III and IV) show a clear separation (with some exceptions) between landraces (CRI) and material that has been conventionally bred or is the product of participatory plant breeding (PPB). Among the landraces, the small red beans that came from Yoro, Francisco Morazán, and Intibucá, were randomly interspersed throughout group III of the phylogenetic tree with no particular pattern, based on genotype origin. In contrast, clustering of genotypes by region of origin is found among the black bean landraces (groups I and II). Group I consists of the black bean landraces, which came from Intibucá, Yoro, Francisco Morazán, and Lempira, and contains the Milpero landraces (HON36, HON47, HON44, and HON39) and 'Negro Opalaca' (HON70), which are all from Intibucá, and 'Negro' (HON43), which is found alone on the next branch is from the neighboring department, Lempira ['Merlot' (HON62) is also found in this region of the tree and is discussed below]. These landraces are the most distantly related genotypes with respect to the rest of the panel. The Milpero landraces were daylight sensitive when grown at Elora in 2014, and their photoperiod sensitivity was likely inherited from a common ancestor. All photoperiod sensitive varieties in the panel may carry the dominant *ppd* gene responsible for control of this trait [48]. The remaining black landraces

from Yoro, Francisco Morazán, and Intibucá are found without any particular pattern throughout the next branch of the tree, and two are found among the small red landrace branches (HON06 and HON30). The STRUCTURE analysis shows that most of the black genotypes are admixed, suggesting a closer genetic relationship to the conventional and PPB germplasm. No small red landraces were found within the black landrace branches; however, two small red landraces, 'Concha Rosada' (HON02) and 'Rosado' (HON22), are displaced and found in group IV among the PPBs. The STRUCTURE analysis shows that Concha Rosada contains only ~10% genetic material from the landrace subpopulation, and Rosado is almost equally admixed between the landrace and PPB subpopulations. In the case of Concha Rosada, this may be explained by the fact that this landrace is widely grown and has been used as a parent in participatory breeding efforts and thus shares ancestry with many PPB varieties. Alternatively, Concha Rosada may not be a traditional landrace but instead a creolized variety derived from a formal-sector variety introduced to the Yorito region in the early 1980s [23]. The reason that Rosado is found among the PPB varieties in the tree is less apparent. Rosado has not been used as a parent for any of the PPB varieties in the panel according to the pedigree information available. Rosado recently arrived in the Yorito region, and a survey of bean farmers revealed that its origin is unknown [29]. According to M. Gomez (FIPAH), the initial population of Rosado showed phenotypic heterogeneity, and some selection has been made to create a uniform line for PPB breeding in collaboration with Zamorano. At Elora 2014, Rosado had uneven maturity, which may indicate that the seed planted that year, and the seed grown for DNA extraction, was not a fixed homogeneous line, and this heterogeneity may have resulted in misplacement of this genotype in the phylogenetic tree.

The organization of the PPB branches of the phylogenetic tree can be explained in part by common ancestry. The Honduran conventional varieties, 'Tio Canela 75' (HON55) and 'Amadeus-77' (HON52), and the landraces, 'Estica' (HON11) and 'Vaina Rosada' (HON34), have been used frequently in generating the PPB genotypes in this panel. For example, five of the seven genotypes in the left-most PPB branch (HON19, HON18, HON15, HON16, and HON14) are derived from crosses with Tio Canela 75 or Estica, or both. The STRUCTURE plot indicates that these genotypes contain a greater proportion of conventional/PPB than landrace genetic material, and therefore, they are found within the PPB branches of the tree. However, six PPB varieties (HON10, HON05, HON12, HON25, HON33, and HON72) are found among the landrace branches of the tree, which indicates that they share greater genetic similarity to their landrace parent than the other genotypes in their pedigrees. For example, FPV 921-65 (HON33) has the landrace Vaina Rosada (HON34) and Amadeus 77 in its pedigree. FPV 921-65 is found in the same branch as Vaina Rosada among the landraces in the tree, and the STRUCTURE analysis shows that FPV 921-65 has more similarity to the landrace genetic subgroup. Marcelino (HON10) and Amilcar (HON05) also contain >50% landrace genetic material (Figure 2B) and are found among the landrace branches. Fourteen PPB genotypes have between 5% and 50% admixture with the landrace subgroup. Alleles favoring agronomic characteristics for local adaptation and culinary traits contributed by landrace parents were likely prioritized under selection among PPB progeny from crosses between landrace and conventional genotypes, resulting in PPB varieties consisting of a large proportion of landrace genetic material.

Common parentage may explain the clustering of the five Honduran conventional genotypes (Amadeus-77, CENTA San Andrés, Tio Canela 75, Aifi Wuriti, and Carrizalito) in the panel. CENTA San Andrés, Amadeus-77, and Carrizalito have common ancestry. Tio Canela 75 is a parent of breeding line EAP 9510-77, which was released in the early 2000s in countries across Central America, including in El Salvador, as CENTA San Andrés, and in Honduras, as Amadeus-77. Tio Canela 75 is also a parent of EAP 9510-1, a sister line of EAP 9510-77, which was released in Honduras as Carrizalito in 2003 and as Telire in Costa Rica in 2004 [49]. Accordingly, these genotypes are closely associated in the phylogenetic tree, although some genetic differentiation has occurred between CENTA San Andres and Amadeus-77. The line 'MD 30-75' (released as Tio Canela 75) was used as a parent in generating four of the five conventional genotypes' mentioned above. This is likely a result of the effort to introduce BGMV resistance into Central American germplasm, as MD 30-75 is a highly resistant line, which carries the

*bgm-1* resistance gene [50]. DEORHO (HON53) was not genotyped in our study; however, it also has MD 30-75 in its pedigree, and it is reasonable to suppose that it would also appear in this region of the dendrogram. Aifi Wuriti was the only black conventional genotype included in the panel. Interestingly, it is most closely related to three small red genotypes (Carrizalito, Conan 33, and Cedrón) rather than the small black genotypes, although Aifi Wuriti does not appear in the pedigrees of any of the small black PPB genotypes in the panel. The final Honduran conventional genotype, 'Dorado', is unusual because it is found among the landrace genotypes in the phylogenetic tree and the STRUCTURE analysis. This may be explained by the lack of a common genotype in its pedigree compared to the other conventional genotypes.

The placement of the North American check varieties (except OAC Rosito) as a separate group within the dendrogram is consistent with the unique alleles that they would be expected to have relative to the genotypes in this panel. The location of the check varieties in the landrace portion of the tree (groups I and II) may reflect the genetic diversity and wide-ranging geographic origins of the germplasm used in the University of Guelph/AAFC and Michigan State University/USDA-ARS bean breeding programs, which are aimed at introducing resistance to abiotic and biotic stresses and improving various agronomic traits. The locations of the North American checks interspersed throughout the admixed portion of the STRUCTURE plot also indicate that genetic diversity has been retained in these modern North American genotypes. OAC Rosito is a special case, since it was derived from El Salvadoran germplasm, and it is found among the PPBs in the phylogenetic analysis where it is most closely related to two PPB genotypes, Campechano (HON57) and Quebradeño (HON26). According to the STRUCTURE analysis, OAC Rosito is more genetically similar to the Honduran Conventional/PPB sub-group than the landraces. This suggests either that El Salvadoran landraces are not similar to Honduran landraces, or more likely that the landrace population from which OAC Rosito was developed was actually a creolized conventional variety. Germplasm provided by Zamorano has been used in El Salvadoran variety development since the early 2000s [26], and this could explain the genetic similarity of OAC Rosito to the PPB varieties in our panel, which were developed in collaboration with Zamorano.

#### *3.2. Optimizing Use of Genetic Diversity of Honduran Landraces and PPB Varieties*

The larger nucleotide diversity among the landraces (π = 3.20 × 10−4) than observed in the PPBs (π = 2.89 × 10−4) in the current study is consistent with general observations that landraces are more diverse than materials that are products of selection [51–54]. However, other studies that compared diversity in wild to domesticated bean accessions found wider diversification between those groups than we found between landraces and PPBs in our panel. Nanni et al. [55] reported that within the Mesoamerican gene pool, nucleotide diversity was 3.2 times higher among wild genotypes (π = 17.34 × 10−3) than domesticated genotypes (π = 5.43 × 10−3). The difference between landraces and PPBs in the current study was only 1.1 times, probably because of the small population size and because these genotype groups do not represent extremes of the diversity continuum that was sampled in the previous study. Both landraces and PPBs are selections from wild accessions, and the PPBs have probably not been as strongly selected as conventionally derived varieties.

The high level of diversity in Honduran landraces suggests they could be a source of novel alleles that could be used in breeding to improve various traits. Landraces are adapted to the environmental conditions of the locations where they were maintained, in some cases, over thousands of years. Landraces that were grown in fluctuating environments and in low-input agricultural systems may be enriched for rare alleles enabling phenotypic plasticity and inherent responsiveness to diverse abiotic and biotic stresses [56]. Landraces, in regions where they are still grown, have often been pushed to marginal production environments where their performance often exceeds that of modern cultivars [29,57,58].

Dry bean landrace germplasm across Mesoamerica is genetically diverse [59]. Soil conditions across this region are poor, and the terrain ranges from low to high altitude with steep slopes, leading to certain trait adaptations in the landraces. For example, 'Common Red Mexican', a red-seeded landrace from Mexico, has been found to be drought resistant [60], while 'Puebla 152', a black-seeded landrace also from Mexico, has superior SNF capacity [61]. Originating in the Andean region, G19833, a 'Chaucha Chuga' landrace from Perú, has tolerance to high levels of soil aluminum and low levels of phosphorus [56] and resistance to a number of bean diseases [62,63]. Our survey of the literature found genes in high landrace π regions associated with abiotic stress tolerance, phosphorus use efficiency, and nitrogen fixation (see Table S1). Conservation of landraces and mobilization of the unique genetic diversity they contain through plant breeding can help address the future need for higher yielding and climate resilient varieties.

#### *3.3. Regions of High Genetic Di*ff*erentiation Indicate Regions Impacted by Selection*

Genetic divergence between the PPB and landrace subpopulations in the panel is indicated by regions of high genetic differentiation (*F*ST), and these genomic regions may contain loci that have been subjected to selection pressure. We identified several regions on chromosomes Pv02, Pv07, Pv09, and Pv11 where *F*ST values exceeded 0.5. Similarly, in a study of genetic diversity of Italian bean landraces, Lioi et al. [64] reported that genomic regions related to domestication were concentrated on Pv02, Pv07, and Pv09 for Mesoamerican types. For comparisons between wild and domesticated bean landraces, Papa et al. [65] also reported significantly larger levels of *F*ST differentiation around genomic regions associated with domestication. While the genetic distance between the Honduran landraces and PPB genotypes included in our study is not likely as wide as that between the wild and landrace genotypes investigated by Papa et al. [65], similar trends towards genetic differentiation between landrace and PPB genotypes developed with modern breeding objectives and germplasm could be expected.

In particular, the genomic regions with large *F*ST differences may be associated with traits that were a focus of selection in PPB breeding. However, an extensive search of the recent bean literature did not reveal any known QTLs associated with agronomic, SNF, or WUE traits that are located within the regions of large *F*ST differentiation identified in this study. This may be because this is the first genomic survey study of Honduran material, and the distinguishing traits between landraces and PPB materials are specific to materials from that region or expressed in that location. In particular, because we are comparing two domesticated groups of genotypes, namely farmer traditional landraces and PPB varieties, the genes underlying the regions of large differentiation found in our study could be those responsible for local adaptation, culinary qualities, and favorable plant traits, rather than traits associated with domestication [5]. Additionally, the conventional germplasm used to generate PPB genotypes, either through crosses with landraces or through varietal selection, is largely limited to material in the Zamorano breeding program, which may have a specific genetic architecture.

In spite of the lack of previous QTL evidence for selection for domestication in the high *FST* regions, several genes in those regions that have been studied for various reasons may be associated with domestication. For example, disease resistance genes, such as those found on Pv02 for Rhizoctonia resistance (Phvul.002G323708, Phvul.002G323712; [40]), Halo blight resistance (Phvul002G326200; [41]), and Anthracnose resistance (Phvul.002G328300; [43]), and on Pv08 for Anthracnose resistance (Phvul.008G019600; [66]), and for the bean-rust interaction (Phvul.008G270500; [67]), have been associated with domestication in several crops [68,69]. Genes controlling agronomic traits have been identified in domestication studies in other crops [70–72]. Additionally, genes that affect survival in diverse growing conditions may have also been favored over the course of domestication. Two such genes are located in a region of high genetic diversity on Pv08; the ethylene-responsive transcription factor (Phvul.008G019600; [73]), an ortholog of an Arabidopsis gene known to be involved in flooding tolerance [74], and the transcription factor IIIA (Phvul.008G270400), which is upregulated in phosphorus-restricted conditions [75]. Soltani et al. [73] suggest that further studies are needed to understand the process of local adaptation and allelic selection using bean landraces and wild

populations. Insight to develop climate-resilient crops can be drawn from the study of crop adaptation under natural selection and domestication [76,77].

#### *3.4. Landraces are Superior Nitrogen Fixers*

Although some genotype- and environment-influenced variability was seen in our study, our examination of symbiotic nitrogen fixation in the Honduran panel revealed a wide range of capacity for this trait. The superiority of the landraces for SNF capacity at all trial locations may be the consequence of the continual selections of these materials under conditions of low soil fertility endemic to Central America. Even today, these materials continue to be grown by small scale farmers who do not have access to fertilizer inputs. Under these conditions, bean genotypes which have developed efficient associations with nitrogen-fixing bacteria would have a larger source of nitrogen for metabolic processes and better phenotypic fitness compared to poor nitrogen-fixing genotypes. Strong nitrogen-fixers would have a competitive advantage in the low input environments and would likely be preferentially selected over time. There may be parallels between the selection pressures during landrace evolution and the selection of heirloom bean varieties, which have also been shown to have strong SNF capacity [18].

There are few studies that have investigated SNF capacity of bean landrace genotypes. Heilig et al. [33] used 'Puebla 152', a black-seeded Mexican landrace known for its nitrogen fixing capacity [78], in a cross with conventional genotype 'Zorro' to create a RIL population to study SNF. The authors found that Puebla 152 fixed between 13.0 to 45.5% of the nitrogen in samples (seed + biomass), which was slightly more than Zorro, which fixed between 5.4 to 44.4% [33]. Many landraces in our study fixed more N than Puebla 152, indicating that Honduran landraces may be a useful source of SNF capacity. The SNF performance of Zorro ranged from 47.7 to 53.0%Ndfa in our study, a mid-range performance among our check genotypes, and overall better than its performance in Heilig's study [33].

The SNF performance of the progeny of the cross between the conventional genotype Zorro and the landrace Puebla 152 may be predictive of the performance of PPB varieties that are crosses between Honduran conventional varieties and landraces. The SNF performance of Honduran PPB varieties ranged between 20.5 to 55.7%Ndfa at Yorito. Although the focus of the participatory breeding program between FIPAH and Zamorano has been to generate higher-yielding genotypes, rather than on improving SNF performance, the SNF capacity seen for the PPB varieties falls within the upper range of that found for the RILs in Heilig's study [33].

For the Honduran panel, the insights gleaned from Yorito are of particular interest because this location has growing conditions representative of small-scale growers across the region, and as much as possible, local growing practices were employed in the trial. At Yorito, there was a range in SNF performance in the landraces, and overall genotypes belonging to this breeding history group performed better than the others studied. In addition, PPB genotypes derived from crosses with the best SNF-performing landraces had strong SNF capacity. Zamorano used the methods of participatory varietal selection to develop these PPB varieties with CIALs, enabling local growers to evaluate genotype performance on their farms. For example, Amilcar and San Jose were tested by various CIALs through the regional adaptation nursery (VIDAC, Vivero de Adaptación Centroamericano) in the mid-2000s. Amilcar has the high-SNF performing landrace Cincuenteño in its pedigree. Generally, native Rhizobia inhabit tropical soils and farmers do not use Rhizobia inoculants, although Zamorano disseminates SNF knowledge through the CIALs, including effective Rhizobia inoculants and protocols for use. In addition, the breeding program has the capacity to test SNF performance in '*bancales*' where soil N levels are low and SNF-related traits, such as nodulation, can be observed. PPB breeding for enhanced SNF capacity could be expanded if grower demand and the threat of climate change and resulting raising input costs warrant it.

In addition, the range in SNF performance among the conventionally bred North American checks and Honduran conventional varieties was wide. The superior SNF capacity Merlot exhibited at Yorito would suggest it has value as a breeding parent for this trait in Honduras; however, it has a larger seed size and a dark red seed coat; traits that are less preferred by Honduran consumers and could be challenging to select against in a breeding program. Of the Honduran conventional genotypes, DEORHO (HON53) fixed the highest amount of N (56.8%), but it performed poorly in Elora fixing 34.4% (2014) and 52.6% (2015) of its N. It is not found in the pedigrees of any PPB genotype included in our panel, but DEORHO has been a popular variety in commercial growing regions of the country. It has disease resistance, high yield, and the preferred light red seed coat color, making it a good candidate for future PPB breeding efforts.

#### *3.5. PPB Genotypes Have Superior Water Use E*ffi*ciency Values*

Plants that have higher water use efficiency (WUE) are more drought tolerant, and WUE can be estimated using carbon differentiation (Δ13C) values measured from plant biomass. During photosynthesis, plants discriminate against the incorporation of the heavy C isotope (13C), depleting 13C in plant biomass and driving lower Δ13C values [79]. Plants with comparatively low biomass Δ13C values can be considered more drought tolerant. WUE has been studied in beans, including landraces. A study by Munoz-Perea et al. [80] of the WUE of 16 dry bean genotypes in drought-stressed and nonstressed environments found that the two landraces included differed in their responses to drought stress, but Common Red Mexican was among the best performers under drought stress conditions. In contrast, in our study, the significantly lower Δ13C values measured in Yorito for the PPB genotypes than the landraces indicates that the PPB varieties in our panel may be more resilient to drought conditions than the landraces.

The drought resistant characteristics of the PPBs were likely contributed by the conventional parents. For example, PM2-Don Rey (HON23) was the most WUE PPB genotype at Yorito, with a Δ13C value of 16.36‰. PM2-Don Rey was developed through PPB methods from a cross between the landrace, Paraísito, and the Honduran conventional variety, Carrizalito. It was released as a drought-resistant variety in 2016 [81], and it is the first variety from the EAP-Zamorano-CIAL PPB collaborations to be released at the national level. A second PPB genotype, Marcelino (HON10), was developed through participatory varietal selection, and it had similar WUE (16.43‰) to PM2-Don Rey. The PPBs FPY-724-43 (HON16; 17.03‰), Cedrón (HON03; 17.41‰), and Amilcar (HON05; 17.42‰) have the next best WUEs.

Three landraces had WUE values below 17.5‰, including Concha Rosada (HON02; 16.85‰), Chapin Rojo (HON27; 16.86‰), and Chirineño (HON67; 17.30‰). Concha Rosada is of particular note because it is favored by poor farmers for its early maturity, which allows it to escape drought conditions late in the growing season [29,30]. Our study indicates that Concha Rosada not only has drought resistance through 'drought escape' but also has WUE characteristics that enable it to survive drought.

Among the conventionally bred genotypes, Carrizalito (HON77) and OAC Rosito (HON63) were the most water use efficient with low Δ13C values of 17.47‰ and 17.54‰ at Yorito. Carrizalito is a commercial Honduran variety, and it was used as a parent contributing disease resistance, agronomic, and likely WUE traits, to PM2-Don Rey. Among the check genotypes, OAC Rosito, recently developed at the University of Guelph from an El Salvadoran plant introduction, had the best WUE performance. This characteristic is likely retained from its domestication in Central America, and this enabled it to outperform the other check genotypes that have been developed for production in the Great Lakes region of North America.

In the coming decades, the effects of climate change are predicted to bring drier conditions to Honduras, and drought-resistant crops will help protect yields through periods of minimal rainfall. It has been proposed that WUE can be improved through selection and breeding. In alfalfa, evaluating genotypes for Δ13C and selecting for lower Δ13C values has been used to improve WUE in this important forage species [79]. The current results indicate that there is variation among the Honduran PPB and conventional bean germplasm in WUE traits, and selecting for lower Δ13C values could be applied to beans in Honduras to generate improved varieties that are more resilient to drought conditions.

#### *3.6. Conventional Genotypes Have Superior Yields*

Releasing varieties with higher yields is the objective of modern breeding programs, and our trial results suggest that improvements impacting yield have been made along the breeding history continuum from landraces to PPBs to conventional varieties. Considering the Honduran germplasm, the landraces were the lowest-yielding group at Yorito, followed by PPBs and the conventional varieties. This result is counter to the findings of early experiments performed by CIALs, where landraces out-yielded conventional materials [23]. However, our trial was conducted at the FIPAH office in Yorito, where soil fertility is less restrictive and the plot is flat, whereas the early CIAL trials were conducted in farmers' fields, which have low-fertility soils and sloped land; conditions for which conventional materials were not developed. The superior performance of the conventional materials in our trial at Yorito is consistent with the aim of modern breeding practices in generating higher yielding varieties. The North American checks were the highest yielding at Elora 2014, as they were bred for performance in this region, whereas the Honduran conventional and PPBs performed poorly at Elora 2014. The landraces also performed well at Elora 2014, and this may be attributed to phenotypic plasticity resulting from retention of useful nucleotide diversity enabling them to perform well in a new environment.

#### *3.7. Utility of Panel Genotypes for Breeding*

The different breeding and/or selection histories for the materials contributing to the phenotypic diversity present in the Honduran panel may provide opportunities for improving different traits in beans in the same way that a number of studies with different crops have found unique benefits from the use of landraces. In wheat, for example, cultivation of landraces in low-input systems has led to the conservation of traits that increase the duration of photosynthesis, which can lead to an increase in grain yield [82]. In a study comparing barley landrace and modern cultivar performances under stress conditions, the landraces were higher yielding and were less likely to fail outright [58]. The advantage of using landraces as parents in breeding programs has also been explored. In a study examining barley yields under drought conditions, progeny from crosses using landrace genotypes were found to be higher yielding than progeny from crosses without landraces in their pedigrees [56]. The authors concluded that breeding crops for vulnerable environments could be enhanced by identifying landrace alleles associated with yield performance and abiotic stress adaptation and employing these in breeding programs [56].

In the current study, landraces, which had superior symbiotic nitrogen fixation characteristics could be excellent sources of novel alleles for this trait. Similarly, PPB materials, which had superior WUE, and cultivated varieties, which had superior yields within their target environments, might be exploited, respectively, for these purposes. In general, all the germplasm types that were tested represent useful resources for breeding for important traits in the face of climate change and increasing production costs/demands.

The diversity for SNF capacity inherent in Honduran bean landraces, and their unique adaptation to the microclimates where they are grown, leads us to conclude that the inclusion of landrace germplasm in breeding for enhanced SNF would produce high fixing genotypes with growth and culinary characteristics already accepted by small-scale bean growers.

#### **4. Materials and Methods**

#### *4.1. Plant Material*

The Honduran Panel was assembled in 2014 at the University of Guelph in collaboration with agronomists at FIPAH. The initial panel contained 27 landraces, 26 PPB varieties, and 5 Honduran conventional checks provided by FIPAH, as well as 6 North American checks sourced from the University of Guelph bean breeding program.

The landraces consisted of traditional inbred varieties unimproved by modern plant breeding, which are grown by subsistence farmers in hillside communities. The PPB varieties were generated either through participatory varietal selection (PVS) or participatory plant breeding (PPB) through a collaboration between the bean breeding program at Zamorano and CIALs associated with FIPAH. The landraces and PPB varieties were sourced by M. Gomez (FIPAH) from six departments in west–central Honduras. The landraces were from Yoro, Intibucá, Francisco Morazán, and Lempira, and the PPB varieties were from Yoro, Francisco Morazán, Santa Bárbara, and Comayagua (Figure 1). Seed was either collected directly from farmers in their communities or sourced from central seed banks maintained by FIPAH and PRR. The five Honduran conventional checks were developed for production in lower to mid-altitude, hillside and valley commercial-production regions of the country. The six North American varieties consisted of Merlot and OAC Rosito as small red market class checks, Zorro as a black market class check, and three navy beans: OAC Mist, a high-nitrogen-fixing genotype, R99, a nonfixing mutant, and its parent line OAC Rico. All genotypes in the panel belong to race Mesoamerica [83].

In the first trial location (Elora 2014), 10 Honduran genotypes were found to be daylight sensitive and were not grown at the subsequent locations. Additionally, seed of 16 genotypes that exhibited uneven maturity in Elora 2014 were sent to Puerto Rico for seed increase over the winter of 2015. For the second trial location (Yorito, 2014–2015), 12 new genotypes (6 landraces, 5 improved, and 1 Honduran conventional check) were added to the panel. For the third trial location (Elora 2015), seed harvested from Elora 2014, from the Puerto Rican seed increase, and from Yorito were used, as available. An additional Honduran conventional variety was grown that year to fill in the experimental design. Overall, a total of 77 genotypes were tested in the Honduran panel, 50 genotypes of which were grown at all 3 locations. A summary of the genotypes included in the Honduran Panel, including trial year, market class, seed source, and pedigree information is provided in Tables 6–8 according to breeding history.

*Plants* **2020**, *9*, 1238

**Table 6.**

Genotypes of the North American check and Honduran

conventional

 breeding categories that were included in the HON panel tested at three field locations,

region

not available.

*Plants* **2020** , *9*, 1238

**Table 8.** *Cont*.

#### *Plants* **2020** , *9*, 1238

**Table 8.** *Cont*.

breeding lines are selected through generations of testing by CIALs; classified by FIPAH. NA

Information not available.

#### *Plants* **2020**, *9*, 1238

#### *4.2. Field Experimental Design and Maintenance*

Field trials were carried out at the University of Guelph Elora research station (ERS) in summer 2014 and 2015 and at Yorito, Honduras in the *Postrera* season (planted December, 2014).

#### 4.2.1. Elora

The Elora fields were selected based on low soil nitrogen levels as measured by preplant soil tests and by site crop rotation histories that indicated that no dry bean crops had been produced in those fields for the previous decade at a minimum. In 2014, nitrates (NO3 −) were low at 7.1 ppm, and ammonium (NH+4) levels were 3.2 ppm. In 2015, nitrates (NO3 −) were very low at 4.8 ppm, and ammonium (NH+4) levels were 4.5 ppm.

A square lattice design (8 × 8) with two replications was used for each trial. At the ERS in 2014, 135 seeds of each genotype were grown in 4-row plots (150 cm × 190 cm, 37.5 cm between rows) with approximately 6 cm between plants. At the ERS in 2015, 60 seeds of each genotype were grown in 4-row plots (150 cm × 90 cm, 37.5 cm between rows) with approximately 5 cm between plants within rows.

Clean seed of each genotype was coated with commercially available Nodulator (Becker-Underwood) *Rhizobium leguminosarum* bv *phaseoli* peat-based inoculant prior to planting. The day before planting, inoculant powder (1/4 teaspoon, approximately 0.4 g in 2014; 1/8 teaspoon, approximately 0.2 g in 2015) was added to each seed envelope and the contents were shaken to coat the seeds. Inoculated seed was stored at the ERS at 4 ◦C until planting to maintain inoculant viability. The entire contents of each envelope (coated seed + loose inoculant powder) was planted. Successful inoculation was confirmed each year by observing pink (active) nodules on a few plants chosen at random throughout the trial at flowering time.

At the ERS, plots were maintained with standard practices throughout the growing season, except no-nitrogen fertilizer was used. Preplant fertilizer (0-20-20) at a rate of 200 kg ha−1 was applied approximately one week prior to planting. Preplant herbicide was applied to control broadleaf weeds, and pesticides were applied as needed throughout the growing season at standard rates to control leafhoppers, anthracnose, and root rot (see details in Wilker et al. [18]). Plots were manually weeded once before canopy closure each year.

At Elora 2014, the harvest was staggered according to maturity. The plots were pulled by hand at maturity and threshed at the side of the field using a Wintersteiger plot combine (Wintersteiger AG, Upper Austria, Austria) with a Classic Seed-Gauge weighing system by Harvest-Master (Juniper Systems Inc., Logan, UT, USA), and plot seed weight and moisture content were recorded. Plots that did not reach reproductive/physiological maturity were not harvested. In 2015, plot harvest took place after all plots reached maturity with an SPC20 Almaco plot combine (ALMACO, Nevada, IO, USA), which automatically recorded moisture and weight (kg ha−1) for each plot at 13% moisture.

#### 4.2.2. Yorito

The Yorito trial site was chosen based on field uniformity, access to irrigation, and proximity to the FIPAH regional office. Soil NO3 − levels were 18 ppm ("moderate" to "high") at Yorito, and the field had been used for bean and maize production in previous seasons. The FIPAH agronomist, M. Gomez, indicated that the trial site conditions were representative of the bean production areas around Yorito.

Seed of each North American check variety were provided by the University of Guelph, and seed of all Honduran bean genotypes in the trial were sourced in Honduras. A square lattice design (8 × 8) with two replications was used for the trial. At Yorito, 100 seeds of each genotype were grown in 2-row plots (100 cm × 500 cm, 50 cm between rows) with approximately 30 cm between plants and 3 seeds sown per hole, as per the traditional planting system.

Inoculant for the trial, a mixture of *Rhizobium etli* (CIAT 632) and *R. tropici* (CIAT 899) strains, was provided by J.C. Rosas (EAP-Zamorano) and was applied according to a protocol provided by *PIF* at a rate of 500 g ha−1 [93]. Briefly, the peat-based inoculant powder was applied to slightly moistened seed to ensure it adhered well to the seeds, and once sown, the seeds were covered with soil to protect the inoculant from the sun. CIAT 899 has high symbiotic stability and efficient N-fixation characteristics [15].

Plots were maintained with standard practices throughout the growing season. A preplant fertilizer of 12-24-12 NPK was used at a rate of 64.81 kg ha−1 as a formulation without N was not available in Honduras. Carbendazim (DEROSAL, Bayer) at a rate of 400 mL ha−1 was used to protect the trial from angular leaf spot (*Pseudoscercospora griseola*), common bacterial blight (*Xanthomonas axonopodis* pv. *phaseoli*), and rust (*Uromyces appendiculatus*).

At Yorito, plots were harvested by hand over a number of days as each plot matured. Yield (kg ha−1) was measured at harvest, the two plots of each genotype were bulked, and a subsample of seed was sent to Guelph for further analyses.

#### *4.3. Phenotyping*

#### 4.3.1. Elora

As plots initiated the reproductive phase, DTF was recorded as the date when 50% of the plants in a plot had one flower open. Relative leaf chlorophyll content was measured twice during the growing season (when the mean number of plots had reached (1) the second trifoliate stage and (2) at 100% flowering) using an SPAD 502 Plus Chlorophyll Meter (Konica-Minolta). The meter was calibrated according to manufacturers' instructions each time the unit was powered-on (https://www.specmeters. com/assets/1/22/2900P\_SPAD\_502.pdf). The middle leaflet in the top-most, fully expanded trifoliate leaf was used for the measurements, and three plants were sampled at random per plot.

As plots reached physiological maturity, DTM was recorded as the date when plots were ready to harvest. Three plants were randomly sampled from mature plots, placed in large paper bags, and dried in a repurposed tobacco kiln (De Cloet Bulk Curing Systems, model TPG-360, Tillsonburg, ON, Canada) at 33 ◦C at the ERS for 24–48 h. Roots were cut from each plant and the above-ground biomass was weighed. Plants were then threshed using an indoor belt thresher (Agriculex SPT-1A, Guelph, ON, Canada), their seed collected, weighed, and counted. Hundred seed weights (HSW) were calculated.

#### 4.3.2. Yorito

As plots initiated the reproductive phase, DTF was recorded as the date when 50% of the plants in a plot had one flower open. Disease ratings and agronomic scores were collected throughout the growing season; however, statistical analyses revealed no significant differences between genotypes, and these traits are not further reported here. DTM was recorded as the date when 90% of the pods in a plot had changed color.

#### *4.4. Isotope Analysis*

Seed from each plot was processed and analyzed as detailed in Wilker et al. [18]. Briefly, seed was finely ground, precisely measured, and the isotope analyses were carried out using mass spectrometry at the Agriculture and Agri-food Canada (AAFC) gas chromatography mass spectrometry facility in Lethbridge, Alberta. Samples were analyzed for δ15N (‰) and δ13C (‰).

To calculate the percent nitrogen derived from the atmosphere (%Ndfa), the natural abundance method was used on seed samples in this study [94]. Seed N represents the total N accumulated by a plant over the course of the growing season, and seed N values are representative of whole-plant N values [21]. Additionally, processing seed samples is more efficient than shoot tissue.

The natural abundance method uses the following equation,

$$\%Ndfa = \frac{\delta^{15}N\_{refrecuce\ plant} - \delta^{15}N\_{fixing\ plant}}{\delta^{15}N\_{refrecuce\ plant} - B}$$

where δ15*Nref. plant* is the rate of δ15*N* in the reference genotype (R99), δ15*Nfixing plant* is the δ15*N* of the N-fixing bean genotype, and B is the average δ15*N* of beans grown in an environment where its entire N source is from fixation [95]. The B-value was obtained for this experiment as described by Farid (2015) [36]. Briefly, δ15*N* was measured and averaged for 20 bean genotypes from both the Andean and Middle American gene pools, which were grown in a growth room in N-free media. Normalized δ15*N* values were used for all genotypes, and an average of δ15*N* values for R99 were used in %Ndfa calculations.

To estimate water use efficiency, δ13C values obtained from GCMS seed analysis were used following the methods proposed by Farquhar et al. [47]. Because the current WUE discussion utilizes Δ13C values, the raw δ13C values were converted to Δ13C using the following equation:

$$
\Delta \mathcal{C} = \frac{\delta^{13} \mathcal{C}\_{air} - \delta^{13} \mathcal{C}\_{p\text{hr}\text{rt}}}{1 + \delta^{13} \mathcal{C}\_{p\text{laut}}}
$$

where δ13*Cair* is the current free atmospheric level of approximately −8‰ and δ13*Cplant* is calculated per plant seed sample using appropriate C isotope standards. For example, a plant with a δ13C value of −28.2‰ yields Δ13C = (−0.008 + 0.0282)/(1 − 0.0282) = 0.0202/0.9718 = 20.7‰.

#### *4.5. Genotyping*

To enable discovery of the genetic structure of the HON panel, 73 genotypes (see Table 1) were genotyped for single nucleotide polymorphisms (SNPs). DNA was extracted following the same methods outlined in Wilker et al. [18]. Briefly, plants were grown in a controlled environment (16 h photoperiod, 22 ◦C) at the University of Guelph, and young-leaf tissue samples were harvested, freeze dried, and the DNA extracted according to manufacturer's instructions using the NucleoSpin Plant II kit (Macherey-Nagel, Dueren, Germany). DNA of adequate quality was sent to the Genome Quebec Innovation Centre (McGill University, Montreal) for SNP genotyping using the Illumina Infinium iSelect Custom Genotyping BeadChip (BARCBEAN6K\_3) containing 5398 SNPs [96]. TASSEL was used to filter the SNP data (MAF > 0.01) to a set representing 72 individuals and containing 4314 polymorphic SNPs [97]. Missing data comprised less than 3% of the data, and these were subsequently imputed using Beagle v4.1 [98] as described by Torkamaneh and Belzile [99].

#### *4.6. Population Structure*

The population structure of the Honduran panel was determined using a number of methods, as follows. First, the population structure was estimated using variational Bayesian inference implemented in fastSTRUCTURE [100]. Five runs were performed for each number of populations (K) set from 1 to 10 using the 4.3K genome-wide SNPs identified in this study. A ChooseK analysis was conducted to determine the number of subpopulations that maximize the marginal likelihood. Then, a principal component analysis (PCA) was conducted in TASSEL and plotted using PCAshiny in R [101]. Finally, the evolutionary relationships among the genotypes of the panel were inferred using the Neighbor-Joining method with the genome-wide SNP data [102]. The taxa were clustered together using the bootstrap test (1000 replicates) [103]. The tree was drawn to scale, with branch lengths (next to the branches) in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the Maximum Composite Likelihood method and the units correspond to the number of base substitutions per site. Evolutionary analyses were conducted in MEGA7 [104,105].

#### *4.7. Genetic Diversity*

The levels of genetic diversity in the landrace and PPB breeding history categories of the HON panel were assessed using the 4.3K SNP dataset and VCFtools [106]. The π statistic provides an indication of polymorphism within a population as measured by nucleotide diversity, and Tajima's D (*D*) provides an indication of selection pressure [107,108]. Both π and *D* were measured in sliding windows of 1 Mb across the genome using the—window-pi and—TajimaD options in VCFtools [106], which resulted in an average of 6 SNPs per window. The pairwise π and *D* values were also calculated among different subpopulations. Genome-wide averages of π and *D* for each breeding history category were generated by taking the average across all windowed calculations. Landrace and PPB π values were compared across the genome, and regions where landrace π exceeded PPB π by more than 3 times were considered highly differentiated, and the regions that were at least 25,000 bp in length were considered significant. To investigate the level of differentiation between the landrace and PPB genotypes, the *F*ST statistic was computed. *F*ST was calculated using the—weir-fst-pop option in VCFtools in sliding windows of 100 bp across the genome [106]. Weighted *F*ST values range from 0 with no genetic differentiation to 1 where fixation of alleles has occurred. *F*ST values exceeding 0.5 were considered significant in our analysis.

#### *4.8. Candidate Gene Investigation*

A literature search was carried out to identify previously reported QTL and genes that colocalize with regions where π values were highly differentiated between the Landrace and PPB categories. JBrowse (https://legumeinfo.org/genomes/jbrowse/) was used to explore the bean genome around SNPs with significant weighted *F*ST values in order to identify candidate genes. A 100 kb region centered on each significant marker was searched.

#### *4.9. Statistical Analysis*

Analysis of variance (ANOVA) tests were performed on the phenotypic data collected from each environment and all environments combined, using the GLIMMIX procedure in SAS (version 9.4, SAS Institute, Cary, NC, USA, 2012). In the combined analysis, genotype, environment, and the genotype-by-environment interaction were considered fixed effects, while all other effects and their interactions were considered random. In the separate environment analyses, genotypes were considered fixed effects, while all other effects and the interaction effects were considered random. The Shapiro–Wilk test was performed on the residuals in the UNIVARIATE procedure to test their normality [109]. Random and independent distributions of the residuals were visually examined by plotting the studentized residuals against the predicted values. Data that generated outlier residuals were removed from the data set. Further, single degree of freedom contrasts were conducted in GLIMMIX between breeding history categories—landraces, PPBs, and conventional and check genotypes—contrasting each category to each of the others. Repeated measures of leaf chlorophyll content (SPAD) were taken, and separate ANOVA tests were used to compare SPAD values at each time point in a combined analysis and by environment. The least squared means (LSmeans) for each trait were computed using the LSMEANS statement in the GLIMMIX procedure for each genotype.

Using the LSmeans calculated above, the pair-wise Pearson's coefficients of correlation were computed for all traits in the CORR procedure in SAS. The PRINCOMP and PRINQUAL procedures were used in SAS to generate the principal component (PC) values, to estimate the proportion of variance accounted for by each PC, and to plot PC1 against PC2 to generate genotype x trait (GT) biplots to determine genotype and trait interactions in each environment [110].

#### **5. Conclusions**

The aim of our study was to evaluate a large set of Honduran landraces and varieties generated through participatory plant breeding, as well as check conventional genotypes, to ascertain their value in future breeding efforts. We used simple genomics and phenotyping to characterize the panel. Our genetic analyses found that the panel is divided into predominantly-landrace and predominantly-PPB groupings, with Honduran conventional genotypes sharing most similarity to the PPBs. Breeding history and pedigrees account for this division. The genetic diversity analysis revealed that landraces have retained a higher level of nucleotide diversity than PPB genotypes, which we attribute to selection pressure imposed by breeding for different production environments/objectives, and the use of a small number of conventional/elite parents in breeding efforts. The nucleotide diversity inherent in landraces can be used to increase the frequency of rare alleles in breeding programs. Beyond genetic characterization, it is important to classify germplasm for trait phenotypic diversity, which could be employed in breeding, because landraces that have evolved in adverse environments contribute adaptive traits to variety development.

Two traits that contribute to climate resiliency, nitrogen fixation capacity, and water use efficiency were evaluated in our study. Genotypes with good nitrogen fixation capacity are an asset for remote hillside growers who have limited funds and limited opportunity to purchase inputs because of poor market access. Landraces were shown to have superior SNF capacity and are already favored by hillside producers. Genotypes with enhanced water use efficiency will also be an asset to hillside growers in a future with drier and more changeable growing conditions. The PPB and conventional varieties in our study show promising characteristics for drought resilience. Further evaluation of PPB varieties under drought conditions is warranted.

Honduran bean production will continue to be carried out predominantly by small-scale hillside producers. The widely grown farmer landraces are locally adapted and accepted by consumers, and future breeding efforts should deliver varieties that maintain the inherent SNF capacity of landraces, while enhancing drought resilience and producing high yields. PPB methods employed in the breeding efforts between EAP-Zamorano and FIPAH- and PRR-supported CIALs have succeeded in generating promising PPB varieties. One variety that combines these characteristics is Amilcar (HON05). It is a small red PPB variety selected among germplasm provided by Zamorano for improvement by a CIAL in Yoro. Ultimately, Amilcar has been widely accepted because of its commercial value, culinary qualities, and disease resistance, but its climate resiliency traits and yield potential are of critical importance. Landrace characterization, conservation, and employment in breeding programs will bring continued benefits.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2223-7747/9/9/1238/s1, **Table S1**: Candidate genes identified in regions of the *P. vulgaris* (2.0) genome where high nucleotide diversity (π) was discovered in landrace genotypes compared to PPB genotypes [5,32,33,35,37–39,46,63,67,111–115]. **Table S2**: Candidate genes identified in regions of the *P. vulgaris* (2.0) genome where SNPs with significantly high weighted FST values (>0.5) were found [40,42–44,46,116,117]. **Table S3**: F-test of fixed and Pearson's χ2 test of random effects in the combined GLIMMIX mixed model analysis of the Honduran panel genotypes tested in multiple field locations in Ontario, Canada and Yorito, Honduras, 2014–2015. **Table S4**: F-test of fixed effect of genotype overall and by breeding history category, and the Pearson's χ2 test of random effects in the GLIMMIX analysis of 63 genotypes tested at Elora, Ontario, Canada, 2014. **Table S5**: F-test of fixed effect of genotype overall and by breeding history category, and the Pearson's χ2 test of random effects in the GLIMMIX analysis of 63 genotypes tested at Elora, Ontario, Canada, 2015. **Table S6**: F-test of fixed effect of genotype overall and by breeding history category, and the Pearson's χ2 test of random effects in the GLIMMIX analysis of 63 genotypes tested at Yorito, Honduras, 2014–2015. **Table S7**: F-test of fixed effects and Pearson's χ2 test of random effects in the GLIMMIX mixed-model analysis of leaf chlorophyll content (SPAD) at early vegetative (SPAD1) and reproductive (SPAD2) stages in the Honduran panel overall (Combined) and at Elora in 2014 and 2015. **Table S8**: Phenotypic (rp) correlations among %Ndfa and other traits estimated in the Honduran panel grown in three locations in 2014–2015. Values shown for traits within environments are significantly correlated at the 5% level. The number of genotypes analyzed in each correlation is presented in brackets for each significant correlation. **Figure S1**: Histograms of Nitrogen derived from the atmosphere (%) values of genotypes comprising the HON panel tested at three field locations from 2014–2015. A. Elora 2014, B. Elora 2015, and C. Yorito. Breeding history category averages with standard errors are presented, followed by individual genotype LSmeans with standard errors. North American check genotypes (CK; purple), Honduran conventional genotypes (CV; blue), landraces (LR; green), and PPB varieties (PB; red). Genotype numbers correspond to those listed in Tables 6–8. **Figure S2**: Histograms of carbon discrimination (Δ) values of genotypes comprising the HON panel tested at three field locations from 2014–2015. A. Elora 2014, B. Elora 2015, and C. Yorito. Breeding history category averages with standard errors are presented, followed by individual genotype LSmeans with standard errors. North American check genotypes (CK; purple), Honduran conventional genotypes (CV; blue), landraces (LR; green), and PPB varieties (PB; red). Genotype numbers correspond to those listed in Tables 6–8. **Figure S3**: Histograms of days to flowering of genotypes comprising the HON panel tested at two field locations from 2014–2015. A. Elora 2014, B. Elora 2015. Breeding history category averages with standard errors are presented, followed by individual genotype LSmeans with standard errors. North American check genotypes (CK; purple), Honduran conventional genotypes (CV; blue), landraces (LR; green), and PPB varieties (PB; red). Genotype numbers correspond to those listed in Tables 6–8. **Figure S4**: Histograms of days to maturity of genotypes comprising the HON panel tested at two field locations from 2014–2015. A. Elora 2014 and B. Elora 2015. Breeding history category averages with standard errors are presented, followed by individual genotype LSmeans with standard errors. North American check genotypes (CK; purple), Honduran conventional genotypes (CV; blue), landraces (LR; green), and PPB varieties (PB; red). Genotype numbers correspond to those listed in Tables 6–8. **Figure S5**: Histograms of yield (kg ha−1) values of genotypes comprising the HON panel tested at three field locations from 2014–2015. A. Elora 2014, B. Elora 2015, and C. Yorito. Breeding history category averages with standard errors are presented, followed by individual genotype LSmeans with standard errors. North American check genotypes (CK; purple), Honduran conventional genotypes (CV; blue), landraces (LR; green), and PPB varieties (PB; red). Genotype numbers correspond to those listed in Tables 6–8. **Data Repository files**: Phenotypic, SNP, and genetic diversity datasets. Available at https://doi.org/10.5683/SP2/DJB7VY

**Author Contributions:** Conceptualization, J.W., A.N., S.H. and M.G.C.; methodology, J.W., K.P.P. and A.N.; software, D.T. and M.E.; validation, J.W. and M.E.; formal analysis, J.W., A.N., M.E. and D.T.; investigation, J.W.; resources, J.W., K.P.P., A.N., M.G.C. and J.C.R.-S.; data curation, J.W. and D.T.; writing—original draft preparation, J.W.; writing—review and editing, K.P.P.; visualization, J.W. and D.T.; supervision, K.P.P.; project administration, J.W.; funding acquisition, K.P.P., and S.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ontario Bean Growers, Agriculture Agri-Food Canada through their Cluster Program, the Ontario Ministry for Food and Rural Affairs, the Ontario Ministry for Research and Innovation, and the Food from Thought Program at the University of Guelph funded by the Canada First Research Excellence Fund. SeedChange/USC Canada, Global Affairs Canada, International Development Research Centre, World Accord, Development Fund, Norway, USAID.

**Acknowledgments:** GCMS analyses were performed by Brett Hill at AAFC, Lethbridge and arranged by Frederic Marsolais, AAFC, London.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

## **DNA Fingerprinting and Species Identification Uncovers the Genetic Diversity of Katsouni Pea in the Greek Islands Amorgos and Schinoussa**

**Evangelia Stavridou 1, Georgios Lagiotis 1, Lefkothea Karapetsi 1, Maslin Osathanunkul 2,3 and Panagiotis Madesis 1,\***

Received: 21 March 2020; Accepted: 5 April 2020; Published: 9 April 2020

**Abstract:** Pea (*P. sativum* L.), one of the most important legume crops worldwide, has been traditionally cultivated in Lesser Cyclades since ancient times. The commonly known traditional pea cultivar, 'Katsouni', is endemic to the islands of Amorgos and Schinoussa and is of great local economic importance. Despite the widespread cultivation of 'Katsouni' in both islands, it is still unknown whether the current Schinoussa and Amorgos pea populations are distinct landraces, and if they have common evolutionary origin. To assist conservation and breeding of the pea crop, the genetic diversity and phylogenetic relationships of 39 pea samples from Amorgos and 86 from Schinoussa were studied using DNA barcoding and ISSR marker analyses. The results indicate that both populations are different landraces with distinct geographical distribution and are more closely related to *P. sativum* subsp. *elatius* than the *P. abyssinicum* and *P. fulvum* species. Further characterization of the 'Katsouni' landraces for functional polymorphisms regarding pathogen resistance, revealed susceptibility to the powdery mildew (*Erysiphe pisi* DC.). This work represents the first investigation on the genetic diversity and population structure of the 'Katsouni' cultivar. Exploiting the local genetic diversity of traditional landraces is fundamental for conservation practices and crop improvement through breeding strategies.

**Keywords:** pea landraces; Amorgos; Schinoussa; DNA Barcoding; ISSR genotyping; HRM analysis; powdery mildew

#### **1. Introduction**

Species in the economically important *Fabaceae* family have been a staple human food for millennia and their use is closely related to human evolution [1]. Legumes, such as Spanish vetchling (*Lathyrus clymenum* L.), lentils (*Lens culinaris* M.) and beans (*Phaseolus vulgaris* L.) are an important plant-based protein source, rich in mineral nutrients, complex starch and fibers, and contain health-promoting antioxidants, such as carotenoids and phenolic compounds [2–4]. The usage of leguminous crops in traditional crop rotation systems, reduces the need for synthetic nitrogen-based fertilizers by forming symbiotic relationships with nitrogen (N)-fixing soil bacteria [5]. Such management practices are of great ecological importance and have high potential for conservation agriculture, considering legumes are functional either as growing crop or as crop residue [6].

Pea (*P. sativum* L.) is among the most important legume crops, such as chickpea (*Cicer arietinum* L.), lentil and faba bean (*Vicia faba* L.), in temperate climates and has a wide geographical distribution,

with field pea being specifically adapted to a wide range of climates and altitudes. The *Pisum* species are of high commercial importance and are cultivated worldwide for dry and fresh consumption. According to the International Legume Database (ILDIS) and to the classification of Maxted and Ambrose (2001) [7], the *Pisum* genus includes three species: i) *Pisum abyssinicum*, ii) *Pisum fulvum* and iii) *Pisum sativum* L., which further includes the wild pea, *Pisum sativum* subsp. *elatius* (M. Bieb. Asch. & Graebn) and the domesticated pea, *Pisum sativum* subsp. *sativum*.

Phylogenetic analyses of various pea taxa with molecular markers indicate that hybridization between wild peas is not an extensive phenomenon [8]. The recently annotated pea genome sequence and the resequencing of data from 42 wild, landrace and cultivar *Pisum* genotypes, provided further insights into legume genome evolution [9]. It has been suggested that the common ancestor of the *Pisum* species was probably cytogenetically similar to *P. sativum* subsp. *elatius*, which evolved across the Mediterranean and Middle East [8] and gave rise to *P. sativum* subsp. *sativum* and *P. fulvum* in the northern Middle East. Regarding *P. abyssinicum*, two main hypotheses exist; it is considered the result of a domestication event from a southern *P. sativum* subsp. *elatius* ancestor [10] followed by a migration to Abyssinia, possibly through ancient human trading routes [11], indicating at least two domestication events independent of *P. sativum* subsp. *sativum* [9,12]. The alternative hypothesis about the origins of *P. abyssinicum* suggests that it derived from a hybridization event between *P. fulvum* and *P. sativum* subsp. *elatius*, which occurred in the western half region of the Fertile Crescent [13] and then a small sample was introduced to north-eastern Africa, where it evolved into the modern *P. abyssinicum* [14]. The very low genetic diversity in the Abyssinian pea suggests that the taxon has recently experienced a severe bottleneck or is a relatively young taxon [10] and the hybridization event has most likely occurred about 4000 years bp [15].

The *Pisum* genus is very diverse, showing the gamut of relatedness that reflect taxonomic identifiers, eco-geography and breeding gene pools [8,12,16,17]. Several phenotypic classification studies on pea germplasm are based on agronomical characteristics and morphological descriptors [18–21], which are unreliable for the evaluation of pea genetic resources and the identification of different cultivars in the Fabaceae family [22], especially considering the environmental effects on the expression of the genotype. Several different molecular methods have been previously employed to assess the genetic diversity in the *Pisum* genus, such as Random Amplification of Polymorphic DNA (RAPD) [23–25], Inter-Single Sequence Repeats (ISSRs) [23,24,26], Single Sequence Repeats (SSR) [27,28], Retrotransposon-Based Insertion Polymorphism (RBIP) markers [13,17,29–31] and Expressed Sequence Tags (EST)-derived genomic markers [32]. Additionally, high-throughput parallel genotyping via genome-wide next generation sequencing techniques have also been used to study the diversity of wild pea [8,33,34].

An alternative method for the simple and accurate authentication of plant species is DNA Barcoding. The *rbcL* and *matK* regions have been recommended as core DNA barcodes for plant identification [35]. In *Fabaceae* species, four coding chloroplast regions (*rpoB*, *rpoC1*, *rbcL*, and *matK*) and two non-coding nuclear regions (*ITS1* and *ITS2*) have been used as barcodes [22,36–38]. DNA barcoding has also been used to reconstruct the phylogenetic relationship of the main Mediterranean leguminous crops [39]. Furthermore, the combination of DNA barcoding with high resolution melting analysis (Bar-HRM) has, thus far, been proved an effective approach for the identification of diverse plant species, their Protected Designation of Origin (PDO) products and quantification of adulterants [40–43].

In Greece, pulses have been traditionally cultivated since the ancient times and is a staple food in the local culinary culture. A popular Greek dish (namely 'Fava') is typically prepared using different legume species, such as yellow- split peas (*P. sativum* L.) or faba beans (*Vicia faba*). However, in the island of Santorini, the authentic PDO 'Fava Santorinis', is exclusively prepared from a local grass pea variety of *L. clymenum*. In other Cycladic islands, especially in Amorgos and Schinoussa, 'Fava' is prepared from the dried peeled and split seeds of an endemic *Pisum* cultivar also known as 'Katsouni', named after the convex shape of a small sickle's lobe used for mowing the crop. 'Katsouni' is a traditional product and a crop of great economic importance for the Lesser Cyclades. Currently, is in the process to be appointed as a PDO EU mark, offering a significant income to the local farmers. 'Katsouni' is fully adapted to the local climatic conditions of the Cyclades, with dry and hot summers and mild winters. It is rich in proteins (over 22%) and can be stored after drying the seeds throughout the year.

Historic records indicate that the 'Katsouni' landrace is the result of long-term selection and evolution from prehistoric times that occurred in Amorgos. Since the mid-19th century 'Katsouni' was transferred to the deserted Schinoussa island by residents of neighbouring Amorgos, who moved to settle there, and has since been cultivated uninterruptedly. However, Schinoussa was not always an abandoned island. Archaeological excavations have revealed findings indicating great activity from the prehistoric times to the Classical and Hellenistic period [44]. In the Byzantine times, trade and commerce were essential components of the island's prosperity [45], which heralded an age of advancement, especially during the Venetian rule (13th–16th century). However, during the Ottoman rule (16th–19th century) the island was deserted as indicated by the famous botanist and traveler de Tournefort (2003) [46]. In the late Middle Ages, in 1537, Cyclades along with Schinoussa were plundered by the Ottoman pirate Haiderin Barbarossa [47] and piracy continued up to the 19th century [48]. With Schinoussa becoming a pirates' den, we hypothesized that either: (i) the current Schinoussa pea population is an independent landrace, which is possibly the result of introgression of *P*. *abyssinicum* (transferred by pirates) into the *P. sativum* subsp. *elatius* germplasm; (ii) the Amorgos and Schinoussa populations belong to the same widely distributed landrace, or (iii) the two landraces emerged from the split of an ancestral population.

The knowledge of genetic relationships and diversity among individual landraces is fundamental for conservation practices and the selection of appropriate parents in breeding programs. Hence, in the present study, we evaluated the application of ISSR marker analysis and DNA barcoding for the molecular characterization of local Amorgos and Schinousssa pea populations. Furthermore, considering that powdery mildew (*E. pisi* DC.) severely affects pea crops worldwide [49], and in the frame of targeting functional polymorphisms, we aimed at characterizing the two landraces for the presence of the powdery mildew resistance gene (*er1-7)* with HRM analysis. The diversity assessment of local landraces may not only provide insights in understanding pea phylogenetics and population genetics, but also broaden pea breeding strategies.

#### **2. Results**

#### *2.1. DNA Barcoding, Sequencing and Tree Analysis*

To identify potential inter- and intra-specific variation between the two pea populations of Amorgos and Schinoussa, 24 pea samples were analyzed in the present study using the *ITS2*, *trnL* and *rpoC*, and 21 for the *psbA*-*trnH* and *matK* barcoding regions. The selection of *matK* and *psbA*-*trnH* was based on the unique SNPs observed in the aligned sequences among the three species and especially between *P. abyssinicum* and *P. sativum* (Figure S1). The amplification of *trnL* and *rpoC* regions was 100% successful, whilst for the *psbA-trnH*, *ITS2* and *matK* sequences the rate was 95.24%, 95.8% and 90.47%, respectively.

The positive amplicons were sequenced, and based on the BLAST results, most of the markers were able to identify the samples at the genus level, but not at the species level, when blasted against the NCBI database. The BLAST entries matched with all the three *Pisum* species (*P. sativum* subsp. *sativum*, *P. sativum* subsp. *elatius*, *P. fulvum* and *P. abyssinicum*), with over 96% similarity in identity and coverage. However, the *ITS2* marker identified all target sequences as *P. sativum* subsp. *elatius* with 100% similarity to the available sequences.

The sequences of each barcode gene were aligned and compared against indicative BLAST entries. In *psbA-trnH* we observed an intel polymorphism among our samples and the NCBI database entries. Most of our samples showed a gap in the barcoding region similar to that observed in the *elatius* subspecies sequences, but not in the *P. abyssinicum* and *P. fulvum* sequences (Figure S1). Other observed variations include the *P. abyssinicum*-specific polymorphism of a thymine (T) in place of a cytocine (C), and an unspecific SNP (G/T) present exclusively in the samples of both islands and the *P. sativum*

subsp. *elatius* NCBI entries (Figure S1). Regarding the *ITS2* barcode, only the brown colored peas from both islands shared a common nucleotide variation of a guanine (G) in place of an adenine (A) (Figure S2). The *matK* barcoding showed a *P. fulvum*-specific polymorphism (C/A) and other unspecific SNPs (Figure S2). However, the observed variations in *trnL* and *rpoC* regions were not consistent among either seed coat color and/or landrace (Figure S2). Furthermore, both *matK* and *psbA-trnH* regions showed higher diversity and pairwise distance values compared to the other markers (Table S1); yet, the observed low genetic variation values across all barcoding regions (Table S1) indicate that the populations are closely related with a recent common ancestor.

To provide a basic illustration of the phylogenetic associations between the landraces, the DNA barcoding data was used to calculate the genetic distances and generate Neighbour-joining dendrograms (Figure 1 and Figure S3). The *psbA-trnH* dendrogram illustrates the clustering of the *P. abyssinicum* sequences supported by 76% bootstrap value, which are distinct from the rest of the samples at 64%, emphasizing that our samples probably do not belong to *P. abyssinicum* (Figure 1). In contrast, both the main structures of the *trnL* and *ITS2* dendrograms (Figure S3) separated the two landraces in two distinct clusters, which also corresponded to the represented geographical regions, although supported by a low bootstrap value (30–60%). The main structure of the *rpoC* dendrogram (Figure S3) presented also two clusters separating the pea taxa from Schinoussa to those of Amorgos island, supported by a higher bootstrap value (59–66%). However, the *matK* dendrogram (Figure S3) did not show any significant patterns. Taking into consideration the results from the barcoding analysis, the two populations probably belong to the *P. sativum* subsp. *elatius*.

**Figure 1.** The *psbA-trnH* dendrogram illustrating the phylogenetic relationships between Amorgos (red circle) and Schinoussa (blue triangle) pea populations. The corresponding NCBI sequences of the *P. sativum* subsp. *elatius* (purple square), *P. sativum* subsp. *sativum* (white square), *P. abyssinicum* (yellow square) and *P. fulvum* (green square) were used as reference taxa.

#### *2.2. ISSR Genotyping*

To further investigate the genetic differences between the two populations, we used fifteen ISSR markers, out of which six were found to be polymorphic. ISSR analysis of the pea populations using the six polymorphic markers yielded 66 bands in total; five unique bands were identified for the Schinousasa population, whereas the remaining 61 bands were shared between the two populations (Table 1). The Schinoussa population presented a significantly greater polymorphism compared to Amorgos by displaying a greater number of different alleles (Na), whilst the number of effective alleles (Ne) did not show significant differences (Table 2).

**Table 1.** Band patterns of the Amorgos and Schinoussa populations resulted from the ISSR analysis.

**Table 2.** Mean value and standard error over loci for Amorgos and Schinoussa populations.

Mean = Mean value, SE = Standard error, N = Number of alleles, Na = Number of different alleles, Ne = Number of effective alleles = 1/(pˆ2 + qˆ2), I = Shannon's Information Index = −1\* (p \* Ln (p) + q \* Ln(q)), h = Diversity = 1 − (pˆ2 + qˆ2), uh = Unbiased diversity = (N/(N − 1)) \* h.

The genetic differentiation between landraces (PhiPT distances) was significantly different (PhiPT = 0.188; P ≤ 0.001), indicating that the two landraces are geographically distinct. The Analysis of Molecular Variance (AMOVA) based on the PhiPT values indicated that most of the genetic diversity occurred within landraces (80%), while the variability among landraces contributed to the 20% of the observed genetic diversity (Table 3).

**Table 3.** AMOVA analysis of the Amorgos and Schinoussa pea populations.

Df = Degrees of freedom, SS = Sum of Squares, MS = Mean Square, Est. Var. = Estimated Variance.

The Principal Coordinates Analysis (PCoA) generated two major clusters, in which samples from Amorgos and Schinoussa were clearly separated (Figure 2A). Additionally, the clustering based on seed coat color supports that the two landraces present genetic differences (Figure 2B), given Amorgos' Black, Brown and Green peas do not overlap with the corresponding colors of Schinoussa peas, supporting the hypothesis that the Amorgos and Schinoussa pea populations are two distinct landraces. The UPGMA dendrogram illustrates that the pea samples originating from the same geographic location are closely clustered (Figure 3).

**Figure 2.** Principal Coordinates Analysis (PCoA) of the Amorgos and Schinoussa pea populations, clustering for: (**A**) population, and (**B**) seed coat color from each landrace. The PCoA analysis shows the separation of the two populations as distinct landraces based on their region of origin.

**Figure 3.** Dendrogram of the Amorgos and Schinousa pea populations based on the UPGMA analysis of the ISSR polymorphisms. Individuals are shape- and color-coded based on their region of origin (Red circle: Amorgos, Blue triangle: Schinoussa).

#### *2.3. Molecular Screening for Powdery Mildew Using HRM Analysis*

HRM analysis using a specific functional marker was employed to characterize the two landraces for the presence of the powdery mildew resistance gene. Regarding the specific InDel marker, all the pea samples tested were found negative to the resistant allele (*er1-7*), across both landraces (Figure 4A). This was also confirmed by sequencing of selected samples, where the resistant allele was absent and therefore the 10-bp sequence TCATGTTATT was present (Figure 4B).

**Figure 4.** HRM analysis coupled with a co-dominant functional marker specific for *er1-7.* (**A**) Normalized fluorescence graph of selected pea samples per seed coat color and landrace. In the x axis deg. indicates temperature in ◦C. (**B**) Sequence alignment of the *er1-7* region from selected pea samples and the corresponding reference sequences obtained from Sun et al. (2016) [50].

#### **3. Discussion**

Exploiting genetic diversity from local traditional landraces is fundamental for conservation practices and breeding programs, especially under the pressure posed for adaptation to climate change worldwide. To promote the local pea landraces of the Cyclades we aimed at identifying the species and understanding the genetic relationship of the local Amorgos and Schinousssa pea populations, using DNA barcoding and ISSR marker analysis. Over the centuries, from the prehistoric times and the Bronze Age (3000–1100 B.C in Greece) to the mid-19th century, the Aegean was a field of pirate action [51,52]. Oral traditions and place names, throughout the history of Schinoussa, suggest that the island was used as a ground and shelter for pirates. One of our hypotheses was that the landrace of Schinoussa belongs to *P. abyssinicum*, however, our results strongly suggest otherwise. Although both landraces showed a geographical clustering according to the PCoA and the UPGMA analysis of the ISSR polymorphisms, both were identified based on DNA Barcoding as *P. sativum* subsp. *elatius*.

Herein, the inability of most of the barcoding markers to discriminate samples at the species level, being not variable enough to resolve phylogeny of the genus, is probably due to the conserved chloroplast sequences and the low mutation rate [53]. The inability to discriminate among the *P. abyssinicum*, *P. sativum* subsp. *elatius* and *P. fulvum* is in agreement with the view that *P. abyssinicum* is an ancient hybrid of the two species [14]. Nevertheless, despite the inability of the *psbA-trnH* spacer on identifying species due to the frequently observed intraspecific inversions [54], in this study it was shown as the most effective marker for separating the two pea landraces from *P. abyssinicum*. Additionally, although the *matK* barcode showed the lowest resolution capacity in our study, it was sensitive enough to discriminate the two pea landraces from *P. fulvum*. Thus, the most informative

barcoding markers to draw conclusions concerning the species identification in the two populations is the combination of *psbA-trnH* and *matK*. Taking into consideration the phylogenetic trees and barcoding results from all the studied markers, the two landraces are more closely related to *P. sativum* subsp. *elatius*, than to the *P. abyssinicum* and *P. fulvum* species.

There is a limited number of studies on *Pisum* germplasm that mainly are focusing on SSR analysis [55]. Herein, the ISSR analysis showed that the populations were distinguished into geographical regions, as separated clusters, indicating the adaptation of these traditional landraces to relatively different agro/climatic conditions (Figures 2 and 3). This geographical isolation could have potentially resulted to the genetic drift of the two landraces. This result is in accordance with the larger molecular differences detected between pea landraces collected in Maritime areas of Spain [26]. Pea is known as self-pollinating with occasional cross-pollination which allows spontaneous hybridization [56]. As a self-pollinated crop, higher genetic diversity is expected among cultivars than within cultivars. However, the AMOVA analysis showed larger genetic diversity within cultivars (80%), which is in accordance with similar findings in other legumes such as chickpea [57]. This may be attributed to the natural interspecific crosses that can occur between *Pisum* species, serving as a source of additional genetic diversity for the selection of common pea [15]. This indicates that despite autogamy, the analysis of genetic diversity on some plants per landrace might be useful in breeding programs [28]. The large genetic diversity might be due to the long-term adaptation of the landraces to the local environment and the diverse agro-ecological systems [55], in combination with putative migration events among the regions, followed by introgression with pre-existing germplasms [28].

Aiming to identify unique traits in the two pea traditional landraces we screened for the powdery mildew resistance gene (*er1*). The er1 is the loss-of-function mutation in the disease susceptibility-related *PsMLO1* gene [58] and the most widespread across resistant pea germplasm, conferring penetration resistance [50,59]. HRM coupled with a functional marker has been previously reported to be highly efficient and cost-effective for routine large-scale screening of pea germplasm for resistance to powdery mildew [60]. This approach allowed for the accurate genotyping of both the homozygous and heterozygous resistant peas [58,60]. In this study, using the *er1-7* functional marker coupled with HRM, we identified that the two landraces of the Cycladic islands are susceptible to powdery mildew. These results indicate the need for further crop improvement of the traditional pea landraces by introducing resistance through molecular-assisted breeding strategies.

#### **4. Materials and Methods**

#### *4.1. Plant Material and DNA Isolation*

The plant material used in this work includes samples from two pea populations from the islands of Amorgos and Schinoussa (Table 4 and Figure 5). The pea samples were grouped based on the geographical region of origin and the seed coat color; groups A-C from Amorgos and groups D-G from Schinoussa (Table 4). Color grouping was based on the three main seed coat colors (black, brown and green) observed in both populations, except for Schinoussa, which had an additional brown-green hue (Table 4). Following, a subsample of seeds from each group was planted in pots containing 2:1 peat:perlite in order to obtain fresh leaf material for DNA extraction. Total genomic DNA was isolated from approximately 0.1 g of fresh leaf material for each sample following the modified CTAB protocol as described by Doyle and Doyle (1987) [61]. After extraction, the DNA samples were re-diluted in 1X TE buffer (10 mM Tris-Cl pH 8.0, 1 mM EDTA) and stored at −20 ◦C. DNA quantity and quality were assessed by regular spectrophotometric procedures using the UV-Vis Spectrophotometer Q5000 (Quawell Technology Inc., U.S.A.) and gel electrophoresis in 1% agarose gel.

**Table 4.** Samples from Amorgos and Schinoussa pea populations used in this work. The samples were grouped in sub-groups (A–G) according to region of origin and seed coat color.

**Figure 5.** Map of the sampling sites in Amorgos and Schinoussa islands.

#### *4.2. DNA Barcoding and Sequencing Analysis*

For the identification of the two landraces, we performed DNA barcoding analysis using the *ITS2* [62], *trnL* [63], *rpoC*, *matK* [35], and *psbA-trnH* [64] barcoding markers. PCR amplification was performed on a Rotor-Gene 6000 real-time 5-Plex HRM PCR Thermocycler (Corbett Research, Sydney, Australia), using the Rotor-Gene Q software version 2.0.2 (Corbett Life Science, Cambridge, UK).  The universal regions were amplified using the following protocol: initial denaturation at 95 ◦C for 4 min, followed by 35 cycles of 95 ◦C for 30 sec, corresponding annealing temperature (Ta) ◦C for 60 sec, and 72 ◦C for 60 sec with a final extension phase at 72 ◦C for 3 min.

After sequencing, the sequences of the five candidate regions were aligned with the MUSCLE algorithm and genetic distances were calculated using Molecular Evolutionary Genetics Analysis X (MEGA X; Version 10.05) based on the K2P-distance model [65] to evaluate divergence between the two populations. The Neighbour-joining clustering method was used to demonstrate the represented differences as an unrooted dendrogram using MEGA X [65]. Statistical support for each constructed tree was provided by two statistical data analysis as bootstrapping (1000 replications) and pairwise distance model.

Species identification based on the sequence similarity approach was performed with the National Center for Biotechnology Information (NCBI) database [66] by basic local alignment search tool (BLAST; setting: blastn, megablast) [67] and all regions of the three various *Pisum* species were used as query sequences. Correct identification was concluded when the best BLAST hit of the query sequence had over 96% query coverage and identity.

#### *4.3. ISSR Genotyping and Data Analysis*

For the distance-based analysis of the two populations we used 15 ISSR markers of which six (UBC811, UBC818, UBC827, UBC841, UBC873, UBC880) were selected for further analysis based on their discrimination efficiency. The total volume of PCR reaction was 25 μL containing 1X PCR buffer, 0.2 mM dNTPs, 10 μM primer, 1 U/μL Taq DNA Polymerase (Kapa Biosystems Ltd.) and 20 ng template DNA. The profile of the PCR reaction program was an initial denaturation for 4 min at 94 ◦C, followed by 35 cycles at 94 ◦C for 30 s, 40 sec annealing at the corresponding Ta ◦C, and 40 sec extension at 72 ◦C, ending with final extension phase at 72 ◦C for 7 min. The amplified PCR products were run on 1.5% agarose gel with 1X TAE buffer at 100 V and visualized with the UV Minibis Pro (DNR Bio-Imaging Systems, Jerusalem, Israel) instrument. Band scoring was performed using the Logger Pro 3.15 software.

DNA fragment profiles were scored in a binary fission with '0 indicating the absence and '1 indicating presence of a band. Using the binary haploid data, a pairwise individual-by-individual genetic distance matrix was constructed using the Jaccard coefficient. The percentage of polymorphic loci (P), Number of alleles (N), Number of different alleles (Na), Number of effective alleles (Ne), gene diversity (expected heterozygosity, He), Shannon's diversity index (I), Diversity (h) and unbiased genetic distances (uh) according to [68] were subsequently calculated. The hierarchical distribution of genetic diversity among and within populations was also characterized by Analysis of Molecular Variance (AMOVA) and Principal Co-ordinate Analysis (PCoA). All of the above analyses were performed using the GenAlex 6.5 software package [69]. The Unweighted Pair Group Method based on Arithmetic Averages (UPGMA) clustering analysis for analyzing the similarity estimates was performed using MEGA X [65] and expressed as a dendrogram.

#### *4.4. Molecular Screening for Powdery Mildew Resistance Using HRM Analysis*

We used a co-dominant functional marker specific for *er1-7*, the InDel111–120, associated with pea resistance to powdery mildew [50,59,60]. Representative samples were selected from each color group of both populations. The total volume of PCR reaction was 20 μL containing 1X PCR buffer, 0.2 mM dNTPs, 10 μM of the er1-7 primer, 1 U/μL Taq DNA Polymerase (Kapa Biosystems Ltd.) and 20 ng template DNA. PCR conditions where: preheating for 5 min and initial denaturation at 95 ◦C, followed by 35 cycles at 94 ◦C for 30 s transition, 30 sec annealing at 56 ◦C. The fluorescence was measured at the end of each extension step during the PCR cycles. The HRM was performed by an initial pre-melt conditioning of the PCR products at 95 ◦C for 5 sec and 50 ◦C for 30 sec, followed by a melt at range of 75–85 ◦C in increments of 0.1 ◦C every 2 sec. Fluorescence was measured at the end of each increment.

#### **5. Conclusions**

This work represents the first investigation focused on the molecular characterization of pea traditional accessions collected from the semi-arid region of Lesser Cyclades with the main goal to evaluate their genetic diversity and their population structure. Based on the DNA Barcoding and the ISSR marker analysis: i) both landraces have been identified as *P. sativum* subsp. *elatius* germplasm and ii) the landraces show a tendency for differentiation, which is in accordance with the geographical distribution of the genetic structure as an underlying evolutionary process. However, further research and phylogenetic analysis is required with larger population sizes for better understanding of the evolutionary processes that led to these differences, as well as for the preservation of existing diversity in *ex-situ* collections.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2223-7747/9/4/479/s1, Table S1: Pairwise distance and diversity analysis within and between populations of the five barcodes, Figure S1: Sequence alignments of the *psbA-trnH* region, Figure S2: Sequence alignments of the *ITS2*, *trnL*, *rpoC* and *matK* regions from selected pea samples and the corresponding Pisum reference sequences available in the NCBI, Figure S3: The *ITS2*, *trnL*, *rpoC* and *matK* dendrograms illustrating the phylogenetic relationships between Amorgos (red circle) and Schinoussa (blue triangle) pea populations.

*Plants* **2020**, *9*, 479

**Author Contributions:** Conceptualization, P.M.; methodology, E.S., G.L., L.K., M.O.; data curation, E.S., G.L.; writing—original draft preparation, E.S.; writing—review and editing, E.S., G.L., M.O.; visualization, G.L., E.S.; supervision, P.M.; funding acquisition, P.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has been co-financed by the European Regional Development Fund of the European Union and Greek national funds through the Operational Program Competitiveness, Entrepreneurship and Innovation, under the call RESEARCH-CREATE-INNOVATE (project code: T1EDK-04448).

**Acknowledgments:** We would like to thank Dimitris Roukas for providing the pea landrace seeds and vital information regarding the origin of the seeds and the history behind its cultivation in the two islands. We would also like to thank the University of Chiang Mai for financial support.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Evaluation of Wild Potato Germplasm for Tuber Starch Content and Nitrogen Utilization E**ffi**ciency**

#### **Silvia Bachmann-Pfabe \* and Klaus J. Dehmer**

Leibniz Institute of Plant Genetics and Crop Plant Research (IPK), Genebank Department, Satellite Collections North, Gross Luesewitz Potato Collections, Parkweg 3a, 18190 Gross Luesewitz, Germany<EMAIL_ADDRESS>
**\*** Correspondence<EMAIL_ADDRESS>Tel.: +49-38209-82-314

Received: 10 June 2020; Accepted: 1 July 2020; Published: 2 July 2020

**Abstract:** Potato wild relatives provide a considerable source of variation for important traits in cultivated potato (*Solanum tuberosum* L.) breeding. This study evaluates the variation of tuber starch content and nitrogen utilization efficiency (NutE) in wild potato germplasm. For the experiments regarding starch content, 28 accessions of ten different tuber-bearing wild *Solanum*-species were chosen, and in vitro plantlets were raised from seeds. Twenty plantlets (= genotypes) per accession were then cultivated in the greenhouse until natural senescence and tuber starch content was determined. The average tuber starch content across all genotypes tested was 21.7% of fresh mass. Contents above 28% of fresh mass were found in 50 genotypes, belonging to the species *S. chacoense*, *S. commersonii*, *S. jamesii*, and *S. pinnatisectum*. Subsequently, 22 wild genotypes revealing high tuber starch contents and four modern varieties of cultivated potato were studied as in vitro plantlets under optimal and low N supply (30 and 7.5 mmol L−1 N). Low N supply lead to a genotype-dependent reduction of shoot dry mass between 13 and 46%. The majority of the wild types also reduced root dry mass by 26 to 62%, while others maintained root growth and even exceeded the NutE of the varieties under low N supply. Thus, wild potato germplasm appears superior to cultivars in terms of tuber starch contents and N utilization efficiency, which should be investigated in further studies.

**Keywords:** genetic resources; *Solanum chacoense*; stress tolerance

#### **1. Introduction**

Besides being one of the most important food crops worldwide, potato (*Solanum tuberosum* L.) plays an important role in industry due to its starchy tubers. Starch is used in bakery products, thickening products, soups and noodles but also for the production of paper, textiles, building materials, pharmaceutical products, chemicals and biodegradable packaging materials [1]. Compared to other starches, potato starch has superior characteristics because it is easily isolated, of high purity and of large granule size, needs low temperatures for gelatinization and produces gels with highest viscosity [2]. Depending on genotype and growing environment, a fresh potato tuber contains about 20% of dry mass, 60–80% of which is starch [3]. Dedicated starch varieties may even reach a starch content of up to 23% of fresh mass [4]. The nitrogen (N) fertilizer regime influences starch yield by positively affecting canopy development and photosynthesis efficiency, dry matter partitioning to the tubers, tuber bulking and tuber yield formation [5,6]. Furthermore, field N availability influences starch quality parameters, such as granule size, viscosity and breakdown [7]. Consistently, Maltas et al. [8] reported a highly significant effect of different N fertilizer rates on total tuber yield, the percentage of large tubers and starch concentration under field conditions in cv. Bintje and Laura.

Depending on environmental conditions and genotype, the potato crop has been found to remove 90 to 190 kg N ha−1 [9]. However, the shallow, less branched and less dense root system of potato does not allow the exploration of a large soil volume or to retrieve nitrogen from deeper soil layers, and hence, potato demands a high level of readily available soil N at the right period of growth [5,9–12]. In a review, Iwama [13] reported that most of the potato roots are present in the upper 30 cm of the soil and only a small fraction extends to 100 cm. In combination with the fact that potatoes are often cultivated on coarse and sandy soils and under irrigation, these areas face an increasing potential for nitrate leaching and contamination of groundwater [14].

Thus, improving the N use efficiency in potato production is not only of economic, but also of environmental concern, and different measures such as split application of N or foliar application of urea are being discussed in order to reduce N leaching [5]. Improving the N uptake and use efficiency of the potato crop itself is also an important approach, and many studies evaluated the N use efficiency in cultivated potatoes [8,15–21]. In contrast, only few studies evaluated the N use efficiency and/or tuber starch contents of native Andean cultivars or wild potato germplasm [22–25]. The wild relatives of the cultivated potato could be an important source of variation for root length and morphology, tuber starch content and N use efficiency. The secondary and tertiary genepool of potato has intensively been studied as a source of disease resistance [26–28], and was used, amongst others, to improve foliar late blight and nematode resistance of *S. tuberosum* [29,30]. In terms of N use efficiency, Errebhi et al. [22] compared 39 wild potato accessions of 23 species with three cultivated varieties under high and zero N in the field. They found some wild potato genotypes which were able to take up significantly more applied N than their cultivated relatives. Genotypes of *S. microdontum* and *S. chacoense* were identified as the ones with the highest N uptake efficiency (NupE) [22]. Selected native Andean cultivars indicated a similar nitrogen use efficiency to commercial cultivars, but showed, despite different environmental conditions, a highly consistent performance across a two-year field study [24]. To our knowledge, the most comprehensive study of tuber starch content and quality in exotic germplasm, was provided by Jansen et al. [25]. Accessions of 46 wild and cultivated potato species showed a high variation in starch contents ranging between 3.8 and 39.6% of fresh mass. Highest starch contents were predominately found in genotypes of species *S. pinnatisectum* and *S. chacoense* [25].

Based on the above-mentioned findings, our study aimed to (I) update and assess the variability of the tuber starch contents in wild potato germplasm and to (II) study the nitrogen use efficiency of genotypes with high tuber starch contents in relation to modern cultivars.

#### **2. Results**

#### *2.1. Variation of Tuber Starch Contents in Wild Potato Germplasm after Greenhouse Cultivation*

In 2013, altogether 28 different wild potato accessions (= populations) representing ten different species were cultivated in the greenhouse to evaluate their tuber starch contents. For each accession 20 different genotypes were cultivated as in vitro plantlets, however, the results only include genotypes which produced sufficient tubers for starch analysis (in total 506 genotypes, Table 1). On average, of all the 506 genotypes analyzed, the starch content in the tubers amounted to 21.7% of fresh mass (FM). The lowest average starch content with 14.2% of FM was measured for accession Gross Luesewitz Potato Collections (GLKS) 31559 (*S. stenotomum*), while accession GLKS 30211 (*S. commersonii*) showed the highest starch content with on average 30.0% of FM. Interestingly, all genotypes of accession GLKS 30211 showed high tuber starch contents ranging between 26.4 and 33.3% of FM, indicating a rather low variation within this population (Table 1). From *S. chacoense,* all the 15 accessions had an average starch content of 22.5% of FM and showed a rather low variation between the populations (CV = 8.71%), and a higher within population variation (CV: 11.5–23.9%, Table 1). Ten accessions of *S. chaoense* encompassed genotypes with a tuber starch content higher than 28% of FM, a target value which was considered as selection criterion for high-starch genotypes in this study. Regarding *S. pinnatisectum*, the tubers of the five accessions studied, had an average starch content of 22.2% of FM, but showed starch contents up to 36.6% of FM. Altogether 11 genotypes of three accessions produced starch contents higher than 28% (Table 1).

**Table 1.** Wild potato accessions of the Gross Luesewitz Potato Collections (GLKS) used in the study ranked according to their tuber starch content. Given are the number of genotypes tested per accession, the mean starch content (%) in the tubers, its minimum and maximum values and the coefficient of variation (CV %).

For subsequent N efficiency experiments, genotypes with starch contents higher than 28% of FM were to be used. In *S. chacoense*, altogether 23 genotypes were identified revealing a tuber starch content above 28% of FM. Overall, 15 genotypes, belonging to accessions GLKS 30135, GLKS 30154, GLKS 30156, GLKS 30159, GLKS 30160, GLKS 30177, GLKS 30181, GLKS 30916, and GLKS 30995 were selected for the N experiments and re-cultivated in 2014 to validate their starch contents (Table 2). In *S. pinnatisectum*, 11 genotypes showed tuber starch contents above 28% of FM, and genotype GLKS 31600\_10 with a starch content of 36.6% was selected for up-coming experiments and re-tested in 2014. Additionally, two genotypes with the highest starch contents of *S. microdontum* (GLKS 30688\_04, GLKS 30688\_12) and *S. stenotomum* (GLKS 31559\_11, GLKS 31559\_14) as well as one genotype of *S. tuberosum* subsp. *andigena* (GLKS 34995\_18) were selected for the N efficiency experiments in order to cover a broader spectrum of *Solanum* species, even if they had starch contents below 28% of FM (Table 2). The tuber starch contents determined in 2013 and after re-testing of selected genotypes in 2014 correlated well (r = 0.72, *p* ≤ 0.01), confirming the high-starch properties of the majority of the selected accessions.

**Table 2.** Selected accessions and genotypes of the Gross Luesewitz Potato Collections used for the N efficiency experiments as well as their respective tuber starch contents (%) in 2013 and 2014.

#### *2.2. Dry Yield of Shoots, Roots and Root-DM:Shoot-DM Ratio in the N Experiments*

Shoot and root DM as well as the root-DM:shoot-DM ratio were predominately affected by the genotype and, to a lesser extent, by the factor treatment. For these traits, the factor genotype explained up to 68% of the variation in the data, while the factor treatment explained between 6 and 22%. The genotype × treatment interaction explained 11 and 15% of the total variation for root DM and root-DM:shoot-DM ratio, respectively, but played only a minor role for shoot DM (3.58%, Table 3).

In the high N (30 mmol L−1) treatment, the shoot DM of the genotypes varied between 214 and 682 mg vessel−1 (Table 4). The lowest shoot DM was observed for the genotypes GLKS 31600\_10, GLKS 30177\_17 and cv. Kiebitz, while the highest biomass was achieved by the genotypes GLKS 30177\_20, GLKS 30181\_06 and GLKS 30160\_15. These genotypes even exceeded the performance of cv. Tomba, which showed the highest shoot DM amongst the standard varieties. Under low N supply (7.5 mmol L−1), shoot DM ranged between 156 and 549 mg vessel−1. Shoot DM decreased under low N supply on average by 115 mg vessel−1 (23%) and the shoot biomass reduction was significant for all the genotypes tested, except for cv. Kiebitz and cv. Eurobravo (Table 4). The strongest reduction (>30%) was observed for genotypes GLKS 30135\_19, GLKS 30995\_18, GLKS 30177\_02 and GLKS 31559\_11. A moderate shoot DM reduction (15 to 20%) at simultaneously high yields in the control was observed for genotypes GLKS 30135\_05, GLKS 30160\_13 and GLKS 30177\_20. Again, cv. Tomba produced the highest shoot DM amongst the standard varieties under low N supply. However, several wild potato genotypes performed as well or even exceeded the shoot DM of cv. Tomba under reduced N conditions (Table 4).

 The letters within one column indicate whether there is a significant difference to the best cultivar Tomba ("b") or not ("a", Dunnett's test *p* ≤ 0.05), underlined are mean significantly higher than that of cv. Tomba; asterisks indicate a significant difference between the treatments within one genotype (pairwise comparisons, Tukey adjustment; \*\*\* *p* ≤ 0.001,\*\**p* ≤ 0.01, \* *p* ≤ 0.05, ns = not significant).

The root DM varied between 30 and 337 mg vessel−1 in the control treatment. The lowest root growth (<80 mg vessel−1) was observed for genotypes GLKS 31600\_10, cv. Kiebitz, GLKS 30181\_18 and GLKS 31559\_14, while the genotypes GLKS 30154\_09, GLKS 30160\_13 and GLKS 30177\_20 produced more than 300 mg vessel−1. Under low N conditions, the root DM varied between 41 and 219 mg vessel−1. The different genotypes either maintained or decreased root DM due to N deficit. The decrease was highest (about 60%) for the genotypes GLKS 30135\_05, GLKS 30177\_02, GLKS 30181\_06, while the genotypes GLKS 30177\_15, GLKS 30181\_18, GLKS 30688\_04, GLKS 30688\_12, GLKS 31559\_14, GLKS 31600\_10 and cv. Kiebitz maintained root mass. Most interestingly, genotype GLKS 30177\_15 produced a high root mass in the control and maintained it under reduced N conditions. Amongst the cultivars, cv. Tomba produced the highest root mass in both treatments. However, several wild potato genotypes produced a significantly higher root biomass under control (GLKS 30160\_13, GLKS 30177\_20) or low N conditions (GLKS 30160\_15, GLKS 30177\_15, GLKS 30177\_20) compared to cv. Tomba (Table 4).

Relating the root to the shoot biomass improves the understanding of genotype-specific reactions to low N supply. The root-DM:shoot-DM ratio varied strongly between the wild types, ranging from 0.14 to 0.61 in the control and 0.10 to 0.48 in the low N treatment, but was generally comparable to that of the standard varieties (Table 5). High values above 0.50 were calculated for genotypes GLKS 30154\_09, GLKS 30177\_01, GLKS 30160\_13 and GLKS 30177\_17, indicating a strong root growth in relation to the shoot. Genotypes reacted differently to low N supply, either by maintaining (13), reducing (11) or increasing (2) root-DM:shoot-DM ratio. The genotypes GLKS\_31559\_14 and GLKS\_31600\_10 showed the lowest root-DM:shoot-DM ratio in the control, but increased root growth at the expense of the shoot under low N.

**Table 5.** Mean root-DM:shoot-DM ratio as well as partitioning of N taken up by the different wild potato genotypes (root-N:shoot-N ratio) cultivated under optimal and reduced N conditions in a climate chamber for 21 days. Results of ANOVA and post hoc comparison of means.

The letters within one column indicate whether there is a significant difference to the best cultivar Tomba ("b") or not ("a", Dunnett's test *p* ≤ 0.05), underlined are mean values significantly higher than that of cv. Tomba; asterisks indicate a significant difference between the treatments within one genotype (pairwise comparisons, Tukey adjustment; \*\*\* *p* ≤ 0.001, \*\* *p* ≤ 0.01, \* *p* ≤ 0.05, ns = not significant).

#### *2.3. N Uptake, N Partitioning and N E*ffi*ciency Parameters in the N Experiments*

Interestingly, the genotypic variation for the shoot and total N uptake was low and explained only about 4% of the total variation. In contrast, the factor genotype explained between 20 and 60% of the total variation for the traits root N uptake, NutE and N partitioning (root-N:shoot-N ratio). Additionally, a clear genotype × treatment interaction was observed for the root-N:shoot-N ratio (Table 3). In the control treatment, the average shoot N uptake was 22.0 mg N vessel−1 and varied only by ± 3.18 mg vessel−1 (CV = 14.5%). In the low N treatment, shoot N uptake was on average reduced by 15.3 mg N vessel−1 (69%), and the reduction was significant for all the genotypes tested (Table 4). Amongst the standard varieties, cv. Tomba had the highest N uptake in the control, and many wild potato genotypes achieved shoot N uptakes as high as cv. Tomba. However, under low N supply, the wild potato genotypes generally lag behind cv. Tomba, except for genotype GLKS 30135\_05.

In both treatments, root N uptake accounted for approximately one quarter of the shoot N uptake but showed a stronger variation between the genotypes. In the control, root N uptake was on average 4.59 mg vessel−1 and varied by ± 1.36 mg vessel−1 (CV = 29.6%). Under low N supply, root N uptake decreased on average by 59% to 1.88 mg vessel−1. The reduction was significant for all the genotypes tested, except for GLKS 31600\_10. Again, cv. Tomba was the best cultivar in terms of root N uptake in both treatments, but, in contrast to shoot N uptake, root N uptake of many wild potato genotypes even exceed that of cv. Tomba under control as well as under reduced N conditions (Table 4).

Total N uptake (shoot + root N uptake) generally ranged from 22.1 to 31.4 mg vessel−1 (except GLKS 31600\_10 with only 12.8) in the control, and from 6.90 to 11.3 mg vessel−1 in the reduced N treatment. Relating the total N uptake of the plant to the total amount of N supplied revealed an average NupE of 102% in the control and 132% in the low N treatment. This indicates that the plants took up all the N provided via the nutrient solution and that an additional amount of N was introduced into the system via the ten shoot tips. The NutE, expressed as the amount of total biomass (shoot + root dry mass) produced per unit N taken up, increased, on average, from 25 units in the control to 57 units in the low N treatment. NutE differed moderately between the genotypes in the control treatment where it ranged from 16 to 35 units. In the reduced N treatment, NutE varied strongly from 28 to 78 units. Although all genotypes significantly increased their NutE under low N conditions, many wild potato genotypes exceeded the NutE of the standard varieties. For example, in the low N treatment 13 of the wild potato genotypes achieved a similar and eight genotypes a significantly better NutE than the best cultivar Tomba (Figure 1). The N partitioning, as the ratio of N taken up by the root and N taken up by the shoots (root-N:shoot-N ratio), gives insight into the distribution of N within the plant. In the control, most of the wild types showed values similar to that of the standard varieties, while nine wild types showed significantly higher values than cv. Tomba. Under reduced N conditions, the majority of the wild types increased the root-N:shoot-N ratio, while it remained constant for all the standard varieties. This indicates a stronger partitioning of N to the root for wild potato genotypes, in particular under N deficit.

To evaluate the stress performance of the different genotypes, two commonly used indices (stress susceptibility index SSI, stress tolerance index STI) were applied based on total plant dry mass (shoot DM + root DM). In Figure 2 the accessions and genotypes are ranked according to their SSI and their STI from the best to the weakest genotype. The SSI, as a measure of the yield stability under stress conditions, identified cv. Kiebitz, GLKS 31559\_14 and GLKS 30181\_18 as the three most stable genotypes across both environments, while the reaction to N stress was most pronounced for GLKS 30177\_17, GLKS 30177\_02 and GLKS 31559\_11. On the other hand, the STI identified GLKS 30177\_20, GLKS 30160\_15 and GLKS 30160\_13 as the most promising genotypes, because they produced high total yields under control as well as under stress conditions, whereas cv. Kiebitz and GLKS 31600\_10 ranked least, due to relatively low yields in both treatments (Figure 2).

**Figure 1.** Mean nitrogen utilization efficiency (NutE) of 22 different wild potato genotypes and the four standard varieties under optimal (**a**) and reduced N conditions (**b**) after 21 days of cultivation in a climate chamber. Error bars indicate standard deviation, the dotted line indicates the mean NutE of the best cv. Tomba. The different letters at the bottom of the bars indicate whether there is a significant difference to the best cultivar Tomba ("b") or not ("a", Dunnett's test *p* ≤ 0.05).

**Figure 2.** Stress susceptibility index (SSI) and stress tolerance index (STI) of different wild potato genotypes and four cultivars based on their total dry mass production (DM) under optimal and reduced N supply during 21 days of cultivation in a climate chamber. Genotypes are ranked according to their performance, with the best being on top of the graph.

#### **3. Discussion**

We assessed the tuber starch content in wild potato accessions and studied the reaction of selected wild types to N deficiency. After greenhouse cultivation, the average starch content of all tested wild potato genotypes amounted to 21.7% of FM, but ranged from minimal 7.1 to maximal 36% of FM. This indicates considerable variation for this trait in the wild relatives of cultivated potato. Furthermore, in 50 out of the 506 genotypes, tuber starch content was higher than 28%, which clearly exceeds the starch contents of modern cultivars. Based on our study, germplasm with high tuber starch content is predominately found in *S. chacoense*, *S. commersonii*, *S. jamesii* and *S. pinnatisectum*. Our results are in line with Jansen et al. [25], where tuber starch contents ranged between 3.8 and 39.6% of FM in wild species cultivated in the greenhouse. Here, genotypes of the species *S. chacoense* and *S. pinnatisectum* revealed the highest starch contents. In comparison, 14 modern starch potato varieties cultivated in pot experiments revealed starch contents between 13.9 and 21.9% of FM [20]. Under field conditions, tuber starch contents ranged from 7 to 20% of FM in a set of 300 potato cultivars, breeding clones, landraces and diploid clones [31]. This points out the potential of wild potato germplasm to increase tuber starch contents in cultivars. However, next to the tuber starch content, also the starch yield, as a result of starch content multiplied by tuber yield, plays an important role for industrial starch production. Whether the high starch contents in the wild species will be maintained when tuber size increases due to breeding, has still to be examined. Studies of Schönhals [31] showed that tuber yield is negatively correlated with tuber starch content, while no significant correlation between tuber yield and starch content was found by Bombik et al. [32].

An increase in tuber starch content and starch yield should be linked with a high resource use efficiency. This holds particularly true for the element N, because potato cultivation bears high risks of N leaching due to its high demand of readily available N in soil and the small root system of the crop [5,6,14]. Based on our results mentioned above, 22 wild potato genotypes with high tuber starch contents and four commercial cultivars were studied for their N use efficiency by cultivating them as in vitro plantlets in 500 mL vessels filled with a nutrient solution containing 30 or 7.5 mmol L−1 N for 21 days in a climate chamber. This system allowed us to screen a high number of plants under low space requirement and highly controlled conditions. Several other reports underline the potential of in vitro cultures for the evaluation of potato germplasm with respect to abiotic and biotic stress or rooting characteristics, because it provides conditions independent from weather conditions, pathogens, N leaching or immobilization events [16,33–35]. Our results revealed a high variation in shoot and root DM development between the wild potato genotypes. N deficiency significantly reduced shoot DM for all wild types and the cultivars (except for cv. Kiebitz and cv. Eurobravo), with the extent of shoot DM reduction being genotype dependent. Most interestingly, genotypes GLKS 30135\_05, GLKS\_30160\_13 and GLKS 30177\_20 of *S. chacoense* showed a moderate shoot DM reduction due to low N supply whilst producing a comparably high shoot biomass under high N conditions. The root DM varied considerably and the wild types GLKS 30160\_13, GLKS 30160\_15, GLKS 30177\_15 and GLKS 30177\_20 clearly exceeded the root growth of the best cultivar Tomba under high and/or low N supply, indicating that wild potato germplasm could considerably contribute to enhance root growth of *S. tuberosum* cultivars. Under N deficit, the genotypes followed different strategies in terms of root development. A significant reduction of root DM was observed for the majority of the genotypes, while seven genotypes maintained root biomass. This was also reflected in the root-DM:shoot-DM ratio which was either maintained (13 genotypes), reduced (11 genotypes) or increased (two genotypes). To sustain or even increase root biomass at the expense of the shoots is a well-known reaction of plants to nutrient deficiency and helps to maintain the nutrient uptake from soil or nutrient solution by exploring a larger (soil) volume [36,37]. In contrast, other wild types seem to preferably invest into shoot growth, probably in order to maintain photosynthetic activity. Different strategies to cope with low N as observed in our experiment are also known from cultivated potatoes. A higher root:total mass ratio under N deficiency was reported for the majority of 17 modern starch and table potato varieties during the course of 18 days of in vitro cultivation [16]. On the other hand, the authors also identified some genotypes which missed the ability to stimulate root growth at the expense of the shoots under N deficiency. This is in accordance with previous studies under climate chamber conditions where some cultivars reduced root FM with increasing N stress, while others showed an increased root

development upon N reduction and even maintained root growth at very low N levels [17]. From their studies, Schum et al. [16,17] observed that genotypes with high biomass production and fast nitrogen uptake under high N supply did not enhance root growth under low N and clearly reduced biomass production. On the other hand, genotypes with comparatively slow growth under high N supply increased root mass under low N supply [16,17]. Besides, it should be considered that other well-known responses of plants to nutrient deficiency stress such as changes in the root architecture, increase of root length, root surface area, root volume or number of root hairs [33], is not necessarily reflected in changes of the total root biomass as measured in our experiment. In relation to the total N applied, the plants took up almost all the N available in both treatments (except for GLKS 31600\_10). This explains the low genotypic variation for total N uptake (Table 3). In some cases, N uptake of the genotypes even exceeded the amount of N provided via the nutrient solution, probably due to the additional N introduced into the system via the transferred shoot tips. Therefore, it is difficult to evaluate the NupE of the different accessions. Nevertheless, on the basis of the uniform N uptakes of the genotypes, the results give insight into genotype-dependent N partitioning and provide a clear picture in terms of NutE.

Most of the N taken up by the different genotypes was partitioned to the shoots under high N supply in our study. Under low N supply, a clear shift towards the roots was observed for many genotypes (Table 5). This confirms previous results, where generally a higher percentage of N was translocated to the roots or tubers under nutrient deficiency [6,36]. The increase in root-N:shoot-N ratio became especially evident for many wild potato genotypes, while this was less pronounced for the standard varieties.

NutE was on average 25.3 units in the control and increased to 57.3 units under N deficiency and the genotypic variation in NutE was particularly high in the low N treatment (Figure 1). This is in line with several studies [16,17,20,24]. In our experiment, many wild potato genotypes exceeded the NutE of the standard varieties under low N supply. Outstanding genotypes were GLKS 30154\_09, GLKS 30916\_08, GLKS 30181\_06, GLKS 30156\_16, GLKS 30177\_15, GLKS 30160\_13, GLKS 30160\_15 and GLKS 30177\_20 of *S. chacoense*. This indicates that these genotypes need considerably less N to produce the same amount of biomass. A high NutE is often related to a good translocation of N from the root to the shoot and/or a reallocation from older leaves to the younger leaves in order to maintain the photosynthetic activity and eventually to the reproductive organs [8,36]. The superior performance of *S. chacoense* genotypes was also found in field studies where, amongst 39 wild potato accessions, genotypes of *S. microdontum* and *S. chacoense* revealed the highest total biomass (tubers + roots + shoots + fruits), a high NupE and N recovery from soil, even exceeding the performance of the control varieties cv. Russet Norkotah and cv. Red Norland [22]. The authors attributed the higher N recovery by the wild species to the deeper penetrating, denser, and more branched root system that is advantageous for nutrient uptake. However, it generally has to be considered that wild potato species form only small tubers, in some cases, produce stolons rather than tubers under long day conditions. Whether the high N recovery rate of the wild type will be maintained after crossing to cultivars needs detailed examination. Hybrids of *S. chacoense* and a haploid *S. tuberosum* line (USW551) were studied in the field with high and zero N supply by Errebhi et al. [23]. Here, hybrids showed highest N use efficiency and produced a total biomass (tubers + roots/stolons + shoots + fruits) higher or similar than that of commercial varieties, but tuber yield was low [23].

For a final assessment of the overall performance of the different genotypes under N deficiency, we studied the two stress indices SSI and STI based on the total DM. Zhao et al. [38] studied different indices to evaluate low N tolerance in maize, and advised to use several indices and not to rely on only one. The SSI for example, proposed by Fischer and Maurer [39] for evaluating the yield stability under stressed and non-stressed environments, does not consider the yield of a respective genotype in relation to the other genotypes tested under control conditions [38]. Cv. Kiebitz and GLKS 31559\_14, for example, exhibited the lowest SSI and could thus be considered as the ones with the lowest N stress susceptibility. That is confirmed by no significant changes in shoot or root mass under low N as

compared to high N. However, these genotypes produced a low total biomass during the three-week in vitro culture compared to the other genotypes tested even in the high N treatment. This might indicate that genotypes with a slow biomass development, and by this a probably rather low internal N demand, react less sensitive to a reduction in N supply than fast growing types with a strong biomass development. By calculating the STI, these genotype-specific growth rates were considered, and here genotypes with a high biomass development under both control and stress conditions rank best. Under this premise, cv. Kiebitz and GLKS 31559\_14 were rather intolerant to N stress, while GLKS 30177\_20, GLKS\_30160\_15, GLKS\_30160\_13 and GLKS 30181\_06, GLKS\_30177\_15 and cv. Tomba were more tolerant. Finally, genotypes being among the best under both indices will be very interesting candidates for further research and pre-breeding. Here, we consider GLKS 30177\_20, GLKS 30177\_15 and GLKS 30160\_15 of *S. chacoense* as the most relevant genotypes, because they combine high shoot and root biomasses in both treatments with a moderate reduction in shoot and root biomass under low N supply. Furthermore, the best performers revealed the highest share of root biomass in relation to total biomass, maintained root-DM:shoot-DM ratio under low N, but partitioned more N to the roots than other genotypes and revealed a high internal N utilization efficiency. *S. chacoense* is a well-known source of pest and disease resistance, resistance to cold-induced sweetening and abiotic stresses such as drought tolerance, but its tubers contain high levels of toxic steroidal glycoalkaloids [40,41]. As a diploid species (2n = 2x = 24, EBN 2), hybridization with tetraploid *S. tuberosum* (2n = 4x = 48, EBN 4) is possible after a ploidy reduction in the *S. tuberosum* parent to the diploid level, followed by backcrossing [42], but it is also an interesting future candidate for diploid breeding programs [43].

Apart from a sole comparison of different genotypes, our study highlights the variation of N use efficiency between genotypes within one population. The five genotypes of GLKS 30177 reacted differently to N stress. While GLKS 30177\_15 and GLKS 30177\_20 belong to the best performing genotypes, GLKS 30177\_02 and GLKS 30177\_17 exhibited a medium to low shoot DM, a strong reduction in root DM and root-DM:shoot-DM ratio under low N, low NutE and a low stability under stress (SSI). Instead, GLKS 30177\_01 can be considered as an intermediate type. These results underline the high diversity of the different genotypes within a wild potato accession which is maintained as a population in gene banks. Furthermore, it highlights the importance to study, describe and maintain individual wild potato genotypes in order to promote the use of wild potato germplasm in breeding and research [44]. Therefore, the tested genotypes in this study are maintained clonally via in vitro propagation at the Gross Luesewitz Potato Collections.

Genotypes GLKS 34995\_18, GLKS 31559\_14 and GLKS 31559\_11 of *S. tuberosum* subsp. *andigena* and *S. stenotomum*, belonging to the cultivated part of series Tuberosa and being most related to *S. tuberosum* according to the taxonomy of Hawkes [45], showed no outstanding performance in the respective experiment. Tuber starch contents as well as shoot and root biomass or N uptakes and efficiencies were on an intermediate to low level. Although revealing the highest tuber starch content, the genotype of *S. pinnatisectum* (GLKS 31600\_10) showed a rather weak performance under the given experimental conditions. This was indicated by the lowest shoot and root biomass as well as the lowest N uptakes compared to the other wild potato genotypes. Since this was the sole genotype of this species, we can only speculate whether this is a generally low-yielding species or if the experimental conditions were unfavorable. However, its reactions to N deficit clearly distinguished it from the other genotypes. Although not statistically significant, absolute root mass increased by about 35% under N deficit and the proportion of root biomass in relation to total biomass increased most at the expense of the shoot mass. In relation to its exceptionally high tuber starch contents, it is worth studying the root parameters of this species in further experiments.

#### **4. Materials and Methods**

#### *4.1. Plant Material*