Characterization and phylogenetic analysis of the complete chloroplast genome sequences of six species of Wikstroemia (Thymelaeaceae) revealed paraphyletic relationship to the monotypic genus Stellera.

Wikstroemia (Thymelaeaceae) is a diverse genus, spanning both Asia and Australia and also recorded from the Hawaiian Islands. However, due to its medicinal properties and resource utilization in pulp production, genetic studies on these important species have been overshadowed. In this study, the chloroplast genome sequences of six species of Wikstroemia were sequenced and analyzed. The chloroplast genomes of the six species ranged between 172,610 bp (W. micrantha) and 173,697 bp (W. alternifolia), and exhibited a typical genome structure consisting of a pair of inverted repeat (IR) regions separated by a large single-copy (LSC) region and a small single-copy (SSC) region. The six chloroplast genomes were similar with a predicted 139 genes that consisted of 92 or 93 protein-coding genes, 38 tRNA genes, and 8 rRNA genes. The overall GC contents were identical (36.7%). Genome comparative analyses were conducted with the inclusion of two additional published species of Wikstroemia in which the sequence divergence and expansion of IRs in the chloroplast genomes were determined. The phylogenetic analysis inferred that Wikstroemia in its current circumscription is paraphyletic to Stellera chamaejasme, while the ITS-based tree analyses could not properly resolve the phylogenetic relationship between Stellera and Wikstroemia. This nding kindled interest in the proposal to synonymize Stellera with Wikstroemia, which was previously brought up, but rejected due to taxonomic conicts. Nevertheless, this study provides valuable genomic information to aid in the taxonomic implications and phylogenomic reconstruction of Thymelaeaceae. visualized using IRscope. To detect the mutational hotspots and divergence regions in the chloroplast genomes of the eight species, sequence alignment of the chloroplast genome sequences was carried out using Geneious Prime


Introduction
Wikstroemia Engl., consisting of about 70 species, is a diverse genus in the family Thymelaeaceae. Members of Wikstroemia are widely distributed in the Asian and Oceanian regions, and also scattered around the Hawaiian Islands [1]. The species are mostly brous trees, shrubs or subshrubs with a woody rhizome. Several species are cultivated as raw material for pulp production [2,3], while a handful of them are reported to have medicinal properties [4,5].
However, studies on Wikstroemia have been con ned to its utilization in pulp production and pharmacological applications; reports on genetic studies of Wikstroemia are scarce.
The only report on the genetic diversity to date is of Wikstroemia ganpi in Korea using inter simple sequence repeat (ISSR) markers [6] and two, published, complete chloroplast genome sequences of Wikstroemia chamaedaphne and Wikstroemia indica [7,8]. Due to the lack of molecular evidence, taxonomic studies on Wikstroemia have relied solely on morphological characteristics [9]. Ironically, the continuous nature of morphological variation in members of Wikstroemia has led to much taxonomic confusion in attempts to distinguish species and has resulted in ambiguities in taxonomic classi cations between Wikstroemia and its sister genera [9,10]. Furthermore, the di culty in detecting natural hybridization among the species of Wikstroemia due to the possibilities in low reproductive isolation and high genetic similarities suggested that Wikstroemia is a large species complex [11].
Chloroplasts are intracellular organelles in plants that contain the entire machinery necessary for the process of photosynthesis [12]. They provides energy through photosynthesis and play an important role in carbon uptake [13]. The chloroplast genome is a circular double-stranded DNA molecule. In plants, the cp genome is speci cally maternally solely inherited and not disturbed by gene recombination [14]. A typical plant cp genome ranges in size from 120 kb to 217 kb [15]. The complete chloroplast genome has a typical quadripartite structure, including a large single-copy (LSC) region, a small single-copy (SSC) region, and two separate inverted regions (IRs) [16]. Due to the ease of availability and advances in next-generation sequencing, the cost of sequencing an organelle genome has become affordable [17]. The trend in plant chloroplast genome sequencing has become increasing popular, thus contributing to studies of taxonomic classi cation and molecular evolution. Owing to its slower evolution rate and ease to sequencing and assembly due to its considerably smaller size, the chloroplast genome has been receiving much attention among biologist and taxonomist in exploring possibilities of species identi cation, species delimitation, phylogenomics, and molecular ecology of plants with challenging taxonomic problems or unique lifestyles [18,19].
The taxonomical placement of Wikstroemia has been controversial. It has experienced a complicated classi cation process when reviewing members of the Thymelaeaceae. However, Stellera chamaejasme of the monotypic genus Stellera was reported to be sister to Wikstroemia based on combined chloroplast DNA sequences (trnT-trnL, trnL-trnF, trnL intron, and rpl16 intron) [20], while Wikstroemia, along with the other 14 sister genera based on palynology ndings, is taxonomically grouped under the Daphne group of the tribe Daphneae [21,22]. Although the phylogenetic work on members of Thymelaeaceae is actively ongoing [23], phylogenetic relationships in Wikstroemia are likely to be understudied. Constituent genera in Thymelaeaceae have experienced similar molecular challenges in which poor phylogenetic resolution is likely to be caused by low genetic variation in selected molecular markers [23]. Such con icts can be overcome by utilizing genome-scale data sets [24]. At the same time, high divergence regions could be identi ed through genome comparisons, which could aid in future phylogenetic studies of such a diverse genus as Wikstroemia.
In this study, we sequenced the complete chloroplast genomes of six species of Wikstroemia, W. alternifolia, W. canescens, W. capitata, W. dolicantha, W. micrantha, and W. scytophylla, to analyze and compare genomes using bioinformatic tools. Our aims were to: (1) characterize the chloroplast genomes of the six species of Wikstroemia; (2) examine variations in sequence repeats and codon usage in the six chloroplast genome sequences; (3) identify highly divergent regions in the chloroplast genome sequences; (4) improve the understanding of the intrageneric/intergeneric phylogeny of Wikstroemia within the Thymelaeaceae based on chloroplast genome sequences and the nuclear ribosomal DNA internal transcribed spacer (ITS) region.

