Sexual Dimorphism and Breed Characterization of Creole Hens through Biometric Canonical Discriminant Analysis across Ecuadorian Agroecological Areas

Simple Summary The first step towards the protection and valorization of the genetic resources of a country is their definition. Although Ecuadorian zootechnical species are very diverse, they are scarcely characterized, and hence the efforts towards their protection are not as fruitful as they could be. The present paper approaches the biometric characterization of the Creole hen population in Ecuador through the study of sexual dimorphism and the differentiation of an agroecologically-based structured population using fourteen zoometric measures as differentiation criteria. Highlands region provinces of Cotopaxi and Tungurahua were the most zoometrically distant from the rest. However, Morona Santiago province population in the Amazonian region slightly differs from those in Guayas, Chimborazo and Bolívar in the Coastal and Highlands regions, respectively. The dual-purpose nature of Ecuadorian Creole hen resources enables the implementation of breeding programs that seek meeting a wider scope of public demands, through the definition of the agroecologically-based breed differentiated production of local hen eggs and meat. Abstract Despite Ecuador having a wide biodiversity of zootechnical species, characterization studies of these genetic resources are scarce. The objective of this research was to perform the biometric characterization of the Creole hen population through 14 quantitative zoometric measures. We evaluated 207 hens and 37 roosters from Ecuador’s three agro-ecological regions: the Sierra (highlands) region (Bolivar, Chimborazo, Tungurahua and Cotopaxi provinces); the Costa (coastal) region (Guayas); and the Oriente Amazonian region (Morona Santiago). Sexual dimorphism was assessed using one-way analysis of variance (ANOVA). Body dimensions were generally significantly higher for males (p < 0.05), especially for length of head, beak, neck, dorsum, tarsus, thigh, leg, and middle finger. Then, individuals were biometrically clustered into populations after a stepwise canonical discriminant analysis (CDA) computing interpopulation Mahalanobis distances. Agroecologically-based structured populations were identified when zoometrical criteria were used to classify the animals. Cotopaxi and Tungurahua provinces were reported to be the most distant from the rest, with a slight differentiation of the Morona Santiago province population from those in Guayas, Chimborazo and Bolívar. Conclusively, Ecuadorian Creole hens were higher than longer contrasting light hen breeds, which favors their dual-purpose aptitude. Hence, the development of selection programs aimed at Ecuadorian differentiated entity of production of eggs and meat across agro-ecological areas is feasible.


Introduction
Genetic diversity of domestic hens existing across Ecuador is not only promoted by climatic stratification but also natural and human driven selection. From the red junglefowl (Gallus gallus), the most likely ancestor of these avian populations [1], the effects of natural selection may have resulted in a high heterogeneity and variability of the morphological characteristics of fowl, with a high potential to adapt to the different environmental conditions [2][3][4][5][6].
Apart from naturally driven selection processes and natural migratory movements, genetic variability in local chicken populations may have been conducted as a result of human action. For instance, human-made migration processes [7] brought about the widespread distribution of poultry genetic material, given the size of animals was convenient and facilitated transport, favoring the expansion of these fowl across the different agroecological levels [3].
These factors led to genetic divergences contributing to poultry production under a family-run backyard system usually developed by each household's women [3]. Husbandry practices characterized by the use of rustic animals in free range conditions with a low capital investment, which enables assuming a relatively low economic risk to implement an efficient productive management to produce high-biological-value protein sources such as meat and eggs [8,9]. Additionally, these products are preferred among consumers because of their pigmentation, taste, and lean quality of meat [10,11]; which translates into acceptable income that returns to each family, closing the cycle [12][13][14][15][16][17].
Breeds originating in the Old World were introduced to Latin American territories by the Spanish colonists and adapted to the different agroecological areas and conditions that they found, forming what has traditionally been addressed as Creole hen populations. For decades, these creole populations occupied local productive niches and evolved towards their current state, but still lacked the necessary characterization actions that may help consolidating and protecting them. In parallel, current breed development and formation until the XVII century provided the basic elements for the directed selection of our days and for the pursuit of concrete characteristics of interest to the farmer or producer. In this context, a new conglomerate of breeds and commercial lines formed in the first world were introduced into developing countries in an attempt to fulfil the growing market demands at a lower cost [18].
This global situation resulted in an alarming loss in the biodiversity of animals of zootechnical interest that the region faces nowadays. According to Food and Agriculture Organization (FAO), the endangerment risk that 81 percent of Latin America and the Caribbean avian breeds are exposed to is unknown [19], as even censuses are not appropriately registered. The increased risk of a population whose endangerment status is unknown is based on the fact that measures towards its protection are not implemented. In this regard, efforts are being made to maintain, conserve and, in turn, benefit from their most profitable or useful traits, such as disease or stress resistance, in commercial breeding plans [18].
Not only local hen breeds face a serious risk of extinction, but there is also a simultaneous loss of the traits that made them survive after the evolutionary process that they followed when they arrived in and adapted to the lands to which they were introduced. Creole hens present a good ability to scavenge and forage, have good maternal qualities, and are hardier than exotic breeds with higher survival rates and minimal care and attention requirements. This rusticity is one of those traits to positively influence avian zootechnical production, given its implication with the adaptation ability of animals to the environment in which they are produced.
After a period characterized by a lack of actions regarding local genetic resources conservation, with policies more likely focusing on intensive production, morphological characterization studies in fowl started being run again in Ecuador. These studies lay the basis for local resources conservation and breeding plans. Zoometric traits have widely been reported to depend on an inherited basis and to be suitable means of prediction for the live weight of the individuals [20][21][22]. Thus, they may play an important role in the subsequent performance of animal carcasses [23]; a relationship that translates to new potential selection criteria, seeking the maximization of the profitability of the products derived from such local genetic resources.
Despite the fact that research projects seeking the zoometrical characterization of Ecuadorian local hen breeds started being implemented using univariate analysis, there is still a patent lack of knowledge regarding the differentiation of such local populations, and hence policies towards the protection of such genetic resources cannot be implemented properly. Therefore, the aim of this study was to perform differentiated zoometric characterization of Creole hens through the application of a canonical discriminant analysis (CDA) to provide insights on the possible clustering patterns described by the population and into which subpopulations can be distinguished using Ecuadorian provinces as the criteria of origin [24]. Conclusively, this approach will enable quantification of the large existing phenotypic variability in the Ecuadorian creole hen population as a strategy to facilitate the rational development of such productively important avian local resources and their use, and to direct the implementation of conservation strategies aimed at ensuring their survival in the competitive world of poultry production and future consolidation as breeds.

