Our taxon sampling scheme covered all major lineages of the subtribe Scorzonerinae as recognized in recent molecular phylogenetic studies (Hatami et al., 2020; Zaika et al., 2020), and had the aim to include a wide range of species of Scorzonera in its former widest sense (Zaika et al., 2020). Consequently, representatives of all described genera based on the taxonomic treatment of Zaika et al. (2020) were included: Epilasia, Gelasia Cass., Geropogon, Koelpina, Lipschitzia Zaika et al., Pseudopodospermum (Lipsch. & Krasch.) Kuth., Pterachaenia, Ramaliella Zaik et al., Scorzonera, Takhtajaniantha Nazarova, and Tragopogon. However, for Tourneuxia variifolia Coss., the only member of Tourneuxia, and Lipschitzia divaricata (Turcz.) Zaika et al., the only member of that genus, wet lab treatment did not yield sufficient reads for analysis. A total number of 152 samples from representatives of subtribe Scorzonerinae were included (Appendix). The outgroup to the Scorzonerinae included eight representatives of selected subtribes of different phylogenetic distances according to Kilian et al. (2009a) and Tremetsberger et al. (2012): Chondrillinae [Willemetia stipitata (Jacq.) Dalla Torre, Chondrilla ramosissima Sm.], Cichoriinae (Cichorium intybus L.), Hypochaeridinae [Hypochaeris achyrophorus L., Leontodon tingitanus Ball., Urospermum dalechampii (L.) F. W. Schmidt], Lactucinae (Lactuca sativa L.), and Scolyminae (Scolymus hispanicus L., Catananche arenaria Coss. & Durieu; Appendix).

The DNA was isolated in most cases directly from herbarium specimens of the selected taxa, or silica-dried leaf material vouchered by a corresponding herbarium specimen (12 samples). Herbarium specimens with the permission for DNA extraction for this study were kindly provided by the following herbaria: Botanic Garden and Botanical Museum Berlin (B), Ferdowsi University of Mashhad (FUMH), Ernst-Moritz-Arndt-Universität Greifswald (GFW), Komarov Botanical Institute of Russian Academy of Science St. Petersburg (LE), Staatliche Naturwissenschaftliche Sammlungen Bayerns, München (M), Shahid Bahonar University of Kerman (MIR), Ludwig-Maximilians-Universität München (MSB), Moscow State University (MW), and Naturhistorisches Museum Wien (W). Appendix also includes the specimen data, links to digitized specimens, as well as the accession numbers of the INSDC (International Nucleotide Sequence Database Collaboration) for the deposited sequence data.

Wet laboratory work was undertaken in the Botanic Garden and Botanical Museum Berlin molecular lab and the Berlin Center for Genomics in Biodiversity Research (BeGenDiv) consortium genomics lab. Genomic DNA was extracted from dried leaf tissues, using the NucleoSpin Plant II kit (Macherey-Nagel GmbH, Düren, Germany) and DNeasy Plant Mini Kit (Qiagen, Hilden, Germany) following the manufacturer’s protocol. The total genomic DNA quantity was measured with a Qubit 2.0 Fluorometer (Life Technologies, Grand Island, New York, United States). The quality of genomic DNA extractions was assessed for level of fragmentation and fragment size using a 0.9% (w/v) agarose gel. A total of 1 μg of genomic DNA in 60 μL was sheared to a target average fragment size of ∼500 bp by sonicating for 55 s using a Covaris S220 (Covaris, Brighton, United Kingdom). Sonication was not carried out for well-fragmented (<600 bp) genomic DNA samples extracted from herbarium specimens.

The DNA libraries were prepared using the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, United States), following the standard protocol provided by the manufacturer. We followed the library preparation wet laboratory method described for the Berlin lab, provided by Jones et al. (2019: Appendix S2). However, in some cases, we used 15, 19, or 20 cycles for PCR amplification. Hybrid capture was performed using MyBaits (Arbor Biosciences, Ann Arbor, MI, United States) and the myBaits COS Compositae 1Kv1 (Mandel et al., 2014), according to the wet laboratory methods described in detail by Jones et al. (2019: Appendix S2), but either 19 or 22 PCR cycles on the last amplification step were done. In preparation for the hybrid capture reactions, usually, around 10 samples were pooled. For one sequencing run (including 76 out of the total of 163 samples), the post-capture reaction was spiked with a pre-capture reaction (ratio 3:1) before sequencing, to yield a higher number of off-target sequences. Further quality checking and sequencing were carried out either at Macrogen Inc. (South Korea) on a Hiseq X platform (300 cycles) in paired-end, high-output mode or the BeGenDiv (Berlin, Germany) on an Illumina NextSeq platform (300 cycles) in paired-end, mid-output mode.

Table 1 and Figure 2 summarize the datasets assembled and their corresponding analyses. All data processing and analyses were done on the high-performance computing system of the Freie Universität Berlin (Bennett et al., 2020). Forward and reverse raw reads of each sample were subject to adapter trimming, quality filtering, and duplicate removal, using the initial data cleaning step of the HybPhyloMaker (Fér and Schmickl, 2018) pipeline, which makes use of Trimmomatic v.0.32 (Bolger et al., 2014) and FastUniq v.1.1 (Xu et al., 2012). HybPiper (Johnson et al., 2016) was used to individually map the cleaned and deduped sequences of each sample to the reference sequences of Carthamus tinctorius, Helianthus annuus, and L. sativa for the 1,061 loci of the COS by Mandel et al. (2014) and assemble them into contigs. To achieve this, HybPiper first searches the reads against and sorts them according to the target sequences (the individual COS loci) using BWA (Li and Durbin, 2009). The appropriate target then guides the read assembly into contigs, using SPAdes (Bankevich et al., 2012). After assembly, the SPAdes contigs are aligned to the targets, scaffolded, and translated. Finally, HybPiper extracts the sequences recovered for the same target of all samples of a giving sampling and generates an unaligned multi-FASTA file for each target. HybPiper flags loci with a paralog warning when coding sequences of more than 85% of the reference length are detected in multiple contigs. Among the competing contigs, the one is selected that has a coverage depth exceeding the other by 10×, or else the one with the greatest percent identity to the reference. Paralog flagging is particularly useful when frequent reticulation events can be suspected, but has the drawback that even simple allelic variation may trigger such warnings. Facing a high percentage of target loci with paralog warnings, we, therefore, assessed the loci with paralog warnings using the script provided by HybPiper. In this way, the inferred paralogous sequences of all loci of the sampling were collected along with all locus sequences without paralogs. The matrices for the individual loci were then aligned using MAFFT v.7.266 (Katoh and Standley, 2013) and gene trees were calculated using FastTree (Price et al., 2010). By visual inspection of these gene trees, loci possessing homologous sequences that were resolved in a sister-group relationship for typically <3 species were considered to indicate allelic variation, selected, and added back to the original loci set (OLS) without paralog warnings. To examine the effect of the additional loci gained through this step on the results of the phylogenetic analysis, the enlarged loci set (ELS) was analyzed and used separately from the OLS. We also used the HybPiper pipeline to assemble the COS exons with flanking introns (the “splash zone” of Weitemier et al., 2014). All loci with paralog warnings were removed from the matrix and the resulting dataset is hereon referred to as COS supercontigs (Table 1 and Figure 2).

