Ribosomal DNA revealed an extensive role of allopolyploidy in the radiation of Ulex L

We studied the phylogeny of Ulex L., a genus of spiny legumes, which its center of diversity in the Iberian Peninsula, using ribosomal DNA markers (rDNA), namely ETS, 5.8S and ITS (45S), and 5S intergenic spacer regions. One of the main findings was the presence of very different haplotypes in 5S-IGS genes and, to a less extent, in ETS and ITS, in seven polyploid taxa. We interpreted these results as an indication of hybrid origins and proposed allopolyploidy for U. argenteus ssp. subsericeus, U. australis ssp. australis, U. australis ssp. welwitschianus, U. densus, U. europaeus ssp. europaeus, U. europaeus ssp. latebracteatus, and U. jussiaei. These results reinforce an early hypothesis which stated that the radiation of Ulex occurred mainly by polyplodization. Phylograms showed two main clades, one grouping the hydrophilic U. gallii, U. breoganii and U. minor, and the other grouping the southern, xerophytic, taxa. The putative allopolyploids showed haplotypes, which grouped in both clades, indicating that allopolyploidy, occurred through hybridization from these hydrophilic and xerophytic lineages. The phylogenetic position of U. micranthus is not certain and it is discussed. The 5S-IGS showed to retain more polymorphisms than ETS gene or ITS markers. This result is compatible with the hypothesis that 5S rDNA region is less vulnerable to inter-loci concerted evolution than 45S, providing a more suitable marker for reconstructing histories of allopolyploid species. We discuss the taxonomic consequences of these results.


Introdution
Polyploidy is widely acknowledged as a major mechanism of speciation in plants. It can lead to speciation in several ways, either within single individuals or following hybridization between closely related populations (autopolyploidy), or from interspecific hybridization events (allopolyploidy) (1).
Nuclear ribosomal DNA has proved to be a valuable instrument to investigate polyploidy events. In eukaryotes, nuclear ribosomal DNA (rnDNA) belongs to two universal gene families: the 5S nrDNA and 45S nrDNA. Both are organized in tandem arrays. In most plants, they have different locations in the genome (2,3). The 45S array is widely utilized in phylogenetic studies, and is composed by an extensive intergenic spacer (IGS), an external transcribed spacer (ETS), three genes encoding the ribosomal units for 18S, 5.8S, and 28S, and two internal transcribed spacers (ITS1 and ITS2) (4). The other tandem repeat codes 5S nrRNA, and includes a small intergenic spacer (5).
After allopolyploidization, ribosomal DNA can undergo two processes, which can lead to the loss of parental haplotypes: concerted evolution and the loss of loci (6,7). Concerted evolution is a process which homogenizes the sequences within a species, making them more similar to each other than to the homologous sequences of other species reducing haplotypic variability. Some authors reported differences between the concerted evolution for the 45S and 5S arrays. Evidence was presented suggesting that, in 45S, the homogenization mechanisms operate within loci as well as between different rnDNA loci, while in 5S concerted evolution only operates within loci or, at least, homogenization mechanisms inter-locus are not very effective (3,8,9).
If concerted evolution is not completed, a hybrid species can retain two or more very divergent haplotypes originating from each parental taxa (10). If these haplotypes are found in other lineages in homozygosity this is often considered a consistent clue of a hybrid origin, and it is used to identify parental lineages (11,12). However, the presence of a single haplotype in a polyploid do not mandatorily indicates autopolyploidy, because this can also arise when the concerted evolution is complete and one of the haplotypes of the parental species is lost (13).
The Ulex L. genus belongs to the Fabaceae family and Genisteae Benth. tribe. It is native of Western Europe and northwest Africa and Iberian Peninsula has been considered as its center of diversity (14). They are spiny shrubs, characterized by the presence of leaves reduced to small spines (15). This genus includes more than 20 taxa (16,17), considering species and subspecies, and polyploidy is frequent. Following Cubas (18) and Cubas et al. (15), the genus includes, at least, 12 polyploid taxa, including tetraploids and hexaploids (Table 1).
Four previous studies tried to clarify the phylogenetic relationships of this genus, using the internal transcribed spacers of nuclear ribosomal DNA (nrDNA) and plastid markers. Ainouche et al. (14) found that Ulex is monophyletic and indicated that U. micranthus Lange is a monospecific lineage, separated from all other species. Similar results were achieved by Pardo et al. (19) with the same markers. However, all these markers showed low variability and generated poorly resolved phylograms (19). Ainouche et al. (14) suggested that the low level of sequence divergence reveals a recent radiation and that polyploidization could have played an important role in the speciation process (14). Cubas et al. (15), using chloroplast microsatellite markers, proposed an evolutionary scenario for the evolution of Moroccan and Andalusian (Spanish) species. According to these authors, Southern Iberian diploid taxa were closely related and the polyploid U. borgiae Rivas Martínez arose from the hybridization of U. baeticus Boiss. ssp. baeticus and U. baeticus ssp. scaber (Kunze) Cubas. Finally, Dias et al. (20) presented sequences of nrDNA of Ulex and advocate that U. europaeus L. subsp. europaeus is an allopolyploid (20). The knowledge on the phylogenetic relationships does not go much further.
The main goal of this paper is to study the phylogeny of Ulex. Considering the high number of polyploid taxa and the low genetic variability reported, we tested the hypothesis of Ainouche et al. (14), which stated that the radiation within this genus resulted mostly from polyploidy events, including allopolyploidy. We used three nuclear fragments of nuclear nrDNA, two from the 45S arrays, namely external transcribed spacer (ETS) and a fragment with two internal transcribed spacers (ITS1 and ITS2), the gene 5.8S (henceforward called ITS), and the 5S intergenic spacer (5S-IGS).

