Genetic admixture and lineage separation in a southern Andean plant

Mountain orogeny has been a major factor in plant evolution in all continents by changing the landscape and climate, creating new habitats and ecological opportunities. In this study we found that diversity in two southern Andean Escallonia species is geographically structured and there is a deep divergence between infraspecific groups that could be associated with ancient evolutionary events like orogeny. We also found evidence of admixture, likely the result of hybridization at the margins of the parental species' distribution range.


Introduction
Mountain orogeny has been a major factor in plant evolution in all continents, and has been linked to recent diversification and speciation events (Hughes and Atchison 2015). The uplift of mountain ranges may change the landscape and climate, creating different environments and microclimates that provides new habitats (Linder 2008) and island-like ecological opportunities (Hughes and Eastwood 2006).
In South America, the Andes mountains have played these roles since their uplift during Pliocene/Miocene . Moreover, the Andes themselves could have played as a North-South corridor allowing the exchange of plant lineages; or as a barrier promoting vicariance. Consequently, these evolutionary processes triggered by the Andean uplift promoted a high speciation rate conducting to great biological diversity in South America (Antonelli et al. 2009;Hartley 2003;Kier et al. 2009;Luebert et al. 2011;Ortiz-Jaureguizar and Cladera 2006;Scherson et al. 2008).
The southern Andes have been identified as the origin of diversification of many groups (e.g. Chuquiraga: Ezcurra 2002, Valeriana: Bell et al. 2012, Oxalis: Heibl and Renner 2012, Heliotropium: Luebert et al. 2011); they are located in Argentina and Chile from 29 S, below Atacama desert to the austral tip of the continent at the Magellan and Fuegian Archipelagos.
During Pleistocene, climate oscillations have greatly influenced biodiversity in all continents (Hewitt 2000(Hewitt , 2004. In Patagonia, glacial cycles generated not only warmth-cold fluctuations in climate but also ice sheet expansions and retreats modelling landscape (Rabassa et al. 2011). As a consequence, plant species have accompanied those changes reacting with population contractions and expansions, long-distance range dispersal, new habitat colonization and in situ survival in refugia; examples of these processes have been studied for many Andean plants (e.g. Austrocedrus chilensis: Pastorino et al. 2004;Souto et al. 2012Souto et al. , 2015 Nothofagus species: Marchelli and Gallo 2006;Premoli et al. 2010;Soliani et al. 2012;Hordeum species: Jakob et al. 2009;Hypochaeris incana: Tremetsberger et al. 2009;Calceolaria polyrhiza: Cosacov et al. 2010, Eucryphia cordifolia;Segovia et al. 2012). Common processes for different lineages are suggested by two shared patterns a) high diversity zones corresponding to putative glacial refugia; and b) latitudinal phylogeographical breaks along the Andes and Patagonian steppe (Sérsic et al. 2011).
Among the plant lineages from South American mountain regions, the endemic genus Escallonia (Sede 2008;Sleumer 1968) is distributed along the whole Andean mountain range where it is highly diverse (Sleumer 1968;Sleumer and Correa 1984), and in different ecosystems in southern Brazil and central Argentina. Escallonia is the most numerous genus of the family Escalloniaceae comprising ca. 40 species; plant morphology displays high variation among species, and many diagnostic characters also show intraspecific variability (e.g. size and shape of leaves and floral organs, petal pigmentation, and presence of hairs and glands in different organs). Moreover, plants with intermediate morphology (IM) between species have been described (Sleumer 1968). Current taxonomy matches this pattern, with some species descriptions containing overlapping characters and ca. 20 species varieties described (Kausel 1953;Sleumer 1968;Sleumer and Correa 1984).
Molecular phylogenetic studies in Escallonia corroborated the monophyly of the genus. Five lineages were strongly supported and geographically structured, suggesting that their evolutionary history might be linked to orogenic processes in South America Zapata 2013). Particularly in the southern Andes, two independent lineages were evident : clade B with only two species (Escallonia virgata, restricted to southern Andes and Escallonia callcottiae, endemic of Juan Fernandez archipelago in Chile) and clade D with the remaining nine species sampled (E. alpina, E. florida, E. illinita, E. leucantha, E. myrtoidea, E. pulverulenta, E. revoluta, E. rosea and E. rubra). Within this major lineage, two species, E. alpina and E. rubra, are highly polymorphic and share the same distribution range along the Andes in Patagonia, although they only differ in maximum elevation: E. alpina reaches higher altitudes (up to ca. 2200 m) than E. rubra, which occurs up to ca. 1700 m. These species are differentiated mainly by the distribution of glands in leaves and flowers and by the type of inflorescence. For both species some varieties have been described on the basis of morphological characters and geographical distribution. There are two varieties described for E. alpina based on vegetative characters: E. alpina var. carmelitana only differs from the type variety in the colour of the stem and leaf size. A preliminary study on the genetic variability, including five populations, supported species boundaries and one population showing intermediate morphological characters was detected and corroborated by AFLP patterns .
The occurrence of individuals with IM between E. alpina and E. rubra [described as hybrids by Kausel (1953);Sleumer (1968) and Sleumer and Correa (1984)] along with the description of a population with genetic admixture ) and the lack of species exclusivity in phylogenetic reconstructions (Zapata 2013), indicate the presence of Escallonia hybrids in Patagonia.
Hybridization is a common phenomenon in plants and there are different views on its role in evolution: from not significant (Mayr 1992) to very relevant in speciation (Rieseberg 1997;Mallet 2007) or in adaptation to new environments (Rieseberg and Wendel 1993). Intermediate morphology is a strong indication of mixed ancestry (Stebbins 1959;Stace 1991); although hybrids are not always morphological intermediates (Rieseberg and Ellstrand 1993) and morphological intermediate individuals are not always hybrids (Morrell and Rieseberg 1998;Park et al. 2003). Hybridization has been frequently documented in the southern Andes (e.g. Fuchsia : Berry 1982;Discaria: Tortosa 1983) but there are few works that explore its evolutionary relevance in this region (Calceolaria: Sérsic et al. 2001;Caiophora: Ackermann et al. 2008;Nothofagus: Acosta and Premoli 2010;Soliani et al. 2012). A better understanding of natural hybrids between Escallonia species will elucidate taxonomic problems in the genus and will provide new hypotheses on its evolution.
Our aims are to characterize the genetic variation between E. rubra and E. alpina, by means of plastid DNA sequences and Amplified Fragment Length Polymorphism (AFLP). Additionally two southernmost populations with intermediate morphological characters between both species are included to investigate the genetic bases of their variation. We further analyze morphological variation in E. alpina (including the broadly distributed var. alpina and var. carmelitana, restricted to north Patagonia) using a geometric morphometrics analysis of leaf shape. We predict that E. rubra, E. alpina var. alpina and E. alpina var. carmelitana will have a high degree of genetic differentiation, while populations with intermediate morphology will show genetic admixture. Finally, leaf shape will be useful to differentiate the two varieties of E. alpina.