Following the approach described by Herrando-Moraira et al. (2019), we assembled plastid genome sequences from the off-target reads using the MIRA v.4.0.2 multi-pass DNA sequence data assembler/mapper (Chevreux et al., 1999) with the MITObim v.1.9 wrappers originally developed for de novo assembly of mitochondrial genomes (Hahn et al., 2013). In a three-step process, the reads are mapped by MIRA to a phylogenetically related reference genome (L. sativa DQ_383816.1) to identify the more conserved regions between the total readpool and this initial reference and then assembled into contigs, yielding a new and gapped reference sequence. In the second step, the MITObim wrapper script uses the gapped reference to fish in the readpool for partly or fully overlapping reads. In the third step, the overlapping reads are mapped to the gapped reference sequence and incorporated, using MIRA again. The two last steps are iteratively repeated until stationarity of the mapped and assembled reads are reached (Hahn et al., 2013). The assembled plastid genomes were annotated with the web application GeSeq (Tillich et al., 2017), using Cichorium intybus (NC_043842.1) and L. sativa (NC_007578.1) as references and, after reorganizing the sequence designation in the resulting GeSeq multi-fasta files to start with the sample designation, the Phyluce package (Faircloth, 2015) script “phyluce_assembly_explode_get_fastas_file” was applied for separating the individual coding regions (CDS) into individual gene files. The resulting dataset is hereon referred to as plastome CDS (Table 1 and Figure 2).

All sequence matrices [COS contigs (OLS and ELS), COS supercontigs, plastome CDS; Table 1 and Figure 2] were aligned individually for each locus or gene, respectively, using MAFFT v.7.266 (Katoh and Standley, 2013) with the parameter –auto. Phyutility vers. 2.2 (Smith and Dunn, 2008) was used to trim alignments by removing sites with a threshold of 50% missing data. Supercontig alignments may represent misassembled contigs, therefore, trimAl v.14 (Capella-Gutiérrez et al., 2009) was used to remove spurious sequences and poorly aligned regions applying the parameters -resoverlap 0.65, -seqoverlap 70, and the gappyout (settings tested in Jones et al., 2019). The same procedure but with the automated1 function was applied to plastome CDS alignments. To avoid obstruction of later tree calculation, locus or gene designations were removed from the sample names in the alignments of the supercontigs and plastid DNA matrices. The tool AMAS (Borowiec, 2016) was then used to retrieve summary statistics for the alignments. Based on these summaries, loci represented in less than 50% of the samples in COS contigs and COS supercontigs alignments were removed. In the plastome CDS alignment 124 samples were represented, for which plastid genomes were successfully assembled. Before the analyses, the plastome CDS alignment was subject to visual inspection and manual cleaning from apparently mis-assembled sequence portions and loci with less than 50% taxa, and less than 10 parsimony-informative (PI) sites were removed from the primary 89 loci dataset, and subsequently, 71 loci were retrieved (plastome CDS71 dataset). For the Bayesian analysis a second, smaller matrix was created from which all genes with less than 50 PI sites (proportion ≤ 0.03%) were removed (plastome CDS26 dataset).

Phylogenetic reconstructions were conducted for three types of datasets: COS contigs (OLS, ELS), COS supercontigs, and plastome CDS (Table 1 and Figure 2). The sampling of both the COS contigs (160 samples) and plastome CDS datasets (124 samples) was based on selecting the representatives of subtribe Scorzonerinae as ingroup and outside of this subtribe as outgroup. To avoid homology assessment problems in the aligning of the flanking intron sequences as would be expected for remote outgroup members in the COS supercontigs analyses, the early diverging Scorzonerinae genus Gelasia was selected as the outgroup for this dataset, and a total of 124 samples across the remainder of the subtribe were included as ingroup (Appendix).

The three types of datasets were analyzed separately under two principal approaches: (1) the multispecies coalescent approach (Liu et al., 2019), in which the species tree is estimated from the individual gene trees resulting from phylogenetic analyses of each locus or gene, and (2) the concatenation approach, using a supermatrix of the concatenated locus or gene alignments for tree inference with ML or Bayesian analyses.

In the multispecies coalescent approach, we first calculated individual gene trees with RAxML-NG using ParGenes (Morel et al., 2019), as a tool for parallel model testing and tree inference of numerous individual loci or genes. ParGenes uses Modeltest-NG (Darriba et al., 2019) to calculate the scores for the best model; we selected BIC as the model test criterium. Parameters for RAxML-NG were set as follows: 50 ML tree searches using 25 parsimony-starting-trees + 25 random-starting-trees; bootstrapping either with the autoMRE option in effect (disables multinodal parallelization of runs), a requested maximum of 3,000 replicates and using the default threshold of 0.03 for the BS convergence assessment (COS contigs and supercontigs datasets), or without autoMRE option and 1,200 replicates requested (plastome CDS dataset); mapping the BS support values onto the best-scoring/best-known ML tree. Recent studies have questioned the traditional assumption of the plastome as a uniform and single locus, indicating that it is rather a mosaic of genes evolving under different constraints (Gonçalves et al., 2019; Walker et al., 2019). Therefore, we also applied the multispecies coalescent approach to the plastome CDS dataset to compare phylogenetic inference with the concatenation approach. Following the recommendations by Mirarab (2019), Newick Utilities (Junier and Zdobnov, 2010) were used before the species’ tree calculation with ASTRAL vers. 5.6.3 (Mirarab and Warnow, 2015), to collapse nodes with very low support (less than 10%), and TreeShrink v. 1.3.1 (Mai and Mirarab, 2018) with standard parameters was used to detect outliers with abnormally long branches and to remove such samples from individual gene trees. ASTRAL was then used to generate the species tree by maximizing the number of quartet trees shared between gene trees and the species tree (Mirarab and Warnow, 2015) and calculating local posterior probabilities (LPP; Sayyari and Mirarab, 2016) as branch support values. The nodes with LPP < 0.5 were considered as not statistically supported and were collapsed in all coalescent species trees. The levels of discordance between the individual nuclear gene trees were assessed with the program PhyParts (Smith et al., 2015). The ELS ASTRAL species tree was used as a reference tree. PhyParts requires rooted trees with the same outgroup. Therefore, 190 individual gene trees and the species tree of the ELS dataset were rooted with C. arenaria as an outgroup using R and the package APE (Paradis and Schliep, 2019) after 19 gene trees missing C. arenaria were excluded. With the script phypartspiecharts.py (Johnson, 2017), the output of PhyParts was visualized by plotting pie charts on the reference tree that show the proportions of concordant and discordant gene trees for each bipartition.

In the concatenation approach, the alignments of each type of dataset were combined into a supermatrix using AMAS (Borowiec, 2016). ML analyses were run using the Multi-Point Interface (MPI) version of RAxML-NG v. 0.8.1 and 0.9.0 (Kozlov et al., 2019). To assess the significance of the partitioning and choice of DNA substitution models for our datasets, two analyses were run: one with non-partitioned datasets and the other with partitioned datasets. The general GTR+G model of sequence evolution was applied to the non-partitioned concatenated datasets. For the partitioned dataset, the best partitioning schemes and substitution models were obtained using PartitionFinder v.2 (Lanfear et al., 2017) with the relaxed clustering algorithm criterion as recommended by Lanfear et al. (2014) for large phylogenomic datasets with the “rcluster” search option (–rcluster-max 100 and –rcluster-percent 0.1), “BIC” model selection parameters, linked branch lengths, and a choice between three substitution models (GTR, GTR+G, and GTR+I+G). Using RAxML-NG, the tree space was explored with 50 (for the plastome CDS71 and COS contigs datasets) or 20 (COS supercontigs datasets) ML tree searches using 25 or 10 random and 25 or 10 parsimony-based starting trees, respectively, followed by standard bootstrapping, which employed the bootstopping test with a maximum replicate number and a bootstrap (BS) convergence requirement with 3% default cutoff for each dataset (Table 1). The BS support values were mapped onto the best-scoring ML tree obtained and the nodes with BS of <50% were collapsed in all ML concatenated trees.

