Phylogeography and palaeodistribution modelling of Nassauvia subgenus Strongyloma (Asteraceae): exploring phylogeographical scenarios in the Patagonian steppe

The Patagonian steppe is an immense, cold, arid region, yet phylogeographically understudied. Nassauvia subgen. Strongyloma is a characteristic element of the steppe, exhibiting a continuum of morphological variation. This taxon provides a relevant phylogeographical model not only to understand how past environmental changes shaped the genetic structure of its populations, but also to explore phylogeographical scenarios at the large geographical scale of the Patagonian steppe. Here, we (1) assess demographic processes and historical events that shaped current geographic patterns of haplotypic diversity; (2) analyze hypotheses of isolation in refugia, fragmentation of populations, and/or colonization of available areas during Pleistocene glaciations; and (3) model extant and palaeoclimatic distributions to support inferred phylogeographical patterns. Chloroplast intergenic spacers, rpl32–trnL and trnQ–5′rps16, were sequenced for 372 individuals from 63 populations. Nested clade analysis, analyses of molecular variance, and neutrality tests were performed to assess genetic structure and range expansion. The present potential distribution was modelled and projected onto a last glacial maximum (LGM) model. Of 41 haplotypes observed, ten were shared among populations associated with different morphological variants. Populations with highest haplotype diversity and private haplotypes were found in central-western and south-eastern Patagonia, consistent with long-term persistence in refugia during Pleistocene. Palaeomodelling suggested a shift toward the palaeoseashore during LGM; new available areas over the exposed Atlantic submarine platform were colonized during glaciations with postglacial retraction of populations. A scenario of fragmentation and posterior range expansion may explain the observed patterns in the center of the steppe, which is supported by palaeomodelling. Northern Patagonian populations were isolated from southern populations by the Chubut and the Deseado river basins during glaciations. Pleistocene glaciations indirectly impacted the distribution, demography, and diversification of subgen. Strongyloma through decreased winter temperatures and water availability in different areas of its range.


