Different patterns of colonization of Oxalis alpina in the Sky Islands of the Sonoran desert via pollen and seed flow

Abstract Historical factors such as climatic oscillations during the Pleistocene epoch have dramatically impacted species distributions. Studies of the patterns of genetic structure in angiosperm species using molecular markers with different modes of inheritance contribute to a better understanding of potential differences in colonization and patterns of gene flow via pollen and seeds. These markers may also provide insights into the evolution of reproductive systems in plants. Oxalis alpina is a tetraploid, herbaceous species inhabiting the Sky Island region of the southwestern United States and northern Mexico. Our main objective in this study was to analyze the influence of climatic oscillations on the genetic structure of O. alpina and the impact of these oscillations on the evolutionary transition from tristylous to distylous reproductive systems. We used microsatellite markers and compared our results to a previous study using chloroplast genetic markers. The phylogeographic structure inferred by both markers was different, suggesting that intrinsic characteristics including the pollination system and seed dispersal have influenced patterns of gene flow. Microsatellites exhibited low genetic structure, showed no significant association between geographic and genetic distances, and all individual genotypes were assigned to two main groups. In contrast, chloroplast markers exhibited a strong association between geographic and genetic distance, had higher levels of genetic differentiation, and were assigned to five groups. Both types of DNA markers showed evidence of a northward expansion as a consequence of climate warming occurring in the last 10,000 years. The data from both types of markers support the hypothesis for several independent transitions from tristyly to distyly.

