Phenotypic Diversity and Multivariate Characterization of Opuntia ficus-indica from the Inter-Andean Dry Valleys of Northern Ecuador

1. Introduction

The genus Opuntia Mill. (Cactaceae) comprises one of the most ecologically and economically significant groups of succulent plants in the Americas. Members of the family Cactaceae are highly specialized dicotyledonous angiosperms adapted to arid and semiarid environments through remarkable physiological strategies, including nocturnal CO₂ fixation via crassulacean acid metabolism (CAM), stem succulence, and efficient water-storage tissues [1,2]. Within this family, Opuntia ficus-indica (L.) Mill., commonly known as prickly pear or cactus pear ("tuna" in Andean Spanish), stands out for its wide geographic distribution, cultural relevance, and its numerous traditional and modern applications as food, fodder, medicine, and raw material for agro-industrial processing [1,3]. The fruit (prickly pear) is commercially exploited in over 30 countries across Latin America, North Africa, the Mediterranean Basin, and South Asia, with global cultivated area concentrated in Mexico (50,000–70,000 ha), Italy, Tunisia, and Morocco [4,5]. This breadth makes O. ficus-indica one of the most globally important xerophytic crops and a flagship species for climate-resilient agriculture under current and projected climate change scenarios [4,5].

Molecular and biogeographic evidence indicates that O. ficus-indica was domesticated from wild Opuntia populations in central Mexico, where the greatest interspecific and intraspecific genetic diversity has been documented [6,7]. Following European colonization, the species was introduced into the Mediterranean Basin during the late 15th and early 16th centuries, where it naturalized rapidly in Sicily, southern Italy, and North Africa owing to favorable climatic conditions [1,8]. Today, O. ficus-indica is deeply embedded in traditional agroecosystems across four continents, contributing to household economies, erosion control, land rehabilitation, and ecosystem services including wildlife habitat provision and carbon sequestration [4,5,8]. In the Andean highlands of South America—particularly in Peru, Ecuador, and Bolivia—the species occupies a central role in subsistence agriculture and local markets, thriving in the dry inter-Andean valleys at altitudes ranging from 1,300 to over 2,600 m above sea level [9].

Beyond its economic significance, O. ficus-indica plays a crucial ecological role in dryland ecosystems, underpinned by a suite of morpho-physiological adaptations. Its broad cladodes with thick cuticles, glochid-bearing areoles, and water-storage parenchyma confer adaptive advantages under extreme water limitation [1,2]. As a CAM plant, O. ficus-indica fixes atmospheric CO₂ at nighttime while keeping stomata closed during the day, dramatically reducing evapotranspiration and increasing water-use efficiency (WUE) compared to C₃ and C₄ species [2,10]. This physiological strategy is a key determinant of the species’ capacity to colonize degraded soils, tolerate prolonged drought, and function as a nurse plant facilitating the establishment of other dryland vegetation in the hypervariable precipitation regimes characteristic of Andean dry valleys [4,5,10].

Despite its ecological and agronomic importance, the taxonomy and morphological classification of Opuntia remain challenging. The genus is characterized by extensive phenotypic plasticity, frequent interspecific hybridization, polyploidy (ploidy levels ranging from diploid 2n = 2x = 22 to octoploid 2n = 8x = 88), and high sensitivity to environmental conditions, which collectively complicate species delimitation and cultivar identification [11,12]. Molecular phylogenetic studies using plastid (atpB-rbcL, matK) and nuclear markers (ITS, ppc) have revealed evidence of reticulate evolution and incomplete lineage sorting, supporting the notion that morphology alone may be insufficient to resolve species boundaries in certain lineages [11]. Nonetheless, morphological descriptors remain indispensable for germplasm documentation, phenotypic characterization, and cultivar differentiation—particularly in regions where molecular facilities are limited or where traditional management systems rely on visual recognition of morphotypes [3,13,14].

The use of standardized morphological descriptor systems—specifically the FAO/CACTUSNET descriptors for Opuntia spp. [13] and the UPOV guidelines for cactus pear [14]—provides a common framework for comparing germplasm across geographic regions and research institutions. Studies employing these protocols in Tunisia [15], Morocco [16], Ethiopia [17], and Algeria [18] have consistently demonstrated that multivariate analytical approaches—particularly principal component analysis (PCA), multiple correspondence analysis (MCA), and hierarchical cluster analysis—are powerful tools for identifying morphotypic groups, detecting discriminant traits, and revealing patterns of phenotypic diversity associated with environmental gradients, management history, and geographic origin [15,16,17,18,19].

In northern Ecuador, O. ficus-indica is widely distributed across a management continuum ranging from cultivated forms under active agronomic management to semi-wild forms (plants established on degraded hillsides or field margins without systematic agronomic care, propagated by vegetative dispersal or abandoned cultivation) and naturalized populations with no apparent human management, across the dry inter-Andean valleys of the provinces of Imbabura, Carchi, and Pichincha. These landscapes encompass pronounced altitudinal gradients (approximately 1,300–2,800 m a.s.l.), annual precipitation ranging from 350 to 800 mm, and high solar radiation—conditions representative of Andean xerophytic ecosystems [9]. Such environmental heterogeneity may influence the phenotypic expression of key morphological traits including cladode dimensions, spine density, fruit pigmentation, and reproductive morphology, as has been reported for altitudinal gradients in other dryland regions [17,20]. Despite the ecological and socioeconomic relevance of the species in this region, systematic studies documenting the morphological diversity of O. ficus-indica in Ecuador remain scarce, geographically fragmented, and largely absent from peer-reviewed literature. This gap limits the development of conservation strategies, genetic improvement programs, and the identification of elite accessions for agricultural or agro-industrial valorization [4,5,21].

To address these gaps, comprehensive morphological characterization is necessary to understand how environmental heterogeneity, local management practices, and potential adaptive differentiation shape phenotypic variation in O. ficus-indica accessions from Andean drylands. Standardized morphological descriptors—combining qualitative traits (e.g., fruit pigmentation, cladode coloration, petal color) with quantitative traits (e.g., cladode dimensions, fruit weight, seed size)—provide a robust framework for evaluating diversity patterns and identifying key discriminant characteristics across germplasm collections [13,14,19]. Such multivariate analyses, when combined with hierarchical clustering approaches, support the detection of unique phenotypes valuable for conservation and the selection of high-performing accessions for breeding programs targeting improved fruit quality, reduced spinescence, or enhanced drought tolerance [4,5,15,22,23].

Therefore, the present study aimed to characterize the morphological diversity of Opuntia ficus-indica accessions collected from the dry inter-Andean valleys of northern Ecuador through the application of 20 standardized qualitative and quantitative morphological descriptors following FAO/CACTUSNET [13] and UPOV [14] protocols, combined with multivariate analytical approaches (PCA, MCA, and Ward–Gower hierarchical clustering). We hypothesized that the pronounced environmental gradients in altitude, temperature, and precipitation regime of this region drive significant phenotypic differentiation among accessions, particularly in traits associated with fruit morphology, cladode structure, and spinescence [6,7,24,25]. The findings of this research contribute to filling a critical knowledge gap regarding Andean cactus pear germplasm, providing a scientific basis for conservation, sustainable management, and future breeding and agro-industrial initiatives [4,5,21,23].

