Phylogeny and biogeography of the northern temperate genus Dracocephalum s.l. (Lamiaceae)

The northern temperate genus Dracocephalum consists of approximately 70 species mainly distributed in the steppe‐desert biomes of Central and West Asia and the alpine region of the Qinghai‐Tibetan Plateau (QTP). Previous work has shown that Dracocephalum is not monophyletic and might include Hyssopus and Lallemantia. This study attempts to clarify the phylogenetic relationships, diversification patterns, and the biogeographical history of the three genera (defined as Dracocephalum s.l.). Based on a sampling of 66 taxa comprising more than 80% from extant species of Dracocephalum s.l., morphological, phylogenetic (maximum parsimony, likelihood, and Bayesian inference based on nuclear ITS and ETS, plastid rpl32‐trnL, trnL‐trnF, ycf1, and ycf1‐rps15, and two low‐copy nuclear markers AT3G09060 and AT1G09680), molecular dating, diversification, and ancestral range estimation analyses were carried out. Our results demonstrate that both Hyssopus and Lallemantia are embedded within Dracocephalum and nine well‐supported clades can be recognized within Dracocephalum s.l. Analyses of divergence times suggest that the genus experienced an early rapid radiation during the middle to late Miocene with major lineages diversifying within a relatively narrow timescale. Ancestral area reconstruction analyses indicate that Dracocephalum s.l. originated in Central and West Asia and southern Siberia, and dispersed from Central and West Asia into the QTP and adjacent areas twice independently during the Pliocene. The aridification of the Asian interior possibly promoted the rapid radiation of Dracocephalum within this region, and the uplift of the QTP appears to have triggered the dispersal and recent rapid diversification of the genus in the QTP and adjacent regions. Combining molecular phylogenetic and morphological evidence, a revised infrageneric classification of Dracocephalum s.l. is proposed, which recognizes nine sections within the genus.


Introduction
Located in the centre of the Eurasian continent and representing the largest arid zone in the Northern Hemisphere, the arid Asian interior (here mainly referred to as Central and West Asia) is characterized by a continental climate with low annual precipitation and highly variable seasonal temperatures (Walter, 1979;Lioubimtseva and Cole, 2006). However, despite the harsh hydrothermal conditions, three regions of the Asia interior, i.e. the Caucasus, Irano-Anatolian region, and mountains of Central Asia, are listed among the 35 biodiversity hotspots defined by Marchese (2015). During the past few decades, numerous studies have made major progress towards understanding the aridification processes and moisture patterns of the Asian interior (Miao et al., 2012, and references therein;Hurka et al., 2019). In contrast, the impact of aridification on the diversification process of plants inhabiting the interior Asia has been poorly documented from a phylogenetic and biogeographic perspective. Meanwhile, a close relationship between the floras of Central and West Asia and the Qinghai-Tibetan Plateau (QTP), another global biodiversity hotspot, has been suggested. Patterns of biogeographical connections between the QTP and other regions (including Central and West Asia) were summarized by Wen et al. (2014). Some lineages were inferred to have originated in the QTP and subsequently migrated to the inland Asia (the "out-of-QTP" hypothesis), whereas other lineages were purported to have originated in Central and West Asia and later dispersed to and diversified within the QTP. However, this latter pattern is still poorly understood and requires additional documentation (Wen et al., 2014). The genus Dracocephalum L. (Lamiaceae, Nepetoideae, Mentheae) is an ideal group to investigate the evolution of the Asian interior arid zone flora and the relationship between the floras of Central and West Asia and the QTP, as the genus is most common in the steppe-desert biomes of the Asian interior and the alpine region of the QTP (Budantsev and Lobova, 1997).
Dracocephalum is the second largest genus in subtribe Nepetinae and comprises approximately 70 species widely distributed in steppes, semi-deserts, deserts, and alpine regions of temperate Eurasia, with disjunct endemic species in North America and North Africa (Wu and Wang, 1977;Budantsev, 1987Budantsev, , 1993Li and Hedge, 1994;Harley et al., 2004;Mabberley, 2017). Often known as 'dragonheads', several species of Dracocephalum have been used as medicinal, ornamental or nectariferous plants (Harley et al., 2004;Mabberley, 2017). The genus is quite variable in leaf, inflorescence, and calyx morphology, as well as in habit (Fig. 1). It can be distinguished from other genera within subtribe Nepetinae by having swollen folds at the base of the sinuses between the calyx teeth, the two anther thecae diverging at 180°relative to each other, and by (generally) having bracts with aristate teeth (Fig. 1;Li and Hedge, 1994;Harley et al., 2004).
Relationships between the genera Dracocephalum and Lallemantia Fisch. & C.A. Mey. have been controversial. Lallemantia consists of five species distributed in Southwest and Central Asia (Harley et al., 2004). It differs from Dracocephalum only in having flattened pedicels and two folds within the upper corolla lip (Wu and Wang, 1977;Harley et al., 2004). Most taxonomists have considered the two genera as distinct entities (Bentham, 1848;Briquet, 1895Briquet, -1897Schischkin, 1954;Huang, 1977;Wu and Wang, 1977;Budantsev, 1987;Harley et al., 2004), but Budantsev (1993) treated Lallemantia as a subgenus of Dracocephalum, subg. Lallemantia (Fisch. & C.A. Mey) A.L. Budantsev. He pointed out that species of Lallemantia possess all the features of Dracocephalum, including aristate-dentate bracts, the sinuses between calyx teeth with a thickened fold, the much broader medium lobe of the posterior calyx lip, and similar nutlet morphology. Moreover, the base chromosome number of Lallemantia (x = 7) is congruent with Dracocephalum (x = 5, 6, 7), and both differ from most other genera of Nepetinae (x = 8,9) (Budantsev, 1985;Budantsev et al., 1992). The base chromosome number of 7 was considered as ancestral within Dracocephalum (Budantsev, 1985;Budantsev et al., 1992). Additionally, Drew and Sytsma (2012), in a study focused on tribe Mentheae, showed that the two genera are closely related and Lallemantia might be embedded within  Dracocephalum. Drew and Sytsma (2012) also showed that another genus, Hyssopus L., may also be embedded within Dracocephalum. Consisting of about six species, Hyssopus is a taxonomically complicated genus distributed in Central and West Asia, northwest Africa, and Europe (Li and Hedge, 1994;Harley et al., 2004). The genus is often used for medicinal, ornamental, and aromatic purposes (Huang, 1977;Mabberley, 2017). Hyssopus and Dracocephalum are morphologically distinct; they have long been considered to be distantly related, and were placed in different subtribes (Huang, 1977;Wu and Wang, 1977;Harley et al., 2004) prior to the treatment of Drew and Sytsma (2012). However, given the limited sampling of Dracocephalum (two species), Hyssopus (one species), and Lallemantia (one species) used by Drew and Sytsma (2012), relationships among these three genera need to be further elucidated using broader taxon sampling.
Here, we define Dracocephalum s.l. to include the three closely related genera Dracocephalum s.s., Hyssopus, and Lallemantia. The phylogenetic framework of Dracocephalum s.l. is reconstructed to clarify relationships within the group based on the most comprehensive geographic and taxonomic sampling to date. We then use this framework to explore the biogeographic history and diversification processes within the genus. In addition, the infrageneric classification of Dracocephalum proposed by Budantsev (1987) is evaluated on the basis of the original morphological and molecular phylogenetic evidence. And finally, an updated taxonomic treatment of Dracocephalum is presented.