Plant material
Our sampling covered most of Ulex taxa, at least those classified as independent species, except the Moroccan U. congestus. Thirty-three populations from twenty-one taxa were examined ( Table  1). The sequences of U. europaeus ssp. europaeus of 45S rnDNA (ETS and ITS) were retrieved from GenBank (Accession numbers: KC806125.1 to KC806140.1, and KC832359.1 to KC832370.1). Genista falcata Brot., Stauracanthus genistoides (Brot.) Samp. and Stauracanthus boivinii (Webb) Samp., were used as outgroups. Sequences were deposited in GenBank (available at www.ncbi.nlm.nih.gov/). Samples consisted of fragments of phyllodes. Vouchers were deposited in the Herbarium of Instituto Superior de Agronomia (Lisbon). We also used preserved material from this herbarium. Criteria for morphological identification followed Cubas (18). In a first scoping, the marker 5S-IGS was amplified with the primers S1 and AS1 (21). At this step, some polyploids presented chromatograms with double peaks or long illegible fragments. In these cases, we design allele-specific primers (Table 2), a technique that uses primers with deliberate mismatches towards one of the template copies to make the reaction specific for one of the templates (22,23). It was impossible to design allele-specific primers for the 45S fragments (ETS and ITS), as well as to split two haplotypes of 5S-IGS of U. argenteus Welw. ex Webb ssp. subsericeus (Cout.) Rothm., mainly due to few variable sites. In these cases, chromatograms with double pics were read by the method described by Sousa-Santos et al. (24), using the indels (24). Table 2 -Primers for PCR amplification. S-1 and AS-1 following Sun et al. (21). The other primers were designed for this study All fragments were amplified with the same PCR protocol. PCR's included an initial denaturation of 2 min at 94º C, followed by 35 cycles of 1 min. denaturation (94º C), 45 s. annealing (53º C), 50 s. elongation at 72º C, and a final extension at 72º C (7 min.). PCR products were purified with the Microclean kit (Microzone Ltd), and sequenced using the same primers (in STABVIDA, Portugal).

