Multi-Environment Yield Components in Advanced Common Bean (Phaseolus vulgaris L.) × Tepary Bean (P. acutifolius A. Gray) Interspecific Lines for Heat and Drought Tolerance

Abstract

Heat and drought are major stresses that significantly reduce seed yield of the common bean (Phaseolus vulgaris L.). In turn, this affects the profitability of the crop in climatic-vulnerable tropical arid regions, which happen to be the poorest and in most need of legume proteins. Therefore, it is imperative to broaden the sources of heat and drought resistance in the common bean by examining closely related species from warmer and drier environments (i.e., Tepary bean, P. acutifolius A. Gray), while harnessing such variation, typically polygenic, throughout advanced interspecific crossing schemes. As part of this study, interspecific congruity backcrosses for high temperature and drought tolerance conditions were characterized across four localities in coastal Colombia. Genotypes with high values of CO₂ assimilation (>24 µmol CO₂ m⁻² s⁻¹), promising yield scores (>19 g/plant), and high seed mineral content (Fe > 100 mg/kg) were identified at the warmest locality, Motilonia. At the driest locality, Caribia, one intercrossed genotype (i.e., 85) and the P. acutifolius G40001 control exhibited sufficient yield for commercial production (17.76 g/plant and 12.76 g/plant, respectively). Meanwhile, at southernmost Turipaná and Carmen de Bolívar localities, two clusters of genotypes exhibited high mean yield scores with 33.31 g/plant and 17.89 g/plant, respectively, and one genotype had an increased Fe content (109.7 mg/kg). Overall, a multi-environment AMMI analysis revealed that genotypes 13, 27, 82, and 84 were environmentally stable with higher yield scores compared to the Tepary control G40001. Ultimately, this study allows us to conclude that advanced common bean × Tepary bean interspecific congruity backcrosses are capable of pyramiding sufficient polygenic tolerance responses for the extreme weather conditions of coastal Colombia, which are likely to worsen due to climate change. Furthermore, some particular recombination events (i.e., genotype 68) show that there may be potential to couple breeding for heat and drought tolerance with Fe mineral biofortification, despite a prevalent trade-off, as a way to fight malnutrition of marginalized communities in tropical regions.

1. Introduction

The common bean (Phaseolus vulgaris L.) is an important food resource and economic income source for farmers around the world [1]. Due to its multiple nutritional components, this species also represents a relevant crop for food security of the poorest. Biofortification efforts have leveraged this potential by offering cultivars with increased Fe content as an alternative to reduce mineral deficiencies [2]. However, an antagonistic effect of the genetic improvement process for drought (and potentially heat) stress tolerance with mineral content has been reported [3], likely jeopardizing long-term intrinsic mineral content properties [4]. In order to minimize this trade-off, germplasm’s screening must reinforce recombinant adaptive phenotypes in concert with high mineral content. Promoting a dual strategy like this promises to diminish the negative impacts of climatic variability on common bean cropping systems, while boosting their mineral content. This would eventually contribute to reduce the production losses in underfed semiarid regions affected by hydric deficit and high temperatures [5].

Drought and heat wave events are becoming more frequent and extreme since the last few decades, particularly in the world’s tropics [6]. In South America, the predictions of temperature increase account for 1.5 °C with several heat waves already recorded [7]. This evidences the need to speed up breeding alternatives to reduce the negative impact of climatic events on agriculture and food security, mainly in regions with strong periods of drought during the year, such as the Caribbean coast of northern South America.

In particular, common bean yield losses due to drought and heat waves in the region may be close to 70% [4]. However, breeding efforts face major bottlenecks to obtain genotypes with good adaptability to adverse abiotic conditions [8] due to the polygenic nature of the tolerant phenotypes. A strategy to develop cultivars with these attributes requires efficiently harnessing genetic resources from other genepools of closely related species typically adapted to heat and drought [9]. The Tepary bean (P. acutifolius A. Gray) is an annual autogamous bean from northwest Mexico [10,11], domesticated near the arid border with the USA, and with a strong preference for hot [12] and dry environments [13,14]. Thanks to this natural adaptation, the Tepary bean is the most heat tolerant among the five cultivated species of the genus Phaseolus. Despite this, the Tepary bean nowadays has limited relevance as a modern crop when compared with the more susceptible, yet commercially accepted, common bean. In order to bridge this gap, we envision using the Tepary bean as donor of alleles to boost drought and heat tolerance in common bean.

Previous works have set the basis to explore the potential of the Tepary bean as an exotic donor by overcoming the phylogenetic distance between the two species, and the naturally low levels of interspecific introgression [15]. For instance, the common bean has been backcrossed with Tepary donors with a relatively good rate of viability [16,17]. Yet, these efforts were unable to pyramid target alleles for heat/drought tolerance, partly due to its complex quantitative genetic inheritance [14]. A promise avenue to overcome this limitation is to rely on advanced interspecific congruity backcross generations [18], as a base population to select for heat tolerance [15]. This is because backcrossing alternately to both parents, as the congruity scheme does, can more efficiently leverage selection of complex polygenic adaptive phenotypes, which is typically more suited for Mendelian traits of simpler genetic architecture [19].

Therefore, the goal of this study was to assess heat and drought tolerance in interspecific congruity backcrosses between common and Tepary beans. First, we aimed quantifying yield components across 86 advanced interspecific lines at four localities in coastal Colombia, where heat and drought waves are frequent and extreme. Second, we envisioned condensing all trials in a multi-locality genotype × environment (G × E) model to examine genotype stability of the best-adapted genotypes. We also addressed the physiological basis of heat stress adaptation in the locality with the most extreme heat events (Motilonia) as a first step to deepen our understanding on the tolerance phenotype.

2. Materials and Methods

2.1. Interspecific Congruity-Backcross Lines

A total of 86 interspecific congruity-backcross lines between the common bean (P. vulgaris) and Tepary bean (P. acutifolius) were considered as part of this study, in addition of the G40001 P. acutifolius control. They correspond to the third generation (and beyond) of congruity backcrosses (detailed pedigree in Table S1). Therefore, it is expected that they partly behave as stable lines with still some potential for further segregation. These genotypes were bred to high temperature and drought conditions by the bean program of the Alliance Bioversity–CIAT (International Center for Tropical Agriculture), and transferred to AGROSAVIA after ATM subscription. All lines were examined to address genotypic differences in phenological and agronomic performance at four localities in coastal Colombia.