Taxon sampling and choice of markers
A total of 125 accessions representing 58 species of Dracocephalum were sampled (Table S1), covering most of the geographical distribution and all three subgenera and seven sections as described by Budantsev (1987). We also included five accessions of four species of Hyssopus and nine accessions of four species of Lallemantia. Outgroups were represented by 15 species from seven genera within subtribe Nepetinae based on Drew and Sytsma (2012).
Nuclear ribosomal internal and external transcribed spacers (ITS and ETS) and four chloroplast DNA markers (rpl32-trnL, trnL-trnF, ycf1, and ycf1-rps15) were selected for the phylogenetic reconstruction of Dracocephalum, based on previous phylogenetic studies of tribe Mentheae Sytsma, 2011, 2012). We also used two low-copy nuclear genes from the pentatricopeptide repeat (PPR) gene family, AT3G09060 and AT1G09680. Previous studies (Yuan et al., 2009(Yuan et al., , 2010Drew et al., 2014;Roy and Lindqvist, 2015) have demonstrated that PPR loci have a rapid rate of evolution and can be useful markers for phylogenetic reconstruction at intergeneric and interspecific levels.

DNA extraction and amplification
Total genomic DNA was extracted either from silica-gel-dried leaf material using the modified CTAB method (Doyle and Doyle, 1987) or from herbarium specimens using the DNeasy Plant Pro Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. Polymerase chain reaction (PCR) mixtures and procedures of the two nuclear ribosomal DNA (nrDNA) and four chloroplast DNA (cpDNA) regions followed those described in Chen et al. (2016). The primer pairs 17SE and 26SE (Sun et al., 1994) were used for the amplification of ITS, and ETS-B (Beardsley and Olmstead, 2002) and 18S-IGS (Baldwin and Markos, 1998) for ETS. Plastid markers were amplified using the primers trnL (UAG) and rpL32-F (Shaw et al., 2007) for rpl32-trnL, "c" and "f" (Taberlet et al., 1991) for trnL-trnF, 4597f and 5778r (Chen et al., 2016) for ycf1, and rps15r and 5711f (Drew and Sytsma, 2011) for ycf1-rps15.
A nested PCR procedure described in Drew and Sytsma (2013) was used to amplify the two low-copy nuclear genes AT3G09060 and AT1G09680, with thermal cycle conditions as described in Yuan et al. (2009). The 25 lL PCR mixtures included 15.25 lL ddH 2 O, 3 lL 109 reaction buffer, 2 lL dNTP mixture (2.5 mM), 0.3 lL bovine serum albumin (BSA, 20 mg/mL), 0.5 lL DMSO, 0.8 lL of each primer (10 lM; Sangon Biotech, Shanghai, China), 1.2 lL 59 Q-solution (Qiagen, Hilden, Germany), 0.15 lL Taq polymerase (2.5 U/lL; Tiangen Biotech, Beijing, China), and 1 lL template DNA. The product from the first round PCR was diluted by 109 with ddH 2 O, and further used as the template for the second PCR. For AT3G09060, the primer pairs 930F and 2080R (Yuan et al., 2009) were used in the initial amplification, and bd1000F and bd2040R (Drew and Sytsma, 2013) were used for the nested PCR amplification. For AT1G09680, the primers 180F and 1760R (Yuan et al., 2009) were used for the initial amplification, and bd244F and bd1683R (Drew, 2011) for the nested amplification.
Sequencing reactions were performed with the dideoxy chain termination method, using the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, CA, USA) and sequenced on an ABI 3730xl DNA Analyzer (Applied Biosystems, CA, USA). The 5 lL PCR sequencing reaction mixtures included 1 lL of BigDye mix. The primers bd1000F and bd2040R were used as the sequencing primers for AT3G09060, and bd244F and bd1683R as the sequencing primers for AT1G09680. The sequencing primers for the remaining six DNA markers were the same as the PCR primers described above.
A total of 1107 sequences were newly generated in present study, and 18 sequences from GenBank were also included. Voucher information and GenBank accession numbers for all sequences are listed in Table S1.