Sample Size and Distribution
The whole sample comprised 281 fowl, 244 hens (84.84%) and 37 roosters (15.16%), evaluated across the three regions and six provinces, as follows: the Sierra Region (Andean highlands), with the provinces of Bolivar (31), Chimborazo (70), Tungurahua (35), Cotopaxi (32); the Amazonian Oriente Region (eastern rainforests) with Morona Santiago (38); and the Costa Region (Pacific coastal lowlands) with Guayas (28), respectively. Stevens [25] provided a very thorough discussion of the sample sizes that should be used to obtain reliable results for canonical discriminant analysis. Strong canonical correlations in the data (R > 0.7), even in cases of relatively small samples (around n = 50) will be detected most of the time. However, to obtain reliable estimates of the canonical factor loadings for interpretation, hence, to be able to draw valid conclusions, Stevens recommends that there should be at least 20 times as many cases as variables included in the analysis, if one wants to interpret the most significant canonical root only, as is the case in our study. To arrive at reliable estimates for two canonical roots, Barcikowski and Stevens [26] recommend, based on a Monte Carlo study, to include 40 to 60 times as many cases as variables.

Study Site Characterization and Sample Animals Management
The study was conducted under field conditions from January to December 2015 and from January 2017 to August 2018. The animals comprising the sample were raised and kept by backyard producers who did not present evidences of crosses with commercial lines among the effectives of their farms, following a randomized design. The map of provinces and climatic floors of Ecuador are shown in Figures 1 and 2. Agroecological zones are detailed in Table 1.      Animals were reared under extensive backyard conditions, and were not vaccinated against viruses or parasites such as coccidia, nor treated against parasites. Chickens were fed on organic corn and were occasionally supplemented with household wastes, vegetables, and other sources of minerals from each area. Antibiotics and multivitamins were not administered.

Biometric Data Collection
Biometrical analysis was performed on each animal, measuring the sixteen quantitative variables proposed by FAO [28]. A summary of the biometric variables measured and the procedure followed is shown in Table 2. Quantitative data was obtained using a digital scale, a gauge with 0.02 mm accuracy, and a tape measure. All the biometric information was collected in a structured file including georeferencing for each producer together with zoometric measurements. Mean and standard deviation of each measurement were computed for each sex and province. one-way analysis of variance (ANOVA) was carried out using the MEANS Statement from the PROC GLM routine of the S.A.S. 9.4 software [31] to determine the existence of differences in the means for the fourteen variables measured between males and females and across provinces. Then, the WALLER option was used to perform post hoc Waller-Duncan k-ratio t test on all main effect to measure specific differences between pairs of means (p < 0.05). Waller and Duncan [32] and Duncan [33] take an approach to multiple comparisons that differs from all the methods previously discussed in minimizing the Bayes risk under additive loss rather than controlling type I error rates. Furthermore, this range test uses the harmonic mean of the sample size, which makes it preferable if the sample sizes are unequal [34].