2. Materials and Methods

2.1. Plant Material

A total of 55 accessions of Opuntia ficus-indica (L.) Mill. were evaluated. These were collected between 2025 and 2026 from three provinces located in the dry inter-Andean valleys of northern Ecuador: Carchi (25 accessions), Imbabura (23 accessions), and Pichincha (7 accessions). Collection sites spanned an altitudinal gradient from 1,381 to 2,758 m a.s.l., with annual precipitation regimes of 350–800 mm and high solar radiation—characteristics representative of Andean xerophytic ecosystems (Figure 1). All accessions were georeferenced using GPS coordinates. Adult plants were selected in situ, prioritizing healthy individuals with no visible symptoms of pests or diseases and with no recent history of agronomic treatments that could alter phenotypic expression. In populations with more than five individuals, three adult plants per accession were randomly selected. All collection activities were carried out on privately owned farmland and cultivated sites with the explicit consent of the landowners. No collection was performed in protected natural areas or wild populations; therefore, no formal biodiversity access permit was required under Ecuador’s Código Orgánico del Ambiente (2017).

2.2. Sampling Strategy and Morphological Evaluation

Morphological characterization was performed following standardized international protocols adapted to the species, prioritizing descriptors with high discriminatory power. A total of 20 morphological descriptors were evaluated—five qualitative and fifteen quantitative—distributed across four plant structures: cladodes (5 descriptors), fruits (9), seeds (4), and flowers (2) (Table 1). Descriptors were drawn from the UPOV guidelines for Opuntia spp. [14] and the FAO/CACTUSNET Minimum List of Descriptors [13], complemented by standardized quantitative measurements established for this species.

Qualitative characters were recorded using ordinal scales and color codes from the Royal Horticultural Society (RHS) color chart, using direct visual comparison with the RHS color table and its RGB equivalents for greater precision and reproducibility. Fruit shape was classified into five nominal categories: 1 = spherical, 2 = elliptical, 3 = lanceolate, 4 = oval, and 5 = rectangular. Quantitative characters were measured using a vernier caliper (precision 0.01 mm), analytical balance (precision 0.01 g), and stereoscopic magnifier. A minimum of 10 mature cladodes, 10 fruits at commercial maturity, and 10 flowers per accession were evaluated from all available representative plants (10–15 individuals in populations exceeding five plants, or all individuals in smaller populations), following UPOV [14] sampling recommendations for uniformity and representativeness.

2.3. Statistical Analysis

2.3.1. Descriptive Statistics and Variable Selection

Descriptive statistics—mean, standard deviation (SD), coefficient of variation (CV), and range (minimum–maximum)—were calculated for the 15 quantitative variables evaluated. Variable selection for multivariate analysis followed three hierarchical criteria: (i) variability: variables with CV ≥ 15% were prioritized, discarding those with low discriminatory capacity; (ii) non-redundancy: Pearson correlations between pairs of variables were evaluated, eliminating one variable from each pair with |r| ≥ 0.70 (fruit length was discarded due to collinearity with fruit weight; seed length was discarded due to collinearity with seed diameter); and (iii) biological relevance: floral scar depth was excluded due to its lowest individual sampling adequacy index (MSA = 0.251; see Section 2.3.3), and flower length and floral scar diameter were excluded for their low individual sampling adequacy (MSA < 0.40) and limited contribution to group separation. Additionally, fruit diameter was excluded due to collinearity with fruit weight (r = 0.807; p < 0.001). The final set comprised 10 representative variables covering vegetative growth, defense structures, fruit characteristics, seed attributes, and floral morphology.

2.3.2. Correction for Multiple Comparisons in Pearson Correlations

Given that 105 pairwise Pearson correlations were simultaneously calculated (15 quantitative variables), p-values were corrected for multiple comparisons using the Bonferroni method (α* = 0.05/105 = 0.000476) and the Benjamini–Hochberg false discovery rate procedure (BH–FDR, α = 0.05). Only associations significant under both methods were considered robust. The exception was the correlation between 100-seed weight and flower length (r = −0.391), significant under BH–FDR (p~BH~ = 0.015) but not under Bonferroni (p~Bonf~ = 0.329), which was interpreted with caution and not used as a variable exclusion criterion.

2.3.3. Validation of the Correlation Matrix for PCA

Prior to Principal Component Analysis (PCA), the adequacy of the correlation matrix was evaluated on the final set of 10 selected variables using Bartlett’s test of sphericity and the Kaiser–Meyer–Olkin (KMO) index. Bartlett’s test contrasts the null hypothesis that the correlation matrix equals the identity matrix; its rejection (p < 0.05) indicates sufficient shared variance to justify PCA. The test yielded χ²(45) = 142.78 (p < 0.001), confirming matrix adequacy. The overall KMO value was 0.59, marginally below the conventional threshold of 0.60 [27] but within the acceptable range (≥ 0.50) established for exploratory factor analyses of morphological data [27,28]; this value was considered adequate given the statistical significance of Bartlett’s test and the biological interpretability of the resulting components.

Analysis of individual sampling adequacy (MSA) values revealed that floral scar depth had the lowest MSA in the matrix (0.251), providing an additional statistical justification—complementary to biological criteria—for its exclusion from the PCA variable set. Cladode diameter also presented a low MSA (0.303); however, it was retained due to its recognized relevance as a primary vegetative descriptor in Opuntia characterization studies [14,27] and its differential contribution to the separation of morphotypic groups observed in the results (see Section 3.3). In exploratory germplasm analyses with mixed variables, retention of variables with MSA < 0.50 is justified when there is a solid biological rationale and the resulting components are interpretable [28,29].

2.3.4. Principal Component Analysis (PCA)

PCA was performed on the correlation matrix of the 10 selected variables, following data standardization (mean = 0, variance = 1), with the objective of identifying the principal axes of phenotypic variation and describing the multivariate structure of the accessions. Interpretation was based on eigenvalues (λ), variance explained per component, and variable loadings on each axis. Components with λ > 1.0 were retained (Kaiser criterion [26]). PCA was performed in InfoStat [30].

2.3.5. Multiple Correspondence Analysis (MCA)

The five qualitative characters (cladode color, fruit shape, peel color, pulp color, and petal color) were analyzed by MCA to explore associations between categories and the multivariate distribution of accessions in qualitative space. MCA was implemented in InfoStat [30]. The significance of each axis was assessed using the χ² test associated with the corresponding eigenvalue. The low cumulative inertia percentage explained by the first two axes (8.19%) is consistent with expected MCA behavior when applied to variables with a high number of categories [31], since total inertia scales with the number of active categories rather than the number of variables [29]. Under these conditions, inertia values below 15–20% in the first two axes are common in germplasm characterization studies using RHS polychromatic descriptors [13,18], and the interpretive relevance of axes should be evaluated based on the biological coherence of the revealed groupings, not the inertia percentage [29,31].

2.3.6. Hierarchical Clustering and Validation of the Optimal Number of Groups