Sequence alignment and phylogenetic analyses
Raw sequences were assembled and edited using Sequencher 4.1.4 (Gene Codes, Ann Arbor, MI, USA), and aligned using MUSCLE (Edgar, 2004) as implemented in MEGA 6.0 (Tamura et al., 2013) based on default settings. The resulting alignments were then manually adjusted in MEGA based on the similarity criterion (Simmons, 2004). Ambiguously aligned regions were removed, gaps were treated as missing data, and indels were not coded. To reduce the influence of alignment uncertainty, gaps in the nrDNA and cpDNA sequences were also removed using Gblocks 0.91b (Castresana, 2000;Talavera and Castresana, 2007) with default parameters (including the "Allowed Gap Positions" = "With Half") prior to phylogenetic analyses. maximum parsimony (MP), Bayesian inference (BI), and maximum likelihood (ML) analyses were used for phylogenetic reconstruction.
The MP analyses were conducted in PAUP* 4.0b10 (Swofford, 2002). All characters were treated as unordered and equally weighted. Branches with a minimum possible optimized length of zero were collapsed (Simmons and Freudenstein, 2011). The heuristic search option was used with tree bisection-reconnection (TBR) branch swapping and 1000 random-addition-sequence replicates, saving up to 10 000 trees per replicate. A strict consensus tree was summarized from all the most-parsimonious trees (MPTs) retained. Bootstrap support (BS) values were calculated via bootstrap analyses (Felsenstein, 1985) from 1000 replicates, each comprising 10 random-addition-sequence replicates with a maximum of 100 trees saved per replicate.
Partitioned BI and ML analyses were performed using MrBayes 3.2.7a (Ronquist et al., 2012) and RAxML-HPC2 8.2.12 (Stamatakis, 2014), respectively, on XSEDE through the web server Cyberinfrastructure for Phylogenetic Research Science (CIPRES) Gateway (http://www.phylo.org/; Miller et al., 2010). Prior to the BI analyses, the nucleotide substitution models for each DNA marker were selected under the Akaike Information Criterion (AIC) using jModelTest 2.1.6 (Darriba et al., 2012) on XSEDE on the CIPRES Gateway. The best-fitting models were TIM3 + Γ for AT3G09060, GTR + I + Γ for AT1G09680 and ITS, GTR + Γ for ETS and ycf1-rps15, and TVM + Γ for rpl32-trnL, trnL-trnF, and ycf1. For those models (TIM3 + Γ, TVM + Γ) that are not implementable in MrBayes, we followed the manual's recommendation (Ronquist et al., 2020) to substitute the preferred model with the next more complex model available in the program (GTR + Γ). Two parallel runs of Markov Chain Monte Carlo (MCMC) analyses were executed in MrBayes for 20 million generations with four chains, each starting with a random tree and sampling one tree every 1000th generation. Convergence of the two runs was accepted when the average standard deviation of split frequencies (ASDSF) dropped below 0.01 and Potential Scale Reduction Factor (PSRF) of Convergence Diagnostic was 1.00. Tracer 1.7.2  was used to inspect the convergence of model parameters and check whether the values of effective sample size (ESS) were ≥200. The first 25% of the resulting trees were discarded as burn-in, and the posterior probabilities (PP) were then calculated using the remaining trees to construct a 50% majority-rule consensus tree. For the ML analyses, a partitioned model (Àq) was selected and 1000 bootstrap iterations (À#|À N) were conducted, followed by a search for the best-scoring ML tree in a single program run (Àf a).
Phylogenetic trees with PP and BS values were visualized using TreeGraph 2 (Stover and M€ uller, 2010). Clades with PP ≥ 0.99 and/ or bootstrap supports (BS) ≥ 85% were considered as strongly supported, 0.95 ≤ PP < 0.99 and/or 75% ≤ BS < 85% as moderately supported, 0.90 ≤ PP < 0.95 and/or 50% ≤ BS < 75% as weakly supported, and PP ≤ 0.90 and/or BS ≤ 50% as poorly supported. The eight loci were initially analyzed separately and topological incongruence among the eight single gene trees was visually inspected. Three combined DNA data sets (PPR, nrDNA, cpDNA) were then reconstructed for phylogenetic analyses, and the incongruence among these data sets was further evaluated based on visual comparison of the resulting trees. The thresholds of support values for hard incongruence were defined as: PP ≥ 0.99 and/or BS ≥ 85%.

Morphological studies
To test the congruence between morphological and molecular data within Dracocephalum, morphological similarities and differences were analyzed based on our previous field investigations and specimen examination. Specimens of Dracocephalum and related genera from 26 herbaria (A, BM, BNU, C, CDBI, E, G, HIB, HIMC, HNWP, IBSC, IFP, K, KUN, LE, MW, NAS, P, PE, PEY, SM, SZ, W, WUK, XJBI, YAK; abbreviations follow Thiers, 2021) were examined. The systematic significance of nutlet morphology was also investigated using light microscopy (LM) and scanning electron microscopy (SEM). Mature nutlets were collected from natural populations and herbarium collections, and also from the Germplasm Bank of Wild Species in Southwest China, Kunming Institute of Botany. A total of 38 species of Dracocephalum, four species of Hyssopus, and four species of Lallemantia were selected, among which 18 species were sampled for the first time (Table S2). Except for some species where only 1-5 nutlets were available, most species were represented by 20 mature nutlets. Measurements and LM analysis were carried out using the Keyence VHX-6000 digital microscope (Keyence Corporation, Osaka, Japan). For SEM, nutlets were directly mounted onto stubs and sputter-coated with gold using the Cressington 108 auto sputter coater (Cressington Scientific Instruments Ltd, Watford, UK) for 90 s at 20 mA. Micromorphological observations were conducted using a Zeiss EVO LS10 scanning electron microscope (Carl Zeiss NTS, Oberkochen, Germany) at 10 kV. Terminologies used for nutlet description followed those of Budantsev and Lobova (1997) and Moon et al. (2009).

Evolution of chromosome numbers
The utility of chromosome numbers in the infrageneric classification of Dracocephalum and related genera was investigated by reconstructing the ancestral chromosome number using ChromEvol 2.0 (Glick and Mayrose, 2014). The chromosome counts were obtained from the Chromosome Counts Database (http://ccdb.tau.ac.il; Rice et al., 2015), and the haploid counts were used for subsequent analysis (Table S3). When several alternative counts were reported for the same species, each alternative count was followed by its frequency (i.e. probability), and the sum of all possible frequencies was equal to 1.0. We examined the data under eight different available models, each including different combinations of the rate parameters gain, loss, duplication, and demi-duplication, with the rates for each treated either as constant or linear. A Bayesian consensus tree reconstructed from the combined PPR matrix and pruned to include only one accession of each species with reported chromosome numbers, was used as an input tree. The two species of Schizonepeta (Benth.) Briq. were retained here as the outgroup, and the remaining outgroup species were removed. The model with the lowest AIC score (i.e. "CONST RATE") was selected and changes in chromosome number across the tree were visualized using FigTree 1.4.4 (Rambaut, 2018).

Divergence time estimation
Divergence times were estimated using BEAST 1.10.4 . Three separate analyses were conducted using the combined PPR data set, the combined nrDNA data set, and the combined cpDNA data set, respectively, and the three data sets were each partitioned by DNA regions in an attempt to accommodate sequence rate heterogeneity. Multiple accessions from a single species were pruned prior to the analyses. Models of nucleotide evolution for each DNA marker were estimated using jModelTest based on the AIC. The suggested best-fitting models were GTR + I + Γ for AT1G09680 and ITS, and GTR + Γ for AT3G09060, ETS, ycf1-rps15, rpl32-trnL, trnL-trnF, and ycf1. Dates were estimated with a lognormal relaxed molecular clock and a Yule process speciation prior. Due to the lack of fossil evidence for the subtribe Nepetinae, the 95% highest posterior density (HPD) from the most recent Lamiaceae-wide dating analyses (Rose et al., 2022) was used here as a constraint for the crown of the Dracocephalum-Schizonepeta clade using a normal distribution, with a mean of 16.5 and a standard deviation (SD) of 2.8. Two separate runs were conducted, each with 1 billion generations and sampled every 5000 generations. The resulting log files and trees from each run were combined with the program LogCombiner 1.10.4 (included in the BEAST package). The combined log file was checked using Tracer to ensure the ESS were all above 200. After assessing the results in Tracer, we discarded the first 25% of the trees as burn-in. The remaining trees were then interpreted by TreeAnnotator 1.10.4 (included in the BEAST package). Finally, the maximum clade credibility (MCC) chronogram was visualized using FigTree.

Diversification analysis
The diversification rate change over time was explored using several programs based on R applications. Only the combined PPR data set was used for this analysis. To display diversification dynamics, a set of 1000 randomly selected trees from the BEAST analysis and the MCC tree with outgroups removed were used to compute semi-logarithmic lineage-through-time (LTT) plots in the R package APE 5.5 (Paradis et al., 2004). Bayesian analysis of macroevolutionary mixtures (BAMM), implemented in BAMM 2.5.0 (Rabosky, 2014), was used to infer diversification rates and rate shift across the tree. The analysis was conducted with the MCC tree after removing all outgroups. Priors were set using the "setBAMMpriors" function in the R package BAMMtools 2.1.7 . The "ExpectedNumberofShifts" was set as 1 and the species sampling percentage ("globalSamplingFraction") was defined as 0.81. Four Markov chains were run for 1 billion generations and sampled every 100 000 generations. The results with ESS >200 were analyzed and plotted in BAMMtools with the first 25% of the sampled data discarded as burn-in. The net diversification rates were plotted using the "plotRateThroughTime" function.

Ancestral area reconstruction
Ancestral area reconstructions (AARs) were conducted using both the Statistical Dispersal-Vicariance Analysis (S-DIVA) (Yu et al., 2010) and Lagrange under the Dispersal-Extinction-Cladogenesis (DEC) model (Ree and Smith, 2008); both were implemented in RASP 4.2 (Yu et al., 2015(Yu et al., , 2019. Based on the extant distribution data of individual species compiled from taxonomic and floristic literature (Schischkin, 1954;Heywood and Richardson, 1972;Wu and Wang, 1977;Davis, 1982;Rechinger et al., 1982) and herbarium records, six geographical regions were delimited: (i) Europe; (ii) North America; (iii) East Asia, including areas of north and northeast China, Korea, Japan, and parts of Russia; (iv) QTP and Table 2 Properties of data sets used for phylogenetic reconstructions Data set Best-fitting substitution The "/" indicates a model was not available. adjacent area, including the Himalayas and Hengduan Mountains; (v) Central and West Asia, and southern Siberia; (vi) North Africa. The AAR analyses were carried out on the MCC tree and a subset of 1000 random Bayesian trees from the BEAST analysis of the PPR data set. The maximum number of ancestral areas at each node was set as two in the analyses. Being the sister group of Dracocephalum, the two species of Schizonepeta were chosen as the outgroups, and the remaining outgroup taxa were pruned prior to the AARs.

Phylogenetic reconstruction
Topological comparison of the eight single gene trees revealed no strongly supported conflict between the two PPR genes, the two nrDNA loci, and among the four cpDNA markers. Thus, the PPR data sets, the nrDNA data sets, and the cpDNA data sets were each concatenated and analyzed 1989. Properties for different data sets are summarized in Table 2. The final combined PPR data set included 2257 positions (921 bp for AT3G09060, 1336 bp for AT1G09680) and the combined nrDNA data set was 1307 bp (818 bp for ITS, 489 bp for ETS). The combined cpDNA data set consisted of 3657 bp (1009 bp for rpl32-trnL, 886 bp for trnL-trnF, 1110 bp for ycf1, 652 bp for ycf1-rps15) after excluding ambiguously aligned regions.
Apart from poorly supported or collapsed nodes, the BI, MP, and ML analyses yielded similar topologies, thus only the Bayesian 50% majority-rule consensus trees of the three combined data sets are presented, the PP and BS support values being superimposed on the nodes (Fig. 2, Figs S1 and S2). Several strongly supported topological differences can be recognized among the resulting trees. Combining any two of the three data sets (PPR, nrDNA, cpDNA) resulted in lower support values or unresolved polytomies at the conflicting nodes (results not shown), so we did not concatenate the three data sets for analysis. The final alignments and all phylogenetic trees were deposited in the Dryad Digital Repository (https://doi.org/10.5061/ dryad.n5tb2rbxm).
All three phylogenetic reconstructions place the species of Hyssopus and Lallemantia deep within Dracocephalum, indicating that Dracocephalum is paraphyletic with respect to these genera, whereas Dracocephalum s.l. is robustly supported (Fig. 2 and Fig. S2, PP = 1.00/MPBS = 100%/MLBS = 100%; all values are reported in this order below; Fig. S1, 1.00/ 75%/97%). Nine strongly-supported clades (Clades A-I) are recognized in the PPR tree (Fig. 2), but only six of them (except Clades E, F, and I) are supported in the nrDNA (Fig. S1) and cpDNA trees (Fig. S2). Relationships among these clades also exhibit incongruences among the three trees, e.g. Clade B is recovered as sister to Clade C in the PPR tree (Fig. 2,  1.00/59%/82%), but is either grouped with D. spinulosum Popov and Clade G in the nrDNA tree (Fig. S1, 1.00/99%/100%) or sister to one of the lineages of Clade E in the cpDNA tree (Fig. S2, 1.00/ 68%/93%). Specific relationships within some clades show conflicts as well, especially between the nuclear trees and the plastid tree. Since the PPR tree (Fig. 2) is the best resolved one and shows the best congruence with morphology (see below), our further discussions of phylogenetic relationships and biogeographic and diversification results will prioritize the PPR data set.

Nutlet morphology
Nutlets of Dracocephalum s.l. are 1.7-4.6 mm long, 0.8-2.4 mm wide, and range from yellowish brown to dark brown or black ( Fig. S3; Table S2). Most species are characterized by trigonous-oblong nutlets (length/ width ≥2) having truncate or rounded apex and triangular or acute base with a lateral bilobed areole, whereas the nutlets of Clade A are usually trigonousovate (length/width ≤1.8), having truncate or rounded base with a straight areole. The nutlet surface is smooth with four major sculpturing patterns: (folded) lumpy, cellular, reticulate, and striate ( Fig. S4A-H; Table S2). The lumpy type is most common in the genus and occurs in all clades except Clade I. Species of Clade I are characterized by either a cellular or reticulate nutlet surface, but the cellular type also occurs in Clades A, C, and E. The striate type can only be found in species of Hyssopus from Clade F. Moreover, nutlets are typically glabrous in Dracocephalum species, whereas those of Hyssopus are covered with glandular or both glandular and eglandular trichomes at the apex (Fig. S4I).

Chromosome number evolution
The ancestral haploid chromosome number of Dracocephalum s.l. is estimated to be n = 7 (Fig. 3), which is also possibly the ancestral haploid chromosome number of Clades A, B, D, G & H. Clades C & I possess the same ancestral type of n = 6, whereas species  Table S1. of Clades E & F have three different base chromosome numbers each (x = 5, 6, 7) and the ancestral base number for these clades is uncertain (Fig. 3).

Divergence time estimates and diversification analyses
According to the recent molecular dating analyses of Lamiaceae (Rose et al., 2022), Dracocephalum s.l. originated during the middle Miocene, with the stem of the genus diverging from Schizonepeta at ca. 16.32 million years ago (Ma) (95% HPD = 10.98-21.99 Ma).
Using this time interval as a secondary calibration point in our BEAST analysis of the PPR data set, our result indicates that Dracocephalum originated at ca. 14.83 Ma ago (95% HPD = 8.97-20.7 Ma), and started to diversify at ca. 11.12 Ma (95% HPD = 6.3-15.95 Ma) ( Fig. 4; Table 3). During the middle to late Miocene (to ca. 7.24 Ma), major lineages of the genus had emerged. Since the relationships among major clades of Dracocephalum show strong conflicts among the PPR, nrDNA, and cpDNA data sets, the time intervals estimated at some nodes vary among the three chronograms (Fig. 4, Figs S5 and S6). However, all three chronograms show that the genus might have experienced an early rapid diversification during 11-7 Ma in the middle to late Miocene. The results of the divergence times for major lineages of Dracocephalum are summarized in Table 3. The LTT plot result shows that the number of lineages of Dracocephalum increased steadily (Fig. 4A), which is consistent with the result from BAMM analysis in which the net diversification rate decreased steadily without any dramatic shift or radiation (Fig. 4B).

Ancestral area reconstruction
The ancestral distributions inferred from the Lagrange (Fig. 5) and S-DIVA (Fig. S7) analyses based on the PPR data set yielded similar results. Central and West Asia and southern Siberia are inferred as the most likely ancestral area of Dracocephalum s.l., as well as of most of the clades A-I therein. The genus later dispersed to Europe and the Sino-Japanese region between 5 and 4 Ma. Two independent migrations from Central and West Asia to the QTP and adjacent regions occurred between 4 and 1 Ma, one in Clade C and the other within the clade formed by Clade G and Clade H.

Discussion
Possible causes of incongruences among the PPR, nrDNA, and cpDNA data sets Topological comparisons of the PPR, nrDNA, and cpDNA phylogenetic trees reveal strongly supported hard incongruences among the three data sets (Fig. 2, Figs S1 and S2), indicating that the three data sets exhibit different histories within Dracocephalum s.l. Therefore, we did not combine these data sets for analyses. Before proceeding to the detailed interpretation of phylogenetic relationships within the genus, we will first consider some possible explanations for the incongruences in our results.

Ancient rapid radiation
The relatively poor resolution at many deep nodes and uncertain relationships among some of the main clades of Dracocephalum s.l., together with the short lengths of the deeper internal branches compared with those of the external branches in all three trees (Fig. 2, Figs S1 and S2), suggests a pattern of "ancient" rapid radiation (Rokas et al., 2005;Whitfield and Lockhart, 2007). Ancient here refers to a high ratio of time since divergences and the time span in which they occurred (Whitfield and Kjer, 2008). This pattern is further supported by our results from BEAST analyses (Fig. 4), which suggest that major lineages of Dracocephalum s.l. emerged during a short timescale at ca. 11-7 Ma. The DNA markers used in this study failed to resolve more relationships at deeper nodes with convincing support; phylogenomic data using high-throughput sequencing might further help to address this issue. As for the strongly supported conflicts in the relationships among clades between the three data sets (e.g. the placements of Clade B and Clade G), incomplete lineage sorting (ILS) accompanied by rapid radiation might be one of the causes. The impact of ILS can be significant when the time elapsed between divergences is short, because the opportunity for the presence of alternative alleles in different lineages increases (Whitfield and Lockhart, 2007;Oliver, 2013). In these situations, different genes would be expected to suggest different phylogenetic relationships.

Chloroplast capture
While comparing the nuclear trees ( Fig. 2 and Fig. S1) and the cpDNA tree (Fig. S2), several hard incongruences can immediately be identified. One of these conflicts involves the D. moldavica L. + D. foetidum Bunge clade, which is recovered in the strongly supported Clade E in the PPR tree (Fig. 2, 1.00/95%/99%), but is grouped with members of Clade F and Clade I in the cpDNA tree (Fig. S2,  1.00/100%/100%). Another instance of incongruence is found in the placement of the clades comprising D. paulsenii Briq. and D. discolor Bunge, which are shown to be members of Clade I (Fig. 2, 1.00/100%/ 100%; Fig. S1, 1.00/88%92%) in the nuclear trees, but form a clade sister to Hyssopus and D. psammophilum C.Y. Wu & W.T. Wang from Clade F in the cpDNA tree (Fig. S2, 1.00/100%100%). Moreover, species represented by multiple accessions tend to be monophyletic in the nuclear trees, but some of these species are polyphyletic in the cpDNA tree, with individuals of the same species recovered in two or more wellsupported clades with other species. Examples of such problematic species are D. integrifolium Bunge and D. tanguticum Maxim. The high support values associated with these relationships suggest that the cytonuclear discrepancies are most likely results of chloroplast capture, i.e. the introgression of a chloroplast genome from one species to another (Rieseberg and Soltis, 1991;Tsitrone et al., 2003). Cases of chloroplast capture have been reported in an increasing number of studies (e.g. Acosta and Premoli, 2010;Folk et al., 2017;Lee-Yaw et al., 2018;Liu et al., 2020). Organellar (both mitochondria and chloroplast) DNA is preferentially (maternally) transmitted across species boundaries, displacing the resident cytoplasmic genome by an alien one through intrageneric or intergeneric hybridization and recurrent backcrossing (Rieseberg and Soltis, 1991;Funk and Omland, 2003;Chan and Levin, 2005;Stegemann et al., 2012). Therefore, chloroplast capture may distort phylogenetic relationships reconstructed using chloroplast genes, resulting in species-level polyphyly or misleadingly close relationships between genetically or morphologically distinct species.

Constraints of nrDNA for phylogenetic reconstruction
The two nuclear trees ( Fig. 2 and Fig. S1) also exhibit topological conflicts with each other, mainly in the placements of five species, i.e. D. spinulosum and D. komarovii Lipsky of Clade E, D. nodulosum Rupr. and D. integrifolium of Clade F, and D. oblongifolium Regel of Clade G. For example, in the nrDNA tree, D. nodulosum is grouped with species of Clade I (Fig. S1), whereas D. integrifolium is sister to the D. heterophyllum Benth. and D. stamineum of Clade E (Fig. S1).
Despite the advantages of nrDNA (e.g. ITS, ETS, and 5S-NTS) for phylogenetic reconstruction, such as biparental inheritance, universality, and intergenomic variability, it may confound evolutionary relationships since it suffers from issues of multiple rDNA arrays, pseudogenes, concerted evolution, secondary structure, and homoplasy due to alignment difficulty ( Alvarez and Wendel, 2003;Feliner and Rossell o, 2007). To exclude potential sample contamination or sequence errors, multiple accessions of the same species were sequenced in this study, and replicate experiments were also carried out for any suspicious sequences. Moreover, our results reveal that matrices with or without gaps removed produce largely consistent topologies (not shown), indicating that the incongruences were not from alignment uncertainty. Although the lowcopy nuclear genes are less subject to concerted evolution and homoplasy, they are not devoid of issues including paralogous copies, ILS, hybridization and introgression. Considering that the PPR tree is most consistent with the morphological evidence, we infer that some of the above-mentioned constraints of nrDNA sequences have possibly contributed to the incongruences. Further studies based on the highthroughput sequencing data may help to clarify of the evolutionary history of Dracocephalum s.l.
Based on a comprehensive sampling of Dracocephalum, our results shed new light on the relationships within the genus. Dracocephalum is shown to be paraphyletic with respect to the two genera Hyssopus and Lallemantia (Fig. 2, Figs S1 and S2), which corroborates the results of Drew and Sytsma (2012). Nine major lineages are recognized based on the PPR tree (Fig. 2), some of which are largely consistent with the sections in Budantsev (1987). Our morphological studies further show that the shape and sculpturing of nutlets, which had been used earlier for the infrageneric classification of Dracocephalum (Budantsev, 1987;Budantsev and Lobova, 1997), are of limited systematic significance. However, these features can be diagnostic characters for some clades/species, e.g. the trigonous-ovate nutlets with truncate bases and straight areoles in most species of Clade A, and the striate nutlet surface in species of Hyssopus (Figs S3 and S4; Table S2).
Our further discussion of the phylogenetic relationships within Dracocephalum s.l. mainly follows the PPR tree (Fig. 2), with the relationships recovered in the nrDNA tree (Fig. S1) and the cpDNA tree (Fig. S2) being compared and discussed. We provide each major clade recovered in the PPR tree with an informal name, which is in most cases based on the section name of Budantsev (1987), with the type of the section being represented in the clade. The descriptions below are based on both molecular phylogenetic and morphological evidence, beginning with the first diverging clade within Dracocephalum s.l.

Ruyschiana Clade (Clade A)
This clade is strongly supported in all three trees ( Fig. 2 and Fig. S1, 1.00/100%100%; Fig. S2, 1.00/ 99%100%). At the same time, it is only recovered as the earliest diverging lineage in the PPR tree with weak support of the successive sister node (Fig. 2, 0.78/-/74%), and the relationship of the two groups is effectively unresolved in the nrDNA and cpDNA trees (Figs S1 and S2). Five species are included in this clade: D. parviflorum Nutt., D. grandiflorum L., and all three species of subg. Ruyschiana. Species of Clade A are widely distributed from North America to the Sino-Japanese region and Siberia, as well as in Central Asia and Europe. The clade is morphologically heterogeneous in both the vegetative and floral features. The North American endemic species D. parviflorum was initially treated within a monotypic series [ser. Cryptodracon (Benth.) A.L. Budantsev] of sect. Dracocephalum by Budantsev (1987), which was later elevated by him to a rank of a (monotypic) section, i.e. sect. Cryptodracon Benth. (Budantsev, 1993). The species is distinct in its dense and spike-like inflorescence with bracts longer than the calyces and corollas barely protruding from calyces. It has been considered to be closely related to D. moldavica and D. foetidum (Budantsev, 1993), but D. parviflorum differs from the two species and other taxa of sect. Dracocephalum by having calyces with the teeth of the posterior lip as long as the lip (vs. half as long as the lip). Placed in sect. Calodracon by Budantsev (1987Budantsev ( , 1993, D. grandiflorum has basal leaves with long petioles as do other members of the section, whereas species treated by him within subg. Ruyschiana are distinct in their pubescent anthers. Despite the morphological differences among the species of Clade A, they share basic chromosome number of x = 7 (Fig. 3) and trigonous-ovate (length/width ≤1.8) nutlets with truncate bases and straight areoles (except for the areoles of the nutlets of D. grandiflorum which are lateral bilobed) (Fig. S3). However, the cellular surface pattern of nutlets of D. grandiflorum is similar to that of the three species of subg. Ruyschiana (Fig. S4E), differing from the folded-lumpy pattern in nutlets of D. parviflorum (Fig. S4A).

Peregrina Clade (Clade B)
Most species of ser. Peregrina Schischk. of sect. Dracocephalum (Budantsev, 1987) are recovered in this robustly supported clade ( Fig. 2 and Fig. S1, 1.00/ 100%/100%; Fig. S2, 1.00/77%/98%), except for D. heterophyllum and D. spinulosum (both in Clade E). This clade is morphologically most similar to the other taxa of sect. Dracocephalum in having a distinctive posterior calyx lip, with the three teeth almost equal in size and only half as long as the posterior lip. However, most of the remaining species of sect. Dracocephalum are grouped in Clade E. The cpDNA tree recovers Clade B as sister to Clade E (Fig. S2, 1.00/ 68%/93%), but Clade B is shown to be more closely related to Clade C in the PPR tree (Fig. 2, 1.00/59%/ 82%), and grouped with D. spinulosum (one of the splits of Clade E) and Clade G in the nrDNA tree (Fig. S1, 1.00/99%/100%).
Dracocephalum bipinnatum Rupr., D. diversifolium Rupr., and D. peregrinum L. form a well-supported clade and are very similar to each other in floral morphology. They also share basic chromosome number of x = 6 (Fig. 3). The major difference between the three species lies in the leaf margin shape, which is pinnatisect in D. bipinnatum, entire to pinnatipartite in D. diversifolium, and spinescent-denticulate in D. peregrinum. Relationships among the three species are not resolved in the nrDNA and cpDNA trees (Figs S1 and S2), whereas the PPR result (Fig. 2) shows that individuals of D. bipinnatum and D. diversifolium are grouped with each other and that the two species are possibly conspecific. They both have nutlets with a lumpy surface, whereas the nutlet surface is reticulate in D. peregrinum. The first diverging species in Clade B, D. nuratavicum Adylov, differs from the above three species in its light yellow (vs. bluish-purple) corolla and crenate-dentate lamina margin, which is more similar to that of the other species of sect. Dracocephalum, especially D. heterophyllum. All four species of Clade B are mainly distributed in Central Asia. Potential morphological synapomorphies for this clade need to be further investigated.
Two major types of sculpturing of nutlet surface can be recognized for the species of Clade C: lumpy and cellular ( Fig. S4; Budantsev and Lobova, 1997). Only two species of Clade C have had chromosome numbers reported, both having a base number of x = 7 (Fig. 3).

Lallemantia Clade (Clade D)
Recent morphological, taxonomic, and molecular phylogenetic studies (Budantsev, 1993;Budantsev and Lobova, 1997;Drew and Sytsma, 2012) have consistently supported the inclusion of Lallemantia in Dracocephalum. Our present study further corroborates that Lallemantia is best treated within Dracocephalum. Although species of Lallemantia can be easily distinguished from the other taxa of Dracocephalum by having flattened pedicels and corollas barely exserted from the calyces [except for L. canescens (L.) Fisch. & C. A. Mey.], the presence of D. scrobiculatum Regel in Clade D ( Fig. 2 and Fig. S2, 1.00/100%/ 100%) renders this clade as morphologically heterogeneous. Species of Lallemantia and D. scrobiculatum also differ in the basic chromosome number, which is x = 7 in Lallemantia but x = 6 in D. scrobiculatum (Fig. 3). At the same time, these species share lumpy nutlet surfaces ( Fig. S4; Budantsev and Lobova, 1997). Lallemantia is mainly distributed in West Asia, with L. royleana (Benth.) Benth. further extending to Central Asia, whereas D. scrobiculatum is native to Central Asia.

Dracocephalum Clade (Clade E)
Clade E is strongly supported in the PPR tree (Fig. 2, 1.00/95%/99%), but is split into four and two separate lineages in the nrDNA tree (Fig. S1) and cpDNA tree (Fig. S2), respectively. This clade largely corresponds to sect. Dracocephalum in Budantsev (1987), but does not include D. parviflorum recovered in Clade A and all the species of Clade B. Dracocephalum stamineum has been treated within the genus Fedtschenkiella by some authors (Kudrjaschev, 1941;Huang, 1977) or in a monotypic subgenus of Dracocephalum (Schischkin, 1954;Budantsev, 1987Budantsev, , 1993, mainly due to its long-exserted stamens. However, D. stamineum possesses two-lipped calyces with the teeth of the posterior lip half as long as the anterior lip, as the other species of Clade E do. They also share similar lumpy nutlet surfaces ( Fig. S4; Budantsev & Lobova, 1997).
Two well-supported subclades can be recognized within Clade E, one of which (Fig. 2, 1.00/95%/98%; Fig. S2, 1.00/100%/100%) is mainly distributed in Central and West Asia to North Africa, and the other (Fig. 2, 1.00/63%/77%) extends from Central to East Asia and Siberia. Except for D. renati, which is endemic to North Africa, the species of the first subclade are mostly restricted to the Irano-Turanian floristic region (Takhtajan, 1986). This subclade (except for D. komarovii and D. spinulosum) corresponds to subsect. Aristodracon A.L. Budantsev in sect. Dracocephalum (Budantsev, 1987), and species within the subclade are similar to each other in both vegetative and floral morphology.
Hyssopus was recently demonstrated to be closely related to Dracocephalum and transferred to subtribe Nepetinae by Drew and Sytsma (2012). It had been previously placed in various subtribes of Mentheae, including subtribe Menthoideae (i.e., Menthinae) sensu Bentham (1876), subtribe Hyssopinae sensu Briquet (1895Briquet ( -1897 [also by Wunderlich (1967) and Li (1977)], and subtribe Menthinae sensu Harley et al. (2004), mainly based on its long-exserted and divergent stamens with the anterior pair longer than the posterior one (species of Nepetinae have parallel stamens with the posterior pair longer). However, Hyssopus shares some diagnostic characters with Dracocephalum, including thickened folds at each sinus between the calyx teeth, anther thecae diverging at 180°relative to one another, and aristate bracts.
Species of sect. Idiodracon are erect herbs or shrubs without basal leaves, having entire laminae that are crenate or serrate, a lanceolate or broadly obovoid medium lobe of the posterior calyx lips, and erect or bent posterior corolla lips (Budantsev, 1987). They also share lumpy nutlet surfaces. Although species of Hyssopus are morphologically quite distinct from the remaining species of Clade F by having the abovementioned long-exserted and divergent stamens with the anterior pair being longer, as well as having nutlets with a striate surface and glandular trichomes (Fig. S4H,I), they exhibit similarities in the bract and calyx morphology. Most species of Hyssopus have multiple stems, aristate bract apices, lanceolate calyx teeth with aristate apices. These characters are similar to those of species from ser. Fruticulosa Schischk. of sect. Idiodracon, especially D. fruticulosum Stephan ex Willd.
Species of Clade F are widely distributed in Central and West Asia, Siberia, and Europe, with H. officinalis L. further extending to northwest Africa.

Confertodracon Clade (Clade G)
This clade includes only one species, D. oblongifolium. Although the three accessions of this species group together in all the three reconstructions (Fig. 2, 0.88/62%/90%; Fig. S1, 1.00/96%/100%; Fig. S2, 1.00/99%/100%), the placement of the clade differs among the three data sets. It is sister to Clade H in the PPR tree (Fig. 2, 1.00/100%/100%), but groups with D. spinulosum and Clade B in the nrDNA tree (Fig. S1, 1.00/99%/100%), and is sister to the large clade formed by Clades H, B, E, I & F in the cpDNA tree (Fig. S2, 1.00/100%/100%). Together with D. scrobiculatum and D. schischkinii Strizhova, D. oblongifolium was placed in sect. Confertodracon by Budantsev (1987). This section is characterized by having dense ascending stems, crenate-dentate to pinnately lobed leaves, lanceolate or triangular-lanceolate calyx teeth, and bent posterior corolla lips. However, D. scrobiculatum has been shown to be a member of Clade D, and it differs from D. oblongifolium by having sessile leaves (vs. with petioles 2-10 mm long), a basic chromosome number of x = 6 (vs. x = 7) (Fig. 3), and lumpy (vs. reticulate) sculpturing of nutlet surfaces (Budantsev and Lobova, 1997). Both species are restricted to Central Asia with a sympatric distribution, and thus their morphological similarity might result from convergent evolution and adaptation to similar habitats. Dracocephalum schischkinii was not sampled in our molecular phylogenetic analyses. Based on our specimen examination and comparison of the three species, we conclude that D. schischkinii is morphologically more similar to D. oblongifolium in the petiolate laminae, and therefore it is possibly also a member of Clade G. Further evidence is needed to confirm the phylogenetic placement of D. schischkinii.

Sinodracon Clade (Clade H)
Section Sinodracon was established by Wu and Wang (1977) and later accepted by Budantsev (1987). Except for D. isabellae Forrest that is not included in our study, all species of the section are recovered in Clade H, which is strongly supported (Fig. 2, 1.00/ 91%/97%; Fig. S1, 1.00/99%100%; Fig. S2, 1.00/ 100%100%). Clade H includes no other species except for those of sect. Sinodracon. As a well-defined lineage, sect. Sinodracon is distributed in the QTP region and can be distinguished from the remaining species of Dracocephalum s.l. by the following combination of characters: plants erect without basal leaves, laminae pinnatipartite to pinnatisect, calyces slightly two-lipped with subequal lanceolate lobes, and posterior lip of corolla erect or bent (Wu and Wang, 1977;Budantsev, 1987Budantsev, , 1993Chen et al., 2021). Budantsev (1985) argued that the pinnatifid laminae with linearlanceolate lobes of the section might indicate its close relationship with subg. Ruyschiana placed in Clade A in our phylogenetic reconstruction; however, Clade H is shown to be sister to Clade G in the PPR tree (Fig. 2, 1.00/100%/100%) and (relatively) distantly related to Clade A. Species of Clade G and Clade H share similar floral morphology and the same basic chromosome number of x = 7 (Fig. 3), but the laminae are crenate-dentate to pinnately lobed (vs. pinnatipartite to pinnatisect) and the stems are ascending (vs. erect) in species from Clade G.

Keimodracon Clade (Clade I)
This clade is strongly supported in the PPR tree (Fig. 2, 1.00/100%/100%), but is split into several groups in the nrDNA and cpDNA trees (Figs S1 and S2). Species of Clade I are mainly distributed in West and Central Asia and Siberia, with a few extending to Russian Far East. They have either cellular or reticulate nutlet surfaces (Fig. S4) and uniform basic chromosome number of x = 6 (Fig. 3). The two poorly supported subclades recovered within Clade I, D. paulsenii~bungeanum 2 (Fig. 2, 0.62/-/-) and D. palmatum~multicolor (Fig. 2, 0.73/-/-), correspond to sect. Keimodracon and sect. Palmata in Budantsev (1987), respectively. The two sections differ in the shape of the posterior corolla lip, which is straight or bent in sect. Keimodracon, but arched in sect. Palmata. Moreover, stems of the species of sect. Keimodracon are usually short and creeping, but are mostly erect in sect. Palmata.
Diversification and biogeographic history of Dracocephalum s.l.
Species of Dracocephalum s.l. are widely distributed in temperate regions of the northern hemisphere. Most species are drought-tolerant and occur in the steppedesert biomes in Central and West Asia and southern Siberia, and biogeographic reconstructions indicate this region to be an ancestral area of the genus and most of its major lineages ( Fig. 5 and Fig. S7). Our BEAST analyses suggest a middle Miocene origin for Dracocephalum s.l. at ca. 14.83 Ma (Fig. 4). As has been discussed above, the genus might have experienced an early rapid radiation during the middle to late Miocene, as major lineages of the genus arose in a short timescale between 11 and 7 Ma, which is consistent with the origin and diversification of other taxa distributed in the Asian interior arid-zone flora such as . Pollen records also demonstrated that herbs first become widespread in Central Asia at ca. 15 Ma, with a further increase at ca. 7 Ma (Barbolini et al., 2020).
The aridification of the Asian interior from the middle Miocene onwards (mainly from ca. 17 to 5 Ma) has considerably influenced the evolution of many temperate lineages in Asia (Shi et al., 2013). The aridification has been shown to be the result of a series of climatic and geological changes, such as global climate cooling and uplift of the QTP (Miao et al., 2012, and references therein;Hurka et al., 2019). Global cooling provides a dominant driving factor for the drying of central Eurasia, because it reduces the moisture load of the westerlies and decreases precipitation in inland Asia (Lu et al., 2010;Miao et al., 2012). The uplift of the QTP and mountain ranges in the vicinity (including the Tianshan and Pamir) at different times and magnitudes further influenced the moisture pattern, causing decreased temperatures and increased continentality throughout Asia's interior (Molnar et al., 2010;Li et al., 2011;Liu & Dong, 2013). The expansion of arid lands due to orogenetic and climatic changes increased the arid niche space, which might have facilitated the rapid diversification of Dracocephalum s.l. during the middle to late Miocene. As summarized in Wu et al. (2018), the rise of dryland floras in different continents was near-synchronous and began in the middle to late Miocene, indicating this period to be pivotal in the assembly and evolution of modern dryland floras.
Two independent dispersal events of Dracocephalum s.l. from Central and West Asia into the QTP and adjacent regions were detected (Fig. 5). One dispersal event occurred within Clade C between ca. 2.57 and 1.46 Ma, while the second occurred within the clade formed by Clade G and Clade H sometime between ca. 3.72 and 2.79 Ma (Figs 4 and 5). The onset of diversification among QTP taxa in Clade C and Clade H is estimated to have occurred in the early Pleistocene at ca. 1.46 Ma (Fig. 4) and in the late Pliocene at ca. 2.79 Ma, respectively. Geological evidence suggests that there were multiple stages of uplift during the formation of the QTP (Harrison et al., 1992;Tapponnier et al., 2001;An et al., 2001An et al., , 2006Mulch and Chamberlain, 2006). The estimated diversification of Dracocephalum within the QTP coincides with the last major period of extensive QTP uplift at ca. 3.5-1.6 Ma (Harrison et al., 1992;Shi et al., 1998). Other genera have also been reported to diversify during this period in response to the uplift of the QTP, such as Solms-laubachia Muschl. (Brassicaceae; Yue et al., 2009) and Phyllolobium Fisch. (Fabaceae;Zhang et al., 2012). It is notable that the phylogenetic relationships within these QTP clades of Dracocephalum are largely unresolved (Fig. 2); the short branch lengths and polytomies indicate that these lineages might have experienced rapid radiations in the QTP and adjacent areas. The recent tectonic events involving the QTP generated myriad novel environmental niches, ostensibly enabling Dracocephalum to migrate out of the arid Asian interior and adapt to the cold alpine terrain of the QTP and facilitating a rapid diversification of the genus.

Taxonomic treatment
Dracocephalum was subdivided into three subgenera, seven sections, and 21 series by Budantsev (1987). Later, he treated Lallemantia as a subgenus of Dracocephalum and in total recognized four subgenera and eight sections within the genus (Budantsev, 1993). Our results reveal nine well-supported clades within Dracocephalum s.l., and most of these clades are largely corroborated with the sections of Budantsev (1987Budantsev ( , 1993. Considering that all four subgenera of Budantsev (1993) are not monophyletic, and Hyssopus is shown to be embedded within Dracocephalum, a revised classification is needed for the genus. Thus, in this study, we suggest subdivision of Dracocephalum s.l. into nine sections, with each section corresponding to a clade recovered in the PPR tree (Fig. 2).
Both Dracocephalum and Hyssopus were established by Linneaus (1753), therefore, the two genera have equal priority. Since Dracocephalum comprises ca. 70 species, far more than Hyssopus, here we treat Hyssopus as a synonym of Dracocephalum and provide new combinations for the species formerly placed in Hyssopus. Dracocephalum

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Fig. S1. Bayesian 50% majority-rule consensus tree of Dracocephalum s.l. based on combined nrDNA (ITS & ETS) data set. PP ≥ 0.50 and MPBS and MLBS ≥ 50% are displayed above the branches following the order PP/MPBS/MLBS ("-" indicates a support value < 50% BS). Scale bar denotes the expected number of substitutions per site in Bayesian analysis. Multiple accessions of the same species are numbered according to Table S1. Fig. S2. Bayesian 50% majority-rule consensus tree of Dracocephalum s.l. based on combined cpDNA (rpl32-trnL, trnL-trnF, ycf1, and ycf1-rps15) data set. PP ≥ 0.50 and MPBS and MLBS ≥ 50% are displayed above the branches following the order PP/MPBS/ MLBS ("-" indicates a support value < 50% BS). Scale bar denotes the expected number of substitutions per site in Bayesian analysis. Multiple accessions of the same species are numbered according to Table S1. Fig. S3. Nutlet morphology of taxa of Dracocephalum, Hyssopus, and Lallemantia.  Fig. S5. Chronogram of Dracocephalum s.l. as inferred from BEAST analysis of the combined nrDNA data set. Blue bars represent 95% HPD of node age. Fig. S6. Chronogram of Dracocephalum s.l. as inferred from BEAST analysis of the combined cpDNA data set. Blue bars represent 95% HPD of node age. Fig. S7. Ancestral area reconstruction of Dracocephalum s.l. based on the chronogram of the combined PPR data set and extant distributions using S-DIVA with maximum areas constrained as two in RASP. Colors of circle in the tips and letters before species names represent the distribution of each species. The relative proportions of colors in the circle of each node represent the probabilities of the inferred alternative ancestral areas of nodes and the most likely ancestral area of each node is also indicated by letter(s) within the circle. White color in the circle indicates the sum of probabilities of all alternative ancestral areas with probability < 15%. The time scale on the bottom axis of the tree is in millions of years. Table S1. Sequence information for all samples used in this study. A "-" indicates missing data, and an "*" indicates a sequence downloaded from GenBank. Herbarium acronyms are listed in bracts after vouchers. Table S2. Nutlet morphology of Dracocephalum s.l. species sampled in present study. An "*" indicates the species for which the nutlet morphology is studied for the first time. Herbarium acronyms are listed in bracts after vouchers. Table S3. Chromosome data of Dracocephalum s.l. and outgroup species used for the analyses of chromosome number evolution in present study.