Introduction
Phylogeographical reconstructions along large, extensive geographic areas are difficult to assess because they ideally need widely distributed, ubiquitous taxa, present at regional scale as a characteristic element of different local communities. Nassauvia Comm. ex Juss. subgen. Strongyloma (DC) Cabrera (Asteraceae) is one of these rare cases; it is in the field: it comprises small shrubs with a tendency to heteroblasty, with 1-5 white flowers per capitulum, and hairy cypselas (Cabrera 1982;Freire et al. 1993). Although Cabrera (1982) suggested wind dispersal, our personal observations in the field suggest that cypselas are dispersed mainly rolling on the ground, because the pappus is quickly deciduous. Distinguishing species within this complex is challenging; however, because variation in morphological traits does not strictly follow species boundaries as traditionally defined (Nicola et al. 2014). Maraner et al. (2012) found that subgen. Strongyloma is polyphyletic based on the ITS region of nuclear ribosomal DNA, but this appears to be an error in sample or sequence identity, or the use of divergent sequence paralogs (Nicola et al. 2014;Nicola, M. V., S. M. Sede, L. A. Johnson, and R. Pozne in preparation). Morphological evidence (Cabrera 1982;Freire et al. 1993;Nicola et al. 2014) and molecular data derived from plastid and nuclear ribosomal DNA regions (Nicola, M. V., S. M. Sede, L. A. Johnson, and R. Pozner in preparation) indicate subgen. Strongyloma is a well-supported monophyletic group with wide and continuous morphological variation coupled with little genetic divergence. We follow here the taxonomical decision of Nicola, M. V., S. M. Sede, L. A. Johnson, and R. Pozner (in preparation) by considering Nassauvia subgen. Strongyloma as a monospecific taxon. The morphological combinations traditionally accepted as species (Cabrera 1982;Freire et al. 1993) are here referred as "morphological variants" (i.e., axillaris, glomerulosa, fuegiana, maeviae, and ulicina morphologies), just to allow some reference points within the morphological, continuous variation observed in this group of plants (Nicola et al. 2014).
Nassauvia subgen. Strongyloma has a remarkably wide distribution, mainly in the Patagonian steppe with a few enclaves in the Andean region, from sea level to 4560 m a.s.l., in southern Bolivia (c. 21°S), in Argentina (from Jujuy Province c. 22°S to northern Tierra del Fuego Province c. 54°S), and in Chile (from Valpara ıso c. 32°S to Magallanes c. 51°S), including a variety of climatic zones but mainly adapted to extremely xeric environments (Cabrera 1982;Nicola et al. 2014). With a particularly abundant distribution, populations with axillaris and glomerulosa morphological variants characterized the Western and the Central Patagonian districts of the Patagonian Phytogeographical Province, respectively (Cabrera and Willink 1980). Sympatric assemblages of different combinations of two or three morphological variants of subgen. Strongyloma are frequent in some areas of Patagonia. With abundant morphologically intermediate individuals, such areas might represent secondary contact zones or populations with incomplete sorting of lineages (Nicola et al. 2014). The fossil record refers the oldest members of subtribe Nassauviinae back to the Miocene, when progenitors of this group were expanding in eastern Patagonia (Barreda et al. 2008). Nassauviinae, together with Mutisieae, are among the oldest, early branching groups within Asteraceae, after Onoseridae and Barnadesioideae (Funk et al. 2009).
The vast, flat, arid Patagonian steppe ecoregion is longitudinally placed from the east of the slopes of the southern Andes to the Atlantic coast, and latitudinally located from central Mendoza Province (c. 32°S) to northern Tierra del Fuego Province (c. 54°S) in Argentina, covering a total area of approximately 787,000 km 2 (Correa 1998). It was shaped during the Middle-Late Miocene (16-5.3 Ma) mainly by the uplift of the Andes. This orogeny blocked and filtered moist winds coming from the Pacific Ocean and produced a substantial rain shadow to the east, resulting in cold and dry conditions throughout the steppe, a contraction of the forest, and the expansion of xerophytic taxa (Barreda and Palazzesi 2007). Major glacial-interglacial cycles in Patagonia began during the Late Miocene (c. 5-6 Ma), with ice reaching its maximum development in the Early Pleistocene (c. 1.15 Ma; Ogg et al. 2008;Rabassa et al. 2011). The major glacial advances in southern South America were the Great Patagonian Glaciation (GPG; c. 1 Ma) and the Last Glacial Maximum (LGM; c. 21 Ka). The GPG resulted in the greatest expansion of ice sheets in the Patagonian steppe compared with LGM and other Pleistocene glaciations (Rabassa et al. 2011). However, a large part of the steppe remained free of ice during the glacial episodes, and climatic conditions were less severe than those occurring in North America (Markgraf et al. 1995).
In Patagonia, Pleistocene glaciations caused several climatic and environmental changes. The sea level lowered, shifting the coastline about four degrees eastward. Climatic continentality of the surrounding areas enhanced, with rising extreme temperatures and decreased precipitation. The Argentinean submarine platform was partially exposed, substantially increasing the space available for animal and plant colonization (Rabassa 2008). Although these climatic and environmental changes are thought to have affected species distribution, the Patagonian steppe remains poorly phylogeographically studied (S ersic et al. 2011). Plant studies have focused mainly on the Andean region (e.g., Muellner et al. 2005;Marchelli and Gallo 2006;Amico and Nickrent 2009;Acosta and Premoli 2010;Vidal-Russell et al. 2011;Segovia et al. 2012;among others), although interest in understanding the evolutionary history and processes of plant diversification in the Patagonian steppe has increased in recent years (Jakob et al. 2009;Cosacov et al. 2010Cosacov et al. , 2013Sede et al. 2012). Considering this background, subgen. Strongyloma is a relevant plant phylogeographical model not only to understand how recent past environmental changes shaped the genetic structure of its populations, but also to explore regional phylogeographical scenarios. Here, we examine the relative importance of probable historical events and/or past demographic processes that can explain the current genetic population structure in subgen. Strongyloma through phylogeographical analyses using plastid DNA sequences. We also investigate how palaeoclimatic conditions have influenced the distribution of this subgenus using distribution modelling (Hijmans and Graham 2006). We analyze three hypotheses: (1) the existence of multiple glacial refugia in flanking zones of the Andes covered by ice followed by postglacial expansion; (2) the expansion of populations over the current marine platform during the Pleistocene and their successive, recent retraction to the surroundings of the current Atlantic coast of Patagonia; and/or (3) the fragmentation of populations that persisted in areas not reached by the ice sheet during glacial periods in the center of the steppe due to past river basins dynamic. Under the first scenario, population expansion from refugia should leave genetic signatures of high diversity and private haplotypes in refugial areas, produced by considerable genome divergence and reorganization. From the perspective of the second hypothesis, in recently reduced areas nearby the Atlantic coast with rapid population retraction, diversity should be low due to successive bottlenecks on the colonizing genomes that lead to a loss of alleles. Lastly, if populations persisted in more or less fragmented areas, then genetic diversity and unique haplotypes within populations should be low: during posterior expansions, haplotypes may spread slowly and more equally (Hewitt 1996). Finally, we search for significant phylogeographical breaks consistent with those found previously for other Patagonian organisms and consider these breaks with prior hypotheses (S ersic et al. 2011;and literature therein;Sede et al. 2012;Cosacov et al. 2013).

