Population genetic analysis and bioclimatic modeling in Agave striata in the Chihuahuan Desert indicate higher genetic variation and lower differentiation in drier and more variable environments
Abstract
PREMISE OF THE STUDY:
Is there an association between bioclimatic variables and genetic variation within species? This question can be approached by a detailed analysis of population genetics parameters along environmental gradients in recently originated species (so genetic drift does not further obscure the patterns). The genus Agave, with more than 200 recent species encompassing a diversity of morphologies and distributional patterns, is an adequate system for such analyses. We studied Agave striata, a widely distributed species from the Chihuahuan Desert, with a distinctive iteroparous reproductive ecology and two recognized subspecies with clear morphological differences. We used population genetic analyses along with bioclimatic studies to understand the effect of environment on the genetic variation and differentiation of this species.
METHODS:
We analyzed six populations of the subspecies A. striata subsp. striata, with a southern distribution, and six populations of A. striata subsp. falcata, with a northern distribution, using 48 ISSR loci and a total of 541 individuals (averaging 45 individuals per population). We assessed correlations between population genetics parameters (the levels of genetic variation and differentiation) and the bioclimatic variables of each population. We modeled each subspecies distribution and used linear correlations and multifactorial analysis of variance.
KEY RESULTS:
Genetic variation (measured as expected heterozygosity) increased at higher latitudes. Higher levels of genetic variation in populations were associated with a higher variation in environmental temperature and lower precipitation. Stronger population differentiation was associated with wetter and more variable precipitation in the southern distribution of the species. The two subspecies have genetic differences, which coincide with their climatic differences and potential distributions.
CONCLUSIONS:
Differences in genetic variation among populations and the genetic differentiation between A. striata subsp. striata and A. striata subsp. falcata is correlated with differences in environmental climatic variables along their distribution. We found two distinct gene pools that suggest active differentiation and perhaps incipient speciation. The detected association between genetic variation and environment variables indicates that climatic variables are playing an important role in the differentiation of A. striata.
Climatic conditions are considered among the most important determinants of different aspects of plant ecology because climate affects plant reproductive ecology (including pollination, inbreeding depression, and seed dispersal), population dynamics (including individual densities and sizes, seed production, germination, and seedling establishment), interactions such as herbivory and predation, local adaptations, and geographic distribution (Foll and Gaggiotti, 2006; Kiambi et al., 2008; Temunović et al., 2012; Hamasha et al., 2013). On the other hand, it has been argued that genetic variation and differentiation of plant populations is modeled—at least in part—by the life histories and ecological traits of the plants and by different evolutionary and historical events; all of these factors are related to climatic conditions (Hamrick and Godt, 1996; Hewitt, 1996).
A relevant question is whether there is an association between environmental variables and genetic differentiation and the possible mechanisms of this association (Temunović et al., 2012; Hamasha et al., 2013; Lira-Noriega and Manthey, 2014). For instance, we would expect that the climatic conditions where the species originally evolved represents the optimal niche of the species and that in these original localities, populations are better adapted and larger. Also, these populations are likely to have been present for a longer period, without going through bottlenecks, thus storing higher levels of genetic variation. Populations in contrasting climatic conditions, however, may be marginal or of recent origin and thus display larger levels of genetic differentiation and less genetic variation (Lira-Noriega and Manthey, 2014). In this way, species living in the same general area could have contrasting patterns of genetic variation and differentiation, depending on where they originally evolved and their later phylogeographic history. Foll and Gaggiotti (2006) detected significant correlations between genetic differentiation (measured either as FST or GST) and average annual temperature in some plant species. Kiambi et al. (2008) found that genetic diversity was a function of annual rainfall in wild rice Oryza longistaminata, while Hamasha et al. (2013) reported higher genetic variation in drier environments in grasses of the genus Stipa.
To study whether there are associations between environmental parameters and genetic diversity and structure and to explore the causes of any associations, we need detailed analyses of genetic variation and differentiation along environmental gradients in common species, mainly because in these species the populations will be large, thus minimizing the potentially confusing factor of differentiation due to genetic drift (and thus not be due to the environmental factors). It will be useful if these species have widespread distributions and thus live in contrasting climatic conditions and are of recent origin; so any patterns are not confused by further genetic differentiation because of genetic drift and by the further accumulation of mutations (Díez et al., 2013; Hamasha et al., 2013; Foll and Gaggiotti, 2006).
The genus Agave L. seems to be an adequate system to explore the associations between environmental parameters and genetic diversity and structure. Agave includes more than 200 species, all endemic to the Americas. Agave is a typical element of arid environments in Mexico (Good-Avila et al., 2006; García-Mendoza, 2011; Eguiarte et al., 2013), and it is of a relatively recent origin, 7.8−10.1 million years ago (Ma) (Good-Avila et al., 2006; Flores-Abreu, 2007). Mexico is its center of diversity, with ca. 159 species, of which 119 are endemic (Garcia-Mendoza, 2011), and several species have a widespread distribution (Gentry, 1982; Eguiarte et al., 2013). The geological, topographical, and climatic conditions are among the main features of the Mexican territory that have been suggested to lead to the diversification of plant species (Rzedowski, 1978; Gentry, 1982; Toledo and Ordonez, 1993). However, in the case of Agave, the factors that have led to its diversification are still debated (Álvarez de Zayas, 1989; Good-Avila et al., 2006; Rocha et al., 2006; Eguiarte et al., 2013). For instance, Silva-Montellano and Eguiarte (2003) found that in A. lechuguilla genetic variation increases while its genetic differentiation decreases in the southern populations (that in general terms are less arid) of the Chihuahuan Desert.
Inter-simple sequence repeats (ISSRs) are dominant PCR-related molecular markers that have been used to study genetic variation and structure in several plant species and have proved useful for measuring genetic diversity in a variety of annual and perennial plants (Tanya et al., 2011; Saini et al., 2013). Using such molecular markers, Wolff and Morgan-Richards (1998) distinguished subspecies in Plantago major. Also, ISSRs have been used to make phylogenetic inferences (Archibald et al., 2006) and in phylogeographic and population genetics analyses (Hess et al., 2000; Eguiarte et al., 2013). They have also been used in several studies of the Agave genus (Rocha, 2006; Vargas-Ponce et al., 2009; Aguirre-Dugua and Eguiarte, 2013; Eguiarte et al., 2013), thus allowing comparisons with a relatively large data set within the genus.
Here we analyzed the genetic variation and structure using ISSRs along the geographical distribution of two subspecies of Agave striata Zucc. in the Chihuahuan Desert and adjacent areas in Mexico. Agave striata is a polycarpic (iteroparous) species with generalist pollination, sometimes involving nectar-feeding bats, and it has a wide geographic distribution through the arid and semiarid environments in Mexico and USA (20°–26° latitude) (Gentry, 1982; García-Mendoza, 2011; Rocha et al., 2005, 2006). We evaluated the environmental similarities and differences within the species to explore the association between bioclimatic variables and genetic variation and differentiation through species distribution modeling and multifactorial analysis of variance.
MATERIALS AND METHODS
Study group
We analyzed genetic and environmental characteristics from populations of the two subspecies of A. striata: A. striata subsp. falcata and A. striata subsp. striata (Fig. 1A and 1B; see Appendix S1 in Supplemental Data with the online version of this article). Agave striata plants are compact, short-stemmed, with a many-leaved rosette that forms a large dense cluster ranging from pale green to red or purplish. The leaves are 0.5−1 cm wide, linear, thick, straight to arching, and convex above. The inflorescences are apicate, with clusters of 50 to 200 flowers. Flowers are long and tubular, with short tepals, petals greenish yellow or red to purple, and anthers bronze to brownish (Gentry, 1982). Agave striata subsp. falcata has fewer leaves, 0.3 to 0.8 cm. thicker than those of A. striata subsp. striata, and its anthers are slate colored to yellow (Gentry, 1982). Agave striata subsp. striata has a wide southern distribution, from the Mexican states of Coahuila to Hidalgo (Fig. 1C), growing in dry to semidry environments (i.e., xerophytic scrub, pine–oak forest, deciduous forest) (Gentry, 1982). In contrast, A. striata subsp. falcata has a more restricted and northern distribution, mainly in the northern Chihuahuan Desert (Fig. 1C), and grows in dry environments (i.e., xerophytic scrub, pine forest).
(A) Agave striata subsp. falcata, from population F9, Coahuila. (B) Agave striata subsp. striata, from population F15, San Luis Potosí. (C) Distribution of Agave striata subsp. striata in green circles (populations F3, F4, F5, F14, F15, and F16) and of subsp. falcata in yellow circles (populations F8, F9, F10, F11, F12, and F13) and the locations of sampling sites used for genetic analysis (circles with black points). The dark lines represent the approximate distribution of the main Sierras in the area: Trans Mexican volcanic belt (1) and the Sierra Madre Oriental (2).
Population sampling
We sampled six populations of each subspecies throughout the range of A. striata, in arid and semiarid regions of central Mexico and the Chihuahuan Desert, over a range of ca. 740 km north to south (Fig. 1C, Table 1). Localities were selected from Gentry (1982). We collected the youngest leaf of each rosette and froze each sample in liquid nitrogen, from an average of 45 individuals per population, reaching a total of 541 individuals. To avoid sampling genetic clones within populations, we sampled plants located at a minimum of 3 m from one another.
| Code/Population | State | Altitude (m a.s.l.) | Latitude N | Longitude W | n | HS (±SD) | %P |
|---|---|---|---|---|---|---|---|
| Agave striata subsp. striata | |||||||
| F3/ Venados | Hidalgo | 1400 | 20°28′00″ | 98°40′45″ | 40 | 0.2375 ± 0.009 | 53.19 |
| F4/ Bucareli | Querétaro | 1600 | 21°05′06″ | 99°36′24″ | 28 | 0.2128 ± 0.014 | 48.93 |
| F5/ Zimapán | Hidalgo | 1780 | 20°43′06″ | 99°39′00″ | 43 | 0.1847 ± 0.008 | 42.55 |
| F14/ Dr. Arroyo | Nuevo León | 1897 | 23°42′20″ | 100°08′32″ | 36 | 0.1943 ± 0.009 | 48.93 |
| F15/ Guadalcazar | San Luis Potosí | 1527 | 22°38′14″ | 100°30′22″ | 50 | 0.2542 ± 0.008 | 55.31 |
| F16/ Matehuala-México | San Luis Potosí | 1899 | 22°17′46″ | 100°48′44″″ | 60 | 0.1949 ± 0.007 | 44.68 |
| Subspecies total | 231 | 0.2131 ± 0.005 | 70.27 | ||||
| Agave striata subsp. falcata | |||||||
| F8/ Saltillo-Monclova | Coahuila | 1126 | 25°41′19″ | 100°59′39″ | 58 | 0.3093 ± 0.009 | 72.34 |
| Km 17–18 | |||||||
| F9/ Saltillo-Monclova | Coahuila | 1296 | 26°19′09″ | 101°21′10″ | 61 | 0.3674 ± 0.012 | 74.46 |
| Km 120–121 | |||||||
| F10/ Ramos Arizpe | Coahuila | 1267 | 25°37′16″ | 100°49′32″ | 50 | 0.3602 ± 0.010 | 74.46 |
| F11/ Saltillo-Matehuala | Coahuila | 1873 | 25°23′20″ | 100°47′47″ | 39 | 0.3635 ± 0.012 | 65.95 |
| Km 24 | |||||||
| F12/ Saltillo-Matehuala | Nuevo León | 1905 | 24°32′18″ | 100°16′44″ | 48 | 0.3457 ± 0.011 | 70.21 |
| Km 112 | |||||||
| F13/28 km al N de Matehuala | Nuevo León | 1868 | 23°50′20″ | 100°30′52″ | 54 | 0.3201 ± 0.010 | 59.57 |
| Subspecies total | 310 | 0.3444 ± 0.007 | 80.85 | ||||
| Agave striata (species total) | 541 | 0.2448 ± 0.005 | 80.85 |
DNA extraction and ISSR protocol
We extracted DNA from leaf samples using a modified CTAB technique (Doyle and Doyle, 1987).
For amplification, we used three different ISSRs primers with the repeated dinucleotides AC, AG, and CT (Rocha, 2006). Primers 811 (5′-GAGAGAGAGAGAGAGAC-3′), 846 (5′-CACACACACACACACART-3′; R = A or G), and 853 (5′-TCTCTCTCTCTCTCTCRT-3′; R = C or T) produced clear and variable banding patterns. PCR reactions were performed in a total volume of 30 µL with final concentrations of 0.2 mM of each of the four dNTPs (Thermo Scientific, Vilnius, Lithuania), 0.4 µM of each primer, 1 U of Taq polymerase (Roche, Basel, Switzerland), and 30–40 ng of template DNA. The amount of MgCl2 varied depending on the primer: 2.4 mM for primer 811, 1.83 mM for primer 846, and 2.7 mM for primer 853. The amplification protocol for primer 811 was 94°C/3 min; followed by 38 cycles of 94°C/40 s, 56°C/20 s, and 70°C/50 s + 4 s each/cycle; and one final step of 70°C/5 min. For primers 846 and 853, the amplification protocol was 95°C/4 min; followed by 37 cycles of 94°C/30 s, 55°C/40 s for primer 846 or 55.5°C/30 s for primer 853, and 72°C/1 min + 4 s each/cycle; and one final step of 70°C/5 min. PCR products were electrophoresed on 1.5% agarose gels and stained with ethidium bromide. We used a DNA ladder of 100 bp (Invitrogen, Carlsbad, California, USA) to visualize the size of the fragments.
Genetic diversity in the subspecies of Agave striata
We generated a matrix of presence (1) or absence (0) of bands in each individual for all the loci identified. With these data, we obtained Bayesian estimates of the expected heterozygosity per population (HS) and standard deviation (±SD) with the software Hickory v. 1.1 (Holsinger et al., 2002). Hickory is a Bayesian-based program that allows direct estimates of structure without assuming previous knowledge of the degree of inbreeding or if the genotypes are in Hardy–Weinberg proportions (Holsinger et al., 2002). The Hickory analyses were run under the f-free model. For each subspecies, we used 10 separate runs with 5000 iterations as burn-in, followed by 25,000 iterations, sampling values every 5th iteration, starting from flat-priors (alpha = 1). The convergence of the different runs was estimated with a Gelman and Rubin (1992) test implemented in Coda V. R package (Plummer et al., 2006). We also calculated the proportion of polymorphic loci with a 95% criterion per population (%P) in TFPGA v. 1.3 (Miller, 1997).
We estimated genetic structure with the F-statistics and their standard deviation following Weir and Cockerham (1984). We also used the Fisher's combined probability (Raymond and Rousset, 1995) to identify loci with significant differences in their allelic frequencies among populations with the software TFPGA v. 1.3 (Miller, 1997), running a Markov chain with 1000 steps of burn-in and a total number of 10,000 steps. The distribution of genetic variation within and among populations of the subspecies was compared with analysis of molecular variance in the program Arlequin v. 1.1 (AMOVA; Excoffier et al., 1992).
The genetic structure of A. striata was analyzed with the program structure v. 2.2 (Pritchard et al., 2000, Falush et al., 2003, 2007). We used the no admixture model implemented in the program, as suggested for dominant markers in the literature (Falush et al., 2007). We ran 50,000 Markov chain Monte Carlo (MCMC) runs with 25,000 burn-in periods for each run, and 30 runs to estimate 1 to 13 populations. We detected the optimal clustering number (K) using the ad hoc statistic ΔK, based on the rate of change in the log probability of data between successive K values, using the program structure-sum (Evanno et al., 2005).
To assess the genetic similarity between pairs of populations, we estimated the linear function FST / (1 – FST) (Rousset, 1997) among all population pairs using the R library gstudio (Dyer, 2014) and visualized the relationships in a neighbor-joining (NJ) dendrogram (Saitou and Nei, 1987) with 300 bootstraps using the program ape v.3.1.0 (Paradis et al., 2004) implemented in R. We also estimated the genetic distance between pairs of populations with Nei's genetic distance (Nei, 1972).
To analyze the existence of a pattern of isolation by distance, we performed a Mantel test (Mantel, 1967) between the geographic distance and the linear function FST / (1 − FST) using the vegan package (Oksanen et al., 2013) of R. Geographic distances were calculated as points on a sphere with the raster library (Hijmans, 2014) of R.
Environmental differentiation of the subspecies of Agave striata
To assess the association between bioclimatic variables and genetic differentiation for each subspecies of A. striata, we used bioclimatic variables obtained from the WorldClim database (http://www.worldclim.org; Hijmans et al., 2005). WorldClim comprises 19 climate layers generated by interpolation of climate data from weather stations around the world. The variables employed are raster grids with a spatial resolution of 30 inches (ca. 1 km2 resolution). To compensate for the possible over-parameterization of the models due to the multicollinearity of the numerous climate layers, we used a principal component analysis (PCA) to select the noncorrelated variables. As a result of this analysis, we selected seven bioclimatic variables. We extracted the climate information from the six analyzed populations of each subspecies throughout the range of A. striata (Table 1). From sampled populations, we obtained climatic information and applied a linear adjustment.
To generate potential distributional models for each subspecies that could explain the patterns of genetic variation and genetic structure, we obtained occurrence points from the MEXU herbarium, the literature, and fieldwork at the Mexican states of Coahuila, Hidalgo, Nuevo León, Querétaro, San Luis Potosí, and Tamaulipas. We compiled a total of 126 locality records, 78 for A. striata subsp. striata and 48 for A. striata subsp. falcata (see Appendix S2 in online Supplemental Data). To reduce the effects of sampling bias, we removed the localities with similar coordinates to obtain unique records, and we selected the localities with 10 or more kilometers of distance among them.
Modeling algorithm
Environmental niche models are typically used to model the potential distribution of a single or multiple species (Illoldi-Rangel et al., 2004; Peterson, 2006). These tools have been very useful to evaluate the role of environmental conditions within species or species complexes (Miller and Knouft, 2006; Leaché et al., 2009; Hufford et al., 2012). In our case, the model for each subspecies helped us to detect potential climate similarities or differences between them as well as to identify the associations between the environmental conditions and the genetic differentiation within the species and between the subspecies.
We built the distribution models using the program MaxEnt v. 3.2.1 (Phillips et al., 2006). It takes into account only the presence records and samples “pseudo-absences” to calculate the model performance measures, an approach that is appropriate for our study because we only included presence data. We used 70% of the records as training data, the remaining 30% for testing the model, and 10 percentile training presence to edit and to define the threshold to validate the models.
To evaluate the accuracy of the models, we used a proportion test (Cruz-Cárdenas et al., 2013) that involves quantifying the test records with prediction probability values higher than the previously defined threshold (10 percentile training presence). We used the default settings for all analyses to measure the degree to which the models differed from random distributions using the AUC, the area under the receiver operating characteristic curve. The AUC scores are interpreted as reflecting the ability of the model to distinguish presence data from background data (Phillips et al., 2006).
MANOVA
We used a multivariate analysis of variance (MANOVA) to test for significant differences between the estimates of the mean for each environmental variable in each subspecies. We extracted climate data for specimen point localities directly from the seven bioclimatic layers (see online Appendix S3). Because the data violated normality and homoscedasticity assumptions, we conducted a nonparametric MANOVA on the log10-transformed data set (Anderson, 2001). We performed a series of ANOVAs and paired comparisons to assess which variable accounted for the overall difference detected by the nonparametric MANOVA with SPSS 15.0 (2006; SPSS, Chicago, Illinois, USA).
RESULTS
Genetic variation
In each of the subspecies of A. striata, we found 47 loci (bands), ranging from 300 to 2000 base pairs (bp) long. For A. striata species, the populations had high levels of genetic diversity (HS = 0.2448 ± 0.005 SD; %P = 80.85) (Table 1). Genetic diversity was higher for the populations of the northern A. striata subsp. falcata (HS = 0.3444 ± 0.007; %P = 80.85) than for the populations of the southern subspecies, A. striata subsp. striata (HS = 0.2131 ± 0.005; %P = 70.28). Including all the populations of both subspecies, genetic variation increased with latitude (Fig. 2A) and temperature seasonality (Fig. 2B), but decreased with annual precipitation (Fig. 2C); nevertheless, the correlations are mainly due to differences between subspecies, as explained later in the MANOVA section.
Correlations between genetic variation and bioclimatic variables. (A) Correlation between expected heterozygosity per population (HS) and latitude. (B) Correlation between heterozygosity per population and temperature seasonality. (C) Correlation between heterozygosity per population and annual precipitation (mm). Diamonds represent subspecies striata and squares represent subspecies falcata.
Agave striata is a species with relatively high population structure (Ѳ = 0.3164 ± 0.011 SD). We found significant differences in the allelic frequencies among populations in 40 of the 47 loci. The southern subspecies, A. striata subsp. striata, had on average higher population differentiation (Ѳ = 0.4085 ± 0.010) than northern subspecies, A. striata subsp. falcata (Ѳ = 0.1362 ± 0.011). The results of the AMOVA showed that 75.47% of the genetic variance of A. striata occurs among populations (ΦST = 0.24529), whereas only 6.99% of the variance is found among subspecies (see online Appendix S4).
The structure analyses showed two main population clusters (K= 2) based on the ∆K statistical test (Evanno et al., 2005) (Fig. 3A), one for each subspecies, with some levels of admixture in some populations. We observed an area of overlap between the distributions of the two subspecies that is located from southern Nuevo León to southern San Luis Potosí. Apparently, in this area there is high gene flow between the two subspecies (i.e., between genetic pools) (Fig. 3A).
Genetic differentiation within Agave striata. (A) Results of an analysis conducted with structure v. 2.2 (Pritchard et al., 2000; Falush et al., 2003, 2007). The green bars indicate that the individual belongs to group 1 (subspecies striata), and red bars indicate that the individual belongs to group 2 (subspecies falcata). The populations are located from left to right and from highest to lowest latitude. (B) Neighbor-joining of genetic distances measured as FST / (1 − FST) (Rousset, 1997) between populations. The number in each node indicates bootstrap value.
The average genetic distance for all A. striata populations was D = 0.0830, but D was higher (0.099) when only subspecies striata was considered and lower (0.039) when only subspecies falcata was considered. The NJ dendrogram (Fig. 3B) shows that the relationship between the two subspecies was resolved. Two admixed populations of A. striata subsp. falcata (F12 and F13) are in a cluster, while A. striata subsp. striata was assigned to another cluster.
Using a Mantel test (Fig. 4; see online Appendices S5 and S6), we found evidence of isolation by distance in the total species sample and also within both subspecies (A. striata: r = 0.556, P = 0.0005; A. striata subsp. striata: r = 0.631, P = 0.003; A. striata subsp. falcata: r = 0.786, P = 0.005).
Mantel test analysis to evaluate isolation by distance within each subspecies of Agave striata.
Modeling subspecies distribution
The validation of the model for each subspecies (see online Appendix S7) indicated that the obtained models were significantly better than the random distributions because the bivariate test and the AUC values were well supported.
The inspection of the models for each subspecies shows a potential differentiation between them and an area of contact. Based on the model and on the biogeographical provinces of Mexico, the distribution model of A. striata subsp. striata is mainly restricted to the southern region of the Chihuahuan Desert and to the Sierra Madre Oriental (Fig. 5). This area corresponds to the biogeographical provinces of Altiplano del Norte, Sierra Madre Oriental, and some portions of the Mexican Gulf and Tamaulipeca provinces. The distribution model of A. striata subsp. striata overlaps in its northwestern distribution with the other subspecies (Fig. 5). In the case of A. striata subsp. falcata, the model indicated a potential distribution including the central and northwestern area of the Chihuahuan Desert (Fig. 5A), corresponding to the provinces of the Altiplano Sur (Zacatecano-Potosino) and Altiplano Norte. The predicted distribution of this subspecies partially overlaps with the subspecies striata in the southeastern portion of its distribution (Fig. 5).
Potential distribution models of the subspecies of Agave striata and the biogeographical provinces of Mexico. Distribution model of both subspecies falcata (yellow) and striata (green) and the zone of overlap in dark green.
MANOVA
The results of the multivariate analysis showed that the subspecies differ in their climatic conditions (Wilks lambda P < 0.0001, F = 210390.1, df = 7). Five climatic variables were significantly different between subspecies (Table 2). These variables are isothermality, temperature seasonality, annual precipitation, precipitation seasonality, and precipitation of the driest quarter. These five variables have higher means in the subspecies striata than in the subspecies falcata, with the exception of temperature seasonality (Table 2), indicating that these are the most relevant parameters accounting for the niche differences between the two subspecies.
| MANOVA | Value | df | F | P |
|---|---|---|---|---|
| Wilks lambda | 0.305 | 7 | 32.556 | <0.001 |
| Univariate contrast |
| Variables | SS | df | F | P |
|---|---|---|---|---|
| Bio 1 = Annual mean temperature (°C) | 0.000 | 1 | 0.049 | n.s. |
| Bio 3 = Isothermality | 0.052 | 1 | 24.386 | <0.001 |
| Bio 4 = Temperature seasonality | 0.650 | 1 | 67.198 | <0.001 |
| Bio 12 = Annual precipitation (mm) | 1.029 | 1 | 53.104 | <0.001 |
| Bio 15 = Precipitation seasonality (coefficient of variation) | 0.042 | 1 | 9.735 | 0.002 |
| Bio 17 = Precipitation of driest quarter (mm) | 0.488 | 1 | 19.516 | <0.001 |
| Elevation (m) | 0.080 | 1 | 1.956 | n.s. |
| Pairwise comparisons | SS | df | F | P |
| Variables | (I) Species | (J) Species | Differences between means (I – J) | |
| Bio 1 = Annual mean temperature (°C) | 1 | 2 | −0.002 | n.s. |
| Bio 3 = Isothermality | 1 | 2 | −0.045 | <0.001 |
| Bio 4 = Temperature seasonality | 1 | 2 | 0.159 | <0.001 |
| Bio 12 = Annual precipitation (mm) | 1 | 2 | −0.200 | <0.001 |
| Bio 15 = Precipitation seasonality (coefficient of variation) | 1 | 2 | −0.041 | 0.002 |
| Bio 17 = Precipitation of driest quarter (mm) | 1 | 2 | −0.138 | <0.001 |
| Elevation (m) | 1 | 2 | 0.056 | n.s. |
- Notes: n.s. = not significant, SS = sum of squares.
DISCUSSION
Originally, the two subspecies of A. striata were described as different species. Agave striata was described by Zuccarini in 1832 (cited by Thiede, 2001, p. 66), while A. falcata was described by Engelmann in 1875 (cited by Thiede, 2001, p. 66). Mainly based on their general similarities and distribution patterns, Gentry (1982, pp. 242–247) proposed that they actually belonged to the same species, each taxa representing a subspecies, A. striata subsp. falcata being found in the North and A. striata subsp. striata in the South. Nevertheless, the two subspecies show differences in the number, the length, and especially the width of the leaves (which are wider in A. striata subsp. falcata than in A. striata subsp. striata). In the contact zone between the two subspecies, the morphological differences between them are less evident (Gentry, 1982). Our genetic differentiation and environmental analyses are consistent with this observation.
Each subspecies of A. striata has contrasting patterns of genetic diversity and structure. The southern subspecies striata has less genetic variation in each population and shows stronger differentiation, as measured by both Ѳ and D, than the northern subspecies falcata. Accordingly, each A. striata subspecies is mainly represented by a single gene pool in our structure analysis. Nevertheless, the same analyses suggest genetic similarity and perhaps gene flow in the contact zone between the two subspecies. This pattern can also be generated by an active and ongoing differentiation (and perhaps speciation) process because the isolation by distance analyses indicate that the more geographically (and in general terms, climatically, see Alvarado-Cárdenas et al., 2013) distant populations have a higher genetic differentiation and that the differentiation is stronger for the southern subspecies striata.
Agave striata and its subspecies have some of the largest genetic distances among populations reported for Agave (Navarro-Quezada et al., 2003; Rocha, 2006; Aguirre-Dugua and Eguiarte, 2013; Eguiarte et al., 2013) and have high genetic diversity, although close to the average found for the genus Agave (Eguiarte et al., 2013). For instance, A. lechuguilla has a very wide distribution in the Chihuahuan Desert and adjacent areas, from central Mexico to the southern United States and, according to isozyme analyses, it has high levels of genetic diversity within populations (HS = 0.394) but low population structure (Ѳ = 0.083) (Silva-Montellano and Eguiarte, 2003), lower than what we found in this study for the population structure of the subspecies falcata (Ѳ = 0.1362). In A. victoriae-reginae, an endemic species with narrow distribution in the Chihuahuan Desert, a study based on isozymes found high genetic diversity (HS = 0.33), and high genetic structure (Ѳ = 0.236) (Martínez-Palacios et al., 1999), but a lower differentiation than the one we found in A. striata subsp. striata (Ѳ = 0.4085). These comparisons suggest that the genetic differentiation in Agave does not depend only on the extent of its distribution and climate, but it also depends on the complex evolutionary history of its populations. For instance, Scheinvar E. et al. (Instituto de Ecología, UNAM, unpublished manuscript), using chloroplast information found four different genetic groups in A. lechuguilla that originated in different glacial–interglacial cycles of the Pleistocene, with a recent expansion to the north from a southern Last Glacial Maximum refuge. The southern refuge of A. lechuguilla contrasts with the possible northern refuge of A. striata, but subsequent historical and phylogeographic analyses are required.
Models of potential distribution showed that each of the subspecies occur in different environmental conditions, although they have an area of overlap in part of their range. This result is consistent with what was observed in the genetic analysis. The MANOVA results also showed that there are differences between the potential distributions of the subspecies in five of the climatic variables. Environmental conditions associated with the heterogeneity of the temperature and annual precipitation apparently play an important role in the genetic differentiation of both taxa, which is consistent with other studies (Foll and Gaggiotti, 2006; Kiambi et al., 2008; Hamasha et al., 2013).
The southern subspecies, A. striata subsp. striata, is mainly distributed in areas with higher precipitation and lower seasonality in temperature, but higher precipitation seasonality (Table 2; online Appendix S3). In consequence, this subspecies grows in different types of vegetation than the other subspecies. The northern subspecies, A. striata subsp. falcata, lives in a similar interval of temperature, but with lower values of annual precipitation, higher temperature seasonality but lower precipitation seasonality, and it is associated with drier vegetation, which characterizes the northern part of the Chihuahuan Desert (Table 2; online Appendix S3). The greater genetic variation in this subspecies correlates with a stressful environment. This genetic variation could be associated with different adaptation strategies, and in this case the genetic variability may be useful to relieve the effects of drastic conditions, as found for organisms growing in similar conditions, but further analyses will be necessary (Nevo, 2001; Ferrer et al., 2004; Hamasha et al., 2013).
In addition to the climatic conditions, the distribution of both subspecies in different portions of the Chihuahuan Desert and in the biogeographic province of the eastern Sierra Madre Oriental, as well as their complicated geography that includes mountain ranges and canyons, could be also acting as a physical barrier to partially separate the two subspecies (Fig. 1C), limiting gene flow among them, as has been observed in other plants (Hensen et al., 2011; Hamasha et al., 2013). These topographic barriers, along with the temperature and precipitation gradient existing from north to south of the Chihuahuan Desert, seem to be potential detonators of the geographic distribution of the subspecies of A. striata. Nevertheless, the patterns are complex because, while some populations around the contact area of the two subspecies clearly show signals of admixture, as is the case for F11, F12, F13 of A. striata subsp. falcata and F15 and F16 of A. striata subsp. striata, other populations of the same area such as F8 of A. striata subsp. falcata or F14 of A. striata subsp. striata show little or no signal of gene flow (Figs. 1 and 2). It will be important to analyze whether the latter populations have differences in ecology, reproductive ecology (phenology, pollinators), or physical environmental conditions that result in a decreased admixture.
Since the proposal of Gentry (1982), both taxa of A. striata have remained within the subspecies category because of their general morphological similarity. However, the detailed evaluation of their genetic and environmental attributes, as well as observations of their morphology, allow us to suggest that both entities may be in an active process of differentiation and perhaps of speciation. The separation of different taxonomic entities based on genetic and environmental attributes have received more attention lately (Leaché et al., 2009). Future detailed morphological (including leaves, inflorescences, flowers, fruits, and seeds) and physiological studies along with detailed gene flow, pollination, and reproductive ecology analyses of both subspecies across their distribution ranges, in particular form the overlap area of the subspecies, will help us to understand whether we truly have an active ongoing speciation process and whether we have a primary or secondary contact zone between diverging lineages. At the moment, we consider that Gentry's (1982) proposal of a single species with two subspecies seems to be consistent with our data and analyses because we found a series of intermediate populations (both in morphology and in the genetic markers) that suggest current gene flow between the two main gene pools. We also found correlations between bioclimatic variables and genetic differentiation in Agave striata. Although these correlations do not prove a causal relationship, they provide a hypothesis to explore the effects of temperature and precipitation in Agave.
ACKNOWLEDGEMENTS
This study was funded by SEMARNAT-CONACYT grant 0246 and CONACYT-SEP grant 46475. We acknowledge the Posgrado en Ciencias Biológicas, UNAM and CONACYT for scholarships to L. Trejo. We thank A. Varela and L. Espinosa for help in the laboratory; J. Gasca and A. González for support in the analysis; A. Silva, N. Flores and R. González for support in the field; and M. E. Olson, G. Castellanos-Morales, and E. Aguirre-Planter and two anonymous reviewers for their comments on the manuscript.