2.2. Multi-Locality Field Trials

The study was conducted in four localities of the Colombian Caribbean region corresponding to the following AGROSAVIA’s stations: Motilonia Research Station (Codazzi, province of Cesar), Caribia Research Station (Sevilla, province of Magdalena), Carmen de Bolívar Research Station (Carmen de Bolívar, province of Bolívar), and Turipaná Research Station (Cereté, province of Córdoba). Motilonia and Caribia belong to the dry Caribbean sub-region, respectively being the warmest and driest localities, while Carmen de Bolívar and Córdoba were representatives of the humid Caribbean sub-region. Soil fertility (i.e., in terms of P, Ca, Mg, and K in mg/kg. N and Fe content in the soil are to be included as part of future screenings), and other properties of the study areas are described in Table S2, as recorded during the production cycle of 2020. The bioclimatic variation at these sub-regions offers contrasts in soil type, precipitation, season durability, maximum and minimum temperature difference, and light quality, all of which are known to affect beans’ agronomical performance. Differences between dry and humid Caribbean sub-regions oscillated (1) between 23 °C and 25 °C (average of 23.7 °C) and 33 °C and 36.3 °C (average of 33 °C) for average minimum and maximum temperatures, between (2) 70% and 80% (average of 80%) for relative humidity, and between (3) 482 mm and 700 mm (average of 591 mm) for precipitation during the rainy season, respectively. Hence, environmental conditions at the localities were differential. Table S3 depicts temperature and precipitation variables for each locality during the crop cycle of July–October 2020.

2.3. Experimental Design

The study employed a completely randomized block design (CRBD) with three repetitions at each locality. The experimental unit per treatment (genotypes) was a plot of four m² with a conventional spatial arrangement of one row spaced at 0.8 m and 0.25 m between plants (13 plants per genotype). A total of 72 genotypes germinated and succeed to growth until scoring in Motilonia, 16 in Caribia, 61 in Turipaná, and 87 in Carmen de Bolívar. Missing genotypes were indicative of mal-adaptation at each locality.

2.4. Phenotyping

The following yield traits, standard for the common bean [20], were measured at the end of the cycle at each locality, NS: average number of seeds per pod, NP: number of pods per plant, YLP: yield per plant (g/plant), PB: filled pod weight (g), SB: seed weight (g), VB: vegetative biomass (g), WS: weight of 100 seeds, and DF: days to flowering. Content of minerals Fe and Zn was quantified for plants with a sufficient number of seeds (minimum 100 g) following the Elemental Analysis by Energy Dispersive X-ray fluorescence (EDXRF) protocol of the Nutritional Quality Laboratory at Alliance Bioversity–CIAT (Palmira, Colombia). Due to the COVID-19 contingency, mineral quantification was only possible at Motilonia, Carmen de Bolívar, and Caribia. For Turipaná, a very high pressure of soil pathogens was also observed, so the amount of seed required for the analysis was insufficient.

Detailed ecophysiological measures were possible at Motilonia for the following variables, (1) E: transpiration, (2) A: net photosynthetic rate, (3) Ca: CO₂ (in the chamber), (4) Ci: internal leaf CO₂, (5) gs: stomata conductance, (6) LT: internal leaf temperature, and (7) VPD: vapor pressure deficit. These measures were gathered using the LI-6800 portable photosynthesis system from 8:00–11:00 am. Additional compound parameters were computed as follows, (8) iWUE: instantaneous water use efficiency was measured as the relation A/E, (9) inWUE: intrinsic efficiency of water use computed as the relation A/gs, and (10) iCE: instantaneous carboxylation efficiency was measured as the relation A/Ci [21]. Measurements were made on three plants per genotype on a fully expanded leaf from the upper third part of the plant. Measurements were carried out between the pre-flowering and flowering of the genotypes. The LI-6800 equipment was calibrated for the local conditions through light and CO₂ curves, obtaining the appropriate range of photosynthetically active radiation (PARi) of 1500 µmol of photons m⁻²s⁻¹, and a reference CO₂ of 400 µmol CO₂ mol⁻¹. Calibration curves were done on the P. acutifolius control genotype G40001.

2.5. Statistical Analyses

Analyses of variance (ANOVA) were performed for the examined traits individually at each locality after checking for normality and considering genotypes as factors in the statistical software SAS enterprise Guide 8.3. A hierarchical grouping analysis was also carried out to determine the diversity in the studied genotypes for yield traits using the hclust function with Ward’s agglomeration method in the R software (R Core Team). The Factoextra [22] and NbClust [23] packages were further used to determine the optimum number of clusters (Figure S1).

In order to assess differences among clusters, Tukey’s comparison of means was performed (p-value < 0.05) after checking for normality using R’s library agricolae. Additionally, principal component (PCA) and Pearson’s correlation analyses were carried out using R’s libraries factoextra and FactoMineR [24]. Missing data were imputed when it did not exceed 10% using the mice library [25] from the same software.

Finally, in order to explore the extent of stability and phenotypic adaptability, previous results were further contrasted via the estimation of phenotypic stability for nine overlapping lines among localities within an AMMI analytical framework. The AMMI model used the routine proposed by [26]. It combined the analysis of variance (ANOVA) to estimate the main effects of G and E, with the principal components analysis (PCA) as a proxy for the G × E interaction, as follows (Equation (1)):

Y_ij = μ + G_i + E_j + Σ_λk α_ik δ_jk + R_ij + ɛ

(1)

where Y_ij corresponded to the value of the ith genotype in the jth environment, μ was the overall mean, G_i was the deviation of the ith genotype from the overall mean, E_j was the deviation of the jth environment from the overall mean, λ_k was the singular value for PC axis k, α_ik and δ_jk, respectively, were the principal component’s scores for axis k of the ith genotype and jth environment, and R_ij and ε were the residual and error terms, respectively.