Chloroplast genome features of six species of Wikstroemia
The total length of the chloroplast genomes of the six species of Wikstroemia analyzed in this study ranged from 172,610 bp (W. micrantha) to 173,697 bp (W. alternifolia). All six chloroplast genomes exhibited the typical quadripartite structure (Figure 1) consisting of a pair of IRs regions (41,850-42,073 bp) separated by an LSC region (86,111-86,7017 bp) and an SSC region (2,799-2,871 bp). All six chloroplast genomes had the same 36.7% GC content. However, the GC content in the chloroplast genome of each species of Wikstroemia was unevenly distributed. The IR region accounted for the highest GC content (38.8--38.9%), followed by the LSC region (34.7-34.9%), while the SSC region was recorded as having the lowest GC content (26.9-29.5%).
All six species of Wikstroemia contained the same number of long repeats (Figure 3a). In general, all of them contained 24 forward repeats and 25 palindromic repeats, except for W. canescens and W. capitata. Long forward repeats that ranged between 30 and 40 bp were found most abundant in W. dolicantha and W. micrantha; while W. alternifolia, W. canescens, W. capitata, and W. scytophylla were recorded with higher number of long forward repeats with the lengths of 41 to 60 bp ( Figure 3b). Long palindromic repeats were equally abundant in W. alternifolia and W. canescens, ranging from 40 to 60 bp and above 60 bp (Figure 3c), while long palindromic repeats were abundant in the range of 30 to 60 bp in W. capitata, W. dolicantha, W. micrantha, and W. scytophylla. Long reverse repeats were only detected in W. canescens and W. capitata, and mostly occurred within the range of 30 to 40 bp ( Figure 3d).

Sequence divergence analysis
The chloroplast genome sequence alignment of eight species of Wikstroemia, using the W. chamaedaphne chloroplast genome as the reference, indicated high sequence conservatism across the chloroplast genomes of eight species of Wikstroemia, but not in the chloroplast genome of W. indica ( Figure 4). As a whole, the size and gene order of the chloroplast genomes in Wikstroemia are well-conserved, but a distinct large gap was observed beginning within the ycf1 gene sequence of the IRa to the 5' region of the trnL-UAG in the IRb of W. indica. Both of the single-copy regions were recorded as having greater sequence divergence than the IR region ( Figure 5). With a Pi-value cut off point of 0.025, eight highly variable gene regions were identi ed; the ndhD-ndhF, ndhF-rpl32, ndhJ, petL-petG, psbI-trnS-GCU, trnG-UCC, trnK-UUU-rps16 and the trnL-UAA-trnF-GAA intergenic spacer regions. Six of these highly variable regions were located in the LSC, while two of them were in the SSC region.

