Studied species

Picris hieracioides is a biennial to perennial, morphologically highly variable species that includes two infraspecific taxa, P. h. subsp. umbellata and P. h. subsp. hieracioides. The subspecies of P. hieracioides are morphologically distinguishable by a combination of several morphological traits on both vegetative and reproductive organs. This was definitely proved also by cultivation experiments under environmentally homogeneous conditions [24]. Beside this, they also differ in longevity, ecological requirements, and geographic distribution [22], [24]. Additionally, each subspecies harbours two highly allopatric genetic (AFLP) and morphologically cryptic groups [24]. Picris hispidissima is a biennial species (M. Slovák, unpublished data) that differs morphologically from P. hieracioides in having a pectinate-ciliate indumentum of involucral bracts and inflation of the peduncle [23]. The last studied species, P. olympica, is a perennial characterized by a low caespitose growth with quite short ascending scapose-like stems bearing only a few capitula. It is considered to be morphologically closely related to P. hieracioides [21], [29]. The investigated species differ conspicuously in their ecological requirements and geographic ranges, with P. hieracioides being the most widespread (distributed across the major part of Europe and Asia Minor) and having wide ecological amplitude [21], [24]. The nominate P. h. subs. hieracioides is mainly biennial, growing on open, dry, xerothermous biotopes with no bedrock preference and thriving in lowlands or lower mountains throughout all of Europe. In contrast, P. h. subsp. umbellata is a short-lived perennial herb inhabiting humid tall herb communities from montane forests to sub-alpine communities across European calcareous mountain ranges. Both P. hispidissima and P. olympica have restricted distributions. The former occurs in dry, calcareous rocky slopes and crevices in the Adriatic coastal mountains from Croatia to Montenegro, while the latter is restricted to a few high mountain ranges in western Anatolia where it inhabits open, craggy grasslands with alkaline bedrocks at alpine levels [21], [29]. All species, but especially P. hieracioides, have spread from their natural biotopes into man-made habitats and migrated along anthropogenic corridors ([22], [24]; M. Slovák unpublished data).

Sampling strategy and plant material

In order to perform the reconstruction of phylogenetic relationships among the studied species and to reveal their genetic patterns, we aimed to include as many populations as possible from the entire geographic ranges of P. hieracioides, P. hispidissima, and P. olympica (Table S1). This, however, was at the expense of sample size at the population level [30], [31], [32]. We also co-analyzed a limited number of samples of the closely related Asian taxa P. japonica (2 populations) and P. nuristanica (1 population). We did not include samples of P. abyssinica since it was proved that it is genetically highly divergent from P. hieracioides and P. nuristanica but is more closely related to other African taxa [26]. Each population of the investigated species had more than several hundreds of individuals, and the taxa under study are neither endangered nor protected and no specific permits were required to collect plant samples at the study sites. We sequenced one to three accessions per population for the nuclear ITS and the plastid rpl32-trnL^(UAG) regions (Table S1). To cover as high proportion of genetic variation of analysed taxa as possible, but simultaneously with respect to sampling limitation at population level, we decided to analyse more than one individual for populations where non-homogenized ITS sequences indicated existence of more than a single allele. The widespread P. hieracioides including populations morphologically corresponding to both subspecies, P. h. subsp. hieracioides and P. h. subsp. umbellata, was collected from 122 sample sites covering both subspecies' entire geographic distributions, and henceforth these sites are referred to as populations (Table S1). Twelve populations of P. hispidissima from its entire distribution area were sampled [23]. Three of the five currently recognized populations of P. olympica were also collected for the purposes of the present study.

Intermediate morphological variation was identified in a few populations of P. hieracioides which were mostly located in zones where the two subspecies or their lineages meet [24]. Since the extent of intermediate morphological variation across these populations was wide, it was impossible to unequivocally sort some of them into discrete categories. Therefore, in order to avoid formation of artificial entities, which might lead to subsequent biasness of results and interpretations, we assigned particular populations a priori not to accepted subspecies, but we considered all of them as P. hieracioides instead.

Selection of appropriate outgroup taxa was based on previously published studies [21], [26] and our preliminary phylogenetic analyses. Thus a priori selected outgroup comprised members of all sections of the genus Picris (P. galilaea (Boiss.) Eig., P. capuligera (Durieu) Walp., P. pauciflora Willd., P. rhagodioloides (L.) Desf., P. scaberrima Ten., P. sinuata (Lam.) Lack, and P. strigosa M. Bieb.) and two members of the closely related genus Helminthotheca (H. aculeata (Vahl) Lack and H. echioides (L.) Holub). Voucher specimens were deposited in the herbarium of the Institute of Botany at the Slovak Academy of Science (SAV).

DNA sequence markers, DNA extraction and PCR amplification

We utilised sequences of the nuclear ribosomal internal transcribed spacer (ITS1-5.8S-ITS2 region) and the non-coding cpDNA intergenic spacer rpl32-trnL^(UAG) [33]. Several single-copy nuclear genes (A25, A28 and B12) [34] and another non-coding chloroplast regions (3′trnV^(UAC)-ndhC, ndhF-rpl32 and trnH^(GUG)-psbA) [33], [35] were tested during preliminary screening. However, these were proved to be unsuitable for the present study due to low levels of variation and/or problems encountered during their PCR amplification.

Total genomic DNA was extracted from silica gel-dried leaf tissue by the DNeasy plant mini kit (Qiagen, Hilden, Germany), following the manufacturer's protocol. The ITS region and the plastid rpl32-trnL^(UAG) spacer (referred as cpDNA henceforth in the text) were amplified using the standard PCR reaction, employing the following universal primer pairs: P1A, P4 and internal ones, ITS2 and ITS3, if necessary [36], [37] and trnL-retF and rpl32-retR [39].