changes combined with information derived from genealogies can contribute to a better understanding the evolutionary history of a reproductive systems of plants. Recent studies have investigated the influence of historical processes such as vicariance events related to climatic changes on the evolutionary dynamics of the plant breeding systems (Dorken & Barrett, 2004;Hodgins & Barrett, 2007;Pérez-Alquicira et al., 2010;Zhou, Barrett, Wang, & Li, 2012). Furthermore, selective pressures associated with the immobility of plants and their reliance on pollen vectors have played an important role in the expression of breeding system variability (Barrett, 2013).
Regions with climatic oscillations that have impacted distribution of organisms represent natural laboratories for investigation of how historical processes such as bottlenecks, isolation, genetic drift, and natural selection have produced complex configurations of population structure and influenced the evolution of reproductive systems (Boyd, 2002;Masta, 2000;Sosenski, Fornoni, Molina-Freaner, Weller, & Dominguez, 2010). An area that has experienced drastic changes in the configuration of the ecosystem is the current desert of northwestern Mexico and southwestern USA. During the last glacial maximum (18,000 years ago), the Sonoran and Arizona deserts experienced an average decrease in temperature of 6°C and increase in humidity (Metcalfe, 2006), fostering an extensive expansion of woodland vegetation. During the last 10,000 years, the climate in these areas became warmer and drier, producing a northward and upward range shift of cool-temperate species. Currently, these species are restricted to the tops of the isolated mountains surrounded by lower elevation desert, known as the Sky Islands of the Sonoran Desert, which serve as habitat for cool-temperate species. On these  Silva-Montellano & Eguiarte, 2002). Several studies have also detected similar latitudinal gradients and western colonization during the Holocene Epoch in species inhabiting Europe (Conord et al., 2012).
In order to have a general overview of the evolutionary history of population structure, molecular markers with both cytoplasmic and nuclear inheritance should be used. Because of the differences in the inheritance of nuclear and cytoplasmic genes (chloroplast and mitochondria), the patterns of genetic structure are frequently dissimilar. On one hand, the dispersion of nuclear genes occurs via pollen and seeds, while the cytoplasmic genes are maternally inherited in most angiosperms, and thus, dispersal of cytoplasmic genes depends exclusively on seed movement (Ennos, 1994;Sears, 1980). Gene flow through seeds is usually more restricted in comparison with pollen flow, and population structure should be compared using plastid and nuclear loci that are expected to differ in the magnitude of gene flow. In this study, we performed a comparison of genetic structure through cpDNA (results previously published in Pérez- Alquicira et al., 2010) and nuclear microsatellite markers of alpine wood sorrel, Oxalis alpina ( Figure 1).
Oxalis alpina (Rose) Knuth (section Ionoxalis, Oxalidaceae) is a heterostylous, tetraploid, perennial herb inhabiting evergreen Madrean woodlands in the Sky Island region (Weller, Dominguez, Molina-Freaner, Fornoni, & LeBuhn, 2007). Heterostyly is a floral polymorphism where two (distyly) or three (tristyly) floral morphs occur in populations. Tristyly includes short-, mid-, and long-styled morphs; the mid-styled morph is absent in the distylous populations. Distyly is the derived breeding system in Oxalis (Gardner, Vaio, Guerra, & Emshwiller, 2012;Weller & Denton, 1976;Weller et al., 2007). Each floral morph is characterized and named by the position of the stigma relative to the two levels of stamens ( Figure 2). Typically, heterostylous incompatibility systems prevent self-fertilization and fertilization between stigmas and anthers that are located at different heights (illegitimate crosses, those crosses that do not produce seeds, sensu Darwin, 1877) (Barrett, 1992). For example, legitimate crosses occur between the short stigma and pollen from the short stamens of the long-and midstyled morphs. Legitimate crosses of the long-styled morph occur when pollen comes from the long stamens of the short-and midstyled morphs, and the midstigma produce seeds after receiving pollen from the mid stamens of the short-and long-styled morphs (Figure 2; Barrett, 1992).
Oxalis alpina exhibits remarkable variation in the reproductive system when viewed in the historical and geographic context of the Sky Island region. First, tristyly occurs mainly in southern ranges while the derived distylous reproductive system occurs in the northern ranges. Great variation in the frequency of the floral morphs has been observed among tristylous populations in the Sky Islands . Second, the incompatibility system includes different degrees of modification, with remarkable modifications in the northern ranges of the distribution. Weller et al. (2007) found a negative relationship between the extent of incompatibility modification in the short-and long-styled morphs and the frequency of the midstyled morph, suggesting that the modifications of the F I G U R E 1 Oxalis alpina (Rose) Knuth growing in the Chiricahua Mts., southeastern Arizona incompatibility system have influenced the loss of midstyled morph in distylous populations (Weber et al., 2013). As the mid-styledmorph becomes less frequent, floral morphology of the remaining short-and long-styled morphs more closely resembles the morphological patterns of distylous populations (Sosenski et al., 2010).
Historical factors associated with range expansion and contraction during climatic oscillations in the Pleistocene Period might also have influenced the tristyly-distyly transition (Pérez-Alquicira et al., 2010). We previously demonstrated that northern populations exhibited the lowest levels of genetic diversity suggesting this area was recently colonized by O. alpina. Thus, genetic drift associated with founder events in northern ranges might have influenced the loss of midmorph in populations that already have low frequencies of the midstyled morph (Pérez-Alquicira et al., 2010). Studying the evolutionary transitions of the breeding system of O. alpina in the historical and geographic context of the Sky Islands provides a unique opportunity to capture the evolutionary steps in the tristyly-distyly transition. In this study, we (1) analyze the influence of climatic oscillations on the genetic structure of O. alpina using microsatellite markers and propose phylogeographic scenarios to aid in understanding of the evolutionary processes behind the variation in the reproductive system of O. alpina and (2) compare the phylogeographic patterns inferred from the microsatellite data with previous findings from cpDNA data to reconstruct the historical scenarios that better explain the current patterns of genetic structure of O. alpina populations.

| Study species
Oxalis alpina occurs from southwestern USA to Guatemala. Based on molecular evidence, the North American Oxalis section Ionoxalis species colonized this area in two events from ancestral southern South America species (Gardner et al., 2012). Furthermore, F I G U R E 2 Reproductive system of Oxalis alpina including three floral morph and their possible genotypes. The genetic system controlling tristyly (remove coma), consists of two linked loci, each with two alleles. The S locus is epistatic over M. The presence of the dominant S allele results in the expression of short-styled morph, while the dominant M allele produces the expression of midstyled morph, and when both loci are homozygous, the long-styled phenotype is expressed. The dotted lines indicate the legitimate crosses, those leading to substantial seed production

| DNA extractions, amplifications, and microsatellite genotyping
Genomic DNA was isolated from 100 mg of leaf tissue following Doyle and Doyle CTAB procedure (1987) and also using the DNeasy