Materials and Methods
Sampling A total of 372 individuals from 63 populations were sampled in Argentina (excepting Tierra del Fuego), covering most of the distribution of subgen. Strongyloma and its wide morphological variation, including axillaris, glomerulosa, fuegiana, and ulicina morphological variants (Table 1, Fig. 1). The maeviae morphological variant is known only from the type collection, and DNA could not be obtained from the herbarium voucher. Fresh leaves were sampled from one to eight individuals at each site, separated by at least 10 m along a linear transect to minimize the potential of sibling relationships among samples. Herbarium vouch-ers were identified following Cabrera (1982) and deposited at SI (Instituto de Bot anica Darwinion) or CORD (Museo Bot anico, Universidad Nacional de C ordoba).
DNA was amplified and sequenced with a profile consisting of 94°C for 3 min followed by 30 cycles of 94°C for 1 min, 50 or 52°C for 1 min for trnQ-5 0 rps16 and rpl32-trnL, respectively, and 72°C for 1 min. Amplification products were purified using PCR 96 cleanup plates (Millipore Corp, Billerica, MA), sequenced with BigDye 3 (Applied Biosystems, Foster City, CA), and purified with Sephadex (GE Healthcare, Piscataway, NJ) before electrophoresis on an AB 3730xl automated sequencer housed in Brigham Young University's DNA Sequencing Centre. Electropherograms were edited and assembled using SEQUENCHER 4.6 (Gene Codes, Ann Arbor, MI). Sequences were aligned with CLUSTALX 1.81 (Thompson et al. 1997) with subsequent manual adjustments using BIOEDIT 5.0.9 (Hall 1999). Indels were coded as binary characters using simple indel coding (Simmons and Ochoterena 2000) as implemented in SEQSTATE 1.4 (M€ uller 2005). The two chloroplast regions were concatenated into a single matrix a priori given their shared ancestry, maternal inheritance, and extremely rare instances of recombination within the chloroplast genome. Unique sequences from this combined matrix were identified as haplotypes. All sequences were deposited in GenBank (rpl32-trnL KM978451-KM978491 and trnQ-5 0 rps16 KM978492-KM978532).