Hierarchical clustering was performed using Ward’s method with Gower distances to simultaneously integrate quantitative and qualitative variables into a unified dissimilarity measure [32]. The Gower distance normalizes each quantitative variable by its observed range and assigns binary dissimilarity (0 = same category, 1 = different category) for qualitative variables, producing values in the interval [0,1]. The resulting distance matrix presented a mean dissimilarity of 0.412 (SD = 0.071; CV = 17.3%; range: 0.186–0.641), with qualitative variables contributing 54.8% of the total mean distance despite representing only 25% of the descriptors. This higher proportional contribution reflects the high categorical richness of RHS chromatic descriptors, whose binary coding in the Gower distance amplifies their relative weight compared to continuous quantitative variables; this behavior justifies the use of Gower distance over exclusively Euclidean metrics [32]. Dendrogram representation quality was assessed using the cophenetic correlation coefficient (r = 0.63; p < 0.001) [33].

The optimal number of groups was determined using two complementary criteria. As the primary criterion, the sequence of fusion distances in the dendrogram was inspected to identify natural gaps between successive fusions [34]: a large gap indicates a real increase in heterogeneity when moving from k to k+1 groups. As supplementary criteria, the Silhouette index [35] and Calinski–Harabasz index [36] were calculated for k = 2–7, acknowledging that their calculation on standardized quantitative variables with Euclidean distance may underestimate the separation achieved by the original Gower distance on mixed data. All supplementary validation calculations were performed in R [37] using the cluster [38] and factoextra [39] packages.

2.3.7. Statistical Validation of Differentiation Between Morphotypes

To statistically confirm morphological differences between the three groups identified by hierarchical clustering, the Kruskal–Wallis test was applied to each quantitative variable, followed by the Dunn post hoc test [40] with Bonferroni correction for pairwise comparisons (α = 0.05). This non-parametric approach was selected due to the non-normal distribution of several morphological variables (Shapiro–Wilk test [41], p < 0.05) and the unequal group sizes (G1 n = 8; G2 n = 26; G3 n = 21). Groups sharing the same superscript letter within a row do not differ significantly.

2.3.8. Independence Tests for Qualitative Characters

Statistical association between qualitative characters and membership in morphotypic groups was evaluated using two complementary approaches. For chromatic characters coded with the RHS scale (cladode color, peel color, pulp color, and petal color), the Kruskal–Wallis test was applied to the ordinal numerical codes, given the high number of categories (15–30) that precludes a valid χ² test due to contingency table sparsity. For fruit shape (five nominal categories), Pearson’s χ² independence test was applied. In all cases, association strength was quantified using Cramér’s V [42], calculated on the complete contingency table—a valid metric regardless of expected frequencies [43]. As a confirmatory check, chromatic variables were collapsed into three tonal groups (pale: RHS 1–7; intermediate: RHS 8–15; intense: RHS 16–30) and subjected to χ² testing. The significance level was α = 0.05.

2.3.9. Phenotypic Diversity Indices

Phenotypic diversity of qualitative characters within each morphotypic group was quantified using the Shannon index [44] (H’ = −Σ p~i~ ln p~i~), where p~i~ is the relative frequency of each category. Pielou’s evenness [45] (J’ = H’/ln S, where S is the observed category richness) was calculated as a complementary measure of distributional uniformity. Both indices were calculated independently for each qualitative character and each group (G1, G2, G3).

2.3.10. Statistical Software

Descriptive, correlation, PCA, and MCA analyses were performed in InfoStat [30]. Clustering validation calculations (Silhouette index [35], Calinski–Harabasz index [36], dendrogram gap analysis [34]), and Gower distance [32] were performed in R [37] using the cluster [38] and factoextra [39] packages. Kruskal–Wallis and Dunn–Bonferroni tests [40] were performed in R [37] using the dunn.test package [46]. The Shapiro–Wilk test [41] was performed in R [37] using the base function shapiro.test. χ² tests, Cramér’s V [42,43], Shannon [44] and Pielou [45] indices were calculated in R [37] using the rcompanion [47] and vegan [48] packages, respectively. Tabular data were read and organized in Microsoft Excel.

3. Results

3.1. Descriptive Statistics and Phenotypic Variability of Quantitative Morphological Traits

Descriptive statistics (Table 2) revealed considerable variation among the morphological traits evaluated, as indicated by the coefficient of variation (CV), which ranged from 15.88% (fruit diameter) to 59.77% (number of seeds per fruit). Traits associated with cladode morphology, such as cladode length (CV = 16.26%) and cladode diameter (CV = 18.65%), showed relatively low variability. Similarly, fruit diameter (15.88%), floral scar diameter (16.17%), and flower length (17.91%) also exhibited limited dispersion.

In contrast, moderate variability (CV between 20% and 40%) was observed in fruit-related traits, including fruit length (21.63%), exocarp thickness (25.98%), fruit weight (37.98%), and 100-seed weight (35.10%).

Several traits displayed high variability (CV > 40%), including the length of the longest spine (47.18%), the length of the shortest spine (47.75%), floral scar depth (59.31%), seed length (48.64%), seed diameter (48.36%), and number of seeds per fruit (59.77%). The high dispersion observed in these traits indicates substantial phenotypic variability among the evaluated accessions. In particular, the number of seeds per fruit showed the highest degree of variation, ranging from 43.67 to 678 seeds.

The analysis of minimum–maximum ranges supports these findings, revealing wide amplitudes in variables such as fruit weight (26.33–166 g) and seed-related traits. Likewise, spine lengths, despite their low absolute values, exhibited high relative variability (CV = 47.18% and 47.75% for the longest and shortest spine, respectively).

Overall, vegetative traits — particularly cladode dimensions — showed greater stability compared to reproductive and defensive traits, which displayed markedly higher variability.

3.2. Pearson Correlation Analysis and Variable Selection for PCA

Pearson correlation analysis (Figure 2) revealed numerous significant associations (p < 0.05) among the evaluated variables. Strong positive correlations were observed among fruit size variables, notably between fruit weight and fruit diameter (r = 0.81; p < 0.0001), as well as between fruit weight and fruit length (r = 0.66; p < 0.0001). Additionally, a near-perfect correlation was found between seed length and seed diameter (r = 0.98; p < 0.0001), indicating redundancy between these two variables.

Conversely, the length of the longest cladode spine exhibited significant negative correlations with several reproductive variables, including cladode length (r = −0.64; p < 0.0001), fruit weight (r = −0.56; p < 0.0001), fruit length (r = −0.53; p < 0.0001), and flower length (r = −0.55; p < 0.0001).

The 100-seed weight also showed moderate negative correlations with the number of seeds per fruit (r = −0.49; p < 0.0001) and with flower length (r = −0.39; p = 0.0031). Collectively, these results enabled the identification of association and redundancy patterns among variables, which were subsequently considered in the variable selection process for the multivariate analysis.