Canonical Discriminant Analysis (CDA)
The multivariate technique involved the use of canonical discriminant analysis on the 14 biometric measurements using the province to which each animal belonged as a labeling classification criteria, to identify the variation provided by the different variables measured under study, and to establish clusters that may identify and outline subpopulations [35][36][37][38]. Hence, we determined the percentage of correctly allocated individuals in their populations of origin in comparison to those animals which were statistically misclassified or attributed to a different province from the one in which they were sampled, to discover a linear combination of quantitative morphological variables that provide maximum separation between the potentially existing different populations when the classification criterion was the province in which the animals were located. CDA was also used to plot pairs of canonical variables to help visually interpret group differences. Variable selection was performed using Forward Stepwise (FSTEP) multinomial logistic regression algorithms. Canonical discriminant analysis was performed using the CANDISC Procedure from the PROC CANDISC routine of the S.A.S. 9.4 software [31].

Canonical Correlation Dimension Determination
Canonical correlation is a form of correlation relating two sets of variables. The maximum number of canonical correlations between two sets of variables is the number of variables in the smaller set. The first canonical correlation is always the one which explains most of the relationship [39]. The canonical correlations are interpreted as Pearson's r, and hence their square is the percentage of variance in one set of variables explained by the other set along the dimension represented by the given canonical correlation (usually the first); that is, Rc-squared is the percentage of shared variance along this dimension [40]. As a rule of thumb, some researchers state that a dimension will be of interest if its canonical correlation is 0.30 or higher, corresponding to about 10% of variance explained. Despite this, some researchers report just the first canonical correlation; it is recommended that all meaningful and interpretable canonical correlations are reported [41].

Canonical Discriminant Analysis Efficiency
Wilks' Lambda test assesses which variables significantly contribute to the discriminant function. As a rule of thumb, the closer Wilks' lambda is to 0, the more the variable contributes to the discriminant function. The significance of Wilk's Lambda can be tested using Chi-Square, and then, if the p-value if less than 0.05, we can conclude that the corresponding function explains group adscription well [42]. For small sample sizes or a small number of treatments, the limiting chi-squared or normal distributions may not adequately describe the actual probability distributions of the test statistics. Here, a finite approximation may be more appropriate than using the limiting distribution. One such method is Fisher's F approximation for Wilks' lambda by Rao [43] as developed in Chávez [44]. According to these authors, under normality conditions, this procedure performs more accurately than a χ 2 approximation [45].

Canonical Discriminant Analysis Model Reliability
Box's M is used to test the assumption of equal covariance matrices in Multivariate Analysis of Variance (MANOVA) and discriminant function analysis (DFA). Box's M has very little power [46] for small sample sizes; hence when we work with a small sample a nonsignificant result may not necessarily indicate that the covariance matrices are equal. In contrast, for large samples a statistically significant result can be reported when it does not actually exist. To address this particular issue, a smaller alpha level (p < 0.001) is recommended [47]. Some authors suggest that Box's M is highly sensitive, hence unless p < 0.001 and sample sizes are unequal, we should ignore it. However, if the results are significant and you have unequal sample sizes, the test is not robust [48].
In multiple regression, another assumption that should be tested for is multicollinearity. The variance inflation factor (VIF) is used as an indicator of multicollinearity. Computationally, it is defined as the reciprocal of tolerance: 1/(1 − R 2 ).
Various recommendations for acceptable levels of VIF have been published in the literature. Perhaps most commonly, a value of 10 has been recommended as the maximum level of VIF [49][50][51][52]. The VIF recommendation of 10 corresponds to the tolerance recommendation of 0.10 (1/0.10 = 10). However, a recommended maximum VIF value of 5 [53] and even 4 [54] can be found in the literature.
For example, a VIF of 8 implies that the standard errors are larger by a factor of 8 than would otherwise be the case, if there were no inter-correlations between the predictor of interest and the remaining predictor variables included in the multiple regression analysis.

Canonical Coefficients and Loading Interpretation and Spatial Representation
A preliminary principal component analysis (PCA) was performed to minimize overall variables into few meaningful variables that contributed most to variations in the populations. As a result, half wing radius ulna length (hwrul) and distal phalanx wing length (dpwl) were discarded, given they reported a component loading lower than |0.5| which suggested their redundant confounding nature, which may base on the fact that they comprise the total length of the wing, defined by proximal humerus wing length (phwl). Discriminant function analysis was used to determine percentage assignment of individuals into their own populations.
The traditional approach to interpreting discriminant functions examines the sign and magnitude of the standardized discriminant weight (also referred to as a discriminant coefficient) assigned to each variable in computing the discriminant functions. Small weights may indicate either that a certain variable is irrelevant in determining a relationship or that it has been discarded because of a high degree of multicollinearity.
Discriminant loadings reflect the variance that the independent variables share with the discriminant function. In this regard, they can be interpreted like factor loadings in assessing the relative contribution of each independent variable to the discriminant function.
In either simultaneous or stepwise discriminant analysis, variables that exhibit a loading of >|0.40| or higher are considered substantive, indicating substantive discriminating variables. With stepwise procedures, this determination is supplemented because the technique prevents nonsignificant variables from entering the function. However, multicollinearity and other factors may preclude a variable from entering the equation, which does not necessarily mean that it does not have a substantial effect. Loadings are considered to have relatively higher validity than weights as a means of interpreting the discriminating power of independent variables because of their correlational nature.
Standardized coefficients allow you to compare variables measured on different scales. Coefficients with large absolute values correspond to variables with greater discriminating ability. Also, discriminant scores can be computed by using the standardized discriminant function coefficients applied to data that have been centered and divided by the pooled within-cell standard deviations for the predictor variables, as discussed in IBM Corp. [55].
The data were standardized following standard procedures of Manly [56] before squared Mahalanobis distances and principal component analysis were computed. Squared Mahalanobis distances were computed between populations using the following formula: where D 2 ij is the distance between population i and j, COV −1 is the inverse of the covariance matrix of measured variable x, and Y i and Y j are the means of variable x in the ith and jth populations, respectively. The squared Mahalanobis distance matrix was converted into a Euclidean distances matrix and used to build a dendrogram using unweighted pair-group method using arithmetic mean (UPGMA) via agglomerative hierarchical cluster procedure with the software DendroUPGMA by Garcia-Vallvé and Pere Puigbo [57]. The Mahalanobis squared distance, defined as the square of the distance between the measures of the standardized values of Z (centroids), was used this way to verify whether there were significant differences between provinces [58].

