Inferring possible population divergence in Espeletia pycnophylla ( Asteraceae ) through morphometric and paleogeographic approaches

The phenotypic structure within and between plant populations is generally influenced by their distribution patterns in space and time; therefore, the study of their divergence is a central issue for the understanding of their microevolutive processes. We boarded the hypothesis that three populations of Espeletia pycnophylla show phenotypic divergence as one of the possible implications of their geographic isolation in the Southern Colombian Andes. We used the Elliptic Fourier Descriptors (leaf shape) and traditional leaf morphometry (leaf size) of 347 leaves to measure inter and intra-population variation and a comparison between a paleogeographic reconstruction with an actual estimate of the distribution areas of E. pycnophylla in order to identify their main changes during the last 14 000 years. The three populations showed significant differences in leaf morphometry and a positive correlation between the matrices of morphometric and geographic dissimilarities, indicating that the inter-population divergence increases between further populations, so that the morphometric structure reflects their spatial distribution. The geographical and paleogeographical estimates evidenced a conspicuous process of reduction and fragmentation of the distribution area of E. pycnophylla since the Late-Glacial until the Holocene. We suggest that these results support possible scenarios of vicariance events, which allow us to approach the divergence of these populations in terms of their historic biogeographic relations. However, genetic analyses are still needed to support these results. Rev. Biol. Trop. 58 (4): 1261-1270. Epub 2010 December 01.

According to Cuatrecasas (1986), E. pycnophylla is the last speciation event of the Espeletia complex, a group of plants which have undergone high diversification rates in the Paramos; high altitude tropical ecosystems in the northern Andes of south America (Monasterio & Sarmiento 1991).The Paramos are well known for their status of biogeographic islands (Carlquist 1974) and their high sensibility to the climate fluctuations during the Andes history (van der Hammen & Cleef 1986), which have been important concomitant factors influencing the Espeletia diversification (Rauscher 2000).However, the paramos have not been always biogeographically isolated; instead, they underwent periodic changes of geographic contraction and expansion during the glacial-interglacial cycles of the Pleistocene (van der Hammen 1974).This possibly influenced the evolution of the group.
The geographic distribution area of E. pycnophylla is actually very fragmented and there are many isolated populations on the Southern Andes of Colombia (Cuatrecasas 1986, Solarte et al. 2007), one in the northern Ecuador and a small population in the central Ecuador (Jorgensen & Leon-Yanes 1999).This isolation is a possible response of the fragmentation of the Paramo continuity during the last 14 000 years.This kind of phenomena induces such inter-population divergence in plant populations due to restricted gene flow, contrasting selecting pressures and derive, the first steps in the microevolutive diversification of plants (Loveless & Hamrick 1984).
To measure this type of divergence some approaches have arisen, such as morphology, citotaxonomy and molecular markers, but other kind of morphometric techniques such as the Elliptic Fourier Analysis of leaf shape has also prove to be a very effective phenotypic based alternative to indirectly estimate genetic variation among and within plant populations, allowing many applications in the study of population evolution (Lönn & Prentice 1995, Olsson et al. 2000, White et al. 1998).This technique accounts for the variation in the shape per se of biological structures (Rohlf & Archie 1984), eliminates the effects of size and recovers a high amount of hereditary variation (Molau & Prentice 1992, Jensen et al. 2002, Rumpunen & Bartish 2002).
Our goal was to test the hypothesis that three geographically isolated populations of E. pycnophylla shows phenotypic divergence correlated with their inter-population spatial distances, as a possible result of their geographic fragmentation, by combining: 1) the variation of leaf shape using the Elliptic Fourier Analysis; 2) a paleogeographic reconstruction of the distribution area of E. pycnophylla in the Southern Colombian Andes during the Late-Glacial, using the model of van der Hammen (1974) and 3) an actual estimate of the distribution area of this species in the same locations.We chose to reconstruct the Late-Glacial distribution of E. pycnophylla because according to paleoecological estimations, it was the time when the Paramo vegetation reached the greatest distribution area and was much more continuous than today.