Since 105 pairwise Pearson correlations were simultaneously computed (15 quantitative variables), p-values were corrected for multiple comparisons using the Bonferroni method and the Benjamini–Hochberg false discovery rate procedure (BH–FDR), both at α = 0.05. Of the 34 nominally significant correlations (uncorrected p < 0.05), 15 remained significant after Bonferroni correction and 27 after BH–FDR correction. All associations highlighted in the text — including the near-perfect correlation between seed length and seed diameter (r = 0.980; p_Bonf < 0.001), the strong positive correlations among fruit size variables (weight vs. diameter: r = 0.807; p_Bonf < 0.001; weight vs. length: r = 0.658; p_Bonf < 0.001), and the negative correlations between the longest spine length and reproductive traits (r = −0.638 to −0.532; all p_Bonf < 0.01) — retained their significance under both correction methods. The sole exception was the association between 100-seed weight and flower length (r = −0.391), which was significant under BH–FDR (p_BH = 0.015) but not under the more conservative Bonferroni correction (p_Bonf = 0.329); this association is therefore interpreted with caution and is not used as a primary criterion for variable exclusion (Table 3).

3.3. Variable Selection for Principal Component Analysis (PCA)

Prior to PCA, the adequacy of the correlation matrix was confirmed through Bartlett’s test of sphericity (χ² = 142.78, df = 45, p < 0.001) and the Kaiser–Meyer–Olkin index (KMO = 0.59). Although marginally below the conventional threshold of 0.60, this value is consistent with those reported in morphological characterization studies of plant genetic resources and does not compromise the descriptive ordination purpose of the analysis. Individual KMO values revealed that floral scar depth (0.42) and seed diameter (0.27) were the descriptors most responsible for reducing the overall index, reflecting their relatively independent variation with respect to the remaining variables. Both descriptors were retained on the basis of their biological relevance, in accordance with standard practice in germplasm characterization.

The first four components jointly explained 67.68% of the total variance, with PC1 accounting for 31.31%, PC2 for 14.30%, PC3 for 12.08%, and PC4 for 9.99% (Table 4). Analysis of the first two components (PC1 + PC2 = 45.61%) revealed the main axes of phenotypic differentiation among accessions (Figure 3).

PC1 (31.31%) was primarily determined by variables related to reproductive output and fruit size. Number of seeds per fruit showed the highest negative loading on this axis (−0.4474), followed by fruit weight (−0.3904) and cladode length (−0.3785), all contributing negatively — that is, accessions with high values in these traits were positioned toward the negative end of PC1. In contrast, longest spine length presented the highest positive loading (0.4776), followed by weight of 100 seeds (0.3457), indicating that accessions with greater spinescence and heavier individual seeds were located toward the positive end of PC1.

PC2 (14.30%) captured an independent axis of variation associated primarily with structural and floral characters. Floral scar depth showed the highest positive loading on this axis (0.6222), followed by cladode diameter (0.5662), indicating that accessions with deeper floral scars and wider cladodes were positioned toward the positive end of PC2. Seed diameter contributed negatively (−0.3223), reflecting an independent source of variation in seed morphology not captured by PC1.

PC3 (12.08%) was dominated by exocarp thickness (0.6022) and fruit weight (−0.4115) and weight of 100 seeds (−0.4289) in opposition, reflecting variation in fruit structural characteristics independent of overall fruit size.

PC4 (9.99%) was primarily associated with shortest spine length (0.5959) and seed diameter (0.5394), capturing variation in secondary defensive and reproductive traits.

The distribution of accessions in the biplot (PC1 × PC2, Figure 3) showed moderate dispersion without complete separation between groups, though clustering tendencies partially coinciding with the Ward–Gower hierarchical classification were evident. Accessions with the most extreme positive PC1 scores (i.e., high spinescence, heavier seeds, lower fruit production) corresponded predominantly to Group 1 and Group 2 accessions from Carchi, while accessions with the most negative PC1 scores (large fruit, high seed counts, lower spinescence) clustered toward Group 3. Overall, the PCA identified a primary functional gradient structuring the morphological diversity of O. ficus-indica in the inter-Andean dry valleys of northern Ecuador, consistent with the morphological gradient described below (Section 3.6).

3.4. Multiple Correspondence Analysis (MCA)

Figure 4 shows the biplot constructed from the qualitative traits of prickly pear (Opuntia ficus-indica) accessions: cladode color, fruit shape, peel color, pulp color, and petal color.

Multiple Correspondence Analysis (MCA) explained a relatively low proportion of total inertia, which is common in this type of analysis when variables with a large number of categories are included. The first dimension (Dim 1) yielded an eigenvalue of 0.89 and accounted for 4.19% of the inertia (χ² = 527.89; p < 0.001), while the second dimension (Dim 2) yielded an eigenvalue of 0.87 and contributed 3.99% (χ² = 502.92; p < 0.001). Together, the first two axes explained 8.19% of the total variability of the qualitative data.

Despite the moderate percentage of explained inertia, the axes revealed a clear and biologically relevant structuring of the accessions. Dim 1 (horizontal axis) captured the greatest source of variability and was primarily associated with extreme or infrequent combinations of chromatic and morphological traits. At the positive end of this axis, the accessions with the highest coordinates were: 15 (4.15), 16 (4.15), 20 (3.54), 28 (3.54), 13 (2.93), 12 (2.59), 18 (2.59), and 23 (2.59). These accessions were characterized by singular qualitative profiles, particularly regarding skin, pulp, and petal colors. Dim 2 (vertical axis) represented a second independent source of variation, with accessions 24 (2.58), 25 (2.54), 3 (2.09), 15 (2.02), and 28 (1.93) standing out at its positive end, and accessions 12 (−2.14) and 24 (−2.14), among others, at its negative end.

The accessions were broadly distributed across the biplot space, with a marked presence in the extreme positive quadrants of both axes, evidencing considerable qualitative diversity within the collection. The accessions furthest from the origin — primarily 15, 16, 20, 24, 25, 28, 13, 12, 18, and 23 — contributed most to total inertia and represented the most differentiated qualitative profiles. The majority of accessions, by contrast, were concentrated near the origin (approximate coordinates between −1.5 and 1.5), indicating the presence of common or intermediate qualitative profiles widely shared within the collection.

Although the biplot did not include group-based coding, the spatial distribution suggests multivariate differences among groups G1, G2, and G3. The more homogeneous groups tended to cluster near the center, while the most heterogeneous group exhibited greater dispersion toward the extremes.

Overall, MCA confirmed the existence of significant qualitative variability among Opuntia ficus-indica accessions. Although the first two axes explained only 8.19% of total inertia, they were sufficient to identify the most distinctive trait combinations and the most singular accessions (particularly 15, 16, 20, 24, 25, 28, and 13). These results complement the univariate frequency analyses by revealing the simultaneous associations among chromatic and morphological traits.

The low cumulative inertia explained by the first two MCA dimensions (Dim 1 = 4.19%, Dim 2 = 3.99%; total = 8.19%) is consistent with the expected behavior of MCA when applied to datasets with a large number of categories, since total inertia scales with the number of active categories rather than with the number of variables [49]. In the present study, the large number of RHS color categories across five qualitative traits (cladode color, fruit shape, peel color, pulp color, and petal color) generates an extensive active category space, which mathematically limits the proportion of inertia that any individual dimension can capture. Despite this, the axes revealed a biologically coherent structure: Dimension 1 separated accessions with rare or extreme chromatic profiles from those with common intermediate profiles, and the most extreme accessions identified by MCA (particularly 15, 16, 20, 24, 25, and 28) coincided with those showing the greatest morphological distinction in the Ward–Gower clustering, thereby confirming the convergent validity of both multivariate approaches.