| Genetic structure
The levels of genetic differentiation among populations measured using the parameters F ST and R ST were estimated using SpaGeDi software version 1.3a (Hardy & Vekemans, 2002). We also obtained the F ST values for distylous and tristylous populations and compared our results with cpDNA data previously published (Pérez-Alquicira et al., 2010). To test for phylogeographic pattern, we used the approach proposed by Hardy, Charbonnel, Fréville, and Heuertz (2003), which is based on the comparison of the observed R ST (before randomization) and an expected value (pR ST ) calculated after 5,000 allele size permutations using SPAGeDi version 1.3a (Hardy & Vekemans, 2002). This test can be interpreted as testing whether If R ST is significantly larger than pR ST, a stepwise mutation model is the most likely explanation for genetic differentiation (Hardy & Vekemans, 2002) and a phylogeographic pattern is suggested. If the difference between R ST and pR ST is not significant, this suggests that allele size is not as important because mutations do not follow a stepwise mutation model, the absence of phylogeographic pattern can be inferred, and F ST should be used instead of R ST .
To determine how many distinct population clusters were supported by the data, we used the software STRUCTURE 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). This program uses a Bayesian method to determine the probability for each individual to be assigned to a particular cluster regardless of its geographic location. We ran the analyses using the admixture model and correlated frequencies with a burn-in and run length of 250,000 and 1,000,000, respectively, k = 1 to 10 (k indicates the number of genetic clusters), and 10 iterations for each k. Further analyses of the substructure within each main genetic cluster were analyzed using the same parameters as for the entire data set (burn-in and run length of 250,000 and 1,000,000, respectively, k = 1 to 10, and 10 iterations for each k).
We did not include the population origin as a "prior" for the analyses. The results were imported into Structure Harvester (Earl & vonHoldt, 2012), which allows the assessment and visualization of the likelihoods for each k value and detects the number of genetic clusters that best fit the data, based on Delta K (Evanno, Regnaut, & Goudet, 2005).
We further carried out a principal coordinate analysis (PCoA) using the Bruvo distances for visualizing the genetic structure and relationships among samples through the R program POLYSAT (Clark & Jasieniuk, 2011). A neighbor-joining (NJ) phenogram was constructed based on Nei's genetic distances to visualize the genetic relationship among populations. We obtained the allelic frequencies through the program POLYSAT (Clark & Jasieniuk, 2011) and then generated 1,000 replicates by a bootstrap resampling method using the seqboot program included in the Phylip package version 3.695 (Felsenstein, 1989 a priori assumptions about population structure. The method is based on the genetic covariance structure among populations analyzed simultaneously (Dyer & Nason, 2004). Populations that exhibit significant genetic matrix correlation will be connected in the network by edges (lines), and the length of the edges is inversely proportional to the genetic covariance between the populations. Therefore, longer edges indicate lower genetic covariance between populations. Populations that are not connected indicate the absence of migration, and the presence of subgraphs (a smaller network within a large network) indicates that a population or group of populations maintain a weak or null genetic connection (Dyer, 2007;Dyer & Nason, 2004;Dyer, Nason, & Garrick, 2010).
In addition, we tested for a correlation between the matrices of genetic distances (Nei, 1972) and the midmorph frequencies using a Mantel test (ADE4 package Dray & Dufour, 2007, in R). The latter test was conducted to determine whether populations with a reduced or zero (distylous population) frequency of the mid morph are genetically more similar. A significant correlation would suggest that distyly evolved from the same tristylous ancestors. The mid morph frequency matrix was constructed using the program Phylip version 3.695 (Felsenstein, 1989).