Because the phylogenetic reconstructions of both the gene trees (as input for the multispecies coalescent analyses) and the species trees based on the concatenated matrix were conducted with ML using RAxML-NG, and tree inference with RAxML-NG can occasionally be misled by non-randomly distributed missing data (Xi et al., 2015), we also applied a Bayesian analysis on the concatenated matrices to test the robustness of the reconstructions. We used PhyloBayes (Lartillot et al., 2013; Lartillot, 2020), which implements a non-parametric approach based on Dirichlet process priors to model nucleotide or amino acid substitutions as site-specific random variables directly inferred from the data, as opposed to being specified a priori. This model classifies amino acid or nucleotide sites into topological categories (therefore “CAT model”). We applied the default CAT-GTR and a discretized gamma distribution with four categories, which are computationally resource-intensive but have shown to perform well in inferring accurate branching patterns in genomic datasets (Whelan and Halanych, 2017; Lartillot et al., 2018). Convergence of the chains can become challenging for larger alignments beyond 20,000 sites (Lartillot et al., 2018), therefore, six chains were independently run for >25,000 cycles, and their stationarity, appropriate burn-in, and convergence were first visually assessed using Tracer vers. 1.7.1 (Rambaut et al., 2018). Subsequently, convergence in tree space and reproducibility of the posterior consensus trees across chains was assessed using the bpcomp and tracecomp commands of PhyloBayes. Chains with an insufficient sampling of parameters [low ESS (=effective sampling size)] and chains stuck in a local optimum (lower loglikelihood values) and, thus, badly converging with the others were discarded. A minimum of two well-converging chains with the highest loglikelihood values was used to calculate the final posterior consensus trees.

The COS contigs and plastome CDS species trees were rooted with the members of the Scolyminae clade (S. hispanicus and C. arenaria) because this is the earliest diverging clade of the tribe in our sampling (Kilian et al., 2009a; Tremetsberger et al., 2012). The COS supercontig tree with its more restricted sampling was rooted with the Gelasia clade (Appendix). TreeGraph vers. 2 (Stöver and Müller, 2010), was used for viewing and displaying the phylogenetic trees. Before, format conversion in FigTree vers. 1.4.4. (Rambaut, 2018) was necessary in some cases.

Addressing the unusual diversity of pollen types in the subtribe as an exemplar, we wanted to test if the resolved phylogenetic relationships help to explain the current distribution of morphological character states in the subtribe. To reconstruct the pollen types at ancestral nodes, we used the backbone of the ELS ASTRAL tree as a phylogenetic hypothesis and treated the seven pollen types identified by Blackmore (1982) as unordered states. The matrix was built and the parsimony ancestral characters state reconstruction was done with Mesquite vers. 3.7 (Maddison and Maddison, 2014; Mesquite Project Team, 2014).

In total, 163 samples were analyzed in our study, including both nuclear and plastome datasets, and for 157 of them, Hyb-Seq data were newly created (Appendix). The DNA was isolated in most cases directly from herbarium specimens of the following age ranges: 3 collected from ≤ 1900, 12 from > 1900 ≤ 1950, 65 from > 1950 ≤ 2000 and 77 (12 of which from silica-dried material) from > 2000.

The average number of reads per sample obtained was 4,952,940, ranging from 147,379 in Scorzonera laciniata DB 44268 to 25,313,374 in H. achyrophorus DB 546. On average, the percentage of mapped reads per sample was 35% (range: 7–58%) and the number of COS loci recovered for each sample ranged from 423 in S. laciniata DB 44268 to 1,041 in Pseudopodospermum picridioides DB 44329, with a mean of 968 loci of the total of 1,061 target loci in the set. A total of 897 COS loci and 832 COS supercontig loci were identified by HybPiper as loci containing potentially paralog sequences and removed from the datasets before analyses. After removing the spurious sequences of the matrices in the COS supercontigs dataset in trimAl and eliminating the loci with less than 50% taxa in both datasets, the OLS dataset and COS supercontigs datasets yielded 147 and 148 loci, respectively. The final concatenated OLS dataset had a length of 44,634 bp, of which 16,315 were parsimony informative (∼35%) and an average of 8% missing data per locus (range: 0–27.5%; Table 1). The final concatenated sequences of the COS supercontigs dataset had a length of 132,646 bp, of which 65,023 were parsimony informative (∼48%) and an average of 7.5% missing data per locus (range: 3.9–11.1%; Table 1). Visual inspection of the 897 gene trees of loci with paralog warnings, suggested that the multiple sequences of 70 loci represented allelic variants, not paralogs. Apart from six loci containing less than 50% taxa, these loci were added to the OLS dataset, forming the enlarged ELS dataset of 211 loci. The final concatenated sequences of the ELS dataset had a length of 64,358 bp, consisting of 23,379 parsimony informative sites (∼35%) and an average of 8.6% missing data per locus (range: 0–27.4%; Table 1). The statistics for the OLS, ELS, and supercontigs datasets before (Supplementary Table 1) and after (Table 1) removing the loci with less than 50% taxa revealed a slight increase in the percentage of variable and PI sites in the condensed datasets, but with no considerable impact on the average percent of missing data per locus. Reducing the number of loci in the plastome dataset from the original 89 to 71 by eliminating those with less than 50% taxa and less than 10 PI sites (plastome CDS71 dataset) caused a reduction in the percentage of variable and PI sites, as the well as the length of concatenated alignments and average sequence length per locus (Supplementary Table 1 vs. Table 1). The concatenated plastome CDS71 sequence was 63,917 bp long, containing 6,321 PI sites (∼8%), and an average of 25.9% missing data per locus (6.6–47.4%). The plastome CDS26 dataset built by removal of all genes with less than 50 PI sites had 26 genes, the concatenated sequence had a length of 48,478 bp, containing 5,300 PI sites (∼11%), and an average of 25% missing data per locus (10–38.8%; Table 1, plastome CDS26 dataset).

In the current study, fourteen phylogenetic analyses based on the nuclear (OLS, ELS, and supercontigs) and plastome CDS datasets were performed under ML and Bayesian inference based on concatenated loci and under the multispecies coalescent model based on the individual gene trees (Table 1; Figures 2, 5–7; and Supplementary Figures 1–5). Trees estimated under the multispecies coalescent model resulted in a lower mean of support values and percentages of nodes with maximum support than those obtained under the concatenation approaches (Table 2). No hard topological incongruences exist between the species trees of both principal approaches (coalescent and concatenation), apart from the single contrary consecutive sister group relationship of the monospecific Scorzonera renzii and S. rupicola clades (Figure 3). Collapsing nodes with very low support and removing outlier samples with abnormally long branches from the gene trees using TreeShrink, before running coalescent analyses, did not cause any notable difference in topology and support values, compared to trees without employing TreeShrink (not shown).

In general, deep tree branches were well-supported and the backbone topology was largely congruent in all inferences of each dataset (Figures 3, 4), while the internal resolution of the lineages was comparatively poor for all datasets (Figures 5–7). However, some small and large groups with congruent topologies and strong support values were also found within major lineages across the trees, and reciprocally, the corresponding tips of some well-supported subclades in shallow nodes generated most cases of topological incongruences between different analyses (Figures 5, 7).