The non-random association between qualitative traits and morphotype groups was formally assessed through independence tests (Table 5). Fruit shape showed a significant chi-square association with group membership (χ²(8) = 24.40, p = 0.002, Cramér’s V = 0.471), confirming that the transition from predominantly rectangular shapes in G1 to spherical shapes in G3 is not attributable to chance. Peel color and pulp color exhibited the strongest group discrimination: both variables yielded highly significant Kruskal–Wallis statistics on ordinal RHS codes (H(2) = 38.55, p < 0.001 and H(2) = 38.91, p < 0.001, respectively) and large effect sizes (V = 0.822 and V = 0.805), corroborated by significant chi-square tests on hue-collapsed categories (skin: χ²(4) = 46.47, p < 0.001; pulp: χ²(4) = 44.15, p < 0.001). In contrast, cladode color and petal color showed no significant differences among groups in either the ordinal analysis (H(2) = 2.88, p = 0.237 and H(2) = 0.32, p = 0.853, respectively) or the collapsed chi-square analysis, despite their high Cramér’s V values (0.547 and 0.764). This apparent paradox reflects high intragroup chromatic heterogeneity rather than intergroup differentiation, and is consistent with the MCA biplot, where these traits did not contribute substantially to axis separation. Collectively, three of the five qualitative traits — fruit shape, peel color, and pulp color — showed statistically significant and practically large associations with morphotype group membership, reinforcing the discriminant validity of the Ward–Gower classification.

3.5. Hierarchical Cluster Analysis and Morphotypic Grouping

The Gower distance matrix — which simultaneously integrates all 20 morphological descriptors into a single dissimilarity measure, as opposed to the univariate coefficients of variation reported in Section 3.1 — yielded a mean dissimilarity of 0.412 (SD = 0.071; range: 0.186–0.641; CV = 17.34%). This mean value indicates that, on average, each pair of accessions differs by approximately 41.2% of the scaled descriptor space. The moderate CV (17.34%) reflects a relatively homogeneous distribution of distances, consistent with a collection encompassing both morphologically similar accessions (low distances, predominantly intragroup comparisons) and pairs with high dissimilarity (large distances, mainly between G1 and G3), which justifies the choice of Gower distance as a metric capable of capturing the full range of phenotypic variation present in the data. The pairwise distance distribution was unimodal and approximately symmetric (skewness = 0.070), with no evidence of bimodal structure or outlying accessions dominating the distance space.

The clustering structure of the 55 Opuntia ficus-indica accessions was determined through hierarchical classification analysis applying Ward’s method with Gower distances (Figure 5). The cophenetic correlation coefficient of the dendrogram was r = 0.63 (p < 0.001), indicating an acceptable fit between the original Gower distance matrix and the dendrogram topology, a value consistent with the application of Ward’s method to non-Euclidean distances (Rohlf, 1970). The dendrogram revealed the formation of three morphologically differentiated groups at a cut-off distance of 0.58 Gower units (48.3% of the maximum fusion distance), with internal fusion levels reflecting the degree of phenotypic cohesion among accessions. Group 1 (n = 8), composed exclusively of accessions from Carchi and Ibarra (Imbabura), exhibited the greatest internal homogeneity. Group 2 (n = 26) was the most numerous and displayed a branched internal structure with two distinguishable subclusters, incorporating accessions from all three sampled provinces. Group 3 (n = 21), with notable representation from Pichincha and Imbabura, showed intermediate fusion distances.

3.5.1. Validation of the Optimal Number of Groups

The selection of k = 3 was based on dendrogram gap analysis. The fusion gap at k = 3 (Δ = 0.106) was 6.6-fold larger than at k = 2 (Δ = 0.016), indicating that the partition into three groups introduces substantially more structural information than the two-group solution. From k = 4 onward, fusion gaps declined sharply and monotonically (Δ = 0.033 at k = 4; Table 6), confirming that no partition with k > 3 introduces additional structure. Although the Silhouette and Calinski–Harabász indices were maximized at k = 2 (0.272 and 15.27, respectively), both showed a sharp decline from k = 4 onward with no recovery, confirming that no partition with k > 3 introduces additional structure. For k = 3, the overall Silhouette was 0.169 with positive values across all three groups (G1 = 0.200; G2 = 0.141; G3 = 0.193) and only 2 out of 55 accessions (3.6%) yielding negative coefficients, indicating that k = 3 represents the partition with the greatest structural support and biological coherence, as it distinguishes a geographically restricted morphotype (G1) that k = 2 subsumes without ecological justification.

3.5.2. Group Establishment

Hierarchical cluster analysis (Ward’s method, Gower distance) enabled the classification of 55 Opuntia ficus-indica accessions into three morphologically differentiated groups, whose composition and geographic distribution are detailed in Table 7. Group 1 (G1) comprised 8 accessions, originating exclusively from the cantons of Bolívar and Mira (Carchi) and Ibarra (Imbabura), reflecting a restricted and relatively homogeneous geographic distribution. Group 2 (G2), the most numerous with 26 accessions, exhibited the broadest geographic range, incorporating materials from all three sampled provinces — Carchi (Bolívar, Mira), Imbabura (Ibarra, Atuntaqui, Urcuquí, Pimampiro), and Pichincha (Quito, Cayambe) — consistent with its greater phenotypic heterogeneity. Group 3 (G3), with 21 accessions, also spanned all three provinces, with notable representation from Pichincha (Quito, Cayambe) and the cantons of Ibarra, Mira, Urcuquí, and Pimampiro. The three morphotypic groups were distributed across the full altitudinal gradient sampled in this study (1,381–2,758 m a.s.l.).

The morphological differentiation among the three groups was statistically confirmed through Kruskal–Wallis tests applied independently to each quantitative variable (Table 8). Significant differences among groups were detected for eight of the ten evaluated variables. The strongest discrimination was observed for longest spine length (H(2) = 28.54, p < 0.001), number of seeds per fruit (H(2) = 24.10, p < 0.001), fruit weight (H(2) = 23.06, p < 0.001), and cladode length (H(2) = 17.06, p < 0.001).

Dunn’s post-hoc comparisons with Bonferroni correction (Table 9) revealed that G1 and G2 did not differ significantly from each other in any of the 10 evaluated variables (p > 0.05 in all cases). G3 differed significantly from both G1 and G2 in fruit weight (p < 0.001), fruit length (p < 0.05), fruit diameter (p < 0.01), cladode length (p < 0.05), and number of seeds per fruit (p < 0.001). For longest spine length, G3 differed significantly from both G1 (p = 0.013) and G2 (p < 0.001), while G1 and G2 were statistically indistinguishable (p = 1.000). For 100-seed weight, G2 showed the highest values and differed significantly from G3 (p = 0.001), while G1 occupied an intermediate position (ab). Cladode diameter and exocarp thickness showed no significant differences among groups (p = 0.791 and p = 0.623, respectively).

3.6. Morphotypic Characterization of Accession Groups