Discriminant Function Cross-Validation
To establish whether the percentage of correctly classified cases is high enough to consider that the discriminant functions issue valid results, as a form of significance we can use leave-one-out cross-validation option. Classification accuracy achieved by discriminant analysis should be at least 25% greater than that obtained by chance.
These results can be supported by Press' Q statistic. This parameter can be used to compare the discriminating power of our function to a model classifying individuals at random (50% of the cases correctly classified), as follows (2) where N is the number of individuals in the sample, n is the number of observations correctly classified (as a coefficient ranging from 0 to 1), and k is the number of groups. The next step is to compute the critical value, which equals the chi-square value at 1 degree of freedom. It is advisable to let alpha equal 0.05. When Q exceeds this critical value, classification can be regarded as significantly better than chance, thereby supporting cross-validation.

Descriptive Statistics, ANOVA and Waller-Duncan Post-Hoc Test
Morphometric analysis indicated highly significantly differences when males were compared to females, as shown in Table 3. 6.54 a 6.10 b a lower mean; b higher mean (p < 0.05). When there is not a significant difference between sexes superindex letters are the same ( a ). Table 4 shows a summary of the significant results of ANOVA and post hoc Waller-Duncan test for the zoometric characteristics of Ecuadorian creole hens across the provinces. The variables analyzed have a very high and variable coefficient of variation across biometric traits.

Canonical Discriminant Analysis Model Reliability
The value of p < 0.05 obtained for Box's M test means the data did not differ significantly from multivariate normal and we could proceed with the analysis. Wilk's lambda statistic was used to assess whether canonical discriminating functions contributed significantly to the separation of treatments, that is, it was used to test the meaning of the discriminating function Table 5. All zoometric variables were included at a preliminary stage of the analysis performed in this study. Tolerance (1/R 2 ) and variance inflation factor (VIF) were analyzed to identify those variables that were responsible for multicollinearity between variables. This analysis revealed that the variables proximal humerus wing length (phwl), half wing radius ulna length (hwrul) and distal phalanx wing length (dpwl), turned out to be highly related (VIF > 4). Therefore, we decided to retain distal phalanx wing length (dpwl) in the analysis, because that measure results from the combination of proximal humerus wing length (phwl) and half wing radius ulna length (hwrul), given their lower VIF. After the removal of redundant variables, the results for tolerance and VIF can be seen in Table 6.