Phylogenetic Analyses
Reconstruction of phylogenetic relationships was conducted using Maximum Parsimony (MP), and Bayesian Inference (BI). The BI was performed in Mr. Bayes 3.1 (27) using MCMC, with two independent runs of four Metropolis-coupled chains of two million generations each, to estimate the posterior probability distribution. Topologies were sampled every 100 generations, and a majority-rule consensus tree was estimated after discarding the first 2000 generations sampled. Sequences were partitioned by region (ITS and ETS), in the concatenated matrix. The MP analysis was conducted in PAUP 4.0b10 (28), using a heuristic search with random stepwise addition (1000 replicates), tree-bisection-reconnection branch swapping. All characters and character-state changes were equally weighted. Bootstrap analysis (1000 replicates) was used to assess the robustness of branches of the MP trees (29). The best-fit model of nucleotide evolution was obtained using Jmodeltest 3.7 (30) with the Akaike Information Criterion (AIC) (31).
Indels shared by two or more haplotypes were coded as binary characters according to the ''simple indel coding'' method of Simmons & Ochoterena (32) and used in parsimony and Bayesian searches (32). The binary gap data and the sequence data were treated as separate partitions.
For most taxa it was possible to concatenate the sequences of ETS and ITS, because we found only one haplotype or very similar ones from each marker, for each plant. This allowed the construction of a concatenated phylogram of 45S fragments with most taxa. However, some specimens of the hexaploid U. australis Clemente showed two very divergent haplotypes in both fragments (in ITS and in ETS), and it was impossible to know which haplotypes of ITS and ETS belong to the same genome, making it impossible to concatenate. The same situation occurred with the 45S sequences of the hexaploid U. europaeus ssp. europaeus, obtained in the GenBank. Therefore, to analyse the phylogeny of these two taxa, we constructed two separated phylograms (non-concatenated) one for ETS and another for ITS.
To analyse the possibility of occurrence of Long Branch Attraction (LBA) in the 45S phylograms, we re-run the matrix of 45S sequences with the Bayesian and Maximum Parsimony approaches. We also calculated the mean for the uncorrected distances, between U. micranthus, suspected to be affected by LBA, and we graphed these distances using a multidimensional scaling, based on the uncorrected pairwise genetic distances.

Results
Sequences from all taxa and studied populations, for 45S and 5S-IGS markers were obtained. The complete sequence data set consisted of 32 concatenated different sequences of ETS and ITS (i.e.: ITS2, 58S and ITS2), 32 different sequences of 5S-IGS, plus 45 and 31 sequences for the separated analysis of non-concatenated ETS and ITS fragments.
Overall, the Bayesian Inference approach resulted in better-solved and strongly supported phylogenetic trees. ETS and ITS resulted in a matrix with 459 and 574 characters, respectively. ETS showed higher variability and phylogenetic utility than ITS, with 44 parsimony informative characters, surpassing the 22 characters from ITS. The 5S-IGS spacer showed 38 parsimony informative characters, although it was the smallest fragment, with only 248 base pairs. Many plants shared the same haplotypes, or very similar ones, even belonging to different taxa. Phylograms are depicted in figures 2 to 5.     Posterior probabilities and bootstrap support values higher than 50 for Maximum Parsimony for each node were separated by a bar, in the left and in the right, respectively.
All diploid taxa showed only one sequence of 5S-IGS and 45S, or two very similar sequences, with a maximum of two double peaks, even when the fragments were amplified with non-allelespecific sequencing primers. When tested with different allele-specific sequencing primers for 5S-IGS haplotypes, diploids did not reveal different haplotypes, independently of the primer used.   (Figures 2 and 3).
In the 5S-IGS phylogram (Figure 2) The position of the diploid U. argenteus is incongruent, as it grouped with U. airensis in the 45S phylogram ( Figure 3, although this subclade is not well supported: Bayesian posterior probabilities and bootstrap are 89 and 57, respectively), but grouped with U. eriocladus in 5S-IGS phylogram (Figure 2).
The phylograms of 5S-IGS and 45S showed an incongruence in the position of U. micranthus. This species grouped as the first divergence of the genus in the 45S phylograms, but in 5S, grouped as a sister species of the Atlantic Clade. As stated before, we rerun the Bayesian and the Parsimony analyses without outgroups of the 45S matrix and found that U. micranthus grouped as an early branching of the Atlantic Clade, as happened in 5S-IGS phylogram (trees not shown). Furthermore, the mean of the uncorrected distance between U. micranthus and the Atlantic Clade (mean distance 2,486 + 0,049 %) is lower than the distance to the Mediterranean Clade (mean distance 2,699+0,001), and also lower than the mean between U. micranthus and all other taxa (mean distance 2,669+0,001), reinforcing the hypothesis of a closer phylogenetic relationship between U. micranthus and the Atlantic Clade. The figure 1 summarizes the relationships among species, using a multidimensional scaling, based on the uncorrected pairwise genetic distances (stress=0.148). The inspection of this figure corroborates the previous analysis and reveals that the 45S haplotypes of U. micranthus are closer to the Atlantic Clade, as expected.