The hierarchical cluster analysis (Ward’s method, Gower distance) applied to the combined quantitative and qualitative morphological descriptors yielded three morphologically distinct groups (G1, n = 8; G2, n = 26; G3, n = 21), which differed consistently in fruit size, vegetative structure, spinescence, and reproductive output (Table 7).

Group 1 (G1) — Small-fruited, spinescent morphotype (Figure 6). This group comprised eight accessions collected predominantly from the dry inter-Andean valleys of Carchi province (Bolívar and Mira cantons) and, to a lesser extent, Ibarra (Imbabura). G1 accessions were characterized by the smallest fruit size among the three groups, with a mean fruit weight of 59.97 ± 28.88 g (range: 28.8–109.0 g), mean fruit length of 6.25 ± 1.8 cm, and fruit diameter of 4.38 ± 0.83 cm. Rectangular fruit shape was the most frequent morphology (50% of accessions). Cladodes were moderate in size (33.9 ± 4.1 cm length; 19.3 ± 3.1 cm diameter) and exhibited comparatively high spinescence, with the longest spine averaging 2.46 ± 1.0 cm. Seed number per fruit was the lowest of the three groups (138 ± 113 seeds/fruit), while the weight of 100 seeds was relatively high (1.44 ± 0.51 g), suggesting fewer but heavier seeds. Pulp coloration was predominantly light, concentrated in categories 1, 3, and 6 of the RHS scale. Petal color was moderately variable across categories 7, 11, 14, 18, 19, and 22.

Group 2 (G2) — Intermediate, morphologically diverse morphotype (Figure 7). G2 was the largest group and the most phenotypically variable, including 26 accessions distributed across all three sampled provinces (Carchi, Imbabura, and Pichincha). Fruit weight was intermediate (74.48 ± 20.99 g; range: 26.3–135.0 g), with mean fruit length of 6.17 ± 1.07 cm and diameter of 4.87 ± 0.61 cm. Fruit shape distribution was the most diverse, with approximately equal representation of oval (n = 8), rectangular (n = 8), and spherical (n = 6) forms. Cladodes showed the shortest mean length (32.98 ± 5.30 cm) but the highest spinescence (2.88 ± 0.70 cm longest spine), consistent with the PCA gradient contrasting defensive and productive traits. Seed number per fruit averaged 155.94 ± 78.00. The weight of 100 seeds was the highest among groups (1.50 ± 0.44 g). Cladode color was the most broadly distributed across RHS categories (codes 1–18), and pulp color spanned a wide chromatic range, reflecting the high phenotypic diversity of this morphotype. G2 included accession 41 (Urcuquí), which recorded the maximum fruit weight of the entire collection (135.0 g).

Group 3 (G3) — Large-fruited, low-spinescence morphotype (Figure 8). Group 3 comprised 21 accessions spanning all three provinces, with notable representation from Pichincha (Quito, Cayambe). This group was defined by the largest fruit dimensions recorded in the study: mean fruit weight of 113.23 ± 29.97 g (range: 68.0–166.0 g), mean fruit length of 7.62 ± 1.5 cm, and fruit diameter of 5.61 ± 0.70 cm. Spherical fruit shape was predominant (9 accessions), followed by lanceolate (6 accessions). Cladodes were the longest of the three groups (39.48 ± 4.84 cm), while spinescence was markedly reduced (mean longest spine: 1.28 ± 0.67 cm), consistent with the negative relationship between spine length and fruit size identified in the Pearson correlation analysis and the PCA. Seed number per fruit was substantially higher in G3 (312.00 ± 122.12 seeds/fruit; range: 98–678), while seed weight was the lowest recorded across the collection (1.15 ± 0.32 g). Pulp coloration was notably diverse, distributed across a broad range of RHS categories (codes 10–23), encompassing tonalities from pale yellow-green to deep red-purple. Accession 14 (Pimampiro, Imbabura) recorded the maximum seed count of the entire collection (678 seeds/fruit), while accession 8 (Bolívar, Carchi) registered the highest fruit weight (166 g).

Taken together, the three morphotypes described a clear morphological gradient: G1 showed the highest spinescence (longest spine: 2.46 ± 1.0 cm) and lowest fruit weight (59.97 ± 28.88 g), G3 the lowest spinescence (1.28 ± 0.67 cm) and highest fruit weight (113.23 ± 29.97 g), and G2 intermediate values across both traits (2.88 ± 0.70 cm; 74.48 ± 20.99 g).

3.7. Qualitative Character Analysis

Qualitative traits showed consistent distributional differences among the three morphotypes (Figure 9; Table 10). Fruit shape was the most discriminating trait: G1 was dominated by rectangular shapes (63%), G2 exhibited the greatest diversity with nearly equal representation of oval, rectangular, spherical, and elliptical categories, and G3 was characterized by the predominance of spherical shape (43%), constituting a clear qualitative gradient aligned with Dimension 1 of the MCA.

Chromatic diversity of fruit traits was greatest in G3, which presented the widest range of peel color (S = 18 RHS categories; H’ = 2.85) and the most heterogeneous pulp coloration (categories 10–23, ranging from pale yellow-green to deep red-purple; H’ = 2.56) — attributes of direct commercial relevance for product diversification. In contrast, G2 showed the greatest diversity in vegetative and ornamental traits, with the broadest ranges of cladode color (H’ = 2.26) and petal color (H’ = 2.81), consistent with its wider geographic distribution across all three sampled provinces. G1 exhibited the lowest chromatic diversity across all five evaluated traits. Pielou’s evenness values were consistently high across all traits and groups (J’ = 0.82–0.98; Table 10), indicating that the observed categories were distributed nearly uniformly within each group, with no dominance by one or two frequent categories.

Cladode color and petal color showed no significant differences among groups in either ordinal or chi-square analyses (p > 0.05; Table 5), reflecting high intragroup chromatic heterogeneity rather than intergroup differentiation, consistent with their limited contribution to MCA axis separation. Complete frequency distributions by RHS category and group are presented in Figure 9.

Overall, the diversity analysis reveals a character-specific pattern: G2 represents the most diverse morphotype with respect to vegetative and ornamental characters (cladode color, fruit shape, petal color), whereas G3 exhibits the greatest chromatic diversity in fruit-related characters (peel and pulp color).

4. Discussion

4.1. Breadth and Structure of Phenotypic Variability

The wide range of phenotypic variation observed across the 55 accessions of Opuntia ficus-indica from the inter-Andean valleys of northern Ecuador (CV: 15.88–59.77%) is consistent with the high morphological plasticity that characterizes the species, as documented in germplasm collections from other biogeographic regions. Studies on Tunisian collections reported coefficients of variation exceeding 40% for reproductive and seed characters assessed using the same FAO/CACTUSNET descriptor set [50], while research on Algerian and Moroccan accessions confirmed that spine- and seed-related traits exhibit the greatest relative variability within the species [16,52]. The greater stability recorded for cladode dimensions (CV: 16–19%) relative to reproductive and defensive characters (CV: 35–60%) mirrors the pattern described by Reyes-Agüero et al. (2005) for wild and cultivated collections from Mexico: primary vegetative organs, subject to stricter structural and physiological constraints, display lower phenotypic variance than traits directly linked to reproductive success or plant–herbivore interactions [7]. Recent studies on Ethiopian populations likewise confirmed that qualitative characters such as pulp color and cladode shape exhibit altitudinally structured distributions, with orange coloration predominating at lower elevations and red at intermediate altitudes, suggesting a differential environmental effect on these chromatic traits [53].