Stable genotypes for each environment were selected by the AMMI biplot, in which scores of the first principal component (CP1) were plotted against the phenotypic means of genotypes and environments. Specifically, CP1 scores close to zero in the y-axis were considered as stable genotypes. Specific genotype–environment adaptability was inferrered from the x-axis, where the means of genotypes and environments were jointly plotted.

This combined analysis was carried out in order to select the best interspecific lines adapted to the target climates at the Caribbean coastal sub-regions of Colombia given two different resolutions: (1) broad adaptation, applicable for genotypes stable across multiple environments, and (2) narrow adaptation, in which specific genotypes outperformed the rest in a single locality but not across all environments.

3. Results

3.1. Phenotypic and Physiological Variability of Interspecific Lines Per Locality

Interspecific lines characterized across four localities (i.e., Motilonia, Caribia, Turipaná, and Carmen de Bolívar) showed a statistical differential response, mainly for parameters associated with yield (

Table 1. One-site analysis of variance for the common bean (P. vulgaris) × Tepary bean (P. acutifolius) interspecific congruity backcross lines.

Yield Trait	Motilonia (72 Genotypes)		Caribia (16 Genotypes)		Turipaná (61 Genotypes)		Carmen (87 Genotypes)
	Rep	Gen	Rep	Gen	Rep	Gen	Rep	Gen
NS	0.4801	<0.0001	0.1	0.52	0.28	0.0177	<0.0001	<0.0001
NP	< 0.0001	<0.0001	0.14	0.38	0.4	0.001	<0.0001	0.005
YLP	< 0.0001	<0.0001	0.04	0.46	0.07	0.04	<0.0001	<0.0001
PB	< 0.0001	<0.0001	0.09	0.38	NA	NA	0.052	0.242
SB	< 0.0001	<0.0001	0.46	NA	NA	NA	<0.0001	<0.0001
VB	0.8546	<0.0001	0.13	0.16	0.03	0.06	<0.0001	<0.0001
WS	NA	NA	0.08	0.4	0.2398	0.42	0.2497	0.001
DF	0.7666	<0.0001	NA	NA	0.08	0.262	0.442	0.122

p-value for repetition (Rep) and genotype (Gen). NS: Average number of seeds per pod, NP: number of pods per plant, YLP: yield per plant (g/plant), PB: pod biomass (g), SB: seed biomass (g), VB: vegetative biomass (g), WS: weight of 100 seeds, and DF: days to flowering.

). These trends are detailed per locality in the next subsections. Since Motilonia and Carmen’s replicates also exhibited differences, uncontrolled variability is discussed in light of local-scale heterogeneity.

3.2. Motilonia (Warmest Dry Caribbean Savanna Sub-Region)

At the Motilonia locality, genotypes had a variable ecophysiological behavior since all ecophysiological parameters shown statistically significant differences, indicating that at least one genotype is different (

Table 2. Analysis of variance for ecophysiological parameters among the common bean × Tepary bean interspecific congruity-backcross lines at Motilonia (warmest dry Caribbean sub-region).

Parameter	Units	Repetitions Effect	Genotype Effect
A	µmol CO2 m−2 s−1	0.8126	<0.0001
Ca	µmol mol−1	0.3461	<0.0001
Ci	µmol CO2 mol−1	0.4298	<0.0001
gs	mol H2O m−2 s−1	0.4597	<0.0001
VPD	kPa	0.2497	<0.0001
E	mmol H2O m−2 s−1	0.9013	<0.0001
Tleaf	°C	0.2881	<0.0001
iWUE	NA	0.6131	<0.001
inWUE	NA	0.184	<0.001
iCE	NA	0.4277	<0.001

p-value for repetition and genotype factors. E: transpiration, A: net photosynthetic rate, Ca: Ambient CO₂ (in the chamber), Ci: internal leaf CO₂, gs: stomatal conductance, E: transpiration, Tleaf: internal leaf temperature, VPD: vapor pressure deficit, iWUE: instantaneous water use efficiency, inWUE: intrinsic efficiency of water use, and iCE: instantaneous carboxylation efficiency.

).

Three clusters were identified according to their productive response (Figure 1a), with cluster 1 being the one with the highest average yield, followed by cluster 2, and finally cluster 3 with the lowest values (Figure 1b). Population variance was mainly explained by dimension 1 (37.5%), 2 (20.4%) and 3 (15.3%). Variables with greater contribution at explaining variability were mostly ecophysiological and performance parameters (Figure 1c, darker vectors). Performance variables like NP, YLP, PB, NS, SB, VB, gs, A, and iCE presented an inverse relationship with variables such as VPD, Ca, LT, and WS.

Pearson correlation analysis corroborated statistically significant correlations among the performance-related parameters (i.e., NS, NP, YLP, PB, SB, VB, and WS), as well as with some ecophysiological parameters (i.e., A, iWUE, gs, and iCE). In short, with higher A, iWUE, gs, and iCE, performance also increased (Figure 1d). In terms of physiological parameters, some correlations were particularly relevant, such as the inverse correlation between E and iWUE (as expected because the E term is the denominator). This result agreed with the ecophysiological behavior of the bean plant since greater transpiration, generated by greater gs, implied having less efficient use of water. Last, a direct correlation between Fe and Zn was found. Zn had a negative correlation with yield parameters (Figure 1d).

Regarding the physiological parameter of greatest interest, A, it was found that some genotypes carried out an efficient assimilation of CO₂ with a lower amount of Ci, involving low VPD and greater water use efficiency. It was also evident through a significant inverse correlation that the higher gs allowed to cool and reduced the temperature of the leaf, implying lower VPD and higher iCE. The inWUE, iWUE, and iCE indices were generally positively correlated (Figure 1).

