Growing Out of the Tropical Forests: Gene Flow of Native Mesoamerican Trees Among Forest and Mayan Homegardens

This work aimed to evaluate domestication effects on the genetic structure of two dioecious species Brosimum alicastrum Sw. (Moraceae) and Spondias purpurea L. (Anacardiaceae), and a heterostylous one Cordia dodecandra A. DC. (Cordiaceae), growing in remnant forests and homegardens within two climatic regions of the Peninsula of Yucatan. The trees of B. alicastrum and C. dodecandra are propagated by seeds in both population types, while those of S. purpurea are propagated asexually in the homegardens. ISSRs genetic markers were amplified from foliar tissue of 18 to 21 plants per population type/region combination for each species. Genetic diversity, genetic differentiation, and genetic structure estimators were obtained and compared among species at the regional and population level. We found higher polymorphism (37.5–41), but lower private alleles (4–4.4) and similar heterozygosity (0.1–0.12) in the species with sexual reproduction compared to S. purpurea (34, 8, and 0.11, respectively). Genetic diversity in B. alicastrum populations varied with the region; in C. dodecandra, to the population type; and in S. purpurea, to both the population type and the region. Unrestricted gene flow among regions was suggested by low ΦRT in C. dodecandra and S. purpurea (−0.006 and 0.002) but not for B. alicastrum (0.1). Gene flow between populations within the regions for the sexually reproducing species was suggested by lower θII (0.005–0.07 and 0.008–0.1) estimates than those of S. purpurea (0.09 and 0.13). Even though the lowest paired FST (0.002–0.05) and ΦST (0.002–0.12) values were found between the northeastern forest and homegarden populations for the three species, the dendrogram, Bayesian assignment, and K-Means analyses suggest that the least differentiated populations are southwestern forest and homegarden populations of B. alicastrum and S. purpurea, and the southwestern forest and northeastern homegarden of C. dodecandra. The sexual reproduction, biotic interactions, and extensive management of B. alicastrum and C. dodecandra in the agroforestry and the urban systems may contribute to connectivity between wild and domesticated populations, while in S. purpurea this connectivity is interrupted by the clonal propagation of the species in the homegardens.