MATERIALS AND METHODS
The species: E. pycnophylla is a long lived giant rosette endemic from the Paramos of Southern Colombia and Northern and Central Ecuador.This species constitute one of the most conspicuous vegetation components of these high altitude ecosystems (Solarte et al. 2007).The main diagnostic characteristics are its sessile and highly pubescent leaves with lanceolate to elliptic shape (Cuatrecasas 1986).The heights of some adult individuals may reach 5m and the stalk diameters range from 6 to 8cm (Laegard 1992).The species are mainly reproduced by a highly auto-incompatible cross pollination (Gonzales et al. 2009) and its seeds lack from pappus, what implies a great restriction for the dispersion over large distances.Moreover, there are no reports of ornitochory, the seeds fall to the ground when the inflorescences are dry enough after the fructification stage, and the seedlings grow up near to the mother plant (Cuatrecasas 1986), producing a clumped pattern of distribution (Benavides et al 2008).
According to Cuatrecasas (1986), E. pycnophylla shows some geographic variation classified in three infraspecific levels: E. pycnophylla var.galerana, E. pycnophylla spp.angelensis and E. pycnophylla spp.llanganatensis (on the central cordillera of Ecuador) but no further evidence has confirmed these groups.Moreover, there are still many isolated populations that remain unstudied and some of them have also such internal patchiness.Here, we considered as a population a group of individuals with a delimited distribution area, spatially isolated from others and clearly identifiable via cartography, according to the classification of Paramo areas in Solarte et al. (2007).
Populations and sampling: We selected three geographically isolated populations of E. pycnophylla in the Southern Colombian Andes, each one being characterized by contrasting climate regimes and geomorphologic characteristics (Solarte et al. 2007).In addition to the spatial distances, the Guaitara river canyon (450m of maximum depth) (Fig. 4) act as a geographical barrier that strengthen the isolation by distance.Within each population we established four sampling areas in each of four orientations (N, S, W, E) and three elevations, which were proportional to the extension of the whole population.Within each sampling area we randomly selected 20 adult individuals (70-100cm tall) and from each individual we took two completely developed and healthy leaves with no damage like hervivory or necrosis (Fig 1).This sampling methodology led us a sample size of 160 leaves per population.Table 1 show some geographic features of each studied population and the definitive sample sizes used for the analyses.

Morphometric analyses:
All the collected leaves were pressed and imaged on their abaxial surface with a Cannon ® Power-Shot ® digital camera at 7.1 megapixels.The images were taken at black and white with standardized conditions (distance, position and angle) upon a white background to increase the image contrast for an optimal acquirement of outlines.We used the software TpsDig to digitalize the outlines and transform those to x,y coordinates.Based on these coordinates, we calculated the Kulh & Giardina's (1982) Elliptical Fourier Descriptors (EFD's) of each leaf outline with the program EFAWin (Isaev 1995).Ten harmonics were enough satisfactory to successfully describe the shape of all leaves.A normalization of the EFD's was applied to maintain invariant the size, location, rotation and starting point of each image, and thereby, recover only the shape per se of the structures (Rohlf & Archie 1984).Additionally, three linear morphometry (LM) variables (total leaf length, maximum leaf width and length-wide proportion) were measured to compare the effectiveness of leaf shape versus leaf size to discriminate between populations.

Statistical analyses:
Within and between groups variance was accounted for understanding the morphometric structuring of populations.Both the normalized EFD's and the LM of leaves were used in canonical variates analysis (CVA) with the software SPSS 15.0, excepting the variable length-width proportion, which was analyzed with a Kruskall-Wallis test due to its departures from normality.
Since every harmonic is described by two sine and cosine components, yielding to four coefficients per harmonic (Kulh & Giardina 1982), 40 normalized EFD's were the variables

