Evidences of delayed size recovery in Araucaria angustifolia populations after post-glacial colonization of highlands in Southeastern Brazil

Up to date, little is known about the relationship between historical demography and the current genetic structure of A. angustifolia. As a first effort towards overcoming this lack, microsatellite data scored in six populations and isozyme allele frequencies published for 11 natural stands of this species were analysed in order to assess molecular signatures of populations’ demographic history. Signatures of genetic bottlenecks were captured in all analysed populations of southeastern Brazil. Among southern populations, signatures of small effective population size were observed in only three out of 13 populations. Southern populations likely experienced faster recovery of population size after migration onto highlands. Accordingly, current genetic diversity of the southern populations gives evidence of fast population size recovery. In general, demographic history of A. angustifolia matches climatic dynamics of southern and southeastern Brazil during the Pleistocene and Holocene. Palynological records and reconstruction of the past climatic dynamics of southeastern and southern Brazil support the hypothesis of different population size recovery dynamics for populations from these regions.


INTRODUCTION
Although the lands in Brazil have never been covered by ice sheets, the remarkable climatic changes during the Pleistocene and Holocene have influenced distribution and evolution of different species.For instance, due to the cold climate with long dry periods during the last glaciation (Ledru et al. 1998), the subtropical highlands of Brazil lacked forest formations and were covered by grassland throughout the Late Pleistocene, about 48000 to 18000 14 C yr BP (uncalibrated radiocarbon years before present).During this adverse period, dispersed stands of Araucaria angustifolia (Bert.)relatively recently, signatures of bottleneck events may be retained in its genome and are likely to be assessed through genetic analyses.Additionally, if such a reduction in effective population size occurred as effect of post-glacial migration from refugia onto highlands, it is expected that climate dynamics match the genetic demographic history of populations.

MATERIALS AND METHODS
Multilocus microsatellite genotypes were scored at five loci [CRCAc2 (Scott et al. 2003), Ag20, Ag45, Ag94 (Salgueiro et al. 2005) and AA01 (Schmidt et al. 2007)] as described by Stefenon et al. (2007).Isozyme allelic frequencies were obtained from their original publications (Auler et al. 2002, Sousa et al. 2004, Mantovani et al. 2006).Founded on the principle that rare alleles are lost in populations which experienced a strong reduction in effective size (Nei et al. 1975), three methods were employed to search for signature of population bottlenecks.The first method consists in testing for heterozygosity excess in comparison to predictable heterozygosity under mutation-drift equilibrium considering the observed number of alleles (Cornuet and Luikart 1996).Heterozygosity excess should not be confused with the excess of heterozygotes (comparison of observed and expected heterozygoties).Bottlenecked populations are expected to reveal a significant excess of heterozygosity than estimated under equilibrium, because rare alleles are lost faster than the heterozygosity.Heterozygosity excess was tested by comparing the expected heterozygosity under mutation-drift equilibrium (H E Q ) with the Hardy-Weinberg heterozygosity (H e ).Estimates of H E Q were obtained by simulating the coalescent process (10000 iterations) for each population following the two-phase mutation model (TPM;Di Rienzo et al. 1994), which assumes that most mutational changes result in an increase or decrease of one repeat unit but also incorporates mutations of larger size.The isozyme data were investigated using the infinite allele model (IAM; Kimura and Crow 1964), which assumes that each arising allele is unique.Statistical significance of the analyses was assessed using the Wilcoxon test (Luikart et al. 1998a).
The second approach is based on the assumption that populations that experienced a recent reduction in   Sousa et al. 2004, c Mantovani et al. 2006.effective size tend to show a distortion of allele frequency distribution (Luikart et al. 1998b).For this mode-shift analysis, alleles were grouped into 10 allele frequency classes and plotted in a frequency histogram.Bottlenecked populations have a propensity to display a shifted distribution, with the incidence of alleles at frequency lower than 0.1 becoming lower than the incidence of alleles in an intermediate allele frequency class.
The third method consists in analysing the ratio of the total number of alleles (k) to the overall range in al-An Acad Bras Cienc (2008) 80 (3) lele size (r ) as M = k/r (Garza and Williamson 2001).For this analysis, alleles of the microsatellite locus Ag20 were binned into classes to fit the step-wise model.To test whether M-values are lower than expected (cases in which a population experienced a bottleneck), 10000 replicates of the coalescent process were simulated assuming a proportion of one-step mutations ( p s ) of 0.2, a mean size of non one-step mutations ( g ) of 3.5 and two different bottleneck population sizes, assuming a mutation rate (μ) of 10 -4 mutations per locus per generation: N e = 50(θ = 0.02) and N e = 500(θ = 0.2), where θ = 4N e μ.Analyses were performed using the programs BOTTLENECK 1.2.02 (Piry et al. 1999) and M_P_VAL (Garza and Williamson 2001).
Further evidences of bottleneck events were investigated assuming that genetic distance values increase rapidly when a bottleneck occurs and, therefore bottlenecked populations tend to present elongated branches in cluster analyses (Takezaki and Nei 1996).Phylograms were constructed for microsatellite data with the program POPULATIONS 2.0 (Langella 2002) using the (δμ) 2 genetic distance (Goldstein et al. 1995) and the neighbour-joining clustering algorithm.