The high variability in seed number per fruit (CV = 59.77%; range: 43.67–678 seeds) exceeds that reported by Barbera et al. (1994) for the Italian cultivars ’Gialla’ and ’Rossa’ (CV ≈ 35–40%) [54], suggesting that Ecuadorian Andean germplasm retains a level of reproductive diversity greater than that of commercial varieties subjected to prolonged selection. This outcome is expected given that the majority of the evaluated accessions correspond to semi-wild or naturalized forms under limited artificial selection pressure, in contrast to commercial cultivars subjected to systematic breeding aimed at reducing seed content [55]. Similarly, the high variability in spine length (CV ≈ 47–48%) is consistent with the absence of consistent agronomic selection toward spinelessness, a trait found exclusively in the most advanced horticultural cultivars [7,56].

4.2. Defensive–Reproductive Trade-Off and the Primary Axis of Differentiation

The first principal component (PC1 = 31.31%) defined the central axis of morphological differentiation as a functional opposition between investment in defensive structures—longest spine length (loading: +0.4776)—and reproductive output—seed number per fruit (loading: −0.4474) and fruit weight (loading: −0.3904). This inverse relationship was confirmed by the significant negative correlation between longest spine length and fruit weight (r = −0.562; p < 0.001), and is phenotypically expressed in the separation between morphotype G1 (high spinescence, small fruit) and morphotype G3 (low spinescence, large fruit). This pattern is consistent with the resource allocation trade-off hypothesis between physical defense and reproductive investment proposed for perennial plants in water-limited environments [57,58]. López-Palacios et al. (2019) documented for Mexican populations of O. ficus-indica that spinier forms consistently bear smaller fruits, in agreement with the results reported here [59]. This pattern of separation between spiny and spineless forms was also documented by Peña-Valdivia et al. (2008) in a multivariate analysis of 46 Mexican accessions of Opuntia spp., in which cladode length and width and spine presence constituted the descriptors with the greatest discriminant power in the first principal component [60], in line with the PC1 axis structure observed in the present study.

The domestication process in the genus Opuntia has historically followed a trajectory favoring spine reduction and increased fruit size—traits that are opposed along the PC1 axis of the present study [7,59]. López-Palacios et al. (2019) characterized 114 samples from 17 Opuntia species from the Mexican highlands and concluded that domestication operates as a continuous gradient of morphological change in fruits and seeds, with spinescence being one of the primary characters under artificial selection pressure [59]. From this perspective, morphotype G3—dominated by accessions from Pichincha and Imbabura bearing larger fruits (113.23 ± 29.97 g) and shorter spines (1.28 ± 0.67 cm)—may represent a more advanced domestication stage than G1, whose materials from Carchi retain greater spinescence (2.46 ± 1.0 cm) and smaller fruits (59.97 ± 28.88 g). This interpretation is consistent with the agroecological history of the Ecuadorian inter-Andean valleys, where Imbabura and Pichincha have a more consolidated tradition of prickly pear cultivation for local markets than the Bolívar region (Carchi), where semi-wild forms on degraded hillsides predominate [7].

4.3. Identified Morphotypes and Their Agronomic Interpretation

Morphotype G3 consistently exhibited the highest values for fruit weight (113.23 ± 29.97 g), cladode length (39.48 ± 4.84 cm), and seed number per fruit (312.00 ± 122.12 seeds), confirming it as the group with the greatest agronomic potential for fresh fruit production. The mean fruit weight of G3 falls within the ’large fruit’ range (≥100 g) established by Mediterranean prickly pear commercial standards [55], and exceeds the ranges reported for unselected accessions from Morocco (63–92 g) [16] and Algeria (42–86 g) [52]. However, the lower 100-seed weight in G3 (1.08 ± 0.38 g) relative to G2 (1.50 ± 0.44 g) is consistent with the size–number trade-off described by Barbera et al. (1994) for Italian cultivars, in which larger fruits harbor more seeds of lower individual mass [54]. This pattern reflects a differential reproductive allocation strategy: G3 invests in a greater number of propagules of lower mass, whereas G2 produces individually heavier seeds. From a plant breeding perspective, G2 accessions with high seed weight (accessions 41, 28, 40) constitute relevant candidates for crosses aimed at improving seed quality within large-fruit genetic backgrounds [55].

Morphotype G1, geographically restricted to the cantons of Bolívar and Mira (Carchi) and Ibarra (Imbabura), represented the most homogeneous group with the greatest relative spinescence. Its limited geographic distribution, combined with its internal phenotypic coherence, suggests that these accessions may correspond to semi-wild forms (sensu Section 1) subject to lower human management pressure and potentially adapted to the more arid conditions along the altitudinal gradient. This interpretation is consistent with the findings of Chougui et al. (2016) in northern Algeria, who documented that spiny accessions of O. ficus-indica with smaller fruit sizes were concentrated at sites of greater aridity, and that spine density and length showed significant correlation with precipitation regimes [52]. Adli et al. (2019) described four morphological groups in naturalized accessions from the Algerian steppe—one spineless and three spiny—whose spinescence was correlated with minimum temperatures and annual precipitation at the collection site [61]. The concentration of G1 in the drier valleys of Carchi (altitude: 1,381–2,100 m a.s.l.; precipitation: 350–500 mm/year) is consistent with this pattern and supports the hypothesis that the environmental gradients of the region contribute to structuring the phenotypic variation of the collection. From a genetic resources management perspective, the retention of these morphotypes in ex situ collections is strategically valuable, given that accessions with greater spinescence have demonstrated higher tolerance to limiting water regimes under field evaluation [62].

Morphotype G2, the most numerous and phenotypically heterogeneous, incorporated accessions from all three sampled provinces and exhibited the highest values for 100-seed weight (1.50 ± 0.44 g), as well as the greatest diversity in cladode color, fruit shape, and petal color. This internal heterogeneity is consistent with its broad geographic distribution and suggests that G2 groups materials of diverse origins subjected to different selection and management regimes. In germplasm characterization studies of prickly pear in ex situ collections from Morocco and Tunisia, the groups with the broadest phenotypic range also corresponded to the most numerous and incorporated accessions from multiple localities, reflecting the interconnected nature of vegetative material exchange networks in traditional agroforestry systems [16].

4.4. Discriminant Power of Morphological Characters