Haplotype network and geographic distribution
Nested clade analysis (NCA) was performed following Templeton et al. (1995) and Templeton (2004) using ANECA 1.2 (Panchal 2007). Although NCA has been criticized (Petit 2007;Beaumont and Panchal 2008;Knowles 2008), it has been validated extensively (Templeton 2008(Templeton , 2009 and is used as an exploratory analysis in conjunction with other population genetic parameters to draw inferences here. A haplotype network was constructed using TCS 1.21 (Clement et al. 2000) with the default 0.95 probability connection limit. Ambiguous connections (loops) were  resolved using predictions from coalescent theory and information about sampling according to three criteria: frequency, network topology, and geography (Crandall and Templeton 1993). The resolved haplotype network was converted into a hierarchical nested design following Templeton et al. (1987) and Templeton and Sing (1993). Clade (Dc) and nested clades (Dn) distances were estimated to assess association between the nested cladogram and geographic distances among sampled localities (Templeton et al. 1995) using GEODIS 2.6 . Null distributions (i.e., under a hypothesis of no geographic association of clades and nested clades) for permutational contingency table test comparisons were generated from 10,000 Monte Carlo replications, with a 95% confidence level. For significant associations, the inference key of Templeton (2011) was used to recognize probable demographic processes and/or historical events of the clades.

Population genetic analyses
Haplotype diversity (h; Nei 1987), nucleotide diversity (p; mean number of pairwise differences per site; Nei 1987), and mean number of pairwise differences (p; Tajima 1983) were calculated for the subgenus as a whole, each sampling location, and each geographic region, using DNASP 5.10.01 (Librado and Rozas 2009).
To investigate hierarchical levels of population structure, two analyses of molecular variance (AMOVA) were performed that consider genetic distances between haplotypes and their frequencies using ARLEQUIN 3.5.1.2 (Excoffier and Lischer 2010). First, populations were grouped according to different latitudinal regions, defined by major rivers: (1) north of Agrio, Neuqu en, and Negro rivers, c. 38°S; (2) between Agrio, Neuqu en, and Negro rivers, and Chubut River, c. 38°-43°S; (3) between Chubut and Deseado rivers, c. 43°-47°S; (4) between Deseado and Chico rivers, c. 47°-50°S; and (5) south of Chico River, c. 50°S (Fig. 1). The aim of the second analysis was to test whether genetic variation was significantly structured among glaciated and nonglaciated sites during GPG in populations south of Chubut River (Rabassa et al. 2011), following a longitudinal transect at c. 71°W. This region was selected given that eight of the 63 populations were located within the ice zone ( Fig. 1). Estimation of molecular diversity indexes and AMOVAs were performed considering only populations with four or more sampled individuals.

Demographic history analyses
To analyze the demographic history of populations in each geographic area, Tajima's D (Tajima 1989), Fu's F S (Fu 1997), and Ramos-Onsins and Rozas' R 2 (Ramos-Onsins and Rozas 2002) statistics of neutrality were calculated to test for evidence of range expansion. Significant negative values of D and F S and small positive values of R 2 indicate an excess of low frequency mutations relative to expectations under the standard neutral model (i.e., strict selective neutrality of variants, constant population size, and lack of subdivision and gene flow). F S and R 2 are the most powerful tests used to detect population growth (Ramos-Onsins and Rozas 2002). Significance was evaluated by comparing observed values with null distributions generated by 10,000 replicates, using the empirical population sample size and observed number of segregating sites implemented by DNASP 5.10.01.

Distribution modelling
To validate the phylogeographical analyses, we modelled the present potential geographic distribution of Nassauvia subgen. Strongyloma and projected it onto a LGM model about 21 Ka using point locality information and environmental data in DIVA-GIS 7.3 (Hijmans et al. 2001) and MAXENT 3.3.3, with the maximum entropy machinelearning algorithm (Phillips et al. 2004).
We recorded latitudinal and longitudinal coordinates from 224 localities of subgen. Strongyloma, covering its entire distribution range primarily using a handheld geographical positioning system (GPS) unit. Only 36 point localities of the total were geo-referenced using Google Earth (http://www.google.com/earth/index.html). Environmental data with a resolution of 2.5 arc minutes (5 km 2 ) for current and past conditions were downloaded from the WorldClim database 1.4 (Hijmans et al. 2005) and were represented by 19 bioclimatic variables derived from the monthly temperature and rainfall values.
Current conditions are interpolations of observed data from climate stations around the world, representing a 50-year period from 1950 to 2000. Past conditions for the LGM are calibrated and statistically downscaled reconstructions based on the WorldClim data for current conditions. We chose the prediction of one of two global climate models tested for the LGM, the community climate system model (CCSM) and the Model for Interdisciplinary Research on Climate (MIROC), by evaluating the MAXENT values of the area under the curve (AUC) and interpreting the resulting maps of each model in DIVA-GIS.
We generated 2000 random points from across Argentina, Bolivia, and Chile using GEO Midpoint (http://www. geomidpoint.com/random/) and extracted environmental data in DIVA-GIS. To avoid overestimation of climatic data that can lead to misleading results (Phillips et al. 2004;Peterson and Nakazawa 2008), we calculated Pearson's correlation coefficients (r ≥ 0.75) and performed Principal Component Analyses in INFOSTAT 2.0 (Di Rienzo et al. 2011) following a similar procedure described by Werneck et al. (2011) and Sede et al. (2012). We identified highly correlated variables and selected ten that we considered more biologically meaningful and directly relevant to subgen. Strongyloma. These were the annual mean temperature (Bio 1), isothermality (Bio 3), temperature annual range (Bio 7), mean temperature of the wettest quarter (Bio 8), mean temperature of the driest quarter (Bio 9), mean temperature of the coldest quarter (Bio 11), annual precipitation (Bio 12), precipitation seasonality (Bio 15), precipitation of the warmest quarter (Bio 18), and precipitation of the coldest quarter (Bio 19).
The bioclimatic layers were trimmed to the surrounding areas of the known geographic distribution of subgen. Strongyloma, and then projected over a wider region from latitude 17°41 0 S to 55°49 0 S and from longitude 57°31 0 W to 75°43 0 O for the present and the LGM in DIVA-GIS.
MAXENT was run using the following settings for current and past models: 10 replicates with auto features, response curves, jackknife tests, logistic output format, random seed, random test percentage = 25% (the test data are a random sample taken from the 100% of the species presence localities, which allows the program to do some simple statistical analysis; therefore, the remaining 75% represents the training data), replicate run type = cross-validate, regularization multiplier = 1, maximum iterations = 500, convergence threshold = 0.00001, and maximum number of background points = 10,000.
To determine the threshold value for each prediction, we used the value of the minimum training presence logistic threshold. Variable importance was determined comparing percent contribution values and jackknife plots.

DNA data sets
DNA sequences of subgen. Strongyloma ranged from 913 to 928 base pairs (bp) for rpl32-trnL (with two indels of seven and eight bp) and 916 to 925 bp for trnQ-5 0 rps16 (with one indel of nine bp). The three indels were coded as present or absent and appended to the end of the matrix. The combined matrix included 1856 characters, of which 37 were variable.

Haplotype network and geographic distribution of haplotypes
Statistical parsimony retrieved a well-resolved network in which three closed loops were resolved unambiguously (Fig. 2). A total of 41 haplotypes (H) were identified. Eleven haplotypes occurred in interior positions and 30 at tip positions. Nine additional haplotypes were not observed in the analyzed individuals and occur only as inferred intermediates in the network. Four frequent, widespread haplotypes were found (H25, H1, H18, and H12; Table 2, Fig. 2).

Haplotype diversity
Considering all individuals from all 63 sampled populations, haplotype diversity (h) ranged from 0 to 0.821 with an average of 0.532, nucleotide diversity (p) ranged from 0 to 0.00218 with an average of 0.00128, and mean number of pairwise differences among haplotypes within populations (p) ranged from 0 and 3.857 with an average of 1.175 (Table 4).
With one exception, the highest values of h were found in populations of central-western Patagonia (populations 5, 9, 30, 31, 36, 37, 42, 44, and 51; Table 1, Fig. 1), located from the flanks of the Andes up to 69°W, and latitudinally from 40°S in southern Neuqu en to 50°S in southern Santa Cruz. Population 57 had the highest haplotype diversity and occurs outside this range in southeastern Santa Cruz. The highest values of p also occur in populations of central-western Patagonia, coinciding with most populations that showed higher haplotype diversity (populations 5, 9, 31, 37, 48, and 57), and in one population (24) located in south-eastern Santa Cruz. The highest values of p similarly occur in populations of central-western Patagonia, coinciding with most populations that showed higher haplotype and nucleotide diversity (populations 5, 9, 31, 37, 48, and 57). Most populations with high genetic diversity also had exclusive haplotypes (populations 5, 30, 31, 36, 37, 42, 44, and 57, Table 4).
The same diversity indices were calculated for each geographic region (Table 5, Fig. 3). The NP2 region showed the highest value of haplotype diversity (h = 0.800), followed by the CP2 region, which also

Population genetic structure
Analyses of molecular variance (AMOVA; Table 6) showed significant differences among all sources of variation tested (P < 0.001), with the exception of regions located inside and outside the boundary of the GPG (Fig. 1). The AMOVA revealed that genetic variation of subgen. Strongyloma was best explained by differences among populations within regions that were glaciated or nonglaciated during GPG (73.6%, P < 0.001). By dividing the regions according to major rivers (Fig. 1), the highest percentage of variation was explained by differences among populations within each of the regions (50.6%, P < 0.001).

Past demographic analyses
Recent demographic expansion was evidenced for the NP1 and the SP1 regions. For the NP1 region, Fu's F S and Ramos-Onsins and Rozas' R 2 statistics showed a negative value of À4.18 and a small positive value of 0.0407, respectively. For the SP1 region, Tajima's D and Fu's F S showed negative values of À1.54473 and À7.324, respectively, while Ramos-Onsins and Rozas's R 2 showed a small positive value of 0.0353 (Table 5).

Potential Pleistocene geographic distribution
Present and LGM species distribution models for subgen. Strongyloma are presented in Figure 4. In predicting the current distribution, the mean value of the AUC for the test data of the model for the present was 0.932 AE 0.008. Because the mean value of the AUC for the test data using the MIROC palaeomodel was 0.836 AE 0.075, compared with 0.789 AE 0.021 using the CCSM palaeomodel, the former was selected for projecting the distribution during the LGM.
The predicted model for the present (Fig. 4A) was a fairly good representation of the known geographic distribution for subgen. Strongyloma, although some inconsistencies between projected and realized distribution areas were found. Low probability (0.21-0.30) of suitable conditions was found in the Islas Malvinas/Falkland Islands, where these plants have never been observed. Table 3. Population processes and/or historical events affecting genetic structure in Nassauvia subgen. Strongyloma based on nested clade analysis (NCA). Geographic region, hierarchically nested clades, results of permutation contingency tests with their associated P value, NCA inference chain and inferred events are shown.
For the training data of the palaeomodel, the highest contribution and the highest permutation importance was for the annual mean temperature (Bio1 = 42.8% and 36.1%, respectively), indicating that this environmental variable contributed most to the construction of the model. In the jackknife test for the training data, the mean temperature of the coldest quarter (Bio11) presented the most useful information when used in isolation to build the model, causing the model to reach the highest gain and allowing a reasonable adjustment to the training data. The mean temperature of the wettest quarter (Bio8) contained a substantial amount of useful information that was not contained in the model when only the remaining variables were used, since excluding Bio8 in the model decreased training gain. Jackknife tests for the gain and for the AUC of the test data yielded similar results to the jackknife test for the gain of the training data, suggesting that mean winter temperature (i.e., the coldest and the wettest quarter) helped MAXENT achieve a good fit to the training data and also to generalize better, achieving comparatively better results on the test data.

Discussion
The phylogeographical and palaeomodelling analyses suggest that the genetic diversity and structure found in Nassauvia subgen. Strongyloma were shaped mainly by three main different historical events: (1) isolation in refugia and restricted gene flow during Pleistocene glacial periods, followed by postglacial range expansions; (2) colonization of new available areas during Pleistocene glaciations and postglacial retraction of populations; and (3) fragmentation of populations in areas not reached by the ice sheet during Pleistocene glacial periods followed by postglacial range expansions.

Pleistocene refugia in central-western and south-eastern Patagonia
Our data suggest that populations of Nassauvia subgen. Strongyloma may have withstood the extreme cold and aridity during the Pleistocene glaciations in two kinds of refugia (Holderegger and Thiel-Egenter 2009): (1) peripheral glacial refugia-at least three along the eastern foothills of the Andes-located in south-western Neuqu en (c. 39°30 0 S 70°30 0 W), north-western Chubut (c. 43°S 71°W), and north-western Santa Cruz (c. 47°S 71°W); and (2) lowland glacial refugia-at least four-located near the Chubut River in central-western Chubut (c. 44°S 69°30 0 W), near the Deseado River in central-northern Santa Cruz (c. 46°30 0 S 69°30 0 W), and near the Chico River in central-western (c. 48°S 71°W), and south-eastern (c. 50°S 69°W) Santa Cruz. The high diversity and several exclusive haplotypes found in these populations of subgen. Strongyloma are consistent with long-term habitation such as would be expected from Pleistocene refugia (Hewitt 1996;Taberlet and Cheddadi 2002;Petit et al. 2003;Provan and Bennett 2008). The finding of exclusive haplotypes suggests there would have been restricted gene flow and that isolation over time allowed the accumulation of genetic differences among populations. Most of these populations were found within the limits of the ice sheet of the GPG. Moreover, restricted gene flow and posterior range expansion were inferred by NCA and corroborated in most cases by demographic analyses, indicating that they might have survived in refugia during glacial periods, as suggested for other Patagonian organisms in the areas mentioned above (S ersic et al. 2011;Sede et al. 2012;Cosacov et al. 2013). The highest probabilities of occurrence projected especially for central-western and southeastern Patagonia, both for the current distribution, as well for the palaeodistribution, reinforce the conclusion that subgen. Strongyloma probably persisted in refugia in these areas, without drastic climatic restrictions.

Colonization of eastern Patagonia during Pleistocene glaciations
During Pleistocene glaciations, the current Atlantic submarine platform was partially exposed as a result of the sea level lowering (c. 100-140 m during full glacial episodes), allowing for plant colonization (Rabassa 2008). Palaeomodelling supports that populations of Nassauvia subgen. Strongyloma could have inhabited that exposed area of the coast. After glaciations, and as the result of the sea level rise, populations were forced to retract to the current Atlantic coast. Low diversity and the absence of exclusive haplotypes in eastern Patagonia are then explained by the expansion of populations over the current Atlantic coast during Pleistocene glaciations, and their recent, postglacial retraction, with the consequent loss of genetic diversity.

Population fragmentation in central Patagonia and phylogeographical breaks
Latitudinal discontinuities have been previously reported for other Patagonian organisms, both plants and animals (S ersic et al. 2011;Sede et al. 2012;Cosacov et al. 2013). Dynamics of the river basins, palaeobasins, and coastline shift during the Pleistocene glacial-interglacial cycles constitute factors that may have converted Patagonia into a highly fragmented landscape (S ersic et al. 2011). In particular, the Chubut and the Deseado river basins may have been important historical barriers that facilitated isolation between populations of different Patagonian species that persisted north and/or south of these river basins (e.g., Phyllotis xanthopygus, Kim et al. 1998;Hordeum spp., Jakob et al. 2009;Mulinum spinosum, Sede et al. 2012;Anarthrophyllum desideratum, Cosacov et al. 2013). However, adverse climatic conditions such as lower winter temperature and precipitation during glacial periods have been suggested as a possible reason for plant population fragmentation in the same geographic area (Sede et al. 2012;Cosacov et al. 2013). Our palaeomodelling results also support a low probability of suitable climatic conditions for Nassauvia subgen. Strongyloma during the LGM in the center of the steppe, between the Chubut (c. 44°S) and the Deseado (c. 47°S) rivers. In fact, the geographic distribution of this subgenus, as in the steppe shrub M. spinosum (Sede et al. 2012), seems to be influenced mainly by the mean temperature of the winter, a biologically important variable that defines the start of the growing season (Paruelo et al. 1998). Our findings also support the traditional water balance hypothesis, which considers aridization as the main factor affecting plant population demography and diversification in Patagonia (Cosacov et al. 2013), related to the enlargement and reduction of the main river basins. Based on past river basin dynamics and climate conditions cited above, a significant difference could be expected between northern and southern haplotypes within subgen. Strongyloma. Nevertheless, there is a close genetic relation between a few northern and centralnorthern lineages with southern lineages distributed in distant areas. However, a general northern-southern haploclade structure can be supported combining the above environmental evidence with NCA inferences and demographic analyses, in two historical steps: (1) both fragmentation between northern and southern haploclades, and isolation with restricted gene flow of one southern haploclade during Pleistocene glaciations, produced the split of the original distribution into two groups of population, one north of the Chubut river basin, and one south of the Deseado river basin; and (2) successive range expansions after the Pleistocene glaciations, mainly toward warmer and wetter conditions in the center of the Patagonian steppe. The most likely posterior demographic expansions would be from both northern and southern populations toward central-eastern Chubut Province. The expansion of northern, central, and southern haploclades could have generated recent contact areas throughout the entire Patagonian steppe, but mainly between the Chubut and the Deseado rivers.
Within a general latitudinal geographic structure of haploclades, other phylogeographical breaks can be distinguished related to river basins. On one hand, the high-Andean haploclade located just north of the Agrio and Neuqu en rivers at c. 38°S was clearly separated from the rest, in accordance with a latitudinal break reported for a lizard species (Liolaemus elongatus, Morando et al. 2003). Low diversity and the presence of only one exclusive haplotype in high-Andean populations may be explained by recent colonization or this was a poor area for the survival of subgen. Strongyloma during the LGM. However, these results should be interpreted with caution because of incomplete sampling of high-Andean region. On the other hand, the northern limit of one central haploclade (CP1) coincides with the break at Limay River basin reported for another lizard species (Liolaemus bibronii, Morando et al. 2007).
Finally, according to haplotype analyses, the four most widespread haplotypes and several other interior and tip haplotypes were shared among populations of different morphological variants, suggesting that there is no strict relation between morphology and genetic diversity. Furthermore, throughout the entire Patagonian steppe, some populations have haplotypes from different clades. This indication of distantly related haplotypes co-occurring suggests the existence of numerous contact zones (Avise 2000) throughout the steppe: within northern lineages in central R ıo Negro; between northern and central clades in central-western R ıo Negro; between northern and southern lineages in southern Neuqu en, north-western Chubut, and north-western Santa Cruz; between central and southern clades in central-northern Santa Cruz; and within southern lineages in central, south-western, and south-eastern Santa Cruz.

Conclusion
Phylogeographical and palaeomodelling patterns found in Nassauvia subgen. Strongyloma appear to reflect climatic fluctuations and landscape changes that occurred during the Pleistocene having an impact on the distribution, demography, and diversification of different lineages within the group. Past events differentially affected areas within the immense geographic range of these shrubs. During glaciations, subgen. Strongyloma persisted in peripheral refugia along the eastern foothills of the Andes, in lowland refugia in central-western and south-eastern Patagonia, over the exposed Atlantic submarine platform toward the east, and in fragmented populations north and south the Chubut and the Deseado river basins. Related to those past scenarios, Pleistocene glaciations indirectly affected the demography of this group of plants through decreased winter temperatures and water availability in different areas of its immense geographic range. Future comparative phylogeographical studies should consider patterns found in subgen. Strongyloma to better comprehend the processes that have affected the Patagonian biota.