Amplifications were carried out in a Mastercycler ep Gradient S thermal cycler (Eppendorf, Hamburg, Germany), using PCR reaction volumes of 25 µL consisting of 0.75 U of Pfu polymerase (Fermentas, St. Leon-Rot, Germany), 0.2 mM of each dNTP, 0.2 µM of each primer, 1 µL of DNA template, and reaction buffer containing 2 mM of MgSO₄ (Fermentas, St. Leon-Rot, Germany). For ITS amplification, we used the following thermocycler program: preheating at 94°C for 3 min, then running it for 35 cycles at 94°C for 30 sec, then at 50°C for 30 sec, and 72°C for 60 sec, with the final extension at 72°C for 10 min, and cooling at 4°C. The cpDNA region was amplified using the “rpl16” program, as described by [38]. PCR products were purified using the NucleoSpin Extract II kit (Macherey-Nagel, Düren, Germany). Cycle sequencing reactions were carried out using the same primers, and sequencing was performed on an ABI PRISM 3130xl sequencer at BITCET Consortium, Comenius University in Bratislava. The sequencing of a few accessions repeatedly failed for the ITS region and these were therefore excluded from further analyses (Table S1).

Sequence alignments and phylogenetic analysis

Sequences were edited and aligned manually using the BioEdit program version 7.0.4.1 [39]. In case of ITS sequences, the electroferograms were carefully inspected for intra-individual polymorphic sites (IPS) having more than one signal (cf [40]). These were labelled with NC-IUPAC ambiguity codes. Polymorphic positions within ITS sequences, in which both bases were detected also separately in different accessions elsewhere in the alignment, were considered additive polymorphic sites (APS) [40].

We believe that the observed intra-individual variation did not arise from PCR errors, because numerous accessions of P. hieracioides and also other taxa possessed homogenized ITS sequences. The geographic pattern of homogeneous and heterogeneous ITS sequences was obviously not random, indicating a genuine variation rather than erroneous signals caused by inaccurate amplification. Although the intra-individual ITS variation might be revealed by cloning procedure, this approach, on the other hand, may result in amplification of numerous ITS sequence types which do not represent relevant variability [41]. Incomplete concerted evolution and/or recombination, processes operating on the multiple-copy regions like ITS, increase the number of unwanted sequence types making the ITS phylogenetic analysis even more complex. Furthermore, with our current sampling involving hundreds of accessions exhibiting intra-individual ITS polymorphism, the cloning procedure would be out of the reasonable solution. As long as most of the sequences were possible to obtain by direct sequencing and polymorphic sites could be designated by IUPAC ambiguity codes, we decided to separately analyse two ITS datasets, one including and one excluding the individuals with APS (see below).

Altogether we generated three alignments: the first one was designated as ITS_1 and involved 117 sequences without APS (Alignment S1. ITS_1). The second one was designated as ITS_2 and comprised 216 sequences coming from all accessions analysed for the ITS region and including all sequences with APS (Alignment S2. ITS_2). The second alignment was constructed after considering that the intra-individual polymorphic sites, especially those with APS, may significantly influence the hierarchical structure of the phylogenetic trees.

Finally, the third cpDNA alignment comprised all plastid sequences, in particular 219 accessions (Alignment S3. cpDNA). In order to detect possible incongruence between particular markers, the incongruence length difference (ILD) test [42], as implemented in the “partition homogeneity” test of PAUP* (1000 replicates, Mulpars off, outgroups included), was employed. Although topologies of the ITS and cpDNA trees were shown to be incongruent, statistical support especially at their basal nodes was comparatively poor, causing their main branching patterns to be unsupported. Therefore, we decided to take the advantage of synergistic effect of combining datasets and generated also concatenated alignments. The concatenated alignment, ITS_cpDNA, contained 117 ITS sequences without APS and corresponding cpDNA sequences (Alignment S4. ITS_cpDNA).

Indels were not coded separately but treated as missing data. The aligned datasets were analysed independently using the following phylogenetic approaches:

1. Maximum parsimony (MP) phylogenetic analysis was performed with the heuristic search option in PAUP* version 4.0b10 [43]. The following settings were utilised: accelerated character transformation (ACCTRAN), gaps treated as missing data, single-site polymorphisms determined uncertainties, tree construction with stepwise addition, 1000 bootstrap replicates with random taxon addition, tree bisection-reconnection (TBR) branch swapping, and retention of multiple trees found during branch swapping (MULTREES option in effect). The identical sequences were merged in McClade version 4.0 PPC [44] to reduce computation time. Clade support was calculated via bootstrap analyses using 10000 re-samplings done with the fast heuristic search in PAUP*. Bootstrap support was categorized according to the following criteria: strong (>85%), moderate (70%–85%), weak (50%–69%), or poor (<50%).
2. Bayesian inference (BI) was run in MrBayes version 3.1.2, using the Markov Chain Monte Carlo algorithm (MCMC) [45]. Bayesian analyses were performed on the CIPRES Portal version 1.15 [46]. Prior to Bayesian analyses, the most appropriate nucleotide substitution models were chosen, using the Akaike Information Criterion (AIC) in jMODELTEST version 0.0.1 [47], [48]. Evolutionary models were calculated for each part of the datasets separately. Specifically, the ITS datasets contained partitions corresponding to the ITS1 and ITS2 spacers and the 5.8S rRNA gene, while concatenated datasets included besides the three ITS region partitions also a fourth partition represented by the cpDNA region. The following models or model combinations were found to be the most appropriate for the datasets studied: (1) ITS_1 and ITS_2 datasets – the SYM + G model for ITS1 and ITS2 sequences and the K80 model for 5.8S rRNA gene sequences; (2) cpDNA dataset – the TVM + I + G model; (3) concatenated ITS_cpDNA – the SYM + G model for ITS1 and ITS2, the K80 model for the 5.8S rRNA gene sequences, and the TVM + I + G model for the cpDNA partition.

All BI analyses were run with four independent Metropolis-coupled MCMC chains (three heated and one cold chain) for ten to twenty five million generations and sampled every 1000^th generation. The first 25% of sampled trees were regarded as ‘burn-in’ trees and were discarded prior to reconstruction of a 50% majority-rule consensus tree. Stationarity was confirmed by checking convergence diagnostic parameters. Specifically, the average standard deviation of split frequencies was lower than 0.01 in all cases; the plots of generations versus log probability of the observed data showed no obvious trends; and the Potential Scale Reduction Factor (PSRF) approached 1. Finally, topologies and node posterior probability values were compared among the runs. The topologies were stabilized among all datasets with only minute differences in branching pattern of terminal clades. Nodes with posterior probability (PP) values of 0.90 and above were regarded as significant and those with PP values below 0.90 regarded as non-significant.