Regarding differences among genotypes, it was found that the cluster with the highest yield (cluster 1) also presented the highest values for other parameters such as NS, NP, PB, SB, A, gs, and iCE. However, it presented the lowest values for WS and VPD. Similarly, the genotypes grouped in cluster 2 also showed high values for NS, A, gs and iCE, and low Ca. Hence, besides presenting high values in variables related to yield and biomass production, the ecophysiological behavior of cluster 1 also indicated a positive relationship between higher yield and CO₂ assimilation, being an attribute of interest for adaptation in the Caribbean conditions of coastal Colombian (

Table 3. Comparisons of means for yield components and physiological traits among clusters of interspecific congruity-backcross lines in Motilonia (dry Caribe).

Trait	Unit	Cluster 1	Cluster 2	Cluster 3
NS	Seeds	2.57 a	2.16 a	1.73 b
NP	Pods	20.27 a	9.21 b	2.78 c
YLP	g/plant	19.72 a	9.05 b	2.73 c
PB	G	5.18 a	2.49 b	0.76 c
SB	G	21.36 a	9.05 b	2.68 c
VB	G	183.49	154.63	130.13
WS	g/100 seeds	25.46 a,b	23.03 b	28.89 a
DF	Days	36.3	34.96	32.92
E	mmol H2O m−2 s −1	0.012	0.013	0.012
A	µmol CO2 m−2 s−1	24.29 a	23.52 a	20.15 b
Ca	µmol mol−1	358.51 b	358.41 b	364.04 a
Ci	µmol CO2 mol−1	284.55	286.45	289.25
gs	mol H2O m−2 s−1	0.73 a	0.72 a,b	0.60 b
LT	°C	35.66	35.99	36.11
VPD	kPa	1.98 b	2.10 a,b	2.29 a
iWUE	NA	2035.99	1819.60	1698.78
inWUE	NA	35.28	33.76	35.33
iCE	NA	0.08 a	0.08 a	0.07 b

NS: average number of seeds per pod, NP: number of pods per plant, YLP: yield per plant (g/plant), PB: pod biomass (g), SB: seed biomass (g), VB: vegetative biomass (g), WS: weight of 100 seeds, DF: days to flowering, E: transpiration, A: Net photosynthetic rate, Ca: Ambient CO₂ (in the chamber), Ci: internal leaf CO₂, gs: stomatal conductance, LT: internal leaf temperature, VPD: vapor pressure deficit, iWUE: instantaneous water use efficiency, inWUE: intrinsic efficiency of water use and iCE: instantaneous carboxylation efficiency. Different letters denote significant differences determined by Tukey’s test (α ≤ 0.05).

).

On the contrary, cluster 3 presented the lowest values for performance traits and biomass production. In addition, regarding the physiological parameters, it was also evident that the VPD value was the highest. This indicated that, the higher the score of this parameter, the greater the impact on other physiological processes, and therefore the overall productive capacity (Table 3). In addition, cluster 3 also presented the lowest values for the variables A, gs, iWUE and iCE, and the highest for Ca and Ci.

Some of the genotypes that experienced interesting behavior at the warmest locality of Motilonia merged key physiological parameters, outstanding mineral content, and good productive responses (

Table 4. The common bean × Tepary bean interspecific congruity-backcross lines selected according to yield and Fe content in the locality of Motilonia (dry Caribbean sub-region). G40001 was used as control from the P. acutifolius genepool, and therefore is presented in the first row.

Selection Criteria	Cluster	Genotype	YLP	WS	Fe	Zn
	1	G40001	31.43	13	66.8	24.3
	1	82	28.96	28.6	56.2	26.7
	1	48	20.55	50	-	-
	1	17	19.03	23.2	94.35	28.8
	1	81	17.7	27	-	-
	1	5	17.04	23.3	82.1	44.9
	1	16	16.78	23.3	83.1	41.5
	1	3	16.47	22.34	87.8	37.4
Yield	1	18	15.07	24.7	79.2	45.4
	1	27	14.18	19.2	79.5	41.2
	2	1	13.34	25.2	80.0	30.4
	2	84	13.26	24.5	-	-
	2	85	12.14	27	61.7	31.6
	2	9	11.69	19.8	103.75	50.1
	2	6	11.69	27	61.9	25.4
	2	39	10.9	38	-	-
	2	23	10.73	23.8	66.1	34.1
	2	13	10.02	25.6	82.5	31.6
	2	9	11.69	19.8	103.75	50.1
	3	12	3.14	27.5	111.25	40.45
	2	26	6.29	19.9	103.5	44.25
Fe content	3	64	1.96	27	102.6	60.25
	3	76	1.32	25.2	118.3	56.15
	3	77	2.08	29.2	120.5	53.6
	3	80	3.79	23	107.7	52.25

YLP: yield (g/plant), WS: weight of 100 seeds (g), and Fe: iron and Zn: Zinc content (mg/kg).

, despite some cases of inverse relationship between yield and Fe content reinforced this well-known trade-off).

Taking into account that seed size is a paramount commercial trait, some other genotypes were also chosen targeting this variation, specifically: genotype 35 (cluster 2) presented a yield of 5.18 g/plant, a WS of 41.72 g/100 seeds, and Fe and mineral content of 72.30 and 30.60 mg/kg; genotype 38 (cluster 3) had YLP of 3.07 g/plant and WS of 50 g/100 seeds; genotype 45 (cluster 3) had a YLP of 1.72 g/plant and a WS of 46 g/100 seeds; and finally genotype 36 (cluster 3) presented a YLP of 2.08 g/plant and a WS of 43.50 g/100 seeds. These multiple selection criteria also targeted Fe content higher than 100 mg/kg given that some recombinant genotypes may offer potential for biofortification, despite the prevailing inverse relationship between Fe content and yield.

3.3. Caribia (Driest Caribbean Sub-Region)

At this locality two clusters were identified for yield components (Figure 2a). Overall, the genotypes presented low yield, possibly due to the environmental conditions of the region (e.g., high relative humidity that may stimulate the presence of soil pathogens). Only two genotypes, the G40001 control from the P. acutifolius genepool and 85, presented high yield value. Still, other genotypes might also be considered of interest in this locality, such as 84, 86, 28, 13 and 49. However, none of the genotypes exceeded the commercially desirable weight of 40 g/100 seeds. The genotype that came closest to this value was genotype 85 with 38 g/plant, which was also the second genotype with the highest yield. None of these genotypes exceeded the value of 100 mg/kg in Fe content.