Evidence of reticulated evolution
As previously mentioned, the data from nrDNA may prove informative to documenting the historical hybridization events, providing information on paternal lineages (10). Thus, since the early 90s, the presence of very divergent nrDNA haplotypes in polyploids is interpreted as an evidence of a hybrid origin, indicating that those taxa are allopolyploids (34,35 It is interesting to note that most of these polyploids register a higher haplotype diversity in 5S-IGS than in 45S. Some previous studies reported that, in 5S, the concerted evolution leads to the homogenization intra-loci, but does not act, or does it very weakly, inter-loci (3,8,9). Based on similar results, Mahelka et al. (2013) advocate that 5S nrDNA may provide a more suitable marker for reconstructing the histories of allopolyploid species than 45S. Therefore, the higher intraspecific diversity of 5S-IGS in Ulex polyploids could be interpreted as a consequence of the higher effectiveness of inter-loci concerted evolution in 45S and the slowness, or absence, of inter-loci concerted evolution in 5S.
Ulex europaeus ssp. europaeus showed also very divergent haplotypes in 45S and 5S-IGS (Figures 2, 4 and 5). We interpret this result as a confirmation of allohexaploidy, as stated previously by Dias et al.(20), based only in the 45S markers. However, differently from the other polyploids, the 45S haplotypes were more diverse (ITS and ETS phylograms, figures 4 and 5), grouping in three main locations, while 5S-IGS haplotypes just groups twice (Figure 2).
Beside concerted evolution, loss of loci may also have contributed to the absence of one parental 5S-IGS haplotypes, as well as 45S haplotypes, in allopolyploids. Furthermore, some authors suggested that this phenomenon is more likely to occur in higher ploidy levels (36,37). This could explain the absence of a third haplotype of 5S-IGS in hexaploid U. europaeus ssp. europaeus. Nevertheless, other hypotheses could be considered, including the mismatch of the allele-specific sequencing primers used to amplify 5S-IGS haplotypes, which may have prevented the amplification of a third 5S-IGS haplotype in U. europaeus ssp. europaeus.
Regarding the identification of parental lineages, we emphasize that U. europaeus ssp. europaeus presented phylogenetic connections with three different lineages (Figures 2, 4 and 5): 1) the lineage of U. minor (Atlantic Clade), connection supported by the data from 5S-IGS, ITS and ETS; 2) an early branching within the Mediterranean Clade, also supported by the data from 5S-IGS and ETS; 3) and a third lineage belonging to a subclade of Mediterranean Clade, supported only by ETS phylogram.
The hexaploid U. argenteus ssp. subsericeus showed three very similar 45S haplotypes, but also three very different in 5S-IGS. Both markers indicate genetic connection with U. erinaceus because these taxa grouped together in the 45S phylogram in a well-supported clade and shared one 5S-IGS haplotype. The other two 5S-IGS haplotypes from U. argenteus ssp. subsericeus grouped in two different clades: together with the diploids U. argenteus ssp. argenteus; and together with U. minor. These data suggest that it is an allohexaploid originated from these three lineages. Nowadays, U. argenteus ssp. argenteus, U. minor and U. erinaceus are not sympatric, but they occur in the extreme south of Iberian Peninsula, which suggests that they may have been in contact in a not-so-distant past.
Regarding U. australis, the presence of haplotypes of 5S-IGS, ITS and ETS, which grouped in three different positions in the phylograms, suggests that this species is an allohexaploid. However, each population was represented by a different combination of haplotypes, which can be explained, at least, in two ways: 1) the concerted evolution, or the loss of loci, fixed different haplotypes in different populations, or 2) under the designation of U. australis there are different unknown taxa with different phylogenetic origins. The former hypothesis is noteworthy because there are few reports on the literature of bidirectional concerted evolution of rnDNA in different populations of the same species (38,39). The second hypothesis seems unlikely because if each the combination of haplotypes represents a different taxon, there would be, at least, five different taxa, with no obvious morphological or ecological differences. However, an extensive sample and a detailed analysis of genetic markers and morphological characters through the geographical range of U. australis is required to fully clarify this question.
The tetraploid U. densus and the hexaploid U. jussiaei showed similar phylogenetic connections with two lineages: the lineage of the diploid U. minor (the Atlantic Clade), supported by one haplotype of 5S-IGS, and with a subclade of the Mediterranean Clade, together with the diploid U. airensis, supported by a second haplotype of 5S-IGS and by 45S. These results suggest that they are allopolyploids, sharing these two parental diploid lineages. Despite the similarity of haplotypes, these taxa have deep morphological and ecological differences (18). These differences could arise from a hypothetical third genome involved in the origin of the hexaploid U. jussiaei, although not detected in our analysis, from different genetic contributions from these two putative parental lineages, or from evolution subsequent to the hybridization events.
The phylogeny of the tetraploid U. europaeus ssp. latebracteatus is harder to clarify with our data. In the 5S-IGS phylogram, three haplotypes grouped in two distant positions: one in the Atlantic Clade and two as an early branching of the Mediterranean Clade. These results suggest that U. europaeus ssp. latebracteatus is an allotetraploid. However, there is an incongruence: in the 45S phylogram, it grouped within the Mediterranean Clade, although not as an early branching ( Figure 3). Despite this, U. europaeus ssp. latebracteatus grouped together with U. europaeus ssp. europaeus in the phylograms of 5S-IGS, ITS and ETS (Figures 2, 4 and 5), which suggests a close phylogenetic connection.
If our hypothesis is correct, the tetraploid U. europaeus ssp. latebracteatus arose from only two different genomes, while the hexaploidy of U. europaeus ssp. europaeus resulted from three different genomes, which means that these two taxa have different phylogenetic origins. Furthermore, Cubas & Pardo (40) showed that U. europaeus ssp. europaeus and U. europaeus ssp. latebracteatus have distinct morphological characters, which are maintained in sympatric populations, and argued that ploidy acts as an effective reproductive barrier (40). Therefore, they do not fit the biological nor the phylogenetic species concepts and should be assign as different species.