LOP: Las Ovejas Population
The geographic distances were measured on a map as straight lines between the centers of the sampling areas of each pair of populations.
used for shape analysis.Wilk's Lambda values (Λ) (being 1 and 0 the minimum and maximum discriminant powers respectively) and the canonical correlation coefficients (being 1 and 0 the maximum and minimum discriminant powers respectively) were used as indicators of the discriminating power of each canonical variate (CV).Post hoc tests of Bonferroni with Hotelling probabilities were applied to inquire differences between population pairs when the CvA showed overall significant discriminations among groups.
Mantel tests (Mantel 1967) were preformed with the software Past 1.94b to prove if the dissimilarity matrices of geographic and morphometric distances among the three populations were linearly correlated.To build these matrices, the inter-population dissimilarities for both the normalized EFD's and the LM variables were calculated as euclidean distances between the centroids of each pair of populations using Past 1.94b.The chosen alpha for all the analyses was 0.05.

Geographical and paleogeographical analysis:
The reconstruction of the distribution area of E. pycnophylla for the Late-Glacial of Pleistocene (14 000 years ago) was made through a cartographic analysis upon the van der Hammen's model of historic altitudinal shifts of Paramo ecosystems in the Colombian Andes (van der Hammen 1974).Based on paleoecological data from fossil pollen grains obtained from different sediment cores, the van der Hammen's model explains how the Paramo ecosystems in Colombia reached their maximum geographic extension in the Late-Glacial, covering an altitudinal range from 2 000 ± 100m to 3 800 ± 100m.After that, the increasing global temperatures that characterized the start of the Holocene produced the uplift of the Paramo to higher elevations.This has been happening as a gradual process until nowadays, when the Paramo sensu stricto covers an average altitudinal range from 3500±100m to 4500±100m.
We obtained the level curves of the Southern Colombian Andes at every 100m and applied the altimetry estimates of Paramo limits suggested by the van der Hammen's model for the Late-Glacial period (inferior limit: 2 000±100m; superior limit: 3 800±100m).We used the SRTM (Shuttle Radar Topography Mission) digital elevation model generated by NASA (National Aeronautics and Spaces Administration) and NGA (National Geospatial Intelligence Agency) at a resolution of 30m.The original image was initially treated via Principal Component Analysis in order to recover the most of the possible information and to facilitate its interpretation in terms of the original variability (Chuvieco 1996).For estimating the current distribution areas of E. pycnophylla we used vegetation coverage images of the Southern Colombian Andes, Landstat TM, October 14 2002 (path 9 row 59, 8 bands, BSQ format, II quadrat, with selection of training fields on screen), which allowed  us to cartographically identify the Espeletia pycnophylla-Calamagrostis effusa and Espeletia pycnophylla-Blechnum loxense associations.Field work made by Solarte et al (2007) demonstrates that E. pycnophylla only occurs in these two associations, which are clearly identifiable via Landstat satellite imagery.Finally, we analyzed the possible relations between the spatial characteristics of the distribution areas by obtaining the respective thematic maps trough a GIS in ArcGis 9.2 platform for both the past and present periods.

RESULTS
Linear Morphometry: Significant differences were found in leaf size among populations.The CvA generated two Cv's (Table 2); the first explained the 99.03% of the original data variance and had a positive perfect correlation with the Total Leaf Length (r=1).The Wilk´s Lambda (Λ) and the canonical correlation values for this function showed highly significant discrimination among populations CCvP from LOP-GvP (p<0.001).The second Cv correlated positively with the Maximum Leaf Width (r=0.98), but this variable didn't show significant differences among populations.Fig. 2 shows these results graphically.Likewise, the Length-Width Proportion of leaves were statistically significant to separate populations CCvP from LOP-GvP (p<0.001).These results revealed that according to leaf size, populations GvP and LOP are a relatively homogeneous group with long leaves and a high Length-Width Proportion, whilst CCvP is a contrasting population with short leaves and low Length-Width Proportion (Table 3).