| Comparison of genetic structure between cpDNA and microsatellite markers
Because one of our objectives was to compare patterns of genetic diversity using nuclear and chloroplast markers, we performed regression analyses between haplotype diversity from cpDNA (chloro- We also estimated the pollen-seed migration ratio (r = mp/ms, where mp and ms are migration values for pollen and seeds, respectively) based on equations from Ennos (1994) and further modified by Petit et al. (2005): G ST corresponds to the genetic differentiation statistic; the subindices b and m correspond to biparental and maternal inheritance markers, respectively, and F IS corresponds to the heterozygote deficit estimated with nuclear codominant markers. Both G ST and F IS were estimated using the SpaGeDi software version 1.3a (Hardy & Vekemans, 2002). Overall, larger values for r indicate that gene flow by pollen is quantitatively more important than gene flow via seeds.

| Genetic variation
Testing for clonal structure revealed the presence of 17 clones distributed in eight populations with the number of plants per clone varying from 2 to 5. For further analyses, we included only genets (one genotype per clone), and therefore, a total of 29 samples were removed from the analyses. The ratio of the number of genets to the number of ramets was from 0.73 to 0.99 (Table 2).
Two distylous (Santa Rita and Sierra Ancha) and one tristylous  (Table 3). For cpDNA sequences, latitude had a significant effect (p = .03) on haplotype diversity. For both genetic markers, longitude and the interaction of latitude and longitude did not show any significant effect (Table 3). The remaining parameters of genetic diversity for microsatellites did not show any significant association with latitude and longitude.

| Genetic structure and genetic relationship among populations
We found F ST = 0.27 and R ST = 0.47 for 17 populations of O. alpina.
The genetic structure for distylous and tristylous populations was similar (F ST = 0.27 and 0.28, respectively). The R ST parameter was marginally significantly higher than the pR ST (0.28) (p = .053). The TA B L E 2 Genetic and clonal diversity characteristics, and bayesian assignment (clusters I and II) based on microsatellite genetic markers for 17 populations of Oxalis alpina from the Sky Island region   (Figures 3a and 4a). The PCoA analysis reflects also the presence of two main clusters, supporting the Bayesian assignment result (Figure 3b). Principal components 1 and 2 were used to construct the PCoA plot because they represent most of the total variation in the genetic data (55% and 6%, respectively). The F ST values for clusters I and II were similar, 0.09 and 0.10, respectively, and for the parameter R ST, the values were 0.09 and 0.059, respectively.
We also examined the genetic structure within each cluster ( Figure 4). Four groups were present in cluster I (Figure 4b).
Approximately 35% of the samples could not be assigned to any group, which indicates the mixed genetic nature of a large number of individuals, while 26%, 17%, 8%, and 13% of the samples were as-   Finally, the correlation between pairwise genetic distances and the frequency of the midmorph was not statistically significant (r = .03, p = .18).

| Comparison of genetic structure between cpDNA and microsatellites markers
Regression analyses did not show any relationship between the microsatellite diversity parameters and cpDNA haplotype diversity

| D ISCUSS I ON
For most angiosperms, nuclear genes are inherited paternally via pollen and maternally via seeds, while cytoplasmic genes found in the chloroplast and mitochondria are maternally inherited (Petit, Kremer, & Wagner, 1993). Complex configurations of gene flow and differences in the distribution of genetic variability within and among populations are expected through nuclear and chloroplast markers (Petit et al., 2005). Evidence from several species of angiosperms and gymnosperms support this prediction (Petit et al., 2005). the GARP algorithm (Genetic Algorithm for Rule-set Production) (Stockwell & Peters, 1999), and the results indicated that during the LGM (Last Glacial Maximum, 18,000 years ago), O. alpina was distributed mainly at lower elevations of the southeastern portion of the Sky Island region. Therefore, northwestern migrations occurred once the climate becomes warmer. Similar patterns of migration have been detected for plant and animal species inhabiting southern Europe and Northern Africa at the end of the Pleistocene period (Davis & Shaw, 2001;Petit et al., 2003;Schmitt, 2007 (Barber, 1999;Downie, 2004;Masta, 2000). However, because eight microsatellite loci were included in our analyses, convergence seems unlikely to explain all of these similarities.
For example, Navascués and Emerson (2005) performed simulations to understand the potential effects of homoplasy on determining genetic relationships among populations. They found that as the number of loci in a study increases, the bias for the estimation of genetic parameters as a consequence of homoplasy decreases. The simulations indicated that four loci produced higher values for the homoplasy index and greater bias for inferring genetic relationship among populations than simulations using nine loci. We included eight loci in our study, and thus, bias in the inferring phylogeographic structure via homoplasy is unlikely but not impossible.
Because dispersion by pollen frequently occurs over larger distances than by seeds, lower levels of genetic structure were expected for nuclear genes in comparison with cytoplasmic genes (Birky et al., 1989;Petit et al., 1993Petit et al., , 2005 (Ennos, 1994;Petit et al., 2005), indicate that gene flow through pollen is quantitatively much more important than seeds. For O. alpina, r = 5.62, which is smaller than the median value of 17 for 93 angiosperms species (Petit et al., 2005). This lower value for the pollen/ seed migration ratio is due to the relatively high G ST value, which is necessary for calculating r (see Petit et al., 2005  Divergent patterns of genetic structure might also be explained by the more rapid coalescence of uniparentally inherited haploid alleles in mitochondrial or chloroplast DNA that have smaller effective population sizes than most nuclear loci (Birky et al., 1989;Palumbi et al., 2001). Smaller effective population size for organellar DNA accelerates the processes of genetic drift in neutral markers, producing more rapid genetic divergence among populations (Palumbi et al., 2001). These processes may in part explain the strong subdivision observed through cpDNA in O. alpina, where five genetic groups of populations were detected.
One of the main results from the popgraph analyses is that the six distylous populations included in this study were assigned to two subgroups and were not all directly connected by gene flow. This is also supported by the results from the Bayesian assignment, neighbor-joining phenogram, and PCoA. According to the population network, three populations (two distylous and one tristylous) from the Chiricahua Mts in cluster I are genetically similar, which is also supported by the results from the Structure and NJ phenogram. These three populations are also located in close geographic proximity to each other, representing a remarkable example of the transition from tristyly to distyly. Weller et al. (2007) found that the Chiricahua Mountains populations harbor a northern-southern gradient of mid-styled morph frequency with some populations having reduced frequency of mid-styled morphs, suggesting that the evolution to distyly is an ongoing process within this geographic range.
Our data for this area seem to provide further support for this hypothesis. However, in the three distylous populations included in the cluster II (Sierra Ancha, Pinal, Pinaleño), distyly probably evolved independently from the transition in the Chiricahua Mts.
We did not find a significant association between genetic distances and mid-styled morph frequency, indicating that distylous populations are genetically dissimilar, and therefore distyly probably originated from different tristylous populations. Additionally, genetic evidence from chloroplast data further suggests that distyly has evolved more than once (Pérez-Alquicira et al., 2010). If distyly is the result of convergent evolution, we expect that deterministic evolutionary forces such as natural selection have played an important role in the evolution of distyly. Modifications of the tristylous incompatibility system support this result. For example, Weller et al. (2007) found that the loss of key elements of the heterostylous incompatibility system in long-and short-styled morphs increases the degree of cross compatibility between these morphs and leads to genic selection against mid alleles. Populations with higher levels of genic selection against the midallele also have higher levels of selffertilization for mid-styled morphs (Weber et al., 2013). Selfing in the mid morph will result in expression of inbreeding depression which may lead to further decline in the frequency of this morph. As these deterministic forces lead to reduced frequency of the midmorph, this morph will be increasingly sensitive to genetic drift, especially with low population sizes. We propose that once populations migrated to northward ranges as a consequence of climatic oscillations, those populations with low frequencies of mid-styledmorph were more likely to lose this reproductive morph due to genetic drift. This process could explain partially the distribution of distylous population in northern ranges of the Sky Island region.
Overall, our results indicate that microsatellite markers provide useful genetic information for tracing the possible routes of colonization of populations of O. alpina in the Sky Island region.
The occurrence of southern and northern populations in cluster II, together with a marginally significant reduction in genetic diversity toward northern areas (H E parameter), support that the route of colonization was probably from south to north, although the genetic signature of northern founder events inferred from microsatellites is not as clear as from cpDNA markers (Pérez-Alquicira et al., 2010). Niche modeling also showed that the distribution of O. alpina during the last glacial maximum (18,000) occurred mainly in the southern ranges compared with the current distribution (Pérez-Alquicira et al., 2010). Northward migration has been a general pattern of colonization in section Ionoxalis because of the South America origin of this section, and species have dispersed on several occasions to North America (Gardner et al., 2012). In view of the reduced genetic differentiation for microsatellites compared to cpDNA genetic markers, which is characteristic of most angiosperms species (Petit et al., 2005), the detection of two genetic lineages via microsatellite markers compared to five genetic clusters via cpDNA is not surprising. Microsatellite markers revealed that distylous populations were associated with different tristylous populations, and this result is also supported by our previous findings using cpDNA suggesting that the tristyly-distyly transition occurred more than once within the Sky Island region.
These results suggest that natural selection has played an important role in the evolution of the reproductive system of O. alpina.
Moreover, because distylous populations are mainly located in northern ranges where founder events were frequent, genetic drift is likely to have influenced the evolutionary transition of tristyly-distyly system.

ACK N OWLED G M ENTS
We extend thanks to Rubén Pérez, Maureen Peters, and Greg Joice for technical assistance, and Andraca-Gómez G. for her support in map editions. JP-A. acknowledges the academic support received from Cátedras-CONACYT. Our project was supported by a National Science Foundation (NSF) grant DEB-0614164 (SGW and A. K. Sakai, co-PIs); a UC MEXUS (University of California Institute for Mexico and The United States) award to SGW and CAD; and NSF subaward and REU supplement RR715-061/4689108 to OVT.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R S CO NTR I B UTI O N S
All authors contributed to the design of the study and sample collection. JP-A and OVT acquired sequencing and microsatellite data, respectively, and conducted data analyses. JP-A, OVT, and SGW wrote the first draft of the manuscript, and all authors contributed critically to its revision and provided approval for the final submission.