Technical details regarding phylogenetic analyses conducted on the COS contigs (OLS and ELS), COS supercontigs, and plastome CDS datasets are given in Table 1. Hereon the support values are provided in parenthesis in the following order: ML analysis on the unpartitioned concatenated dataset [bootstrap percentage (BS-NP)], ML analysis of partitioned concatenated dataset [bootstrap percentage (BS-P)], Bayesian inference [posterior probability (PP)], and coalescent-based approach (LPP).

Overall, the phylogenetic inferences in this study based on the nuclear on-target and plastome off-target Hyb-Seq data revealed the same major clades as those inferred by Zaika et al. (2020). Clade designations, therefore, follow Zaika et al. (2020: Figures 1, 2).

The exemplary ancestral character state reconstruction for the seven pollen types on the nuclear backbone tree is provided in Figure 4C and Supplementary Data Sheet 1. The distribution of the pollen types in the lineages of the Scorzonerinae is characterized by multiple state dimorphisms and the reconstruction at ancestral notes shows many homoplasies.

The phylogenetic backbones of the four analyses based on the plastome CDS matrix were largely congruent with each other (Figures 4, 7 and Supplementary Figure 5). Analyses of the plastome dataset resolved the same major lineages as in the analyses of the nuclear (COS contigs and COS supercontigs) datasets, but their relationships showed several conspicuous differences (Figures 3, 4). This concerns the relationships of Gelasia, Geropogon, Takhtajaniantha, the S. rupicola, and S. renzii clades, and the clade comprising Koelpinia, Pterachaenia, and Ramaliella.

Subtribe Scorzonerinae was resolved as monophyletic with full support. Within the subtribe, Geropogon was resolved as sister to the remainder of the subtribe, also with full support in all analyses; Tragopogon was not represented in the plastome datasets. ML and Bayesian concatenated analyses inferred a sister relationship between a clade containing Takhtajaniantha as sister to Scorzonera (90 BS-NP, 92 BS-P, 0.97 PP) and another clade including Epilasia, Gelasia, Koelpinia, Pseudopodospermum, Pterachaenia, and Ramaliella (99 BS-NP, 99 BS-P, 1 PP), while in the coalescent analysis, the relationship between Takhtajaniantha and the Scorzonera lineage was not statistically supported and the clade, therefore, collapsed, resulting in a trichotomy (Figures 4A,B). Koelpinia, Ramaliella, and Pterachaenia formed a clade with full support in all analyses, while their internal relationship was not statistically supported by the ML trees and incongruent in the other two (Figures 4A,B). Gelasia, Epilasia, and Pseudopodospermum formed a second clade, which was supported only in the Bayesian analysis (0.99 PP), with Gelasia as sister to Epilasia and Pseudopodospermum; in the ML and Bayesian trees, the sister group relationship between Epilasia and Pseudopodospermum was resolved (98 BS-NP, 98 BS-P, 1 PP), while in the ASTRAL tree, this relationship was not statistically supported (Figures 4A,B). Within the Scorzonera lineage, all inferences indicated similar backbone relationships among major clades. A sister-group relationship between the S. angustifolia and S. albicaulis clades was supported in all analyses (100 BS-NP, 100 BS-P, 1 PP, 0.99 LPP). The relationship of the Scorzonera s.typ. clade, as sister to the clade including the Podospermum, S. purpurea, and S. renzii clades, was resolved in both the Bayesian and ML trees (99 BS-NP, 96 BS-P, 1 PP) but not statistically supported in the ASTRAL tree, and the relationship of the S. renzii clade, as sister to the clade formed by Podospermum and the S. purpurea clades was resolved in all analyses (99 BS-NP, 99 BS-P, 1 P, 0.96 LPP). However, the sister relationship of Podospermum and the S. purpurea clade was recovered with low support (71 BS-NP, 68 BS-P, 0.55 P, 0.94 LPP; Figures 4A,B).

Similar to the results for the nuclear datasets, congruence in the phylogenetic backbones of the different analyses of the plastome datasets contrasts with frequent cases of topological incongruences and differences in statistical support values in shallow nodes (Figure 7 and Supplementary Figure 5). A closer consideration follows in the section “Discussion.”

Our study presents the first phylogenomic analysis of the Scorzonerinae based on the myBaits COS Compositae 1Kv1 (Mandel et al., 2014). Most of the data were generated from herbarium specimens; 38 (24%) of the samples newly sequenced were 50 years or older, which demonstrates the power of Hyb-Seq in the valorization of the wealth of material preserved in herbaria around the world.

In principle, species trees estimated under the multispecies coalescent model and trees from the concatenation approach were congruent. However, there were lower means of support values and percentages of nodes with maximum support in coalescent model trees compared with those from the concatenated approach. Removal of outlier samples with abnormally long branches from the gene trees with TreeShrink did not cause any notable difference in topology and support values of the coalescent model trees. This indicates that our datasets were largely free from significant biases caused by contamination, mistaken orthology, and misalignment, which are common reasons for long branches (Mai and Mirarab, 2018). Low support values and statistically unsupported nodes in species trees of the multispecies coalescent approach are due to a lack of support in individual gene trees or gene tree discordance. The PhyParts output (Figure 4A) demonstrates massive discordance among the gene trees in the Scorzonerinae in that all backbone nodes were supported only by a minority of gene trees and the vast majority of shallower nodes have similar levels of discordance (Supplementary Figure 6). In contrast, the root of the Scorzonerinae was found in 60% of the gene trees, markedly corroborating them as a separate lineage. Important sources of gene tree discordance are hybridization and incomplete lineage sorting (Liu et al., 2015) because of ancient rapid radiation (Whitfield and Lockhart, 2007; Oliver, 2013). Trees based on the multispecies coalescent approach better reflect the phylogeny than such based on concatenated matrices, because the multispecies coalescent model takes incomplete lineage sorting into account in that it allows conflicting phylogenetic signals among individual loci and provides a more realistic measure of support (Liu et al., 2015, 2019; Sayyari and Mirarab, 2016). We, therefore, conclude that other biological causes such as reticulate or convergent evolution (Liu et al., 2019) remain the main causes of gene tree discordance in the Scorzonerinae.

Among the analyses of concatenated datasets, partitioning led to a slightly higher percentage of nodes with full support and typically a higher mean of the support values, compared to the trees estimated from unpartitioned datasets applying the generic nucleotide substitution model GTRGAMMA (Table 2). Log-likelihood scores were also higher for ML trees with partitioning (Table 1). Our findings, thus, confirmed the positive impact of data partitioning on tree resolution and support values, concordant with previous studies (Xi et al., 2012; Kainer and Lanfear, 2015; Warnow, 2015; Lanfear et al., 2017; Jones et al., 2019), although the impact is not striking in our case.

Adding flanking introns of the COS exons to generate supercontig datasets increased the alignment length nearly threefold, the number of variable and PI sites more than three- and four-fold, respectively, and correspondingly, the percentage of PI sites increased (c. 37%). This was expected because the non-coding intron regions possess a great amount of variation and are less sensitive to selective pressures, such as selection-driven convergence (Weitemier et al., 2014; Folk et al., 2015; Johnson et al., 2016; Gernandt et al., 2018; Kates et al., 2018). The resulting increase in the mean of support values and percentage of nodes with full support in terms of LPP in ASTRAL species trees and BP support levels in ML concatenated analyses, compared to OLS trees, was, nevertheless, altogether moderate but more visible in the species trees of the multispecies coalescent analyses (Table 2). The corresponding full trees (Figures 5–7; Supplementary Figures 1–5) show increased support for branches at deeper nodes and otherwise in particular an improved resolution and support for branches at shallow nodes, leading, in several cases, to a substantial improvement in the reconstruction of species relationships. This agrees with previous studies that highlighted the potential of supercontigs for improved statistical support (Jones et al., 2019; Bagley et al., 2020; Gardner et al., 2021). Signal saturation at most of the deeper nodes and the opposing impact of discordant phylogenetic signals due to reticulate evolution at deeper and, in particular, at shallow nodes may explain the rather moderate impact of supercontigs on support values.