1. 3. Net-like approaches were used to identify and display potential contradictory signals in the datasets. All alignments including ITS_2 containing intra-individual polymorphic sites were analysed using the neighbour-net analysis of [49] in SplitsTree version 4.10 with uncorrected P-distance and default settings. To visualize the relationships among cpDNA haplotypes and to detect possible ancestral polymorphism, the cpDNA dataset was subjected to haplotype network analysis based on the parsimony method of [50] using TCS version 1.21 [51] limited to 30 steps of parsimonious connection in creating the network.
2. 4. In order to precisely specify the amount of phylogenetically informative signal in the datasets, we analysed spectrum of supporting nucleotide positions. There are three groups of supporting positions recognized by Wägelle and Rödding [52]: (1) symmetrical or binary positions have two different character states in functional outgroup and ingroup and thus support both group of a split equally; (2) asymmetrical positions support only one group which possesses the same nucleotide at particular position, while the other group harbours different and more than one character state at this position; (3) noisy positions include same character states present in all sequences of the functional ingroup but also at least in one sequence of the functional out-group and thus represent convergences or chance similarities between ingroup and outgroup, or alternatively ingroup autapomorphies.
3. 5. We also apply a Bayesian coalescent-based approach to estimate a species tree employing *BEAST as implemented in the program BEAST version 1.7.4 [18]. Two input files, one with 117 homogeneous ITS sequences (Alignment S1. ITS_1) and the second one containing corresponding cpDNA sequences, were used for the BEAST analysis. Populations assigned to particular taxa were used as OTU's. The input file for *BEAST was created in BEAUti version 1.7.4, with the following settings: two data partitions (corresponding to the two loci), the best-fit evolutionary model for each partition as determined by jMODELTEST, uncorrelated lognormal clock, a Yule process model for the species tree prior, and other parameters as default. Four independent MCMC analyses were run each for 120 million generations, sampling every 1000^th generation. Another MCMC analysis was run with settings suitable for a subsequent JML analysis (see below), i.e. with piecewise constant population size model and 40 million generations sampling every 1000^th generation. The computer program Tracer version 1.5 [53] was used to check convergence of all parameters to the stationary distribution in each run and TreeAnnotater version 1.7.4 was employed to set the burn-in (discarding the first 30000 trees) and to calculate the maximum clade credibility tree.
4. 6. We performed statistical tree topology tests on ML gene trees inferred from the ITS_1, cpDNA and ITS_cpDNA alignments to find out whether discrepancies between topologies shown in particular gene trees and coalescent species tree are statistically significant. To this end, looking at the topology of the coalescent species tree, we enforced the following constrains on gene trees: basal position of P. olympica within the subsection Hieracioides; sister relationship of P. nuristanica and P. japonica; and monophyly of P. hispidissima and P. hieracioides (Table 1). Constrained trees were built under the evolutionary substitution model as specified for each alignment above, using the maximum likelihood (ML) criterion and heuristic search with TBR swapping algorithm and 10 random sequence addition replicates. The site-wise likelihoods for the best unconstrained ML tree and all constrained trees were calculated in the computer program raxmlGUI version 1.3 [54] and consequently were compared using the approximately unbiased, weighted Shimodaira-Hasegawa, and weighted Kishino-Hasegawa tests as implemented in the computer package CONSEL version 0.1j [55]–[57]. A p-value of <0.05 was chosen for rejection of the null hypothesis that the log likelihoods of the constrained and best unconstrained trees are not significantly different.
5. 7. To test whether hybridization influenced species relationships and could be the source of gene tree incongruence, we employed the program JML [20]. This software uses a posterior distribution of species trees, population sizes and branch lengths to simulate replicate sequence datasets under the coalescent with no migration. The minimum pairwise sequence distance between sequences of two species is evaluated on the simulated datasets and compared to the one estimated from the original data (i.e. from the ITS or cpDNA dataset). This procedure, the posterior predictive test, is a good predictor of hybridization events that disturb the bifurcating species tree model. Two separate JML analyses were run to simulate sequence replicates in the ITS and cpDNA datasets. For these analyses, 40,000 species trees resulting from the JML-specified *BEAST analysis were used. Settings for particular simulations involved: (1) relative mutation rate as inferred from the log file generated during the *BEAST analyses (set to 1.016 and 0.385 for the ITS and the cpDNA simulation, respectively); (2) heredity scalar (2 and 1 for the ITS and cpDNA simulation, respectively); and (3) appropriate model of sequence evolution for both markers. In each analysis, 9,000 trees were removed as burn-in and every 10^th tree was used for simulations. Based on the original sequence data files, minimum pairwise sequence distances between all pairs of species and exact probabilities of observing these distances in simulations under assumption of no migration were calculated. All pairwise sequence distances with p-value <0.05 were recorded as potential cases of hybridization.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Log likelihoods and p-values of AU, WSH, and WKH tests for tree comparisons considering different topological scenarios.

https://doi.org/10.1371/journal.pone.0104929.t001

The ITS region

In addition to homogeneous sequences, we identified numerous accessions containing intra-individual single nucleotide polymorphic sites in the ITS dataset (Table S2). Out of 653 nucleotide positions, 148 were detected to be polymorphic in at least one individual of the subsection Hieracioides, with up to 16 intra-individual polymorphic sites per sequence. Forty two polymorphic sites showed an additive pattern, ranging between 1 and 16 APS per individual sequence (Table S2). The heterogeneous sequences varied conspicuously in terms of number of APS with respect to potentially parental ITS variants, ranging from those with fully additive patterns to those comprising almost completely homogenized contigs.