Canonical Coefficients and Loading Interpretation and Spatial Representation
The canonical discriminant analysis identified five discriminating canonical functions. The first had a high discriminatory power, as denoted by the eigenvalue of 5.422. The results are presented in Table 7. The first function obtained explains 91.9% of total variance. The fifth function contributes to the explanation of variance with 49.9% of the information to the analysis, that is, relatively low. The results for the tests of equality of group means to test for differences across provinces once redundant variables have been removed are shown in Table 8. The greater the value of F and the lower the value for Wilks' Lambda, the better the discriminating power a certain variable has and the lower the rank position it presents. Those variables presenting equal values of lambda and F had equivalent discriminatory power, as shown by beak and dorsal length. When this happens, it is necessary check whether these similarities are based on a multicollinearity problem or are because the variables, indeed, have a similar discriminant power. Once F and Wilks' Lambda had been assessed, we evaluated the magnitude of standardized and non-standardized coefficients, reported in Table 9, to determine whether there had been a reduction in the discriminant power of individual variables as a result of multicollinearity between pairs, which implies a reduction in the separate discriminant power of each of the two variables involved in the multicollinear relationship.
As shown in Table 8 for the variables neck and dorsal length, we observe that standardized coefficients fell below 0.4, hence, hence there was a decrease in the discriminating power of the non-individual function as a result of the effect of multicollinearity, that is, neck and dorsal length variables being related and explaining a somehow redundant fraction of variability. The greater the reduction in the standardized coefficient, the more important the multicollinearity problem between variables holding similar Wilks' lambda and F values. Standardized coefficients are shown in Table 9. According to Hair Jr [59], absolute values below |0.3| are indicative of multicollinearity problems when F and Wilks' Lambda have been previously checked to be similar for a certain pair of variables.
Unstandardized coefficients, calculated on raw scores for each variable, are of most use when the investigator seeks to cross-validate or replicate the results of a discriminant analysis or to assign previously unclassified subjects or elements to a group. As we are assessing the potential misclassification of individuals belonging to previously defined populations as a way to define such populations themselves, we must interpret standardized coefficients, and hence unstandardized coefficients were discarded [60]. Furthermore, the unstandardized coefficients cannot be used to compare variables or to determine what variables play the greatest role in group discrimination because the scaling for each of the discriminator variables (i.e., their means and standard deviations) usually differ. The maximum number of canonical discriminant functions generated is equal to the number of groups minus one. In the present study, the number of canonical discriminant functions was 5 for each series, as we used the six provinces as a labelling criterion. After the evaluation of standardized coefficients, the resulting discriminant functions were as follows: To determine which is the variable that we have to discard out of each pair for which a multicollinearity problem has been detected, we checked discriminant loadings, which are presented in Table 10. Discriminant loadings measure the existing linear correlation between each independent variable and the discriminant function, reflecting the variance that the independent variables share with the discriminant function. In this regard, they can be interpreted like factor loadings in assessing the relative contribution of each independent variable to the discriminant function. A graphical representation of discriminant loadings is shown in Figure 3, with those variables whose vector extends further apart from the origin being the most representative discriminating ones.  A territorial map was created by plotting the discriminating values for each observation (Z) for the first function on the x axis and those values for the second discriminant function on the y axis. Figure 4 graphically depicts the canonical discriminant analysis of individuals across the six sampling provinces. Cotopaxi and Tungurahua provinces are evidently independent, while the populations from the provinces of Guayas, Chimborazo and Bolívar displayed significant overlap. Further away from the latter three, we find the population from the province of Morona Santiago. Centroids designed the central observation for each province group. The probability that an unknown case belongs to a particular group was calculated by measuring the relative distance of Mahalanobis to the centroid of a population. To compute discriminant scores or centroids, we substituted the mean for each possible province in the three first dimensions [61]. Then, to calculate the optimal cut-off point, that is, the probability of classification we followed the procedures in Hair, Black, Babin and Anderson [52]. Then we could determine whether a certain case was appropriately classified. A territorial map was created by plotting the discriminating values for each observation (Z) for the first function on the x axis and those values for the second discriminant function on the y axis. Figure 4 graphically depicts the canonical discriminant analysis of individuals across the six sampling provinces. Cotopaxi and Tungurahua provinces are evidently independent, while the populations from the provinces of Guayas, Chimborazo and Bolívar displayed significant overlap. Further away from the latter three, we find the population from the province of Morona Santiago. A territorial map was created by plotting the discriminating values for each observation (Z) for the first function on the x axis and those values for the second discriminant function on the y axis. Figure 4 graphically depicts the canonical discriminant analysis of individuals across the six sampling provinces. Cotopaxi and Tungurahua provinces are evidently independent, while the populations from the provinces of Guayas, Chimborazo and Bolívar displayed significant overlap. Further away from the latter three, we find the population from the province of Morona Santiago. Centroids designed the central observation for each province group. The probability that an unknown case belongs to a particular group was calculated by measuring the relative distance of Mahalanobis to the centroid of a population. To compute discriminant scores or centroids, we substituted the mean for each possible province in the three first dimensions [61]. Then, to calculate the optimal cut-off point, that is, the probability of classification we followed the procedures in Hair, Black, Babin and Anderson [52]. Then we could determine whether a certain case was appropriately classified. Centroids designed the central observation for each province group. The probability that an unknown case belongs to a particular group was calculated by measuring the relative distance of Mahalanobis to the centroid of a population. To compute discriminant scores or centroids, we substituted the mean for each possible province in the three first dimensions [61]. Then, to calculate the optimal cut-off point, that is, the probability of classification we followed the procedures in Hair, Black, Babin and Anderson [52]. Then we could determine whether a certain case was appropriately classified.
It can be observed that the provinces of Cotopaxi and Tungurahua are located in different places on the Cartesian plane, that is, remote from each other. The opposite situation happened with the provinces of Bolívar, Guayas and Chimborazo, and it can also be observed that Morona Santiago is slightly separated from the three latter provinces.