Contraction and expansion in IR region
The genes adjacent to the IR borders were consistent across members of Wikstroemia, except in W. indica, which varied for its adjacent genes at the IRb/SSC (JSB) and IRa/SSC (JSA) border ( Figure 6). Instead of the presence of rpl32 and ndhF genes in the SSC region, adjacent to JSB and JSA respectively, the ycf1 gene was located across both the JSA and JSB in the chloroplast genome of W. indica. The trnL-UAG gene was also placed adjacent to the JSA, in the SSC region of the W. indica chloroplast genome. On the other hand, six species (W. alternifolia, W. chamaedaphne, W. dolicantha, W. indica, W. micrantha, and W. scytophylla) had their rps19 gene crossing the IRa/LSC (JLA) border.

Phylogenetic analysis
The ML and BI trees based on the complete chloroplast genome sequences revealed that all the branch nodes for eight species of Wikstroemia included in the phylogenetic tree were supported with high bootstrap values and Bayesian posterior probabilities (ML: ≥90%; BI: ≥ 95%) (Figure 7). In addition, it was suggested a paraphyletic relationship was present in the genus Wikstroemia. Two species, W. alternifolia and W. canescens, were clustered with Stellera chamaejasme; while six species of Wikstroemia (W. capitata, W. chamaedaphne, W. dolicantha, W. indica, W. micrantha, and W. scytophylla) formed a monophyletic group.
For the ITS sequences in the ML tree revealed a paraphyletic relationship between Wikstroemia and S. chamaejasme; while most of the branch nodes within the Wikstroemia clade were not highly supported (Figure 8a). Strong bootstrap supports were recorded for the sistership between W. alternifolia and W. canescens, and W. micrantha and W. stenophylla. Weakly supported sisterships were present between W. dolicantha and W. scytophylla, and W. capitata and W. ligustrina. Contrarily, the BI analysis displayed a monophyletic relationship within the Wikstroemia clade (Figure 8b). Similar to the ML tree, sisterships were strongly supported between W. alternifolia and W. canescens; and W. micrantha and W. stenophylla, but not between W. dolicantha and W. scytophylla, and W. capitata and W. ligustrina in the BI tree.