When analyzing the relationship between the performance variables, it was found that the variability of the population was mainly explained by the first two components with 49% and 19.2%, respectively. A positive relationship was observed between the variables associated with yield, such as PB, SB, NP, YLP, and NS. However, the contents of Fe and Zn had an inverse behavior against those parameters, and to a lesser extent the variables WS and WB (Figure 2c). Pearson’s correlations made it possible to corroborate these relationships between variables, showing that the yield components had a positive correlation among them. The correlations between Fe and Zn content with some yield parameters were negative, indicating an inverse behavior (Figure 2d).

In the clustering analysis, it was found that cluster 1 was made up of 75% of the genotypes, and was statistically different from cluster 2, which was the one that presented the highest values for all variables, as detailed in

Table 5. Comparisons of means between clusters at locality of Caribia (dry Caribe).

Parameter	Unit	Cluster 1	Cluster 2
NP	Pods	10.29 b	22.37 a
YLP	g/plant	3.17 b	9.53 a
PB	G	9.37 b	24.32 a
SB	G	4.98 b	11.82 a
VB	G	21.78 b	45.28 a

NP: number of pods per plant, YLP: yield per plant (g/plant), PB: pod biomass (g), SB: seed biomass (g), VB: vegetative biomass (g). Different letters denote significant differences determined by Tukey’s test (α ≤ 0.05).

.

Given the selection criteria applied to this locality, only two genotypes presented yield scores that were above the desired threshold: control genotype G40001 from the P. acutifolius genepool (cluster 2) with YLP of 17.76 g/plant, WS of 12.28 g/100 seeds, Fe of 63.80 mg/kg, and Zn of 22.40 mg/kg; and genotype 85 (cluster 2), with YLP of 12.76 g/plant, WS of 38.41 g/100 seeds, Fe of 71.30 mg/kg, and Zn of 35.20 mg/kg.

3.4. Turipaná (Humid Caribbean Sub-Region)

In the Turipaná locality, four clusters were found when considering the main performance parameters, as indicated in Figure 3. The largest cluster was the second, including up to 62.29% of the characterized interspecific genotypes, followed by clusters 4, 3, and 1, with 26.22%, 8.19%, and 3.27%, respectively (Figure 3a). The variability of the population was mostly explained by component 1 (46.7%) and component 2 (19.4%, Figure 3c). The yield (YLP) and its components (NP and NS) were positively related and provided the greatest contribution to explain the variability of the population.

Relationships among performance parameters were corroborated with Pearson’s linear correlations, in which these traits presented direct and significant correlations (p-value < −0.05, Figure 3d). Finally, when comparing the clusters, it was found that the highest yield was obtained within cluster 1, followed by cluster 3, for the variables NS, NP, and YLP. On the contrary, cluster 2 presented the lowest values for these same parameters (

Table 6. Comparisons of means for yield components among clusters of interspecific congruity-backcross lines in Turipaná (humid Caribbean sub-region).

Cluster	Units	Cluster 1	Cluster 3	Cluster 4	Cluster 2
NS	Seeds	6.20 a	4.55 b	3.07 c	2.09 d
NP	Pods	103.17 a	35.73 b	3.07 c	2.09 d
YLP	g/plant	33.31 a	22.59 b	9.22 c	2.81 d

NS: average number of seeds per pod, NP: number of pods per plant, YLP: yield per plant (g/plant). Different letters denote significant differences determined by Tukey’s test (α ≤ 0.05).

).

In this way, the genotypes that presented the highest yield scores, and that belonged to clusters 1, 3, and 4, were chosen and ranked, as shown in

Table 7. Interspecific congruity-backcross lines selected according to yield per plant above 10 g in Turipaná (humid Caribbean sub-region). G40001 was used as control from the P. acutifolius genepool.

Cluster	Genotype	NS	NP	YLP	VB	WS	DF
1	G40001	5.90	106.35	44.06	100.27	11.32	52.67
3	85	4.50	25.50	30.65	91.00	20.10	48.00
3	17	4.20	50.00	23.50	33.67	16.50	52.67
1	13	6.50	100.00	22.57	72.00	17.77	48.00
3	83	4.10	38.50	21.80	43.00	25.76	55.00
3	27	4.90	40.00	20.00	15.00	19.50	44.50
3	82	5.07	24.67	17.00	37.00	21.56	48.00
4	86	2.78	32.33	16.04	37.67	21.49	60.67
4	19	3.11	16.33	13.67	17.67	18.39	62.00
4	34	3.07	31.67	12.70	22.33	17.09	57.33
4	84	2.57	16.00	12.50	11.00	15.00	48.00
4	72	4.30	25.00	11.00	33.00	18.61	55.00
4	75	2.90	27.00	10.67	126.33	17.52	48.00
4	6	3.47	38.33	10.30	42.00	15.39	48.00

NS: average number of seeds per pod, NP: number of pods per plant, YLP: yield per plant (g/plant), VB: Vegetative biomass (g), WS: weight of 100 seeds (g), DF: days to flowering.

.

On the other hand, given the importance of grain weight, genotype 68 presented a good WS of 43.82 g/100 seeds, although in terms of NS had a modest score of 2.90 (besides a NP of 22.50, YLP of 9.00 g/plant, GV of 39.50, and DF of 51.50).

3.5. Carmen de Bolívar (Humid Caribbean Sub-Region)

In Carmen de Bolívar cluster 1 was the most numerous with 44.82% genotypes, cluster 2 with 36.78%, and cluster 3 with 18.39% (Figure 4a). Population variability was mainly explained by component 1 (46.4%) and 2 (17.46%). It was also found that most of the yield components were related, except for DM, DF, WS, Fe and Zn, which presented an opposite response to the rest of the parameters (Figure 4c). This was corroborated with Pearson’s correlations (Figure 4d), in which a high inverse correlation was found between YLP and WS (p-value < −0.01), indicating that the genotypes with the highest yield had low grain weight (i.e., higher yield was obtained at the expense of WS).

In addition, other statistically significant correlations of interest were found, although with lower Pearson’s correlation scores. Specifically, positive correlations between Zn and VB, and between Fe and VB, suggested that plants with higher VB might accumulate higher content of Zn and Fe in the grain.

On the contrary, a low negative correlation was found between WS and Fe. This trend suggested that the grains with the highest weight had the lowest Fe content, as expected due to the prevailing trade-off between yield and nutrition. A high and positive correlation was also found between Fe and Zn.