Results of BI, MP and split decomposition analyses based on the ITS_1 dataset (117 sequences without APS) did not show the analysed Picris taxa from the subsection Hieracioides to form a monophyletic group (Figures 1 and 2A). Specifically, accessions of the presumed outgroup species, P. scaberrima and P. strigosa, were clustered within the ingroup taxa (Figures 1 and 2A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Majority-rule consensus tree of Bayesian inference based on the reduced dataset_ITS_1.

The numbers above branches refer to posterior probability values of Bayesian inference/the bootstrap support as inferred for the maximum parsimony analyses (values ≥50% are shown). Support values for terminal branches are not shown. Dashed lines represent branches collapsed in strict consensus tree of the maximum parsimony analyses. Colour line below the tree indicates affiliation of grouping to the ITS genetic lineages as mentioned in the text. Each accession label includes the population code following Table S1.

https://doi.org/10.1371/journal.pone.0104929.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Neighbour-Net diagram (A) and split support spectrum (B) for the ITS_1 dataset.

The Neighbour-Net diagram is based on uncorrected P-distances. Colour symbols indicate affiliation of accessions to the ITS genetic lineages as mentioned in the text while the shape of symbols indicate affiliation of accessions to taxa. Bootstrap supports of selected important splits are indicated above edges. Column height in the split spectrum represents the number of clade-supporting positions, i.e., putative primary homologies. Column parts above the y-axis represent the ingroup partition while those below the axis correspond to the outgroup partition.

https://doi.org/10.1371/journal.pone.0104929.g002

The Bayesian 50% majority-rule consensus tree and the strict consensus MP tree based on the ITS_1 dataset (117 sequences without APS) displayed essentially identical topologies with only small differences in terminal positions and/or in clade supports (Figure 1). All ingroup taxa including both aforementioned presumed outgroup species formed a large strongly supported clade (BS = 90%, PP = 1.00), with shallow hierarchical structure comprising numerous sub-clades with various levels of support (Figure 1). All non-hieracioides species, P. hispidissima, P. japonica, P. olympica, P. nuristanica, P. scaberrima, and P. strigosa, formed their own strongly supported sub-clades, with P. olympica and P. scaberrima being the most divergent (Figure 1). Picris hieracioides is not monophyletic and forms four ITS groups with allopatric or parapatric spatial distributions (Figures 1, 2A and 3): ITS_1 composed of two sub-clades: ITS_1A – mostly lowland populations morphologically corresponding to the P. h. subsp. hieracioides morphotype from the Apennine Peninsula, the Balkan Peninsula, western Turkey, rarely from Central and northwestern Europe and a single population morphologically assignable to the P. h. subsp. umbellata morphotype from northwestern Europe (blue in Figures 1, 2A and 3); ITS_1B – mostly lowland populations of the P. h. subsp. hieracioides morphotype from Central, northwestern Europe and Balkan Peninsula and a single population of the P. h. subsp. umbellata morphotype from the Western Carpathians (pale blue in Figures 1, 2A and 3); ITS_2 – mountain populations of the P. h. subsp. umbellata morphotype from Alps and the Western Carpathians (red in Figures 1, 2A and 3); ITS_3 – mostly mountain populations corresponding to the P. h. subsp. umbellata morphotype from the Sierra Nevada, Pyrenees, Jura Mts., Alps, Germany and Sweden (orange in Figures 1, 2A and 3); and ITS_4 – an unresolved group of populations assignable to both of P. hieracioides morphotypes from the Apennine and the Iberian Peninsulas, Belgium and Poland (brown in Figures 1, 2A and 3). We also endeavoured to analyse the dataset ITS_2 (all 216 ITS sequences including those with APS) using MP and BI. Both approaches, however, resulted in a completely unresolved clade with all accessions from the subsection Hieracioides placed in a basal polytomy (data not shown).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Map displaying geographic distribution of genetic variation in the analysed Picris taxa, as inferred from the nuclear ITS sequence data.

Coloured symbols indicate affiliation of accessions to the ITS genetic lineages as mentioned in the text. Black symbols indicate populations including individuals with unhomogenised sequences containing additive polymorphic sites. Symbol shapes indicate taxa as follows: P. hieracioides – circle, P. hispidissima – triangle, and P. olympica – square. The circle symbols with empty centres refer to populations morphologically corresponding to P. hieracioides subsp. hieracioides, while those with solid symbols refer to populations morphologically corresponding to P. hieracioides subsp. umbellata. The populations of P. nuristanica, P. japonica and outgroup taxa are not shown. Numbers refer to population codes as denoted in Table S1.

https://doi.org/10.1371/journal.pone.0104929.g003

The split decomposition diagram based on the ITS_1 dataset (117 ITS sequences) showed four main splits having various levels of bootstrap support (Figure 2A): (1) P. olympica separated from the rest by a set of long parallel and strongly supported edges; (2) lowland populations of P. hieracioides, morphologically corresponding to P. h. subsp. hieracioides, are clearly separated from all ingroup and outgroup taxa (ITS 1A and 1B). The exception represented few populations of the P. h. subsp. hieracioides morphotype that clustered together with populations attributable to the P. h. subsp. umbellata morphotype in an unresolved ITS4 group; (3) a large heterogeneous group associating several smaller well supported splits corresponding to P. hispidissima, P. japonica, P. nuristanica, and two lineages composed of mountainous populations of the P. h. subsp. umbellata morphotype (ITS_2 and ITS_3); and (4) a group of both morphotypes of P. hieracioides (ITS_4) and P. strigosa, a presumably outgroup species. Picris scaberrima formed a long edge linked with outgroup taxa by short parallel splits. Importantly, the second neighbour-net analysis based on the ITS_2 dataset (206 ITS sequences) revealed essentially the same pattern as the previous one (figure not shown). This apparently indicates that heterogeneous sequences did not significantly affect the overall split pattern. Individuals with un-homogenized ITS sequences and numerous APS appeared mostly in the central position of the network, while those possessing more homogenized sequences, with lower number of APS, were preferentially incorporated in terminal positions (figure not shown).

Analyses of spectrum of supporting positions revealed that majority of phylogenetically informative positions are shared by outgroup taxa (Figure 2B). The entire functional ingroup, involving all accessions of the subsection Hieraciodes plus P. strigosa and P. scaberrima, are supported only by 7 nucleotide positions. Within the functional ingroup, the highest support obtained P. olympica (9 positions), P. japonica (4 positions), the P. hieracioides clade with lowland populations corresponding to the P. h. hieracioides morphotype (3 positions) and P. nuristanica (2 positions). All other above mentioned splits/subclades were supported at maximum by a single position.

The cpDNA spacer

MP and BI analyses yielded trees with different topologies to a certain extent. However, there was no strongly supported incongruence between them. In order to maintain consistency in the presentation of our results we discuss only the topology based on the BI phylogenetic tree.

In contrast to the ITS trees, P. scaberrima and P. strigosa appeared clustered together with other outgroup taxa. All accessions of the Hieracioides group appeared in a largely unresolved clade with five essentially supported cpDNA subclades and two unsupported groups (cpDNA_A–G) placed in a basal polytomy (Figure 4). The topology and variation of non-hieracioides taxa in the cpDNA spacer tree were in some aspects congruent and in other discordant with respect to the ITS tree (see Figure 1). Although not forming distinct subclades, Picris japonica and P. nuristanica possessed their own exclusive haplotypes, similarly as they did in the ITS data (see Figure 5). Contrastingly, P. olympica and P. hispidissima were heterogeneous and shared haplotypes with P. hieracioides (Figures 4, 5 and 6). Repeatedly, P. hieracioides was heterogeneous and was found in seven genetic cpDNA geographically correlated groups (Figures 4, 5 and 6): unresolved group cpDNA_A – comprising populations of P. hieracioides morphologically corresponding to both subspecies and originating from the entire distribution area, P. hispidissima from Croatia and Montenegro, and P. olympica (cyan in Figures 4, 5 and 6); cpDNA_B – mountain populations of the P. h. subsp. umbellata morphotype from the Apennines, Alps, and Carpathians, while populations from lowland of northwestern Europe corresponding morphologically rather to P. h. subsp. hieracioides and finally P. hispidissima from Croatia (violet in Figures 4, 5 and 6); cpDNA_C – P. hispidissima and lowland populations of the P. h. subsp. hieracioides morphotype from Croatia and Montenegro (olive in Figures 4, 5 and 6); unresolved group cpDNA_D – including mountain and lowland populations of P. hieracioides from Sierra Nevada, the Carpathians and north-western Europe morphologically corresponding to its both subspecies, P. japonica and P. nuristanica (pink in Figures 4, 5 and 6); cpDNA_E – P. hispidissima from Montenegro (grey in Figures 4, 5 and 6); cpDNA_F – populations of both morphotypes of P. hieracioides from the Iberian mountain ranges, the western Alps, the Jura Mts, and northwestern Europe (green in Figures 4, 5 and 6); cpDNA_G – lowland populations of the P. h. subsp. hieracioides morphotype from the Apennine and the Balkan Peninsula, mountain populations of the P. h. subsp. umbellata morphotype from the Alps and P. hispidissima from Croatia (light brown in Figures 4, 5 and 6). The highest density of rare haplotypes was observed on the Balkan Peninsula (nine haplotypes), the Iberian Peninsula mountains (seven haplotypes), and the Apennine Peninsula (five haplotypes) (Figure 6).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Bayesian majority-rule consensus tree based on the cpDNA intergenic spacer dataset.

The numbers above branches refer to posterior probability values of Bayesian inference/the bootstrap support as inferred for the maximum parsimony analyses (values ≥50% are shown). Support values for terminal branches are not shown. Dashed lines represent branches collapsed in strict consensus tree of the maximum parsimony analyses. Colour line below the tree indicates affiliation of groupings to the cpDNA genetic lineages as mentioned in the text. Each accession label includes the population code following Table S1.

https://doi.org/10.1371/journal.pone.0104929.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Maximum parsimony network of the cpDNA haplotypes of Picris populations.

The symbol sizes are proportional to the haplotype frequencies, the lines represent mutational steps, and black dots are unsampled haplotypes. Colour symbols indicate affiliation of accessions to the cpDNA genetic lineages as mentioned in the text.

https://doi.org/10.1371/journal.pone.0104929.g005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 6. Map displaying geographic distribution of genetic variation in the analysed Picris taxa, as inferred from the cpDNA intergenic spacer dataset.

Coloured symbols indicate affiliation of accessions to the cpDNA genetic lineages as mentioned in the text. Bi-coloured symbols denote populations possessing two different haplotypes. Symbol shapes indicate taxa as follows: P. hieracioides – circle, P. hispidissima – triangle, and P. olympica – square. The circle symbols with empty centres refer to populations morphologically corresponding to P. hieracioides subsp. hieracioides, while those with solid symbols refer to populations morphologically corresponding to P. hieracioides subsp. umbellata. Population of P. nuristanica and P. japonica and outgroup taxa are not shown. Numbers refer to population codes as denoted in Table S1. The symbols and their numbers highlighted in red indicate populations harbouring rare haplotypes.

https://doi.org/10.1371/journal.pone.0104929.g006

The parsimony-based haplotype network revealed 34 very closely related ingroup haplotypes, mostly connected by one or two steps (Figure 5) and the split decomposition diagram displayed a stellar-like structure with numerous haplotypes in an unresolved basal position (data not shown). Similarly as in the ITS dataset, the majority of supporting positions was harboured by outgroup species (data not shown). Both hieracioides clades as well as all their sub-clades were supported extremely weakly, in particular, by two positions at the maximum (data not shown).

Concatenated phylogeny and coalescence based species tree

Although the ILD test revealed significant (p<0.001) incongruence between both gene trees, we performed concatenated analyses because our split spectrum analyses show that the observed discordances are not strongly supported. Thus, incongruence between both genetic regions could rather reflect paucity of phylogenetically informative signal (Figures 1, 2AB and 4). Concatenated phylogeny also did not unambiguously resolve evolutionary relationships among ingroup taxa, but still depicted several variably supported (sub)clades. Position of two presumed outgroup taxa, P. scaberrima and P. strigosa, varied among phylogenetic analyses and remained uncertain (see Figure 7AB and MP/BI phylogenies not shown). Non-hieracioides ingroup species, namely P. japonica, P. olympica and P. nuristanica, each clustered in strongly supported sub-clades, repeatedly with P. olympica being the most divergent (Figures 7AB). Importantly, although weakly supported, the majority of accessions assigned to morphotypes of both subspecies of P. hieracioides as well as to P. hispidissima tended to form taxon specific groupings (Figure 7A).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 7. Neighbour-Net diagram (A) and split support spectrum (B) for the concatenated ITS_cpDNA dataset.

Coloured symbols indicate affiliation of accessions to taxa. Bootstrap supports of selected important splits are indicated above edges. Column height in the split spectrum represents the number of clade-supporting positions, i.e., putative primary homologies. Column parts above the y-axis represent the ingroup partition while those below the axis correspond to the outgroup partition.

https://doi.org/10.1371/journal.pone.0104929.g007

The coalescence-based species tree (Figure 8) showed different topology than the concatenated gene tree and brought more light into the phylogenetic relationships of the analysed taxa: the presumed outgroup taxa, P. scaberrima and P. strigosa, appeared together in a sister position to the subsection Hieracioides; P. olympica appeared in a basal position to all taxa from the subsection Hieracioides; the Asian P. japonica and P. nuristanica were clustered together in a sister position to the P. hispidissima and P. hieracioides clade. Relationships among P. hispidissima and P. hieracioides were, however, unsupported.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 8. The maximum clade credibility species tree obtained from Bayesian inference in *Beast based on the ITS and cpDNA loci.

The posterior probability values are indicated above the branches.

https://doi.org/10.1371/journal.pone.0104929.g008

Statistical tree topology tests

Although all tested evolutionary scenarios found in the coalescent species tree (see above and Material and Methods) were not observed in the present gene trees, they cannot be rejected at the significance level of 0.05 by any of the three statistical tree topology tests performed (Table 1). The only exception is the monophyly of P. hispidissima and P. hieracioides, which can be excluded for the ITS_1 and cpDNA alignment at the significance level of 0.05 only by the AU test but not by the WSH and WKH tests. Thus, the phylogenetic relationships within the subsection Hieracioides as unravelled by the Bayesian coalescent-based analysis cannot be very likely excluded also for the gene trees.

Test of hybridization and/or introgression

Only a single case of potential hybridization within the ITS dataset was detected, in between a pair of outgroup taxa, Helminthotheca aculeata and H. echioides (p<0.001). By contrast, 31 cases with p<0.05 were found in the cpDNA dataset involving accessions of the corresponding pairs of the following taxa: P. olympica vs. P. hieracioides, P. olympica vs. P. rhagodioloides; P. sinuata vs. P.capuligera; and H. aculeata vs. H. echioides (Table 3).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. List of the distances among couples of taxa with p-values <0.05, as inferred by posterior predictive distributions.

https://doi.org/10.1371/journal.pone.0104929.t003

Intricate phylogenetic relationships within the subsect. Hieracioides

Tree-building and network analyses of both the nuclear and plastid markers of the studied Picris taxa revealed phylogenies with rather low resolutions and topological discrepancies (see Results), although utilization of these markers has been proposed as one of the most informative approaches even at lower taxonomic levels [35], [58]. More light on the evolutionary relationships within the subsection Hieracioides was cast by the Bayesian coalescence analysis which generates species instead of gene trees. Although relationships in the present species tree were not strongly statistically supported and were not detectable in all gene trees, they could not be rejected by statistical tree topology tests performed on the concatenated dataset (see Results). Thus, in spite of the complex patterns, we might assume the following: (1) the analysed members of the subsection Hieracioides have a monophyletic origin, with P. scaberrima and P. strigosa being their sister taxa; (2) P. olympica is basal within the subsection Hieracioides, i.e. it is sister to all other members of this subsection; (3) all Asian species are sister to the European ones; and (4) P. japonica + P. nuristanica are sister to P. hieracioides + P. hispidissima. Interrelationships among populations of P. hieracioides and P. hispidissima are left unresolved, however.

Which factor(s) might be responsible for low resolution, varying supports and discordances between gene tree topologies? In the present study, we have shown by the split spectrum analysis that one of the primary and crucial factors responsible for this pattern is a paucity of parsimony informative and node supporting nucleotide positions. Deficiency of phylogenetically informative signal in the datasets analysed was unambiguously evidenced by the low number of primary nucleotide homologies (Figures 2B and 7B). This was indicated also by the presence of short internal branches in all phylogenetic trees (Figures 1 and 4) and the high number of rare plastid haplotypes in derived positions connected mostly with ancestral haplotypes by only one or two mutational steps (Figure 5). Low level of phylogenetically informative positions in the used markers, which are among the most informative [35], [58], strongly indicates a recent diversification of the studied taxa and genetic groups [1], [59], [60], [61], [62]. It would be desirable to confirm this by molecular dating. However, in case of insufficient amount of phylogenetic information, such an analysis is assumed to be vague and statistically unsupported.

Furthermore, incomplete lineage sorting represents another evolutionary force that might be responsible for the discordances and complex genetic patterns in the six Picris taxa [1], [3], [4], [5]. The Bayesian coalescent-based approach considers the incomplete lineage sorting as a main source of uncertainties in inferred species trees [18]. Thus, the unresolved relationships among genetic lineages of P. hieracioides and P. hispidissima as well as the low nodal support could signalise either influence of this phenomenon or violation of prior assumption by hybridisation or other processes. Likewise, the disjunctive occurrence of the ancestral plastid haplotypes in populations of P. olympica and the spatially distant P. hieracioides might be the consequence of incomplete lineage sorting. By contrast, the mosaic pattern detected in the ITS dataset may repeatedly suggest influence of hybridization.

Hybridisation might substantially influence phylogenetic reconstructions at even lower hierarchical levels. Decreased phylogenetic tree resolution, incongruent patterns in different gene trees, and occurrence of alleles transcending taxa or genetic groups are commonly observed in plant phylogenies in which extensive gene exchange between taxa has been documented [2], [3], [17], [41], [60], [61], [62], [63], [64], [65], [66]. Presence of extensive and recurrent gene exchange among investigated Picris taxa might be indicated by the following facts. (1) Considerable intra-individual and intra-population ITS sequence heterogeneity and allele transcending clades and taxon borders in P. hieracioides, P. hispidissima, and P. olympica. On the other hand, maintenance of two or more ITS variants within a single genome might be attributed to other molecular-genetic processes [11], [67]. (2) Occurrences of heterogeneous ITS sequences predominantly in the contact zones of the studied species or genetic groups, while locations of homogeneous sequences with little or no genetic admixture mostly along the sampled area's margins inhabited mainly by single lineages or species (see Results, Figure 3). (3) Strict self-incompatibility in P. hieracioides and P. hispidissima, which has been recently documented ([27]; M. Slovák unpublished data). (4) Last but not least, successful crossing between populations of morphotypes and genetic lineages of P. hieracioides as well as with P. hispidissima obtained during our preliminary field experiments (M. Slovák, unpublished data).

Posterior predictive checking revealed hybridization/introgression at least in five pairs of the studied species. Four pairs were recorded exclusively by simulations on the cpDNA dataset and one pair was detected by simulations on both the plastid and the nuclear dataset. Since the method of posterior predictive checking is powerful unless the hybridization event has occurred very rapidly after the speciation event [9], it might be hypothesised that hybridizations among the Picris taxa are rather recent or have occurred later after speciation events. This hypothesis corresponds well mostly in the case of comparatively distantly related species, P. olympica with P. hieracioides, P. olympica with P. rhagodioloides, and P. capuligera with P. sinuata. Such hybridisation events might take part during the Pleistocene glacial/interglacial range shifts, especially, at sites where lowland taxa could have had recurrent contact with mountainous ones [2], [68]. Posterior predictive checking, however, did not support hybridisation scenarios among the pairs of closest relative taxa where it would be expected: among subspecies and genetic lineages of P. hieracioides as well as between P. hieracioides and P. hispidissima. One of the most plausible explanations might be that the hybridization events among those taxa happened just after their diversification. Especially in case of the two subspecies of P. hieracioides, the divergence of genetic lineages might be shallow causing hybridization to be undetectable [9].

The overall failure to detect hybridization events by simulations on the ITS dataset might be attributed to the presence of recombination and/or concerted evolution within the ITS locus. As stated in [9], some patterns resulting from concerted evolution can potentially bias the estimates of population sizes toward lower values. This causes the loci to coalesce faster and therefore hybridization might be harder to detect. Only one hybridization event was found by both genetic markers, viz., in H. aculeata and H. echioides. The sample size might be the limitation factor in the posterior predictive checking approach. Since only a single individual per each aforementioned species (outgroup) was sampled, their population sizes cannot be properly estimated for species tree reconstruction. Therefore we assume this result might be a methodological artefact.

Finally, it remains arguable whether some of the studied taxa, but especially P. hispidissima, fit into the reticulate evolution scenario or, alternatively, its large genetic heterogeneity arose from genetic erosion. The latter hypothesis might be explained by recurrent introgression with adjacent populations of P. hieracioides after their diversification. We are more inclined to the latter scenario because all taxa are morphologically rather well delimited, inhabit ecologically specific biotopes, and trans-taxon alleles are concentrated predominantly in areas of sympatry or parapatry.

Although methods to potentially distinguish between hybridisation and incomplete lineage sorting have been developed, there is still a high risk of confusing these two phenomena or they cannot be discerned in some cases at all [9], [12], [13]. In the investigated closely related Picris taxa from the subsection Heiracioides, it is highly tricky or even impossible to discern between these two processes. Reasons for this ambiguity might come from weak genetic divergence of analysed taxa, and importantly from presence of gene transfer via introgression which, moreover, most probably occurred repeatedly during their evolutionary history.

Comparison of genetic patterns with results inferred from AFLP and morphology

Spatial distribution of genetic variation in populations of P. hieracioides and P. hispidissima varied considerably across phylogenetic trees, with both concordant and discordant patterns being detected (Figures 1 and 4). The ITS data were highly consistent with taxonomic/morphological delimitation of the studied taxa. The genetic ITS groups of all studied species as well as within P. hieracioides exhibited high levels of concordance with the results of the AFLP analyses [24]. The congruence of genetic patterns resolved by ITS and AFLP, especially, with regard to the phylogenetic signal, has been already demonstrated for various plant groups [69]. The spatial distribution of genetic variation in populations and genetic lineages of P. hieracioides and in P. hispidissima was, however, clearly blurred by populations harbouring individuals with non-homogenized ITS sequences. Such a diffuse pattern of genetic variation disrupting geographic structure is not unexpected in the nuclear biparentally inherited genetic markers. In fact, this reflects both seed and pollen-mediated gene dispersal [70], especially, in widespread taxa with unrestricted gene flow [71]–[73]. Numerous heterogeneous ITS sequences did not, however, display APS for some positions variable in the potentially parental haplotypes (Table S2). This pattern can be attributed to in vivo recombination between ITS ribotypes, or to the eventual combination of minor, undetected alleles with more widespread ones [11], [67], [74]. Detection of such minor ITS ribotypes within taxa/populations might be prevented by preferential amplification of the major ITS copies or by quantitative masking of the rare ribotypes [74]. Another highly plausible explanation on the origin of incomplete homogenization of APS involves concerted evolution [2], [40], [58], [63], [64], [67]. Its presence in our ITS dataset might be one of the crucial reasons responsible for the failure to detect mutual hybridisation among populations and genetic lineages of P. hieracioides and P. hispidissima (see discussion above).

In contrast, genetic variation inferred from the plastid cpDNA spacer was, however, only partly congruent with that captured by ITS, AFLPs, and morphology. Further, numerous discrepancies in the number and composition of the genetic cpDNA groups were noted (see ‘Results’). Thus, the plastid pattern exhibited an apparently mosaic-like spatial structure and two or more species/genetic groups shared same plastid haplotypes. This was, especially, true for the sympatric/parapatric populations (Figure 6). Complex and incongruent genetic plastid patterns might, beside of the aforementioned hybridisation and incomplete lineage sorting, arise from different attributes of the employed genetic markers such as their evolutionary rate, inheritance mode (bi-parental vs. uni-parental), or number of copies (multiple vs. single copies), which reflects different time horizons of evolution [11]. The nuclear ITS marker evolves faster than plastid ones and therefore enable detection of more recent diversification events [11], [75].

Evolution of the studied Picris taxa and taxonomic implications

Our results suggested that evolution of the studied Picris species was shaped by the interplay of several factors and evolutionary processes. Their diversification was most likely rapid, and may have occurred in an allopatric manner. Since all of these taxa, but especially, P. hispidissima and P. olympica, have unique habitat and ecology preferences (see Material and Methods), we presume that ecological factors also played an important role in their speciation.

The centre of the diversity of the studied species complex, as indicated by the distributions of P. olympica and P. hispidissima as well as by the accumulation of rare haplotypes, lies in the southern European peninsulas and partly also the Alps and the Carpathians (Figure 6). These regions have been shown to be the most important glacial refugia for the majority of European flora and fauna [76], [77]. This might suggest that the evolutionary history of the investigated Picris taxa was connected with glacial/interglacial cycles in the Quaternary. This is in concordance with conclusions concerning the origin of the AFLP lineages, as detected within P. hieracioides [24].

Genetic data presented herein together with recent distribution [21] indicate that P. olympica very likely evolved in alpine levels of high mountain ranges of Anatolia in Asia Minor. Here it very likely survived isolated for a long period, although remote secondary contacts and interaction with other closely related taxa cannot be excluded. Likewise, autecological characteristics and distribution of P. hispidissima suggest that evolutionary history of this taxon has undoubtedly been confined with coastal, calcareous mountain ranges of the western part of the Balkan Peninsula. This region has also been repeatedly proven as important evolutionary hotspot and glacial refugium for other plant groups [32], [68].

Furthermore, genetic lineages detected within P. hieracioides most likely diversified in allopatry or parapatry in southern European peninsulas and adjacent mountainous regions of the Alps and the Carpathians. The mountain-dwelling populations morphologically corresponding to the P. h. subsp. umbellata most likely evolved in European mountains or their proximity (the Alps, central Apennines, the Carpathians, Iberian mountain ranges and the Jura Mts.) [24]. If it had existed there already during Pleistocene climatic oscillations, it could potentially survive in the peripheral glacial areas together with other mesophylous herbs [78], [79]. As proved by the highest cpDNA haplotype richness, diversification centre of populations morphologically assignable to the lowland P. h. subsp. hieracioides was located more southwards, especially, in lowland and hilly regions of the Apennine Peninsula and plausibly also in the Balkan Peninsula. Subsequently, they colonised other habitats of their present day occurrence in central and north-western Europe as well as in Asia Minor from these southern refugia.

An important taxonomic question is whether the analysed taxa should be recognised as distinct species/subspecies, or whether they rather represent informal phylogeographic lineages within a single heterogeneous taxon. The latter possibility has been favoured in case of several other complex plant species groups [30], [32], [80], [81]. Picris hieracioides has been here evidenced to be highly heterogeneous across both markers. Actually, both plastid and ITS nuclear alleles have transcended borders of both subspecies of P. hieracioides, as recently delimited in [24]. Moreover, several populations of P. hieracioides displayed transitional morphological variation. In addition, from ecological point of view, P. hieracioides is characterised by an exceptionally large range of habitat preferences [22], [24]. All these facts might indicate that both currently delimited subspecies of P. hieracioides are artificial biological units, representing rather ecotypes evolved recurrently in different evolutionary lineages under similar environmental influences than real taxa.

Ecotype hypothesis, however, strongly contradicts results of our previous study focused not only on investigation of wild populations of P. hieracioides, but also on cultivation experiments realised under environmentally homogeneous conditions [22], [23], [24]. These experiments included more than 300 individuals from 32 populations that belonged to both subspecies of P. hieracioides, originating from almost the entire studied area [24]. Subsequent comparative morphological investigations revealed that both subspecies retained morphologically stable under homogeneous environmental conditions. Importantly, these analyses showed that environmental conditions have had only a minute impact on the morphological characters identified as taxonomically diagnostic. Likewise, P. hispidissima retained its morphological distinctiveness under environmentally homogenous conditions (M. Slovák, unpublished data). Thus, taxonomically important morphological features of these entities are not a consequence of environmental plasticity, but are genetically inherited instead. All three taxa are well defined morphologically, ecologically and partially also genetically (AFLP's) at least within ‘pure’ populations untouched by their recent mutual secondary contacts.

Both, subspecies of P. hieracioides as well as P. hispidissima are recently diversified and thus genetically weakly separated taxa. Additionally, their genetic variation has been significantly affected by their reciprocal multiple secondary contacts, taking place after anthropogenically mediated range expansions. Due to the lack of habitat and breeding isolation barriers, hybridisation and introgression among these taxa occur frequently in contact zones. This process leads to the genetic and morphological homogenisation of affected populations, erasing their distinctness in contact zones. Similar cases of taxon fusion caused by hybridisation preceded by migration of populations enhanced by human activities have been reported for instance in members of the genus Knautia [80] and Cerasus [82]. For the latter case the term “anthropohybridisation” was coined. Establishment of sound taxonomic concept in such difficult cases is very problematic and usually it is a matter of subjective choice. Therefore, at the present state of knowledge, we prefer to leave both subspecies of P. hieracioides at currently accepted taxonomic ranks [24]. In contrast, P. hispidissima should be included under highly variable P. hieracioides as separated subspecies: P. h. subsp. hispidissima (Bartl.) Slovák and Kučera, comb. nova, hoc loco (basionym: Crepis hispidissima Bartl. in Bartling & Wendland, Beitr. Bot. 2: 125, 1825).

Although this study does not provide comprehensive and firm taxonomic conclusions, it represents an important background for future investigations on this intricate and dynamically evolving species complex. Last but not least, the present investigation documents rapid process of taxon erosion, happening due to the anthropogenic disturbance of natural biotopes.