Other evolutionary insights
Five polyploids, namely: U. gallii ssp. gallii, U. gallii ssp. breoganii, U. eriocladus, U. erinaceus and U. borgiae, display low intraspecific diversity both in 45S and 5S-IGS haplotypes, and they do not show topological incongruences between the 45S and 5S-IGS phylograms. They can be autopolyploids, or they can be allopolyploids in which the homogenization process of nrDNA is completed, and the alleles of 45S and 5S-IGS from same parental species were lost. Therefore, the autopolyploid origin of these taxa needs confirmation by additional genetic markers. Ulex borgiae is probably an example of this, as Cubas et al. (18) (Figure 6), suggesting that the speciation of diploid lineages, in Ulex, occurred mainly by allopatric or parapatric speciation. Ulex parviflorus, both subspecies of U. baeticus, and U. canescens grouped in the same subclade, within the Mediterranean Clade, in both 5S-IGS and 45S phylograms. These results support a close phylogenetic relationship between these three taxa.
As stated before, in the phylograms of both markers, U. erinaceus groups as an early branching of Mediterranean Clade, not in this subclade of U. canescens. Accordingly, the morphological similarity between both stenoendemisms: U. canescens from Cape of Gata (Spain) and U. erinaceus, endemic from Cape of Sagres (Portugal), namely the presence of a dense branching and a cushion-like form (18), is probably a homoplasy, resulting from exposure to strong winds and dry habitats.
The incongruence between the positions of U. micranthus in the trees of 5S and 45S could result from Long Branch Attraction. This is a known methodological bias which leads to the erroneous grouping of two or more long branches as sister groups, in which a ingroup taxa can be 'pulled' to the base of the tree by long-branched outgroups (41,42). As far as we know, there is no method to prove that a particular topology results from Long Branch Attraction. Bergsten (2005) argued that the reconstruction of the trees without outgroups (unrooted trees) can indicate the corrected position of the taxon suspected of erroneous grouping. As pointed by Borowiec et al. (42), if the outgroups were to influence the branching order, we would expect the internal topology to be changed, if they were excluded. Indeed, the construction of unrooted trees of Parsimony and Bayesian approaches confirmed that U. micranthus grouped within the Atlantic Clade as an early branching, as it happened in the phylogram of 5S-IGS. Furthermore, the uncorrected distance between the U. micranthus and the Atlantic Clade is lower than the distance to the Mediterranean Clade, as well between U. micranthus and all other taxa. Future studies using data from other markers are needed to corroborate this hypothesis.
This is an interesting hypothesis, because, if confirmed, the first split within the genus Ulex would divide hydrophytes and xerophytes diploids. Allopolyploid taxa seem to arise from hybridation of these hydrophytes and xerophytes diploid lineages. We emphasize that U. micranthus, U. gallii ssp. gallii and U. gallii ssp. breoganii occur almost exclusively in the Atlantic Bioregion, although U. minor exists in both Atlantic and Mediterranean Bioregions. Nevertheless, under Mediterranean climate, U. minor always occurs in marshes or moist soils, usually with edaphic water compensation.