Plant material
Collection trips were undertaken on the eastern side of the Andes, in the Argentinean patagonian region (  (Thiers 2015).
For population genetic analyses we collected eight individuals from each location. Individuals were collected at least 15 m apart in order to avoid collecting close relatives. When eight individuals were not available, at least four were collected for each population. Fresh leaves were separated and dried with silica gel.
For morphological analyses we used pressed leaves of E. alpina including individuals from both varieties used in the population genetic analyses and additional herbarium material identified either as var. alpina or var. carmelitana [See Supporting Information 1a].

Total DNA isolation
Genomic DNA was extracted from silica dried leaves following a cetyltrimethylammonium bromide (CTAB) protocol (Doyle and Doyle 1987) with some modifications. Twenty mg of dried leaves were cooled by immersing in liquid nitrogen and then ground into fine powder. The samples were transferred to 1.5mL tubes with 700 mL warm CTAB buffer [2 % CTAB, 100 mM Tris-HCl pH 8.0, 1.4 M NaCl, 20 mM EDTA, 2 % polyvinylpyrrolidone] and kept in a water bath at 65 C for an hour with continuous shaking and mixing by inversion. Tubes were then removed from the water bath and left to cool at room temperature for 10 min, then 700 mL of chloroform: isoamyl alcohol (24:1) was added and the contents were mixed by vortexing. Centrifugation was carried out at 12 000 rpm for 1 min and the aqueous (top) layer transferred into a new 1.5 mL tube. A solution of NaCl (5M; 300 mL) and ice-cold isopropanol (600 mL) were added and mixed gently by inversion; tubes were incubated for 30 min at À20 C. Centrifugation was carried out at 12 000 rpm for 1 min, then the supernatant was discarded and the pellet was washed twice with ice-cold ethanol (70 %; 700 lL). After drying at room temperature, the pellet was suspended in 50 lL of 10:1 TE (10 mM Tris: 1 mM EDTA) buffer.

Plastid DNA
Amplification and sequencing. The plastid intergenic spacers trnS-trnG, 3 0 trnV-ndhC and the ndhF gene [see Supporting Information 2] were amplified with a profile consisting of 94 C for 3 min, followed by 30 cycles of 94 C for 1 min, 52 C for 1 min and 72 C for 1 min. Polymerase chain reactions (PCRs) were performed in a final volume of 25 lL with 50 ng of DNA template, 0.2 lM of each primer, 25 lM deoxynucleotide triphosphates, 5 mM MgCl2, 1Â Taq buffer and 1.5 units of Taq polymerase provided by Invitrogen-Life Technologies (Brazil). Automated sequencing was performed by Macrogen Inc. (South Korea). We edited and assembled electropherograms in BioEdit 7.0.9 (Hall 1999 Haplotype network. To understand the relation between the species of Escallonia distributed in the southern Andes, an haplotype network was constructed using statistical parsimony (0.95 probability connection limit) and the genealogical reconstruction algorithms of Templeton et al. (1992) as implemented in the R package pegas 0.8-1 (Paradis 2010 AFLP protocol was performed essentially as described by Vos et al. (1995), with fluorescent labelled primers that allowed automatic detection of the amplified fragments. Genomic DNA samples (50-100 ng) were digested to completion with EcoRI and MseI and the fragment ends were ligated to EcoRI-and MseI-specific adaptors [Supporting Information 2] in a single reaction for three hours at 37 C. The digestion-ligation products were diluted 20-fold into 10 mM Tris-HCl, 0.1 mM EDTA (pH 8.0) and amplified using EcoRI-A and MseI-C as pre-selective primers. The resulting template was diluted 20-fold prior to amplification with selective primers: EcoRI (FAM)-ACT and MseI-CAC [Supporting Information 2]. The fluorescence-labelled selected amplification products were separated on a sequencer with an internal size standard at Macrogen Inc. (South Korea). Nine random samples (9.6 % of individuals) were duplicated in order to assess reproducibility. Fragment-size fluorescent data was visualized using Peak Scanner (Applied Biosystems, available at http:// products.invitrogen.com/ivgn/product/4381867, last accessed 31 May 2016) and automatically scored using the R package RawGeno v2.0-1 (Arrigo et al. 2012 available at http://sourceforge.net/projects/rawgeno, last accessed 31 May 2016). Peaks between 70 and 500 bp with fluorescence higher than 120 relative fluorescent units (RFU) were retained and non reproducible fragments were removed. Presence or absence of each fragment was recorded for each individual obtaining a binary matrix [Supporting Information 3]. Pearson correlation coefficients between each fragment size and its frequency were calculated in order to assess the possibility of homoplasy (Vekemans et al. 2002).

Genetic structure
The AFLP matrix was used to infer genetic structure from Escallonia populations in order to assess the identity of IM populations and also to investigate the relation between E. alpina varieties. Differentiation due to genetic structure was tested using Wright's fixation index (Fst) for all populations and for separate species. Additionally a matrix of pairwise Fst genetic distances between all populations was constructed. Principal Coordinate Analysis (PCoA, Gower 1966) was performed using Jaccard genetic distances between individuals calculated from the binary matrix using software FAMD v1.108 (Schlü ter and Harris 2006); this multivariate analysis projects pairwise genetic distances upon a set of derived orthogonal axes, this reduced dimensional space allows the recognition of genetically similar individuals. Additionally, a bayesian analysis of population structure was performed using software STRUCTURE v2.3 (Falush et al. 2003(Falush et al. , 2007Pritchard et al. 2000) to infer the distribution and number of genetic clusters (K). The default model was used which assumes correlated frequencies between clusters and allows individuals to have a mixed ancestry. The log-likelihood probability of the data was calculated for each possible K value from 1 to 15 using 20 runs of 1 000 000 MCMC iterations following a burn in period of 100 000 iterations. Convergence of parameters between different runs was analyzed using Tracer v1.4 (Rambaut et al. 2013). The best fit number of clusters was calculated according to Evanno et al. (2005) using STRUCTURE HARVESTER (Earl and vonHoldt 2011).

Leaf shape analysis of E. alpina varieties
In order to characterize both varieties of E. alpina, leaf shape of var. alpina and var. carmelitana was described by traditional methods, using specimens from our collections (Table 1) and herbarium material [Supporting Information 1b]. Additionally, leaf shape was studied by performing an Elliptic Fourier Analysis (Kuhl and Giardina 1982). A total of 88 leaves from 35 collections of E. alpina including var. alpina and var. carmelitana (see Table 1 and Supporting Information 1b) were cut, hydrated, pressed and photographed using a Cannon PowerShot G11 digital camera over a clear background. The digital files were manipulated to obtain a binary image using Fiji (Schindelin et al. 2012). Final digitized images are available as Supporting Information 4. The R package Momocs V0.9 (Bonhomme et al. 2014) was used to extract the outlines as coordinates and to perform succeeding analysis. Each leaf shape was reconstructed using the first 15 harmonics, this number was selected by comparing the original outline and several reconstructed shapes using an increasing harmonics number.
A principal component analysis (PCA, Hotelling 1933) was performed on the Fourier coefficients of the leaves shape, and ordination of samples was then plotted on the first two principal components (PCs) axes.
The PCs explaining >99 % of variance were used as shape variables in subsequent analyses. A Linear Discriminant Analysis (LDA, equivalent to Canonical Variate Analysis) was carried to show maximum separation among E. alpina varieties. The LDA matrix was tested using a leave-one-out cross validation procedure to compute the percentage of correctly assigned individuals for each variety (Martens and Dardenne 1998). Additionally, PCs were subjected to multi-variate analyses of variance (MANOVA) to evaluate the significance of the separation between the two groups.

Sequence variation and haplotype network
The PCR amplification of intergenic spacers resulted in fragments of 740 bp (trnS-trnG) and 602 bp (trnV-ndhC), and the amplification of ndhF gene (incomplete sequence) resulted in a 1803 bp-fragment. All accessions, including multiple accessions for one species, showed exclusive haplotypes (Fig. 3).
In the network (Fig. 3), one major group is composed of E. leucantha, E. illinita, E. rosea, E. revoluta, E. myrtoidea, E. alpina (including var. carmelitana), E. rubra, and E. florida, which are separated by few mutational steps. The remaining species are separated by a long chain of mutational steps from the main group.
In particular, new sequences of E. alpina var. alpina and var. carmelitana were included in the major group together with those downloaded from the GenBank database. Varieties were grouped together, although neither of them shares a unique haplotype. Moreover, all accessions of E. alpina var. alpina are more closely related to other species like E. rubra and E. florida than they are to E. alpina var. carmelitana accessions; both varieties were separated by at least 7 mutational steps.

AFLP
Non replicated peaks (4.8 %) were eliminated and a matrix with 174 fragments for the 94 individuals sampled was obtained. Linear regression of fragment size and fragment frequency was not significant (R 2 ¼ À0.153; P ¼ 0.07). Escallonia rubra and E. alpina (including both varieties) had, respectively, 7 and 19 exclusive fragments. At the same time E. alpina var. alpina and E. alpina var. carmelitana had seven exclusive fragments each.

Genetic structure
When comparing all populations, Wright's fixation index showed a significant degree of genetic structure (Fst ¼ 0.17, P < 0.001). As this Fst value encompassed differences among taxonomic units we performed separate analyses for E. rubra (Fst ¼ 0.13, P < < 0.001) and E. alpina var. alpina (Fst ¼ 0.03, P ¼ 0.002). Escallonia rubra showed a much higher level of genetic structure, although both were significantly different from zero. Pairwise differences among populations are shown in Table 2. Fst distances were generally low when comparing populations from the same taxonomic unit (Fst ¼ 0.01-0.11) and higher when comparing populations from different taxonomic units (Fst ¼ 0.10-0.39). The notable exception was population 224 that showed high Fst values (0.2-0.39) when compared with any other population. The two populations with IM showed low levels of genetic differentiation when compared to populations of E. alpina var. alpina and E. rubra (Fst ¼ 0.01-0.13) but higher values when compared to E. alpina var. carmelitana (0.17-0.22).
The distribution of individuals among the first and second principal coordinates of the PCoA analysis is shown in Figure 4A; the first 2 eigenvalues are 0.19 and 0.10, and both account for 28.65 % of total variation. Escallonia. rubra and E. alpina individuals are separated along the first axis, while individuals that share morphological diagnostic characters of both species are spread in an intermediate position (Fig. 4). Along the second axis, E. alpina var. carmelitana populations are discretely grouped and well separated from the rest of E. alpina and E. rubra populations.
STRUCTURE results are shown in Figure 4B. Individuals were assigned to three genetic groups, as suggested by Evanno's method, mostly corresponding with taxonomy:   (Table 3) revealed a significant level of differentiation between E. rubra, E. alpina var. alpina and E. alpina var. carmelitana (25.67 % P < 0.001). The highest percentage of variation was found within populations (63.65 %, P < 0.001).

Leaf shape
Morphological examination allowed us to characterize leaves of E. alpina var. carmelitana as narrowly obovate to spatulate, lanceolate, tapering gradually to the base, with the apex acute to obtuse, rarely truncate (see Fig.  2C), while those of the type variety are mostly obovate, sometimes suborbicular-cuneate or spatulatelanceolate, with the apex typically obtuse, sometimes subacute, rarely truncate or invaginated (see Fig. 2B).
Geometric morphometry further reinforces these observations ( Fig. 5A and B); the PCA of the Elliptic Fourier descriptors shows two groups of leaves that correspond to E. alpina var. alpina and E. alpina var. carmelitana, although they are partially overlapped (Fig. 5A). The first two PCs retained 89.7 % of total variation. As seen from the reconstruction of the shapes along the first axis, the major source of leaf outline variation is anisotropy (length to width ratio). The mean leaf shape for each group is shown in Figure 5B. Cross-validation performed over LDA values (using the first nine PCs) was highly successful for both varieties, as 93 % of all leaves was well classified (E. alpina var. alpina 94.4 % and E. alpina var. carmelitana 91.7 %). A MANOVA test also showed significant differences between the two varieties (P < 0.001).

Discussion
Southern populations show genetic admixture and could be the result of hybridization between E. alpina var. alpina and E. rubra The results of our analyses at the population level support preliminary evidence of genetic admixture ) and the morphological hypothesis of Sleumer (1968) Fig. 2B) and E. rubra (stipitate glands on hypanthium, pedicels and young stems, and resinous spots on the adaxial leaf surface, Fig. 2A). We were able to identify herbarium collections with IM from the same geographic location [Supporting Information 1a]; moreover, most of the collections identified by Sleumer (1968) as natural hybrids are located in southwestern Santa Cruz in the proximity of Lago Argentino. Accidental interspecific cross pollination is likely to occur where more than one species are in sympatry and share pollinators. Floral biology could help elucidate the origin of the hybrid populations. To date, there are few studies on the floral biology of Escallonia species; Valdivia and Niemeyer (2006) found that E. myrtoidea is commonly visited by generalist pollinators such as bees and this is in accordance with our observations in the field. Although the parental species have also been recorded in the same region (Sleumer 1968), we only observed and collected populations with IM.
New questions arise if we assume these southern populations are indeed of hybrid origin. Why are hybrids more common in the southernmost locations? And why are these populations mostly composed by hybrids? One plausible answer is that hybrid fitness could be dependent on habitat: in most of the parental distribution range, hybrid fitness may be significantly lower than their parental fitness, promoting isolation. On the other hand, in marginal habitats selective pressures are different, and hybrids could be just as fit or even favoured by natural selection. In our study, IM populations were found in the southern limit of both parental species distribution range. The role for hybridization in evolution in marginal or altered habitats has been discussed in previous works Burke and Arnold (2001), and Rieseberg et al. (2003). A similar process could have occurred if parental species had shared glacial refugia. Reproductive isolation mechanisms that may have arisen during and after speciation could have been lost later as a consequence of reduced or fragmented habitats and climatic changes. This scenario has been proposed for plant species in Europe (e.g. Heuertz et al. 2006;Palmé et al. 2004 The fact that current hybrid populations are composed mostly by admixtured individuals suggests that gene flow was not an isolated event in the region and it may have played a significant role during the evolution of Escallonia species; this could be relevant for reconstructing processes as demographic contractions and expansions, colonization routes and possible survival in refugia. A more comprehensive phylogeographycal analysis is necessary to answer questions on the history and origin of these putative hybrid populations, and an experimental design would be useful to test hypotheses concerning hybridization mechanisms. According to the protologue (Kausel 1953), E. alpina var. carmelitana was differentiated from the type variety by the leaf size (which is highly variable even within individuals) and stem colour (which is affected by age and state of preservation). During the morphological study of collections of E. alpina var. carmelitana we observed differences in leaf shape rather than size. Leaf shape differences were strongly corroborated by geometric morphometrics analysis and this new character will be useful for species identification purposes.
E. alpina var. carmelitana is a clearly distinct group, with genetic and morphological differences and with a restricted geographic distribution, which deserves to be considered an independent entity. Populations of E. alpina var. carmelitana are located between 35 and 38 S, 70 W: this region is characterized by two high mountain ranges (the Andes and Cordillera del Viento), volcanoes, e.g. Copahue, and valleys; the zone is well irrigated, with several tributaries of the Neuquén river, although there are only small lakes, e.g. Varvarco Campos and Caviahue (Bran et al. 2002). A high level of endemism has historically been reported for this zone of northern Patagonia (Cabrera and Willink 1973;Simpson and Neff 1985). Moreover, it has been proposed that high genetic diversity of populations in northern Patagonia may be a consequence of in situ survival in refugia during glacial cycles Nevertheless, the strong pattern of genetic, morphological and geographical differentiation found in E. alpina var. alpina vs E. alpina var. carmelitana do not seem to be a consequence of recent periods of isolation and expansion due to climatic events during the Quaternary. High mountains may have acted as dispersal barriers (Struwe et al. 2009) promoting isolation and vicariance, even generating deep divergence among species and high levels of biological diversity (Hughes and Eastwood 2006;Pennington et al. 2010). Sérsic et al. (2011 and literature herein) found shared patterns of genetic distribution along the Andes, both for plants and animals: particularly, a high diversity zone around 35 -39 S. It has been hypothesized that many Patagonian plants have diversified in this region as a consequence of Miocene-Pliocene orogeny and associated tectonic processes; e.g. Calceolaria. polyrhiza (Cosacov et al. 2010), Hordeum (Jakob et al. 2009) and Hypochaeris. incana (Tremetsberger et al. 2009). Our results, combining morphology, plastid sequences and AFLP, together with the geographically structured phylogenies of Escallonia support the hypothesis that Andean orogeny has played an important role in the diversification of the genus.

Conclusions
New evidence of genetic admixture in Escallonia populations might be a result of interspecific hybridization. Further studies on ecology, pollination and floral biology could help to understand the role of interspecific gene flow in the evolution of the genus. Moreover, a comprehensive phylogeographic work, using co-dominant markers, will give hindsight in the recent evolution of Escallonia species and their interaction over time.
Additionally, we conclude that E. alpina var. alpina and E. alpina var. carmelitana are distinct lineages, and taxonomy should be revised to reflect their separation. Escallonia species, as currently circumscribed, could be concealing a richer diversity. We expect that new studies, combining morphology and genetics, will improve our understanding of biological diversity in the genus and general trends of plant evolution in Patagonia.

Accesion Numbers
All new sequences were deposited in GenBank with accession numbers KU759574-KU759579. performed the geometric morphometric analysis; S.S. and S.M. analysed the data and led the writing.