Vanilla planifolia Andrews (Orchidaceae): Labellum Variation and Potential Distribution in Hidalgo, Mexico

: Vanilla planifolia is a species of commercial importance. However, vanilla presents gene erosion problems due to its clonal reproduction. In the Huasteca of Hidalgo, there is no information on vanilla populations. Therefore, the objectives of this study were to identify the current populations and the potential distribution of, and the morphological variation in, the labellum of V. planifolia in the Huasteca of Hidalgo. Twenty-two accessions were located and selected. Based on 21 environmental variables, the niche modeling of the potential distribution was carried out with the MaxEnt program; with the Jackknife test being used to identify the variables that contributed to the model. Flowers from 22 accessions were collected and the labellum of each ﬂower was dissected. Subsequently, 64 morphological variables were obtained and various multivariate analyses were performed. The results showed three regions, deﬁned by the highest to the lowest probability that V. planifolia was distributed. The precipitation of the driest month, altitude, and vegetation cover delimited the distribution. Five different morphotypes were distinguished, and the main differences were associated with the middle part of the labellum as well as the entrance of pollinators to the ﬂower; therefore, the characterization of the labellum showed an infraspeciﬁc variation in V. planifolia in populations of the Huasteca of Hidalgo.


Introduction
Vanilla planifolia Andrews is a species of economic and ecological importance that is distributed from Mexico to Costa Rica [1,2]. Vanilla is native to Oaxaca and the crop was developed in the north of Veracruz, Mexico [1,3,4], although its cultivation has spread to different regions of the world. V. planifolia is a hemi-epiphytic or rupicolous plant that develops in evergreen or almost evergreen tropical forests, in primary or secondary vegetation at a height between 150 and 900 m above sea level (masl) [3,5]. It grows in evergreen or sub-evergreen forests with year-round rains on calcareous soils. In wetter areas, it can be found in young secondary forests. The flowers appear from March to April, and flowering is activated by low winter temperatures followed by warm temperatures in early spring [3]. Vanilla is subject to Special Protection (Pr) by the Mexican Government under NOM-059-SEMARNAT-2010, since there are only 30 registered collects in their natural environment [3]. Because Vanilla planifolia is propagated by cuttings, there are problems of genetic erosion in the crops [2,[6][7][8] and susceptibility to diseases (fungi and bacteria) [9].
Through the modeling of ecological niches and potential distribution, the environmental and anthropogenic variables that affect the distribution of a species can be identified [10][11][12][13], in addition to determining if there is gene flow between populations [14][15][16]. These analyses

Species Distribution Modeling
The potential distribution of Vanilla planifolia Andrews was modeled with 21 environmental variables: 20 variables of 30 s of the spatial resolution were obtained from the WorldClim database (www.worldclim.org, accessed on 18 May 2022) [56,57] and a vegetation cover variable was obtained from the CONABIO database [58] ( Table 1). The potential distribution was modeled using MaxEnt v. 3.4.1 [23,31,59]. The logistic output format was used to visualize the ideal habitat (probability of presence) of V. planifolia based on the different environmental variables [28,60]. The combined pixels in the model were recorded as the possible maximum entropy distribution space given by MaxEnt. Therefore, each cell on the map gives an estimate of the value of the suitability of the habitat on a scale that goes from 0 (least suitable) to 1 (most suitable) [23,31,61]. Table 1. Environmental variables used to obtain the potential distribution of V planifolia in the Huasteca of Hidalgo, Mexico.

Code
Environmental Variables Units Bio1 Annual mean temperature °C Bio2 Mean diurnal range °C Bio3 Isothermality Dimensionless Bio4 Temperature seasonality CV Bio5 Max temperature of the warmest month °C Bio6 Min temperature of the coldest month °C Bio7 Temperature annual range °C Bio8 Mean temperature of the wettest quarter °C Bio9 Mean temperature of the driest quarter °C Bio10 Mean temperature of the warmest quarter °C

Species Distribution Modeling
The potential distribution of Vanilla planifolia Andrews was modeled with 21 environmental variables: 20 variables of 30 s of the spatial resolution were obtained from the WorldClim database (www.worldclim.org, accessed on 18 May 2022) [56,57] and a vegetation cover variable was obtained from the CONABIO database [58] ( Table 1). The potential distribution was modeled using MaxEnt v. 3.4.1 [23,31,59]. The logistic output format was used to visualize the ideal habitat (probability of presence) of V. planifolia based on the different environmental variables [28,60]. The combined pixels in the model were recorded as the possible maximum entropy distribution space given by MaxEnt. Therefore, each cell on the map gives an estimate of the value of the suitability of the habitat on a scale that goes from 0 (least suitable) to 1 (most suitable) [23,31,61]. Table 1. Environmental variables used to obtain the potential distribution of V. planifolia in the Huasteca of Hidalgo, Mexico.

Bio1
Annual mean temperature • C Bio2 Mean diurnal range • C Bio3 Isothermality Dimensionless Bio4 Temperature seasonality CV Bio5 Max temperature of the warmest month • C Bio6 Min temperature of the coldest month • C Bio7 Temperature annual range • C Bio8 Mean temperature of the wettest quarter • C Bio9 Mean temperature of the driest quarter • C Bio10 Mean temperature of the warmest quarter • C Bio11 Mean temperature of the coldest quarter • C Bio12 Annual precipitation mm Bio13 Precipitation of the wettest month mm Bio14 Precipitation of the driest month mm Bio15 Precipitation seasonality CV Table 1. Cont.

Bio16
Precipitation of the wettest quarter mm Bio17 Precipitation of the driest quarter mm Bio18 Precipitation of the warmest quarter mm Bio19 Precipitation of the coldest quarter mm Cover Vegetation cover 16 types Alt Altitude m The accuracy of the model was evaluated by calculating the area under the curve (AUC) ROC (Receiver Operating Characteristic), which considers each value of the prediction result as a possible discrimination threshold and then calculates the corresponding sensitivity and specificity of each value. Sensitivity is the proportion of test localities that are present which were correctly predicted (1-extrinsic omission rate). The quantity (1-specificity) is the proportion of all of the pixels predicted to have suitable conditions for the species [23]; therefore, based on the AUC value, the model can be considered as poor (AUC < 0.8), fair (0.8 ≤ AUC < 0.9), good (0.9 ≤ AUC < 0.95), or very good (0.95 ≤ AUC < 1.0) [25].
Subsequently, a Jackknife test [62] was carried out, which allows for analyzing the contribution of each environmental variable individually and jointly, to form the distribution model of V. planifolia [63]. Therefore, through this test and the percentage of contribution of the species, the variables that locate the potential distribution of V. planifolia in the Huasteca of Hidalgo, Mexico were determined.

Morphological Characterization of the Flower
In April 2013, during the flowering season, 328 turgid flowers with pollinia and no apparent damage to the floral structure were collected from 22 Vanilla planifolia Andrews populations (22 accessions) ( Figure 2A). The petals, sepals, and labellum were dissected and stored in a preservative solution (27% ethanol, 4% lactic acid, 3% benzoic acid, 3% glycerin, and 63% distilled water) inside 220 mL bottles with their respective collection label.

Bio11
Mean temperature of the coldest quarter °C Bio12 Annual precipitation mm Bio13 Precipitation of the wettest month mm Bio14 Precipitation of the driest month mm Bio15 Precipitation seasonality CV Bio16 Precipitation of the wettest quarter mm Bio17 Precipitation of the driest quarter mm Bio18 Precipitation of the warmest quarter mm Bio19 Precipitation of the coldest quarter mm Cover Vegetation cover 16 types Alt Altitude m The accuracy of the model was evaluated by calculating the area under the curve (AUC) ROC (Receiver Operating Characteristic), which considers each value of the prediction result as a possible discrimination threshold and then calculates the corresponding sensitivity and specificity of each value. Sensitivity is the proportion of test localities that are present which were correctly predicted (1-extrinsic omission rate). The quantity (1specificity) is the proportion of all of the pixels predicted to have suitable conditions for the species [23]; therefore, based on the AUC value, the model can be considered as poor (AUC < 0.8), fair (0.8 ≤ AUC < 0.9), good (0.9 ≤ AUC < 0.95), or very good (0.95 ≤ AUC < 1.0) [25].
Subsequently, a Jackknife test [62] was carried out, which allows for analyzing the contribution of each environmental variable individually and jointly, to form the distribution model of V. planifolia [63]. Therefore, through this test and the percentage of contribution of the species, the variables that locate the potential distribution of V. planifolia in the Huasteca of Hidalgo, Mexico were determined.

Morphological Characterization of the Flower
In April 2013, during the flowering season, 328 turgid flowers with pollinia and no apparent damage to the floral structure were collected from 22 Vanilla planifolia Andrews populations (22 accessions) ( Figure 2A). The petals, sepals, and labellum were dissected and stored in a preservative solution (27% ethanol, 4% lactic acid, 3% benzoic acid, 3% glycerin, and 63% distilled water) inside 220 mL bottles with their respective collection label. To identify morphometric variation, the work was based on the geometric morphometry of contours that are used in the analysis of anatomical structures. The shape of a To identify morphometric variation, the work was based on the geometric morphometry of contours that are used in the analysis of anatomical structures. The shape of a structure is analyzed as a whole and not in fragments [1], in addition to characterizing the shapes through multivariate analysis [64,65].
The procedure to characterize the labellum was based on Hernández-Ruíz et al. [39] and Lima-Morales et al. [42]. First, the labellum stored in a preservative solution was stained with methylene blue (0.08%) ( Figure 2B,C). Photographs were taken at a distance of 30 cm with a Sony digital camera (SONY α, DSLR-SLT-A55) equipped with a macro lens (Sony DT 30 mm F/2.8 SAM). Once the digital images of all of the flowers were obtained, the initial landmark points were placed ( Figure 3A). In the curved regions, extra points were added without overloading the contour edges so as not to generate redundant information [66]. With the first landmark points, the labellum was separated into seven regions: A, B, C, D, E, F, and G ( Figure 3B). Then, the secondary lines were placed to record the entire shape of the labellum; thus, a total of 57 lines and 7 angles (morphological variables) were obtained ( Figure 3C).
The procedure to characterize the labellum was based on Hernández-Ruíz et al. [39] and Lima-Morales et al. [42]. First, the labellum stored in a preservative solution was stained with methylene blue (0.08%) ( Figure 2B,C). Photographs were taken at a distance of 30 cm with a Sony digital camera (SONY α, DSLR-SLT-A55) equipped with a macro lens (Sony DT 30 mm F/2.8 SAM). Once the digital images of all of the flowers were obtained, the initial landmark points were placed ( Figure 3A). In the curved regions, extra points were added without overloading the contour edges so as not to generate redundant information [66]. With the first landmark points, the labellum was separated into seven regions: A, B, C, D, E, F, and G ( Figure 3B). Then, the secondary lines were placed to record the entire shape of the labellum; thus, a total of 57 lines and 7 angles (morphological variables) were obtained ( Figure 3C).

Statistical and Numerical Analysis
For all of the labellum lines and angles, the mean and the coefficient of variation were obtained. Subsequently, an analysis of variance under a completely random unbalanced scheme was performed to determine if there were significant differences between the accessions. The 22 accessions were considered to be the source of variation; therefore, each collection had a different number of replicates (Table 2). Subsequently, a multivariate analysis of Principal Components and a hierarchical cluster analysis based on the Euclidean distance of each mean were performed to identify infraspecific variation in the labellum of Vanilla planifolia Andrews in the Huasteca of Hidalgo using the Software SAS 9.1 (SAS Institute, Cary, NC, USA) and the JMP 10.0.2 (SAS Institute, Cary, NC, USA).

Statistical and Numerical Analysis
For all of the labellum lines and angles, the mean and the coefficient of variation were obtained. Subsequently, an analysis of variance under a completely random unbalanced scheme was performed to determine if there were significant differences between the accessions. The 22 accessions were considered to be the source of variation; therefore, each collection had a different number of replicates (Table 2). Subsequently, a multivariate analysis of Principal Components and a hierarchical cluster analysis based on the Euclidean distance of each mean were performed to identify infraspecific variation in the labellum of Vanilla planifolia Andrews in the Huasteca of Hidalgo using the Software SAS 9.1 (SAS Institute, Cary, NC, USA) and the JMP 10.0.2 (SAS Institute, Cary, NC, USA).  In the Huasteca of Hidalgo, 22 accessions of Vanilla planifolia Andrews were located in the municipalities of Atlapexco, Jaltocán, and Huejutla de Reyes (Table 3). These sites presented the conditions of vanilla plantations that were more than 20 years old and in traditional systems of acahuales and Coffea arabica Benth plantations under the shade of native trees (for example Pimenta dioica (L.), Bursera Jacq. ex L. spp., and Ceiba pentandra (L.) Gaertn.) and minimal management. Sites with intensive management and the recent acquisition of cuttings were excluded. Agricultural use S12 367 S13 331 Pahuatlán S14 381 Ichcatepec S15 545 The municipality of Huejutla presented the largest number of V. planifolia populations with 63.6% of the total, while Atlapexco had 27.2% and Jaltocán had 9.2% ( Figure 4). The vanilla populations were located between 273 and 545 masl; 31.8% of the accessions were in a warm humid climate, 45.4% in a humid semi-warm climate, 9.1% in a humid warm climate with the coldest month less than 18 • C, and 13.7% in a humid semi-warm climate of group C (Table 3).

Potential Distribution
The MaxEnt model predicted the potential distribution of Vanilla planifolia Andrews with a training area under the curve (AUC) of 0.985 ( Figure 5A). The red curve (indicating the degree of fit of the sampling data) and the blue one (indicating the fit of the model)

Potential Distribution
The MaxEnt model predicted the potential distribution of Vanilla planifolia Andrews with a training area under the curve (AUC) of 0.985 ( Figure 5A). The red curve (indicating the degree of fit of the sampling data) and the blue one (indicating the fit of the model) were identical and the values were considered to be acceptable. Figure 5B shows the omission rate calculated on both the training presence records and the test records. In the omission rate, a small part fell below the predictions, and another remained above the predictions because the sample used and the training samples were dependent. Once the modeling carried out by MaxEnt was validated, the potential distribution of V. planifolia was obtained ( Figure 6). The 22 accessions were located only in the northern and northwestern parts of Hidalgo that belong to the region of the Huasteca of Hidalgo. The areas in red showed the highest probability of presenting populations of V. planifolia. By contrast, the green areas are the places where there is little probability of finding a population of V. planifolia ( Figure 6). Once the modeling carried out by MaxEnt was validated, the potential distribution of V. planifolia was obtained ( Figure 6). The 22 accessions were located only in the northern and northwestern parts of Hidalgo that belong to the region of the Huasteca of Hidalgo. The areas in red showed the highest probability of presenting populations of V. planifolia. By contrast, the green areas are the places where there is little probability of finding a population of V. planifolia ( Figure 6).
The populations marked by GPS were divided into three groups that were differentiated by the probability of finding V. planifolia. In Group I (GI), the largest number of populations was present (red area, 67-100% probability); therefore, it is the area with adequate environmental conditions for the development of V. planifolia in the Huasteca of Hidalgo. Group II (GII) was located in the orange area, with a 51-66% probability that the V. planifolia populations were distributed there. Finally, in Group III (GIII), the probability of finding populations of V. planifolia was 34-50% ( Figure 6). The areas in gray were the areas where V. planifolia is not distributed and cannot be cultivated. The populations marked by GPS were divided into three groups that were differentiated by the probability of finding V. planifolia. In Group I (GI), the largest number of populations was present (red area, 67-100% probability); therefore, it is the area with adequate environmental conditions for the development of V. planifolia in the Huasteca of Hidalgo. Group II (GII) was located in the orange area, with a 51-66% probability that the V. planifolia populations were distributed there. Finally, in Group III (GIII), the probability of finding populations of V. planifolia was 34-50% ( Figure 6). The areas in gray were the areas where V. planifolia is not distributed and cannot be cultivated.

Effect of Environmental Variables
The variables that contributed the most to the potential distribution model generated by MaxEnt were the precipitation of the driest month (Bio14) (43%) and vegetation cover (Cover) (14.9%). Therefore, the two variables Bio14 and Cover were the determinants to generate the potential distribution model of Vanilla planifolia Andrews in the Huasteca of Hidalgo (Table 4).

Effect of Environmental Variables
The variables that contributed the most to the potential distribution model generated by MaxEnt were the precipitation of the driest month (Bio14) (43%) and vegetation cover (Cover) (14.9%). Therefore, the two variables Bio14 and Cover were the determinants to generate the potential distribution model of Vanilla planifolia Andrews in the Huasteca of Hidalgo (Table 4). The environmental variables that were individually most important for the potential distribution of V. planifolia were precipitation of the driest month, precipitation seasonality, precipitation of the coldest quarter, altitude, precipitation of the driest quarter, and temperature seasonality (Figure 7). The least important variables, individually, for the potential distribution of V. planifolia were temperature annual range and vegetation cover, which individually do not present a direct effect but, if they are eliminated, affect the distribution of the model when analyzed with the other variables together, as well as the mean diurnal temperature range (Figure 7).
Annual mean temperature (Bio1) 4.9 Mean diurnal range (Bio2) 2.7 Altitude (Alt) 1.3 Precipitation of the wettest quarter (Bio16) 0.7 Temperature annual range (Bio7) 0.6 The environmental variables that were individually most important for the potential distribution of V. planifolia were precipitation of the driest month, precipitation seasonality, precipitation of the coldest quarter, altitude, precipitation of the driest quarter, and temperature seasonality (Figure 7). The least important variables, individually, for the potential distribution of V. planifolia were temperature annual range and vegetation cover, which individually do not present a direct effect but, if they are eliminated, affect the distribution of the model when analyzed with the other variables together, as well as the mean diurnal temperature range (Figure 7).

Labellum Characterization
The 64 morphological variables which were analyzed presented low coefficients of variation (3-10%). In addition, highly significant differences were observed between the accessions for each of the variables (Table 5).

Diversity Clustering
The multivariate analysis of the clustering showed that, on a Euclidean distance of 0.831 in the dendrogram of Figure 9, the five morphotypes were confirmed for Vanilla planifolia Andrews accessions from the Huasteca region of Hidalgo. Similar tones mean variables with similar values, in addition to the fact that intense blue tones show the highest values while white tones represent the lowest values. Morphotypes I and II differed from Morphotypes III, IV, and V at a distance of 1.008 because the variables that represented the shape of the mid-apical and mid-basal region were the ones that presented the most information and mainly served to differentiate the morphotypes in the dendrogram. Subsequently, Morphotype I was separated from Morphotype II by the angles of the labellum (aA, aB, aD, aE, and aG). The Morphotypes III, IV, and V were separated by the angles of the labellum and the basal region (aA, aB, aD, aE, aG, A2, A3, A4, and A5) ( Figure 9A).

Diversity Clustering
The multivariate analysis of the clustering showed that, on a Euclidean distance of 0.831 in the dendrogram of Figure 9, the five morphotypes were confirmed for Vanilla planifolia Andrews accessions from the Huasteca region of Hidalgo. Similar tones mean variables with similar values, in addition to the fact that intense blue tones show the highest values while white tones represent the lowest values. Morphotypes I and II differed from Morphotypes III, IV, and V at a distance of 1.008 because the variables that represented the shape of the mid-apical and mid-basal region were the ones that presented the most information and mainly served to differentiate the morphotypes in the dendrogram. Subsequently, Morphotype I was separated from Morphotype II by the angles of the labellum (aA, aB, aD, aE, and aG). The Morphotypes III, IV, and V were separated by the angles of the labellum and the basal region (aA, aB, aD, aE, aG, A2, A3, A4, and A5) ( Figure 9A).
Based on PC1 (Figure 8), the morphotypes located on the positive side of the graph were more elongated and broader in the mid-basal and basal region of the labellum (Morphotype I and Morphotype II), while those on the negative side were narrower in this region (Morphtype III, Morphtype IV, and Morphtype V) (Figure 8). In PC2, on the negative side, the morphotypes were thin in the left part of the middle region of the labellum (D2) and the basal region (B1, aA, and aB), but long in the middle basal region and the left part of the middle region of the labellum (B6, B7, and D1) (Morphotypes II, III, and IV). Morphotype I and V were located in the middle part of PC2 (Figure 8).
For PC3, the labellum on the positive side was longer in the right part of the middle region (D, D4, and D8) and the left part of the apical middle region of the labellum (E1). The labellum was narrower on the right side of the mid-region (D3 and aD) and short on the right side of the apical mid-region (E4). However, only Morphotype III was found in the positive part, while the other morphotypes were in the middle of PC3 (Figure 8). The morphological expression profile shows the behavior of each variable for each accession ( Figure 9B). The variables with similar behaviors were grouped so that the five morphotypes were separated based on the differences of each accession. The height of the peaks corresponds to the value of each variable, and high peaks represent high values; therefore, Morphotype I is the one with the highest peaks for each accession, while it is the opposite for Morphotype V ( Figure 9B).

Potential Distribution Model
The potential distribution model of Vanilla planifolia Andrews identified three regions where there was a higher probability of finding this species and that could be considered as conservation areas for the germplasm present in the Huasteca region of Hidalgo [19,22]. This measure will prevent the disappearance of vegetation and changes in land use that reduce the potential distribution areas, as has been reported for V. planifolia in Oaxaca [28], San Luis Potosí [29,30], and Mexico in general [31].
The value obtained to validate the distribution model of Vanilla planifolia was 0.994, which means that the model prediction of the potential distribution of V. planifolia was acceptable (0.95 ≤ AUC < 1.0) because the current test data are adjusted with the training data [25]. The AUC value is high because V. planifolia is localized to specific environmental conditions [5], while in species that are located in different environments, the AUC value tends to be lower [20,67,68].
The evaluation of the model allows us to know its usefulness; therefore, it must be validated to know if the results are significant. To this end, the omission (or commission) rate is used, which is a binomial test that is dependent on a threshold based on omission and predicted area [23,69,70]. The omission rate is the fraction of test locations that fall into pixels that are not expected to be suitable for V. planifolia; the predicted area is the fraction of all of the pixels that are predicted to be suitable for the species [25]. In Figure 5B, the omission in the test examples was adjusted to the predicted omission rate, which is the omission rate for the test data modeled from the distribution given by MaxEnt (the omission rate predicted is a straight line due to the cumulative output format). Thus, the potential distribution modeled by MaxEnt was validated since the omission of the test data is close to the predicted omission [23,25,69].

Environmental Variables That Define the Potential Distribution
The main environmental variable that defined the potential distribution of Vanilla planifolia Andrews in the Huasteca of Hidalgo was the precipitation of the driest month. The abundance of rain in the month of April is a determinant for the establishment of vanilla populations ( Figure 10B), similar to what Armenta-Montero et al. [31] report. Trinidad-García et al. [30] and Reyes-Hernández et al. [29] mention that the total precipitation in the Huasteca Potosina is one of the main factors that influence the distribution of vanilla; for Hernández-Ruíz et al. [28] it was the month with the highest rainfall in Oaxaca. However, because it is a species under cultivation, Soto-Arenas and Cribb [5] found that V. planifolia is established in dry conditions in the spring months for Veracruz.
In the case of altitude, there is a greater probability of finding populations of Vanilla planifolia at lower altitudes than at higher altitudes for the Huasteca region of Hidalgo ( Figure 10A). These values fall within the range that has been reported for other populations of V. planifolia that are distributed from 150 to 800 masl [3,5]; in Oaxaca they are located from 200 to 1190 masl [28], and in the Huasteca Potosina, from 61 to 678 masl [29].
The climate in which the populations of the Huasteca of Hidalgo occur is similar to the Totonacapan region conditions (Puebla-Veracruz), where there are populations in warm humid and warm sub-humid conditions [71]. The vegetation cover did not turn out to be a decisive factor because, individually, it does not affect the distribution. However, based on the Jackknife test together with the other variables, it provides information in the construction of the model. abundance of rain in the month of April is a determinant for the establishment of vanilla populations ( Figure 10B), similar to what Armenta-Montero et al. [31] report. Trinidad-García et al. [30] and Reyes-Hernández et al. [29] mention that the total precipitation in the Huasteca Potosina is one of the main factors that influence the distribution of vanilla; for Hernández-Ruíz et al. [28] it was the month with the highest rainfall in Oaxaca. However, because it is a species under cultivation, Soto-Arenas and Cribb [5] found that V. planifolia is established in dry conditions in the spring months for Veracruz. In the case of altitude, there is a greater probability of finding populations of Vanilla planifolia at lower altitudes than at higher altitudes for the Huasteca region of Hidalgo (Figure 10A). These values fall within the range that has been reported for other populations of V. planifolia that are distributed from 150 to 800 masl [3,5]; in Oaxaca they are located from 200 to 1190 masl [28], and in the Huasteca Potosina, from 61 to 678 masl [29].
The climate in which the populations of the Huasteca of Hidalgo occur is similar to the Totonacapan region conditions (Puebla-Veracruz), where there are populations in warm humid and warm sub-humid conditions [71]. The vegetation cover did not turn out to be a decisive factor because, individually, it does not affect the distribution. However, based on the Jackknife test together with the other variables, it provides information in the construction of the model. The largest number of populations (59.1%) was in the type of agricultural use cover. Vanilla planifolia was located in secondary vegetation, made up of acahuales for timber or Coffea arabica Benth plantations, and similar to that reported in the areas of Veracruz and part of Puebla [4,34,71]. Further, 40.9% of the accessions were associated with tropical or subtropical evergreen broadleaf forests, which are mainly used for Coffea arabica plantation, and the cultivation of Pimenta dioica (L.) Merr., Ceiba pentandra (L.) Gaertn, and Pouteria sapota (Jacq.) H. E. Moore and Stearn (direct observation in the field). These conditions are similar to some acahuales in Puebla and Veracruz [4,34,71], which serve as reservoirs for native species [72][73][74].
The use of computational predictive models has allowed for the identification of the distribution of species through the analysis of the environmental conditions of the sites where they are collected [75,76]. Geographic Information Systems, together with predictive algorithms, allow for the modeling of ecological niches. These models constitute an important technique in analytical biology which is oriented mainly to the conservation and management of species [10,24,63]. The distribution and geographical area of the plants are influenced by two main factors: altitude and climate [77,78]. However, plant species adapt to variations in environmental conditions [79,80] with significant changes in the composition and structure of populations; therefore, there may be species with a wide or very restricted distribution (endemic) [35,36]. In the case of Vanilla planifolia, the main factors were the precipitation of the driest month and altitude. In addition, because its distribution is restricted, it is considered to be an endemic species to Mexico [81,82].

Labellum Morphotypes
Generally, biotic and abiotic factors influence the morphological variation in vegetative and reproductive characters [83]. Within the reproductive characters, the flowers present quantitative variation in the populations of the same species. This variation represents the basis of natural selection that can eventually result in diversification and speciation [46,84]. The size of the flower in some species is modified due to environmental variation; therefore, they become larger at high altitudes, cold temperatures, and high humidity, and shorter at low altitudes, warm temperatures, and dry conditions [46,85,86]. However, these variations occur due to the plasticity that individuals present under different environmental conditions [84,87], as reported for Arabidopsis thaliana (L.) Heynh [45], Narcissus triandrus L. [85], and Campanula rotundifolia Pall. Ex Roem. and Schult [86].
In many orchids, there is a high degree of pollinator specificity in flower shape; therefore, as they are specialized in pollination, the variation is minimal within the same species [48].
The Vanilla planifolia populations analyzed in this work had significant differences in all of the morphological variables. However, through the PCA, the variables that allowed for the separation of the populations into five different morphotypes were identified.
The traits that separated the groups for PC1 corresponded to the basal regions and the middle region of the labellum, similar to those reported for Vanilla planifolia Andrews in Oaxaca [40] and San Luis Potosí [42], while for Vanilla. pompona Schiede, they were from the middle and basal region [41].
In the case of PC2, the morphological variables that influence the separation of the groups corresponded to the middle regions and left part of the callus region, similar to what Hernández-Ruíz et al. [40] reported for V. planifolia from Oaxaca, but differing from what Lima-Morales et al. [2] reported in San Luis Potosí. In PC3, the morphological variables that separated the groups in this study corresponded exclusively to the callus region, a situation that coincides with V. planifolia from San Luis Potosí [2], and partially coincides with what they reported for Oaxaca for V. planifolia [40]. Considering the three PCs and comparing them with V. planifolia from Oaxaca [40], San Luis Potosí [42], and V. pompona [41], the relevant regions are the middle part and the callus. These regions define the entrance of the pollinating insect to the flower as suggested by Hernández-Ruíz et al. [41] for V. pompona and confirmed by Hernández-Ruiz et al. [40] and Lima-Morales et al. [42] for V. planifolia.
Although the variables of each PC are independent between each PC, as observed in Figure 8, the important variables are located in the same regions and are closer to each other; therefore, even though they are independent in the multivariate analysis, they are directly related to the structure of the labellum, a situation that was not observed in the published articles on Vanilla from San Luis Potosí and Oaxaca [40,42].
In addition, through the hierarchical clustering heatmap and the morphological expression profiles of the labellum of V. planifolia (Figure 9), the areas that varied between the five morphotypes were obtained: Morphotype I (MI). Represented by nine accessions, MI had a wide labellum in the basal and mid-basal region, elongated on the right side, broad on the left side of the middlemiddle apical region, and larger in the region of the apical lobes. It is the largest labellum compared to the other morphotypes ( Figure 9B, Table A1).
Morphotype II (MII). Represented by three accessions. The main characteristic of MII was the larger basal region where the labellum joins the base of the column. For the mid-basal region, the structure was more elongated in a similar way to the median and apical median region of the labellum ( Figure 9B, Table A1).
Morphotype III (MIII). With only two accessions, MIII had the third largest labellum size and was intermediate between Morphotypes I and II (larger) and Morphotypes IV and V (smaller labellum size) ( Figure 9B, Table A1).
Morphotype IV (MIV). Represented by seven accessions, MIV had a small labellum (only surpassed by Morphotype V) in the basal region and the middle region ( Figure 9B, Table A1).
Morphotype V (MV). This morphotype had only one accession, and presented the smallest labellum size, mainly in the mid-basal region, apical mean, and apical lobes ( Figure 9B, Table A1).
Morphotypes I and II were those with the largest labellum sizes, followed by Morphotype III, then Morphotype IV, and the one with the smallest labellum was Morphotype V.
These five morphotypes would represent the vanilla populations that develop in the Hausteca of Hidalgo because, as previously mentioned, the 22 accessions analyzed come from acahuales and coffee plantations with little or no management, the age of the plants is more than 20 years old, and they have not been recently acquired from other regions such as Papantla, Veracruz.

Geographic Distribution of Morphotypes
Reproductive characters generally show a certain degree of morphological variation, a product of the genetic variability of each species and on which natural selection acts as the main force of speciation [40]. In situations where morphological variation is associated with environmental factors, it has been documented that it is generally expressed as gradual or mosaic patterns across a landscape or geographic area [83]. When the phenotype of a plant is affected by any of these factors, environmental patterns can be treated as geographic patterns of phenotypic variation [83,88]. Particularly, this type of variation is related to species with a wide geographical distribution, and which occupy discontinuous territories in the form of mosaics [89].
However, in the case of Vanilla planifolia Andrews, the distribution of the five morphotypes was not associated with abiotic factors: the five morphotypes were distributed within the same soil moisture regime (Udic from 270 to 330 days of moisture), at the same elevation that goes from 250 to 556 masl, and in areas with a total annual rainfall of 1000-2000 mm and an average temperature of 21-23 • C. Therefore, vanilla had no climatic pattern since the same morphotype was distributed in several types of climates, as reported by Soto-Arenas and Solano-Gómez [82], and as seen in the type of vegetation of the Huasteca of Hidalgo ( Figure 11) [58]. 2000 mm and an average temperature of 21-23 °C. Therefore, vanilla had no climatic pattern since the same morphotype was distributed in several types of climates, as reported by Soto-Arenas and Solano-Gómez [82], and as seen in the type of vegetation of the Huasteca of Hidalgo ( Figure 11) [58].

Final Considerations on the Labellum Variation
Since the presence of five morphotypes cannot be explained by environmental variables, other factors such as biotic factors could be considered. McCormirck and Jacquemyn [90] suggest that micro factors such as mycorrhizae, tutors, and pollinators are factors that can affect and modify the spatial distribution of orchids in general. Damon et al. [91] sug-

Final Considerations on the Labellum Variation
Since the presence of five morphotypes cannot be explained by environmental variables, other factors such as biotic factors could be considered. McCormirck and Jacquemyn [90] suggest that micro factors such as mycorrhizae, tutors, and pollinators are factors that can affect and modify the spatial distribution of orchids in general. Damon et al. [91] suggest that the distribution and abundance of euglossine bees (Euglossini Latreille) in agroecosystems and forest fragments in southern Mexico is associated with relict forests and coffee plantations due to light and humidity, conditions that occur in the Huasteca of Hidalgo.
Shipunov and Bateman [49] pointed out that the size and shape of the labellum are important factors for the attraction of pollinators. Benítez-Vieyra et al. [92] reported the same situation for the orchid Geoblasta pennicillata (Rchb. f.) Hoehne ex M.N Correa, which attracts its pollinator by having a labellum shaped like the female wasps of the species Campsomeris bistrimacula Lepeletier. In Cryptostylis R. Br. orchids, the larger labellum functions as a stimulant for pollinating wasps; therefore, in some orchid species, the labellum is under constant selection pressure from pollinators [93].
The differences between Vanilla planifolia morphotypes were concentrated in the shape of the labellum attraction zone and were exposed to selection by pollinators ( Figure 12). In the Huasteca region of Hidalgo, the vanilla plantations in acahual and with little management of the crop present natural pollination due to the presence of wasps and bees of the Euglossa Latreille and Eulaema Lepeletier genera [32,94]. Pollinators influence the variation in the size, shape, and color of floral structures in some plant species [44,95,96]; therefore, they are one of the main causes of floral evolution [45,46]. The morphological variation in the flowers depends on the level of specialization with the pollinator (the case of some orchids) [97]; for this reason, the variation in size of the flower is minimal because there is a strong relationship between the pollinator and the flower that is stable in climatic variations [44,93,98,99]. In addition, the labellum is a very important organ, not only to identify and differentiate highly related taxonomic entities [100], but also to study the processes and mechanisms that generate variation and adaptation within and between populations of V. planifolia, as has been reported for other regions of Mexico [40,42].  [32,94]. Pollinators influence the variation in the size, shape, and color of floral structures in some plant species [44,95,96]; therefore, they are one of the main causes of floral evolution [45,46]. The morphological variation in the flowers depends on the level of specialization with the pollinator (the case of some orchids) [97]; for this reason, the variation in size of the flower is minimal because there is a strong relationship between the pollinator and the flower that is stable in climatic variations [44,93,98,99]. In addition, the labellum is a very important organ, not only to identify and differentiate highly related taxonomic entities [100], but also to study the processes and mechanisms that generate variation and adaptation within and between populations of V. planifolia, as has been reported for other regions of Mexico [40,42]. Possibly, the variation in the morphology of the labellum of the populations of the Huasteca Hidalguense are mainly related to the size of the pollinator; however, it is necessary to carry out studies on the biological interactions between plant and pollinator [101], Figure 12. Labellum morphological variables exposed to selection by pollinators of V. planifolia.
Possibly, the variation in the morphology of the labellum of the populations of the Huasteca Hidalguense are mainly related to the size of the pollinator; however, it is necessary to carry out studies on the biological interactions between plant and pollinator [101], as in other species of Vanilla [102], to determine if the shape between the morphotypes is related to the size of the pollinators present in the region or if the variation corresponds to other environmental factors as a product of the plasticity or accommodation of the plants to the environment in which they develop. In addition, to confirm that the morphological variation reflects the genetic diversity of V. planifolia, analyses with molecular markers must be carried out to characterize the germplasm [103,104] and propose generic improvement programs (new hybrids or new varieties) [105] and conservation programs in the Huasteca de Hidalgo, Mexico.

Conclusions
In the Huasteca region of Hidalgo, Mexico, 22 accessions of V. planifolia were located in acahuales and Coffea arabica Benth plantations with native vegetation and minimal management. The potential distribution map shows that, based on the probability of presence, the populations of V. planifolia were located in three groups from the highest to lowest probability of the presence of vanilla. The main environmental variables that delimit the potential distribution of V. planifolia in the Huasteca of Hidalgo were precipitation and altitude. In addition, five different labellum morphotypes which were possibly related to plant-pollinator interaction were obtained. However, it is necessary to deepen the study of the morphology associated with the floral ecology of the germplasm of V. planifolia in the Huasteca of Hidalgo.  Acknowledgments: Thanks to the farmers of Huasteca of Hidalgo, especially to Manuel Ambrocio for his help.

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