RESULTS
The accentuated genetic differentiation reported among populations from southeastern and southern regions (e.g.Sousa et al. 2004, Stefenon et al. 2007) has been explained as effect of post-glacial migration from different refugia, which in addition to the geographical isolation, resulted in temporal isolation of the southeastern populations (Stefenon et al. 2007).Four populations analysed in this study belong to the southeastern region (CJ, CJm, CJs1 and CJs2 from São Paulo State; Fig. 1) and revealed signatures of bottleneck in at least one of the analyses.Among the southern populations (from Paraná, Santa Catarina and Rio Grande do Sul States), molecular signature of low population size was detected just in three out of 13 stands.
In the analysis of heterozygosity excess (Tables I  and II) signature of bottleneck was detected only for population CJm, with excess of heterozygosity observed in six out of seven loci (Table II).Three out of these six loci revealed significant excess of heterozygosity at 5%  I and II.
level.The Wilcoxon test for the multilocus analysis was highly significant (P = 0.008).
In the mode-shift analysis of isozyme data, populations CJm, CJs1, CJs2 and FTB revealed a shifted distribution of alleles (Fig. 2A), characteristic of bottlenecked populations.All populations analysed with microsatellite markers revealed a normal distribution of alleles in this analysis (Fig. 2B), suggesting absence of bottleneck signature.
In the M-ratio analysis, populations BJ, RG and CJ revealed signatures of bottleneck (Table I) when considering an effective population size of 50 individuals for the post-migration population (θ = 0.02).When θ was set to 0.2 in the M-ratio analysis (N e = 500; Table I), signature of low effective population size was retained just in population CJ, emphasizing the occurrence of a bottleneck in this stand.Additional evidence of a bottleneck event in population CJ is given by the elongated branch revealed by this population in the neighbour-joining phylogram (Fig. 2C).Although displaying a relatively low bootstrap support (55%), such elongated branches are usually observed for bottlenecked populations (Takezaki and Nei 1996).

DISCUSSION
In this study, the occurrence of genetic bottlenecks in A. angustifolia populations was tested with complementary methods, providing a more comprehensive view of the demographic events related with historical low effective population size.The graphical method and the heterozygosity excess analysis are more powerful in detecting molecular signature when pre-bottleneck θ is small, the bottleneck was severe and the population recovered sample size quickly, while the M-ratio analysis is more efficient in opposite circumstances (Busch et al. 2007).Moreover, the heterozygosity excess and the shift-mode analyses are more powerful in capturing signatures of bottleneck from loci evolving under the IAM model (Cornuet and Luikart 1996), while the M-ratio was specifically designed to the analysis of allele distributions of loci evolving under the SMM and TPM models (Garza and Williamson 2001).The population size sampled for microsatellite analysis is also an important factor, given that (due to its high polymorphism) a higher number of alleles is expected to be observed if a larger sample size is investigate, fact that is not observed in isozyme analysis (see the similar number of observed alleles in populations with 10-fold difference in population size in Table II).The number of sampled individuals for the microsatellite analysis was three times the number of alleles of the most variable locus (AA01 with 20 alleles), ensuring that the majority of the alleles in the population were sampled (Garza and Williamson 2001).Concerning the timeframe since the occurrence of the bottleneck, the methods cover different and complementary time windows: 0.2N e to 4N e generations for the heterozygosity excess method (Cornuet and Luikart 1996), 2N e to 4N e generations for the model shift method (Luikart et al. 1998b) and a comparative larger timeframe for the M-ratio analysis (Garza and Williamson 2001).Assuming a very low population size in the colonization of highlands by A. angustifolia may be unrealistic, given the high amount of pollen and seeds produced (one single tree can produce more than 1000 seeds per year; Mantovani et al. 2004).Thus, N e = 50 individuals should be a reasonable approximation of the effective population size of the post-glacial founder populations of A. angustifolia.Assuming this post-bottleneck effective population size, the timeframe ranges from 10 to 200 generations for the heterozygosity excess analysis and from 100 to 200 generations for the shift-mode analysis.For the M-ratio analysis, this timeframe goes up to more than 200 generations.We intended to detect bottleneck events occurred as effect of the post-glacial migration of the species.Assuming a generation interval ranging from 15 to 30 years for A. angustifolia, the methods applied will detect signatures of bottleneck events occurred between 150 and more than 6000 years ago, covering the period of the migration onto highlands (1000 to 3500 14 C yr BP).