Of the 1,061 COS loci, a significant number were flagged as potentially paralogous by HybPiper (897 COS contigs and 832 COS supercontigs). Jones et al. (2019) discussed the relationship between whole-genome duplications (WGDs) in the evolution of a lineage and the number of paralogous loci. WGDs have played a major role also in the diversification of the Compositae, and the Cichorieae are among the tribes, which are known to have experienced WGDs (Huang et al., 2016). Jones et al. (2019) found 721 loci flagged as paralogous in a small sampling across the tribe, which was the highest value among all Compositae tribes included. At the same time, they reported distinctly smaller numbers for individual subsets of this sampling (Jones et al., 2019: Table 1), confirming that a wider systematic range of a sampling potentially increases the number of putative paralogs. The Vernonieae, for which no WGDs, so far, have been identified, had, however, the second-highest number of putative paralogs, and Jones et al. (2019) assumed that this is due to lineage-specific or even species-specific WGDs within the Vernonieae. The Scorzonerinae may be similar in this respect: the number of putative paralogous loci for the Scorzonerinae here is even higher than in Jones et al. (2019). Moreover, in our study, the number of putative paralogs detected was only slightly increased (6.6%) when comparing the Scorzonerinae only vs. Scorzonerinae plus eight samples from more distantly related Compositae subtribes (supercontigs sampling vs. COS contigs sampling). Therefore, the Scorzonerinae alone appear to be a major source of potential paralogs. WGD, as allo- or autopolyploidy is not apparent in the subtribe with an only exceptional occurrence of polyploidy in single species of Koelpinia and Takhtajaniantha with a basic number of x = 6 or 7 (Zaika et al., 2020). However, lineage diversification and evolution of innovations are not directly initiated by WGD but through processes of post-WGD genome rearrangements, including massive loss of duplicated genes, resulting in post-polyploid genome diploidization (Mandákova and Lysak, 2018). Therefore, the occurrence of WGD in the deeper history of a lineage is rarely manifested by persisting polyploidy. We suggest that the high paralogy found within Scorzonerinae may be evidence for extensive reticulate evolution at different evolutionary timescales.

The extension of our original COS dataset to include loci with paralog warnings that were found to be based on allelic variation only, allowed us to investigate the impact of an increase from 147 loci (in the original set, OLS) by 64 loci or c. 43% to 211 loci (in the extended data set ELS) on the phylogenetic reconstruction. The higher number of loci brought about a corresponding increase in variable and PI sites, whereas their percentages remained almost constant (Table 1). ML concatenation analyses of the ELS dataset resulted in a somewhat improved resolution, indicated by the slightly higher percentage of nodes with full support whereas the topologies obtained from both datasets were largely congruent with each other (Table 2 and Figures 3, 4). Also, in the coalescent analyses, the ELS tree shows improved resolution and support, notably in the case of the successive sister group relationships of Pseudopodospermum, Takhtajaniantha, and Scorzonera to the Tragopogon-Geropogon-Epilasia clade.

The first plastid DNA phylogeny of the Scorzonerinae by Zaika et al. (2020: Figure 2) confirmed the major lineages found in their nrITS phylogeny that was also corroborated by our analyses but failed to resolve relationships between them (with a few exceptions). Most of the relationships inferred from our plastome CDS dataset are therefore new results and show that cytonuclear discordance (Lee-Yaw et al., 2019) is a major issue already in the backbone of the species trees (Figure 4, discussed in more detail below).

The success in using the off-target reads for assembling plastid genomes has been shown to depend significantly on the spiking of the post-capture library with the unenriched pre-capture library (Jones et al., 2019). Among the 124 samples, of which we could successfully assemble the plastome and build the plastome CDS matrix, more than half were spiked with the unenriched pre-capture library. Considering the effort to assemble and annotate complete plastomes in relation to the cost of a sequencing job, it seems questionable if spiking is a recommendable strategy compared to sequencing unenriched libraries alone. MITObim v.1.9 wrappers (Hahn et al., 2013) for plastome assembly, which makes use of the MIRA v.4.0.2 multi-pass DNA sequence data assembler/mapper (Chevreux et al., 1999) in an automated iterative process, was encountered as a fast approach but with a severe drawback. Numerous stretches in the individual CDS sequences, where the read coverage was insufficient, were found crudely misassembled and needed manual cleaning. We, therefore, also restricted our use of the plastome to the protein-coding CDS and refrained from attempting to use the entire plastome, thus also including non-coding regions, for the phylogenetic analysis. This restriction must be seen as critical in the light of recently reported gene tree conflicts in phylogenies based on plastome data (e.g., Gonçalves et al., 2019; Walker et al., 2019; Cho et al., 2020; Köhler et al., 2020), which questions the assumption of the plastome as a single locus. Cho et al. (2020) recovered different topologies for deep branches in Sonchus comparing phylogenetic reconstruction based on the protein-coding CDS only, with one based on the entire plastome, and Köhler et al. (2020) even recorded conflicting topologies among major clades of Opuntioideae (Cactaceae), when exploring different assemblies of top-informative CDS markers. However, it cannot in all cases be excluded so far that non-biological causes are responsible for such conflicts. The non-coding portions of the plastome are characterized not only by a higher frequency of substitutions compared to the more conserved coding regions but also by microstructural mutations ranging from length-variable mononucleotide and other hypervariable stretches to inversions, providing a challenge for the correct assessment of the positional homology of the nucleotides. Employing correctly aligned non-coding regions, thus, becomes more mandatory and crucial, the lower the overall genetic distance (Escobari et al., 2021).

It seems at least a precautionary measure to extend the application of the multispecies coalescent approach also to the plastid DNA matrix. However, our reconstructions based on the protein-coding CDS show no hard topological incongruences that would indicate the presence of significantly discordant plastid gene trees, apart from generally lower support of the branches and a consequently lower resolution of the ASTRAL tree (Figures 4A,B, 7, Supplementary Figure 5, and Table 2). This result agrees with the general finding that in the Asteraceae the CDS regions are strongly conserved (e.g., Kim et al., 2005; Pascual-Díaz et al., 2021). Our exclusive use of the protein-coding CDS sequences limited the phylogenetic resolution, such as that shallow nodes, and, thus, relationships at the species level, frequently remained unresolved, apparently because relatively few genetic differences have accumulated among taxa of more recent origin.

Our phylogenetic reconstructions corroborate the monophyly of the Scorzonerinae and fully resolved its backbone phylogeny (Figure 3). Moreover, the taxon composition of the major lineages identified is highly consistent in all our nuclear and plastid DNA analyses (Figures 5–7); whereas the analyses of Sanger sequenced data by Zaika et al. (2020: Figures 1, 2) successfully identified phylogenetic lineages within Scorzonerinae, the relationships among lineages of the subtribe were mostly unresolved. In so far as our study corroborates the major lineages identified by Zaika et al. (2020) and resolved their relationships, it confirms the monophyly of the genera recognized by the authors in their revised generic taxonomy of the subtribe.

