Degradation of grasslands results in significant loss of the products and services that these ecosystems provide to humans, such as livestock fodder, water capture, soil retention and carbon sequestration¹. Some of the main problems affecting grasslands are fragmentation, expansion of invasive species, conversion to agricultural uses, and human population growth, among others². Grassland rehabilitation is becoming more frequent, and in recent years research has focused on selecting the grass species genotypes most apt for this application^3,4,5. However, genetic material selection has concentrated mainly on agronomic traits, with minimal importance given genetic structure and potential adaptation to climatic conditions.

A grass species’ genetic structure is largely determined by its adaptation to environmental conditions within its distribution area. For this reason, understanding population genetic structure can help delimit genetic type distributions and their potential use in ecological restoration programs^6,7. In northern Mexico, one of the species most used in grasslands restoration is sideoats grama [Bouteloua curtipendula (Michx.) Torr.] because it adapts to a wide range of climates and has excellent forage value. Due to its potential, recent research has emphasized selection of outstanding B. curtipendula genotypes for use in restoring grasslands in the state of Chihuahua, Mexico^8,9,10. Identifying ideal genotypes requires analysis of the genetic structure of B. curtipendula populations in the state and of their environmental aptitude.

Mathematical models, such as the maximum entropy method applied in the MaxEnt software¹¹, are an effective way of estimating areas with potential environmental suitability for the distribution of a species. Applying these models via MaxEnt provides advantages such as that only presence data is required, reliable results are produced with limited data, continuous and/or categorical environmental data are used, and are these graphically displayed via maps, facilitating interpretation. MaxEnt also employs a Jackknife test to identify the relevance of each environmental variable in species distribution^11,12. Models generated with this program have been used widely to estimate climatic niche for flora and fauna species^13,14. It has been used to identify areas with environmental aptitude with potential B. curtipendula distribution in Mexico and the United States^15,16, a prerequisite for ecosystem rehabilitation¹⁷. However, MaxEnt has only been applied at the species level. It has not been used to estimate areas with environmental suitability for genetic types within B. curtipendula, which could make ecosystem restoration programs more efficient. The present study objective was to evaluate B. curtipendula genetic structure in the state of Chihuahua, Mexico, by applying MaxEnt models at the genetic group level and identifying those areas with environmental suitability for the resulting genetic groups.

Sampling was done in 51 populations, distributed in 29 municipalities of the state of Chihuahua, in northern Mexico (32° and 25° N; -103° and -109° W). To include as much genetic diversity as possible, sampling sites were located in the state’s arid and semi-arid ecoregions (Figure 1).

Figure 1 Geographic location of collected genetic material, corresponding to 51 sideoats grama (Bouteloua curtipendula) populations in Chihuahua, Mexico 

The genetic structure of B. curtipendula was analyzed using Amplified Fragment Length Polymorphisms (AFLP) molecular markers. Plant leaves were collected at each site and DNA extracted from them based on the method proposed by Doyle & Doyle¹⁸. The AFLP analysis was done using the method proposed by Vos et al.¹⁹. First, 2 µl diluted DNA were digested by the EcoRI and MseI restriction enzymes and the digested DNA fragments ligated with EcoRI and MseI adapters. For pre-amplification, an extra nucleotide was added to the primers (EcoRI + A and MseI + A). Selective amplification was done using four combinations of fluorescence-labeled primers: MseI + CTG - EcoRI + AAG, MseI + CTG - EcoRI + ACT, MseI + CAG - EcoRI + AGG, MseI + CAG - EcoRI + AAC. The polymerase chain reaction (PCR) was done in a thermal cycler (Applied Biosystems Veriti 2720), with the following program: one cycle of 94 °C for 30 sec, 65 °C for 30 sec, 72 °C for 1 min; 12 cycles of 94 °C for 30 sec, 65 °C for 30 sec, 72 °C for 1 min; and 23 cycles of 94 °C for 30 sec, 56 °C for 30 sec, 72 °C for 1 min. The selective amplification products (2 µL) were mixed with 8 µl formamide and 1 µL Eco 700 GeneScan label (Applied Biosystems). Separation of the amplified fragments was done in a LI-COR DNA analyzer, loading 0.8 µL sample per well. Fluorescence-marked oligos or primers at different wavelengths (700 nm and 800 nm) were used. With the band pattern produced by the AFLP analysis, a binary band presence (1) / absence (0) matrix was constructed.