Discriminant Function Cross-Validation
When classification and leave-one-out cross-validation matrices are evaluated, it can be observed that 97.14% has been estimated to be correctly classified for Bolivar, with 88.57% being validated for the same province. For Chimborazo, 93.33% has been correctly classified, with 89.33% validated. For Cotopaxi, Guayas, Morona Santiago and Tungurahua the total of observation was appropriately classified with 90.63%, 50%, 80% and 94.44% being validated, respectively.
Cross-validation reported a result for Press Q parameter of 1064.1418 (Press Q = [N − (nK)] 2 /N(K − 1) = 227 − (221 × 6) 2 /227(6 − 1)). Hence, Q was above 6.63 (significance level of 0.01), Chi 2 s critical value for a degree of freedom at a chosen confidence level. Predictions were significantly better than chance, according which it would have a correct classification rate of 50% The absolute values of Mahalanobis' distances between the local populations of the six provinces involved in the analysis are shown in Table 11. The shortest distance is found between Bolivar and Chimborazo, while the longest distances are those found between the Province of Cotopaxi and the rest. Contrastingly, the distances of Tungurahua with Bolívar, Chimborazo, Cotopaxi are similar. The level of statistical significance for all Mahalanobis distances was high and similar (p < 0.0001). Table 11. Distance of Mahalanobis between locations (above the diagonal) and F statistics (numerator degrees of freedom (dfn) = 6, denominator degrees of freedom (dfd) = 214) for the square distances between locations (below the diagonal). Mahalanobis' distance analysis is based on the analysis of generalized squared Euclidean distances adjusted for unequal variances. The Mahalanobis distance (D 2 ), defined as the square of the distance between the measures of the standardized values of Z, was used to verify whether there were significant differences between the provinces. Thus, the greater the value of the distance, the greater the distance between the means of the provinces considered as well [58]. As can be seen in the dendrogram Figure 5, the provinces of Bolivar-Chimborazo, Guayas-Morona Santiago, and Cotopaxi-Tungurahua are populations represented as subpopulations.

Items
The shortest Euclidean distance was observed between the provinces of Bolivar and Chimborazo; whereas the opposite happened between the Province of Cotopaxi and the others. The distances between Tungurahua and Bolívar, Chimborazo, or Cotopaxi are similar. In contrast, Morona Santiago is slightly far from the provinces of Guayas, Chimborazo, Bolívar.