Cluster 3 presented the highest values for NS, NP, YLP, SB, and VB, and the lowest values for WS and DM. Regarding these latter two traits, cluster 1 was the one with the highest scores (

Table 8. Comparisons of means between clusters of interspecific congruity-backcross lines, as assessed at Carmen de Bolívar (humid Caribbean sub-region).

Cluster	NS	NP	YLP	SB	VB	WS	DM
1	3.46 c	18.84 c	5.98 c	5.99 c	57.60 b	36.18 a	76.11 a
2	4.51 b	23.34 b	11.61 b	11.61 b	67.13 b	30.35 b	72.70 b
3	5.25 a	34.39 a	17.89 a	16.91 a	87.86 a	28.11 b	72.37 b

NS: average number of seeds per pod, NP: number of pods per plant, YLP: yield per plant (g/plant), SB: seed biomass (g), VB: vegetative biomass (g), WS: weight of 100 seeds (g), DM: days to maturation. Different letters denote significant differences as by the Tukey’s test (α ≤ 0.05).

). It was expected that genotypes with the lowest yield could reach the highest grain weight, possibly due to more photo-assimilate accumulation on less grains.

The genotypes selected by yield mainly belonged to clusters 3 and 2, in decreasing order, as indicated in

Table 9. The common bean × Tepary bean interspecific congruity-backcross lines selected for yield, seed weight, and seed Fe in the locality of Carmen de Bolívar (humid Caribbean sub-region). G40001 was used as control from the P. acutifolius genepool.

Criteria	Cluster	Genotype	YLP	WS	Fe	Zn
	3	G40001	27.12	29.47
	3	21	21.12	26	71.5	36.9
	3	86	20.96	27.33	64.4	36.8
	3	82	20.15	28.33	53.1	34.2
	3	84	20.07	28.33	58	30.9
	3	81	19.77	30.33	70.3	44.4
	3	5	19.11	27.67	62	38
	3	72	18.7	28.67	83.6	41
	3	83	17.43	26	67.9	32.9
	3	9	16	26	79.2	45.5
	3	63	15.79	32	92.25	48.4
	3	67	15.76	33	81.8	51.2
	3	7	15	29	77.5	42
	3	18	14.9	24.67	66.4	33.6
	3	25	13.37	24	56.3	31.6
	3	11	11	29	73.6	43.9
	2	43	16.48	35.67	70.2	28.4
	2	62	15.13	31	76.6	42.4
	2	80	14.78	32.67	90.55	41.7
	2	77	14.75	29.42	109.7	58.4
	2	40	14.61	44	66.6	27
Yield	2	79	13.96	26.25	99.9	56.3
	2	29	13.91	40.67	63.1	29.7
	2	17	13.52	24.89	70.3	35.2
	2	85	12.72	28	53.5	35.9
	2	73	12.57	30.33	82.3	46.3
	2	26	12.23	25.33	74.9	35
	2	69	12.03	33.67	92	46.8
	2	41	11.99	48.33	64.4	34.1
	2	16	11.95	28.67	64.6	35.8
	2	13	11.9	26	75.1	45
	2	36	11.81	47.02	73.5	33.5
	2	28	11.7	26.03	56.4	35.4
	2	24	11.35	25	61.1	32.1
	2	65	11.02	29.33	93.15	58.7
	2	23	10.92	25.62	62.6	35.5
	2	19	10.72	25	72.3	40.4
	2	22	10.62	24.33	82.2	39
	2	64	10.34	30	101.1	50.4
	2	3	10.14	27.7	71.4	38.1
	2	71	10.09	30.67	92.65	46.4
	1	45	9.42	41	65.4	39.2
	1	56	9.03	47.87	74	38.4
	1	35	8.24	44	74.3	39.5
	1	47	8.15	46.75	74.3	37.4
WS	1	37	7.36	40.04	61.7	41.9
	1	49	7.09	39.78	74.7	33.1
	1	50	7.02	41.42	96.5	33.9
	1	48	6.76	40.89	65.8	37.9
	1	53	6.58	39.91	57.8	37.5
	1	14	7.68	23.49	101.4	50.2
Fe content	1	78	7.52	25.33	104.4	57
	2	6	8.4	28.66	99.45	41.6
	2	76	7.44	27	106.05	57.9

YLP: yield per plant (g/plant), WS: weight of 100 seeds (g), Fe: iron (mg/kg), Zn: Zinc (mg/kg).

. Of these groups, only genotypes 77 and 64 presented high values for minerals. Some additional genotypes of interest were identified considering weight of seeds.

Finally, mineral content was also accounted as selection criterion for some genotypes in order to advance future biofortification outcomes in the presence of increasing heat and drought constrains.

3.6. Phenotypic Stability

The interspecific congruity-backcross lines with the highest yield and the best stability across localities were identified in a panel of nine genotypes that grew and produced in all four environments. Remaining genotypes were interpreted as mal-adapted.

The selection criterion for phenotypic stability was based on distance from the origin to the first component (ASV) in the standard AMMI analysis, which meant that the genotypes with the lowest stability index value (YSI) were discarded (

Table 10. Stability index with AMMI statistics for yield per plant in nine common bean × Tepary bean interspecific congruity-backcross lines evaluated across four localities in the Caribbean coast of Colombia. G40001 (P. acutifolius) was used as control.

Genotypes	ASV	YSI	rASV	rYSI	Square Root	Yield (g/plant)
13	1.08	10	5	5	3.26	12.14
27	0.94	10	4	6	3.17	11.54
28	1.1	14	6	8	2.61	7.43
68	0.80	12	3	9	2.10	5.62
82	1.43	11	8	3	3.83	17.32
84	0.74	6	2	4	3.51	13.22
85	1.69	11	9	2	3.88	17.06
G40001	0.12	2	1	1	5.13	28.16
9	1.36	14	7	7	2.76	8.89

ASV: AMMI stability value, YSI: stability index value, rASV: ranking ASV, rYSI: ranking.

, Figure 5).