Our trees obtained from the nuclear datasets show robust support for the early divergence of two principal clades within the Scorzonerinae, one consisting of Gelasia, Koelpinia, Pterachaenia, and Ramaliella and the other of Pseudopodospermum, Takhtajaninatha, Epilasia, Tragopogon, Geropogon, and Scorzonera (Figure 3). The monospecific NW African Tourneuxia, for which sequences could not be generated, can be expected to be either sister to both clades, or sister to the members of the first clade, according to the findings by Zaika et al. (2020). Both principal clades are also present in NW Africa, which is estimated to belong to the ancestral area of the tribe, of which subtribe Scorzonerinae (Zaika et al., 2020) is one of the early diverging clades (Tremetsberger et al., 2012; Kilian et al., in prep.).

The Scorzonerinae have, however, developed their highest diversity in the E Mediterranean-SW to Middle Asian region. The further diversification of the first principal clade, with Gelasia and Koelpinia as successive sisters to a clade of Pterachaenia and Ramaliella (Figure 3), is centered in the E Mediterranean-SW to Middle Asian region, where the Pterachaenia and Ramaliella clades are exclusively found. The same holds for the second principal clade, where Epilasia, Pseudopodospermum, Takhtajaniantha, and Scorzonera are successive sisters to a Tragopogon-Geropogon clade (Figure 3). Epilasia and Takhtajaniantha are restricted to that area and the other lineages have their highest diversity there. The phylogeny of Tragopogon (Mavrodiev et al., 2004) clearly shows that the N African species evolved from later diversification and migration.

Tremetsberger et al. (2012) estimated the origin of the Scorzonerinae clade in the late Oligocene (c. 25 mya) and the diversification between Gelasia (represented in their study by G. hirsuta) and Scorzonera and Tragopogon, corresponding to the origin of the two principal clades, in the late Early Miocene (c. 17–18 mya). This age estimation for the diversification of the subtribe was corroborated by Fernández-Mazuecos et al. (2016) and coincides with the onset of a series of tectonic, orographic, and climatic changes in the Miocene. These led to an expansion of open vegetation and the formation of diversified mountain habitats and triggered diversification also in other subtribes (Tremetsberger et al., 2012; Kilian et al., 2017).

Scorzonera and Tragopogon are estimated by Tremetsberger et al. (2012) to have split in the transition between the Middle and Late Miocene (c. 11 mya). However, Bell et al. (2012), who analyzed diversification and diversification rates of Tragopogon with its >150 species, estimated the origin of that genus later in the Late Miocene c. 7.4 mya (3.7–11.6 HPD). The highest species diversity of Tragopogon is centered in the mountainous habitats around the Paratethys basin, where the Black, Caspian, and Aral Seas represent relics of the former Paratethys Sea. They estimated its diversification to have taken place rather late and rapidly from 2.6 mya onward, thus, well after the Messinian Salinity Crisis in the late Miocene between 5.3 and 5.96 mya, which brought about a dramatic aridization of the entire region (Bell et al., 2012). Regarding the distribution of Tragopogon, parallels can be drawn to some extent also for the other large genera of the subtribe (Gelasia, Pseudopodospermum, and Scorzonera). Furthermore, a rather recent diversification seems likely given the shallow genetic differences within these major lineages (see below).

The evolutionary history of the subtribe, however, must have been somewhat more complicated, considering the backbone tree based on the uniparentally inherited plastome CDS (Figures 4A,B). Lacking a sample of Tourneuxia, which was inferred as the first diverging clade of the subtribe in the plastid DNA tree by Zaika et al. (2020), our analyses resolved Geropogon, represented by two samples, as sister to the remainder of the subtribe (Figure 4). We assume that the Tragopogon-Geropogon clade is another early diverging clade of the subtribe in the plastome phylogeny because in both nuclear and plastid DNA phylogenies (Mavrodiev et al., 2012; Zaika et al., 2020) Geropogon is the well-supported sister group to Tragopogon, for which we could not assemble a plastome. The position of Geropogon in our plastid DNA tree is, therefore, in striking contrast to the topology of the nuclear phylogeny with the Tragopogon-Geropogon clade deeply nested in the principal clade with Scorzonera (Figure 3). We propose two potential explanations for this discrepancy: the plastome phylogeny may present the species phylogeny, encapsulating the early diverging origin of the Geropogon-Tragopogon lineage before it underwent reticulation with one lineage of the remainder of the subtribe. Alternatively, the nuclear phylogeny may present the species phylogeny and the Geropogon-Tragopogon lineage acquired the plastome from an early diverging lineage of the subtribe, e.g., Gelasia, by hybridization with extensive backcrossing, thus through chloroplast capture (Rieseberg and Soltis, 1991). The hybridization left no apparent morphological traces in the Geropogon-Tragopogon lineage due to the absence of significant nuclear gene flow. The clear morphological distinction of the Geropogon-Tragopogon lineage, compared to the remainder of the subtribe, which is evident by its uniseriate involucre, may provide an argument for the former scenario. The pollen type of Tragopogon shared with Koelpinia, and that of Geropogon, shared with a few members of Gelasia (G. lanata type; Blackmore, 1982, 1986) do not help to decide between the scenarios. The fact that Lipschitzia, the monotypic Central-E Asian genus accommodating the former Scorzonera divaricata and not represented in our analyses, is resolved, as sister to the Tragopogon-Geropogon clade in the plastid DNA tree by Zaika et al. (2020), may support the cytoplasmatic gene flow hypothesis. However, that Lipschitzia is not an early diverging lineage in the nuclear phylogeny by Zaika et al. (2020) further complicates the picture.

Within the large clade sister to Geropogon in the plastid DNA tree (Figures 4A,B), two main clades are resolved: one including Takhtajaniantha and Scorzonera, their sister group relationship is, however, not statistically supported in the coalescent analyses, while the other includes the remainder of Scorzonerinae. Zaika et al. (2020), found the same clades with mostly strong statistical support, and also the Koelpinia-Pterachaenia-Ramaliella and the Epilasia-Pseudopodospermum subclades in the second clade. In the nuclear trees, in contrast, Scorzonera, Takhtajaniantha, Pseudopodospermum, and Epilasia form consecutive sister groups to the Geropogon-Tragopogon clade (Figures 4A,B). The Koelpinia-Pterachaenia-Ramaliella clade is remarkable as the only major clade besides Scorzonera, the composition of which is congruent in both phylogenies.

Within the Scorzonera clade, an odd incongruence occurs with the monospecific S. renzii clade: in the plastid DNA tree it groups with the S. purpurea-Podospermum clade (ASTRAL tree) or the S. purpurea-Podospermum-Scorzonera s.str. clade (concatenation trees; also, in Zaika et al., 2020), but with the S. rupicola-S. angustifolia-S. albicaulis clade in the nuclear trees (Figures 4A,B). Therefore, the relationship of S. renzii remains an open question.

Different processes can cause cytonuclear discordance: ancestral polymorphism and population splitting resulting in incomplete lineage sorting, fixation of different organellar genomes from existing variation by the selection, or hybridization with cytoplasmic introgression (chloroplast capture; Lee-Yaw et al., 2019). Considering that the coalescent method rules out incomplete lineage sorting, our findings of massive nuclear gene tree conflicts and the high rate of paralogs may support the hypothesis that hybridization with cytoplasmic and/or nuclear introgression led to multiple events of reticulate evolution during both the early and later diversification within the Scorzonerinae. This is similar but distinctly more pronounced than what was concluded for the subtribe Lactucinae (Kilian et al., 2017).