Discussion
The plastid genomes of the species of Wikstroemia in this study were rather conserved, which is similar to other angiosperms [25]. The cp genome lengths of the six species of Wikstroemia did not vary much and had cp genome sizes similar to typical angiosperms [26]. The same number and contents of the genes were predicted in this study, suggesting that evolution of the gene sequences was consistent across the six species. Similar to most angiosperms, sequence repeats for A/T were detected to be more abundant than those of G/C in the Wikstroemia cp genomes and such probability was affected by the stability between the nucleotides A/T and G/C in the genome [27,28].
The expansion and contraction of the IR region are major evolutionary events that in uence the length of the cp genomes [29]. Our study indicated that the contractions and expansions of the IR regions exhibited relatively stable patterns within Wikstroemia, with slight variation in which gene recombination between the repetitive sequence or poly-A structure and tRNA could be one of the reasons for the change in length in the IR region [30]. However, W. indica indicated dissimilarity in its IR borders, which differed from most angiosperms [31]. We suspect that the chloroplast genome IR contraction and expansion in W. indica is severe and may be due to extensive gene transfer and larger IR expansion due to the results of the double strand break repair mechanism [32][33][34]. Although W. indica had a smaller cp genome size (151,731 bp) when compared to other species of Wikstroemia sequenced in this study, a higher GC content was detected in its cp genome [8](37.4%). Upon comparison, we found that the chloroplast genome of W. indica had a shorter IR region and larger SSC region when compared to other species of Wikstroemia. Changes in the placement of IR the borders in the chloroplast genome of W. indica was likely due to the contraction of its IR region, causing a loss in number and content of the genes. Among the genes that were not found in W. indica, but present in other species of Wikstroemia were a pair of ndhA, ndhG, and ndhI that were supposed to be present in the IR region; genes such as ccsA, ndhD, ndhE, ndhH, psaC, rps15, and trnL-UAG that were commonly duplicated in the IR regions were reduced to only one copy, and were transferred to the SSC region; while the ndhF and rpl32 genes that were common genes in the SSC region were not detected. Therefore, it can be concluded that the contraction of the IR region that caused gene loss has contributed to the difference in cp genome content between W. indica and the other seven species of Wikstroemia.
Molecular evidence based on chloroplast genome sequences revealed a non-monophyletic relationship between the species of Wikstroemia due to W. alternifolia and W. canescens clustering with S. chamaedaphne. Information on the phylogenetic relationships of the species of Wikstroemia is scarce. Although taxonomic work is tedious in a genus with diverse species, continuous efforts among taxonomists working on the members of Thymelaeaceae have provided some insights on the taxonomic status of Wikstroemia. To provide better insight on the phylogenetic relationships at the nuclear level, we used ITS sequences to perform ML and BI analyses. Unlike phylogenomic tree analyses on the complete chloroplast genome sequences, low bootstrap support and Bayesian posterior probabilities were observed at species level in the genus Wikstroemia. However, the molecular placement of the species of Wikstroemia were identical in both the ML and BI trees, while the most distinct difference between both phylogenetic trees was the placement of S. chamaejasme. In the ML tree based on the ITS sequences, S. chamaejasme was clustered within the Wikstroemia clade; while S. chamaejasme was sister to Wikstroemia in the BI tree. The discordance between the chloroplast and nuclear phylogenies in this study may be due to phylogenetic sorting, convergence, unequal rates of evolution, long branch attraction, and introgression [35]. However, low branch node supports in both the ITS-based ML and BI trees suggested that either the inclusion of additional nuclear gene sequences, or the application of the restriction site-associated DNA sequence (RAD-Seq) technique that integrates up to 10% of the nuclear genome [36], could be helpful in resolving the phylogenetic relationship within Wikstroemia. Evidently, in this study, the use of a single nuclear gene sequence, i.e. ITS, which was thought to be useful in delimitation of many plants at species level [37], is insu cient for resolving the phylogenetic relationships between Stellera and Wikstroemia.
Today, members of Wikstroemia comprise species previously placed under Capura L., Daphne L., Diplomorpha Meisn., Daphnimorpha Nakai, Lonicera L., Passerina L., Restella Pobed., and Stellera L. [1,38]. Eventually, the monotypic genus Stellera, which exhibits highly similar morphological characteristics, had troubled some taxonomists when compared to Wikstroemia. At least ve species were placed under Stellera before they were transferred to Wikstroemia; while many were transferred to allied genera, such as Daphne, Diarthron and Thymelaea, in the tribe Daphneae [38]. This is understandable as Stellera has a longer taxonomic history¸ i.e. back to 1747, when compared to other genera in the Daphneae. Further division of the tribe into several genera emphasized the identity of the species. As a result, the sole species of Stellera, S. chamaejasme, as the type species, is only species left in the genus. Based on the literatures, we found that Wikstroemia has an interesting nomenclatural history, in which two genera, Diplomorpha and Daphnimorpha, were synonymized and were excused. For Stellera to combination with Wikstroemia was proposed before by transferring the type species, S. chamaejasme, under the monotypic subgenus Chamaejasme [39,40]. However, the proposal was rejected as Stellera has priority over Wikstroemia [41] and based on the Rules of Nomenclature, the combination can only be accepted if Stellera is proposed as a nomen genus rejiciendum [42](nov. gen. rejic.). Therefore, we do not exclude the possibility that Stellera could be synonym to Wikstroemia, however, based on the phylogenetic trees in this study, Wikstroemia is polyphyletic and Stellera is unquestionably a sister to Wikstroemia.
A comprehensive history of Wikstroemia was well-reviewed [43]. Wikstroemia was named after the Swedish botanist Johan Emanuel Wikström to commemorate his contributions in botany. However, the use of his name, Wikstroemia, was rst proposed in 1821 by Heinrich Schrader for a genus in the family Theaceae [44]. In the same year, Kurt Sprengel described a new species with a new genus name, Wikströmia glandulosa, which was later corrected to a species of Eupatorium [43]. Neither of those names applied to a genus in the Thymelaeaceae. Eventually, members of Thymelaeaceae were once grouped under the genus Capura, which was established since 1771. Yet, in 1833, Endlicher proposed the genus name Wickstroemia for these species to honor Wikström, causing several nomenclature con icts to occur. Thus the genus Capura was rejected and synonymized under Wickstroemia despite being earlier [45]. The spelling Wickstroemia was corrected to Wikstroemia in 1841, and is retained until today. However, changes in the genus name still occurred after that. Heller did not agree with the use of Wikstroemia Endl. to be practical and he treated it as a homonym of Wikstroemia Schrad [43]. As a consequence, a new genus, Diplomorpha, was established to replace Wikstroemia in Thymelaeaceae in 1841 [46]. Unfortunately, the use of Diplomorpha did not last long as members of Diplomorpha were either corrected to Sauropus or synonymized under Wikstroemia, including W. canescens [1]. Although many species of Diplomorpha were treated as unresolved due to incomplete description of their origin [38,47], based on their epithets, we speculate that these unresolved names could be members of Wikstroemia. It is noteworthy that the genus Daphnimorpha was similarly treated and is synonymized under Wikstroemia due to the high similarity in morphological characteristics when compared to the latter [22], while W. australis is still recognized as the type species for Wikstroemia. Subgeneric classi cation in Wikstroemia was rst proposed to include only two sections, Euwikstroemia and Diplomorpha [48]. The two sections were then included as subgenus Wikstroemia and subgenus Diplomorpha, along with the proposal of a new subgenus, Chamaejasme. However, subgenus Chamaejasme was then rejected due to taxonomic con ict [39][40][41]. Alongside the proposal of a novel genus, Daphnimorpha, which consisted of two species, Daphnimorpha capitellata (synonym Diplomorpha capitellata), and D. kudoi [49], the suggestion to retain subgenus Diplomorpha as an independent genus instead of a subgenus was brought up [50]. However, the suggestion was not heeded; Daphnimorpha was synonymized under Wikstroemia [22] and the subgeneric classi cation of Wikstroemia, consisting only of the subgenera Wikstroemia and Diplomorpha, is generally accepted currently [42,51].
One of the key morphological characteristics proposed to differentiate Wikstroemia from allied genera is the presence of petaloid scales in the ower [39]. However, the presence of disc in the owers was not emphasized in Wikstroemia [10]. Owing to that, some species of Daphne may be misidenti ed as Wikstroemia due to incomplete morphological descriptions in both genera, causing an overlap in features of classi cation. A study conducted on the variations in morphological characteristics for members of Thymelaeaceae revealed that species of Wikstroemia and Diplomorpha have a scaly, subulate or clapper-shaped disc, while species of Daphnimorpha have a fan-shaped or half-cylindrical disc [49]. The three genera share only the same morphological feature of the petaloid scales instead of a disc, which occurs in Daphne, Edgeworthia, and Eriosolena. It was emphasized that discs and petaloid scales should not be treated as homologous structures, but to be considered as mutually exclusive in members of the Thymelaeaceae [52].