GENETIC DYNAMICS OF BOTTLENECKED POPULATIONS
After a bottleneck, populations enter a recovery phase and new alleles are created by mutation.This population expansion and arise of new alleles generates a heterozygosity deficiency which can erase the signature of population decline (Cornuet and Luikart 1996).Consequently, signature of population reduction will be lost if population size recovery is fast.Recent forest exploitation should not be a source of bias in comparing A. angustifolia populations from southeastern and southern Brazilian regions, given that both share the same recent events (comprising the last 100 years).Assuming that the southeastern region was colonized earlier than the southern region, the preserved signature of a bottleneck suggests that southeastern stands recovered slowly after post-glacial migration onto highlands.On the other hand, southern populations may have displayed somewhat quicker effective size recovery, even if founded with low population size.Additionally, migration events among the southern stands might have overturned the effect of a bottleneck event on the allelic diversity of these populations.In small populations, the movement of just a few migrants can erase the signature of bottleneck in two to three generations (Busch et al. 2007).
Based on analyses of allelic structure of microsatellite loci, molecular signatures of genetic bottleneck were detected in an isolated stand but not in a large population of Pinus taeda in central Texas/USA (Al-Rabab'ah and Williams 2004).The occurrence of long dry periods during the Holocene was the reason postulated for the persistent low population size of the isolated population.Similar to P. taeda, prolonged drought is a likely factor to have intensified the effect of reduced effective size of A. angustifolia populations in southeastern Brazil.Ledru et al. (1998) propose that the modern climate in central Brazil was reached just about 2500 14 C yr BP.If A. angustifolia colonized the southeastern Brazil about 3000 14 C yr BP as suggested by the reconstruction of vegetation and polar advections trajectory (Ledru et al. 1998) and by the palynological record (Behling 1998(Behling , 2002)), relatively dry periods followed the post-glacial migration for at least 500 years (these dating should be considered with caution due to imprecision of the radiocarbon method) and may have delayed population recovery in this region.At present day, the region of Campos do Jordão, southeastern Brazil, displays a low rainfall period of four months between May and August (Behling 1997b).
A shifted frequency distribution of alleles at five isozyme loci in the endemic P. maximartinezii was also interpreted as effect of an extreme bottleneck (Ledig et al. 1999).The authors evoked a rapid post-bottleneck expansion as the reasoning for the diminished effects of genetic drift observed.Similarly, rapid population expansion from a bottleneck event was suggested to explain the unimodal distribution of pairwise differences among individuals (chloroplast microsatellites analysis; Echt et al. 1998), the small estimated effective population size and the high selfing rate (nuclear microsatellites analysis; Boys et al. 2005) assessed in the widely distributed P. resinosa.
In the present study the mean number of alleles per locus observed was similar among bottlenecked and non-bottlenecked populations for microsatellites (44.7 and 44.3 respectively; p > 0.05; t-test = 0.61; d.f.= 4).Among the populations analysed with microsatel-An Acad Bras Cienc ( 2008) 80 (3) lite markers, population CJ revealed the lowest number of alleles and gene diversity (H e ).On the other hand, population RG revealed the highest number of alleles and gene diversity among the six investigated populations, although signature of bottleneck was captured in this stand in the M-ratio analysis with N e = 50(θ = 0.02).The high diversity observed in this population supports a quick effective size recovery in the southern region.As for microsatellites, the mean number of alleles per locus at isozyme markers was statistically not different between bottlenecked and non-bottlenecked populations (2.2 and 2.3 respectively; p > 0.05; t-test = 0.59; d.f.= 9).Comparing the eight populations analysed by Auler et al. (2002), population FTB revealed the third highest gene diversity (H e ) and mean number of alleles per locus.Considering the signature of bottleneck captured in this population in the mode-shift analysis, this comparatively high diversity is also evidence of a quick effective size recovery of the southern populations.
The general assumptions of the analyses employed in this study are based on multilocus investigations and on the mutation-drift model, assuming that a bottleneck event generated the observed genetic signatures (i.e.loss of heterozygosity and alleles and change in allele frequencies).Selection could also quickly change allele frequencies and generate similar loss of alleles (strong balancing selection or heterozygosity advantage could maintain alleles at intermediate frequencies and cause populations to appear to have been recently bottlenecked).However, it is unlikely that selection would reduce diversity simultaneously at many independent loci (Luikart et al. 1998a, b).For instance, no evidence of mode-shift distortion was observed in multilocus analysis of eight non-bottlenecked populations reported to have evidences of balancing selection at isozyme loci (Luikart et al. 1998b).

MATCHING OF PAST CLIMATIC DYNAMICS AND SIGNATURE OF BOTTLENECK EVENTS
If the observed signatures of past population demography is an effect of the post-glacial migration of A. angustifolia from refugia onto highlands, it is expected that reconstructed climate dynamics match the genetic inferences.The presence of small populations of A. angustifolia in low elevated areas along rivers in southeastern Brazil and in deep protected valleys and/or on the Atlantic facing slops at lower elevations in southern Brazil during the glacial period is indicated by pollen analytical studies (Behling and Lichte 1997, Behling et al. 2002, 2004, 2007).Expansion of A. angustifolia from refugia as gallery forest at low elevations onto the higher mountains in southeastern Brazil (e.g.Campos do Jordão) started already during the late glacial period (Behling 1997b), while in southern Brazil a significant expansion onto highlands occurred when climate conditions were more suitable, about 3500 14 C yr BP by the expansion of gallery forests and somewhat latter since about 1000 14 C yr BP into the grassland (Behling 1998, 2002, Behling and Pillar 2007).
Based on the pollen record from four peat bogs in Santa Catarina, Behling (1995) suggested that A. angustifolia might have persisted in protected highland valleys during the Late Pleistocene (> 10000 years ago).Just minor expansion as gallery forest along rivers occurred at the end of this period.The first major expansion of A. angustifolia onto highlands in Santa Catarina occurred as a result of a very moist climate, around 1000 14 C yr BP (Behling 1995).Palynological data from Serra dos Campos Gerais in Paraná State (Behling 1997a) revealed a Late Quaternary vegetation and climate history very similar to that of Santa Catarina with a long dry season until the beginning of the Holocene.The expansion of A. angustifolia onto highlands in the Late Holocene is evidence of a shorter annual dry period around 2850 14 C yr BP.The first broad expansion of A. angustifolia onto highlands in Paraná State occurred about 1500 14 C yr BP (Behling 1997a).Past environmental changes reconstructed on basis of the pollen record from Cambará do Sul, in the Rio Grande do Sul State (Behling et al. 2004) corroborated the results from Santa Catarina and Paraná.A. angustifolia was likely found just in deep protected valleys and/or wetter coastal slopes.Replacement of grassland vegetation by A. angustifolia started around 1100 14 C yr BP, reflecting the arrival of a wetter period without marked annual dry season.The pollen record from the southeastern highlands is more restricted, but a rather comprehensive reconstruction of the past climatic dynamics is available by Behling (1997b).The late Quaternary period (from about 35000 to 17000 14 C yr BP) was represented by a cold and dry climate in this region.Behling (2002) and Behling et al. (2004).
A comparatively warmer but still dry climate followed until around 2600 14 C yr BP, forcing A. angustifolia to continue in moist refugia.After this period, a cool and moist climate allowed the expansion of A. angustifolia from refugia onto São Paulo highlands.
A comparative analysis of palynological records from southeastern and southern Brazilian highlands (Fig. 3) corroborates the occurrence of an expressive expansion of A. angustifolia populations in the southern States starting about 1500 to 1000 14 C yr BP (Cambará do Sul, Rio do Rastro and Campos Gerais in Fig. 3), but a lack of such a major expansion in the southeastern region (São Paulo State; Itapeva in Fig. 3).This fact matches the molecular evidences of a post-glacial bottleneck in Campos do Jordão region followed by a slow recovery of effective population size, while southern populations experienced a fast size recovery after the post-glacial migration.

CONCLUSIONS AND OUTLOOK
Despite the evidences of genetic bottlenecks revealed in this study, Busch et al. (2007) demonstrated that some biological features may lead to violations of the assump-tions made for each of the methods employed here, obscuring the molecular signature of genetic bottlenecks.In addition, a low number of markers (both isozymes and microsatellites) were analysed.Therefore, corroboration of the demographic patterns observed should be obtained through genome sequencing which allows the use of more powerful tests for departure from mutationdrift equilibrium.
The action of indigenous people modifying the environment (e.g.cleaning the forest understory in the plantation of crops or in the excavation of pit houses; Bitencourt and Krauspenhar 2006) and transporting seeds during migrations likely was an important factor in shapping the current occurrence area and genetic structure of A. angustifolia in Brazil.The peak expansion of Araucaria forest during the Holocene coincides with the period of highlands occupation by indigenous groups of Taquara/Itararé Tradition (Bitencourt andKrauspenhar 2006, Iriarte andBehling 2007).Because of this close relationship between A. angustifolia and the indigenous people, the anthropogenic transport of seeds may not be ignored as contributing to the establishment of new A. angustifolia populations.Such new populations could be established with a small number of seeds and generate a kind of bottlenecked population, resulting in similar genetic signatures.
Here, signature of low effective size during populations' establishment and subsequent generations was assessed.However, consequences of recent population reduction and fragmentation in patterns of reproduction and species adaptedness will be observed just in future generations.Considering the increasing concern in conservation of A. angustifolia genetic resources, highlighting the species demographic history may aid to foresee and minimize unwanted events related to decreasing demographic and genetic population size.For instance, since isolation likely delays genetic recovery from a small population size (as suggested by the molecular signature of the southeastern populations), promoting connectivity among fragments may be a fundamental issue in scheduling genetic conservation of the extant remnants of A. angustifolia.

ACKNOWLEDGMENTS
We would like to thank Klabin S/A, Acauan Family and the staff of Parque Estadual de Campos do Jordão for supplying plant material for the study and Dr. L. Leinemann and Dr. VA Sousa for helpful comments in a previous version of this paper.VM Stefenon is supported with a grant from Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Brasília, Brazil.The comments of two anonymous referees aided the improvement of this paper.

Fig. 1 -
Fig. 1 -Map showing the natural distribution range of A. angustifolia (grey areas) in the States of Rio Grande do Sul, Santa Catarina, Paraná (southern region) and São Paulo (southeastern region).The State to which each population belongs is indicated by the arrows.Geographic coordinates of individual populations are given in Tables I and II.

Fig. 2 -
Fig. 2 -Molecular signatures of bottleneck events in A. angustifolia populations.(A) Plotting of the frequency distribution of allele classes for isozymes.Populations CJm, CJs1, CJs2 and FTB revealed a shifted distribution, indicating historical reduction of the effective population size.(B) Plotting of the frequency distribution of allele classes for microsatellites.All populations revealed a non-shifted alleles distribution, meaning absence of bottleneck signatures.(C) Neighbour-joining phylogram computed for microsatellite data based on (δμ) 2 genetic distance.Evidence of a bottleneck event in population CJ is given by the elongated branch length.Numbers at nodes in the phylogram are bootstrap values after 1000 replicates.

TABLE I Geographic location of populations and summary of the coalescent-based analyses of microsatellite data. Statistically significant values are highlighted in bold.
† RS: Rio Grande do Sul; SC: Santa Catarina; PR: Paraná; SP: São Paulo.§ Number of loci revealing heterozygosity excess.Number of loci with significant heterozygosity excess is given within parenthesis.