The AMMI analysis showed that the first component explained 59.3% of the variance due to G × E interaction for total yield, while the second component explained 33.3%. Therefore, the selection of a single multiplicative term of the AMMI model was sufficient to explain a large proportion of the relevant data [27]. The genotypes closest to the origin point in the AMMI biplot were those with little contribution to the interaction effects, and therefore more stable. Among those, genotypes 13, 27, 82, and 84 were more stable (i.e., broad adaptation) with higher yields, compared to the G40001 P. acutifolius control. Genotypes 68, 84, 85, and the P. acutifolius G40001 control, respectively located closer to localities Caribia, Motilonia, Turipaná and Carmen de Bolívar, suggesting narrow adaptation [28] to such specific environments in the dry and humid Caribbean sub-regions of Colombia.

Heat and drought events were likely more extreme in the former (i.e., Motilonia, Turipaná), while severe drought was presumably preponderant in the latter (i.e., Turipaná and Carmen de Bolívar). In short, this AMMI approach allowed selecting four stable genotypes across environments, and four additional genotypes with narrow adaptations to specific sites [29,30,31,32]. Furthermore, the AMMI analysis made possible a straightforward interpretation of the results using a biplot. This enabled improving the precision of the G × E interaction, since it eliminated the error around phenotypic stability parameters generated by the effect of some environments on particular genotypes [30].

4. Discussion

Across all four localities, the P. acutifolius control genotype G40001 presents the best performance, evidencing its ability to adapt to the different conditions in coastal Colombia. This behavior is consistent with extensive reports that mention the potential of the species P. acutifolius to grow and produce in warm arid and semi-arid conditions [13,14]. In addition, this study reveals that genotypes from advanced interspecific congruity backcrosses with this species exhibit promising responses to extreme conditions in Coastal Colombia, suggesting the inheritance of polygenic adaptation across interspecific boundaries. Yet, selection of elite genotypes may also be guided at each locality with the aim to identify genotypes with particular narrow adaptability.

Looking back into the pedigree of the elite interspecific lines, key donors include the SER, SEF, HTA, GGR, and SMG genotypes (i.e., general combining ability due to the additive parental breeding values), as well as the progenies from the crosses among G 40119, G 40264, G 40036, DAA 9, SAR, SAB, SEF, RAB, SMR 180, ALB 74, INB 841, SMC 199, and DAB 295 (i.e., specific combining ability due to relative non-additive performance of particular crosses). Hence, genotypes resulting from the Alliance Bioversity–CIAT bean breeding program show potential as parental donors for interspecific crosses [33,34,35,36], or even as cultivars themselves, in the face of heat and drought stresses in coastal Colombia [37]. In this way, this works offers novel perspectives to harness, via advanced interspecific schemes, the highly contrasting diversity present at beans’ germplasm, while bridging major bottlenecks to breed traits of interest in the face of climate change.

Currently, the Tepary bean represents an important genetic resource to improve resistance to heat and drought stresses, which are the most limiting factors in common bean [12]. Throughout this research, we have demonstrated the potential of this species to contribute outstanding trait scores via interspecific congruity backcrosses with the common bean. This approach offers new scopes to utilize the conserved bean variation to strength the sustainability of the common bean crop in hot semiarid and arid conditions in tropical regions [38]. Such hybrid breeding strategy is necessary to leverage the over-dominance of genetic variation previously locked within species because the genomic architecture [9,39] conferring thermotolerance in the Tepary bean [40] likely differs from the one in the common bean [41]. Specifically, GWAS (Genome Wide Association Studies) and candidate gene approaches have unveiled complex genomic architectures for drought tolerance, subjected to rampant divergent selection and limited cross-species standing variation [42,43,44,45]. Similarly, the genetic basis of heat tolerance, as inferred in common bean from GEA (Genome–Environment Association) scans and candidate genes [46,47], only holds partially in the Tepary bean [9]. In short, the Tepary bean could help bridging long-standing gaps in the breeding of common bean for drought and heat tolerance, due to the interspecies mosaic in genetic determination mechanisms for these abiotic stresses. Specifically, the Tepary bean may offer promising [47], accessible, and novel genetic pathway sources of resistance [48], targeting tropical geographical regions with prolonged periods of intense heat and drought events [49].

Additionally, our research also invites consideration of the feasibility to simultaneously breeding for drought tolerance and adaptation to low soil fertility, such as via line 68 (a.k.a., SMG8). This interspecific backcross recombinant is particularly outstanding because it breaks a prevalent inverse relationship between Fe content and yield, pervasively observed across the common bean × Tepary bean congruity backcross families. The fact that some interspecific recombination events break intrinsic trait correlations evidences the potential to mitigate phenotypic trade-offs by taking advantage of cross-species variability. In order to achieve a more integrative view on the interaction among various trait types, oncoming trials must select promissory interspecific lines after considering further localities in the Caribbean region of Colombia experiencing combined restrictive environments (e.g., high water limitation and low soil quality).

Moreover, this study reinforces the potential of interspecific congruity backcrosses between exotic donors of adaptation and more susceptible elite genepools. Congruity schemes perform multiple alternating backcrosses [8,15,18,50], allowing for a better retention of the exotic phenotype while overcoming the hybrid sterility, genotype incompatibility, and embryo abortion typically found in simpler backcrosses [9]. As shown here, congruity backcross pedigrees behave equally well for the introgressive breeding of polygenic quantitative adaptation [51] (i.e., heat/drought stress tolerance [17]). Therefore, this scheme brings a realistic pathway to efficiently unlock polygenic adaptation, and effectively utilize hidden variation from tolerant landraces, related crop species, and well-adapted crop wild genepools [9].

Yet, we must recognize the pilot nature of the present study when it comes to the physiological basis of the heat and drought tolerance phenotypes. Even though we carried out ambitious enviromics’ multi-locality trials of interspecific crosses, typically regarded as of low viability, more efforts are required to bridge the gap between yield and physiological components during stress responses. In addition, unexplained variation among replicates was evidenced for Motilonia and Carmen de Bolívar in the single-locality analysis. This is possibly due to some intrinsic factors such as intra-family segregation, spatial variability of the soil, or variation of the edaphic microbiota due to high temperature and low organic matter content (i.e., below 2%, which has been reported as a condition that can affect the site-specific response of bean yield [52]).