A genetic structure analysis, based on the Bayesian clustering algorithm, was applied to the binary data using the STRUCTURE version 2.3.4 software^20,21. The STRUCTURE program was run 30 times for each K number of genetic clusters and analyzed from K=1 to K=10. Ten thousand (10,000) Markov-Monte Carlo chain (CMMC) repetitions and 100,000 burn-in periods were done in each run. This analysis was performed using a correlated allele frequency and admixture model. The optimal number of K clusters was considered to be that which attained the highest value for the average posterior probability (LK) and ΔK, according to the criteria proposed by Evanno et al²². The LK and ΔK values were obtained from the Structure Harvester website²³. An analysis of the association between climatic zones and the distribution of the B. curtipendula populations in the genetic groups was done with a χ² test of independence (α=0.05).

An analysis of molecular variance (AMOVA) was applied to compare the groups formed in the genetic structure analysis²⁴, using the GenAIEx ver. 6 software²⁵. Using the F (Φ_ST) statistics produced with the AMOVA, the inter-group gene flow index was calculated with the formula N_m= [0.25 (1- Φ_ST)/(Φ_ST)]²⁶. The genetic data were also analyzed using the Monmonier algorithm to detect possible eco-geographical barriers affecting interpopulational gene flow. This analysis was run with the Barrier ver. 2.2 software²⁷, in which the Bootstrap values of each barrier were calculated with 100 Dice coefficient distance matrices.

For the genetic groups formed by the STRUCTURE analysis, the GenAIEx ver. 6 software²⁵ was used to calculate diversity statistics, polymorphic loci percentages, the average number of alleles per locus, the number of effective alleles, and the Shannon information (I) and Nei diversity (He) indices. Both these indices were estimated based on the assumption that each locus represents a pair of alleles when the presence or absence of an AFLP fragment is identified in a band. Diversity statistics for each population were compared using the Wilcoxon test with the Bonferroni correction (α=0.05).

Genetic group environmental fitness was identified using the MaxEnt algorithm in the MaxEnt ver. 3.3.3 software¹¹. For each genetic group, the model was run separately with the coordinates of the genetically analyzed populations. Of the total data set, 75 % was used to test the models and the remaining 25 % to validate the models using the bootstrap test with 50 replicates. The generated environmental suitability models were evaluated using receiver operating characteristic (ROC) curve and area under the curve (AUC) analysis. The AUC score is useful for measuring model performance; the higher (closer to 1) the AUC value the better the model estimates species presence probability.

A total of 22 climatic variables were used as predictors to model the potential distribution of B. curtipendula. Nineteen were bioclimatic variables: annual mean temperature (Bio1); mean diurnal range (Bio2); isothermality (Bio3); temperature seasonality (Bio4); max temperature of warmest month (Bio5); min temperature of coldest month (Bio6); temperature annual range (Bio7); mean temperature of wettest quarter (Bio8); mean temperature of driest quarter (Bio9); mean temperature of warmest quarter (Bio10); mean temperature of coldest quarter (Bio11); annual precipitation (Bio12); precipitation of wettest month (Bio13); precipitation of driest month (Bio14); precipitation seasonality (Bio15); precipitation of wettest quarter (Bio16); precipitation of driest quarter (Bio17); precipitation of warmest quarter (Bio18); and precipitation of coldest quarter (Bio19)²⁸. Three additional variables were included: average annual solar radiation (Rad); average annual evapotranspiration (Vapr); and average annual wind speed (Wind). These variables were obtained from the WorldClim database (https://www.worldclim.org) and limited to the geographic space of the state of Chihuahua using the ArcMap ver. 10.3 software. Climate data are interpolation estimates for 1950-2000, with a 30 arc-seconds spatial resolution. Applying MaxEnt generated a logistic map showing potential B. curtipendula distribution within a 0 (inadequate) to 1 (optimal) value range.

The genetic structure analysis (STRUCTURE) divided the sampled B. curtipendula populations into two genetic groups since K=2 produced the highest ΔK and LK values (Figure 2).

Values range from K=1 to K=10; K=2 is the optimum number of groups.

Figure 2 Delta K (ΔK) and average posterior probability (LK) values for the genetic structure of 51 sideoats grama (Bouteloua curtipendula) populations in Chihuahua, Mexico (grupos = groups). 

Group genetic divergence appears to have been generated by adaptations to Chihuahua’s climatic conditions since there was a high association between the state’s ecoregions and group formation (χ² = 32.9; P<0.0001). Most (89 %) of the Group 1 populations are in the semi-arid region and 11 % are in the arid region. In contrast, 90 % of the Group 2 populations are in the arid region and only 10 % in the semi-arid region (Figure 3B). These results agree with those of the BARRIER analysis, which identified genetic barriers coinciding with the border between the arid and semi-arid regions (Figure 3C).

A) STRUCTURE analysis using K= 2 run and 186 AFLP fragments; the colors represent the proportion of probability of belonging to each genetic group. B) Group structure in geographic context of state’s climatic zones. C) Genetic barriers identified by BARRIER analysis; yellow lines represent barriers and numbers are Bootstrap values (1000 bootstraps). Pie charts represent percentage of populations in each genetic group within each region.