Taxonomic consequences
Our data can help to clarify some issues in Ulex taxonomy. Allopolyploidy provides a high degree of reproductive isolation from other ploidies, even with the progenitor species (43). For this reason, allopolyploids should not be assigned as subspecies, once there are not conspecific with lower ploidies, and should be assigned as species. This is the case of U. argenteus ssp. subsericeus, U. europaeus ssp. europaeus and U. europaeus ssp. latebracteatus.
Both markers (5S-IGS and 45S) confirm the splitting between Ulex airensis and U. parviflorus, as proposed by Espírito-Santo et al (33), and the validity of Ulex airensis as a different species. Lastly, Rothmaler (44) synonymized the stenoendemisms U. canescens and U. erinaceus, being followed by Guinea & Webb (16). Our results showed U. canescens is phylogenetically closer to U. parviflorus and U. baeticus, supporting the splitting between both stenoendemisms.

Conclusions
The main conclusion of this study is the confirmation of the hypothesis brought by Ainouche et al. (2003), that Ulex radiation occurred mainly by polyploidization, especially allopolyploidy. Indeed, allopolyploidy, as a speciation mechanism, is supported by ribosomal DNA markers, explaining the origin of, at least, seven taxa: the hexaploids U. australis ssp. australis, U. australis ssp. welwitschianus, U. argenteus ssp. subsericeus, U. jussiaei, U. europaeus ssp. europaeus, and the tetraploids U. densus and U. europaeus ssp. latebracteatus. This number is a third of the 21 studied taxa. However, our data are less informative on the origin of the polyploids U. gallii ssp. gallii, U. gallii ssp. breoganii, U. eriocladus, U. erinaceus and U. borgiae, which can be allopolyploids or autopolyploids. This approach should be complemented with additional genetic markers. This set of new data have important consequences in the taxonomy of Ulex, namely the assignment of specific rank to the allopolyploids, the confirmation of the validity of U. airensis as a species, and the splitting between U. canescens and U. erinaceus.
Our data also suggests, although do not prove, that the first splitting in this genus separate xerophytic and hydrophilic diploid. All the putative allopolyploids seem to be originated from the hybridization of these two early diverging lineages.
Other important result is the confirmation, for the genus Ulex, that 5S-IGS provided a more suitable marker for reconstructing histories of allopolyploid species than 45S, as stated by Mahelka et al. (3), probably because the inter-loci homogenization is less effective in 5S than in 45S repeats.