Leaf shape morphometry with EFD's:
Significant differences were found in leaf shape among the three populations.The CvA generated two significant Cv's (Table 2), which together explained the 100% of the original EFD's variance.The first Cv recovered the  76.2% of the variation and significantly differentiated the three populations (p<0.001).However, this function showed that the CCvP population differs to a greater extent from the rest.The second function recovered the 23.8% of the total variation and mainly differentiated LOP from GvP. Fig. 3 shows these results graphically.
Relation between geographical distances and morphometric differences: High and significant positive linear correlations were found between geographic and leaf size (LM) (r=0.90, p<0.05) and geographic leaf shape (EFD's) (r=0.99)morphometric matrices, according to the Mantel correlation test.This showed that the more geographical distance between any pair of populations, the more morphometric divergence between them.However, this correlation was stronger with the EFD's.This last proves that leaf shape variation among populations is more linearly sensible to the geographic distances than leaf size.

Geographic and paleogeographic analyses:
The map of the paleogeographic reconstruction of the distribution area of E. pycnophylla 14.000 years ago (Fig. 4A) showed a large and uninterrupted central population with some small marginal populations at the south west and north of the distribution, covering a total estimated area of 8 351.4km 2 .This ancient population was first split into two parts; the CCvP and a population formed by the conjunction of the current GvP and LOP areas.After that, the fragmentation of the last two populations occurred.The map of the current distribution (Fig 4B ) shows a high degree of spatial disjunction, which covers a total estimated area of 669.90km 2 .These results demonstrates that the 92% of the geographic distribution area of E. pycnophylla has not only fragmented, but conspicuously reduced since the Late-Glacial through all the Holocene, resulting in a set of geographically isolated populations.

DISCUSSION
Variation in leaf shape and size: We think that the results of this research support the hypothesis of an allopatric population divergence of E. pycnophylla in the Paramos of the Southern Colombian Andes.The EFD's of leaf shape had greater statistical effectiveness and larger resolution than LM of leaf size in distinguishing inter-population differences.This is in agreement with a general pattern shown by similar studies (Lönn & Prentice 1995, McLellan & Endler 1998, Jensen et al. 2002, Rumpunen & Bartish 2002).However, the results of each individual method deserve a particular explanation.The differences in leaf  size could be related to plasticity or to different environmental pressures influencing each population.The significantly smaller leaf size in population CCvP may be a possible adaptation to the sandy soils, lower precipitation rates and lower water availability that govern its habitat (Solarte et al 2007).This is a typical response of plants to reduce the potential area for evapotranspiration and tissue dehydration.
On the other side, although the differences in leaf shape shown by the EFD's cannot be attributed to conspicuous biologic features or adaptations, making difficult its interpretation, Rohlf & Archie (1984) states that it indeed allows detecting differences in the overall shape per se of biological structures, which usually contain a hereditary control (Dickinson et al. 1987).In fact, some researchers have proved that whilst the LM accounts for a high percentage of environmental and developmental variation and a small percentage of genetic variation, the EFD's applied to leaf shapes do exactly the opposite, serving as very accurate estimators of genetic variation within and between plant populations.For example, Molau & Prentice (1992) and Lönn & Prentice (1995) found that for three Saxifraga species and Hippocrepis emerus respectively, the EFD's only accounted for a very small percentage of variation explained by environmental plasticity.However, the measurement of genetic variation trough molecular markers is still needed to prove if this applies to E. pycnophylla.
Identifying the mechanisms behind the inter-population differences in leaf shape is beyond the scope of this work because we still need to apply genetic measurements and other tools to solve this question.However, we suggest as a hypothesis that the restrictions in seed dispersion and pollination over long distances among isolated populations could be responsible factors.

Relation between geographic distances and morphometric divergence:
The significant correlation found between morphometric and geographic distances matrices are also supported by the geographic and paleogeographic analyses.In case of having an important genetic

A B
Gu ai ta ra Ri ve r GVP LOP CCVP Distribution area of E. pycnophylla Glacier control of leaf shape, this result could indicate that the processes involved with gene flow among the populations of E. pycnophylla are proportionally related with their geographical distances.These kinds of tendencies have been detected in other plant species with truncated gene flow due to restrictions in pollination (Lönn & Prentice 1995, Brunell & Whitkus 1993).In the case of E. pycnophylla, this result may have historic and present implications, which respectively are: 1) Closer populations were joint in a more immediate past than distant ones did, directly or indirectly by intermediate populations; 2) Closer populations as GPv and LOP, which have a minimal linear distance of 12km, have a higher probability of undergoing inbreeding contact than distant ones.This is because the spatial distances between them may not represent striking restrictions for pollination vectors.In a combined perspective of both hypotheses, current nearer populations may undergo gene flow although they were split in the past.¿Evidence of a vicariance event?: Nowadays, an increasing number of palaeoecological evidences shows that the periodic phenomenon of climate change during the Pleistocene (glacial-interglacial cycles) and their repercussion on the dynamics of the spatial distribution of Paramos in the northern Andes were a strong force that influenced many episodes of extinction, speciation, range contractions and expansions for many groups of plants, implying continuous events of vicariance for the biota (Simpson 1971, van der Hammen 1974, Simpson & Todzia 1990).Applying these ideas to the case of E. pycnophylla and based on our geographic and paleogeographic modeling, we suggest that the three studied populations are the result of a vicariance scenario.
Although our data of leaf shape variation, geographic and paleogeographic analysis of the distribution area of E. pycnophylla supports the mentioned explanations, a greater survey is still needed to test if this pattern is consistent when adding variation from marginal populations and from different locations within each population.
Due to the difficult access and mobilization in the Paramos, we only sampled one localized area in each of the three populations; this may have induced such pseudoreplication and makes necessary the addition of new sampling areas in contrasting microhabitats to average unwanted factors.But on the whole, the findings of this study offer a good opportunity to keep searching for biogeographic processes in response to climate change and geographical dynamics in plant species of Paramo ecosystems.