Figure 3 Genetic structure of 51 sideoats grama (Bouteloua curtipendula) populations in Chihuahua 

The two groups generated by the STRUCTURE analysis were different (P<0.0001) according to the AMOVA. However, these differences only explained 7 % of overall variation, indicating the presence of substantial intragroup genetic variability. In grasslands it is common for intergroup differences to explain small proportions of overall variation due to generally high interpopulation genetic flow^29,30. The low interpopulation genetic differentiation (7 %) identified by the AMOVA apparently contrasts with the marked association between the state’s climatic zones and group formation. This can occur because formation of a single locus can be closely linked to adaptation to climatic conditions. For example, a study of Bouteloua gracilis populations identified a close relationship between the frequency of a single locus and precipitation in the driest quarter of the year (R²=0.84) and precipitation seasonality (R²= 0.77)⁷.

Compared to previous research using AFLP markers, the observed intergroup genetic exchange (F _st = 3.41) can be considered high. A study of Stipa pulcherrima grass distributed throughout Europe and Asia found a 0.76 gene flow among 30 populations³¹, while one of Microlaena stipoides found an F _st of 0.02 among 85 Australian populations²⁹.

The STRUCTURE analysis groups exhibited differences (P<0.05) in the evaluated diversity parameters. Group 2, mainly distributed in the arid region, had the highest values (P<0.05) in all the diversity parameters (Table 1). This coincides with the higher genetic diversity levels observed in Festuca ovina in arid zones³². In this study a positive correlation was identified between the Nei genetic diversity index and mean annual temperature (r= 0.56), while a negative correlation (r= -0.60) was found between this index and mean annual precipitation. Similar results have been reported for Dactylis glomerata³³, and many studies indicate that plant populations tend to have greater diversity in adverse environments^34,35. Plant populations in extremely arid environments tend to develop greater genetic diversity as a mechanism of adaptation to drought³⁶.

Based on the Shannon information index (I), both STRUCTURE groups presented high intragroup genetic diversity (Group 1= 0.302, Group 2= 0.427). By comparison, an evaluation of 56 accessions of Panicum virgatum produced an I value of 0.317³⁷, and one of 281 cultivars of Pennisetum purpureum produced I values ranging from 0.12 to 0.34³⁸. Both these studies included large numbers of genotypes, further emphasizing that the present genetic diversity results for B. curtipendula are relatively high in both groups.

The MaxEnt model was run using the coordinates of 20 populations for Group 1 and 32 for Group 2. Average AUC value for the Group 1 climatic niche was 0.91 (SD = 0.031) while that for Group 2 was 0.93 (SD = 0.015). Both values indicate the estimated environmental suitability of both genetic groups is highly reliable^11,39. The respective climatic niches of the groups clearly differ (Figure 4). Group 1 is distributed mainly in the center and south of the state (Figure 4A), while Group 2 is distributed in the center, west and northeast (Figure 4B). This suggests that the genetic groups diverged by evolutionary adaptation and that each genotype is adapted to regional climatic conditions.