Conclusion
To the best of our knowledge, this study presents the rst genome-scale analysis on species of Wikstroemia. The ndings revealed high conservation of genes in the chloroplast genomes. The identi cation of highly variable gene regions in the chloroplast genome sequences of Wikstroemia could potentially be useful in resolving phylogenetic relationship in the genus. A strong sistership between Wikstroemia and the monotypic genus Stellera was present. The ML and BI trees based on the complete chloroplast genome sequences revealed that all the branch nodes for eight species of Wikstroemia included in the phylogenetic tree were supported with high bootstrap values and Bayesian posterior probabilities (ML: ≥90%; BI: ≥ 95%), while the ITS-based tree analyses could not properly resolve the phylogenetic relationship between Stellera and Wikstroemia. Nevertheless, the molecular data obtained in this study will serve as a valuable resource for providing greater insights into the taxonomy and phylogeny of Thymelaeaceae.

Plant Materials and DNA Extraction
Fresh leaf materials of six species of Wikstroemia, W. alternifolia, W. canescens, W. capitata, W. dolicantha, W. micrantha and W. scytophylla, were collected from botanical gardens and natural populations in China (Table 1). Voucher specimens were deposited in the Herbarium of Yunnan Normal University [53]. The total genomic DNA was extracted using the Axygen® AxyPrep Multisource Genomic Miniprep DNA kit (Corning, USA), following the manufacturer's protocol.
Chloroplast Genome Sequencing, Assembly and Annotation A sequence library was constructed and sequencing was performed on an Illumina HiSeq 2500-PE150 platform (Illumina, USA). All raw reads were ltered using NGS QC Toolkit version 2.3.3 with default parameters to obtain clean reads [54]. The plastome was de novo assembled using NOVOPlasty [55] with the rbcL gene sequence of Daphne kiusiana (GenBank accession KY991380) as the seed sequence. Gene annotation was performed in Geneious Prime (Biomatters, New Zealand) using the complete chloroplast genome sequence of W. chamaedaphne (GenBank accession MN563132) as the reference genome. The circular physical map of the chloroplast genome was visualized using OGDRAW [56].

Repeat analyses
Short sequence repeats (SSRs) were identi ed using MISA-web [57], in which parameters for identi cation of perfect mono-, di-, tri-, tetra-, penta-, and hexanucleotide motifs were set for a minimum of 10, 5, 4, 3, 3, and 3 repeats, respectively. Long repeats, including forward, palindrome, reverse and complement repeats, were determined using REPuter [58] with a hamming distance of 3 and a minimal repeat size of 30 bp.