AKNOWLEDGMENTS
We are grateful to Belisario Cepeda and Maria Elena Solarte for contributing the seminal ideas of the project, to Jorge Bedoya and the Reserva Natural Pueblo viejo for their authorizations and logistic support during field work; to Germán Narváez and the TERRA group who supplied clue information for the geographical analyses; to Jorge Benavides for his useful revision and to IDEA WILD and the vIPRI-UDENAR which provided the financial support.

Fig. 1 .
Fig. 1.Schematic representation of an average adult of E.pycnophylla.The sampled leaves for shape analysis were those disposed between 45º and 90º with respect to the vertical axis of the plant.This guaranteed the selection of mature leaves.Image modified fromBenavides et al (2007).

Fig. 2 .
Fig. 2. Scatterplot of the CvA for the comparison of the leaf size among the three populations of Espeletia pycnophylla using linear morphometry.This graph represents the morphometric location of each leaf on the first and second canonical functions, showing the inter and intra population variation.The biplot projections represent the weights of the Total Leaf Length (TLL) and the Maximum Leaf Width (MLW) within each function.

Fig. 3 .
Fig. 3. Scatterplot of the CvA for the comparison of the leaf shape among the three populations of Espeletia pycnophylla using the Elliptic Fourier Descriptor.This graph represents the morphometric location of each leaf on the first and second canonical functions, showing inter and intra-population variation.

Fig. 4 .
Fig. 4. (A) Reconstruction of the distribution area of E. pycnophylla during the late glacial period (14 000 years ago) on the Southern Colombian Andes, according to van der Hammen's paleogeographic model (1974).(B) Estimation of the current distribution area of E. pycnophylla on the Southern Colombian Andes.Only the three studied populations are show in the map.

TABLE 1
Geographic characteristics of the studied populations CCvP: Chiles Cumbal volcanoes Population GvP: Galeras volcano Population.The respective research authorization for this population was obtained according to the permission DEC309/2000 issued by Parques Nacionales Naturales de Colombia.

TABLE 2
Discrimination statistics for the comparison of the leaf size and leaf shape among the three populations of Espeletia pycnophylla using Linear Morphometry