Group 1 (A) and Group 2 (B). Red color represents areas of greatest environmental aptitude and blue those of least aptitude.

Figure 4 MaxEnt model maps for two sideoats grama (Bouteloua curtipendula) genetic groups identified by AFLP markers and STRUCTURE analysis  

Restoration programs involving B. curtipendula would probably benefit most by using local genotypes from each ecoregion. However, grassland restoration programs carried out to date in Mexico have employed genetic material from outside the region, mainly due to low seed availability. Bouteloua curtipendula varieties from the United States have been used; for example, the El Reno variety had unfavorable performance due to its low adaptability to the Mexican climate. A comparison of the El Reno variety with 277 B. curtipendula genotypes from different states in Mexico found that more than half of the genotypes provided better productive potential than the El Reno variety⁴⁰. Other studies⁴ have reported that, compared to native material, the El Reno variety has low establishment and forage production capacities. Local materials are most effective in grassland restoration programs since they guarantee a greater probability of success and preserve local genetic structure^41,42. This agrees with the present results, emphasizing that B. curtipendula revegetation programs in Chihuahua should use genotypes specific to each of the state’s ecoregions.

The high diversity identified within each genetic group suggests the possibility of selecting outstanding genotypes for each ecoregion. Previous studies have addressed genotype selection^9,43, but focused mainly on productive characteristics and gave little weight to the environmental adaptation of each genotype. In contrast, the present results provide basic information on the potential environmental suitability of B. curtipendula populations in Chihuahua, which could be valuable in selection programs for productive genotypes.

The variables that contributed most to the Group 1 environmental suitability model were precipitation seasonality (79.5 %) and precipitation in the coldest quarter (7.4 %). The Group 2 environmental suitability model was mainly influenced by precipitation in the coldest quarter (36.2 %), annual thermal oscillation (13.2 %), average annual solar radiation (8.1 %) and average temperature of the driest quarter (7.4 %; Table 2). These results agree with previous studies modeling the B. curtipendula climatic niche in Mexico and the United States. In one, the ecological descriptors of thermal oscillation and precipitation contributed most to potential B. curtipendula distribution in Mexico¹⁶, while another found average annual temperature contributed most to potential distribution in the United States¹⁵.

A detailed analysis using response curves of the most influential variables showed Group 1 genotypes to have a higher probability of developing in areas with a 120+ precipitation seasonality coefficient and that they develop best in areas with 0 to 20 mm precipitation in the coldest quarter (December, January and February) (Figure 5). Group 2 genotypes also develop best in areas with 0 to 20 mm precipitation in the coldest quarter. In addition, they prefer areas with a 49 to 53 annual thermal oscillation, 18,500 to 19,000 w m^-2 annual average solar radiation and a 6 to 10 °C average temperature of driest quarter. Group 1 genotypes apparently do not resist long periods of drought and needs occasional precipitation during the warmer months (June-October). Group 2 genotypes can resist less thermal oscillation, higher temperatures and greater solar radiation, but need precipitation concentrated in the warmer months. These response curve results highlight the evolutionary differences between the two genetic groups and the importance of using them in their source ecoregion.

Figure 5 Response curves for two genetic groups of sideoats grama (Bouteloua curtipendula) based on the variables with greatest influence on environmental aptitude.  

In conclusion, the state of Chihuahua, Mexico, sideoats grama Bouteloua curtipendula populations are divided into two genetic groups. Their distributions are highly influenced by ecoregion climatic conditions such that each genetic group has a different climatic niche. Restoration programs involving Bouteloua curtipendula could benefit from using local genotypes from specific ecoregions in environmentally suitable areas. These genotypes could also be used in edaphoclimatic conditions similar to those of their point of origin. The genetic diversity identified within each gene pool provides an opportunity for developing outstanding genotypes for use in ecoregion-specific grassland restoration programs. However, climatological projections are still needed to consider how climate change may affect Bouteloua curtipendula genetic types and what this could mean for future grassland restoration projects.