The shallow morphological differentiation between the major Scorzonerinae lineages may also be of some significance in this context. As already noted by Zaika et al. (2020), morphology badly reflects the major lineages and genera, making their diagnoses burdensome. Most major lineages and clades lack non-homoplastic synapomorphies (Mavrodiev et al., 2004; Hatami et al., 2020; Zaika et al., 2020). Likewise, reliable diagnostic morphological characters seem largely absent for the two principal clades and their subclades. In contrast to the case of the Lactucinae, where similar sets of features evolved multiple times in different lineages (Kilian et al., 2017, in prep.), certain patterns of character state distribution across, more or less, distant clades in the Scorzonerinae may find a more parsimonious explanation through reticulation than convergent evolution. Pollen features (Blackmore, 1982, 1986) may illustrate this thought. The Scorzonerinae are unique among the tribe in this respect. Not only does the pollen possess colpori with only two lacunae in contrast to three in all other subtribes, but their pollen also exhibits seven distinctive morphological types. According to the hypothesis of character evolution provided by Blackmore (1982), they form four main developmental branches. The distribution of these pollen types in the nuclear backbone tree and their parsimony ancestral character reconstruction (Figure 4C) shows a conspicuous mixture of multiple state dimorphisms and dual or plural homoplasies in the terminals. Similar patterns could be elaborated for other morphological character sets and their comparisons could help to infer the parentage of hybrid lineages.

In principle, our study confirmed that the Hyb-Seq approach is useful to resolve close relationships at interspecific levels. However, branch support strongly decreased for internal nodes closer to the tips compared to backbone relationships. We interpret this because of the presumed young age of the speciation within the major lineages, which was inferred for Tragopogon by Bell et al. (2012) and can likely be generalized for the other larger lineages of the subtribe, Scorzonera, Pseudopodospermum, and Gelasia. In addition, we suspect that reticulation events between and within the major lineages have contributed to blurring interspecific relationships, as can be deduced from the comparatively low percentage of loci confirmed as orthologous across the subtribe (see above). Consequently, in many cases, infrageneric and interspecific relationships received no or low statistical support and showed incongruences among analyses, which, considering their lack of support, are not meaningful. Some species groups, in contrast, were consistently well-supported as monophyletic. In the following, we focus on these well-supported species groups, discussing their internal relationships, morphological synapomorphies, and possible occurrence of cytonuclear discordance.

Gelasia: Both nuclear and plastid DNA analyses strongly support the Gelasia lanata group composed of G. flaccida (S. flaccida Rech. f.), G. lanata (L.) Zaika et al. (Figure 1D), and G. psychrophila (Boiss. & Hausskn.) Zaika et al. (BS ≥ 99, 1 PP, 1 LPP), although the relationships between these species remained unresolved (Figures 5–7 and Supplementary Figures 1–5). Morphologically, they possess a combination of the tuberous root, scape-like stems, and achenes with densely white lanate hairs. The Gelasia persica group, including G. persica (S. persica Boiss. & Buhse), G. wendelboi (S. wendelboi Rech. f.), and G. xylobasis (S. xylobasis Rech. f.) is also monophyletic in both nuclear and plastid analyses (BS ≥ 89, PP ≥ 0.95, LPP ≥ 0.9, Figures 5–7, and Supplementary Figures 1–5), except for the OLS and ELS coalescent and supercontig ML-unpartitioned-concatenation analyses. These three species are endemic to Iran and possess lanceolate or oblong leaves with three to five main veins and glabrous achenes. They share this combination of morphological characters with G. tomentosa (L.) Zaika et al. and G. cinerea (Boiss.) Zaika et al., which is consistent with the nuclear DNA analyses that placed the latter two species with the G. persica group in a moderately to the highly supported clade (BS ≥ 98, 1 PP, LPP ≥ 0.98, 73 BS-NP, 84 BS-P, Figures 5, 6 and Supplementary Figures 1–4). However, in the plastid DNA analyses, the G. persica group in its narrow sense formed instead a well-supported clade with G. seidlitzii (Boiss.) Zaika et al. (89 BS-NP, 95 BS-P, 0.95 PP, 1 LPP, Figure 7, and Supplementary Figure 5), thus, represent a clear case of cytonuclear discordance. A larger Gelasia villosa group, besides G. villosa (Scop.) Cass., which provides the type of the generic name, including G. doriae (Degen & Bald.) Zaika et al., G. ensifolia (M.Bieb.) Zaika et al., G. biebersteinii (Lipsch.) Zaika et al., G. caespitosa (Pomel) Zaika et al., and G. callosa (Moris) Zaika et al. were usually resolved with low to moderate support only but fairly consistently in both nuclear (67 BS-P, 55 BS-NP, 1 PP, 0.9 LPP, incompletely so in the supercontig analyses and ML-unpartitioned-concatenation analysis of OLS, Figures 5, 6 and Supplementary Figures 1–4) and plastid DNA analyses (88 BS-P, 59 BS-NP, 0.7 LPP, incompletely so in the Bayesian analysis, Figure 7, and Supplementary Figure 5). We have, however, not noticed morphological synapomorphies for this group.

Pseudopodospermum: All nuclear and plastid DNA analyses resolved the Pseudopodospermum papposum group, including P. ovatum (Trautv.) Zaika et al., P. ferganicum (S. ferganica Krasch.), P. papposum (DC.) Zaika et al., P. picridioides (Boiss.) Hatami, and P. limnophilum (S. limnophila Boiss.), as monophyletic with high support (98 BS, 0.99 PP, 1 LPP, Figures 5–7, and Supplementary Figures 1–5). Monophyly of this group is also supported by morphology, as these species share a tuberous root, tuberculate achenes, and a pappus with entirely plumose bristles (Hatami and Mirtadzadini, 2022). Also, well-supported is the Pseudopodospermum undulatum group composed of P. undulatum (Vahl) Zaika et al., P. brevicaule (Vahl) Zaika et al., P. baeticum (DC.) Zaika et al., and P. reverchonii (Debeaux & Hervier) Zaika et al. (97 BS, 1 PP, 0.98 LPP, Figures 5, 6 and Supplementary Figures 1–4). Their internal relationships remain largely unresolved. The plastid DNA analyses resolved P. reverchonii as sister to P. baeticum (100 BS, 0.99 PP, 1 LPP, Figure 7, and Supplementary Figure 5) and P. brevicaule as sister to P. undulatum (100 BS, 1 PP, 0.98 LPP, Figure 7, and Supplementary Figure 5), but both clades are part of a larger polytomy, and their relationship is, thus, unresolved. The distribution of the species group is centered in the W Mediterranean, with P. undulatum extending through N Africa to the Arabian Peninsula, a geographic pattern that is exclusive to this species group in Pseudopodospermum. A third well-supported group in the nuclear DNA analyses is that of Pseudopodospermum incisum (DC.) Zaika et al., which also includes, apart from the name-giving species, P. calyculatum (Boiss.) Zaika et al., P. bicolor (S. bicolor Freyn & Sint.), and P. violaceum (D. F. Chamb.) Zaika et al. (BS ≥ 0.96, 1 PP, 1 LPP, Figures 5, 6 and Supplementary Figures 1–4); the plastid DNA analyses did not resolve this clade. The four species share dentate to pinnatisect leaves as exclusive synapomorphy in contrast to the entire leaves in all other species of Pseudopodospermum.