Local-scale spatial variability is undeniably regarded as a major driving force in various cropping systems [53,54,55]. Fortunately, the effect of soil heterogeneity in the statistical analysis of bean trials is starting to be recognized by the workflow that the Alliance Bioversity–CIAT has carried out as part of the bean genetic improvement program. Specifically, a novel platform called Mr. Bean is undergoing testing in order to better handle local-scale spatial information [56], and will certainly prove useful for future trials. Meanwhile, as part of future projects, we aim extending time-seriated physiological screening from the warmest Motilonia locality to more mild climates in coastal Colombia such as Turipaná and Carmen de Bolívar, where heat is predicted to become more intense and frequent in the oncoming years. Beyond this, it would also be desirable to better control for the effect of specific microbiological activity on the performance parameters. Envisioning novel spatial analysis and experimental designs would help minimizing the impact of uncontrolled intra-genotype variability, even at the G × E interface. These trans-disciplinary approaches are ultimately needed to gather a more mechanistic understating of the tolerance response across the environmental continuum, specifically concerning traits’ trade-offs.

5. Perspectives

After having assessed adaptation in advanced common bean × Tepary bean interspecific congruity backcrosses through four localities in coastal northwest South America, the study of G × E interaction would require to better account for the inter-annual variability. We aim addressing this perspective in the years to come by validating and extending in time the key predictions from the present study. Future field trials will have to be explicitly contrasted with the current analysis of the 2020 production cycle, which was typically heat and warm for coastal Colombia (therefore serving as a standard condition to capture adaptive variation for heat and drought).

Furthermore, the use of parental Tepary beans may also confer resistance to soil pathogens as part of future introgressive breeding schemes. This expectation is due to the ability of common bean × Tepary bean interspecific lines to survive in localities with presumably high incidence of biotic stresses, such as Caribia. After all, some common bean parental genotypes are derived from selection schemes carried out by the Alliance Bioversity–CIAT at the locality of Santander de Quilichao (province of Cauca), where the pressure of soil pathogens is high and similar to the one observed in some regions of coastal Colombia.

Other interspecific lines examine here might also exhibit future potential to introgress resistance to bean common mosaic virus (BCMV). This is the case of lines 84 and 86 (a.k.a., SEF 16 and SEF 54), both with red seed coat color, as well as the pure 118 line (a.k.a., SER118), developed by the breeding program of Escuela Agrícola Panamericana Zamorano (EAP) in Honduras. Yet, these expectations remain speculative at the studied localities in Coastal Colombia. Hence, future studies must carry out detailed characterizations of in situ biotic pressures at these localities, as a key step to deepen our understanding on the role of hybrid breeding for biotic stress resistance in bean species. These studies must be envisioned within the framework of the plant disease triangle (PDT), which postulates that any plant disease is the result of the interaction between a host’s genotype, the biotic stress, and their unique environmental combination.

Another necessary next step will involve exploring the underlying genomic basis of drought and heat tolerance across the species boundary via high-throughput sequencing [9]. First, genomic diversity and signatures of selection for heat and drought tolerance may be explored by tracing back changes in the allelic and haplotype frequencies across the available interspecific congruity-backcross genealogies that made part of this study [15]. Theoretical expectations are that selection signatures of introgressive breeding will be enriched in (i) key underlying regions of the tolerance adaptive phenotype [18], as well as on (ii) low recombining regions in high linkage disequilibrium that are more prone to lineage sorting (i.e., genetic drift) [57] and linked selection due to smaller effective population size [58,59].

Second, more explicit additive models such as GWAS mapping [60], and polygenetic genomic prediction (GP) will respectively enlighten the genomic architecture of drought and heat tolerance in Phaseolus beans, while improving predictions of the genomic estimated breeding values (GEBVs) in yet to be established congruity-backcross hybrid seedlings (i.e., genomic-assisted backcrossing–GABC) [60]. Given the complex nature of the heat and drought tolerances, it is foreseen that the GP and GABC approaches would outperform GWAS-type modeling [61], without meaning that the latter will be incapable to capture variants with moderate effects segregating at medium frequencies, still useful to boost more traditional and scalable marker-assisted backcrossing (MABC) initiatives [62]. Merging these approaches will eventually insight into molecular evolution and pre-breeding under heat and drought stresses [12], both at the yield and ecophysiological levels [63].

6. Conclusions

This research explored the yield, ecophysiology, and mineral content of Tepary and common bean interspecific congruity-backcross lines in four localities with contrasting conditions at the Caribbean region of Colombia. Mineral content of the grain in terms of Fe and Zn was carried out completely in three localities, while ecophysiological responses were recorded in one location as a pilot proxy to intermingle yield and physiological responses at extreme heat and drought conditions. Genotypes obtained by interspecific congruity backcrosses emphasize the importance of conserving and utilizing agrobiodiversity as part of the pre-breeding for adaptation to heat and drought conditions. Inter-specific backcrosses 13, 27, 82, and 84 showed desirable adaptation, yield scores, and stability across all four localities, close to the Tepary bean’s G40001 control genotype. Yet, selection for narrow adaptation is also feasible because some interspecific lines exhibited good responses in terms of yield and mineral content per locality.

In the same way, as part of the characterization of ecophysiological parameters, it was possible to identify how the cluster of genotypes with the best performance also presents the highest scores for assimilation of CO₂ and WUE. Therefore, parameters related to gas exchange may be used as indirect selection indices (i.e., proxy traits to speed-up breeding [64,65]). In order to unveil this potential, we recommend continuing with high-throughput phenotypic characterizations of target yield, ecophysiological and mineral content parameters in an extended panel of interspecific congruity-backcross lines.

Last but not least, expanded field trials across contrasting agroecological regions must procure including localities beyond dry and humid coastal northwest South America, as part of an ambitious enviromic [66] approach [67]. A multi-climate enviromic perspective offers opportunities to better optimize G × E arrangements for yield components and biofortification initiatives in the face of novel antagonistic abiotic and biotic interactions unleashed by climate change. This initiative will ultimately enable selecting for heat/drought tolerant [68,69] parents and new recombinants [70], hopefully by breaking the long-standing genetic trade-off between yield and Fe biofortification.