INTRODUCTION
Domesticated species diverge from wild populations because of isolation, selection, and the founder effect. Selection and preservation of a few individuals in their native habitats can lead to bottleneck events and the establishment of a few individuals in traditional agroecosystems lead to founder events (Casas et al., 2007;Miller and Gross, 2011). The founder effect reduces the genetic diversity of domesticated populations, which in turn can increase the genetic differentiation between domesticated populations and wild populations (Miller and Gross, 2011). Besides, conscious or unconscious domestication can cause changes in domesticated plants phenology, in the communities of pollinators and/or fruit dispersors, which interrupts genetic flow between populations under cultivation and wild populations (Zohary, 2004;Ross-Ibarra et al., 2007). However, in tropical species in which wild populations have decreased drastically as a consequence of deforestation, cultivated populations may have higher genetic diversity and the native trees that are kept in the agroecosystems may even form bridges that maintain gene flow between the remnant fragments of the forests (Duminil et al., 2007).
Gene flow and genetic diversity depend primarily on the reproductive system of the species and the dispersal ability of pollen, and to a lesser extent, on seeds dispersal (Loveless, 1992;Hamrick and Got, 1996;Nybom, 2004;Duminil et al., 2007;Gamba and Muchhala, 2020). Sexual reproduction of plants, when accompanied by mechanisms that promote crosspollination, such as dioecy and self-incompatibility, increase genetic diversity and gene flow among populations, mainly by its effect on the isolation of the population, the effective population size (N e ), and recombination rates (reviewed by Glémin et al., 2019). In plants with clonal propagation, genetic diversity can be high and gene flow low because of increased isolation and lower N e , but with high heterozygosity, because the emergence of new mutations produces a heterozygous genotype in diploid organisms (Glémin et al., 2019). Plant species pollinated by the wind, vertebrates, and large insects disperse pollen at larger distances than those pollinated by small insects. Hence, the latter often have a higher genetic population structure (Gamba and Muchhala, 2020); in plants with seeds dispersal by gravity, a stronger genetic structure is expected (Loveless, 1992;Hamrick and Got, 1996).
Genetic population structure in plants is associated with lifehistory traits that affect gene exchange among individuals and populations (Loveless, 1992;Hamrick and Got, 1996;Nybom, 2004;Duminil et al., 2007). The expected genetic structure of tropical trees populations derives from the effect of the species low densities species which produce a smaller effective population size (N e ); the dependency on small insects as pollination vectors; and the restricted seed dispersal that results in lower gene flow and can promote self-pollination and biparental inbreeding at small rates significantly different from those of the trees of temperate systems (Dick et al., 2008;Duminil et al., 2009).
At the Yucatan Peninsula, the fragments of the deciduous low height forest, medium height sub-deciduous and semi-evergreen forests that made up the original vegetation coexist within a vegetation matrix for agricultural and livestock production (Porter-Bolland et al., 2015). In this vegetation matrix, people in rural and urban settlements keep homegardens, traditional agroecosystems-also know and as solares Mayas-around the family house in which native species are preserved, cultured, and used (Ramírez- Barajas and Calmé, 2015). The solares Mayas are considered sites of domestication of many plants, that could have been managed by the Yucatec Maya since pre-Hispanic times (Colunga-García Marín and Zizumbo-Villarreal, 2004;Zizumbo-Villarreal and Colunga-García Marín, 2010).
Peninsular homegardens are comprised of a group of species, generally fruit trees and shrubs. This group of species (that has between 12 and 32 species) form a pattern in the floristic composition of homegardens of the peninsula and have been called "structural species" (Jimenez-Osornio et al., 2005;Montañez-Escalante et al., 2014). Three of the structural tree species are Brosimum alicastrum Sw. (Moraceae; known in Spanish and Maya as ramón and ox, respectively), Cordia dodecandra A. DC. (Cordiaceae; known in Spanish and Maya as ciricote, and K'oopté, respectively) and Spondias purpurea L. (Anacardiaceae; known in Spanish and Maya as ciruela and abal, respectively). Brosimum alicastrum is known to have been domesticated by the Mayas who used it as a substitute for corn when maize was scarce (Puleston, 1971;Peters, 1983). The species might be one of the main sources of food for the Mayas of the classic period (years 250-900 CE) (Puleston, 1971). The diversified use of C. dodecandra (Rico-Gray et al., 1991), its use in ancestral ceremonies such as Hanal Pixan (López, 2012) and as a building material for Mayan houses (Sánchez Suárez, 2006) suggests pre-Hispanic management and domestication of the species (Colunga-García Marín and Zizumbo-Villarreal, 2004). Evidence of the domestication of S. purpurea 2,000 years ago has been reported in Mexico (Miller and Schaal, 2006;Miller, 2008;Piperno and Smith, 2012), and the Yucatan Peninsula was an important center of domestication (Ruenes-Morales et al., 2010;Fortuny-Fernández et al., 2017). The continued management of those species in the forest and the homegarden populations since pre-Hispanic times allow the domestication and conservation of those native trees in the Peninsula of Yucatan (Gómez-Pompa, 1987;Colunga-García Marín and Zizumbo-Villarreal, 2004).
In this work, we evaluated the genetic diversity of the three species in the forest and the homegarden populations belonging to two different climate regions of the state of Yucatan, as a part of a project to evaluate the effects of domestication on the diversity of the multipurpose native trees of Yucatán. Previous analyses with ISSRs genetic markers and morphological traits of the leaf, flower, and fruit, suggest that B. alicastrum, C. dodecandra, and S. purpurea display genetic lineages with larger fruits are more frequent in the homegarden than in the forest populations, as a consequence of the ongoing domestication process (Ferrer et al., 2020). Our aim, in the present study, was to find out if the process of domestication causes a differential genetic structure among the species with contrasting reproductive systems. Specifically, we expected that the gene flow between forest and homegardens populations of each region would be maintained in B. alicastrum and C. dodecandra, while in S. purpurea gene flow would be interrupted by clonal propagation in homegardens, even if genetic diversity is similar between forests and homegarden populations. We used ISSRs genetic markers to evaluate genetic diversity, differentiation, and structure of the three species in the northeastern and southwestern regions of the state of Yucatan where forest and homegarden populations coexist.

Study Species
Brosimum alicastrum is an evergreen tree that measures up to 45 m in height and 150 cm in diameter at breast height (DBH). Its flowers are globular, arranged in dense axillary clusters; malee flowers are grouped in globose catkins, have no corolla, and are yellow; female flowers are arranged in oblong heads and are green (Pennington and Sarukhán, 2005). It is a dioecious tree, as most of the 15 species of the Neotropical genus are (Berg, 2001), but there are reports of changes in sexual expression upon reaching maturity, either as monoecious or gynodioecious patterns (Lopez-Mata, 1987). In Yucatan's forests, it blossoms from January to June and bears fruit between April and September, having no defined or marked periods (Vázquez-Yanes et al., 2001;Hernández-Cruz, 2018). Pollination is anemophilous and dispersal is ornitho-chiroptera (Vázquez-Yanes et al., 2001), although there are reports of pollination by bees (Murawski and Hamrick, 1991). In the Yucatan it is managed in homegardens to obtain forage, with one or two prunings of 60% of the canopy per year, affecting the phenophases and in some cases, preventing fruition (Mex-Mex, 2017;Hernández-Cruz, 2018). In the forest, populations present low densities and receive no management, but they are logged for their wood (Mex-Mex, 2017).
Cordia dodecandra is a deciduous hermaphrodite tree that reaches heights of up to 30 m and diameter at breast height (DBH) up to 70 cm. Its flowers are arranged in a cymosepaniculated inflorescence, the calix is a 1-1.3 cm tubularcampanulate structure, the orange corolla is funnel-shaped 3-5 cm × 2.5-3.2 cm, with 12-16 lobes; and the fruits are conic fleshy drupes 3-4 cm, with a sole white seed (Pennington and Sarukhán, 2005). The species is heterostylous with two floral morphos and self-sterility (Canché-Collí and Canto, 2014), a trait that has been proposed as basal for the genus (Opler et al., 1975) and recorded for species of this Neotropical genus (Gasparino and Barros, 2009). In Yucatan, it blooms from March to June and bears fruit from May to June in the forests, and in homegardens, it blooms almost the whole year and bears fruit from March to November (Grajales-Puc, 2018). The leaf fall depends on the percentage of edaphic humidity, and with constant irrigation, its leaves remain and produce few flowers or fruits (Campos-Bobadilla et al., 2016). Bees and hummingbirds are potential pollinators and seed dispersal is carried out by birds and mammals (Vázquez-Yanes et al., 2001;Canché-Collí and Canto, 2014). In Yucatan, forest populations have very low density, receive no management and their flowers and fruits are not used. However, adult trees are logged for their wood, very much appreciated in the forestry industry; in homegardens, the trees are cultivated as ornamentals for their flowers, and their fruits are used for food (Ferrer et al., 2020;Hurtado-Torres et al., 2020).
Spondias purpurea is a deciduous tree that grows to a height of 15 m and diameter at breast height (DBH) up to 80 cm. It has small flowers of different colors arranged in short bunches measuring 0.5-8 cm, long, pink to red pentamerous calyx about 1 mm long, with five red to purple petals 3 mm long, with an annular fleshy disk divided into five adjoining nectaries. The fruit is an oblong-cylindrical drupe measuring up to 3.5 cm × 1.5 cm in forests and up to 9 cm × 4 cm in homegardens, with different colors from yellow to purple; the endocarp is yellow, fibrous and has one or more seeds (Pennington and Sarukhán, 2005;Ruenes-Morales et al., 2010). The species is polygamodioecious, among the 11 Neotropical species of the genus that are hermaphrodite (Mitchell and Daly, 2015), and the absence of seeds, in cultivated varieties, is associated with parthenocarpy (Avitia- García, 1997). In Yucatan's forests, flowering and fruition, take place from April to March and from March to April, respectively; in homegardens, both occur from January to July (Ruenes-Morales et al., 2010;Grajales-Puc, 2018). Pollinators for these species have not been identified, but it is believed that bees and small-sized wasps could be potential pollinators (Miller, 2011) and fruit dispersers are small and medium-sized vertebrates (Mandujano et al., 1994). In the forest, populations can be scarce and have low densities, are not managed nor are their flowers and fruits harnessed, but they are an important source of food for wild fauna (González Iturbe-Ruenes, 2019). In homegardens, fruits are food, and bark and leaves are used in traditional Mayan medicine (Ruenes-Morales et al., 2010).

Study Sites
The state of Yucatan is located in the biogeographic province of the Peninsula of Yucatan (Espinosa-Organista et al., 2008) in a gradient ranging from 0 to 400 m a.s.l. with a predominant humid subtropical climate (A), and dry climate (B) at the northeastern coast, according to Köppen climate classification modified for Mexico (1973). Data were analyzed according to the distribution of the localities in two regions: southwestern region and northeastern region. The southwestern region presents the tropical climate variant Aw0 sub-humid with rains in summer, the driest of the sub-humid climates (P/T ≥ 43.2 mm/ • C). The northeastern region has the variant Aw1(x') subhumid tropical climate with rains in summer, intermediate between the subhumid climates (P/T < 43.2 mm/ • C and less than 55.3 mm/ • C) and 10.2% of winter rain, according to the Köppen classification modified by García (1973).

Sampling
Sampling was carried out in the fragments of low height deciduous and medium-height semi-deciduous and semievergreen forests of two municipalities of the state of Yucatan: Tzucacab at the southwestern region and Tizimín in the northeastern region (Figures 1-3). In each region, a locality was selected in which fragments of the forest had been identified for each of the species, and between five to fifteen homegardens with the species were also selected (Figures 1-3). For each population type (forest and homegarden) twenty-two adult individuals with a DBH larger than 15 cm were selected. Each of the adult trees was tagged, georeferenced, and foliar tissue and young leaves were collected for DNA extraction.

DNA Purification and Amplification
DNA was extracted with liquid nitrogen from 0.25 g of lyophilized tissue following the CTAB extraction protocol modified by Falcón and Valera (2007). Amplification, visualization, and recording of the ISSRs markers for the three species were performed following the protocols described by Ferrer et al. (2020). Each of the identified bands was considered a locus. From the total of 84 individuals selected for each species, DNA amplification was achieved in 74 B. alicastrum, 81 C. dodecandra, and 83 S. purpurea and the number of loci for each species was 36, 40, and 68, respectively ( Table 1).

Statistical Analysis
The data on the presence/absence of each locus was recorded for each individual and they were grouped by population type (forest and homegarden) for each region (northeastern and southwestern). Statistics on the genetic diversity, genetic differentiation, and genetic relationships were estimated for each combination: region/population type (northeast homegarden, northeast forest, southwest homegarden, and southwest forest).
Because of the dominant inheritance nature of the ISSRs markers, estimations of the allelic frequencies cannot be made directly, so they were calculated with a Bayesian approach, without a priori information on the allelic frequencies (Zhivotovsky, 1999) to estimate statistics on genetic diversity, percentage of polymorphic loci (PPL) and expected heterozygosity (He), assuming Hardy-Weinberg balance, with the aid of AFLP SURV v. 1.0. (Vekemans et al., 2002). The frequency down-weighted marker values (DW) were also estimated (Schönswetter and Tribsch, 2005), as indicators of the proportion of rare alleles in the population, using the equations described by Winkler et al. (2010). The values of PPL, He, and DW were standardized to 0, because of the differences in the size of the samples, with the formula (x i − x)/x, where x i is the average value for each combination (region/population type) and x is the mean of the 4 combinations. These values were represented in a map elaborated in ArcMap 10.0 (ESRI-United States), together with the centroids of the geographic location of the individuals for each combination of region/population type. The base layers used were from the municipalities, the Yucatan Peninsula and Mexico's, all of them obtained from INEGI (National Institute of Statistics, Geography, and Informatics).
Genetic differentiation was estimated with the global estimator θ II , the estimators between the regions, between types of populations-inside the regions-and inside the populations ( RT , ST , and PT , respectively), and the statistical paired F ST and ST between populations. The statistical θ II was estimated with a Bayesian approach that allows comparing four related models in the program Hickory v. 1.1 (Holsinger andLewis, (2001-2007)). The four models are: (1) complete (includes the estimates of genetic diversity, inbreeding coefficient, and genetic differentiation index); (2) f = 0 (includes only the estimates of genetic diversity and genetic differentiation index); (3) θ = 0 (includes only the estimates of the genetic differentiation index and inbreeding coefficient); and (4) f free (the estimation of genetic diversity and genetic differentiation index is done choosing several inbreeding coefficients as priors). The election of the model that fits better the data is based on the value of the deviance information criterion (DIC). The model that has the smallest value and that differs in at least six units from others is chosen, when two models have similar values, the decision is made considering the model that have the smallest fitting values to the data of allelic frequencies (Dbar) and for the number of parameters used in each model (pD). The values of RT , ST , and PT were estimated with Arlequin v 3.5 using AMOVA (Excoffier et al., 2005). The paired F ST values for the four combinations of the region and population type were estimated with a random permutation of individuals with AFLPSURV v. 1.0. (Vekemans et al., 2002). The paired ST for the populations were estimated with Arlequin v 3.5 using an AMOVA approach, the significance of the estimated values was tested comparing the value obtained against the ordered values of 999 permutations (Excoffier et al., 2005).
The genetic relations were inferred in a heuristic manner, analyzing the consensus Neighbor-joining dendrograms built with Phylip v. 3.698 (Felsenstein, 2019) for 1,000 matrices of Frontiers in Ecology and Evolution | www.frontiersin.org FIGURE 1 | Distribution of forest (square) and homegarden (triangle) individuals of Brosimum alicastrum in two regions of Yucatán, México. The graphic depicted the standardized mean values for PPL-horizontal lines, DW -solid color, and He-diamond pattern [ x i − x /x where x i is the mean value per region/population type and x the overall mean value]; negative and positive values indicate values below and above the overall mean. The dendrogram for the paired F ST between each combination of the region (northeast -NE-dark green, and southwest-SW-light green) and population type (forests -squares-and homegardens -diamonds-) suggests almost no differentiation between southwestern forest and homegarden populations with a 100% bootstrap support, and this group and the homegarden of northeastern population. the paired F ST values for the four combinations of the region and population type were estimated with a random permutation of individuals with AFLPSURV v. 1.0. (Vekemans et al., 2002). The final dendrograms were edited with MEGA X (Kumar et al., 2018). The genetic structure of the populations was analyzed with the Bayesian assignment of the individuals to K clusters, using as priors the geographic location of the individuals with the mixture model at group level (combinations of region/population type) and of individuals. The values of the posterior probabilities were obtained with 500 simulations of the allelic frequencies for 10 partitions considering three to 20 clusters (K). The results of these analyses were used to analyze the model of admixture, using the running values of the program, and the best partition was chosen based on the logarithm of the value the marginal verisimilitude using BAPS v. 6.0 (Corander et al., 2008). The frequency of the posterior probabilities of belonging to a genetic group was clustered for the combinations of the region and population type using distruct v. 1.1 (Rosenberg, 2004). Besides, the belonging of the individuals to K clusters was estimated using the AMOVA K-means approach described by Meirmans (2012), considering pseudo F to determine the clusters in B. alicastrum and C. dodecandra; while the Bayesian Information Criteria was used in S. purpurea.

RESULTS
Mean PPL values varied 34.5% for S. purpurea and 41.0% for B. alicastrum. Mean DW values varied between 4.045 for C. dodecandra and 8.540 for S. purpurea. Mean He values varied between 0.096 for C. dodecandra and 0.116 for B. alicastrum (Table 1). Maximum PPL values for B. alicastrum and C. dodecandra were higher than for S. purpurea. However, minimum and maximum DW values were around three times higher in S. purpurea compared with the species that propagate sexually, particularly in homegarden ( Table 1). Genotypic diversity estimated as the He had minimum values which were slightly smaller in B. alicastrum and C. dodecandra-species with sexual propagation-compared with those of S. purpurea.
The full and f = 0 models had the smallest values for the deviance information criterion (DIC) and deviance (Dbar, Dhat); the differences of six units between the DIC values for the two models suggest that both are equally plausible ( Table 2). The θ II values were different than 0 for the three species. The values for the full model and the model f = 0 suggest that the genetic structure is weaker in C. dodecandra than that of B. alicastrum and S. purpurea, with values between 12 and 15 times lower in the full model and between 14 and 18 times lower in model f = 0, in the respective comparisons with the two species ( Table 3).
The genetic differentiation between regions RT was significant for B. alicastrum and S. purpurea, but not for C. dodecandra, while genetic differentiation between the population type-within the regions-ST and within populations PT were significant for the three species ( Table 4). The variance explained that the level of region ( RT ) was ca. three times larger in B. alicastrum than in S. purpurea, while at the level of population type-within the region-( ST ) was ca. three times smaller in B. alicastrum than in C. dodecandra and S. purpurea ( Table 4). The values of PT explain between 85 and 88% of genetic variance, in the three species ( Table 4).
On the other hand, almost no differentiation was found between the southwestern forest and homegarden populations, and this group and the homegarden of the northeastern population with a 100% bootstrap support in B. alicastrum and a 79% bootstrap support in S. purpurea; while in C. dodecandra the smallest differentiation was found between the southwestern FIGURE 3 | Distribution of forest (square) and homegarden (triangle) individuals of Spondias purpurea in two regions of Yucatán, México. The graphic depicted the standardized mean values for PPL-horizontal lines, DW -solid color, and He-diamond pattern x i −x /x where x i is the mean value per region/population type and x the overall mean value]; negative and positive values indicate values below and above the overall mean. The dendrogram for the paired F ST between each combination of the region (northeast -NE-dark yellow, and southwest-SW-light yellow) and population type (forests -squares-and homegardens -diamonds-) suggests almost no differentiation between southwestern forest and homegarden populations with a 79% bootstrap support, and this group and the homegarden of northeastern population.
forest and the northeastern homegarden population and this group and the northeastern forest population with a 100% bootstrap support according to the dendrograms for the paired F ST values (Figures 1-3).
In C. dodecandra two genetic groups were found and in B. alicastrum and S. purpurea three, according to the Bayesian assignment and the K-means. Genetic groups were shared among forest and homegarden populations within each region for the three species, although their frequencies varied among population type-region combination (Figure 4). In B. alicastrum, the populations of the northeastern region present two of the three groups (three in the northeastern forest according to the K-Means) and in the southwestern region the three groups. C. dodecandra forest populations in both regions and northeastern homegarden populations present predominantly one genetic group, while in the population from southwestern homegarden the second group predominates. S. purpurea forest populations presented two of the three groups; in the northeastern population the predominant group was the least represented in the southwestern region; and the homegarden populations also presented two groups, those from the Northeast to the third group, along with the group that was more frequent in the forest of the same region and in the southwestern populations to the same groups as those of the forest of that same region (Figure 4).

DISCUSSION
Genetic diversity in B. alicastrum and C. dodecandra-species with sexual reproduction-was higher than in S. purpurea, species which is clonally propagated, when estimated as polymorphism, but the DW and He values were higher in S. purpurea. The genetic diversity of the populations varies with the region in B. alicastrum; with the population type in C. dodecandra, and with both the population type and region in S. purpurea. The species have a very low genetic differentiation in C. dodecandra, the heterostylous species, and low in the dioecious species B. alicastrum and S. purpurea, according to the θestimates, moderate; null genetic differentiation between regions for C. dodecandra and low for B. alicastrum and S. purpurea according to the RT estimates; and with very low genetic differentiation between forest and homegarden of the same region for B. alicastrum and low genetic differentiation between those of C. dodecandra and S. purpurea according to the ST estimates. Even though the lowest paired F ST and ST values were found between the northeastern forest and homegarden for the three species, the dendrograms based in the paired F ST , the analyses of Bayesian assignment, and the K-Means suggest that the least differentiated populations are southwestern forest and homegarden populations of B. alicastrum and S. purpurea, and the southwestern forest and northeastern homegarden of C. dodecandra.
Domestication of trees, when accompanied by a shift in the reproductive system, may affect the genetic diversity of the species, a fact that is prevented in trees that remain as predominantly outcrossing species under cultivation (Miller and Gross, 2011). This fact was of particular importance for B. alicastrum and C. dodecandra, which rely on sexual reproduction in the homegardens, and remains as dioecious and distylous in the homegardens. Polymorphism of these species was slightly higher, but the number of single alleles and heterozygosity was smaller than that of S. purpurea, species that is clonally propagated in homegardens. This result agrees with observations of species with sexual and asexual reproduction having a larger diversity of effective alleles and heterozygotes compared with those that only reproduce asexually (Hamrick et al., 1992). The above is related to the ability of the lines that propagated clonally to protect polymorphisms through the fixation of the heterozygotes. Besides, any novel mutation in a clonal line will be preserved, increasing the number of rare or exclusive alleles (Balloux et al., 2003;Glémin et al., 2019). The domestication of S. purpurea in homegardens may be considered as a process that increases the genetic diversity of the species in this agroecosystems because heterozygotes variants of the species are kept in them; and that, according to Zohary (2004), turn them into "frozen cultivars" because of the clonal propagation.
The cultivation of B. alicastrum, C. dodecandra and S. purpurea in the Mayan solar (homegarden), as part of the ongoing process of domestication, allows their genetic conservation; which is of particular interest due to local and regional threatening of tropical forests of the Yucatan Peninsula. The regions and the population types of the species affected genetic differentially the diversity of the three species. Standardized estimators of genetic diversity of populations can function as indicators of isolation of the populations (Schönswetter and Tribsch, 2005). For instance, PPL and DW values are expected to be higher than the mean in isolated populations because of the accumulation of mutations, while in populations established recently the values are expected to be lower than the mean. According to this interpretation, it can be inferred that the populations of B. alicastrum in the forest and homegarden within the same regions were not completely isolated; however, those of the northeastern region, irrespectively of belonging to the forest or homegarden population, were established more recently than those of the southwestern region. At a regional level, this observation agrees with the results of different phylogeographic studies of the species at continental and regional (Central America) levels (Poelchau and Hamrick, 2013;Lander and Monro, 2015), which suggest that the species migrated from South America, the center of the origin of the species. Cordia dodecandra populations were found to be more isolated in the homegarden and/or alternatively, the populations of the forest may have been established more recently. Considering that the adult individuals in homegardens have a larger diameter at breast height, and therefore are older than those found in the forest (Hurtado-Torres et al., 2020), the difference observed in genetic diversity would be reflecting the presence of older lineages in homegardens. The decrease of the genetic diversity of the forest populations of C. dodecandra is a (2) f = 0 (includes only the estimates of genetic diversity and genetic differentiation); (3) θ = 0 (includes only the estimates of genetic diversity and inbreeding indexes); and, (4) f free (includes estimates of genetic diversity and genetic differentiation based on variable inbreeding indexes values used as priors).

Brosimum alicastrum
Full model consequence of the bottlenecks caused by the logging of adult trees to sell their wood and of the high levels of deforestation in the Yucatan Peninsula. The populations of S. purpurea in the northeastern homegardens and southwestern forests are the ones that could have been more isolated, given the higher genetic diversity that they hold. This result was not expected because in homegardens from both the northeast and southwest, the species is propagated clonally (Ruenes-Morales et al., 2010), which reproductively isolates it from the forest populations. However, the higher genetic diversity of the population can be associated with the exclusive presence of a genetic group in northeastern homegardens (Figure 3). The above results highlight the importance of homegardens as in circa situm reservoirs of the three species, be it because they maintain similar or higher levels of genetic diversity than in the forests. In the case of B. alicastrum and C. dodecandra the maintenance of gene flow between forest and homegardens populations promotes panmixia, whereas for S. purpurea, homegardens safeguard unique alleles and genotypes, even when the gene flow is restricted. The gene flow, among the forest and domesticated populations, is guaranteed in trees because they are characterized by having long life cycles and recurrent reproductive events (Miller and Gross, 2011). As a consequence, a higher genetic diversity, low genetic differentiation between populations, and a weak or inexistent genetic structure is expected (Loveless, 1992;Hamrick and Got, 1996). Even though in tropical regions the low population density and the asynchrony in the reproductive phenology of trees can cause a genetic differentiation by disrupting the gene flow mediated by pollen, it is expected that in species with reproductive systems that foster cross-pollination, the genetic structure will be weaker at a population than at regional level (Dick et al., 2008). Null or very low genetic differentiation in C. dodecandra, the heterostylous species, suggests that the gene flow is higher in this species, even at the regional level, while low genetic differentiation in B. alicastrum and S. purpurea, the two dioecious species, suggests that the gene flow is maintained between forest and homegarden populations, but is limited within each region. Both the dioecious species and the heterostylous species have outcrossing rates of 1 because by definition none of them can self-fertilize, and, in this sense, cross-pollination not only allows a larger gene exchange among individuals, it also promotes higher gene flow between populations (Glémin et al., 2006(Glémin et al., , 2019. Hence, the differences that are observed between the species respond to other life-history traits that can promote the gene flow of the species such as flowering synchrony, pollen dispersal vectors, and seeds dispersal vectors (Hamrick et al., 1992;Duminil et al., 2007Duminil et al., , 2009Dick et al., 2008;Gamba and Muchhala, 2020). Our results suggest that pollen and seed dispersal among forest remnants and homegardens is not interrupted by the management of the species in the solar Maya, and that biotic interaction of the species largely contribute to the genetic conservation of these species under domestication. The management practices that accompany the domestication process may affect the reproductive phenology, and when asynchrony blooming occurs under cultivation, isolation of the wild and domesticated populations is expected (Zohary, 2004). In B. alicastrum and C. dodecandra, species with sexual reproduction, cross-pollination, recurrent reproduction events, pollen dispersal by wind, large bees and vertebrates, and the zoochory can promote a high gene flow and explain the low genetic differentiation observed between the types of populations and between the regions. Asynchronous reproductive cycles are expected in the trees living in the tropics, as an effect of climate heterogeneity, with an interruption of gene flow promoting structuring genetics (Bawa and Hadley, 1990;Martin et al., 2009). We know that management affects the reproductive phenology, be it as a direct effect of removing the foliage, flowers, and fruits of B. alicastrum to feed domestic animals (Mex-Mex, 2017), or indirectly by the irrigation that increases soils humidity in homegardens, prolonging the presence of leaves and decreasing the frequency of bloom and fruition in C. dodecandra (Campos-Bobadilla et al., 2016). However, these changes related to unconscious domestication of the species (sensu Zohary, 2004), do not cause an asynchrony of the phenophases that impedes the blooming overlap of forest and homegarden populations (Grajales-Puc, 2018;Hernández-Cruz, 2018) therefore it is expected that the gene flow among populations from the forest and homegarden will persist as long the blooming overlap remains. In a recent review (Gamba and Muchhala, 2020), the reproductive system, lifecycle, and distribution, along with the pollination vectors and seed dispersals were found as the main factors in explaining genetic differentiation between populations. Our results are consistent with those of Gamba and Muchhala (2020) who found that the species pollinated by wind (B. alicastrum in our study) have higher genetic differentiation than those pollinated by vertebrates (C. dodecandra in our study). In the case of B. alicastrum, anemophilous pollination can be important to preserve cross-pollination within the populations by promoting genetic recombination and gene exchange. However, fruit dispersal by bats and their ability to travel large distances may be more determinant in the gene flow between populations within the same region. Brosimum alicastrum forest and homegardens populations are far apart from one another (8 and 25.1 km in the Northeast and Southwest regions (estimation based on the centroids of all the individuals). Artibeus jamaicensis, the most common bat frugivorous species in the Peninsula, may travel between 0.6 and 8 km each day foraging for food (Ortega and Castro-Arellano, 2001), so dispersal between both populations can be a common event, especially between the northeastern populations where the lowest geographical distances and genetic differentiation indexes values were found. Potential C. dodecandra pollinators include  Apis mellifera, native bees, and the hummingbirds Amazilia rutila, Archilochus colubris, and Chlorostilbon canivetii (Canché-Collí and Canto, 2014). The two latter species are migratory (Arizmendi and Ornelas, 1990), and it is possible that the movement of these birds may promote extensive gene flow because Archilochus colubris is known to travel up to 32 km per day in spring (Courter et al., 2013), a time that coincides with the flowering of the species from March to June. As it happens in B. alicastrum, C. dodecandra fruits are consumed by bats (Vargas-Contreras et al., 2009), and may disperse their seeds, accounting for the gene flow maintained between populations within the same region. The distances between B. alicastrum and C. dodecandra populations of the northeast and southwest range from 188 to 222.9 km (estimates based on the centroids of all the individuals), therefore neither hummingbirds nor bats can overcome these distances. To explain the low genetic differentiation between populations types and regions of B. alicastrum and C. dodecandra, it is necessary to consider that besides being present in homegardens, the species have been used in urban reforestation plans, in ranches to shadow and feed the cattle, in land boundaries as living fences (Vázquez-Yanes et al., 2001;Campos-Bobadilla et al., 2016;Mex-Mex, 2017;Hurtado-Torres et al., 2020). The isolated trees in these landscapes, can be promoting gene flow between remnant forest fragments and homegardens, as has been observed in another neotropical species, Swietenia humilis (White et al., 2002) and in Cordia africana (Derero et al., 2011). As a consequence of the mobility of the pollinating vectors and seed dispersals of B. alicastrum and C. dodecandra, we found that the traditional agroecosystems, mainly the homegardens, play an important role in maintaining connectivity of the remnant populations of the forest within the climate regions, and may even promote connectivity at a regional biogeographic scale.
On the other hand, the higher genetic differentiation of the populations of S. purpurea between forest and homegarden and between regions may be explained, as a consequence of a unidirectional gene flow. The results of this work and those of a previous paper (Ferrer et al., 2020) indicate that the populations of the homegardens host genotypes that came from the forest. Previous works on the spatial genetic structure of S. purpurea suggest that the pollen dispersal allow cross-pollination at low densities, ranging from 27.8 and 854 m and that seeds dispersal takes place between 7.2 and 833 m, which guarantees that the seeds are scattered far away from the mother tree (Cristóbal-Pérez et al., 2020). However, these dispersal ranges would not allow the pollen and seed exchange between populations from the forest and homegarden in our study, which have distances between them of 10 and 23.22 km in the northeast and southwest regions (estimations based on the centroids of all individuals). It is important to highlight that, unlike B. alicastrum and C. dodecandra, S. purpurea gene flow (via pollen and seeds) can be interrupted. The directional selection upon the fruits of the species (Ferrer et al., 2020) promotes that only female or even parthenocarpic individuals are maintained in the homegarden, and these individuals are incapable of contributing to the gene pool because they do not produce pollen nor seeds. In this sense, a conscious selection of the species (sensu Zohary, 2004), together with clonal propagation leads to the interruption of the gene flow of the species. This genetic differentiation of the populations is also evident at a biogeographic level for the known distribution of the species (Miller and Schaal, 2006;Miller, 2008;Fortuny-Fernández et al., 2017).
Because of the genetic differentiation of the species, a genetic structuring was found in the three species. The grouping patterns of the dioecious species are similar, suggesting that the populations of the southwest are more similar than those of the northeast, which is corroborated by the presence of similar genetic groups, according to the Bayesian assignments and the K-means, in these regions and the presence of exclusive groups in the northeastern populations (Figures 1, 3). This genetic structure is congruent with the one obtained for the three species using STRUCTURE for the Bayesian assignment of the genetic K groups reported by Ferrer et al. (2020). The South-North dispersal patterns that have been reported for B. alicastrum and S. purpurea in phylogeographic studies at the intercontinental level (Miller and Schaal, 2006;Miller, 2008;Poelchau and Hamrick, 2013;Lander and Monro, 2015;Fortuny-Fernández et al., 2017) also match with what we observed at the regional level. The population genetic structure of C. dodecandra has a different pattern that is related to the predominance of a genetic group, according to the Bayesian assignments and the K-means, in all populations except in the southwestern forest where this group is less frequent (Figure 2). In this case, structuring of the forest and homegarden populations may respond to larger connectivity between populations managed in other agroforestry systems or urban areas, or even to a differential selection of the trees that are kept in homegardens at the northeast region as compared with those of the southwest region. The role that other agroforestry systems and urban areas have in the gene flow of the species requires further study.
Domestication, as an ongoing process, that begun in prehispanic times and continues nowadays, has preserved the genetic diversity of B. alicastrum, C. dodecandra, and S. purpurea. The conservation of the species may rely on the Mayan homegardens and other agroforestry systems, either because the connection of wild and domesticated populations is allowed in species such as the dioecious B. alicastrum and the heterostylous C. dodecandra; or because, unique genotypes of species such as the dioecious S. purpurea are maintained in forest and homegarden, even when the gene flow is interrupted by clonal propagation. Overall, with the genetic diversity data, these results suggest that the dynamics of colonization and extinction of the populations in the forest and homegardens in Yucatan are intimately linked to the domestication process of the species because Maya Yucatec manages this species as forest resources.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
MF, HE-M, PM-E, MR-M, and JJ-O conceived and designed the work and collected samples from the field. MF and CT-R collected and analyzed the data in the laboratory. MF conducted statistical analyses and data interpretation. All authors drafted the article and gave their final approval of the version to be published.

FUNDING
CONACyT provided a grant to the project CB-2014 236428 and the scholarship 24904 to CATG.