Codon Usage
Protein-coding sequences of each chloroplast genome were extracted and the relative synonymous codon usage (RSCU) was analyzed using MEGA 7 [59].

Comparative genome and divergence analyses
The complete chloroplast genome sequences of two species of Wikstroemia, W. chamaedaphne (GenBank accession MN563132) and W. indica (GenBank accession MN453832), that were available in the NCBI GenBank, were downloaded and included in subsequent analyses. By using the chloroplast genome sequences of W. chamaedaphne as the reference genome, nucleotide variation in the chloroplast genome sequence alignment of the eight species of Wikstroemia were visualized using mVISTA [60] in Shu e-LAGAN mode. To detect the expansion and contraction of the IR region in the chloroplast genomes across the eight species, the IR/SC boundaries of the chloroplast genomes were visualized using IRscope. To detect the mutational hotspots and divergence regions in the chloroplast genomes of the eight species, sequence alignment of the chloroplast genome sequences was carried out using Geneious Prime (Biomatters, New Zealand). Calculations of the nucleotide variability (Pi) among the eight chloroplast genomes were performed using DnaSP v5 [61] with window length of 1,000 bp and a step size of 500 bp.

Polymerase chain reaction and Sanger sequencing
Polymerase chain reaction (PCR) ampli cation was carried out in a 20 µL volume reaction using the ITS universal primer set: 5F: 5'-GGAAGTAAAAGTCGTAA-CAAGG-3' (forward) and 4R: 5'-TCCTCCGCTTATTGATATGC-3'(reverse). The PCR reactions for the nuclear ribosomal DNA ITS region contained 10 µL of 2× Taq PCR Starmix with loading dye (Genstar Biosolutions, China), 0.4 µM of each primer and 20 ng of genomic DNA as a template. PCR ampli cations was conducted on a T100™ Thermal Cycler (Bio-Rad, USA), with an initial denaturation at 93°C for 5 min; 40 cycles of denaturation at 93°C for 30 s, annealing at 60°C for 30 s, extension at 72°C for 30 s; and a nal extension at 72°C for 5 min. PCR products were sent for direct Sanger sequencing at both ends using an ABI 3730 DNA Analyzer (Applied Biosystems, USA).

Phylogenetic Analyses
The phylogenetic analysis was conducted on 17 complete chloroplast genome sequences of the Thymelaeaceae. Two species, Psidium guajava (Myrtaceae; GenBank accession KY635879) and Gossypium gossypioides (Malvaceae; GenBank accession HQ901195) were included as outgroups. Sequence alignment was carried out using MAFFT [62]. Maximum-likelihood (ML) tree was constructed using RAxML 8.2.11 [63], in which the general-time-reversible (GTR) and gamma distributed (+G) (+GTR+G) DNA substitution model was selected and all branch nodes were calculated under 1,000 bootstrap replicates; while the Bayesian inference (BI) was conducted using MrBayes [64], in which the Markov Chain Monte Carlo (MCMC) was conducted with 2,000,000 generations and sampling was collected at every 100 cycles. Both tree analyses were conducted through the online programs available in the CIPRESS Science Gateway web portal The ITS sequences were aligned and manually trimmed for their primer sequences to obtain clean sequences. A total of 26 additional ITS sequences derived from members of the Thymelaeaceae were downloaded from the NCBI GenBank and MUSCLE-aligned against the ITS sequences of the six species of Wikstroemia used in this study using MEGA 7 [59]. Two species, P. guajava (Myrtaceae; GenBank accession MN295360) and Gossypium australe (Malvaceae; GenBank accession AF057763), were included as outgroups. The alignment was trimmed using trimAL v1. 2[67] with the gappyout method in order to reduce the systematic errors produced by poor alignment. The optimal DNA substitution model for the ML analysis using the "Find Best DNA/ Protein Model (ML)" function embedded in MEGA 7 [59] was calculated to be the Kimura two-parameter (K2P) with discrete Gamma model (+G4) and invariant included (+I) (=K2P+G+I). ML analysis was performed using MEGA 7 [59] with 1,000 bootstrap replicates. The BI analysis was conducted with the method previously described [64].

Additional Information
Supplementary Materials Table S1: Information of introns and exons for protein-coding genes in six specie of Wikstroemia used in this study.