Takhtajaniantha: All nuclear analyses strongly support a Takhtajaniantha austriaca group composed of T. crispa (M. Bieb.) Zaika et al., T. grubovii (S. grubovii Lipsch.) and T. austriaca (Willd.) Zaika et al. (100 BS, 1 PP, 1 LPP, Figures 5, 6 and Supplementary Figures 1–4), their internal relationships, however, were not well-supported. Also, well-supported in all nuclear analyses is the sister group relationship of T. mongolica (Maxim.) Zaika et al. and T. pseudodivaricata (Lipsch.) Zaika et al. (100 BS, 1 PP, 0.99 LPP, Figures 5, 6 and Supplementary Figures 1–4), having cauline stem leaves with the opposite arrangement as a synapomorphy. In all nuclear analyses, except the supercontig ones, the T. austriaca group clade is sister to a clade including the remainder of the Takhtajaniantha taxa of our analysis, thus, in addition to the latter two species, these are T. capito (Maxim.) Zaika et al., T. ikonnikovii (Krasch. & Lipsch.) Zaika et al., T. pusilla (Pall.) Nazarova, and T. tau-saghyz (Lipsch. & G. G. Bosse) Zaika et al.

Scorzonera s. typ. clade: Among the four species of this clade, all nuclear (except the supercontig coalescent analysis) and plastid analyses support the sister group relationship of S. radiata Fisch. ex Ledeb. and S. humilis L. well (100 BS, 1 PP, LPP ≥ 0.99, Figures 5–7, and Supplementary Figures 1–5). Morphologically, this correlates with the red spot on the apex of phyllaries in both species. Scorzonera aristata stands also morphologically apart from tuberculate-ribbed achenes instead of the otherwise common smooth-ribbed achenes; S. parviflora is odd as a halophytic species with creeping rootstock.

Scorzonera albicaulis clade: All nuclear analyses support a Scorzonera bracteosa group well, consisting of S. bracteosa C. Winkl., S. crassicaulis Rech. f., and S. tragopogonoides Regel & Schmalh. (BS ≥ 90, 1 PP, 1 LPP, Figures 5, 6 and Supplementary Figures 1–4). OLS Bayesian and ML-concatenation analyses of ELS also resolved the interspecific relationships of this group, with S. bracteosa as sister to S. crassicaulis and both as sister to S. tragopogonoides, while in the plastid analyses S. tragopogonoides was not included. The three species possess very large capitula (5–7 cm long at fruiting) and bracteal leaves below the capitula.

Podospermum clade: In all nuclear and plastid analyses, the species relationships within the Podospermum clade are largely unresolved. The only congruent interspecific topology in all nuclear analyses and both plastid ML-concatenation analyses is the sister group relationship of S. meshhedensis (Rech. f.) Rech. f. and S. songorica (Kar. & Kir.) Lipsch. & Vassilcz. (100 BS, 1 PP, 1 LPP, Figures 5–7, and Supplementary Figures 1–5), already found in previous studies (Hatami et al., 2020; Zaika et al., 2020). Both species share leaf heteromorphy (entire and pinnatisect) and an easily detachable pappus.

All in all, the Hyb-Seq approach helped to clarify infrageneric and interspecific relationships in the Scorzonerinae but only in a few cases. Infrageneric classification of the major lineages, thus, remains a challenge, both as such and concerning a suitable methodological approach.

New Combinations

The following new combinations have become necessary for species confirmed by our analyses as members of phylogenetic lineages that are classified as separate genera of the Scorzonerinae following Zaika et al. (2020):

Revised infrageneric classification of Scorzonera

Maintaining the circumscription of Scorzonera by Zaika et al. (2020), which is confirmed by our reconstruction, we propose the following revised infrageneric classification.

1. Scorzonera sect. Rupicolae E. Hatami, N. Kilian & K.E. Jones, sect. nov. – Type: S. rupicola Hausskn.

Diagnosis. Suffruticose, dense pulvinate-caespitose perennial; stem with few or reduced leaves; capitula few on a single stem, easily falling; achenes attenuated at the apex.

Note. S. rupicola was resolved already by Zaika et al. (2020) in a clade of its own. This suffruticose perennial has reduced leaves, easily falling capitula, and beaked achenes (Figure 1F). The combination of these morphological characters is unique in Scorzonera. Moreover, S. rupicola has no close allies in the genus.

2. Scorzonera sect. Renzianae E. Hatami, N. Kilian & K.E. Jones, sect. nov. – Type: S. renzii Rech. f.

Diagnosis. Uppermost leaves filiform, spreading, subtending the capitula; capitula with short peduncles or subsessile in a racemiform synflorescence; phyllaries thick membraneous, subkeeled by ± prominent midrib.

Note. S. renzii was resolved as the only species of the “S. renzii clade” in Zaika et al. (2020) and our analyses. Morphologically, the combination of racemiform synflorescence, capitula with subtending leaves, and thick membraneous phyllaries with ± prominent subkeeled midrib are exclusive to this species in Scorzonera. Additionally, our analyses show this species to have no closer allies.

3. Scorzonera sect. Piptopogon C. A. Mey. ex Turcz. in Bull. Soc. Imp. Naturalistes Moscou 21(3): 97. 1848 ≡ Scorzonera subg. Piptopogon (C. A. Mey. ex Turcz.) C. Díaz & Blanca in Anales Jard. Bot. Madrid 43: 330. 1987. – Type: Scorzonera macrosperma Turcz. ex DC. [= S. albicaulis]

Diagnosis. Perennial herbs or subshrubs; leaves graminoid; capitula with 4–12 flowers only; achenes without carpopodium, beaked or at least attenuate at apex; pappus dirty yellow, caducous; pollen with 24 lacunae.

Note. All nuclear and plastid analyses support the sister relationship of the S. angustifolia and S. albicaulis clades and the members of these clades also share morphological similarities, particularly having the same pollen type and apically attenuate to beaked achenes. The section has been recognized for a long time. It was also resolved as a separate clade by Zaika et al. (2020) with some 13 species. It is here confirmed to also contain S. angustifolia L.

4. Scorzonera sect. Scorzonera

Diagnosis. Leaves entire; capitula solitary or by a few only; inner phyllaries with apical red or blackish spots; achenes without carpopodium.

Note. The typical section includes only four species: Scorzonera aristata, S. humilis, S. parviflora, and S. radiata (inclusion assumed by Zaika et al. (2020) and confirmed by our study).

5. Scorzonera sect. Podospermum (DC.) Benth. in Bentham & Hooker, Gen. Pl. 2: 532. 1873 ≡ Podospermum DC. in Lamarck & Candolle, Fl. Franç., ed. 3, 4: 61. 1805, nom. cons. ≡ Scorzonera subg. Podospermum (DC.) Lipsch., Fragm. Monogr. Gen. Scorzonera 1: 7. 1935. – Type: Scorzonera laciniata L.

Diagnosis. Leaves (at least some) pinnately divided, more rarely graminoid; phyllaries subapically often corniculate; achene with conspicuous tubular carpopodium; achene surface mostly glabrous or somewhat hairy.

Note. The sister-group relationship of the S. purpurea and Podospermum clades was resolved with high support in all nuclear analyses and received at least weak to moderate support also in the plastid DNA analyses, and the presence of a tubular carpopodium in the achenes is a synapomorphy for the Podospermum-S. purpurea clade. We take this as justification for treating them in a single section. Within this section, the members of the S. purpurea clade are, nevertheless, set apart by their caudex with (instead of lacking) blackish-brown fibrous leaf sheath residues, graminoid leaves (instead of having at least some pinnately divided leaves), and pink or purplish (instead of yellow) flowers. These morphological differences may, however, be taken as justification to treat both clades alternatively as separate sections.