Among the qualitative descriptors, peel color and pulp color showed the greatest discriminant capacity between morphotypes (Cramér’s V: 0.822 and 0.805, respectively), whereas cladode color and petal color showed no significant inter-group differentiation (p > 0.05). Díaz-Delgado et al. (2024) obtained similar results in Opuntia accessions from Tenerife (Canary Islands), where canonical discriminant analysis correctly classified accessions primarily based on pulp color, while cladode color showed a low contribution to group separation [63]. Hadjkouider et al. (2017) identified, across five Opuntia species from the Algerian steppe, that of 49 UPOV descriptors, only eight—including flower length and longest spine length—provided genuine discriminant power in the PCA [18]. Similarly, Gallegos-Vázquez et al. (2011) differentiated three groups among 29 Mexican commercial varieties using 24 UPOV quantitative characters, with fruit and seed descriptors showing the highest classificatory capacity [64], which converges with the high Cramér’s V values for peel and pulp color observed in the present study. The high intra-group heterogeneity of cladode color and petal color—evidenced in the MCA biplot and in the high per-group Shannon H’ values (H’ = 1.67–2.81)—may reflect that these characters are under greater environmental control, whereas fruit colors, linked to betalain synthesis, would be under stronger genetic control and thus more diagnostically informative of the genotype [65,66].

Fruit shape was the qualitative character most strongly associated with morphotype membership (χ²(8) = 24.40; p = 0.002; V = 0.471), with a transition from rectangular forms in G1 to predominantly spherical forms in G3. This gradient is consistent with the morphological inventory of the FAO/CACTUSNET collection, in which spherical and elliptical forms predominate among selected horticultural cultivars, while rectangular and lanceolate forms are associated with more primitive or naturalized types [13]. The predominance of spherical forms in G3 is also consistent with the geometric relationship between shape and size: at equal volume, the spherical form minimizes the surface-to-volume ratio, a property relevant to reducing postharvest water loss [55].

4.5. Chromatic Diversity and Agro-Industrial Potential

The high chromatic diversity of pulp color detected in G3 (S = 14 RHS categories; H’ = 2.56), spanning from pale yellow-green to intense red-purple, has direct implications for the O. ficus-indica value chain. Betalains—pigments responsible for the red, purple, and orange hues of the pulp—are bioactive compounds of growing interest to the food and cosmetic industries for their antioxidant and antimicrobial properties and their potential as natural colorants [65,66]. Recent studies have demonstrated that betacyanin and betaxanthin concentrations in the pulp and peel of O. ficus-indica are closely linked to fruit color, with red and orange cultivars exhibiting the highest contents of bioactive compounds and antioxidant activity [67,68]. Accordingly, G3 accessions with deep red to intense purple pulp coloration (RHS categories 18–23) represent priority genetic materials for breeding programs aimed at developing cultivars with high functional and nutraceutical value [66,68]. This potential has been directly quantified: Castellanos-Santiago & Yahia (2008) documented across 10 Mexican cultivars that betacyanin concentrations in red- and purple-colored pulp exceeded those of white and yellow cultivars by up to fourfold, thereby establishing a direct correlation between chromatic category and agro-industrial value [69].

The chromatic diversity observed in petal color (H’ = 1.73–2.81 depending on the group) and cladode color (H’ = 1.67–2.26) may represent an asset for the ornamental and landscaping sector, which is increasingly interested in cacti for their low water requirements and adaptability to degraded soils. However, these characters showed no significant differentiation between morphotypes, suggesting that their variation is relatively independent of selection for productive characters and could be exploited through intra-group crosses.

4.6. Implications for Conservation and Genetic Improvement

The identification of three morphotypes with differentiated phenotypic profiles across the dry valleys of Imbabura, Carchi, and Pichincha provides a structured basis for designing ex situ conservation strategies. G1, with its geographic distribution restricted to the most arid valleys of Carchi and its high internal phenotypic homogeneity, constitutes a genetic pool of priority conservation interest, given that geographically restricted morphotypes are more vulnerable to loss through local disturbances such as land-use changes or extreme climatic events [70]. G2 represents the largest proportion of the total phenotypic diversity of the collection and should be the most broadly represented group in any core collection.

The extreme accessions identified by PCA and MCA—particularly accessions 15, 16, 20, 24, 25, and 28—merit special attention as genetic resources with the potential to broaden the phenotypic base of breeding programs. Accession 41 (Urcuquí, Imbabura; G2), which recorded the highest fruit weight in the entire collection (135 g), and accession 14 (Pimampiro, Imbabura; G3), with the highest seed number per fruit (678), represent extreme phenotypes relevant for direct selection or as parental lines in crosses. Although reducing seed number is a commercially desirable breeding objective [55], the germplasm bank must preserve variability in this character for future reproductive biology and genetic studies of the species. In sum, the results of the present study contribute to filling the information gap on the morphological diversity of O. ficus-indica in Ecuador and lay the groundwork for molecular and physiological studies that will elucidate the genetic and environmental mechanisms governing phenotypic differentiation along the altitudinal and climatic gradient of the inter-Andean valleys of northern Ecuador [59,71]. The future projection of these materials in the context of climate change adds urgency to their conservation: Acharya et al. (2019) noted that O. ficus-indica accessions with morphology closer to wild forms exhibit greater resilience to water deficits in arid production systems [62], thus conferring on morphotypes such as G1 an additional functional value beyond their taxonomic relevance.

5. Conclusions

This study provides the first systematic multivariate characterization of Opuntia ficus-indica (L.) Mill. germplasm from the dry inter-Andean valleys of northern Ecuador, documenting substantial morphological diversity across 55 accessions through the integrated application of 20 FAO/CACTUSNET and UPOV descriptors, principal component analysis, multiple correspondence analysis, and Ward–Gower hierarchical clustering. The breadth of phenotypic variation recorded (CV: 15.88–59.77%) is comparable to that of long-established germplasm collections from the Mediterranean Basin and sub-Saharan Africa, underscoring the relevance of inter-Andean dry valleys as reservoirs of intraspecific diversity for this globally important xerophytic crop.

Hierarchical clustering resolved three morphologically and statistically distinct groups: a small-fruited, highly spinescent morphotype (G1) restricted to the most arid sites of Carchi, likely representing semi-wild populations that may retain adaptive traits for drought-tolerance breeding; a phenotypically diverse intermediate morphotype (G2) distributed across the full altitudinal gradient; and a large-fruited, low-spinescence morphotype (G3; mean fruit weight: 113.23 ± 29.97 g) consistent with a more advanced domestication trajectory and of greatest agronomic interest for fresh-market production. The primary axis of morphological differentiation (PC1 = 31.31%) revealed a significant negative association between spinescence and reproductive output — a morphological gradient consistent with documented domestication patterns in Opuntia, although its mechanistic basis requires ecophysiological confirmation. The strong discriminant power of fruit peel and pulp coloration (Cramér’s V > 0.80) identifies betacyanin-rich accessions as priority materials for nutraceutical applications.

Collectively, these results establish a morphological baseline for O. ficus-indica in Andean agroecosystems and identify priority accessions for ex situ conservation and genetic improvement. Inherent limitations include the absence of molecular markers — which precludes distinguishing genetic differentiation from phenotypic plasticity — the lack of formal environment–morphology modeling, and a sample size (55 accessions, three provinces) that does not fully cover the distribution range. Future work should prioritize: (i) SSR or SNP genotyping of the identified morphotypes to assess genetic structure; (ii) controlled water-stress experiments to quantify phenotypic plasticity; and (iii) multi-site agronomic trials evaluating yield, post-harvest quality, and drought tolerance of G3 accessions under contrasting altitudinal conditions.