Discussion
The morphometric measurements show highly significant differences in relation to sex, as reported by Yakubu and Salako [62] in indigenous chickens in Nigeria, which reported such differences to be based on the hormonal effects of sex that condition growth. These results were consistent with similar results in the literature [6,63].
The high coefficient of variation observed in the results is similar to that reported in different populations of chickens in Mexico [64], and also in indigenous chickens in Nigeria [61], which demonstrates the variability of the morphometry in the birds studied, which may be due to genetic divergence processes followed by the populations, such as migration [65], which resulted in the morphological modification of the populations to adapt to the characteristics of the different environments and the orography to which the birds were introduced [63, [66][67][68].
Measurements for head length are higher than those measured in Batsi Alak Hens of Mexico whose mean for males and females varied from 4.16 to 4.6 cm [67], and lower than those reported for Yoruba ecotypes of Nigeria with an average of 9.90 cm [68].
In terms of crest length, the values are lower than those measured in indigenous Nigerian roosters and higher than those indicated for hens of the same country [62]. Likewise, Yoruba and Fulani ecotypes [67] reported a similar value to that reported for the province of Cotopaxi when comparing fowl in general, without considering males and females separately [69]. However, in autochthonous Catalonian chicken breeds, such as Patridged Penedesenca and Blonde Empordanesa, we observe crest sizes double those measured in the present study [70].
For beak length, the values are analogous to those found in Botswana hens. This could be supported by the fact that both studies were conducted across three agroecological regions [71]. However, higher values were reported by Yakubu and Salako [62] in indigenous Nigerian fowl for males and females, as well as in native Catalonian breeds Partridged Penedesenca and Blonde Empordanesa. The opposite situation was described by Batsi Alak hens from Mexico [67] and Fulani

Discussion
The morphometric measurements show highly significant differences in relation to sex, as reported by Yakubu and Salako [62] in indigenous chickens in Nigeria, which reported such differences to be based on the hormonal effects of sex that condition growth. These results were consistent with similar results in the literature [6,63].
The high coefficient of variation observed in the results is similar to that reported in different populations of chickens in Mexico [64], and also in indigenous chickens in Nigeria [61], which demonstrates the variability of the morphometry in the birds studied, which may be due to genetic divergence processes followed by the populations, such as migration [65], which resulted in the morphological modification of the populations to adapt to the characteristics of the different environments and the orography to which the birds were introduced [63, [66][67][68].
Measurements for head length are higher than those measured in Batsi Alak Hens of Mexico whose mean for males and females varied from 4.16 to 4.6 cm [67], and lower than those reported for Yoruba ecotypes of Nigeria with an average of 9.90 cm [68].
In terms of crest length, the values are lower than those measured in indigenous Nigerian roosters and higher than those indicated for hens of the same country [62]. Likewise, Yoruba and Fulani ecotypes [67] reported a similar value to that reported for the province of Cotopaxi when comparing fowl in general, without considering males and females separately [69]. However, in autochthonous Catalonian chicken breeds, such as Patridged Penedesenca and Blonde Empordanesa, we observe crest sizes double those measured in the present study [70].
For beak length, the values are analogous to those found in Botswana hens. This could be supported by the fact that both studies were conducted across three agroecological regions [71]. However, higher values were reported by Yakubu and Salako [62] in indigenous Nigerian fowl for males and females, as well as in native Catalonian breeds Partridged Penedesenca and Blonde Empordanesa. The opposite situation was described by Batsi Alak hens from Mexico [67] and Fulani ecotypes of Nigeria [68], which reported similar mean values to those of Bolivar, Guayas, Morona Santiago. This could be ascribed to the similarity between the climates of the locations in which the study took place.
For the neck length trait, the values found were equivalent to those found in Partridged Penedesenca and Blonde Empordanesa [70] and Yoruba and Fulani Nigerian ecotypes [68]; while lower values were found in indigenous Nigerian hens for both males and females [62]. The highest values in literature were reported for Batsi Alak Hens from Mexico with measures of 19 to 17 cm in males and females, respectively [67].
Body and dorsal lengths along with head length have been related in literature to the potential of animals for egg production [68]. When these data are compared with those reported by Moazami-Goudarzi [5] and their studies of local Tanzanian chicken ecotypes, the values for males of the Singamagazi ecotype were slightly higher than the average reported for Morona Santiago hens, but comparable to the males of the Kuchi ecotype. However, for the male of Mbeya ecotype and female of the Singamagazi ecotypes [5], these values were similar to those reported for the Patridged Penedesenca hens and for the Blonde Empordanesa [70] and equivalent to those found in Bolívar, Chimborazo, and Guayas. Studies conducted on Botswanan hens across three different agroecological areas found similar values for females and males [71] to those measured in Cotopaxi and Tungurahua. In turn, Yakubu and Salako [62] reported higher average measurements than those of Ecuadorian creole hens found in this study or those found for the Fulani and Yoruba.
Higher values for the ventral length variable are found in native breeds of Partridged Penedesenca, Blonde Empordanesa hens [5], and Yarubi and Fulani ecotypes [68]. The thoracic perimeter variable is a good indicator of meat yield in most species of poultry [68]. Higher values were found than those obtained for Nigerian males and females [62]. However, native Nigerian birds bred for research purposes belonging to the Anak Titan ecotype [72] and the Yoruba and Fulani ecotypes [68] reported similar values to those found for Chimborazo individuals.
The Batsi Alak hens in Mexico reached values similar to those reported in Bolívar, Chimborazo and Cotopaxi [67]. These values were common in backyard hens in Mexico [63]. The measurements for half wing radius ulna length (hwrul) were shorter than those for the males and females of the Batsi Alak hen breed from Mexico [67]. However, the values for distal phalanx wing length (dpwl), were similar to those observed in Batsi Alak hens from Mexico for both males and females. This may be due to the fact that the study was carried out at an altitude of 1200 to 2760 m above mean sea level, similar conditions found at the location where our study took place [67]. Thigh length in Partridged Penedesenca and Blonde Empordanesa hens [70] were similar to those from Ecuadorian fowl from Cotopaxi and Tungurahua provinces, but shorter than the Nigerian Yoruba ecotype hens [68].
Regarding the circumference of the leg, the indigenous hens of Nigeria reported similar values to those measured in Bolívar and Chimborazo and slightly similar ones to those reported for the hens of Cotopaxi. The dimensions of the leg have been related in literature with the type of production, with those animals presenting higher dimensions (both in width and length), being more appropriate to suit the requirements for meat production and characteristic of carnic breeds [68]. Tanzanian local chicken ecotypes presented higher tarso-metatarsal lengths [5] for Singamagazi and Kuchi than those of Ecuadorian Creole hens. However, females of the same ecotypes reported similar average values [73] to those measured in Morona Santiago hens.
Similarly, males of the Ching'wekwe ecotype and females of the same ecotype and female of the Morogo ecotype presented similar values to those of the provinces of Guayas and Tungurahua, but much higher than those of the province of Cotopaxi. Likewise, this value was similar to that reported by Nigerian birds [62], which may support the fact that ecotype, may be strongly conditioned by the agroecological conditions in the area in which these avian population are based.
Similarly, larger sizes were recorded for breeds such as the Partridged Penedesenca and Blonde Empordanesa with an average of 8 cm, which may be because the birds studied also come from four different eco-climates, homologous to those of the areas considered in this research.
For birds from India [74], average values were similar to those measured in Bolívar and Chimborazo. Long tarsi have been associated with dry regions and flat topographies, as they allow birds to travel long distances in search for food, unlike birds with short tarsi, which could be attributed to the effects of natural selection [72].
Short tarsi have been identified with a greater ability to escape from predators [74], hence, they have been directly related to processes of adaptation and improvement of survival. Functionally, from a productive point of view, tall animals tend to be destined for meat production and small animals for egg production [68]. In addition, the length of tarsi may be related to the prediction of live weight in the field, as reported by some authors [75][76][77]. Ecuadorian creole hens present morphological traits which would made them more prone to produce eggs. However, the dimensions of certain morphological variables could make them suitable for meat production, so it is essential to implement breeding programs to select and direct crosses that allow ecotypes to be obtained by classifying individuals depending on to their productive potential.
Canonical discriminant analysis suggests that the variables perimeter of the leg (legc), metatarsus tarsus length (metl) and middle finger phalanx length (mfpl) were the ones that had the greatest discriminatory capacity between provinces. The results revealed the presence of wide ranges of variation within and among Creole hens in Ecuador. However, four large population blocs could be identified, namely Cotopaxi, Tungurahua, Morona Santiago, and that comprising the populations of Bolivar, Chimborazo and Guayas, a fact that could be attributed to the different conditions found across the various agroecological zones in the country, ethnic groups handling these resources and cultural implications that they have, along with the huge migration events suffered by these resources when facing natural and/or man-made challenges.
Toalombo, et al. [78] identified a common pattern of haplogroups of mitochondrial DNA for Ecuadorian chickens reared across Ecuadorian agroecological systems, which suggests these animals may belong to the same maternal lineages. In fact, the maternal origin for these populations could presumably be attributed to pre-Columbian Asiatic matrilines or Iberian matrilines arriving during the Spanish colonization. Furthermore, mitochondrial findings support the fact that current Ecuadorian local chickens do not show maternal influences from commercial lines and maintain high levels of genetic diversity without evidence of genetic drift and/or population bottlenecks. This high diversity may be owing to internal heterogeneity, which, as suggested by the results in the present study, may be promoted by the breeding policies that are carried out. Additionally, the patterns found by canonical discriminant analysis are supported by the results reported by Toalombo et al. [78], as a certain internal substructure can be found. However, the absence of any breeding program, registers, or zootechnic management produces a high fragmentation of the potential Ecuadorian breeds.
The most likely reason for which Morona Santiago, Bolívar, Chimborazo, and Guayas populations present a mixed structure which does not permit their complete segregation from one another may have its basis on the fact that poultry farmers in this province are prone to preserve their birds, avoiding the introduction of individuals from external populations, despite the fact that poultry production has already been developed in their areas.
The provinces of Morona Santiago, Bolívar, and Chimborazo are nearby to each other, sin contrast to the province of Guayas, located in the Costa Region. However, there is a provincial road that joins Chimborazo (Sierra) and Guayas (Costa). It should be noted that agricultural fairs are likely to take place along the road connecting both provinces. These events act as exchange centers of genetic material, which is mainly performed with minor species such as birds, given the considerable ease to transport such resources. The same would happen with Morona Santiago, which is located in the Ecuadorian Amazon, but near Chimborazo, which is geographically located in the center of the country. The provinces in question are characterized by their agricultural and livestock background, such that 13% of the population is engaged in poultry production activity; hence, it could be assumed that the populations maintain the genetics of their fowl over time.
The connection with Guayas, as stated by the General Secretariat of the Andean Community in 2009, lays on the fact that during the seventies, there was a reduction in domestic agricultural-livestock production, which led to the migration of the inhabitants of Nabuzo-Penipe (Chimborazo) to the coast of Ecuador. This social movement implied people carried easily-used animal species such as hens along with them. Complementarily, as a result of the eruptions of the Tungurahua volcano, constant since the late 1990s [79], the Canton Penipe experienced immigration from the surrounding populations, which was also noticeable in Nabuzo (Chimborazo), the area least affected by volcanic ash. This confirms that the genetic material did not suffer such migratory pressure and hence its resources remain intact, which in turn explains the clustering revealed by the canonical discriminant analysis and which may outline the same genetic structure.

Conclusions
The results revealed the presence of wide ranges of variation within and among Creole hens in Ecuador. However, four large population blocs could be identified, namely Cotopaxi, Tungurahua, Morona Santiago, and that comprising the populations of Bolivar, Chimborazo, and Guayas, a fact that could be attributed to the different conditions found across the various agroecological zones in the country, the ethnic groups handling these resources and cultural implications that they have, along with the huge migration events suffered by these resources when facing natural and/or man-made challenges. This addresses the high number of opportunities to implement programs for the promotion, conservation, and genetic improvement of these local resources, through their selection and crossing after their definition and characterization, since their inherent resistance and adaptation to the different environmental conditions allow the definition of technical and scientific strategies to exploit their productive potential. In addition, the local nature of these resources should be highlighted as intangible ancestral value in the field of food sovereignty, and should be inserted within the productive policies at a governmental level, as has already been done with other species, with the aim of improving rural livelihoods and meeting the growing demand for poultry products. Funding: This study received funding from "Instituto de Investigaciónes (IDI) -Vicerrectorado de Investigaciones -Escuela Superior Politcnica de Chimborazo" by the project proposal: "Caracterización morfológica y genética de la gallina criolla del Ecuador".